Magnetostatics inconsistent sign for current and field

I downloaded the 2D example of the magnetostatic problem at 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,

Shouldn’t it read Bx=-dy(Az), or By=dx(Az)?

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,

real Mu0 = 4.pi1.e-7; //Vacuum magnetic permeability
real MuC = 1.25e-6; //Copper magnetic permeability
real J0 = 500.e4; //Current density

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., BoxRadiusnn/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
border b4(t=0., 1.){x=MagnetExternalRadius-(MagnetExternalRadius-MagnetInternalRadius)t; y=MagnetHeight/2.;};
border b5(t=0., 1.){x=MagnetInternalRadius; y=MagnetHeight/2.-MagnetHeight

border b6(t=0., 1.){x=-MagnetInternalRadius-(MagnetExternalRadius-MagnetInternalRadius)t; y=-MagnetHeight/2.;};
border b7(t=0., 1.){x=-MagnetExternalRadius; y=-MagnetHeight/2.+MagnetHeight
border b8(t=0., 1.){x=-MagnetExternalRadius+(MagnetExternalRadius-MagnetInternalRadius)t; y=MagnetHeight/2.;};
border b9(t=0., 1.){x=-MagnetInternalRadius; y=MagnetHeight/2.-MagnetHeight

int nnH = max(2., 10.nnMagnetHeight);
int nnR = max(2., 10.nn(MagnetExternalRadius-MagnetInternalRadius));
mesh Th = buildmesh(
+ 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;

func Pk = P2;
fespace Ah(Th, Pk);
Ah Az;

macro Curl(A) [dy(A), -dx(A)] //

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

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

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(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) ;

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).