|------------------------------------------------------------------------------ | | Peter Richter (prichter@princeton.edu) | Princeton Plasma Physics Laboratory | August 2000 ! ! Updated Mar02, Rob Andre PPPL: For Impurity Reaction Rates | | documentation.txt - Documentation for PREACT module. | |------------------------------------------------------------------------------ Contents: I. Abstract/Introduction II. Installation III. Usage IV. Error Handling V. Programming Interface VI. Available Data VII. Simple Examples VIII. Known Problems IX. References See README for instructions on installing the module. ------------------------------------------------------------------------------- I. ABSTRACT/Introduction ------------------------------------------------------------------------------- The PREACT module performs lookups and interpolation of tables of the rate (weighted product of cross-section and velocity "sigma*v") of various charge exchange, impact ionization, and fusion reactions. It has independent interfaces for use in both C++ and Fortran-77 codes. The C++ interface is fully object-oriented, and the Fortran-77 interface provides a layer above this C++ interface so that PREACT may be embedded in legacy Fortran-77 transport codes. The module generates data tables with real*8 precision, and for convenience supports both a real and a real*8 programming interface. The interface is "vector oriented", returning a list of reaction rate coefficients given a list of input parameters appropriate to the type of data sought; taking advantage of the vector feature can significantly benefit applications performance. PREACT is now used for all atomic and nuclear rate data in TRANSP. Installing PREACT and exploiting its vector capabilities resulted in a 10-15% speed-up of TRANSP's serial Monte Carlo fast ion calculation. Available reaction table groups are: CX -- (neutralizing) charge exchange between light neutrals and light fully stripped ions (H,He,Li) (Li data is not complete) II -- impact ionization + non-neutralizing charge exchange ionization of light neutrals by light ions (H,He,Li) (Li data is not complete) FS -- various fusion reactions involving D,T,He3 EI -- electron impact ionization and H+,e- recombination SV -- stopping on fully stripped light impurities (C+6,O+8, scaled by Z for other impurities) A partial list of tables: CX II FS EI SV [sigma*v](Erel/A) x x x [sigma*v](Erel) x (Te) x (Eb/Ab,Ti/Ai) x x (Eb,Ti) x <>(Ti) x [sigma*v] indicates a mono-energetic interaction; indicates an average over a Maxwellian; <> indicates a bimaxwellian convolution (i.e. thermonuclear rates) all energies and temperatures are given in KeV or KeV/AMU; output rate coefficients are in m**3/sec. Also available: CX -- tables giving expectation values for energy and normalized velocity (relative to incident beam) for beam-target interactions FS -- gyro-averaged sigma*v table, useful for fast ion - fast ion fusion rate calculations in the presence of a magnetic field Author: Peter Richter (prichter@princeton.edu) Contact: Doug McCune (transp_support@pppl.gov) Update RGA (1/2002): Rates for impurity ionization, recombination and radiation have been added to the preact C++ and fpreact fortran interface. The interface is similar to the other reaction cross sections except for the following aspects: 1) The rates are returned for all of the charge states of an atomic species so the vector interfaces return results in a two dimensional array. 2) The electron density is an input along with the electron temperature with the resulting rates in the form in units of (sec-1) for ionization and recombination or (joules/sec) for radiation. 3) The f77 layer exposes only a small fraction of the underlying C++ code. 4) The initial source of impurity rate data has been taken from the ADPAK module which is now a dependency of preact. ADPAK is Russell Hulse's Atomic Data package extracted from the MIST code. Other sources are expected to follow. 5) The structure of the code is separate from the other reaction code and data table architecture. In addition, there is a routine for calculating the coronal equilibrium using the ionization and recombination rate results. ------------------------------------------------------------------------------- II. INSTALLATION ------------------------------------------------------------------------------- See README for instructions on building, testing and installing the module. ------------------------------------------------------------------------------- III. USAGE ------------------------------------------------------------------------------- The environment variable PREACTDIR must be set to be the full pathname of the directory where the Plasma Reaction Tables are stored. A "high level" fortran interface, fpreact, exists and should be used in preference to direct preact calls in most cases, when calling preact from fortran code. From C++ code, the user interfaces with the PREACT module via the six C++ classes -- CXReaction (for charge exchange reactions) -- IIReaction (for ion impact ionization reactions) -- FSReaction (for fusion reactions) -- EIReaction (for electron impact ionization and recombination reactions (vs. Te only)) -- SVReaction (for impurity stopping cross section tables (vs. Erel only)) -- ImpReaction (for impurity ionization, recombination and radiation) The header files for these classes are located in the directory include/preact, and the file preact/testpreact.cpp shows how the classes are used in a client program. When compiling a client program that uses the PREACT module, the user must link his code with the PREACT libraries (libpreact.a, libsigsub.a, libaladdinsub.a, libr8slatec.a, libadpak.a, libcppsub.a, libcomput.a and libportlib.a) and the Fortran run-time libraries (e.g. libfor.a, -lUfor, and -lFutil on a DEC Alpha). Note that, as always, the order in which the libraries are given is important. Thus, the compile command for a C++ client program 'client.cpp' that uses the PREACT module would look like % g++ -I$PREFIX/include -o client client.cpp \ -L$PREFIX/lib/ -lpreact -lsigsub -laladdinsub -lr8slatec \ -ladpak -lcppsub -lcomput -lportlib /usr/lib/libfor.a -lUfor -lFutil on a DEC Alpha (PREFIX points to the NTCC software root directory). Although preact itself contains a fortran interface, this has been extended for convenience in the library "fpreact.a". "fpreact" uses an f90 static module to store (and hide from the user) the integer "handles" needed by fortran to reference the underlying c++ objects. The test program "fpreact_test.f90" demonstrates use of both the real and real*8 interfaces. The f90 interface module "fpreact_calls" can be used to have the compiler select the right fpreact routine depending on the precision of calling arguments. For the impurity rate data, the f77 and f90 interfaces use the atomic number of the impurity element to reference the underlying c++ object instead of using an integer handle as is used by the other reactions in fpreact. In order to run a fortran client program on an Alpha, one would first type something like % f90 -c client.for to create object files for each of the client program source files, and would then type something like % f90 -I$PREFIX/include -o client client.o \ -L$PREFIX/lib -lfpreact -lpreact -lsigsub -laladdinsub \ -lr8slatec -ladpak -lcppsub -lcomput -lportlib \ -lstdc++ to build the client program (note the trailing reference to the c++ runtime library may take different forms on different systems -- also not for gcc you might need to include libgcc.a). To view the current parameters for a particular reaction's tables, type the command: % preact_list The 'preact_list' program will prompt the user for the reaction and then show the parameters for the given reaction's tables. To change the parameters for a particular table (which will require recomputing the table unless one with those same parameters already exists), type: % preact_change The 'preact_change' program will prompt the user for the necessary information and then change the table parameters for the appropriate reaction. To remove all existing computed tables from the database (thus re-initializing the database), type: % preactinit These are the only commands that the user need be concerned with. (The 'preact_init' program in the directory $PREFIX/bin is not to be called explicitly by the user; it is run by 'gmake install' and 'preactinit' script.) ------------------------------------------------------------------------------- IV. Error handling ------------------------------------------------------------------------------- The main type of error handled by the preact/fpreact software is an out-of-range request for table interpolation. The default error handling is: warn and use low end table value for out-of-range low; write message, set status code and exit for out-of-range high. The error handling "policy" is user controllable, with the following policies available and indepently controllable on each dimension of each table (from fpreact/hatom_mod.f90): ! "policy options" for range errors on table lookups. ! for each initialization a "policy" must be set for two ! conditions: out-of-range low, and out-of-range high ! -- can also modify independently for each dimension of each table... integer, parameter :: r_ignore=0,r_warn=1,r_halt=2 integer, parameter :: r_ignore_clear=3,r_warn_clear=4 ! ! These policies have effect when an energy or temperature argument ! is given beyond the range of the table data: ! ! r_ignore means: ignore range error; use nearest table value ! r_warn means: use nearest table value, but issue a warning message ! r_halt means: issue a warning message and exit ! r_ignore_clear means: set result to zero *silently* ! r_warn_clear means: set result to zero with warning message ! ! actual policies: ! default: ignore if out-of-range low ! halt with message if out-of-range high The fortran interface fpreact sets "reasonable" error handling policies; the code for this is in the private module fpreact/hatom_mod.f90. The C++ interface allows either the default policy or a single policy for all table dimensions to be defined at object instantiation time; policy for individual dimensions can then be modified by subsequent public "policy" method invocations. ------------------------------------------------------------------------------- V. Programming Interfaces (C++ and Fortran) ------------------------------------------------------------------------------- The C++ interface for the PREACT module consists of the following publiclly accessible methods: ***** CXReaction class methods: ************************************ #include "grid_1d.h" #include "grid_2d.h" class CXReaction { public: // constructor, destructor... CXReaction(char tag[], int errlo = WARN, int errhi = HALT); ~CXReaction(); // sigma, vrel as functions of Erel... fp_64 sigma(fp_64 Erel); fp_64 vrel(fp_64 Erel); // sigma*v for beam to cold (T=0) target... void coldTarget_sigmaTimesV(fp_64 E[], fp_64 f[], int n); void coldTarget_sigmaTimesV_32(fp_32 E[], fp_32 f[], int n); range coldTarget_ERange() const; void coldTarget_EPolicy(int errlo, int errhi); // for beam to warm (T>0) target... void warmTarget_sigmaTimesV(fp_64 E[], fp_64 T[], fp_64 f[], int n); void warmTarget_sigmaTimesV_32(fp_32 E[], fp_32 T[], fp_32 f[], int n); range warmTarget_ERange() const; range warmTarget_TRange() const; void warmTarget_EPolicy(int errlo, int errhi); void warmTarget_TPolicy(int errlo, int errhi); // of target reagent... void energyAvg_perAMU(fp_64 E[], fp_64 T[], fp_64 f[], int n); void energyAvg_perAMU_32(fp_32 E[], fp_32 T[], fp_32 f[], int n); range energyAvg_perAMU_ERange() const; range energyAvg_perAMU_TRange() const; void energyAvg_perAMU_EPolicy(int errlo, int errhi); void energyAvg_perAMU_TPolicy(int errlo, int errhi); // ... void average_vDotNorm(fp_64 E[], fp_64 T[], fp_64 f[], int n); void average_vDotNorm_32(fp_32 E[], fp_32 T[], fp_32 f[], int n); range average_vDotNorm_ERange() const; range average_vDotNorm_TRange() const; void average_vDotNorm_EPolicy(int errlo, int errhi); void average_vDotNorm_TPolicy(int errlo, int errhi); ...} ***** IIReaction class methods: ************************************ #include "grid_1d.h" #include "grid_2d.h" class IIReaction { public: // constructor, destructor... IIReaction(char tag[], int errlo = WARN, int errhi = HALT); ~IIReaction(); // sigma, vrel as functions of Erel... fp_64 sigma(fp_64 Erel); fp_64 vrel(fp_64 Erel); // sigma*v for beam to cold (T=0) target... void coldTarget_sigmaTimesV(fp_64 E[], fp_64 f[], int n); void coldTarget_sigmaTimesV_32(fp_32 E[], fp_32 f[], int n); range coldTarget_ERange() const; void coldTarget_EPolicy(int errlo, int errhi); // for beam to warm (T>0) target... void warmTarget_sigmaTimesV(fp_64 E[], fp_64 T[], fp_64 f[], int n); void warmTarget_sigmaTimesV_32(fp_32 E[], fp_32 T[], fp_32 f[], int n); range warmTarget_ERange() const; range warmTarget_TRange() const; void warmTarget_EPolicy(int errlo, int errhi); void warmTarget_TPolicy(int errlo, int errhi); ...} ***** FSReaction class methods: ************************************ #include "grid_1d.h" #include "grid_2d.h" class FSReaction { public: // constructor, destructor... FSReaction(char tag[], int errlo = WARN, int errhi = HALT); ~FSReaction(); // sigma, vrel as functions of Erel... fp_64 sigma(fp_64 Erel); fp_64 vrel(fp_64 Erel); // sigma*v for beam to cold (T=0) target... void coldTarget_sigmaTimesV(fp_64 E[], fp_64 f[], int n); void coldTarget_sigmaTimesV_32(fp_32 E[], fp_32 f[], int n); range coldTarget_ERange() const; void coldTarget_EPolicy(int errlo, int errhi); // for beam to warm (T>0) target... void warmTarget_sigmaTimesV(fp_64 E[], fp_64 T[], fp_64 f[], int n); void warmTarget_sigmaTimesV_32(fp_32 E[], fp_32 T[], fp_32 f[], int n); range warmTarget_ERange() const; range warmTarget_TRange() const; void warmTarget_EPolicy(int errlo, int errhi); void warmTarget_TPolicy(int errlo, int errhi); // <> thermonuclear... void thermonuclear_sigmaTimesV(fp_64 T[], fp_64 f[], int n); void thermonuclear_sigmaTimesV_32(fp_32 T[], fp_32 f[], int n); range thermonuclear_TRange() const; void thermonuclear_TPolicy(int errlo, int errhi); // gyro average... void gyroAvg_sigmaTimesV(fp_64 E[], fp_64 Var[], fp_64 f[], int n); void gyroAvg_sigmaTimesV_32(fp_32 E[], fp_32 Var[], fp_32 f[], int n); range gyroAvg_ERange() const; range gyroAvg_VarRange() const; void gyroAvg_EPolicy(int errlo, int errhi); void gyroAvg_VPolicy(int errlo, int errhi); ...} ***** EIReaction class methods: ************************************ #include "grid_1d.h" #include "grid_2d.h" class EIReaction { public: // constructor, destructor... EIReaction(char tag[], int errlo = WARN, int errhi = HALT); ~EIReaction(); // sigma, vrel as functions of Erel... fp_64 sigmav(fp_64 Te); // sigma*v for beam to cold (T=0) target... void sigvTe_sigmaTimesV(fp_64 Te[], fp_64 f[], int n); void sigvTe_sigmaTimesV_32(fp_32 Te[], fp_32 f[], int n); range sigvTe_TRange() const; void sigvTe_TPolicy(int errlo, int errhi); ...} ***** SVReaction class methods: ************************************ #include "grid_1d.h" #include "grid_2d.h" class SVReaction { public: // constructor, destructor... SVReaction(char tag[], int errlo = WARN, int errhi = HALT); ~SVReaction(); // sigma, vrel as functions of Erel... fp_64 sigmav(fp_64 Erel); // sigma*v for beam to cold (T=0) target... void coldTarget_sigmaTimesV(fp_64 E[], fp_64 f[], int n); void coldTarget_sigmaTimesV_32(fp_32 E[], fp_32 f[], int n); range coldTarget_ERange() const; void coldTarget_EPolicy(int errlo, int errhi); ...} ***** ImpReaction class methods: ************************************ Two classes are relavent for impurity reactions. A React class defines a single quantity (ionization, recombination or radiation). Impurity rates are returned by operating on the React objects. An ImpReaction class holds the four Reacts needed to describe the rates of a single impurity species (ionization, recombination, radiation (with or without bremsstrahlung), bremsstrahlung radiation). Also relevant is the fwrap template class which provides access to memory like a multidimensional fortran array and is typedefed to the names ap_32 and ap_64 depending on the precision of the result. Two static methods exist in ImpReaction to support the calculation of the coronal equilibrium and the radiation for a coronal equilibrium. These are template functions which will select the precision according to the type of the rate arguments. // // ------------------------ ImpReaction ----------------------- // Encapsulates the reactions associated with impurity rates // Currently includes: // - electron collisional ionization rate (sec-1) // - recombination rate due to radiative and dielectronic // recombination (sec-1) // - radiation in (joules/sec) due to collisional excitation // and recombination and optionally brem...(default is no brem...) // - brem... radiation only in (joules/sec) // // For all three rates: // arguments: 0-> Te (keV) electron temperature // 1-> Ne (m-3) electron density // // result: rank 1 with shape [Z+1] // [i]-> rate for charge state i, i=0...Z // // Note: If brem... is set to be included in radiation rad() React, then // trying to call brem() will throw an exception. This is only to // prevent double counting of brem... in total radiation. // // example: // ImpReaction r(6) ; // build carbon // ap_32 xarg(2), result(7) ; // xarg(0) = 10. ; // Te in keV // xarg(1) = 1.e19 ; // Ne in m-3 // r.ioniz()(xarg,result) ; // fill result array // class ImpReaction { public: // // --- constructor --- // The option value is used to select the objects used to implement // the impurity rates. // enum option {ADPAK_DIRECT=0} ; ImpReaction(int nucz, option r=ADPAK_DIRECT) ; ~ImpReaction() ; // --- reaction data --- int nucz() const ; // atomic number of this element string name() const ; // name of this reaction string info() const ; // informative // // -- modify the out of range policy for all three Reacts -- // Only the policy is modifiable here because the minimum and maximum // should be intrinsic to the implementation method. The clear_minimum // and clear_maximum values are set to the default of 0.0 by these calls. // Changing the policy could result in rebuilding an internally cached table. // i=0 -> Te // =1 -> Ne void r_set_min_policy(size_t i, Range::policy pmin) ; // set policy for axis Range i void r_set_max_policy(size_t i, Range::policy pmax) ; // --- Reacts --- React ioniz() ; // Ionization (sec-1) React recom() ; // Recombination (sec-1) React rad() ; // Radiation (joules/sec) React brem() ; // Brem (joules/sec) // --- static --- static string table_name(int knucz) ; // return name in the form "Helium (He)" // // ----- static coronal equilibrium ----- // These functions can operate in either single evaluation or vector // form depending on the rank of the input arrays // // coronal: calculates the coronal equilibrium given the ionization // and recombination rate // single mode, // ioniz[Z+1] -> input -- rank 1 array of ionization rate over // charge states 0...Z // recom[Z+1] -> input -- rank 1 array of recombination rate over // charge states // coron[Z+1] -> output -- rank 1 array of fraction of particles in // each charge state. This is normalized // to a sum of 1. // vector mode, // ioniz[Z+1,N] -- input -- same but rank 2 for N different sets of data // recom[Z+1,N] -- input // coron[Z+1,N] -- output // // coronal_rad: returns the total radiation in (joules/sec) based on the // supplied coronal_equilibrium and radiation rate // single mode, // coron[Z+1] -> input -- rank 1 array of fraction of particles in // each charge state as returned by coronal(), need not // be normalized to a sum of 1. // rad[Z+1] -> input -- rank 1 array of radiation rate in (joules/sec) // total[] -> output -- rank 0 array of total radiation in (joules/sec) // // vector mode, // coron[Z+1,N] -- input -- same but rank 2 for N different sets of data // rad[Z+1,N] -- input -- same but rank 2 for N different sets of data // total[N] -- output -- same but rank 1 for N different sets of data // template static void coronal(const fwrap& ioniz, const fwrap& recom, fwrap& coron) ; template static void coronal_rad(const fwrap& coron, const fwrap& rad, fwrap& total) ; } ; // // ---------------------------- React ------------------------ // Describes a single function in a reaction. This object has // reference semantics so copying it will refer to the same // underlying reaction (aReact). The function it describes has a fixed // number of arguments, nargs(), and can return an array as // a result of its invocation. // class React { public: // --- object stuff --- React(aReact* q) ; React(const React& q) ; React& operator=(const React& q) ; ~React() ; // --- arguments --- size_t nargs() const ; // number of arguments of function fp_64 min(size_t i) const ; // return the minimum of this argument range 0<=i