The Edge03d element is used for numerical solution in freefem

[Dainy-Jia]

I have a three-dimensional model：
E+u\times B-\curl B=g in \Omega, E \times n=0 on \partal \Omega
here, u represents the velocity field, E represents the electric field, and B represents the magnetic field.
The discrete format of this model is
(E_h^n,F)+(u_h^n \times B_h^{n-1},F )-(\curl B_h^{n},F)=(g,F),
In the code, I do the decoupling operation, which involves the exact solution with respect to the values of u and B. But unfortunately, I didn’t get the convergence order right.
load “Element_Mixte3d”
load “msh3”
load “medit”
int MaxIter = 2;
real[int] errL2(MaxIter), orderL2(MaxIter);

macro C(t) 1./picos(pit)-1.5t //EOM
macro D(t) sin(pit)+1.5 //EOM

real cpu = clock();
real mu=1.;
real muc=2.;
real Rm=1.;

func real[int] uExact(real t){
real[int] uExactvec(3);
uExactvec(0) = (exp(-t)+1.)sin(pix)sin(piy)sin(piz);
uExactvec(1) = 0.;
uExactvec(2) = (exp(-t)+1.)(z^2-z)sin(pix)sin(piy);
return uExactvec;
};
func real[int] udxExact(real t){
real[int] udxExactvec(3);
udxExactvec(0) = (exp(-t)+1.)picos(pix)sin(piy)sin(piz);
udxExactvec(1) = 0.;
udxExactvec(2) = (exp(-t)+1.)picos(pix)(z^2-z)sin(piy);
return udxExactvec;
};
func real[int] udyExact(real t){
real[int] udyExactvec(3);
udyExactvec(0) = (exp(-t)+1.)picos(piy)sin(pix)sin(piz);
udyExactvec(1) = 0.;
udyExactvec(2) = (exp(-t)+1.)picos(piy)*(z^2-z)sin(pix);
return udyExactvec;
};

func real[int] udzExact(real t){
real[int] udzExactvec(3);
udzExactvec(0) = (exp(-t)+1.)picos(piz)sin(pix)sin(piy);
udzExactvec(1) = 0.;
udzExactvec(2) = (exp(-t)+1.)(2.*z-1.)sin(pix)sin(piy);
return udzExactvec;
};

func real[int] E1Exact(real t){
real[int] E1Exactvec(3);
E1Exactvec(0) = (sin(pit)+1.5)(y^2-y)sin(piz);
E1Exactvec(1) = (sin(pit)+1.5)(x^2-x)sin(piz);
E1Exactvec(2) = (sin(pit)+1.5)(x^2-x)sin(piy);
return E1Exactvec;
};
macro EExact(t) [(sin(pit)+1.5)(y^2-y)sin(piz), (sin(pit)+1.5)(x^2-x)sin(piz),(sin(pit)+1.5)(x^2-x)sin(piy)] //EOM
func real[int] EdxExact(real t){
real[int] EdxExactvec(3);
EdxExactvec(0) = 0.;
EdxExactvec(1) = (sin(pit)+1.5)(2.x-1.)sin(piz);
EdxExactvec(2) = (sin(pit)+1.5)(2.x-1.)sin(piy) ;
return EdxExactvec;
};
func real[int] EdyExact(real t){
real[int] EdyExactvec(3);
EdyExactvec(0) = (sin(pit)+1.5)(2.y-1.)sin(piz);
EdyExactvec(1) = 0.;
EdyExactvec(2) = (sin(pit)+1.5)(x^2-x)picos(piy);
return EdyExactvec;
};

func real[int] EdzExact(real t){
real[int] EdzExactvec(3);
EdzExactvec(0) = (sin(pit)+1.5)(y^2-y)picos(piz);
EdzExactvec(1) = (sin(pit)+1.5)(x^2-x)picos(piz);
EdzExactvec(2) = 0.;
return EdzExactvec;
};

func real[int] B1Exact(real t){
real[int] B1Exactvec(3);
B1Exactvec(0) = (1./picos(pit)-1.5t)pi(x^2-x)(cos(piy)-cos(piz));
B1Exactvec(1) = (1./picos(pit)-1.5t)(pi*(y^2-y)cos(piz)-(2.x-1.)sin(piy));
B1Exactvec(2) = (1./picos(pit)-1.5t)(2.x-2.y)sin(piz);
return B1Exactvec;
};
macro BExact(t) [(1./picos(pit)-1.5t)pi(x^2-x)(cos(piy)-cos(piz)), (1./picos(pit)-1.5t)(pi(y^2-y)cos(piz)-(2.x-1.)sin(piy)),(1./picos(pit)-1.5t)(2.x-2.y)sin(piz) ] // EOM
func real[int] BdxExact(real t){
real[int] BdxExactvec(3);
BdxExactvec(0) = (1./picos(pit)-1.5t)pi(2.x-1.)(cos(piy)-cos(piz));
BdxExactvec(1) = (1./picos(pit)-1.5t)(-2.)sin(piy);
BdxExactvec(2) = (1./picos(pit)-1.5t)(2.)sin(piz);
return BdxExactvec;
};
func real[int] BdyExact(real t){
real[int] BdyExactvec(3);
BdyExactvec(0) = (1./picos(pit)-1.5t)(-pi^2)(x^2-x)sin(piy);
BdyExactvec(1) = (1./picos(pit)-1.5t)(pi(2.y-1.)cos(piz)-pi(2.x-1.)cos(piy));
BdyExactvec(2) = (1./picos(pit)-1.5t)*(-2.sin(piz));
return BdyExactvec;
};

func real[int] BdzExact(real t){
real[int] BdzExactvec(3);
BdzExactvec(0) = (1./picos(pit)-1.5t)pi^2(x^2-x)sin(piz);
BdzExactvec(1) = (1./picos(pit)-1.5t)(-pi^2)(y^2-y)sin(piz);
BdzExactvec(2) = (1./picos(pit)-1.5*t)pi(2.*x-2.*y)cos(piz);
return BdzExactvec;
};

func real[int] BBcd(real t){
real[int] BBcdvec(3);
BBcdvec(0) = (1./picos(pit)-1.5t)pi(x^2-x)(cos(piy)-cos(piz));
BBcdvec(1) = (1./picos(pit)-1.5t)(pi*(y^2-y)cos(piz)-(2.x-1.)sin(piy));
BBcdvec(2) = (1./picos(pit)-1.5t)*(2.*x-2.*y)sin(piz);
return BBcdvec;
};

func real[int] EBcd(real t){
real[int] EBcdvec(3);
EBcdvec(0) = (sin(pit)+1.5)(y^2-y)sin(piz);
EBcdvec(1) = (sin(pit)+1.5)(x^2-x)sin(piz);
EBcdvec(2) = (sin(pit)+1.5)(x^2-x)sin(piy);
return EBcdvec;
};

func real[int] g(real t){
real[int] gvec(3);
gvec(0) = E1Exact(t)(0) + uExact(t)(1) * B1Exact(t)(2) - uExact(t)(2) * B1Exact(t)(1)
- 1. / Rm * (BdyExact(t)(2) - BdzExact(t)(1));
gvec(1) = E1Exact(t)(1) + uExact(t)(2) * B1Exact(t)(0) - uExact(t)(0) * B1Exact(t)(2)
- 1. / Rm * (BdzExact(t)(0) - BdxExact(t)(2));
gvec(2) = E1Exact(t)(2) + uExact(t)(0) * B1Exact(t)(1) - uExact(t)(1) * B1Exact(t)(0)
- 1. / Rm * (BdxExact(t)(1) - BdyExact(t)(0));
return gvec;
};

//real cpu = clock();
// mesh iteration
for (int i = 0; i < MaxIter; i++){
int n=pow(2, i + 2);
mesh3 Th = cube(n, n, n, [x, y, z]);
fespace Vh(Th, P1);
fespace Ph(Th, RT03d);
fespace Wh(Th, Edge03d);
Wh [E1, E2,E3],[F1,F2,F3];
Ph [B1,B2,B3], [B1old,B2old,B3old];
Vh u1old, u2old, u3old, u1,u2,u3;

int nt = pow(2*(i+2),2) / 4;
real tfinal = 1./4.;
real dt = tfinal / nt, alpha = 1. / dt;
real tnow = 0;
u1old = uExact(0.)(0);
u2old = uExact(0.)(1);
u3old = uExact(0.)(2);
[B1old,B2old,B3old] = BExact(0.);
problem solverE([E1,E2,E3],[F1,F2,F3], solver=CG) = 
                int3d(Th)(E1*F1+E2*F2+E3*F3)  
               + int3d(Th)((u2*B3old-u3*B2old)*F1+(u3*B1old-u1*B3old)*F2+(u1*B2old-u2*B1old)*F3) 
               - int3d(Th)((1./Rm)*(F1*(dy(B3)-dz(B2))+F2*(dz(B1)-dx(B3))+F3*(dx(B2)-dy(B1))))     
                - int3d(Th)(g(tnow+dt)(0)*F1+g(tnow+dt)(1)*F2+g(tnow+dt)(2)*F3)         
                + on(1, 2, 3, 4, 5, 6, E1=0.,E2=0.,E3=0.);

for (int i = 1; i <= nt; i++) {
real tnext = tnow + dt;  
u1=uExact(tnext)(0);
u2=uExact(tnext)(1);
u3=uExact(tnext)(2);
[B1,B2,B3]=BExact(tnext);
solverE;                         
tnow = tnext;            
[B1old,B2old,B3old]=BExact(tnow);

}

// E_error
errL2[i] = sqrt(int3d(Th)((E1 - E1Exact(tfinal)(0)) ^ 2+ (E2 - E1Exact(tfinal)(1)) ^ 2+ (E3 - E1Exact(tfinal)(2)) ^ 2));

};

for (int i = 1; i < MaxIter; i++){
orderL2[i] = log(errL2[i - 1] / errL2[i]) / log(2.);
};

cout << “================OUTPUT===============” << endl;
cout << "CPU time = " << clock() - cpu << endl;
cout << "errors and convergence rates: " << endl;
cout << "E = " << endl;
cout << “h\t L2 error\t order” << endl;
for (int i = 0; i < MaxIter; i++){
cout << “1/” << 2 ^ (i + 2) << “:\t”
<< errL2[i] << “\t” << orderL2[i] << “\t\t” << endl;
};

[Monirul25 (Monirul)]

Hello, dear Dainy-Jia, Can you please provide your code as your_code.edp format so that it can run directly.

[Dainy-Jia]

Sure. Thank you for your reply!

Blockquote
solver_E_3D(1).edp (6.8 KB)

[fb77 (François Bouchut)]

Hello,
Your space Ph for the magnetic field chosen as RT03d seems inappropriate. I have taken P0 instead.
Then the curl B term needs to be integrated by parts, and the boundary condition E\times n=0 needs to be set correctly.
Another thing is that all macros need parenthesis () around the argument because they do literal replacement of the argument. When the argument is tnow+dt for example, without parenthesis it would not distribute the factor to each term tnow and dt.
solver_E_3D.edp (7.1 KB)

[Dainy-Jia]

Thank you very much for your reply！
Why does the magnetic field here use [P0,P0,P0] elements instead of RT03d elements?

In the model I presented above, these two equations are solved simultaneously. E belongs to the H(curl) space, and B belongs to the H(div) space.
The discrete format is as follows：

So naturally, the electric field (E_h^n) is approximated using Edge03d(the lowest-order edge element), and the magnetic field (B_h^n) is approximated using RT03d(the lowest-order face element).
By decoupling, I separately calculated the discrete format regarding the magnetic field as follows：

The code is

Blockquote
slover_B_3D.edp (4.8 KB)
The calculation result is:
T=1/4; solver_B_3D
================OUTPUT===============
CPU time = 645.957
errors and convergence rates:
B =
h L2 error order divBh
1/4: 0.0448658 0 6.77796e-17
1/8: 0.0223948 1.00245 3.09493e-16
1/16: 0.0111921 1.00068 6.56293e-16
This is a good result, but why does the order of convergence get wrong when we decouple the electric field? Looking forward to your reply！

[fb77 (François Bouchut)]

Dear Dainy-Jia,
I don’t know why RT03d does not work for B in the equation on E. Maybe one has to take finer meshes to see the convergence? Or maybe \nabla\times B, which is involved in the equation on E, is not well approximated when B in RT03d?
For the equation on B, RT03d works, and also P0, and Edge03d (I have tried with your code).
Thus I think that P0 for B will work for the coupled problem. Indeed with P0, the variational formulation on B is equivalent to setting directly
\frac{B^n_h-B_h^{n-1}}{\tau}+\nabla\times E_h^n=0, because when E is in Edge03d, \nabla\times E is in P0. Thus you can replace B_h^n in the equation on E by this expression and find a variational formulation for E_h^n alone, involving (\nabla\times E_h^n,\nabla\times F).
This formulation is good for E in Edge03d, and you have a single problem on E instead of a coupled problem in (E,B) to solve.
You could also do the converse: B in Edge03d and E in P0, and replace the value of E from the first equation into the second one (but then you have probably to make a correction with a “pressure type” Lagrange multiplier to get \nabla\cdot B=0 weakly, and a regularization term (\nabla\cdot B,\nabla\cdot C) ).

[Dainy-Jia]

Thank you for your reply and patient explanation. However, I still have one question. Why is the boundary condition chosen as:
+int2d(Th,1,3)(1.e30*(E1F1+E3F3))
+int2d(Th,2,4)(1.e30*(E2F2+E3F3))
+int2d(Th,5,6)(1.e30*(E1F1+E2F2))
Looking forward to your reply！

[fb77 (François Bouchut)]

This is for E\times n=0, knowing that
labels 1,3 correspond to n=(0,\pm 1,0)
labels 2,4 correspond to n=(\pm 1,0,0)
labels 5,6 correspond to n=(0,0,\pm1)
thus
for labels 1,3 we have to set E1=0 and E3=0,
for labels 2,4 we have to set E2=0 and E3=0,
for labels 5,6 we have to set E1=0 and E2=0
This is done by penalty.

[frederichecht (Frédéric Hecht)]

Please do not put 1e30 as a penalty term i integral because you have term off diagonal and it is true penalisation method,
1e6 is OK , I think. Remark on(1,3)(E1=0,E2=0,E3=0) impose E.tangent = 0 and tangent is (±1,0,0)

[fb77 (François Bouchut)]

Dear Frédéric,
If I understand well, you say that because of the particular structure of the space Edge03d, setting E1=0,E2=0,E3=0 indeed only imposes E.tangent=0.
This is good news, it simplifies the boundary statements!