For the magnetic field, the boundary conditions I chose are: \mathbf{B}\times\mathbf{n}|_{\partial \Omega}=0,\quad \nabla\cdot \mathbf{B}|_{\partial \Omega}=0
Regarding the boundary conditions \mathbf{B}\times\mathbf{n}|_{\partial \Omega}=0, I know the correct approach is:
However, I applied a divergence-free condition to the boundary \nabla\cdot \mathbf{B}|_{\partial \Omega}=0. How is this specifically implemented in the FreeFem++ code?
The issue is indeed to have an appropriate variational formulation, and to choose a good finite element space.
Here since you already set {\mathrm{div} B}=0 in \Omega, you don’t really need to impose {\mathrm{div} B}=0 in \partial\Omega because the latter follows from the former.
If you want nevertheless to enforce such condition you could add a term -\nabla\mathrm{div} B in the B equation. Then you have to take a space of the type P1 for B (you need \nabla B in L^2).
On the contrary if you don’t put such term you can take a space with only \nabla\times B\in L^2, such as RT03d or Edge03d, this should be more efficient.
I notice that you don’t have to really impose {\mathrm{div} B}=0 in \Omega because if it holds at t=0, the B equation implies that it will hold for all t>0 (but maybe not exactly at the discrete level).
Thank you very much for your answer, Professor. I understand what you mean. For the magnetic field, I chose a second-order Nedelec element \subset H_0(\text{curl},\Omega).
Regarding why I set \nabla\cdot \mathbf{B}=0 on \partial\Omega, the specific reasons are as follows (using a first-order time scheme as an example).
If set \mathbf{H}_{h} =\nabla(\nabla\cdot \mathbf{B}_{h}^{n+1}), under the boundary conditions \nabla\cdot\mathbf{B}_{h}^{n+1}|_{\partial \Omega}=0, we can obtain
Your equation on B is second-order, and you would like to impose a condition on \nabla\cdot B on the boundary. This is a condition on the first-order derivatives of B. Such condition is usually imposed in the weak sense by a suitable variational formulation (think of the case of the Laplace equation). Boundary conditions in the strong sense are imposed only when no derivative of the unknown is involved (as your condition B\times n=0).
Your scheme at the continuous level is nice, but looks impossible to be converted in a finite element framework. The reason is that you choose a space that is curl conforming, but no div conforming. As a consequence there are some jumps of B\cdot n between elements. Even worse, there are some jumps of \nabla\cdot B between elements, and thus you cannot take \nabla(\nabla\cdot B) as test function. To make your estimate work you would need a space that is H^2 conforming, which is out of range.
The inequality you would like to achieve would imply that \nabla\cdot B=0 if it holds initially.
This property is the only thing that is relevant, but looks unrealizable for the same reason as above: the finite element space is not div conforming.
I see only these possibilities:
– Use a H^1 conforming space like P1, and include the -\eta\nabla(\nabla\cdot B) term, which means put a \eta(\nabla\cdot B)(\nabla\cdot H) in the variational formulation. This will ensure a dissipation of \nabla\cdot B.
– Use a discontinuous Galerkin formulation to take into account the jumps (but looks really difficult)
– Use a curl conforming space as is your initial choice, and do nothing for the discrete div constraint.
