- Acknowledgements
- Introduction
- System requirements
- Quick-start
- Modifying the impurity database
- SD1D Integration
The majority of this code is based on the atomic1D code, which is in turn based on the excellent OpenADAS analysis tool provided at cfe316/atomic.
The data used for modelling is entirely supplied by the OpenADAS project, a public-access database managed by the University of Strathclyde.
This project provides a C++
library for the analysis of atomic processes in a SOL/divertor plasma. It is primarily intended as a library for integration into the BOUT++ SD1D project, although the testing suites provided may be generally useful.
The project consists of three main parts:
- The
atomicpp
library -- provides all analysis code as well the Cython.pyx
wrapper; Prad.cpp
testing suite inC++
-- tests basic functionality ;Prad.py
verification suite inPython3
-- uses a Python wrapper to theatomicpp
library to check the results from the library.
The project is built and run using a Makefile
architecture. All commands are localised to the project directory -- to uninstall simply delete the project folder.
It is recommended that a package manager (apt-get
(Linux) or MacPorts
(Mac), etc) be used to install the required packages.
For the core atomicpp
library and the Prad.cpp
C++ test suite
- A C++ compiler which supports C++11 standard (tested with GNU compiler
gcc/g++ version 6.3.0
) gmake
to run theMakefile
(tested withGNU Make 3.81
)
For the Prad.py
Python3 verification suite
- A Python3 installation with an installed SciPy stack and Cython. Both the SciPy stack and Cython are available in the Anaconda distribution (tested with
Python 3.6.1 |Anaconda 4.4.0 (x86_64)
)
To extend the json_database
of OpenADAS rate-coefficients
- The OpenADAS_to_JSON project. See Modifying the impurity database and the OpenADAS_to_JSON
README
for more details. - A Fortran compiler (tested with GNU compiler
gcc/gfortran version 6.3.0
)
- The
atomicpp
module folder containingRateEquations
(.cpp
and.hpp
): for evaluation of the population, power and momentum balance derivatives.ImpuritySpecies
(.cpp
and.hpp
): Called byRateEquations
, for storing the data of a specified plasma impurity.RateCoefficient
(.cpp
and.hpp
): called byImpuritySpecies
, for providing a callable interface to a set of OpenADAS density- and temperature-resolved rate coefficients.sharedFunctions
(.cpp
and.hpp
): for functions which are useful to all other modules - currently a JSON reader and a file checker.json.hpp
: header-only JSON reader from the nlohmann::json project. A copy of the header file is included, although it is recommended that users download an up-to-date version here.SD1DData
(.cpp
and.hpp
): not required for SD1D integration. For storing and processing data from an SD1Dcollect
, which is then saved as a json file. Run theatomic1D/reference/data_dict_export.py
program in the SD1D output folder to produce the json file, and then copy this to theatomic++
directory (wherePrad.cpp
is being executed). This data is used to train and verify theatomic++
code.
impurity_user_input.json
(same directory asPrad.cpp
): a plain-text file for providing hard-coded data relating to the impurity (such as year for OpenADAS data, etc). To modify, see the JSON project for nomenclature and follow the same style as the template.sd1d-case-*.json
(same directory asPrad.cpp
): a json file created byatomic1D/reference/data_dict_export.py
, for use withSD1DData
for training the program (not needed if integrating into SD1D).
The Makefile
provides most of the desired functionality of the project. It has 5 main commands;
make cpp
: checks if the required source files have been built since they were last modified -- if not, compiles theatomicpp
library and thePrad.cpp
program into.o
files and then links the.o
files into an executable (Prad
).make cpp_run
: (default formake
) performs themake cpp
functionality and also runs thePrad
executable via./Prad
make py
: checks if the required source files have been built since they were last modified -- if not, compiles theatomicpp
library and generates a Python moduleatomicpy
which is declared inatomicpy.pyx
(Cython) and built withsetup.py build_ext --inplace
.make py_run
: performs themake py
functionality and also runs thePrad.py
script viapython Prad.py
make clean
: reverts the project to a fresh install state.
A separate project is supplied at OpenADAS_to_JSON. This project downloads ADF11 files for the specified impurity from OpenADAS, uses Fortran helper functions to read the data and exports the data as JSON files. To change which species is being considered you'll need to do the following;
- Modify the
elements
tag of themakefile
- Run
make json_update
- Copy the
OpenADAS_to_JSON/json_database
file from that project ontoatomicpp/json_database
(overwrite) - Update
impurity_user_input.json
to include the impurity data required by the program - Update the
impurity_symbol_supplied
variable of Python and C++ if using the supplied testing programs.
Sub-contents
- Impurity rate equations
- Hydrogen rate equations
- The
DerivStruct
data structure - Print method for the data structure
- The
computeDerivs
function - Shared interpolation
- Kahan-Neumaier summation
N.b. the SD1DData (.hpp and .cpp) files of the atomicpp module directory are not required
The principal purpose of this code is to extend the radiation model of the SD1D SOL/divertor plasma simulation code, which is built on the BOUT++ project.
The derivative evaluator is provided in the RateEquations
class. This class is initialised from a ImpuritySpecies
object, shown here for a Carbon (c
) impurity. The derivative evaluator is supplied as computeDerivs
-- see the RateEquations.hpp
header for identity and units of the arguments.
#include "atomicpp/ImpuritySpecies.hpp"
#include "atomicpp/RateEquations.hpp"
std::string impurity_symbol="c";
atomicpp::ImpuritySpecies impurity(impurity_symbol);
atomicpp::RateEquations impurity_derivatives(impurity); //Organised as a RateEquations object for cleanliness
impurity_derivatives.setThresholdDensity(1e9); //Density threshold - ignore ionisation stages which don't have at least this density
impurity_derivatives.setDominantIonMass(1.0); //Dominant ion mass in amu, for the stopping time calculation
atomicpp::DerivStruct derivative_struct = impurity_derivatives.computeDerivs(Te, Ne, Vi, Nn, Vn, Nzk, Vzk);
The 'impurity' methods for evaluation of the OpenADAS rate-coefficient data can be equally applied to the dominant ion. A modified method computeDerivsHydrogen
is supplied to avoid having to supply the same variables as the 'impurity' and 'dominant ion' values, although computationally the method is almost identical (ion-ion drag and charge-exchange are neglected). The equivalent case is
#include "atomicpp/ImpuritySpecies.hpp"
#include "atomicpp/RateEquations.hpp"
std::string hydrogen_symbol="h";
atomicpp::ImpuritySpecies hydrogen(hydrogen_symbol);
impurity_derivatives.setThresholdDensity(1e9); //Density threshold - ignore ionisation stages which don't have at least this density
atomicpp::DerivStruct derivative_struct_H = hydrogen_derivatives.computeDerivsHydrogen(Te, Ne, Nhk, Vhk);
The computeDerivs
functions return a DerivStruct
data structure. This is defined as
struct DerivStruct{
double Pcool;
double Prad;
std::vector<double> dNzk;
std::vector<double> F_zk;
double dNe;
double F_i;
double dNn;
double F_n;
};
which may be unpacked as
double Pcool = derivative_struct.Pcool; //Electron-cooling power, in J m^-3 s^-1 (needed for electron power balance)
double Prad = derivative_struct.Prad; //Radiated power, in in J m^-3 s^-1 (for comparing to diagnostic signal)
std::vector<double> dNzk = derivative_struct.dNzk; //Change in each ionisation stage of the impurity population, in particles m^-3 s^-1
std::vector<double> F_zk = derivative_struct.F_zk; //Force on each particle of ionisation stage k of the impurity population, in N
double dNe = derivative_struct.dNe; //Perturbation change in the electron density (in particles m^-3 s^-1) and
double F_i = derivative_struct.F_i; // perturbation force (in N) on the electron population due to atomic processes
double dNn = derivative_struct.dNn; //Perturbation change in the neutral density (in particles m^-3 s^-1) and
double F_n = derivative_struct.F_n; // perturbation force (in N) on the neutral population due to atomic processes
The RateEquations
class also includes a printing method for the DerivStruct
data structure
impurity_derivatives.printDerivativeStruct(derivative_struct);
which returned a formatted result
The core functionality of the RateEquations
object is provided by the computeDerivs
function. This is given as
DerivStruct RateEquations::computeDerivs(
const double Te,
const double Ne,
const double Vi,
const double Nn,
const double Vn,
const std::vector<double>& Nzk,
const std::vector<double>& Vzk){
resetDerivatives(); //Reset all the derivatives to zero, since the object has a memory of the previous step
// Perform 'sharedInterpolation' - find the lower-bound gridpoint and fraction into the grid for both Te and Ne
std::pair<int, double> Te_interp, Ne_interp;
if (use_shared_interpolation){
Te_interp = findSharedInterpolation(rate_coefficients["blank"]->get_log_temperature(), Te);
Ne_interp = findSharedInterpolation(rate_coefficients["blank"]->get_log_density(), Ne);
} else {
throw std::runtime_error("Non-shared interpolation method requires switching of method. Declare Te_interp and Ne_interp as doubles instead of <int, double> pairs.");
// //Pass Te_interp and Ne_interp as doubles instead of pairs and the program will auto-switch to the full interpolation method.
// double Te_interp = Te;
// double Ne_interp = Ne;
}
calculateElectronImpactPopulationEquation(Ne, Nzk, Vzk, Te_interp, Ne_interp);
calculateChargeExchangePopulationEquation(Nn, Nzk, Vzk, Te_interp, Ne_interp);
//Apply neumairSum corrections
for(int k=0; k<=Z; ++k){
dNzk[k] += dNzk_correction[k];
F_zk[k] += F_zk_correction[k];
}
verifyNeumaierSummation();
calculateIonIonDrag(Ne, Te, Vi, Nzk, Vzk);
calculateElectronImpactPowerEquation(Ne, Nzk, Te_interp, Ne_interp);
calculateChargeExchangePowerEquation(Nn, Nzk, Te_interp, Ne_interp);
// auto derivative_tuple = makeDerivativeTuple();
// return derivative_tuple;
DerivStruct derivative_struct = makeDerivativeStruct();
return derivative_struct;
};
computeDerivs
relies 'shared interpolation' to calculate the scaling factors for interpolation, which requires that the underlying grids are identical
/**
* @brief find the lower-bound gridpoint and fraction within the grid for the given point at which to interpolate
* @details Using bilinear interpolation, the scaling factors for interpolating the rate coefficients are the same
* regardless of which process is called (since the underlying log_temperature and log_density grids are the same).
* Therefore, the grid-point and fraction pair may be shared by any rate-coefficient. Have overloaded call0D such that,
* if a <int, double> pair is supplied as an argument then the shared interpolation method will be called
*
* @param log_grid grid-points for which data is given
* @param eval point at which the interpolation should be performed
*
* @return <int, double> pair where int is the lower-bound grid-point and fraction is the scaling factor (fractional distance
* between the lower and upper-bound gridpoints)
*/
std::pair<int, double> findSharedInterpolation(const std::vector<double>& log_grid, const double eval){
// Perform a basic interpolation based on linear distance
// values to search for
...
std::pair<int, double> interp_pair(interp_gridpoint, interp_fraction);
return interp_pair;
}
To add the elements of a list with significantly varying orders of magnitude to very high precision (avoids floating point rounding error) the Kahan-Neumaier algorithm is implemented as follows;
/**
* @brief Uses Neumaier algorithm to add the elements of a list
* @details Extension on Kahan summation algorithm for an unsorted list
* Uses a compensated sum to improve precision when summing numbers of
* significantly different magnitude
*
* @param list_to_sum The list of numbers to sum
* @return neumaier_pair the uncompensated sum and the compensation
* Compensated sum is sum + correction - however, this is left external
* in case summation is restarted
*/
std::pair<double, double> neumaierSum(const std::vector<double>& list_to_sum, const double previous_correction = 0.0){
...
std::pair<double, double> neumaier_pair(sum, correction);
return neumaier_pair;
}