I need to solve for the toroidal flux of the magnetic field above an accretion disc.

For this purpose, I define the region over which the flux is non zero,

a= 2;
f[R_] = a (1 - R^2)^(1/2 + 1/100);
Ω = ImplicitRegion[z <= f[R], {{R, 0, 1}, {z, 0, f[0]}}];

Mathematica graphics

Then I solve for the force-free Grad-Shafranov equation:

eqn0 =R D[P[R, z], {R, 2}] + R D[P[R, z], {z, 2}] - D[P[R, z], R] == -R/2;

as follows

P0 = NDSolveValue[{eqn0, DirichletCondition[P[R, z] == 0, R == 0], 
    DirichletCondition[P[R, z] == 0, z == f[R]]},  P, {R, z} ∈ Ω];

I then plot the resulting (normalized) flux map:

np = NMaximize[P0[R, z], {R, z} ∈ Ω][[1]];
ContourPlot[ P0[R, z]/np, {R, z} ∈ Ω,PlotPoints -> 50, 
ImageSize -> Small, AspectRatio -> Sqrt[f[0]]]

Mathematica graphics

and look at the value of the flux on the boundary

Plot[P0[R, z]/np /. z -> f[R], {R, 0, 1}, PlotRange -> All]

Mathematica graphics

it seems to satisfy the boundary condition.

If I now look at the pressure above the cap:

grad2 = Grad[P0[R, z]/np, {R, z}] // (#.#/R^2 &);
Plot[grad2 /. z -> f[R], {R, 0, 1}]

Mathematica graphics

It is not smooth enough…

On top of that

If I decide to extend the height of the column over which the flux is defined, to say


Then the accuracy of he map deteriorates considerably for the map

Mathematica graphics

and even more for the pressure:

Mathematica graphics


How can improve the accuracy of the solution found by NDSolveValue?

I understand that there are options such as PrecisionGoal or Method, but I guess I am trying to ask a more general question:

More generally, what is the best learning strategy within mathematica to be able to find such improvement? (a.k.a how not to get lost in the documentation?). The idea being, next time I can figure this myself :-)

Mathematica gives the following warning, which is undoubtedly a hint

NDSolveValue::femcscd: The PDE is convection dominated and the result may not be stable. Adding artificial diffusion may help.

but I do not know how to follow it up.

PS: If, instead of

 f[R_] = a (1 - R^2)^(1/2 + 1/100);

I have

f[R_] = a (1 - R^2)^(1/2);

the integrator also fails miserably, which is rather odd.


This is not a full answer but what I found so far…

Following @user21's link

If I use

 Method -> {"PDEDiscretization" -> {"FiniteElement", 
"MeshOptions" -> {"MaxCellMeasure" -> 0.0001}, 
"IntegrationOrder" -> 3}}

as in

P0 = NDSolveValue[{eqn0, DirichletCondition[P[R, z] == 0, R == 0], 
    DirichletCondition[P[R, z] == 0, z == f[R]]}, 
   P, {R, z} ∈ Ω
   , Method -> {"PDEDiscretization" -> {"FiniteElement", 
       "MeshOptions" -> {"MaxCellMeasure" -> 0.0001}, 
       "IntegrationOrder" -> 3}}];

I get

Mathematica graphics

which is clearly better.

It still fails to get a very accurate pressure for a=15 though.

Mathematica graphics

One should not consider this an answer because it falls short in defining a strategy to identify what the proper options should be. Here I am just fishing randomly.

  • $\begingroup$ I can not look at this right away, but have a look at this desciption of convection dominated PDEs. $\endgroup$ – user21 Jul 20 '15 at 17:36
  • $\begingroup$ @user21 thanks. Seems quite relevant indeed. $\endgroup$ – chris Jul 20 '15 at 17:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.