Functions
kernel_cram.hpp File Reference

Header for kernel_cram.cpp, kernel_cram_coefficients.cpp and kernel_cram_from_fortran.cpp. More...

#include "ScaleUtils/SuperLU/SRC/slu_dcomplex.h"
#include "Standard/Interface/Reporter.h"

Functions

int kernel_cram (double *n, double *n0, double **source_term, int source_order, double delta_t, int itot, int non, int *non0, int *kd, double *diag, double *a, 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 *a_n, double **p_a_val, int **p_a_ci, int **p_a_rp, int **p_a_diag_idx, int *a_nnz, int **p_newIdx, int **p_oldIdx, double **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)
 
int kernel_cram_from_fortran (double *n, double *n0, double *source_term, int source_order, double delta_t, int itot, int non, int *non0, int *kd, double *diag, double *a, int *loc, int zero_flux_step, int is_adjoint, int cram_order, int internal_steps, int remove_negatives, double cutoff, int *ZAI)
 
int kernel_cram_coefficients (doublecomplex *alpha0, doublecomplex *alpha, doublecomplex *theta, int cram_order)
 

Detailed Description

Header for kernel_cram.cpp, kernel_cram_coefficients.cpp and kernel_cram_from_fortran.cpp.

Author
Aarno Isotalo

Function Documentation

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_orderOrder of the source term polynomial
[in]delta_tLength of the time step in seconds
[in]itotThe number of nuclides in the system
[in]nonNumber 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_stepLogical integer, is the flux zero?
[in]is_adjointLogical integer, is this adjoint calculation?
[in]cram_orderOrder of the CRAM approximation (just use 16).
[in]internal_stepsNumber of internal substeps to use
[in]remove_negativesHow 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]cutoffConcentrations 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.

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.

References sortCRSnuclides().

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

int kernel_cram_from_fortran ( double *  n,
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 
)

Interface function for calling kernel_cram() from Fortran. Converts the source_term to the form expected by kernel_cram().

Note
All arguments are the same as for kernel_cram(), except that source_term is a 1D array whose elements must be ordered so that the coefficients for given nuclide are consequtive in memory. This only converts source_term, loc indexing must be converted by the caller.
Returns
EXIT_SUCCESS or EXIT_FAILURE. On failed exit the results (i.e., n) are unspecified.

References kernel_cram().

int kernel_cram_coefficients ( doublecomplex *  alpha0,
doublecomplex *  alpha,
doublecomplex *  theta,
int  cram_order 
)

Set CRAM coefficients for partial fraction decomposition form.

Sets the poles and residues of chebyshev rational approximation on the negative real axis. Available orders are 4, 6, 8, 10, 12, 14, and 16.

Order 14 and 16 coefficients are from: Maria Pusa, "Correction to partial fraction decomposition coefficients for Chebyshev rational approximation on the negative real axis." arXiv:1206.2880v1 [math.NA] (2013). Others have been provided by Maria Pusa directly.

Parameters
[out]alpha0,alpha,thetaThe coefficients. alpha and theta must be preallocated to hold at least order/2 values
[in]cram_orderThe desired order of the CRAM approximation.
Returns
EXIT_SUCCESS or EXIT_FAILURE. On failed exit the results are unspecified.

References NAME.

Referenced by kernel_cram().

Generated by   doxygen 1.8.10