“Problem” and “varf” format

General Discussion
1

Hi, everyone!
I have a three-dimensional model：

E+u\times B-(1/Rm)* \curl B=g             in \Omega
(\partial B)/(\partial t)+\curl E=0           in \Omega

The boundary condition is

B\cdot n=0,     E\times n=0    On \partial 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 )-(B_h^n,\curl F)=(g,F),
(1/\tau)*(B_h^n-B_h^{n-1},C)+(\curl E_h^n,C)=0.

In the code, I do the decoupling operation, which involves the exact solution with respect to the values of u. I used two types of code to solve these two models, but the convergence order obtained using “Problem” was correct, while it was incorrect when changed to “varf”.

Blockquote
solver_EB_3D.edp (6.2 KB)
solver_EB3D_save.edp (6.3 KB)

2

Hi,
In your right-hand side you miss the boundary values for B. You can include them by

     real[int] bon = a2(0,WhPh);
     b=b+bon;

When you set [E1[], B1[]] = sol, the dof are not in the correct order. A correct way is

    WhPh [Q1,Q2,Q3,Z1,Z2,Z3];
     Q1[]=sol;
    [E1,E2,E3]=[Q1,Q2,Q3];[B1,B2,B3]=[Z1,Z2,Z3];

Another difference is that in one file you have [u1,u2,u3]=(Et(0.)+1.)*u1Exact; and in the other [u1,u2,u3]=(Et(tnext)+1.)*u1Exact;
full file with the corrections (not the last one with Et(0.) or Et(tnext))
solver_EB3D_save.edp (6.5 KB)

3

Thank you very much for your reply!
Based on the code you returned, I have been able to successfully run the correct result. Since the velocity u, the electric field E and the magnetic field B are solved in coupling, I need to define to combine the space ‘Uh’

  fespace Uh(Th,[P1,P1,P1], Edge03d, RT03d);
  Uh [Q1,Q2,Q3,Z1,Z2,Z3,N1,N2,N3];
...
real[int] sol = A^-1 * b;
      Q1[]=sol;
    [u1,u2,u3]=[Q1,Q2,Q3];[E1,E2,E3]=[Z1,Z2,Z3]; [B1,B2,B3]=[N1,N2,N3];

However, defining ‘Uh’ this way will result in an error. Do you have any good methods？
Looking forward to your reply！

4

You can simply set
fespace Uh(Th, [P1,P1,P1,Edge03d, RT03d]);
The important thing is to put Q1[]=sol where [Q1,...] is declared in the space Uh involved in the matrix and rhs definitions matrix A1 = a1(Uh,Uh); real[int] b = rhs(0, Uh); so that the degrees of freedom are correctly interpreted.

5

Thank you for your suggestion. I can run it correctly now！