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;
};