CRAM depletion solver. More...

#include "Origen/Solver/cram/src/kernel_cram.hpp"
#include <stdio.h>
#include <cmath>
#include <ctime>
#include <iostream>
#include <limits>
#include <vector>
#include "ScaleUtils/SuperLU/SRC/slu_zdefs.h"
#include "Standard/Interface/jdebug.h"

Classes
struct	LUstore

Macros
#define	SORT_NUCLIDES 1

#define	MAKE_EXPM_INPUT 0

#define	MAKE_MATLAB_INPUT 0

#define	NAME "Origen/Manager/Wrapper/kernel_cram.cpp"
	Identifier for error messages. More...

Functions
void	CopyLU (SuperMatrix L, SuperMatrix U, LUstore *LU)

void	makeExpmInput (int itot, int a_n, int nnz, doublecomplex a_val_complex, int a_ci, int rp, double delta_t, doublecomplex Bvals, double *source_term, int source_order, int ZAI)

void	makeMatlabInput (int itot, int a_n, doublecomplex a_val_complex, int a_ci, int a_rp, int ZAI)

void	sortCRSnuclides (int a_n, int a_nnz, double a_val, int a_ci, int *a_rp, int newIdx, int oldIdx, int sort_nuclides, int ZAI)

int	kernel_cram (double n_out, double n0, double *source_term, int source_order, double delta_t, int itot, int non, int non0, int kd, double diag, double offdiag, int loc, int zero_flux_step, int is_adjoint, int cram_order, int internal_steps, int remove_negatives, double cutoff, int *ZAI)

void	OrigenToCRSMatrix (int p_a_n, double p_a_val, int p_a_ci, int p_a_rp, int p_a_diag_idx, int p_a_nnz, int p_newIdx, int p_oldIdx, double *p_src_dummy_n0, int itot, int non, int non0, int kd, double diag, double offdiag, int loc, double *source_term, int source_order, int zero_flux_step, int is_adjoint, int sort_nuclides, int ZAI)

Detailed Description

CRAM depletion solver.

Author: Aarno Isotalo

Macro Definition Documentation

#define SORT_NUCLIDES 1

Sort nuclides if ZAI is provided (logical int). Sorting improves performance by factor of 2. Unless problems arise this should always be on.

Referenced by kernel_cram().

#define MAKE_EXPM_INPUT 0

Controls printing of input for an external test program. For internal testing purposes only. 0 = not printed, 1 = print with separate source, 2 = print with source merged to the matrix with dummy elements

Referenced by makeExpmInput().

#define MAKE_MATLAB_INPUT 0

Controls printing of the matrix in matlab form (logical int). For internal testing purposes only.

Referenced by makeMatlabInput().

#define NAME "Origen/Manager/Wrapper/kernel_cram.cpp"

Identifier for error messages.

Referenced by kernel_cram().

Function Documentation

void CopyLU	(	SuperMatrix *	L,
		SuperMatrix *	U,
		LUstore *	LU
	)

Special fun fun copy the L and U matrices in SuperLU format.

This copies the pointers, not their contents. New Store is allocated. nzval (the non-zero values) pointers are also copied but then the nzval pointers in L U are set to point to new, freshly allocated memory. This way the data from L, U is now in LU without copying any memory. nzval in L, U are effectively reset, but they are not needed anyway (L,U are used to store the next LU decomposition for different pole in CRAM so nzval would be overwritten anyway.

Because the pointers, rather than the memory they point to is copied, all copies (except nzval) end up pointing to the same space, i.e., the "work" array in the solver. This works because Every A and thus every L, U has the same sparsity pattern.

The first pair of nzval (for term k=0) will point to work while the rest will point to the new arrays allocated on the previous k in this function.

Parameters

References LUstore::L, and LUstore::U.

Referenced by kernel_cram().

void makeExpmInput	(	int	itot,
		int	a_n,
		int	a_nnz,
		doublecomplex *	a_val_complex,
		int *	a_ci,
		int *	a_rp,
		double	delta_t,
		doublecomplex *	Bvals,
		double **	source_term,
		int	source_order,
		int *	ZAI
	)

Write the CRS matrix to a format used by an external test program "expm". This function exists for development purposes only. The input is written to file "expmInput" in current directory.

References MAKE_EXPM_INPUT, and r.

Referenced by kernel_cram().

void makeMatlabInput	(	int	itot,
		int	a_n,
		doublecomplex *	a_val_complex,
		int *	a_ci,
		int *	a_rp,
		int *	ZAI
	)

Write the matrix and ZAI to a .m file

References MAKE_MATLAB_INPUT, and r.

Referenced by kernel_cram().

void sortCRSnuclides	(	int	a_n,
		int	a_nnz,
		double **	p_a_val,
		int **	p_a_ci,
		int **	p_a_rp,
		int *	newIdx,
		int *	oldIdx,
		int	sort_nuclides,
		int *	ZAI
	)

Sorts nuclides in a CRS matrix decay/transmutation system by ascending ZAI.

Input is a CRS matrix and a ZAI array. ZAI can be nullptr in which case sorting is not performed. The oldIdx and newIdx arrays are generated either way (they are trivial if no sorting was done).

Parameters

[in]	a_n	Size of the CRS matrix A
[in]	a_nnz	Number of elements in CRS matrix A
[in,out]	p_a_val[a_nnz]	Non-zero elements of CRS matrix A
[in,out]	p_a_ci[a_nnz]	Column indices of CRS matrix A
[in,out]	p_a_rp[itot+1]	Row pointers of CRS matrix A
[out]	oldIdx[itot]	oldIdx[i] is the old index of the nuclide that has new index i
[out]	newIdx[itot]	newIdx[i] is the new index of the nuclide that had old index i
[in]	sort_nuclides	Logical, whether to actually sort the nuclides
[in]	ZAI[itot]	ZAI identifiers of the nuclides. Can be nullptr in which case sorting is not actually done. (Still needs to be called to set the index arrays.)

Referenced by OrigenToCRSMatrix().

int kernel_cram	(	double *	n_out,
		double *	n0,
		double **	source_term,
		int	source_order,
		double	delta_t,
		int	itot,
		int	non,
		int *	non0,
		int *	kd,
		double *	diag,
		double *	offdiag,
		int *	loc,
		int	zero_flux_step,
		int	is_adjoint,
		int	cram_order,
		int	internal_steps,
		int	remove_negatives,
		double	cutoff,
		int *	ZAI
	)

CRAM depletion solver with source term and adjoint capabilities.

Based on the publication: Pusa, M. and Leppänen, J. "Computing the matrix exponential in burnup calculations." Nucl. Sci. Eng., 164 (2010) 140-150.

The required linear system solutions are done with the SuperLU library: Xiaoye S. Li. "An Overview of {SuperLU}: Algorithms, Implementation, and User Interface." ACM Trans.Math. Softw., 31;3 (2005) 302-325.

Fortran should call kernel_cram_from_fortran() instead.

Source term can have polynomial time dependency. Due to numerical issues source order is limited to cram_order/2-2. Non constant source term will cause some deterioration in accuracy, but below this limit it should be negligible.

Adjoint calculation is also possible by passing is_adjoint=1. In adjoint mode n0, source_term and delta_t>0 must be those for the adjoint problem while everything else must have values defining the forward system.

Internal substeps: Dividin steps to pieces greatly improved accuracy. This is particularly significant for decay problems where accuracy without might not be good enough for all purposes. The solver can do this division internally, which is more efficient than calling the solver multiple times as the LU decompositions only need to be formed once. internal_steps=1 means means that there is no substepping. Setting internal_steps=2 increases running time by roughly 25%. Each substep above 2 increases it by another 10-15% of the running time without substepping. See *** for details.

Parameters

[out]	n_out[itot]	Final concentrations. Must be Pre-allocated.
[in]	n0[itot]	Initial concentrations
[in]	source_term[itot][source_order+1]	Polynomial coefficients s_i,j for the source term
[in]	source_order	Order of the source term polynomial
[in]	delta_t	Length of the time step in seconds
[in]	itot	The number of nuclides in the system
[in]	non	Number of reactions (inc. decay) in the system
[in]	non0[itot]	Number of parents for each nuclide
[in]	kd[itot]	Number of decay parents for each nuclide
[in]	diag[itot]	Diagonal elements of A, i.e., (-lambda_eff)
[in]	offdiag[non]	Reaction specific decay constant of each reaction.
[in]	loc[non]	Index of the parent in each reaction.
[in]	zero_flux_step	Logical integer, is the flux zero?
[in]	is_adjoint	Logical integer, is this adjoint calculation?
[in]	cram_order	Order of the CRAM approximation (just use 16).
[in]	internal_steps	Number of internal substeps to use
[in]	remove_negatives	How to treat negative vaues in results: 0: Do nothing. 1: Remove all negative values. 2: Remove all negative values in forward mode. Remove all negative values in adjoint mode if source term is constant and both source and initial concentrations are non-negative.
[in]	cutoff	Concentrations whose ABSOLUTE value is below cutoff times sum(abs(n_out)) are set to zero. This is applied after remove_negatives.
[in]	ZAI[itot]	ZAI identifiers of the nuclides.

Note

The indexing in loc must be in c style, i.e., starts from 0. Be avare of this when calling from Fortran.

offdiag (which is called "a" in most of Origen source code) and loc must be sorted in the special Origen way:

decay parents of first nuclide,
activation parents of first nuclide,
decay parents of second nuclide,
activation parents of the second,
and so on.

The culling feature from earlier versions was removed. Use internal substeps instead.

Returns: zero on successfull exit or nonzero on failed exit. On failed exit the results (i.e., n) are unspecified.

Examples:: tstCRAM.cpp.

References LUstore::C, CopyLU(), LUstore::equed, Origen::info(), kernel_cram_coefficients(), makeExpmInput(), makeMatlabInput(), NAME, OrigenToCRSMatrix(), r, LUstore::R, and SORT_NUCLIDES.

Referenced by kernel_cram_from_fortran(), Solver_cram::solve_impl(), TEST(), and TEST_F().

void OrigenToCRSMatrix	(	int *	p_a_n,
		double **	p_a_val,
		int **	p_a_ci,
		int **	p_a_rp,
		int **	p_a_diag_idx,
		int *	p_a_nnz,
		int **	p_newIdx,
		int **	p_oldIdx,
		double **	p_src_dummy_n0,
		int	itot,
		int	non,
		int *	non0,
		int *	kd,
		double *	diag,
		double *	offdiag,
		int *	loc,
		double **	source_term,
		int	source_order,
		int	zero_flux_step,
		int	is_adjoint,
		int	sort_nuclides,
		int *	ZAI
	)

Convert ORIGEN form input to regular compressed row storage (CRS) matrix representing the equivalent homogeneous equations.

This has two main functions:

Build the actual coefficient matrix
Homogenize the equations (i.e., dn/dt = An +s ==> dn'/dt = A'n')

which are done simultaneously. In addition this:

Allocates all diagonal elements (even zeros)
Sorts nuclides by ZAI
Sorts each row

For details about the homogenization see:

Note: The first nine arguments are used for output. They are pointers to the variables p_<name in kernel_cram()> The other arguments are same as for kernel_cram(). The output, except for p_a_n and p_a_nnz which are scalars, will be allocated by this function.

Examples:: tstCRAM.cpp.

References sortCRSnuclides().

Referenced by kernel_cram(), TEST(), and tstOrigenToCRSMatrix::~tstOrigenToCRSMatrix().

Classes

Macros

Functions

Detailed Description

Macro Definition Documentation

Function Documentation