# The boundary condition of Hall-MHD model
**URL:** https://community.freefem.org/t/the-boundary-condition-of-hall-mhd-model/3484
**Category:** General Discussion
**Created:** [September 13, 2024, 1:59am UTC](https://community.freefem.org/t/the-boundary-condition-of-hall-mhd-model/3484 "2024-09-13T01:59:38Z")
**Posts on this page:** 12
**Page:** 1
### Author:  [@weixj](https://community.freefem.org/u/weixj)
#### Post date: [September 13, 2024, 1:59am UTC](https://community.freefem.org/t/the-boundary-condition-of-hall-mhd-model/3484/1 "2024-09-13T01:59:38Z")
Hello, everyone.
Now I’m using freefem++ to solve the incompressible hall-MHD equations, but I’ve encountered a problem regarding the boundary conditions of the magnetic field.
\mu\nabla\times (\nabla\times\mathbf{B})+(\mathbf{u}^n\cdot\nabla)\mathbf{B}^n-(\mathbf{B}^n\cdot\nabla)\mathbf{u}^n+\eta\nabla\times\left((\nabla\times\mathbf{B}^n)\times\mathbf{B}^n\right)=\mathbf{g}
\nabla\cdot\mathbf{B}=0
Before solving, assume that the velocity field \mathbf{u}^n,\mathbf{B}^n is already known (I decoupled the equation).
I use the boundary condition of magnetic
\mathbf{B}\times \mathbf{n} = \mathbf{0} or \mathbf{B}\times \mathbf{n} = \mathbf{B}\_e , where \mathbf{B}\_e is the analytical solution that I set, and \mathbf{n} is the unit external normal vector of the boundary. The domain is [0,1]^3.
```auto
solve Magnetic(Bx1, By1, Bz1, wx, wy, wz, solver = GMRES)
= int3d(Th)( mu*curl(Bx1, By1, Bz1)'*curl(wx, wy, wz) + mu*div(Bx1, By1, Bz1)*div(wx, wy, wz) )
+ int3d(Th)( UgradV(ux0, uy0, uz0, Bx0, By0, Bz0)'*[wx, wy, wz] - UgradV(Bx0, By0, Bz0, ux0, uy0, uz0)'*[wx, wy, wz]
+ eta*crossProduct(dy(Bz0)-dz(By0), dz(Bx0)-dx(Bz0), dx(By0)-dy(Bx0), Bx0, By0, Bz0)'*curl(wx, wy, wz)
- (g1t*wx + g2t*wy + g3t*wz))
;
```
with
```auto
macro Grad(ux, uy, uz) [(dx(ux)), (dy(ux)), (dz(ux)), (dx(uy)), (dy(uy)), (dz(uy)), (dx(uz)), (dy(uz)), (dz(uz))] //
macro div(ux, uy, uz) ( dx(ux) + dy(uy) + dz(uz) ) //
macro UgradV(ux, uy, uz, vx, vy, vz) [(ux)*(dx(vx)) + (uy)*(dy(vx)) + (uz)*(dz(vx)), (ux)*(dx(vy)) + (uy)*(dy(vy)) + (uz)*(dz(vy)), (ux)*(dx(vz)) + (uy)*(dy(vz)) + (uz)*(dz(vz))] //
macro curl(Bx, By, Bz) [dy(Bz) - dz(By), dz(Bx) - dx(Bz), dx(By) - dy(Bx)] //
macro crossProduct(ux, uy, uz, vx, vy, vz) [(uy)*(vz) - (uz)*(vy), (uz)*(vx) - (ux)*(vz), (ux)*(vy) - (uy)*(vx)] //
```
But, the result is wrong. Maybe I should add a stabilization term, but I do not how to do
---
### Author:  [@fb77](https://community.freefem.org/u/fb77)
#### Post date: [September 13, 2024, 8:35am UTC](https://community.freefem.org/t/the-boundary-condition-of-hall-mhd-model/3484/2 "2024-09-13T08:35:43Z")
To be consistent with the boundary condition B\times n=0 you should add a penalty boundary term
`+int2d(Th)(1.e15*crossProduct(Bx1, By1, Bz1, N.x, N.y, N.z)'*crossProduct(wx, wy, wz, N.x, N.y, N.z))`
For a non homogeneous condition B\times n = B\_e you should add also
`-int2d(Th)(1.e15*[Bex,Bey,Bez]'*crossProduct(wx, wy, wz, N.x, N.y, N.z))`
---
### Author:  [@weixj](https://community.freefem.org/u/weixj)
#### Post date: [September 13, 2024, 8:49am UTC](https://community.freefem.org/t/the-boundary-condition-of-hall-mhd-model/3484/3 "2024-09-13T08:49:49Z")
I understand. But is there any requirement for the penalty parameter? I see you gave `1.e15`
---
### Author:  [@fb77](https://community.freefem.org/u/fb77)
#### Post date: [September 13, 2024, 1:08pm UTC](https://community.freefem.org/t/the-boundary-condition-of-hall-mhd-model/3484/4 "2024-09-13T13:08:37Z")
The value is not important, it needs just to be much larger than the other terms, from 1.e8 to 1.e30 should do the job. The standards of FFEM uses tgv=1.e30.
---
### Author:  [@weixj](https://community.freefem.org/u/weixj)
#### Post date: [September 14, 2024, 7:35am UTC](https://community.freefem.org/t/the-boundary-condition-of-hall-mhd-model/3484/5 "2024-09-14T07:35:20Z")
OK, I understan.
Finally, I need to solve a Poisson equation for pressure
```auto
solve poisson(p, pp, solver = CG) = int3d(Th)( dx(p)*dx(pp) + dy(p)*dy(pp) + dz(p)*dz(pp) )
- int3d(Th)( (dx(p0)*dx(pp) + dy(p0)*dy(pp) + dz(p0)*dz(pp)) - div(ux, uy, uz)*pp/dt )
+ int2d(Th)( 1.e-8*p*pp );
```
Can I set the parameter to 1.e-8?
---
### Author:  [@fb77](https://community.freefem.org/u/fb77)
#### Post date: [September 14, 2024, 10:00am UTC](https://community.freefem.org/t/the-boundary-condition-of-hall-mhd-model/3484/6 "2024-09-14T10:00:46Z")
Yes it should work with 1.e8 (not 1.e-8). Alternatively you can use the built-in command
`+on(1,2,3,4,5,6,p=0)`
1,2,3,4,5,6 refers to the boundary labels. This built-in command is fully equivalent to
`+int2d(Th,1,2,3,4,5,6)(tgv*p*pp)`
---
### Author:  [@weixj](https://community.freefem.org/u/weixj)
#### Post date: [September 14, 2024, 1:07pm UTC](https://community.freefem.org/t/the-boundary-condition-of-hall-mhd-model/3484/7 "2024-09-14T13:07:01Z")
It should be 1e-8 (I try to set 1.e8, but it is wrong)
This is my equation, which has the artificial Newman boundary condition.
Now there is another problem.The convergence order of the L2 norm of the magnetic field I calculated is not correct (as the space is divided finer, the order decreases), but the convergence order of the H1 semi norm is no problem. I don’t know if this is due to the code.
```auto
h = [1./2, 1./4, 1./8, 1./12, 1./14, 1./16, 1./18, 1./20, 1./22, 1./24];
Convergence B_L2: 0 2.937716417 3.003305649 2.935914756 2.751860597
2.933047063 2.438786889 1.90737733 4.622726824 0.3392563845
Convergence B_H1_semi: 0 1.885937569 1.974759866 1.985809507 1.985307756
1.983228133 1.998874224 1.943471302 2.039907908 2.000346527
```
---
### Author:  [@fb77](https://community.freefem.org/u/fb77)
#### Post date: [September 15, 2024, 11:43am UTC](https://community.freefem.org/t/the-boundary-condition-of-hall-mhd-model/3484/8 "2024-09-15T11:43:53Z")
The 1.e8 is for homogeneous Dirichlet boundary condition. If you want Neumann, you just don’t put the penalty term!
About the order of convergence I don’t know, but your results look pretty good.
---
### Author:  [@Alice](https://community.freefem.org/u/Alice)
#### Post date: [December 16, 2024, 4:48am UTC](https://community.freefem.org/t/the-boundary-condition-of-hall-mhd-model/3484/9 "2024-12-16T04:48:21Z")
Dear Professor,
For the case non homegenous case B \times n=B\_e, why should add both
`+int2d(Th)(1.e15*crossProduct(Bx1, By1, Bz1, N.x, N.y, N.z)'*crossProduct(wx, wy, wz, N.x, N.y, N.z))`
`-int2d(Th)(1.e15*[Bex,Bey,Bez]'*crossProduct(wx, wy, wz, N.x, N.y, N.z)) ` ?
Why is adding only `-int2d(Th)(1.e15*[Bex,Bey,Bez]'*crossProduct(wx, wy, wz, N.x, N.y, N.z)) ` not enough?
---
### Author:  [@fb77](https://community.freefem.org/u/fb77)
#### Post date: [December 16, 2024, 10:46am UTC](https://community.freefem.org/t/the-boundary-condition-of-hall-mhd-model/3484/10 "2024-12-16T10:46:20Z")
Because both terms together give
10^{15}\*\int\_\Gamma(B\times N-B\_e)\cdot(w\times N)
Hence putting this in the variational formulation implies that
\int\_\Gamma(B\times N-B\_e)\cdot(w\times N)\simeq 0
for all test w.
We have to assume that B\_e\cdot N=0 otherwise the boundary condition is impossible to satisfy. Then since w\times N is an arbitrary vector on the boundary, orthogonal to N, it gives that B\times N-B\_e\simeq 0 on \Gamma.
After that it remains the variational formulation (interior terms) for all test w
such that w\times N=0 on \Gamma.
---
### Author:  [@Alice](https://community.freefem.org/u/Alice)
#### Post date: [December 17, 2024, 4:52am UTC](https://community.freefem.org/t/the-boundary-condition-of-hall-mhd-model/3484/11 "2024-12-17T04:52:56Z")
Dear Professor,
What should we do for the case n \times B =n \times B\_e?
---
### Author:  [@fb77](https://community.freefem.org/u/fb77)
#### Post date: [December 17, 2024, 7:48am UTC](https://community.freefem.org/t/the-boundary-condition-of-hall-mhd-model/3484/12 "2024-12-17T07:48:59Z")
Just replace B\_e by B\_e\times N in the code!