SPEC 3.20
Stepped Pressure Equilibrium Code
Functions/Subroutines
Free-Boundary Computation

Functions/Subroutines

subroutine bnorml (mn, Ntz, efmn, ofmn)
 Computes \({\bf B}_{Plasma} \cdot {\bf e}_\theta \times {\bf e}_\zeta \;\) on the computational boundary, \(\partial {\cal D}\). More...
 
subroutine casing (teta, zeta, gBn, icasing)
 Constructs the field created by the plasma currents, at an arbitrary, external location using virtual casing. More...
 
subroutine casinggrid (xyz, ntz, Pbxyz, Jxyz, vcstride, gBn)
 Compute the field produced by the plasma currents, at an arbitrary, external location using virtual casing. More...
 
subroutine surfacecurrent (teta, zeta, pxyz, jj)
 Compute the position and surface current on the plasma boundary for the virtual casingCalculation of integrand More...
 
subroutine dvcfield (Ndim, tz, Nfun, vcintegrand)
 Differential virtual casing integrand. More...
 

Detailed Description

Function/Subroutine Documentation

◆ bnorml()

subroutine bnorml ( integer, intent(in)  mn,
integer, intent(in)  Ntz,
real, dimension(1:mn), intent(out)  efmn,
real, dimension(1:mn), intent(out)  ofmn 
)

Computes \({\bf B}_{Plasma} \cdot {\bf e}_\theta \times {\bf e}_\zeta \;\) on the computational boundary, \(\partial {\cal D}\).

free-boundary constraint

  • The normal field at the computational boundary, \(\partial {\cal D}\), should be equal to \(\left({\bf B}_P + {\bf B}_C\right)\cdot {\bf e}_\theta \times {\bf e}_\zeta\), where \({\bf B}_P\) is the "plasma" field (produced by internal plasma currents) and is computed using virtual casing, and \({\bf B}_C\) is the "vacuum" field (produced by the external coils) and is given on input.
  • The plasma field, \({\bf B}_P\), can only be computed after the equilibrium is determined, but this information is required to compute the equilibrium to begin with; and so there is an iteration involved [7] .
  • Suggested values of the vacuum field can be self generated; see xspech() for more documentation on this.

compute the normal field on a regular grid on the computational boundary

  • For each point on the compuational boundary, casing() is called to compute the normal field produced by the plasma currents.
  • An FFT gives the required Fourier harmonics.
See also
casing.f90
Parameters
[in]mntotal number of Fourier harmonics
[in]Ntztotal number of grid points in \(\theta\) and \(zeta\)
[out]efmneven Fourier coefficients
[out]ofmnodd Fouier coefficients

The surface currents get precomputed at all points, but maybe we don't have to integrate over the whole grid. Start with a large stride over the plasmaboundary (= how many points to skip at each step), then get progressively finer. Do at least vcasingits iterations (sqrt(vcasingits) resolution per dimension) and halve the stride until the accuracy is reached.

At least one step in the resolution cascade is required to determine an estimate of the current accuracy.

References allglobal::ate, allglobal::ato, allglobal::aze, allglobal::azo, bnorml(), casing(), casinggrid(), allglobal::cfmn, allglobal::cpus, allglobal::dxyz, allglobal::globaljk, allglobal::gteta, allglobal::guvij, allglobal::gzeta, constants::half, inputlist::igeometry, allglobal::ijimag, allglobal::ijreal, allglobal::im, allglobal::in, allglobal::jiimag, allglobal::jireal, allglobal::jxyz, inputlist::lcheck, allglobal::lcoordinatesingularity, inputlist::lrad, fileunits::lunit, inputlist::lvcgrid, allglobal::mpi_comm_spec, allglobal::mvol, allglobal::myid, allglobal::ncpu, inputlist::nfp, allglobal::notstellsym, allglobal::nt, allglobal::nxyz, allglobal::nz, constants::one, fileunits::ounit, allglobal::pbxyz, constants::pi, constants::pi2, allglobal::pi2nfp, allglobal::prevcstride, allglobal::rij, allglobal::sfmn, allglobal::sg, numerical::small, surfacecurrent(), constants::ten, allglobal::tt, constants::two, inputlist::vcasingits, inputlist::vcasingper, inputlist::vcasingtol, inputlist::vcnt, inputlist::vcnz, allglobal::virtualcasingfactor, numerical::vsmall, inputlist::wmacros, allglobal::yesstellsym, constants::zero, and allglobal::zij.

Referenced by bnorml(), preset(), and spec().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ casing()

subroutine casing ( real, intent(in)  teta,
real, intent(in)  zeta,
real, intent(out)  gBn,
integer, intent(out)  icasing 
)

Constructs the field created by the plasma currents, at an arbitrary, external location using virtual casing.

Compute the external magnetic field using virtual casing.

Theory and numerics

  • Required inputs to this subroutine are the geometry of the plasma boundary,

    \begin{eqnarray} {\bf x}(\theta,\zeta) \equiv x(\theta,\zeta) {\bf i} + y(\theta,\zeta) {\bf j} + z(\theta,\zeta) {\bf k}, \end{eqnarray}

    and the tangential field on this boundary,

    \begin{eqnarray} {\bf B}_s=B^\theta {\bf e}_\theta + B^\zeta {\bf e}_\zeta, \end{eqnarray}

    where \(\theta\) and \(\zeta\) are arbitrary poloidal and toroidal angles, and \({\bf e}_\theta \equiv \partial {\bf x}/\partial \theta\), \({\bf e}_\zeta \equiv \partial {\bf x}/\partial \zeta\). This routine assumes that the plasma boundary is a flux surface, i.e. \({\bf B} \cdot {\bf e}_\theta \times {\bf e}_\zeta = 0\).
  • The virtual casing principle (Shafranov & Zakharov (1972) [19], Lazerson (2012) [9] and Hanson (2015) [3]) shows that the field outside/inside the plasma arising from plasma currents inside/outside the boundary is equivalent to the field generated by a surface current,

    \begin{eqnarray} {\bf j} = {\bf B}_s \times {\bf n}, \end{eqnarray}

    where \({\bf n}\) is normal to the surface.
  • The field at some arbitrary point, \(\bar {\bf x}\), created by this surface current is given by

    \begin{eqnarray} {\bf B}({\bf \bar x}) = -\frac{1}{4\pi} \int_{\cal S} \frac{\left( {\bf B}_s \times d{\bf s} \right) \times {\bf \hat r}}{r^2}, \end{eqnarray}

    where \(d{\bf s} \equiv {\bf e}_\theta \times {\bf e}_\zeta \; d\theta d\zeta\).
  • For ease of notation introduce

    \begin{eqnarray} {\bf J} & \equiv & {\bf B}_s \times d{\bf s} = \alpha \; {\bf e}_\theta - \beta \; {\bf e}_\zeta, \end{eqnarray}

    where \(\alpha \equiv B_\zeta = B^\theta g_{\theta\zeta} + B^\zeta g_{\zeta\zeta}\) and \(\beta \equiv B_\theta = B^\theta g_{\theta\theta} + B^\zeta g_{\theta\zeta}\).
  • We may write in Cartesian coordinates \({\bf J} = j_x \; {\bf i} + j_y \; {\bf j} + j_z \; {\bf k}\), where

    \begin{eqnarray} j_x & = & \alpha \; x_\theta - \beta \; x_\zeta \\ j_y & = & \alpha \; y_\theta - \beta \; y_\zeta \\ j_z & = & \alpha \; z_\theta - \beta \; z_\zeta . \end{eqnarray}

  • Requiring that the current,

    \begin{eqnarray} {\bf j} & \equiv & \nabla \times {\bf B} = {\sqrt g}^{-1}(\partial_\theta B_\zeta -\partial_\zeta B_\theta) \; {\bf e}_s + {\sqrt g}^{-1}(\partial_\zeta B_s -\partial_s B_\zeta ) \; {\bf e}_\theta + {\sqrt g}^{-1}(\partial_s B_\theta-\partial_\theta B_s ) \; {\bf e}_\zeta, \end{eqnarray}

    has no normal component to the surface, i.e. \({\bf j}\cdot \nabla s=0\), we obtain the condition \(\partial_\theta B_\zeta = \partial_\zeta B_\theta\), or \(\partial_\theta \alpha = \partial_\zeta \beta\). In axisymmetric configurations, where \(\partial_\zeta \beta=0\), we must have \(\partial_\theta \alpha=0\).
  • The displacement from an arbitrary point, \((X,Y,Z)\), to a point, \((x,y,z)\), that lies on the surface is given

    \begin{eqnarray} {\bf r} \equiv r_x \; {\bf i} + r_y \; {\bf j} + r_z \; {\bf k} = (X-x) \; {\bf i} + (Y-y) \; {\bf j} + (Z-z) \; {\bf k}. \end{eqnarray}

  • The components of the magnetic field produced by the surface current are then

    \begin{eqnarray} B^x &=& \oint\!\!\!\oint \!d\theta d\zeta \,\,\, (j_y r_z - j_z r_y)/r^3,\\ B^y &=& \oint\!\!\!\oint \!d\theta d\zeta \,\,\, (j_z r_x - j_x r_z)/r^3,\\ B^z &=& \oint\!\!\!\oint \!d\theta d\zeta \,\,\, (j_x r_y - j_y r_x)/r^3 \end{eqnarray}

    up to a scaling factor virtualcasingfactor \(=-1/4\pi\) that is taken into account at the end.
  • When all is said and done, this routine calculates

    \begin{eqnarray} \int_0^{2\pi} \int_0^{2\pi} \verb+vcintegrand+ \;\; d\theta d\zeta \end{eqnarray}

    for a given \((X,Y,Z)\), where vcintegrand is given in Eqn. \((\ref{eq:integrand_casing})\).
  • The surface integral is performed using DCUHRE , which uses an adaptive subdivision strategy and also computes absolute error estimates. The absolute and relative accuracy required are provided by the inputvar vcasingtol . The minimum number of function evaluations is provided by the inputvar vcasingits.
Parameters
[in]teta\(\theta\)
[in]zeta\(\zeta\)
[out]gBn\( \sqrt g {\bf B} \cdot {\bf n}\)
[out]icasingreturn flag from dcuhre()

References casing(), allglobal::cpus, dvcfield(), allglobal::dxyz, allglobal::globaljk, allglobal::mpi_comm_spec, allglobal::myid, allglobal::ncpu, allglobal::nxyz, fileunits::ounit, constants::pi, constants::pi2, inputlist::vcasingits, inputlist::vcasingper, inputlist::vcasingtol, fileunits::vunit, inputlist::wmacros, and constants::zero.

Referenced by bnorml(), casing(), dvcfield(), preset(), and surfacecurrent().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ casinggrid()

subroutine casinggrid ( real, dimension(1:3), intent(in)  xyz,
real, dimension(1:3), intent(in)  ntz,
real, dimension(1:vcnz*nfp*vcnt, 1:3), intent(in)  Pbxyz,
real, dimension(1:vcnz*nfp*vcnt, 1:3), intent(in)  Jxyz,
integer, intent(in)  vcstride,
real, intent(out)  gBn 
)

Compute the field produced by the plasma currents, at an arbitrary, external location using virtual casing.

casinggrid() computes the virtual casing field using a fixed resolution grid on the plasma boundary, in contrast to casing() which uses an adaptive integration routine. Because the evaluation positions are known in advance, the most expensive terms of the computation (surfacecurrents, or in the casing() routine dvcfield) are precomputed. The integral computed by casinggrid() is

\begin{eqnarray} \int_0^{2\pi} \int_0^{2\pi} \frac{{\bf j}(\theta, \zeta) \times {\bf D}(\theta, \zeta)}{\|{\bf D}(\theta', \zeta')\|^3}+ \;\; d\theta d\zeta \end{eqnarray}


with \({\bf r}\) approximately the distance vector from an evaluation point of the current surface \({\bf x}(\theta, \zeta)\) to the external point \((x,y,z)\), and the exact formula includes a regularization factor \(\epsilon\) Eqn. \((\ref{eq:vcasing_distance})\).

Numerically it uses a Kahan Summation algorithm to accumulate the contributions of all inner points, because truncation error quickly dominates the computation when thousands of points are used in both dimensions.

Parameters
[in]xyz\((x,y,z)\) cartesian coordinates on the computational boundary
[in]ntzcartesian normal vector on the computational boundary
[in]Pbxyzarray of cartesian coordinates on the plasma boundary
[in]Jxyzarray of surface currents \({\bf B}_{Plasma} \cdot {\bf e}_\theta \times {\bf e}_\zeta \;\) on the plasma boundary
[in]vcstrideinteger stride for the fixed resolution grid. vcstride = 1 computes the integral at full resolution, vcstride = 2 computes the integral at half resolution, etc.
[out]gBnnormal field \({\bf B}_{Plasma} \cdot n \;\) on the computational boundary

References allglobal::cpus, allglobal::dxyz, allglobal::globaljk, allglobal::mpi_comm_spec, allglobal::myid, allglobal::ncpu, inputlist::nfp, allglobal::nxyz, constants::one, fileunits::ounit, constants::pi2, allglobal::pi2nfp, constants::three, inputlist::vcasingeps, inputlist::vcnt, inputlist::vcnz, fileunits::vunit, and constants::zero.

Referenced by bnorml().

Here is the caller graph for this function:

◆ surfacecurrent()

subroutine surfacecurrent ( real, intent(in)  teta,
real, intent(in)  zeta,
real, dimension(1:3), intent(out)  pxyz,
real, dimension(1:3), intent(out)  jj 
)

Compute the position and surface current on the plasma boundary for the virtual casingCalculation of integrand

  • An adaptive integration is used to compute the integrals. Consequently, the magnetic field tangential to the plasma boundary is required at an arbitrary point. This is computed, as always, from \({\bf B} = \nabla \times {\bf A}\), and this provides \({\bf B} = B^\theta {\bf e}_\theta + B^\zeta {\bf e}_\zeta\). Recall that \(B^s=0\) by construction on the plasma boundary.

  • Then, the metric elements \(g_{\theta\theta}\), \(g_{\theta\zeta}\) and \(g_{\zeta\zeta}\) are computed. These are used to "lower" the components of the magnetic field, \({\bf B} = B_\theta \nabla \theta + B_\zeta \nabla \zeta\).

  • The distance between the "evaluate" point, \((X,Y,Z)\), and the given point on the surface, \((x,y,z)\) is computed.
  • If the computational boundary becomes too close to the plasma boundary, the distance is small and this causes problems for the numerics. I have tried to regularize this problem by introducing \(\epsilon \equiv \)inputvar vcasingeps. Let the "distance" be

    \begin{eqnarray} D \equiv \sqrt{(X-x)^2 + (Y-y)^2 + (Z-Z)^2} + \epsilon^2. \label{eq:vcasing_distance} \end{eqnarray}

  • On taking the limit that \(\epsilon \rightarrow 0\), the virtual casing integrand is

    \begin{eqnarray} \verb+vcintegrand+ \equiv ( B_x n_x + B_y n_y + B_z n_z ) ( 1 + 3 \epsilon^2 / D^2 ) / D^3, \label{eq:integrand_casing} \end{eqnarray}

    where the normal vector is \({\bf n} \equiv n_x {\bf i} + n_y {\bf j} + n_z {\bf k}\). The normal vector, Nxyz , to the computational boundary (which does not change) is computed in preset().

    Todo:
    This needs to be revised.

Parameters
[in]teta\(\theta\)
[in]zeta\(\zeta\)
[out]pxyz\((x,y,z)\) cartesian coordinates on the plasma boundary
[out]jjSurface current \({\bf B}_{Plasma} \cdot {\bf e}_\theta \times {\bf e}_\zeta \;\) on the plasma boundary

References allglobal::ate, allglobal::ato, allglobal::aze, allglobal::azo, casing(), allglobal::cpus, allglobal::first_free_bound, constants::half, inputlist::igeometry, allglobal::im, allglobal::in, allglobal::irbc, allglobal::irbs, allglobal::izbc, allglobal::izbs, inputlist::lrad, allglobal::mn, allglobal::mpi_comm_spec, allglobal::mvol, allglobal::myid, allglobal::ncpu, inputlist::nvol, constants::one, allglobal::tt, allglobal::yesstellsym, and constants::zero.

Referenced by bnorml(), and dvcfield().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ dvcfield()

subroutine dvcfield ( integer, intent(in)  Ndim,
real, dimension(1:ndim), intent(in)  tz,
integer, intent(in)  Nfun,
real, dimension(1:nfun), intent(out)  vcintegrand 
)

Differential virtual casing integrand.

Differential virtual casing integrand with the calling convention required by NAG for the adaptive integration routine.
This function wraps surfacecurrent and uses a global variable globaljk to determine the position on the computational boundary (not thread safe!).

Parameters
[in]Ndimnumber of parameters (==2)
[in]tz\(\theta\) and \(\zeta\)
[in]Nfunnumber of function values (==3)
[out]vcintegrandcartesian components of magnetic field

References casing(), allglobal::cpus, allglobal::dxyz, allglobal::first_free_bound, allglobal::globaljk, constants::half, allglobal::mpi_comm_spec, allglobal::myid, allglobal::ncpu, allglobal::nxyz, constants::one, surfacecurrent(), constants::three, inputlist::vcasingeps, and constants::zero.

Referenced by casing().

Here is the call graph for this function:
Here is the caller graph for this function: