Re: [ff3d-users] magnetic modeling

[Stephane Del Pino]

Le lundi 25 septembre 2006 09:50, Schmidt Beat a écrit :
> Hi,
Hello.

> I'm trying to use ff3d for 3D magnetic modeling. I've started out with a
> simple cylindrical magnet (current in a cylindrical sheet) and would like
> to calculate the vector potential A, as a first step. The first problem I
> have encountered was that the solver does not seem to accept the following
> form of the Robin boundary condition :
>
> C*Ax+dnu(Ax)=0
>
> Only Ax+dnu(Ax)=0 is accepted.
In fact there is a stupid "bug" here: C must be a function (note that this 
does not change the performances) ... 

> The second problem is that even when the boundary condition is given in
> this second form, the solver stops in the middle of the calculation and
> returns an error.
Which version of ff3d are you using?
Which error do you get?

I tried your data and fixed a few mistakes, and it worked fine for me ... with 
the cvs release.

Here are the changed for the ff3d file:
> double pi = 3.1415926536;
> double mu0 = 4 * pi *10e-7;
> double R=0.5; //rayon de l'aimant
> double C=0.25; //inverse du rayon de la sphère limitante R=4
function C = 0.25;
> double j0 = 1; //normalisation du courant d'aimantation
>
> scene S = pov("test.pov");
> mesh M = structured(<80,80,80>,(-4,-4,-4),(4,4,4));
You do not need such a fine mesh, but you nee your mesh to have vertices on 
the (Ox) axis. I propose <21,21,21> ...
> domain O = domain(S, inside(<1,1,1>));
> function jx = - one(<1,0,0>)*z/0.5;
> function jz = one(<1,0,0>)*x/0.5;
>
> solve(Ax,Ay,Az)in O by M method(type = penalty)
> {
>     pde(Ax)
>         dx(dx (Ax)) +dy(dy (Ax)) + dz(dz (Ax)) = mu0*jx;
Here: 
1) you'd better use the div(grad(Ax)) operator
2) you should write -div(grad(Ax)) = -mu0*jx;
If you don't, the operator you are solving will not be positive definite so 
conjugate gradient will not work and you will have to use bicg for 
instance ...
>         Ax+dnu(Ax)=0 on <1,1,1>;
C*Ax+dnu(Ax)=0 on <1,1,1>;

And the same for the following:
>     pde(Ay)
>         dx(dx (Ay)) +dy(dy (Ay)) + dz(dz (Ay)) = 0;
>         Ay+dnu(Ay)=0 on <1,1,1>;
>
>     pde(Az)
>         dx(dx (Az)) +dy(dy (Az)) + dz(dz (Az)) = mu0*jz;
>         Az+dnu(Az)=0 on <1,1,1>;
>
> }

Best regards,
Stephane.