diff --git a/src/mam4xx/gas_chem.hpp b/src/mam4xx/gas_chem.hpp index d1f8242ce..40df3dfdd 100644 --- a/src/mam4xx/gas_chem.hpp +++ b/src/mam4xx/gas_chem.hpp @@ -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 @@ -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 // ----------------------------------------------------- @@ -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]); @@ -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]; @@ -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; @@ -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] @@ -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; @@ -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 @@ -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) { @@ -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 @@ -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 @@ -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 { diff --git a/src/mam4xx/gas_chem_mechanism.hpp b/src/mam4xx/gas_chem_mechanism.hpp index e6ed56d88..d1c979790 100644 --- a/src/mam4xx/gas_chem_mechanism.hpp +++ b/src/mam4xx/gas_chem_mechanism.hpp @@ -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]) { @@ -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) { @@ -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; @@ -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]; @@ -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]) { @@ -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;