The MORTY C/Fortran ORIGEN API

This section describes the available features and use of a limited-scale API called MORTY (the MELCOR ORIGEN TidY interface). The MORTY API was originally developed to link ORIGEN directly into the MELCOR code for severe accident consequence analysis; however, it can be used in a broad set of applications for which the direct use of the ORIGEN solver through embedded calls to ORIGEN methods is desired.

The actual design of MORTY is a C++ “shim” to a core set of ORIGEN capabilities relating to specifying materials, moving materials between logical units, and obtaining quantities such as elemental masses (in kilograms, kg) and decay heat (in watts, W). By utilizing primitive data types such as integers, floating-point real numbers, and strings over ORIGEN’s native data types, the MORTY API streamlines user interaction with ORIGEN, prioritizing user-friendliness and accessibility. This design choice facilitates seamless interoperability between C++ and Fortran, enabling users to efficiently utilize ORIGEN’s capabilities. With this approach, users can effectively operate ORIGEN without requiring a deep understanding of its intricate methods and data structures.

This section is divided into three components: the first section describes the functions available in the MORTY API and their usage; the second section briefly describes linking ORIGEN via shared libraries from external applications (to call the MORTY API without the need to recompile the SCALE source code); and the third section provides demonstrative examples of how to use the MORTY interface to perform decay calculations.

MORTY API Functions and Uses

The MORTY API features a unique global singleton instance (referred to as the global MORTY) which maintains a persistent state across function calls. This instance holds information such as the current step index and the set of depletion materials being tracked in the problem. However, the user need not ever interact with the global MORTY singleton, as all interfaces are handled through the provided functions.

Caution

As a global singleton object, the MORTY API is by design not thread-safe and thus should not be used in multi-threaded applications.

Conventions

The MORTY API uses two important conventions for returning status information to the user: the data types STATUS and MAT_ID, both of which are mapped to C++ integer types.

STATUS is an integer used to convey the return status of a function; it is typically associated with the current step index of the current operation or zero. Negative numbers for status always indicate an error.

MAT_ID is an integer used to convey the material ID corresponding to an operation; for example, the material ID of the newly created material is returned when initializing a material; it is also used when setting properties on materials (e.g., setting initial mass, transferring mass between materials).

Description of MORTY Functions (C/C++)

void morty_clear()

Clear the global MORTY singleton object (remove all materials / timesteps)

bool morty_has_material(int mat_num)

Returns true if the MORTY global singleton contains the requested material

MAT_ID morty_init_material(const std::string &filename)

Initialize a new material using the decay library at the provided filepath. Note that ORIGEN aliases (e.g., decaylib.data) could be used in place of a fully qualified filepath.

Returns the material ID of the created material if successful, 0 if unsuccessful. If the user-provided library could not be found but the backup (default ORIGEN decay library) was successfully loaded, returns -(material ID).

int morty_set_default_substeps(int nsteps)

Set the default number of substeps to be used for calculating inventories for each major “step.” (Default is 10.)

std::string morty_get_error()

Get the most recent error message stored on the global MORTY singleton

Because MORTY is intended for embedding into separate applications, direct output of errors to the standard error / out streams is not necessarily desirable. Thus, users can retrieve the most recent error from the stack via the morty_get_error() function.

Returns an empty string if the error buffer is empty.

int morty_nuclides_size(MAT_ID mat_num)

Get the size of the nuclide vector for material mat_num.

STATUS morty_get_nuclides(MAT_ID mat_num, int *zaids, int zaids_size)

Gets the full list of nuclide IDs in IZZZAAA format for material number mat_num. Returns 0 upon success.

STATUS morty_load_initial_mass_from_file(MAT_ID mat_num, const std::string &response, const std::string &filename, int pos)

Load the initial isotopic vector for material mat_num from an external ii.json file. response is the name of the response to be loaded, filename is the fully qualified path to the ii.json file, and pos is the position of response (i.e., in time) to load.

Returns:

• -mat_num if the material does not exist on the global MORTY singleton or if loading concentrations from the ii.json object fails

• -1 if the interface file is not found.

• 0 on success.

STATUS morty_move_to_s(double dt)

Move to cumulative time dt in the calculation, in seconds. Note that MORTY always works in units of _absolute_ time; ergo, to move backward in time one need only specify a smaller value of dt than the current maximum time.

Returns -1 if unsuccessful, otherwise returns the step index of the current timestep.

STATUS morty_total_time_s(double &t)

Get the total (cumulative) time across all timesteps. Returns the _total_ number of timesteps.

STATUS morty_move_mass_kg(MAT_ID source, MAT_ID target, const int *z_list, int zlist_size, double *masses_kg, int masses_size)

Move masses (in kg) from material source to target. Elements are specified by element number (e.g., 92 for uranium). Two parallel arrays specify the element and mass to be transferred. Isotopes are proportionally moved according to their mass fraction by element on the source material.

Returns:

• -source if the source material does not exist or the array sizes do not match

• -target if the target material does not exist

• z if the element is not a valid element ID

• -index if the user requests to move more mass of an element from source than is present, the index of the element in the list

• 0 if successful

STATUS morty_get_decay_heat_w(MAT_ID mat_num, const int *z_list, int zlist_size, double *dh, int dh_size)

Get the decay heat by element for material number mat_num at the current time. Decay heats are calculated per-element based on the elements provided (by atomic number) in the z_list array. Returns -mat_num if there is a problem (material does not exist, array sizes do not match, other problems) and 0 upon success.

MORTY function signatures (Fortran)

The MORTY API is primarily oriented for generating useful and easy-to-use function signatures translated into Fortran bindings via SWIG. Certain conventions that may appear unwieldy or cumbersome in C++ (i.e., pointers with array dimensions) are intended to be wrapped into convenient Fortran wrapper functions that accept Fortran-native types with ISO_C bindings (e.g., C_INT, C_DOUBLE.).

The following list describes each of the function signatures as they are accessed in Fortran.

morty_init_material

Arguments:

• libname: Character array containing the fully qualified library path

morty_set_default_substeps

Arguments:

• nsteps: C-type integer with the number of substeps to use per step

morty_get_error

Arguments:

• None

Returns variable-length, null-terminated C-type character array.

(Note: error message is automatically truncated to the length of the character buffer provided.)

morty_nuclides_size

Arguments:

• mat_num: C-type integer corresponding to the material number

Returns C-type integer of the number of nuclides tracked on the material

morty_load_initial_mass_from_file

Arguments:

• mat_num: C-type integer corresponding to the material number to load inventories into

• response: Character array corresponding to the response name to load

• filename: Character array containing the fully qualified filepath of the ii.json file to load

• pos: C-type integer corresponding to the position to load

• norm: (Optional) Total mass (in kg) to normalize loaded concentrations to; otherwise initial mass is inferred from the composition on the ii.json file.

morty_save_mass_to_file

Arguments:

• mat_num: C-type integer corresponding to the material number from which to save inventories

• response: Character array label for the response to be saved

• filename: Character array of a fully qualified path of the ii.json file to save to

morty_set_initial_mass_kg

Arguments:

• mat_num: C-type integer corresponding to the material number to set initial inventories for

• zaid_list: Assumed-size array of IZZZAAA nuclide IDs (C-type ints)

• masses_kg: Masses for each isotope in the zaid_list, in units of kg (C-double type reals)

morty_move_to_s

Arguments:

• dt: Time to move to, in seconds (C-type double)

morty_total_time_s

Arguments:

• dt: Total time at current step, in seconds (C-type double, intent INOUT)

morty_move_mass_kg

Arguments:

• source: Source material to move mass from (C-type integer)

• target: Target material to move mass to (C-type integer)

• z_list: Assumed-size array of atomic numbers of elements to move (C-type integers)

• masses_kg: Assumed-size array of masses of each element to move, in kg (C-type doubles)

morty_get_mass_kg

Arguments:

• mat_num: Material number to get elemental masses for (C-type integer)

• z_list: Atomic numbers of elements to get masses for (C-type integers)

• masses_kg: Masses (in kg) for each element (C-type doubles)

morty_get_decay_heat_w

Arguments:

• mat_num: Material number to get elemental masses for (C-type integer)

• z_list: Atomic numbers of elements to get masses for (C-type integers)

• dh: Decay heat (in watts) for each element specified (C-type doubles)

Linking to ORIGEN via external applications

To link to ORIGEN via an external application, one needs to include the following SCALE libraries (found in ${SCALE}/lib):

• liborigen.*

• libampx.*

• libs5lib.*

• libs6inp.*

• lib6slib.*

• libnemesis.*

• libnemesis_fortran.*

Where the extension is the appropriate OS-specific shared library extension (e.g., .so, .dylib, etc.).

Note that CMake users can locate these libraries automatically by using the command:

find_package(SCALE)

Additionally, ORIGEN depends upon the following third-party modules (not included with SCALE):

• Trilinos

• HDF5

All of the MORTY C/C++ functions are defined in the MelcorInterface.hh header, found in the top-level include directory, i.e. include/origen/Manager/wrap/MelcorInterface.hh.

MORTY use case examples

Setting and retrieving mass on a material

The following example demonstrates setting a mass of 1 kg of Nd-148 on material 1 and then subsequently retrieving the mass of material 1.

const std::string libname = "decaylib.def";

// return value of 1 => created material #1
STATUS initialized = morty_init_material(libname);

std::vector<int>    zaids    = {60148};
std::vector<int>    z_list   = {60};
std::vector<double> ref_vals = {1.0};

// 1.0 kg Nd-148
int status = morty_set_initial_mass_kg(
    1, zaids.data(), zaids.size(), ref_vals.data(), ref_vals.size());


// Since Nd-148 is stable, we should get back 1.0 kg
int                 mat_vec_size;
std::vector<double> dec_conc_morty(zaids.size(), 0.0);

STATUS success = morty_get_mass_kg(1,
                                   z_list.data(),
                                   z_list.size(),
                                   dec_conc_morty.data(),
                                   dec_conc_morty.size());

The equivalent FORTRAN code would be:

integer(C_INT), allocatable :: zaids(:), zlist(:)
integer(C_INT) :: m_status
real(C_DOUBLE), allocatable :: ref_vals(:), dec_conc(:)

m_status = morty_init_material("decaylib.def")
EXPECT_EQ(1, m_status)

! 1.0 kg of Nd-148
allocate(zaids(1), zlist(1), ref_vals(1), dec_conc(1))

zaids = 60148
zlist = 60
ref_vals = 1.0
m_status = morty_set_initial_mass_kg(1, zaids, ref_vals)

m_status = morty_move_to_s(1000.0d0)

! Should get back 1.0 kg Nd
m_status = morty_get_mass_kg(1, zlist, dec_conc)

Performing a basic decay calculation and retrieving mass / decay heat

The following example sets a series of radioactive materials (Sr-90, Cs-137, and Np-239) on material 1 and then decays them for a period of 3600 seconds (1 hour). The user can then retrieve the element masses (in kilograms) from material 1.

auto gm = global_morty();
gm->clear();
gm->init();

const std::string libname = "decaylib.def";

// sr-90, cs-137, np-239
std::vector<int>    ref_ids  = {38090, 55137, 93239};
std::vector<int>    ref_z     = {38, 55, 93};
std::vector<double> ref_vals = {5.2, 2.1, 10.0};

// return value of 1 => created material #1
STATUS initialized = morty_init_material(libname);

int mat_num = 1;
int status  = morty_set_initial_mass_kg(mat_num,
                                       ref_ids.data(),
                                       ref_ids.size(),
                                       ref_vals.data(),
                                       ref_vals.size());

const double target_time = 3600.0; // 1 hour
morty_move_to_s(target_time);


std::vector<double> dec_conc_morty(ref_vals.size(), 0.0);
std::vector<double> dec_heat_morty(ref_vals.size(), 0.0);

morty_get_mass_kg(1,
                  ref_z.data(),
                  ref_z.size(),
                  dec_conc_morty.data(),
                  dec_conc_morty.size());

morty_get_decay_heat_w(1,
                  ref_z.data(),
                  ref_z.size(),
                  dec_heat_morty.data(),
                  dec_heat_morty.size());

The FORTRAN equivalent is:

integer(C_INT), allocatable :: ref_zaids(:), ref_zaids2(:), ref_z(:), ref_z2(:)
integer(C_INT) :: m_status
real(C_DOUBLE), allocatable :: ref_vals(:), ref_vals2(:), dec_conc(:), dec_consts(:)
real(C_DOUBLE) :: ref_dec, dec_val

allocate(ref_zaids(4), ref_z(4), ref_vals(4))
ref_zaids = (/ 92235, 92238, 53133, 55137 /)
ref_z = (/ 92, 92, 53, 55 /)
ref_vals = (/ 30.0, 970.0, 10.0, 1.0 /) ! masses in kg

call morty_clear()

! Initialize first material
m_status = morty_init_material("decaylib.def")
m_status = morty_set_initial_mass_kg(1, ref_zaids, ref_vals)

allocate(ref_zaids2(3), ref_z2(3), ref_vals2(3), dec_consts(3))
ref_zaids2 = (/38090, 55137, 93239/)
ref_z2 = (/ 38, 55, 93 /)
ref_vals2 = (/ 5.2, 2.1, 10.0 /)
dec_consts = (/ 7.60018E-10, 7.30203E-10, 3.40515E-6 /) ! decay constants in 1/s

! Initialize second material
m_status = morty_init_material("decaylib.def")
m_status = morty_set_initial_mass_kg(2, ref_zaids2, ref_vals2)

m_status = morty_move_to_s(1800.0d0) ! decay for 0.5 hours
m_status = morty_move_to_s(3600.0d0) ! decay for 1 hour (total)

allocate(dec_conc(3))
m_status = morty_get_mass_kg(2, ref_z2, dec_conc)

Loading initial concentrations from an II.json file

The following example demonstrates how to load initial problem concentrations onto a material from an external II.json file. Here we search for a specific response ("I'm a pickle!") from the JSON file ("pickle_rick.json").

const std::string  libname = "decaylib.def";
STATUS initialized = morty_init_material(libname); // should be 1

std::string response_name = "I'm a pickle!";
std::string morty_json    = "pickle_rick.json";

// loaded == 0 => successful load
STATUS loaded = morty_load_initial_mass_from_file(
                1, response_name.c_str(), morty_json.c_str(), 0);

The equivalent method in FORTRAN would be:

integer(C_INT) :: m_status
real(C_DOUBLE), allocatable :: in_vals(:)
integer(C_INT), parameter :: ref_z(3) = (/ 92, 53, 55 /)
integer(C_INT), parameter :: ref_zaids(5) = (/ 92235, 92238, 53133, 55133, 8016 /)
real(C_DOUBLE), parameter :: ref_vals(3)  = (/ 10000.0d0, 9.672249d0, 1.000880d0 /) ! masses in kg
real(C_DOUBLE), parameter :: init_masses(5) = (/ 350.0d0, 9650.0d0, 10.0d0, 1.0d0, 1000.d0 /)
real(C_DOUBLE), parameter :: scaled_vals(3) = (/9.081827d2, 8.784170d-1, 9.089823d-2 /)
call morty_clear()

m_status = morty_init_material("decaylib.def")
m_status = morty_set_initial_mass_kg(1,ref_zaids, init_masses)

m_status = morty_move_to_s(3600.0d0)

m_status = morty_save_mass_to_file(1,"I'm a pickle!", "pickle_rick.json")

m_status = morty_init_material("decaylib.def")
m_status = morty_load_initial_mass_from_file(2, "I'm a pickle!", "pickle_rick.json", 0)

allocate(in_vals(3))
m_status = morty_get_mass_kg(2, ref_z, in_vals)

To load from a JSON normalized to a total mass of 1000 kg, one would simply call the morty_load_initial_mass_from_file function with an additional argument:

m_status = morty_init_material("decaylib.def")
m_status = morty_load_initial_mass_from_file(3, "I'm a pickle!", "pickle_rick.json", 0, 1000.d0)

Moving mass between two different materials

The following example illustrates moving an unstable nuclide from one material to another. The case starts with 1.0 kg of I-133 (unstable) and 1.0 kg of Ba-134 (stable) on material 1 (the source). The material is initially decayed for 20.80011 hours (one half-life of I-133) such that 50% of the original I-133 remains. After this, 0.3 kg of I is moved from material 1 to material 2; afterward, the two mixtures are decayed for an additional 20.8 hours (i.e., one half-life of I-133). Note that the separation occurs only after the timestep forward (i.e., separation occurs at 20.8 hours but doesn’t “register” until the next step forward in the simulation).

const std::string libname = "decaylib.def";

std::vector<int>    ids    = {53133, 56134}; // I-133, Ba-134
std::vector<int>    ele_z  = {53, 56};
std::vector<double> masses = {1.0, 1.0};

// Now initialize material, etc.
Origen::global_morty()->clear();
MAT_ID source_mat = morty_init_material(libname);

morty_set_initial_mass_kg(
    source_mat, ids.data(), ids.size(), masses.data(), masses.size());

// Initial decay, such that 50% of I-133 (0.5 kg) remains
const double decay_time_i133_hl = 20.80011 * 3600;

std::vector<double> final_mass(2, 0.0);

// Source mat should have { 0.5, 1.0 } kg of I / Ba
morty_get_mass_kg(source_mat,
                  ele_z.data(),
                  ele_z.size(),
                  final_mass.data(),
                  final_mass.size());



std::vector<double> target_masses = {0.0, 0.0};
// target_mat = 2 => created material #2
MAT_ID              target_mat    = morty_init_material(libname);

morty_set_initial_mass_kg(target_mat,
                          ids.data(),
                          ids.size(),
                          target_masses.data(),
                          target_masses.size());

morty_move_to_s(decay_time_i133_hl);

// Now move 0.3 kg from material 1 => 2 and decay for 20.8 h
// (one half-life)
std::vector<int>    elements_to_move = {53};
std::vector<double> masses_to_move   = {0.3};
morty_move_mass_kg(source_mat,
                   target_mat,
                   elements_to_move.data(),
                   elements_to_move.size(),
                   masses_to_move.data(),
                   masses_to_move.size());

morty_move_to_s(decay_time_i133_hl * 2);

// Separation only takes effect after we move forward in time
morty_get_mass_kg(source_mat,
                  ele_z.data(),
                  ele_z.size(),
                  actual_source_masses.data(),
                  actual_source_masses.size());
morty_get_mass_kg(target_mat,
                  ele_z.data(),
                  ele_z.size(),
                  actual_target_masses.data(),
                  actual_target_masses.size());

// Source material should now have { 0.1, 1.0 } kg I / Ba
// Target material should now have { 0.15, 0.0 } kg I / Ba

The equivalent FORTRAN call would be as-follows:

integer(C_INT) :: m_status
real(C_DOUBLE), allocatable:: source_masses(:), target_masses(:), mass_moved(:)
real(C_DOUBLE) :: decay_time, cumulative_time
integer(C_INT), allocatable :: zaids(:), ele_moved(:)
character(80,kind=C_CHAR) :: libname

integer(C_INT), parameter :: ref_ids(2) = (/ 53133, 56134 /) ! I-133, Ba-134
integer(C_INT), parameter :: ref_z(2) = (/ 53, 56 /)
real(C_DOUBLE), parameter :: ref_masses = (/ 1.0, 1.0 /)
real(C_DOUBLE), parameter :: decay_time_i_133_hl = 74880.396 ! 20.80011 h * 3600 s/h

call morty_clear()

allocate(source_masses(2), target_masses(2))

! Source material (material 1)
libname = 'decaylib.def'
m_status = morty_init_material(trim(libname))
m_status = morty_set_initial_mass_kg(1,ref_ids, ref_masses)

! Target material (material 2)
m_status = morty_init_material(trim(libname))
target_masses(:) = 0.0
m_status = morty_set_initial_mass_kg(2, ref_ids, target_masses)

! Initial decay such that 50% of I-133 (0.5 kg) remains
morty_move_to_s(decay_time_i133_hl)

! Now move 0.3 kg from material 1 to material 2 and decay for 20.8 h
! (i.e., one half-life)
allocate(ele_moved(1), mass_moved(1))
ele_moved(:) = 53
mass_moved(:) = 0.3

m_status = morty_move_mass_kg(1,2,ele_to_move,mass_moved)
m_status = morty_move_to_s(decay_time_i133_hl*2)

m_status = morty_get_mass_kg(1, ref_z, source_masses)
m_status = morty_get_mass_kg(2, ref_z, target_masses)

! source_masses should now have (0.10, 1.0) kg I / Ba
! target_masses should now have (0.15, 0.0) kg I / Ba