xspech.f90 File Reference

Main program. More...

Functions/Subroutines
program	spec_main
	Main program of SPEC. More...
subroutine	xspech
	Main subroutine of SPEC. More...
subroutine	read_command_args
	Read command-line arguments; in particular, determine input file (name or extension). More...
subroutine	spec
	This is the main "driver" for the physics part of SPEC. More...
subroutine	final_diagnostics
	Final diagnostics. More...
subroutine	ending
	Closes output files, writes screen summary.

Detailed Description

Main program.

Function/Subroutine Documentation

spec_main()

program spec_main

Main program of SPEC.

This only calls the xpech() subroutine to do a stand-alone SPEC run.

Returns

none

References xspech().

Here is the call graph for this function:

xspech()

subroutine xspech

Main subroutine of SPEC.

This orchestrates a stand-alone SPEC run:

• read the input file
• solve the MRxMHD equilibrium (see spec() )
• run some diagnostics on the results
• write the output file(s)

reading input, allocating global variables

• The input namelists and geometry are read in via a call to readin() . A full description of the required input is given in global.f90 .
• Most internal variables, global memory etc., are allocated in preset() .
• All quantities in the input file are mirrored into the output file's group /input .

preparing output file group iterations

• The group /iterations is created in the output file. This group contains the interface geometry at each iteration, which is useful for constructing movies illustrating the convergence. The data structure in use is an unlimited array of the following compound datatype:

  DATATYPE  H5T_COMPOUND {
        H5T_NATIVE_INTEGER "nDcalls";
        H5T_NATIVE_DOUBLE "Energy";
        H5T_NATIVE_DOUBLE "ForceErr";
        H5T_ARRAY { [Mvol+1][mn] H5T_NATIVE_DOUBLE } "iRbc";
        H5T_ARRAY { [Mvol+1][mn] H5T_NATIVE_DOUBLE } "iZbs";
        H5T_ARRAY { [Mvol+1][mn] H5T_NATIVE_DOUBLE } "iRbs";
        H5T_ARRAY { [Mvol+1][mn] H5T_NATIVE_DOUBLE } "iZbc";
  }

restart files

• wrtend() is called to write the restart files.

References allglobal::broadcast_inputs(), allglobal::check_inputs(), allglobal::cpus, ending(), final_diagnostics(), sphdf5::finish_outfile(), sphdf5::hdfint(), sphdf5::init_convergence_output(), sphdf5::init_outfile(), numerical::machprec, sphdf5::mirror_input_to_outfile(), allglobal::mpi_comm_spec, allglobal::myid, allglobal::ncpu, fileunits::ounit, preset(), read_command_args(), numerical::small, spec(), numerical::vsmall, sphdf5::write_grid(), allglobal::wrtend(), and xspech().

Referenced by final_diagnostics(), spec(), spec_main(), and xspech().

Here is the call graph for this function:

Here is the caller graph for this function:

read_command_args()

subroutine read_command_args

Read command-line arguments; in particular, determine input file (name or extension).

• The input file name, ext , is given as the first command line input, and the input file itself is then ext.sp .
• Alternatively, you can directly specify the input file itself as ext.sp .
• You can also generate a template input file using xspec -i .
• Or print help information using xspec -h .
• Additional command line inputs recognized are:
  • -readin will immediately set Wreadin=T ; this may be over-ruled when the namelist screenlist is read

References allglobal::cpus, allglobal::mpi_comm_spec, allglobal::myid, fileunits::ounit, and inputlist::wreadin.

Referenced by xspech().

Here is the caller graph for this function:

spec()

subroutine spec

This is the main "driver" for the physics part of SPEC.

Picard iterations are performed (if in free-boundary mode) and within each Picard iteration, the fixed-boundary problem is solved (also iteratively).

packing geometrical degrees-of-freedom into vector

• If NGdof.gt.0 , where NGdof counts the geometrical degrees-of-freedom, i.e. the \(R_{bc}\), \(Z_{bs}\), etc., then packxi() is called to "pack" the geometrical degrees-of-freedom into position(0:NGdof) .

initialize adiabatic constants

• If Ladiabatic.eq.0 , then the "adiabatic constants" in each region, \(P_v\), are calculated as

  \begin{eqnarray} P_v \equiv p_v V_v^\gamma, \label{eq:adiabatic_xspech} \end{eqnarray}

  where \(p_v\equiv\,\)pressure(vvol) , the volume \(V_v\) of each region is computed by volume() , and the adiabatic index \(\gamma\equiv\,\)gamma .

solving force-balance

• If there are geometrical degress of freedom, i.e. if NGdof.gt.0 , then

  • Todo:

    If Lminimize.eq.1 , call pc00aa() to find minimum of energy functional using quasi-Newton, preconditioned conjugate gradient method, E04DGF

  • If Lfindzero.gt.0 , call newton() to find extremum of constrained energy functional using a Newton method, C05PDF .

post diagnostics

• The pressure is computed from the adiabatic constants from Eqn. \((\ref{eq:adiabatic_xspech})\), i.e. \(p=P/V^\gamma\).
• The Beltrami/vacuum fields in each region are re-calculated using dforce() .
• If Lcheck.eq.5 .or. LHevalues .or. LHevectors .or. Lperturbed.eq.1 , then the force-gradient matrix is examined using hesian() .

free-boundary: re-computing normal field

• If Lfreebound.eq.1 and Lfindzero.gt.0 and mfreeits.ne.0 , then the magnetic field at the computational boundary produced by the plasma currents is computed using bnorml() .
• The "new" magnetic field at the computational boundary produced by the plasma currents is updated using a Picard scheme:

  \begin{eqnarray} \texttt{Bns}_i^{j} = \lambda \, \texttt{Bns}_i^{j-1} + (1-\lambda) \texttt{Bns}_i, \label{eq:blending_xspech} \end{eqnarray}

  where \(j\) labels free-boundary iterations, the "blending parameter" is \(\lambda\equiv\,\)gBnbld , and Bns \(_i\) is computed by virtual casing. The subscript "$i$" labels Fourier harmonics.
• If the new (unblended) normal field is not sufficiently close to the old normal field, as quantified by gBntol , then the free-boundary iterations continue. This is quantified by

  \begin{eqnarray} \sum_i | \texttt{Bns}_i^{j-1} - \texttt{Bns}_i | / N, \label{eq:gBntol_xspech} \end{eqnarray}

  where \(N\) is the total number of Fourier harmonics.
• There are several choices that are available:
  • if mfreeits=-2 : the vacuum magnetic field (really, the normal component of the field produced by the external currents at the computational boundary) required to hold the given equlibrium is written to file. This information is required as input by FOCUS [8] for example. (This option probably needs to revised.)
  • if mfreeits=-1 : after the plasma field is computed by virtual casing, the vacuum magnetic field is set to exactly balance the plasma field (again, we are really talking about the normal component at the computational boundary.) This will ensure that the computational boundary itself if a flux surface of the total magnetic field.
  • if mfreeits=0 : the plasma field at the computational boundary is not updated; no "free-boundary" iterations take place.
  • if mfreeits>0 : the plasma field at the computational boundary is updated according to the above blending Eqn. \((\ref{eq:blending_xspech})\), and the free-boundary iterations will continue until either the tolerance condition is met (see gBntol and Eqn. \((\ref{eq:gBntol_xspech})\)) or the maximum number of free-boundary iterations, namely mfreeits , is reached. For this case, Lzerovac is relevant: if Lzerovac=1 , then the vacuum field is set equal to the normal field at every iteration, which results in the computational boundary being a flux surface. (I am not sure if this is identical to setting mfreeits=-1 ; the logic etc. needs to be revised.)

output files: vector potential

• The vector potential is written to file using ra00aa() .

References inputlist::adiabatic, allglobal::ate, allglobal::ato, allglobal::aze, allglobal::azo, allglobal::bbe, allglobal::bbo, allglobal::beltramierror, bnorml(), allglobal::bnserr, allglobal::cfmn, allglobal::cpus, dforce(), allglobal::dma, allglobal::dmb, allglobal::dmd, allglobal::dmg, allglobal::dpflux, allglobal::dtflux, allglobal::efmn, allglobal::first_free_bound, allglobal::forceerr, inputlist::gamma, inputlist::gbnbld, inputlist::gbntol, inputlist::helicity, hesian(), allglobal::ibnc, allglobal::ibns, inputlist::igeometry, allglobal::iie, allglobal::iio, allglobal::im, allglobal::imagneticok, allglobal::in, allglobal::irbc, allglobal::irbs, inputlist::isurf, allglobal::ivnc, allglobal::ivns, inputlist::ivolume, allglobal::izbc, allglobal::izbs, inputlist::ladiabatic, inputlist::lautoinitbn, inputlist::lcheck, inputlist::lconstraint, allglobal::lcoordinatesingularity, inputlist::lfindzero, inputlist::lfreebound, allglobal::lgdof, inputlist::lhevalues, inputlist::lhevectors, inputlist::lhmatrix, numerical::logtolerance, inputlist::lperturbed, allglobal::lplasmaregion, inputlist::lrad, fileunits::lunit, allglobal::lvacuumregion, inputlist::lzerovac, allglobal::mbpsi, inputlist::mfreeits, allglobal::mn, allglobal::mpi_comm_spec, inputlist::mu, constants::mu0, allglobal::myid, allglobal::ncpu, newton(), inputlist::nfp, allglobal::nfreeboundaryiterations, allglobal::ngdof, allglobal::notstellsym, inputlist::nppts, inputlist::nptrj, allglobal::ntz, inputlist::nvol, inputlist::odetol, allglobal::ofmn, constants::one, fileunits::ounit, packxi(), inputlist::pflux, inputlist::phiedge, constants::pi2, inputlist::pressure, inputlist::pscale, ra00aa(), inputlist::rbc, inputlist::rbs, allglobal::sfmn, allglobal::solution, inputlist::tflux, inputlist::vcasingtol, constants::version, volume(), numerical::vsmall, allglobal::vvolume, inputlist::wmacros, allglobal::wrtend(), xspech(), allglobal::yesstellsym, inputlist::zbc, inputlist::zbs, and constants::zero.

Referenced by dforce(), sphdf5::init_outfile(), and xspech().

Here is the call graph for this function:

Here is the caller graph for this function:

final_diagnostics()

subroutine final_diagnostics

Final diagnostics.

• sc00aa() is called to compute the covariant components of the magnetic field at the interfaces; these are related to the singular currents
• if Lcheck=1 , jo00aa() is called to compute the error in the Beltrami equation
• pp00aa() is called to construct the Poincare plot by field-line following.

References allglobal::beltramierror, brcast(), allglobal::btemn, allglobal::btomn, allglobal::bzemn, allglobal::bzomn, allglobal::cfmn, allglobal::cpus, allglobal::diotadxup, allglobal::dtflux, allglobal::efmn, inputlist::igeometry, allglobal::imagneticok, allglobal::iquad, allglobal::ismyvolume(), allglobal::ismyvolumevalue, inputlist::isurf, inputlist::ivolume, jo00aa(), inputlist::lcheck, allglobal::lcoordinatesingularity, allglobal::lmns, allglobal::lplasmaregion, inputlist::lsvdiota, inputlist::ltransform, allglobal::lvacuumregion, allglobal::mn, allglobal::mpi_comm_spec, inputlist::mu, allglobal::myid, allglobal::ncpu, inputlist::nppts, inputlist::nptrj, allglobal::nt, allglobal::ntz, inputlist::nvol, allglobal::nz, inputlist::odetol, allglobal::ofmn, constants::one, fileunits::ounit, constants::pi2, pp00aa(), allglobal::sfmn, tr00ab(), volume(), allglobal::whichcpuid(), inputlist::wmacros, xspech(), and constants::zero.

Referenced by xspech().

Here is the call graph for this function:

Here is the caller graph for this function: