Solver/Driver

Functions/Subroutines
subroutine	ma02aa (lvol, NN)
	Constructs Beltrami field in given volume consistent with flux, helicity, rotational-transform and/or parallel-current constraints. More...

Detailed Description

Function/Subroutine Documentation

ma02aa()

subroutine ma02aa	(	integer, intent(in)	lvol,
		integer, intent(in)	NN
	)

Constructs Beltrami field in given volume consistent with flux, helicity, rotational-transform and/or parallel-current constraints.

Parameters

[in]	lvol	index of nested volume for which to run this
[in]	NN	number of degrees of freedom in the (packed format) vector potential;

sequential quadratic programming

• Only relevant if LBsequad=T . See LBeltrami for details.
• Documentation on the implementation of E04UFF is under construction.

Newton method

• Only relevant if LBnewton=T . See LBeltrami for details.

linear method

• Only relevant if LBlinear=T . See LBeltrami for details.
• The quantity \(\mu\) is not not treated as a "magnetic" degree-of-freedom equivalent to in the degrees-of-freedom in the magnetic vector potential (as it strictly should be, because it is a Lagrange multiplier introduced to enforce the helicity constraint).
• In this case, the Beltrami equation, \(\nabla \times {\bf B} = \mu {\bf B}\), is linear in the magnetic degrees-of-freedom.

• The algorithm proceeds as follows:

  plasma volumes

  • In addition to the enclosed toroidal flux, \(\Delta \psi_t\), which is held constant in the plasma volumes, the Beltrami field in a given volume is assumed to be parameterized by \(\mu\) and \(\Delta \psi_p\). (Note that \(\Delta \psi_p\) is not defined in a torus.)
  • These are "packed" into an array, e.g. \(\boldsymbol{\mu} \equiv (\mu, \Delta\psi_p)^T\), so that standard library routines , e.g. C05PCF, can be used to (iteratively) find the appropriately-constrained Beltrami solution, i.e. \({\bf f}(\boldsymbol{\mu})=0\).

  • The function \({\bf f}(\boldsymbol{\mu})\), which is computed by mp00ac(), is defined by the input parameter Lconstraint:

    • If Lconstraint = -1, 0, then \(\boldsymbol{\mu}\) is not varied and Nxdof=0.
    • If Lconstraint = 1, then \(\boldsymbol{\mu}\) is varied to satisfy the transform constraints; and Nxdof=1 in the simple torus and Nxdof=2 in the annular regions. (Note that in the "simple-torus" region, the enclosed poloidal flux \(\Delta\psi_p\) is not well-defined, and only \(\mu=\boldsymbol{\mu}_1\) is varied in order to satisfy the transform constraint on the "outer" interface of that volume.)

    • Todo:

      If Lconstraint = 2, then \(\mu=\boldsymbol{\mu}_1\) is varied in order to satisfy the helicity constraint, and \(\Delta\psi_p=\boldsymbol{\mu}_2\) is not varied, and Nxdof=1. (under re-construction)

  vacuum volume

  • In the vacuum, \(\mu=0\), and the enclosed fluxes, \(\Delta \psi_t\) and \(\Delta \psi_p\), are considered to parameterize the family of solutions. (These quantities may not be well-defined if \({\bf B}\cdot{\bf n}\ne 0\) on the computational boundary.)
  • These are "packed" into an array, \(\boldsymbol{\mu} \equiv (\Delta\psi_t,\Delta\psi_p)^T\), so that, as above, standard routines can be used to iteratively find the appropriately constrained solution, i.e. \({\bf f}(\boldsymbol{\mu})=0\).

  • The function \({\bf f}(\boldsymbol{\mu})\), which is computed by mp00ac(), is defined by the input parameter Lconstraint:

    • If Lconstraint = -1, then \(\boldsymbol{\mu}\) is not varied and Nxdof=0.
    • If Lconstraint = 0,2, then \(\boldsymbol{\mu}\) is varied to satisfy the enclosed current constraints, and Nxdof=2.
    • If Lconstraint = 1, then \(\boldsymbol{\mu}\) is varied to satisfy the constraint on the transform on the inner boundary \(\equiv\) plasma boundary and the "linking" current, and Nxdof=2.

• The Beltrami fields, and the rotational-transform and helicity etc. as required to determine the function \({\bf f}(\boldsymbol{\mu})\) are calculated in mp00ac().
• This routine, mp00ac(), is called iteratively if Nxdof>1 via C05PCF to determine the appropriately constrained Beltrami field, \({\bf B}_{\boldsymbol{\mu}}\), so that \({\bf f}(\boldsymbol{\mu})=0\).
• The input variables mupftol and mupfits control the required accuracy and maximum number of iterations.

• If Nxdof=1, then mp00ac() is called only once to provide the Beltrami fields with the given value of \(\boldsymbol{\mu}\).

debugging: finite-difference confirmation of the derivatives of the rotational-transform

• Note that the rotational-transform (if required) is calculated by tr00ab(), which is called by mp00ac().
• If Lconstraint=1, then mp00ac() will ask tr00ab() to compute the derivatives of the transform with respect to variations in the helicity-multiplier, \(\mu\), and the enclosed poloidal-flux, \(\Delta\psi_p\), so that C05PCF may more efficiently find the solution.
• The required derivatives are

  \begin{eqnarray} \frac{\partial {{\,\iota\!\!\!}-}}{\partial \mu}\\ \frac{\partial {{\,\iota\!\!\!}-}}{\partial \Delta \psi_p} \end{eqnarray}

  to improve the efficiency of the iterative search. A finite difference estimate of these derivatives is available; need DEBUG, Lcheck=2 and Lconstraint=1.

References allglobal::ate, allglobal::cpus, df00ab(), allglobal::dma, allglobal::dmb, allglobal::dmd, allglobal::dpflux, allglobal::dtflux, constants::half, inputlist::helicity, allglobal::im, allglobal::imagneticok, allglobal::in, allglobal::ivol, allglobal::labintegral, allglobal::lbbintegral, allglobal::lblinear, allglobal::lbnewton, allglobal::lbsequad, inputlist::lcheck, inputlist::lconstraint, allglobal::lcoordinatesingularity, allglobal::localconstraint, allglobal::lplasmaregion, inputlist::lrad, allglobal::lvacuumregion, ma02aa(), allglobal::mbpsi, allglobal::mn, mp00ac(), allglobal::mpi_comm_spec, inputlist::mu, inputlist::mupfits, inputlist::mupftol, allglobal::myid, allglobal::ncpu, allglobal::nfielddof, constants::one, fileunits::ounit, packab(), numerical::small, allglobal::solution, constants::ten, numerical::vsmall, inputlist::wmacros, allglobal::xoffset, and constants::zero.

Referenced by dfp100(), and ma02aa().

Here is the call graph for this function:

Here is the caller graph for this function: