# Magnetostatics inconsistent sign for current and field **URL:** \ **Category:** General Discussion\ **Created:** [September 6, 2019, 7:29am UTC](https://community.freefem.org/t/magnetostatics-inconsistent-sign-for-current-and-field/125 "2019-09-06T07:29:05Z")\ **Posts on this page:** 4\ **Page:** 1
**Author:** ![alberto](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/alberto/32/255_2.png) [@alberto](https://community.freefem.org/u/alberto)\ **Post date:** [September 6, 2019, 7:29am UTC](https://community.freefem.org/t/magnetostatics-inconsistent-sign-for-current-and-field/125/1 "2019-09-06T07:29:06Z")
Hello, I downloaded the 2D example of the magnetostatic problem at [https://modules.freefem.org/modules/magnetostatic/](https://modules.freefem.org/modules/magnetostatic/) and I have added a few lines to compute By = dy(Az) and plot its values in 2D. According to the sign of the currents the value of By at the center should be positive, but it turns out to be negative when i run the script. Also, if I plot the current obtained from J = dx(Hy) - dy(Hx) to check for self-consistency, i see that the currents are exactly opposite with respect to the ones set at the beginning. Also, I think the signs of the weak form of the equation are correct. Could somebody explain this behavior to me? Thank you very much, Adele ---
**Author:** ![denis.crete](https://avatars.discourse-cdn.com/v4/letter/d/9d8465/32.png) [@denis.crete](https://community.freefem.org/u/denis.crete)\ **Post date:** [January 4, 2020, 10:02pm UTC](https://community.freefem.org/t/magnetostatics-inconsistent-sign-for-current-and-field/125/2 "2020-01-04T22:02:16Z")
> [@alberto](#): > > By = dy(Az) Shouldn’t it read Bx=-dy(Az), or By=dx(Az)? Denis ---
**Author:** ![alberto](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/alberto/32/255_2.png) [@alberto](https://community.freefem.org/u/alberto)\ **Post date:** [January 7, 2020, 8:38am UTC](https://community.freefem.org/t/magnetostatics-inconsistent-sign-for-current-and-field/125/3 "2020-01-07T08:38:51Z")
Dear Denis, thank you for pointing out the typo in my message. Because B = curl(A) the correct way of writing is Bx = dy(Az). The other component is By = - dx(Az). I copied my script and I pasted it below. If you run it, you will se that the currents computed to check self-consistency are opposite with respect to the currents imposed as boundary conditions. I don’t know why it comes out this way, they should be exactly the same. Thank you, Alberto //Parameters real Mu0 = 4._pi_1.e-7; //Vacuum magnetic permeability real MuC = 1.25e-6; //Copper magnetic permeability real J0 = 500.e4; //Current density //Mesh int nn = 25; //Mesh quality real BoxRadius = 2.; //“Inifite” boundary real MagnetHeight = 1.; //Magnet height real MagnetInternalRadius = 0.1; //Magnet internal radius real MagnetExternalRadius = 0.3; //Magnet external radius int BoxWall = 10; int nnL = max(2., BoxRadius_nn/2.); border b1(t=0., 2.pi){x=BoxRadiuscos(t); y=BoxRadius_sin(t); label=BoxWall;}; border b2(t=0., 1.){x=MagnetInternalRadius+(MagnetExternalRadius-MagnetInternalRadius)_t; y=-MagnetHeight/2.;}; border b3(t=0., 1.){x=MagnetExternalRadius; y=-MagnetHeight/2.+MagnetHeight_t;}; border b4(t=0., 1.){x=MagnetExternalRadius-(MagnetExternalRadius-MagnetInternalRadius)_t; y=MagnetHeight/2.;}; border b5(t=0., 1.){x=MagnetInternalRadius; y=MagnetHeight/2.-MagnetHeight_t;}; border b6(t=0., 1.){x=-MagnetInternalRadius-(MagnetExternalRadius-MagnetInternalRadius)_t; y=-MagnetHeight/2.;}; border b7(t=0., 1.){x=-MagnetExternalRadius; y=-MagnetHeight/2.+MagnetHeight_t;}; border b8(t=0., 1.){x=-MagnetExternalRadius+(MagnetExternalRadius-MagnetInternalRadius)_t; y=MagnetHeight/2.;}; border b9(t=0., 1.){x=-MagnetInternalRadius; y=MagnetHeight/2.-MagnetHeight_t;}; int nnH = max(2., 10._nn_MagnetHeight); int nnR = max(2., 10._nn_(MagnetExternalRadius-MagnetInternalRadius)); mesh Th = buildmesh( b1(nnL) + b2(nnR) + b3(nnH) + b4(nnR) + b5(nnH) + b6(-nnR) + b7(-nnH) + b8(-nnR) + b9(-nnH) ); int Coil1 = Th((MagnetExternalRadius+MagnetInternalRadius)/2., 0.).region; int Coil2 = Th(-(MagnetExternalRadius+MagnetInternalRadius)/2., 0.).region; int Box = Th(0., 0.).region; //Fespace func Pk = P2; fespace Ah(Th, Pk); Ah Az; //Macro macro Curl(A) [dy(A), -dx(A)] // //Current Ah Jz = J0\*(-1._(Th(x,y).region == Coil1) + 1._(Th(x,y).region == Coil2)); plot(Jz, nbiso=30, fill=true, value=true, cmm=“J”, wait=true); //Problem func Mu = Mu0 + (MuC-Mu0)\*((region==Coil1)+(region==Coil2)); func Nu = 1./Mu; varf vLaplacian (Az, AAz) = int2d(Th)( Nu \* Curl(Az)’ \* Curl(AAz) ) - int2d(Th)( Jz \* AAz ) + on(BoxWall, Az=0) ; matrix Laplacian = vLaplacian(Ah, Ah, solver=sparsesolver); real[int] LaplacianBoundary = vLaplacian(0, Ah); Az[] = Laplacian^-1 \* LaplacianBoundary; //Magnetic induction Ah Bx, By; Bx = dy(Az); By = -dx(Az); Ah B = sqrt(Bx^2 + By^2); //Magnetic field Ah Hx, Hy; Hx = Nu \* Bx; Hy = Nu \* By; Ah H = sqrt(Hx^2 + Hy^2); //Current Ah J = dx(Hy) - dy(Hx); //Lorentz force Ah Lx, Ly; Lx = -J \* By; Ly = J \* Bx; Ah L = sqrt(Lx^2 + Ly^2); //Plot plot(Az, nbiso=20, fill=true, value=true, cmm=“A”, wait=true); plot(B, nbiso=20, fill=true, value=true, cmm=“B”, wait=true); plot(Bx, nbiso=20, fill=false, value=true, cmm=“Bx”, wait=true) ; plot(By, nbiso=20, fill=true, value=true, cmm=“By”, wait=true) ; plot(H, nbiso=20, fill=true, value=true, cmm=“H”, wait=true); plot(L, nbiso=20, fill=true, value=true, cmm=“L”, wait=true); plot(J, nbiso=20, fill=true, value=true, cmm=“Jpost”, wait=true) ; ---
**Author:** ![denis.crete](https://avatars.discourse-cdn.com/v4/letter/d/9d8465/32.png) [@denis.crete](https://community.freefem.org/u/denis.crete)\ **Post date:** [January 14, 2020, 11:22pm UTC](https://community.freefem.org/t/magnetostatics-inconsistent-sign-for-current-and-field/125/4 "2020-01-14T23:22:20Z")
Hello, I suspect that there is a wrong sign in the definition of vLaplacian, as Rot(Rot A)=Rot(B)= mu0.J = Grad(Div A)-Nabla^2(A); with the (usual) gauge div A=0, the result is Laplacian(A)+mu0\*J=0. If I am right, - int2d(Th)(Jz \* AAz) should be +int2d(Th)(Jz \* AAz). HTH Denis