I want to compute flux coordinates $\{\psi,\theta,\chi\}$ as functions of cylindrical coordinates $\{r,\theta,z\}$ in the problem of ballooning mode instability in mirror traps (also called open traps) for plasma confinement and heating. Magnetic field in a mirror trap is expressed through the flux potentials as follows. In the cylindrical coordinates
$\vec{B}(r,z)= \left\{\frac{1}{r}\frac{\partial \psi (r,z)}{\partial z},0,-\frac{1}{r}\frac{\partial \psi (r,z)}{\partial r}\right\}$
or
$\vec{B}(r,z)= \left\{G(r,z)\frac{\partial \chi (r,z)}{\partial r},0,G(r,z)\frac{\partial \chi (r,z)}{\partial z}\right\}$
Function $\psi(r,z)$ can be found as a solution of the Grad-Shafranov equation, however method of searching of functions $G$ and $\chi$ is not widely known. Moreover, for investigating plasma stability it is better to express $r$, $z$ and $G$ as functions of $\psi$ and $\chi$. Appropriate equations were derived in 2011 by Arsenin and Terekhin:
(Eq.9)
${\partial_\psi}\left(G(\psi ,\chi ) {\partial_\psi r^2(\psi ,\chi )}\right) +{\partial_\chi}\left(\frac{{\partial_\chi r^2(\psi,\chi)}}{G(\psi ,\chi ) r^2(\psi ,\chi )}\right)=0$,
(Eq.15)
$2 \pi \frac{\partial p(\psi )}{\partial \psi } \left(\left(G(\psi ,\chi ) {\partial_\psi {r^2}(\psi ,\chi )}\right)^2 +\frac{\left({\partial_\chi {r^2}(\psi ,\chi )}\right)^2}{{r^2}(\psi ,\chi )}\right)+\frac{\partial G(\psi ,\chi )^2}{\partial \psi }=0$
Here $r^2(\psi,\chi)$ is the square of $r(\psi,\chi)$ and $p(\psi)$ is radial profile of the plasma pressure. These equations should be solved on a rectangle area
$\Omega=\text{Rectangle}[\{0,0\},\{\psi_m,\chi_m\}]$
with the boundary conditions
$r^2(\psi,\chi)==0$,
$G(\psi,\chi)==1$,
$G(\psi,\chi)\ \partial_\psi r^2(\psi, \chi)=2/B_0(\chi)$
at the axis $\psi=0$ of the axial symmetry, where$B_0(\chi)$ is the magnetic field on the axis. Instead of the last (Neuman) BC we can try posing Dirichlet BC at the opposite edge of the rectangle:
$r^2(\psi_m, \chi)=2\psi_m/B_0(\chi)$
At the boundaries $\chi=0$ and $\chi=\chi_m$ we have the following conditions:
$\partial_\chi r^2(\psi,\chi)=0$.
Formal solution of Eq.15 can be found with DSolve. It expresses $G(\psi,\chi)$ through a double integral of $r^2(\psi,\chi)$. Calculating the double integral might be time-consuming so finding both $r2=r^2$ and $G$ with the use of NDSolve seems preferable. Unfortunately, straightforward attaсk on the problem failed with the error message:
"The PDE coefficient ... does not evaluate to a numeric scalar
at the coordinate {0.0025000000000000005`,2}; it evaluated to
Indeterminate instead"
EDIT #1
As it is said in FEMDocumentation/tutorial/FiniteElementProgramming, to solve nonlinear PDE it should be put in an inactive form. Here is my second attempt to solve the above formulated problem.
Here is inactive version of my code. Note that $G2$ and $r2$ are squares of $G$ and $r$ respectively.
Introduce some functions and parameters:
ClearAll[\[Psi], \[Chi], \[Theta], B, B0, G, r, z, r2]
(* set functions p[\[Psi]] and B0[\[Chi]]*)
funs = {B0 ->
Function[\[Chi],
B00 (1 + (K - 1) Sin[(\[Pi] \[Chi])/(2 \[Chi]m)]^2)]
, p -> Function[\[Psi], p0 (1 - \[Psi]/\[Psi]m)]}
(* set input parameters *)params = {\[Chi]m ->
B00 Sqrt[K] L, \[Psi]m -> B00 a^2/2, B00 -> 1, L -> 1, K -> 16,
a -> 0.1, p0 -> A/(4 \[Pi]), A -> 1}
(* define region where we look for a solution *)
\[CapitalOmega] = Rectangle[{0, 0}, {\[Psi]m, \[Chi]m}] //. params
Prepare functions r2[0] and G2[0] for setting Initial Seeding an BCsNote that we put into initial seeding and BCs a so called paraxial (i.e. near-axis solution of our equations):
(* Define functions for Initial Seeding and BCs *)
r2[0][\[Psi]_, \[Chi]_] = (2 \[Psi])/
B0[\[Chi]] - (
B0[\[Chi]] Derivative[1][B0][\[Chi]]^2 -
8 \[Pi] Derivative[1][p][0] -
B0[\[Chi]]^2 (B0^\[Prime]\[Prime])[\[Chi]])/(
2 B0[\[Chi]]^3) \[Psi]^2 /. funs //. params
G2[0][\[Psi]_, \[Chi]_] =
1 - (8 \[Pi] \[Psi] Derivative[1][p][0])/B0[\[Chi]]^2 /. funs //.
params
Set Eq. 9 in active and inactive forms and check that they are equivalent:
eq09 = \!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Psi]\)]\((
\*SqrtBox[\(G2[\[Psi], \[Chi]]\)]
\*SubscriptBox[\(\[PartialD]\), \(\[Psi]\)]r2[\[Psi], \[Chi]])\)\) + \
\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Chi]\)]\((
\*FractionBox[\(
\*SubscriptBox[\(\[PartialD]\), \(\[Chi]\)]r2[\[Psi], \[Chi]]\), \(
\*SqrtBox[\(G2[\[Psi], \[Chi]]\)] r2[\[Psi], \[Chi]]\)])\)\) - 0;
\[DoubleStruckCapitalC]09 = {{Sqrt[G2[\[Psi], \[Chi]]], 0}, {0, 1/(
Sqrt[G2[\[Psi], \[Chi]]] r2[\[Psi], \[Chi]])}};
\[DoubleStruckCapitalC]09 // MatrixForm
Inactive[Grad][r2[\[Psi], \[Chi]], {\[Psi], \[Chi]}];
\[DoubleStruckCapitalC]09 .
Inactivate[Grad[r2[\[Psi], \[Chi]], {\[Psi], \[Chi]}]];
Eq09 = Inactive[Div][\[DoubleStruckCapitalC]09 .
Inactive[Grad][
r2[\[Psi], \[Chi]], {\[Psi], \[Chi]}], {\[Psi], \[Chi]}]
Activate[Eq09]
% == eq09
Set Eq. 15 in active and inactive forms and check that they are equivalent:
(* Eq. 15 from the above cited paper *)
eq15 = \!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Psi]\)]\(G2[\[Psi], \[Chi]]\)\) +
2 \[Pi] (\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Psi]\)]\(p[\[Psi]]\)\)) \
(G2[\[Psi], \[Chi]] (\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Psi]\)]\(r2[\[Psi], \
\[Chi]]\)\))^2 + (\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Chi]\)]\(r2[\[Psi], \
\[Chi]]\)\))^2/r2[\[Psi], \[Chi]]) - 0
\[DoubleStruckCapitalC]15 = {{G2[\[Psi], \[Chi]], 0}, {0, 1/
r2[\[Psi], \[Chi]]}};
\[DoubleStruckCapitalC]15 // MatrixForm;
(* This version of Eq. 15 begin to work somehow *)
Eq15 = D[G2[\[Psi], \[Chi]], \[Psi]] + 2 \[Pi] \!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Psi]\)]\(p[\[Psi]]\)\) Inactive[
Grad][r2[\[Psi], \[Chi]], {\[Psi], \[Chi]}] . \
\[DoubleStruckCapitalC]15 .
Inactive[Grad][r2[\[Psi], \[Chi]], {\[Psi], \[Chi]}];
Eq15 = D[G2[\[Psi], \[Chi]], \[Psi]] + 2 \[Pi] \!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Psi]\)]\(p[\[Psi]]\)\) Inactive[
Grad][
r2[\[Psi], \[Chi]], {\[Psi], \[Chi]}] . \
\[DoubleStruckCapitalC]15 .
Inactive[Grad][r2[\[Psi], \[Chi]], {\[Psi], \[Chi]}];
Activate[Eq15]
% == eq15 // Simplify
Set boundary conditions:
BCD17 = DirichletCondition[r2[\[Psi], \[Chi]] == 0, \[Psi] == 0];
BCD20 = DirichletCondition[
r2[\[Psi], \[Chi]] == r2[0][\[Psi], \[Chi]], \[Psi] == 0];
BCD21 = DirichletCondition[
G2[\[Psi], \[Chi]] == G2[0][\[Psi], \[Chi]], \[Psi] == 0];
BCN22 = NeumannValue[(* G[\[Psi],\[Chi]]\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Psi]\)]\(r2[\[Psi], \[Chi]]\)\)\
\[Equal] *)2/B0[\[Chi]], \[Psi] == 0];
Combine a system of PDEs:
(* Сombine system of PDEs *)
(* Version with Neuman and Dirichlet BCs on lower edge of
\[CapitalOmega] *)
PDEs = {Eq09 == BCN22, Eq15 == 0, BCD17, BCD21} /. funs //. params;
(* Version with Dirichlet BCs on lower and uppers edges of
\[CapitalOmega] *)
PDEs = {Eq09 == 0, Eq15 == 0, BCD17, BCD20, BCD21} /. funs //. params
Solve it:
"Solution"
NDSolveValue[PDEs, {r2,
G2}, {\[Psi], \[Chi]} \[Element] \[CapitalOmega],
InitialSeeding -> {r2[\[Psi], \[Chi]] == r2[0][\[Psi], \[Chi]],
G2[\[Psi], \[Chi]] == G2[0][\[Psi], \[Chi]]}]
Unfortunately this story ends with the following error:
FindRoot::dfmin: The minimal damping factor of 1/10000 has been reached.
NDSolveValue::fempsf: PDESolve could not find a solution.
P.S.
In my other Question How to pose Dirichlet and Neumann BCs on same boundary? I noted that NDSolve cannot solve Laplace equation if Dirichlet and Neuman BCs are posed on same edge of the rectangle region $\Omega$. This is also the case with the above problem. Therefor instead of Neuman BC
$G(\psi,\chi)\partial_\psi r^2(\psi, \chi)=2/B_0(\chi)$
at $\psi=0$ it is better to pose Dirichlet BC
$r^2(\psi,\chi)=2\psi/B_0(\chi)$
at $\psi=\psi_m$.
Any help would be greatly appreciated.
NDSolveValue::femnlmdor::"The maximum derivative order of the nonlinear PDE coefficients for \ the Finite Element Method is larger than 1. It may help to rewrite \ the PDE in inactive form."$\endgroup$FindRoot::dfmin: The minimal damping factor of 1/10000 has been reached.I don't understand how to bypass it. $\endgroup$