[ff3d-users] Problem in Solving Laplace with Neumann BC's

[Juzar Thingna]

Hello Everyone,
I'm trying to solve the laplace equation on a rectangular cube. On y min I
set the Dirichlet BC's whereas on y max I want to set the value of the
co-normal derivative. Since the equation is a laplace equation the solution
along the y-direction should be a straight line. Since on the other surfaces
the normal derivative = 0 automatically. Here is a copy of my input file:

vector n=(50,70,50);
vector a=(0,-37,0);
vector b=(25,-7,25);
mesh M=structured(n,a,b);
scene S=pov("ferromagnet_1.pov");
domain O =domain(S,inside(<1,0,0>));
function dv=0.130;
solve(v) in O by M
{
pde(v)
div(grad(v))=0;
v=10.0 on M ymin;
dnu(v)= dv on M ymax;
};
mesh Z = tetrahedrize (O,M);
save(medit,"v_fm_1",v,Z);
save(medit,"v_fm_1",Z);
save(raw,"v_fm_1.dat",v,M);
vector n=(50,70,50);
vector a=(0,-37,0);
vector b=(25,-7,25);
mesh M=structured(n,a,b);
scene S=pov("ferromagnet_1.pov");
domain O =domain(S,inside(<1,0,0>));
function dv=0.130;
solve(v) in O by M
{
pde(v)
div(grad(v))=0;
v=10.0 on M ymin;
dnu(v)= dv on M ymax;
};
mesh Z = tetrahedrize (O,M);
save(medit,"v_fm_1",v,Z);
save(medit,"v_fm_1",Z);
save(raw,"v_fm_1.dat",v,M);

But when I plot my solution as a function of y I notice that there are many
discrepancies. I have attaced my plot here. In the vicinity of ymax I get
the derivative to be = dv but the solution is not a straight line as
expected. I earlier thought it could be because of the roughness of my grid,
but even when my grid is finer I have the same problems.
Is there a way to get around this?

Regards,
Juzar Thingna
Center for Computation Science and Engineering
National University of Singapore

[Stephane Del Pino]

Hello.

Default setting for ff3d to solve linear systems is Conjugate Gradient with and
epsilon of convergence set to 1E-5.

Here, to get a better solution, you need to get a smaller epsilon. Taking 1E-8
seems enough for your case:

Use the conjugate gradient option in your solver block:

...
solve(v) in O by M
cg(epsilon=1E-8)
{
...

Beware that your are solving your problem in a box which corresponds exactly
to your fictitious domain geometry. This can eventually lead to strange
behaviors. If your mesh fits the geometry I recommend not to use the fictitious
domain but directly the conformal FEM. It is very simple in your case, write
...
solve(v) in M
cg(epsilon=1E-8)
{
...

Best regards,
Stéphane.