Solver for Beltrami (linear) system

Functions/Subroutines
subroutine	mp00ac (Ndof, Xdof, Fdof, Ddof, Ldfjac, iflag)
	Solves Beltrami/vacuum (linear) system, given matrices.unpacking fluxes, helicity multiplier More...

Detailed Description

Function/Subroutine Documentation

mp00ac()

subroutine mp00ac	(	integer, intent(in)	Ndof,
		real, dimension(1:ndof), intent(in)	Xdof,
		real, dimension(1:ndof)	Fdof,
		real, dimension(1:ldfjac,1:ndof)	Ddof,
		integer, intent(in)	Ldfjac,
		integer	iflag
	)

Solves Beltrami/vacuum (linear) system, given matrices.unpacking fluxes, helicity multiplier

• The vector of "parameters", \(\boldsymbol{\mu}\), is unpacked. (Recall that \(\boldsymbol{\mu}\) was "packed" in ma02aa() .) In the following, \(\boldsymbol{\psi} \equiv (\Delta \psi_t, \Delta \psi_p)^T\).

construction of linear system

• The equation \(\nabla \times {\bf B} = \mu {\bf B}\) is cast as a matrix equation,

  \begin{eqnarray} {\cal M} \cdot {\bf a} = {\cal R},\end{eqnarray}

  where \({\bf a}\) represents the degrees-of-freedom in the magnetic vector potential, \({\bf a} \equiv \{ {\color{red} {A_{\theta,e,i,l}}}, {\color{blue}{A_{\zeta, e,i,l}}}, \dots \}\).
• The matrix \({\cal M}\) is constructed from \({\cal A}\equiv\,\)dMA and \({\cal D}\equiv\,\)dMD, which were constructed in matrix() , according to

  \begin{eqnarray} {\cal M} \equiv {\cal A} - \mu {\cal D}. \end{eqnarray}

  Note that in the vacuum region, \(\mu=0\), so \({\cal M}\) reduces to \({\cal M} \equiv {\cal A}\).
• The construction of the vector \({\cal R}\) is as follows:
  • if Lcoordinatesingularity=T , then

    \begin{eqnarray} {\cal R} \equiv - \left( {\cal B} - \mu {\cal E} \right) \cdot \boldsymbol{\psi} \end{eqnarray}

  • if Lcoordinatesingularity=F and Lplasmaregion=T , then

    \begin{eqnarray} {\cal R} \equiv - {\cal B} \cdot \boldsymbol{\psi} \end{eqnarray}

  • if Lcoordinatesingularity=F and Lvacuumregion=T , then

    \begin{eqnarray} {\cal R} \equiv - {\cal G} - {\cal B} \cdot \boldsymbol{\psi} \end{eqnarray}

  The quantities \({\cal B}\equiv\,\)dMB, \({\cal E}\equiv\,\)dME and \({\cal G}\equiv\,\)dMG are constructed in matrix() .

solving linear system

It is not assumed that the linear system is positive definite. The LAPACK routine DSYSVX is used to solve the linear system.

unpacking, ...

• The magnetic degrees-of-freedom are unpacked by packab() .
• The error flag, ImagneticOK , is set that indicates if the Beltrami fields were successfully constructed.

construction of "constraint" function

• The construction of the function \({\bf f}(\boldsymbol{\mu})\) is required so that iterative methods can be used to construct the Beltrami field consistent with the required constraints (e.g. on the enclosed fluxes, helicity, rotational-transform, ...).

  See also

  ma02aa() for additional details.

  plasma region

  • For Lcoordinatesingularity=T , the returned function is:

    \begin{eqnarray} {\bf f}(\mu,\Delta\psi_p) \equiv \left\{ \begin{array}{cccccccr} (& 0&,& 0&)^T, & \textrm{if }\texttt{Lconstraint} &=& -1 \\ (& 0&,& 0&)^T, & \textrm{if }\texttt{Lconstraint} &=& 0 \\ (&{{\,\iota\!\!\!}-}(+1)-\texttt{iota (lvol )}&,& 0&)^T, & \textrm{if }\texttt{Lconstraint} &=& 1 \\ (& ?&,& ?&)^T, & \textrm{if }\texttt{Lconstraint} &=& 2 \end{array}\right. \end{eqnarray}

  • For Lcoordinatesingularity=F , the returned function is:

    \begin{eqnarray} {\bf f}(\mu,\Delta\psi_p) \equiv \left\{ \begin{array}{cccccccr} (& 0&,& 0&)^T, & \textrm{if }\texttt{Lconstraint} &=& -1 \\ (& 0&,& 0&)^T, & \textrm{if }\texttt{Lconstraint} &=& 0 \\ (&{{\,\iota\!\!\!}-}(-1)-\texttt{oita(lvol-1)}&,&{{\,\iota\!\!\!}-}(+1)-\texttt{iota(lvol)}&)^T, & \textrm{if }\texttt{Lconstraint} &=& 1 \\ (& ?&,& ?&)^T, & \textrm{if }\texttt{Lconstraint} &=& 2 \end{array}\right. \end{eqnarray}

  vacuum region

  • For the vacuum region, the returned function is:

    \begin{eqnarray} {\bf f}(\Delta\psi_t,\Delta\psi_p) \equiv \left\{ \begin{array}{cccccccr} (& 0&,& 0&)^T, & \textrm{if }\texttt{Lconstraint} &=& -1 \\ (& I-\texttt{curtor} &,& G-\texttt{curpol} &)^T, & \textrm{if }\texttt{Lconstraint} &=& 0 \\ (&{{\,\iota\!\!\!}-}(-1)-\texttt{oita(lvol-1)}&,& G-\texttt{curpol} &)^T, & \textrm{if }\texttt{Lconstraint} &=& 1 \\ (& ?&,& ?&)^T, & \textrm{if }\texttt{Lconstraint} &=& 2 \end{array}\right. \end{eqnarray}

• The rotational-transform, \({{\,\iota\!\!\!}-}\), is computed by tr00ab() ; and the enclosed currents, \(I\) and \(G\), are computed by curent() .

early termination

• If \(|{\bf f}| < \) mupftol , then early termination is enforced (i.e., iflag is set to a negative integer). (See ma02aa() for details of how mp00ac() is called iteratively.)

Parameters

[in]	Ndof
[in]	Xdof
	Fdof
	Ddof
[in]	Ldfjac
	iflag	indicates whether (i) iflag=1: "function" values are required; or (ii) iflag=2: "derivative" values are required

References allglobal::adotx, allglobal::ate, allglobal::ato, allglobal::aze, allglobal::azo, allglobal::cpus, curent(), inputlist::curpol, inputlist::curtor, allglobal::ddotx, allglobal::diotadxup, allglobal::ditgpdxtp, allglobal::dma, allglobal::dmas, allglobal::dmb, allglobal::dmd, allglobal::dmds, allglobal::dmg, allglobal::dpflux, allglobal::dtflux, inputlist::epsgmres, inputlist::epsilu, allglobal::gmreslastsolution, constants::half, inputlist::helicity, allglobal::idmas, allglobal::im, allglobal::imagneticok, allglobal::in, intghs(), inputlist::iota, allglobal::iquad, allglobal::ivol, allglobal::jdmas, allglobal::labintegral, allglobal::lbbintegral, inputlist::lconstraint, allglobal::lcoordinatesingularity, inputlist::lgmresprec, allglobal::liluprecond, inputlist::lmatsolver, allglobal::lplasmaregion, inputlist::lrad, allglobal::lvacuumregion, numerical::machprec, allglobal::mn, allglobal::mns, mp00ac(), allglobal::mpi_comm_spec, mtrxhs(), inputlist::mu, inputlist::mupftol, allglobal::myid, allglobal::nadof, allglobal::ncpu, allglobal::ndmas, allglobal::ndmasmax, inputlist::nitergmres, allglobal::notmatrixfree, allglobal::notstellsym, allglobal::nt, inputlist::ntor, allglobal::nz, inputlist::oita, constants::one, fileunits::ounit, packab(), rungmres(), numerical::small, allglobal::solution, tr00ab(), inputlist::wmacros, allglobal::xoffset, allglobal::yesstellsym, and constants::zero.

Referenced by ma02aa(), and mp00ac().

Here is the call graph for this function:

Here is the caller graph for this function: