Hi everyone,
I’m confused with the results I’m getting with the RT1 space for the mixed formulation of the Poisson Equation. I adapted the LaplaceRT.edp code to have an analytical solution.
When I use the RT0&P0 pair, I get div u ~ 0, as expected. However, the RT1&P1 pair gives a huge value for div u, but when I project the div u into P1, the projection has a small error.
I’ve got a similar problem with [P2,P2]&P1 pair as well.
// Solve
// u + grad p = g
// div u = 0
load "Element_Mixte" // for RT1
mesh Th=square(50,50);
fespace Vh(Th,RT1);
fespace Ph(Th,P1);
// Exact solution //
real omega = 2.0*pi;
func u1exact = sin(omega*x)*cos(omega*y);
func u2exact = -cos(omega*x)*sin(omega*y);
func pexact = sin(omega*x)*sin(omega*y);
func dxpexact = omega*cos(omega*x)*sin(omega*y);
func dypexact = omega*sin(omega*x)*cos(omega*y);
func g1 = u1exact + dxpexact;
func g2 = u2exact + dypexact;
// Exact solution //
Vh [u1,u2],[v1,v2];
Ph p,q,qq;
solve laplaceMixte([u1,u2,p],[v1,v2,q],solver=UMFPACK) =
int2d(Th)( p*q*1e-10 + u1*v1 + u2*v2 - p*(dx(v1)+dy(v2)) + (dx(u1)+dy(u2))*q )
- int2d(Th) ( [g1,g2]'*[v1,v2])
+ int1d(Th,1,2,3,4)(pexact*(v1*N.x+v2*N.y))
// + on(1,2,3,4,u1=u1exact,u2=u2exact)
;
real avg = int2d(Th)(p)/Th.area;
p = p - avg; // mean value zero
plot(p,fill=1,wait=1,value=true);
q = dx(u1)+dy(u2);
plot(q,coef=0.05,wait=1,value=true);
real errorDivUL2 = sqrt(int2d(Th)( (dx(u1)+dy(u2))^2.0 ));
cout << "Before postprocessing: errorDivUL2 = " << errorDivUL2 << endl;
// Project divu
solve L2Proj(q,qq,solver=CG) =
int2d(Th)(q*qq) - int2d(Th)(qq*(dx(u1)+dy(u2)) );
plot(q,fill=1,wait=1,value=1,cmm="div u");
// Errors //
errorDivUL2 = sqrt(int2d(Th)( q^2.0 ));
cout << "After postprocessing: errorDivUL2 = " << errorDivUL2 << endl;
real errorUL2 = sqrt(int2d(Th)( (u1-u1exact)^2.0 + (u2-u2exact)^2.0 ));
real errorPL2 = sqrt(int2d(Th)( (p-pexact)^2.0 ));
cout << "errorUL2 = " << errorUL2 << endl;
cout << "errorPL2 = " << errorPL2 << endl;
If anyone could point out the issue, that would be very helpful.