Skip to content

Commit

Permalink
Merge pull request #229 from eagles-project/mjs/gas_chem_deepDive
Browse files Browse the repository at this point in the history 
Add Comments and Concerns to gas_chem, plus one minor fix
  • Loading branch information
@mjs271
mjs271 authored Aug 25, 2023
2 parents 30e2bed + 0b0ce15 commit 403291a
Show file tree
Hide file tree
Showing 2 changed files with 90 additions and 11 deletions.
70 changes: 59 additions & 11 deletions src/mam4xx/gas_chem.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,8 @@ namespace gas_chemistry {
// BAD CONSTANTs
const int itermax = 11;
const Real rel_err = 1.0e-3;
const Real high_rel_err = 1.0e-4;
// NOTE: high_rel_err is unused currently
// const Real high_rel_err = 1.0e-4;
const int max_time_steps = 1000;

KOKKOS_INLINE_FUNCTION
Expand Down Expand Up @@ -92,6 +93,27 @@ void newton_raphson_iter(const Real dti, const Real lin_jac[nzcnt],
// work array
Real epsilon[clscnt4]) {

// dti := 1 / dt
// lrxt := reaction rates in 1D array [1/cm^3/s]
// lhet := washout rates [1/s]
// iter_invariant := dti * solution + ind_prd
// factor := boolean controlling whether to do LU factorization
// lsol is local solution--appears to be identical to 'solution'
// looks like 'solution' is local to imp_sol() and 'lsol' is local to
// newton_raphson_iter(), and the final solutions is 'base_sol'
// these are volume mixing ratios [kmol species/kmol dry air]
// lsol := base_sol (at initialization), then when converged base_sol = lsol
// solution := array from imp_sol that holds the intermediate solutions and
// also holds the solution after converged [kmol species/kmol dry air]
// converged := array for entrywise convergence bools
// convergence := overall bool flag for convergence
// prod/loss := chemical production/loss rates [1/cm^3/s]
// NOTE: max_delta doesn't appear to be used for anything within gas_chem.hpp
// however, it looks like it's written to output in MAM4
// max_delta := abs(forcing / solution) if abs(solution) > 1.0e-20 and
// 0 otherwise
// epsilon := rel_err = 1.0e-3 (hardcoded above)

// -----------------------------------------------------
// the newton-raphson iteration for f(y) = 0
// -----------------------------------------------------
Expand Down Expand Up @@ -124,6 +146,9 @@ void newton_raphson_iter(const Real dti, const Real lin_jac[nzcnt],
imp_prod_loss(prod, loss, // out
lsol, lrxt, lhet); // in

// the units are internally consistent here, providing that
// iter_invariant, prod, loss all have units [1/s] to match up with
// solution (vmr) [-] and dti [1/s]. however, there could be other answers
for (int mm = 0; mm < clscnt4; ++mm) {
forcing[mm] =
solution[mm] * dti - (iter_invariant[mm] + prod[mm] - loss[mm]);
Expand All @@ -141,6 +166,9 @@ void newton_raphson_iter(const Real dti, const Real lin_jac[nzcnt],
// ... convergence measures
// -----------------------------------------------------------------------

// NOTE: is there a particular reason we don't check on the first iteration?
// seems like it'd be better to avoid the if on every iteration loop.
// same deal below
if (nr_iter > 0) {
for (int kk = 0; kk < clscnt4; ++kk) {
int mm = permute_4[kk];
Expand Down Expand Up @@ -184,13 +212,21 @@ void newton_raphson_iter(const Real dti, const Real lin_jac[nzcnt],
converged[kk] = true;

int mm = permute_4[kk];
// TODO: is there a computational reason this needs to happen?
// I suspect not, given that epsilon is hard-coded to 1e-3, meaning that
// all of this logic surrounding 'converged[kk] = ...' is unnecessary
bool frc_mask = haero::abs(forcing[mm]) > small;
if (frc_mask) {
// this ends up effectively being:
// if (small < abs(forcing) <= eps * abs(sol))
// => converged
// so the lower bound appears unnecessary
converged[kk] =
haero::abs(forcing[mm]) <= epsilon[kk] * haero::abs(solution[mm]);
} else {
// and this is just; if (abs(forcing) <= small <= eps) => converged
// and the implicit comparison of small and eps is not helpful
converged[kk] = true;

} // frc_mask
if (!converged[kk]) {
convergence = false;
Expand All @@ -200,9 +236,9 @@ void newton_raphson_iter(const Real dti, const Real lin_jac[nzcnt],
if (convergence) {
return;
}
} // end nr_iter
} // end nr_iter
} // newton_raphson_iter
} // end if (nr_iter > 0)
} // end nr_iter loop
} // newton_raphson_iter() function

KOKKOS_INLINE_FUNCTION
void imp_sol(Real base_sol[gas_pcnst], // inout - species mixing ratios [vmr]
Expand All @@ -221,6 +257,10 @@ void imp_sol(Real base_sol[gas_pcnst], // inout - species mixing ratios [vmr]
// this source is meant for small l1 cache machines such as
// the intel pentium and itanium cpus
// ---------------------------------------------------------------------------

// NOTE:
// extfrc := external in-situ forcing [1/cm^3/s]

const Real zero = 0;
const Real half = 0.5;
const Real one = 1;
Expand All @@ -239,6 +279,8 @@ void imp_sol(Real base_sol[gas_pcnst], // inout - species mixing ratios [vmr]
// -----------------------------------------------------------------------
// ... class independent forcing
// -----------------------------------------------------------------------
// FIXME: BAD CONSTANT
// what does this 4 represent, and would it ever be different?
indprd(4, // in
ind_prd, // inout
reaction_rates, extfrc); // in
Expand All @@ -253,6 +295,8 @@ void imp_sol(Real base_sol[gas_pcnst], // inout - species mixing ratios [vmr]
int cut_cnt = 0;
int fail_cnt = 0;
int stp_con_cnt = 0;
// track how much of the outer time step = delt (interval) has been completed
// during Newton-Raphson iteration
Real interval_done = zero;
// time_step_loop
for (int i = 0; i < max_time_steps; ++i) {
Expand All @@ -275,6 +319,13 @@ void imp_sol(Real base_sol[gas_pcnst], // inout - species mixing ratios [vmr]
// ... set the iteration invariant part of the function f(y)
// -----------------------------------------------------------------------

// TODO: the units seem wrong here--could these arrays hold quantities
// with different units?
// ind_prd has units [1/cm^3/s] (for the entries that are nonzero)
// dti units are [1/s], and
// solution is a volume mixing ratio [kmol species/kmol dry air]
// NOTE: this could be correct if solution had units [1/cm^3]
// which would line up with a number concentration
for (int mm = 0; mm < clscnt4; ++mm) {
iter_invariant[mm] = dti * solution[mm] + ind_prd[mm];
} // mm
Expand Down Expand Up @@ -307,19 +358,16 @@ void imp_sol(Real base_sol[gas_pcnst], // inout - species mixing ratios [vmr]
stp_con_cnt = 0;

if (cut_cnt < cut_limit) {
cut_cnt += cut_cnt;

cut_cnt += 1;
if (cut_cnt < cut_limit) {
dt *= half;

} else {
dt *= 0.1;
} // cut_cnt < cut_limit
// FIXME: is this part of the below FIXME?
// FIXME: figure out how we want to do error handling/logging
// break;
// cycle time_step_loop
} else {
// FIXME
// write(iulog,'('' imp_sol: Failed to converge @
// (lchnk,lev,col,nstep,dt,time) = '',4i6,1p,2e21.13)') &
// lchnk,lev,icol,nstep,dt,interval_done+dt
Expand All @@ -342,7 +390,7 @@ void imp_sol(Real base_sol[gas_pcnst], // inout - species mixing ratios [vmr]
if (haero::abs(delt - interval_done) <= 0.0001) {
if (fail_cnt > 0) {
// FIXME: probably handle this more gracefully via error logging?
printf("'imp_sol : @ (lchnk,lev,col) = \n");
printf("ERROR: imp_sol failure @ (lchnk,lev,col) = \n");
}
break;
} else {
Expand Down
31 changes: 31 additions & 0 deletions src/mam4xx/gas_chem_mechanism.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,14 @@ void adjrxt(Real rate[rxntot], Real inv[nfs], Real m) {
rate[1] *= inv[6] * inv[6] / m;
} // adjrxt

// TODO: unless rxt[0:6] and/or het_rates[1:4] have different units than the
// rest of the arrays the below additions seem fishy
// Units:
// rxt := reaction rates in 1D array [1/cm^3/s]
// het_rates := washout rates [1/s]
// TODO: the lines of concern *kind of* bear resemblance to the similarly
// concerning lines in linmat(), though it's difficult to tell if that results
// in consistent units
KOKKOS_INLINE_FUNCTION
void imp_prod_loss(Real prod[clscnt4], Real loss[clscnt4], Real y[gas_pcnst],
const Real rxt[rxntot], const Real het_rates[gas_pcnst]) {
Expand Down Expand Up @@ -121,7 +129,10 @@ void imp_prod_loss(Real prod[clscnt4], Real loss[clscnt4], Real y[gas_pcnst],
KOKKOS_INLINE_FUNCTION
void indprd(const int class_id, Real prod[clscnt4], const Real rxt[rxntot],
const Real extfrc[extcnt]) {
// extfrc := external in-situ forcing [1/cm^3/s]
// thus, prod must have units [1/cm^3/s]
const Real zero = 0;
// this is hard-coded to 4 outside of this function
if (class_id == 1) {
prod[0] = zero;
} else if (class_id == 4) {
Expand Down Expand Up @@ -157,6 +168,9 @@ void indprd(const int class_id, Real prod[clscnt4], const Real rxt[rxntot],
prod[29] = +extfrc[7];
} // indprd
}

// NOTE: at this point we are taking RHS units of (maybe) [1/s] to [s]
// and the units are internally consistent
KOKKOS_INLINE_FUNCTION
void lu_fac(Real lu[nzcnt]) {
const Real one = 1;
Expand Down Expand Up @@ -192,6 +206,9 @@ void lu_fac(Real lu[nzcnt]) {
lu[31] = one / lu[31];
} // lu_fac

// TODO: again, the units for b[1:2] look like they could be inconsistent
// lu = sys_jac [s]--mostly, maybe?; when passed in within gas_chem.hpp
// b = forcing [1/s]--maybe; when passed in from gas_chem.hpp
KOKKOS_INLINE_FUNCTION
void lu_slv(Real lu[nzcnt], Real b[clscnt4]) {
b[29] *= lu[31];
Expand Down Expand Up @@ -223,11 +240,21 @@ void lu_slv(Real lu[nzcnt], Real b[clscnt4]) {
b[3] *= lu[5];
b[2] -= lu[4] * b[3];
b[2] *= lu[3];
// the above 2 lines are equivalent to b[2] = (b[2] - (lu[4] * b[3])) * lu[3]
b[1] -= lu[2] * b[2];
b[1] *= lu[1];
// the above 2 lines are equivalent to b[1] = (b[1] - (lu[2] * b[2])) * lu[1]
b[0] *= lu[0];
} // lu_slv

// TODO: as in imp_prod_loss() unless rxt[0:6] and/or het_rates[1:4] have
// different units than the rest of the arrays the below additions seem suspect
// Units:
// rxt := reaction rates in 1D array [1/cm^3/s]
// het_rates := washout rates [1/s]
// TODO: the lines of concern *kind of* bear resemblance to the similarly
// concerning lines in imp_prod_loss(), though it's difficult to tell if that
// results in consistent units
KOKKOS_INLINE_FUNCTION
void linmat(Real mat[nzcnt], const Real rxt[rxntot],
const Real het_rates[gas_pcnst]) {
Expand Down Expand Up @@ -265,6 +292,10 @@ void linmat(Real mat[nzcnt], const Real rxt[rxntot],
mat[31] = -(+het_rates[30]);
} // linmat

// TODO: the below calculations appear to have inconsistent units for
// mat[0, 2, 3, 5]
// NOTE: it's *possible* we could be ok, if the odd-looking sums in linmat()
// turn out to be ok
KOKKOS_INLINE_FUNCTION
void nlnmat(Real mat[nzcnt], const Real lmat[nzcnt], const Real dti) {
mat[0] = lmat[0] - dti;
Expand Down

0 comments on commit 403291a

Please sign in to comment.