From 21b1876cf9b1f563de52530b57a17755d34e02f6 Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Tue, 19 Aug 2025 09:50:36 -0500 Subject: [PATCH 1/4] basic implementation without fixes thermo consts --- src/core/impl/TChem_Impl_MOSAIC.hpp | 755 ++++++++++++++++++++++++++++ 1 file changed, 755 insertions(+) diff --git a/src/core/impl/TChem_Impl_MOSAIC.hpp b/src/core/impl/TChem_Impl_MOSAIC.hpp index 98e04a3d..b5e764d7 100644 --- a/src/core/impl/TChem_Impl_MOSAIC.hpp +++ b/src/core/impl/TChem_Impl_MOSAIC.hpp @@ -2421,6 +2421,761 @@ struct MOSAIC{ aer_liquid(mosaic.ilim2_a) = 0.0; } // adjust_liquid_aerosol + KOKKOS_INLINE_FUNCTION static + void compute_activities(const MosaicModelData& mosaic, + const real_type_1d_view_type& aer_liquid, + const real_type_1d_view_type& ma, + const real_type_1d_view_type& mc, + const real_type_1d_view_type& electrolyte_solid, + const real_type_1d_view_type& electrolyte_liquid, + const real_type_1d_view_type& electrolyte_total, + const real_type_1d_view_type& log_gam, + const real_type_2d_view_type& log_gamZ, + const real_type_1d_view_type& gam, + const real_type_1d_view_type& activity, + real_type& T_K) { + + // local variables + real_type_1d_view_type xmol("xmol", mosaic.nelectrolyte); + real_type a_c, Keq_ll, sum_elec, gam_ratio, water_a, mSULF, gam_ratio; + ordinal_type jA; + // FIXME: these need to be set beforehand + ordinal_type jaerosolstate, jphase, jhyst_leg; + + // get aerosol water activity + aerosol_water(mosaic, electrolyte_liquid, water_a, jaerosolstate, jphase, jhyst_leg); + + if (water_a == 0.0) { + return; + } + + // get sulfate ratio to determine regime + real_type XT = 0.0; + calcuate_XT(aer_liquid, mosaic, XT); + + if (XT > 2.0 || XT < 0.0) { + // SULFATE POOR: fully dissociated electrolytes + + // anion molalities (mol / kg water) + ma(mosaic.ja_so4) = 1.e-9 * aer_liquid(mosaic.iso4_a) / water_a; + ma(mosaic.ja_hso4) = 0.0; + ma(mosaic.ja_no3) = 1.e-9 * aer_liquid(mosaic.ino3_a) / water_a; + ma(mosaic.ja_cl) = 1.e-9 * aer_liquid(mosaic.icl_a) / water_a; + ma(mosaic.ja_msa) = 1.e-9 * aer_liquid(mosaic.imsa_a) / water_a; + + // cation molalities (mol / kg water) + mc(mosaic.jc_ca) = 1.e-9 * aer_liquid(mosaic.ica_a) / water_a; + mc(mosaic.jc_nh4) = 1.e-9 * aer_liquid(mosaic.inh4_a) / water_a; + mc(mosaic.jc_na) = 1.e-9 * aer_liquid(mosaic.na_a) / water_a; + a_c = ( + (2. * ma(mosaic.ja_so4) + + ma(mosaic.ja_no3) + + ma(mosaic.ja_cl) + + ma(mosaic.ja_msa)) - + (2. * mc(mosaic.jc_ca) + + mc(mosaic.jc_nh4) + + mc(mosaic.jc_na)) ); + // FIXME: consider adding update_thermodynamic_constants instead + // of computing equil. constants directly. + // Make sure temp. is handled + fn_Keq(mosaic.Keq_ll_298(2), mosaic.Keq_a_ll(2), mosiac.Keq_b_ll(2), T_K, Keq_ll); + mc(mosaic.jc_h) = 0.5 * ( (a_c) + + (ats::sqrt(a_c*a_c + 4. * Keq_ll)) ); + + if (mc(mosaic.jc_h) == 0.0) { + mc(mosaic.jc_h) = 1.e-10; + } + + sum_elec = 2. * electrolyte_liquid(mosaic.jnh4no3) + + 2. * electrolyte_liquid(mosaic.jnh4cl) + + 3. * electrolyte_liquid(mosaic.jnh4so4) + + 3. * electrolyte_liquid(mosaic.jna2so4) + + 2. * electrolyte_liquid(mosaic.jnano3) + + 2. * electrolyte_liquid(mosaic.jnacl) + + 3. * electrolyte_liquid(mosaic.jcano3) + + 3. * electrolyte_liquid(mosaic.jcacl2) + + 2. * electrolyte_liquid(mosaic.jhno3) + + 2. * electrolyte_liquid(mosaic.jhcl); + + if (sum_elec == 0.0) { + for (ordinal_type jA = 0; jA < mosaic.nelectrolyte; jA++) { + gam(jA) = 1.0; + } + // FIXME: duplicated code to avoid goto statement + gam(mosaic.jlvcite) = 1.0; + gam(mosaic.jnh4hso4) = 1.0; + gam(mosaic.jnh4msa) = 1.0; + gam(mosaic.jna3hso4) = 1.0; + gam(mosaic.jnamsa) = 1.0; + gam(mosaic.jcamsa2) = 1.0; + + activity(mosaic.jlvcite) = 0.0; + activity(mosaic.jnh4hso4) = 0.0; + + activity(mosaic.jnh4msa) = mc(mosaic.jc_nh4) * ma(mosaic.ja_msa) * + gam(mosaic.jnh4msa) * gam(mosaic.jnh4msa); + + activity(mosaic.jna3hso4) = 0.0; + activity(mosaic.jnahso4) = 0.0; + + activity(mosaic.jnamsa) = mc(mosaic.jc_na) * ma(mosaic.ja_msa) * + gam(mosaic.jnamsa) * gam(mosaic.jnamsa); + activity(mosaic.jcamsa2) = mc(mosaic.jc_ca) * ma(mosaic.ja_msa) * + gam(mosaic.jcamsa2) * gam(mosaic.jcamsa2) * gam(mosaic.jcamsa2); + + gam_ratio = gam(mosaic.jnh4no3) * gam(mosaic.jnh4no3) / + gam(mosaic.jhno3) * gam(mosaic.jhno3); + + return; + } + + // ionic mole fractions + xmol(mosaic.jnh4no3) = 2. * electrolyte_liquid(mosaic.jnh4no3) / sum_elec; + xmol(mosaic.jnh4cl) = 2. * electrolyte_liquid(mosaic.jnh4cl) / sum_elec; + xmol(mosaic.jnh4so4) = 3. * electrolyte_liquid(mosaic.jnh4so4) / sum_elec; + xmol(mosaic.jna2so4) = 3. * electrolyte_liquid(mosaic.jna2so4) / sum_elec; + xmol(mosaic.jnano3) = 2. * electrolyte_liquid(mosaic.jnano3) / sum_elec; + xmol(mosaic.jnacl) = 2. * electrolyte_liquid(mosaic.jnacl) / sum_elec; + xmol(mosaic.jcano3) = 3. * electrolyte_liquid(mosaic.jcano3) / sum_elec; + xmol(mosaic.jcacl2) = 3. * electrolyte_liquid(mosaic.jcacl2) / sum_elec; + xmol(mosaic.jhno3) = 2. * electrolyte_liquid(mosaic.jhno3) / sum_elec; + xmol(mosaic.jhcl) = 2. * electrolyte_liquid(mosaic.jhcl) / sum_elec; + + jA = mosaic.nh4so4; + if (xmol(jA) > 0.0) { + log_gam(jA) = xmol(mosaic.jnh4no3) * log_gamZ(jA,mosaic.jnh4no3) + + xmol(mosaic.jnh4cl) * log_gamZ(jA,mosaic.jnh4cl) + + xmol(mosaic.jnh4so4) * log_gamZ(jA,mosaic.jnh4so4) + + xmol(mosaic.jna2so4) * log_gamZ(jA,mosaic.jna2so4) + + xmol(mosaic.jnano3) * log_gamZ(jA,mosaic.jnano3) + + xmol(mosaic.jnacl) * log_gamZ(jA,mosaic.jnacl) + + xmol(mosaic.jcano3) * log_gamZ(jA,mosaic.jcano3) + + xmol(mosaic.jcacl2) * log_gamZ(jA,mosaic.jcacl2) + + xmol(mosaic.jhno3) * log_gamZ(jA,mosaic.jhno3) + + xmol(mosaic.jhcl) * log_gamZ(jA,mosaic.jhcl); + gam(jA) = ats::pow(10., log_gam(jA)); + activity(jA) = ats::pow(mc(mosaic.jc_nh4),2.) * + ma(mosaic.ja_so4) * + ats::pow(gam(jA,3.)); + } + + jA = mosaic.jnh4no3; + if (xmol(jA) > 0.0) { + log_gam(jA) = xmol(mosaic.jnh4no3) * log_gamZ(jA,mosaic.jnh4no3) + + xmol(mosaic.jnh4cl) * log_gamZ(jA,mosaic.jnh4cl) + + xmol(mosaic.jnh4so4) * log_gamZ(jA,mosaic.jnh4so4) + + xmol(mosaic.jna2so4) * log_gamZ(jA,mosaic.jna2so4) + + xmol(mosaic.jnano3) * log_gamZ(jA,mosaic.jnano3) + + xmol(mosaic.jnacl) * log_gamZ(jA,mosaic.jnacl) + + xmol(mosaic.jcano3) * log_gamZ(jA,mosaic.jcano3) + + xmol(mosaic.jcacl2) * log_gamZ(jA,mosaic.jcacl2) + + xmol(mosaic.jhno3) * log_gamZ(jA,mosaic.jhno3) + + xmol(mosaic.jhcl) * log_gamZ(jA,mosaic.jhcl); + gam(jA) = ats::pow(10., log_gam(jA)); + activity(jA) = mc(mosaic.jc_nh4) * + ma(mosaic.ja_no3) * + ats::pow(gam(jA,2.)); + } + + jA = mosaic.jnh4cl; + if (xmol(jA) > 0.0) { + log_gam(jA) = xmol(mosaic.jnh4no3) * log_gamZ(jA,mosaic.jnh4no3) + + xmol(mosaic.jnh4cl) * log_gamZ(jA,mosaic.jnh4cl) + + xmol(mosaic.jnh4so4) * log_gamZ(jA,mosaic.jnh4so4) + + xmol(mosaic.jna2so4) * log_gamZ(jA,mosaic.jna2so4) + + xmol(mosaic.jnano3) * log_gamZ(jA,mosaic.jnano3) + + xmol(mosaic.jnacl) * log_gamZ(jA,mosaic.jnacl) + + xmol(mosaic.jcano3) * log_gamZ(jA,mosaic.jcano3) + + xmol(mosaic.jcacl2) * log_gamZ(jA,mosaic.jcacl2) + + xmol(mosaic.jhno3) * log_gamZ(jA,mosaic.jhno3) + + xmol(mosaic.jhcl) * log_gamZ(jA,mosaic.jhcl); + gam(jA) = ats::pow(10., log_gam(jA)); + activity(jA) = mc(mosaic.jc_nh4) * + ma(mosaic.ja_cl) * + ats::pow(gam(jA,2.)); + } + + jA = mosaic.jna2so4; + if (xmol(jA) > 0.0) { + log_gam(jA) = xmol(mosaic.jnh4no3) * log_gamZ(jA,mosaic.jnh4no3) + + xmol(mosaic.jnh4cl) * log_gamZ(jA,mosaic.jnh4cl) + + xmol(mosaic.jnh4so4) * log_gamZ(jA,mosaic.jnh4so4) + + xmol(mosaic.jna2so4) * log_gamZ(jA,mosaic.jna2so4) + + xmol(mosaic.jnano3) * log_gamZ(jA,mosaic.jnano3) + + xmol(mosaic.jnacl) * log_gamZ(jA,mosaic.jnacl) + + xmol(mosaic.jcano3) * log_gamZ(jA,mosaic.jcano3) + + xmol(mosaic.jcacl2) * log_gamZ(jA,mosaic.jcacl2) + + xmol(mosaic.jhno3) * log_gamZ(jA,mosaic.jhno3) + + xmol(mosaic.jhcl) * log_gamZ(jA,mosaic.jhcl); + gam(jA) = ats::pow(10., log_gam(jA)); + activity(jA) = ats::pow(mc(mosaic.jc_na),2.) * + ma(mosaic.ja_so4) * + ats::pow(gam(jA,3.)); + } + + jA = mosaic.jnano3; + if (xmol(jA) > 0.0) { + log_gam(jA) = xmol(mosaic.jnh4no3) * log_gamZ(jA,mosaic.jnh4no3) + + xmol(mosaic.jnh4cl) * log_gamZ(jA,mosaic.jnh4cl) + + xmol(mosaic.jnh4so4) * log_gamZ(jA,mosaic.jnh4so4) + + xmol(mosaic.jna2so4) * log_gamZ(jA,mosaic.jna2so4) + + xmol(mosaic.jnano3) * log_gamZ(jA,mosaic.jnano3) + + xmol(mosaic.jnacl) * log_gamZ(jA,mosaic.jnacl) + + xmol(mosaic.jcano3) * log_gamZ(jA,mosaic.jcano3) + + xmol(mosaic.jcacl2) * log_gamZ(jA,mosaic.jcacl2) + + xmol(mosaic.jhno3) * log_gamZ(jA,mosaic.jhno3) + + xmol(mosaic.jhcl) * log_gamZ(jA,mosaic.jhcl); + gam(jA) = ats::pow(10., log_gam(jA)); + activity(jA) = mc(mosaic.jc_na) * + ma(mosaic.ja_no3) * + ats::pow(gam(jA,2.)); + } + + jA = mosaic.jnacl; + if (xmol(jA) > 0.0) { + log_gam(jA) = xmol(mosaic.jnh4no3) * log_gamZ(jA,mosaic.jnh4no3) + + xmol(mosaic.jnh4cl) * log_gamZ(jA,mosaic.jnh4cl) + + xmol(mosaic.jnh4so4) * log_gamZ(jA,mosaic.jnh4so4) + + xmol(mosaic.jna2so4) * log_gamZ(jA,mosaic.jna2so4) + + xmol(mosaic.jnano3) * log_gamZ(jA,mosaic.jnano3) + + xmol(mosaic.jnacl) * log_gamZ(jA,mosaic.jnacl) + + xmol(mosaic.jcano3) * log_gamZ(jA,mosaic.jcano3) + + xmol(mosaic.jcacl2) * log_gamZ(jA,mosaic.jcacl2) + + xmol(mosaic.jhno3) * log_gamZ(jA,mosaic.jhno3) + + xmol(mosaic.jhcl) * log_gamZ(jA,mosaic.jhcl); + gam(jA) = ats::pow(10., log_gam(jA)); + activity(jA) = mc(mosaic.jc_na) * + ma(mosaic.ja_cl) * + ats::pow(gam(jA,2.)); + } + + // Note: these are commented out in MOSAIC also. + // jA = mosaic.jcano3; + // if (xmol(jA) > 0.) { + // gam(jA) = 1.; + // activity(jA) = 1.; + // } + + // jA = mosaic.jacl2; + // if (xmol(jA) > 0.) { + // gam(jA) = 1.; + // activity(jA) = 1.; + // } + + jA = mosaic.jcano3; + if (xmol(jA) > 0.0) { + log_gam(jA) = xmol(mosaic.jnh4no3) * log_gamZ(jA,mosaic.jnh4no3) + + xmol(mosaic.jnh4cl) * log_gamZ(jA,mosaic.jnh4cl) + + xmol(mosaic.jnh4so4) * log_gamZ(jA,mosaic.jnh4so4) + + xmol(mosaic.jna2so4) * log_gamZ(jA,mosaic.jna2so4) + + xmol(mosaic.jnano3) * log_gamZ(jA,mosaic.jnano3) + + xmol(mosaic.jnacl) * log_gamZ(jA,mosaic.jnacl) + + xmol(mosaic.jcano3) * log_gamZ(jA,mosaic.jcano3) + + xmol(mosaic.jcacl2) * log_gamZ(jA,mosaic.jcacl2) + + xmol(mosaic.jhno3) * log_gamZ(jA,mosaic.jhno3) + + xmol(mosaic.jhcl) * log_gamZ(jA,mosaic.jhcl); + gam(jA) = ats::pow(10., log_gam(jA)); + activity(jA) = mc(mosaic.jc_ca) * + ats::pow(ma(mosaic.ja_no3),2.) * + ats::pow(gam(jA,3.)); + } + + jA = mosaic.jcacl2; + if (xmol(jA) > 0.0) { + log_gam(jA) = xmol(mosaic.jnh4no3) * log_gamZ(jA,mosaic.jnh4no3) + + xmol(mosaic.jnh4cl) * log_gamZ(jA,mosaic.jnh4cl) + + xmol(mosaic.jnh4so4) * log_gamZ(jA,mosaic.jnh4so4) + + xmol(mosaic.jna2so4) * log_gamZ(jA,mosaic.jna2so4) + + xmol(mosaic.jnano3) * log_gamZ(jA,mosaic.jnano3) + + xmol(mosaic.jnacl) * log_gamZ(jA,mosaic.jnacl) + + xmol(mosaic.jcano3) * log_gamZ(jA,mosaic.jcano3) + + xmol(mosaic.jcacl2) * log_gamZ(jA,mosaic.jcacl2) + + xmol(mosaic.jhno3) * log_gamZ(jA,mosaic.jhno3) + + xmol(mosaic.jhcl) * log_gamZ(jA,mosaic.jhcl); + gam(jA) = ats::pow(10., log_gam(jA)); + activity(jA) = mc(mosaic.jc_ca) * + ats::pow(ma(mosaic.ja_cl),2.) * + ats::pow(gam(jA,3.)); + } + + jA = mosaic.jhno3; + if (xmol(jA) > 0.0) { + log_gam(jA) = xmol(mosaic.jnh4no3) * log_gamZ(jA,mosaic.jnh4no3) + + xmol(mosaic.jnh4cl) * log_gamZ(jA,mosaic.jnh4cl) + + xmol(mosaic.jnh4so4) * log_gamZ(jA,mosaic.jnh4so4) + + xmol(mosaic.jna2so4) * log_gamZ(jA,mosaic.jna2so4) + + xmol(mosaic.jnano3) * log_gamZ(jA,mosaic.jnano3) + + xmol(mosaic.jnacl) * log_gamZ(jA,mosaic.jnacl) + + xmol(mosaic.jcano3) * log_gamZ(jA,mosaic.jcano3) + + xmol(mosaic.jcacl2) * log_gamZ(jA,mosaic.jcacl2) + + xmol(mosaic.jhno3) * log_gamZ(jA,mosaic.jhno3) + + xmol(mosaic.jhcl) * log_gamZ(jA,mosaic.jhcl); + gam(jA) = ats::pow(10., log_gam(jA)); + activity(jA) = mc(mosaic.jc_h) * + ma(mosaic.ja_no3) * + ats::pow(gam(jA,2.)); + } + + jA = mosaic.jhcl; + if (xmol(jA) > 0.0) { + log_gam(jA) = xmol(mosaic.jnh4no3) * log_gamZ(jA,mosaic.jnh4no3) + + xmol(mosaic.jnh4cl) * log_gamZ(jA,mosaic.jnh4cl) + + xmol(mosaic.jnh4so4) * log_gamZ(jA,mosaic.jnh4so4) + + xmol(mosaic.jna2so4) * log_gamZ(jA,mosaic.jna2so4) + + xmol(mosaic.jnano3) * log_gamZ(jA,mosaic.jnano3) + + xmol(mosaic.jnacl) * log_gamZ(jA,mosaic.jnacl) + + xmol(mosaic.jcano3) * log_gamZ(jA,mosaic.jcano3) + + xmol(mosaic.jcacl2) * log_gamZ(jA,mosaic.jcacl2) + + xmol(mosaic.jhno3) * log_gamZ(jA,mosaic.jhno3) + + xmol(mosaic.jhcl) * log_gamZ(jA,mosaic.jhcl); + gam(jA) = ats::pow(10., log_gam(jA)); + activity(jA) = mc(mosaic.jc_h) * + ma(mosaic.ja_cl) * + ats::pow(gam(jA,2.)); + } + + // FIXME: duplicated code to avoid goto statement + gam(mosaic.jlvcite) = 1.0; + gam(mosaic.jnh4hso4) = 1.0; + gam(mosaic.jnh4msa) = 1.0; + gam(mosaic.jna3hso4) = 1.0; + gam(mosaic.jnamsa) = 1.0; + gam(mosaic.jcamsa2) = 1.0; + + activity(mosaic.jlvcite) = 0.0; + activity(mosaic.jnh4hso4) = 0.0; + + activity(mosaic.jnh4msa) = mc(mosaic.jc_nh4) * ma(mosaic.ja_msa) * + gam(mosaic.jnh4msa) * gam(mosaic.jnh4msa); + + activity(mosaic.jna3hso4) = 0.0; + activity(mosaic.jnahso4) = 0.0; + + activity(mosaic.jnamsa) = mc(mosaic.jc_na) * ma(mosaic.ja_msa) * + gam(mosaic.jnamsa) * gam(mosaic.jnamsa); + activity(mosaic.jcamsa2) = mc(mosaic.jc_ca) * ma(mosaic.ja_msa) * + gam(mosaic.jcamsa2) * gam(mosaic.jcamsa2) * gam(mosaic.jcamsa2); + + gam_ratio = gam(mosaic.jnh4no3) * gam(mosaic.jnh4no3) / + gam(mosaic.jhno3) * gam(mosaic.jhno3); + + // end SULFATE POOR regime + } else { + // SULFATE RICH: solve for SO4= and HSO4- ions + + sum_elec = 3. * electrolyte_liquid(mosaic.jh2so4) + + 2. * electrolyte_liquid(mosaic.jnh4hso4) + + 5. * electrolyte_liquid(mosaic.jlvcite) + + 3. * electrolyte_liquid(mosaic.jnh4so4) + + 2. * electrolyte_liquid(mosaic.jnahso4) + + 5. * electrolyte_liquid(mosaic.jna3hso4) + + 3. * electrolyte_liquid(mosaic.jna2so4) + + 2. * electrolyte_liquid(mosaic.jhno3) + + 2. * electrolyte_liquid(mosaic.jhcl); + + if (sum_elec == 0.0) { + for (ordinal_type jA = 0; jA < mosaic.nelectrolyte; j++) { + gam(jA) = 1.0; + } + // FIXME: duplicated code to avoid goto statement + gam(mosaic.jnh4no3) = 1.0; + gam(mosaic.jnh4cl) = 1.0; + gam(mosaic.jnano3) = 1.0; + gam(mosaic.jnacl) = 1.0; + gam(mosaic.jcano3) = 1.0; + gam(mosaic.jcacl2) = 1.0; + gam(mosaic.jnh4msa) = 1.0; + gam(mosaic.jnamsa) = 1.0; + gam(mosaic.jcamsa2) = 1.0; + + // compute equilibrium pH + // cation molalities (mol / kg water) + mc(mosaic.jc_ca) = 1.e-9 * aer_liquid(mosaic.ica_a) / water_a; + mc(mosaic.jc_nh4) = 1.e-9 * aer_liquid(mosaic.inh4_a) / water_a; + mc(mosaic.jc_na) = 1.e-9 * aer_liquid(mosaic.ina_a) / water_a; + + // anion molalities (mol / kg water) + mSULF = 1.e-9 * aer_liquid(mosaic.iso4_a) / water_a; + ma(mosaic.ja_hso4) = 0.0; + ma(mosaic.ja_so4) = 0.0; + ma(mosaic.ja_no3) = 1.e-9 * aer_liquid(mosaic.ino3_a) / water_a; + ma(mosaic.ja_cl) = 1.e-9 * aer_liquid(mosaic.icl_a) / water_a; + ma(mosaic.ja_msa) = 1.e-9 * aer_liquid(mosaic.imsa_a) / water_a; + + real_type dumK, c_bal, aq, bq, cq; + + gam_ratio = ats::pow(gam(mosaic.jnh4hso4),2.) / + ats::pow(gam(mosaic.jhhso4),2.); + fn_Keq(mosaic.Keq_ll_298(0), mosaic.Keq_a_ll(0), mosaic.Keq_b_ll(0), T_K, dumK); + dumK = dumK * ats::pow(gam(mosaic.jhhso4),2.) / + ats::pow(gam(mosaic.jh2so4),3.); + + + c_bal = mc(mosaic.jc_nh4) + mc(mosaic.jc_na) + 2. * mc(mosaic.jc_ca) - + ma(mosaic.ja_no3) - ma(mosaic.ja_cl) - mSULF - ma(mosaic.ja_msa); + + aq = 1.0; + bq = dumK + c_bal; + cq = dumK * (c_bal -mSULF); + + //--quadratic solution + if (bq != 0.0) { + real_type xq = 4. * (1./bq) * (cq/bq); + + } else { + real_type xq = 1.e6; + } + + if(atsabs(xq) < 1.e-6) { + real_type dum = xq * (0.5 + xq*(0.125 + xq*0.0625)); + real_type quad = (-0.5*bq/aq)*dum; + if (quad < 0.0) { + quad = -bq/aq - quad; + } + } else { + quad = 0.5 * (-bq + atssqrt(bq*bq - 4.*cq)); + } + + //--end of quadratic solution + mc(mosaic.jc_h) = atsmax(quad, 1.e-7); + ma(mosaic.ja_so4) = mSULF * dumK / (mc(mosaic.jc_h) + dumK); + ma(mosaic.ja_hso4) = mSULF - ma(mosaic.ja_so4); + + activity(mosaic.jcamsa2) = mc(mosaic.jc_ca) * + ats::pow(ma(mosaic.ja_msa),2.) * + ats::pow(gam(mosaic.jcamsa2),3.); + + activity(mosaic.jnh4so4) = ats::pow(mc(mosaic.jc_nh4),2.) * + ma(mosaic.ja_so4) * + ats::pow(gam(mosaic.jnh4so4),3.); + + activity(mosaic.jlvcite) = ats::pow(mc(mosaic.jc_nh4),3.) * + ma(mosaic.ja_hso4) * + ma(mosaic.ja_so4) * + ats::pow(gam(mosaic.lvcite),5.); + + activity(mosaic.jnh4hso4) = mc(mosaic.jc_nh4) * + ma(mosaic.ja_hso4) * + ats::pow(gam(mosaic.jnh4hso4),2.); + + activity(mosaic.jnh4msa) = mc(mosaic.jc_nh4) * + ma(mosaic.ja_msa) * + ats::pow(gam(mosaic.jnh4msa),2.); + + activity(mosaic.jna2so4) = ats::pow(mc(mosaic.jc_na),2.) * + ma(mosaic.ja_so4) * + ats::pow(gam(mosaic.jna2so4),3.); + + activity(mosaic.jnahso4) = mc(mosaic.jc_na) * + ma(mosaic.ja_hso4) * + ats::pow(gam(mosaic.jnahso4),2.); + + activity(mosaic.jnamsa) = mc(mosaic.jc_na) * + ma(mosaic.ja_msa) * + ats::pow(gam(mosaic.jnamsa),2.); + + // Note: these lines are also commented out in MOSAIC + // activity(jna3hso4,ibin)= mc(jc_na,ibin)**3 * ma(ja_hso4,ibin) * + // & ma(ja_so4,ibin) * gam(jna3hso4,ibin)**5 + + activity(mosaic.jna3hso4) = 0.0; + + activity(mosaic.jhno3) = mc(mosaic.jc_h) * + ma(mosaic.ja_no3) * + ats::pow(gam(mosaic.jhno3),2.); + + activity(mosaic.jhcl) = mc(mosaic.jc_h) * + ma(mosaic.ja_cl) * + ats::pow(gam(mosaic.jhcl),2.); + + activity(mosaic.jmsa) = mc(mosaic.jc_h) * + ma(mosaic.ja_msa) * + ats::pow(gam(mosaic.jmsa),2.); + + // sulfate-poor species + activity(mosaic.jnh4no3) = 0.0; + activity(mosaic.jnh4cl) = 0.0; + activity(mosaic.jnano3) = 0.0; + activity(mosaic.jnacl) = 0.0; + activity(mosaic.jcano3) = 0.0; + activity(mosaic.jcacl2) = 0.0; + } + + xmol(mosaic.jh2so4) = 3. * electrolyte_liquid(mosaic.jh2so4) / sum_elec; + xmol(mosaic.jnh4hso4) = 2. * electrolyte_liquid(mosaic.jnh4hso4) / sum_elec; + xmol(mosaic.jlvcite) = 5. * electrolyte_liquid(mosaic.jlvcite) / sum_elec; + xmol(mosaic.jnh4so4) = 3. * electrolyte_liquid(mosaic.jnh4so4) / sum_elec; + xmol(mosaic.jnahso4) = 2. * electrolyte_liquid(mosaic.jnahso4) / sum_elec; + xmol(mosaic.jna3hso4) = 5. * electrolyte_liquid(mosaic.jna3hso4) / sum_elec; + xmol(mosaic.jna2so4) = 3. * electrolyte_liquid(mosaic.jna2so4) / sum_elec; + xmol(mosaic.jhno3) = 2. * electrolyte_liquid(mosaic.jhno3) / sum_elec; + xmol(mosaic.jhcl) = 2. * electrolyte_liquid(mosaic.jhcl) / sum_elec; + + // 2H.SO4 + jA = mosaic.jh2so4; + log_gam(jA) = xmol(mosaic.jh2so4) * log_gamZ(jA,mosaic.jh2so4) + + xmol(mosaic.jnh4hso4) * log_gamZ(jA,mosaic.jnh4hso4) + + xmol(mosaic.jlvcite) * log_gamZ(jA,mosaic.jlvcite) + + xmol(mosaic.jnh4so4) * log_gamZ(jA,mosaic.jnh4so4) + + xmol(mosaic.jnahso4) * log_gamZ(jA,mosaic.jnahso4) + + xmol(mosaic.jna3hso4) * log_gamZ(jA,mosaic.jna3hso4) + + xmol(mosaic.jna2so4) * log_gamZ(jA,mosaic.jna2so4) + + xmol(mosaic.jhno3) * log_gamZ(jA,mosaic.jhno3) + + xmol(mosaic.jhcl) * log_gamZ(jA,mosaic.jhcl); + gam(jA) = ats::pow(10.,log_gam(jA)); + + // H.HSO4 + jA = mosaic.jhhso4; + log_gam(jA) = xmol(mosaic.jh2so4) * log_gamZ(jA,mosaic.jh2so4) + + xmol(mosaic.jnh4hso4) * log_gamZ(jA,mosaic.jnh4hso4) + + xmol(mosaic.jlvcite) * log_gamZ(jA,mosaic.jlvcite) + + xmol(mosaic.jnh4so4) * log_gamZ(jA,mosaic.jnh4so4) + + xmol(mosaic.jnahso4) * log_gamZ(jA,mosaic.jnahso4) + + xmol(mosaic.jna3hso4) * log_gamZ(jA,mosaic.jna3hso4) + + xmol(mosaic.jna2so4) * log_gamZ(jA,mosaic.jna2so4) + + xmol(mosaic.jhno3) * log_gamZ(jA,mosaic.jhno3) + + xmol(mosaic.jhcl) * log_gamZ(jA,mosaic.jhcl); + gam(jA) = ats::pow(10.,log_gam(jA)); + + + // NH4HSO4 + jA = mosaic.jnh4hso4; + log_gam(jA) = xmol(mosaic.jh2so4) * log_gamZ(jA,mosaic.jh2so4) + + xmol(mosaic.jnh4hso4) * log_gamZ(jA,mosaic.jnh4hso4) + + xmol(mosaic.jlvcite) * log_gamZ(jA,mosaic.jlvcite) + + xmol(mosaic.jnh4so4) * log_gamZ(jA,mosaic.jnh4so4) + + xmol(mosaic.jnahso4) * log_gamZ(jA,mosaic.jnahso4) + + xmol(mosaic.jna3hso4) * log_gamZ(jA,mosaic.jna3hso4) + + xmol(mosaic.jna2so4) * log_gamZ(jA,mosaic.jna2so4) + + xmol(mosaic.jhno3) * log_gamZ(jA,mosaic.jhno3) + + xmol(mosaic.jhcl) * log_gamZ(jA,mosaic.jhcl); + gam(jA) = ats::pow(10.,log_gam(jA)); + + + // LETOVICITE + jA = mosaic.jlvcite; + log_gam(jA) = xmol(mosaic.jh2so4) * log_gamZ(jA,mosiac.jh2so4) + + xmol(mosaic.jnh4hso4) * log_gamZ(jA,mosaic.jnh4hso4) + + xmol(mosaic.jlvcite) * log_gamZ(jA,mosaic.jlvcite) + + xmol(mosaic.jnh4so4) * log_gamZ(jA,mosaic.jnh4so4) + + xmol(mosaic.jnahso4) * log_gamZ(jA,mosaic.jnahso4) + + xmol(mosaic.jna3hso4) * log_gamZ(jA,mosaic.jna3hso4) + + xmol(mosaic.jna2so4) * log_gamZ(jA,mosaic.jna2so4) + + xmol(mosaic.jhno3) * log_gamZ(jA,mosaic.jhno3) + + xmol(mosaic.jhcl) * log_gamZ(jA,mosaic.jhcl); + gam(jA) = ats::pow(10.,log_gam(jA)); + + + // (NH4)2SO4 + jA = mosaic.jnh4so4; + log_gam(jA) = xmol(mosaic.jh2so4) * log_gamZ(jA,mosaic.jh2so4) + + xmol(mosaic.jnh4hso4) * log_gamZ(jA,mosaic.jnh4hso4) + + xmol(mosaic.jlvcite) * log_gamZ(jA,mosaic.jlvcite) + + xmol(mosaic.jnh4so4) * log_gamZ(jA,mosaic.jnh4so4) + + xmol(mosaic.jnahso4) * log_gamZ(jA,mosaic.jnahso4) + + xmol(mosaic.jna3hso4) * log_gamZ(jA,mosaic.jna3hso4) + + xmol(mosaic.jna2so4) * log_gamZ(jA,mosaic.jna2so4) + + xmol(mosaic.jhno3) * log_gamZ(jA,mosaic.jhno3) + + xmol(mosaic.jhcl) * log_gamZ(jA,mosaic.jhcl); + gam(jA) = ats::pow(10.,log_gam(jA)); + + + // NaHSO4 + jA = mosaic.jnahso4; + log_gam(jA) = xmol(mosaic.jh2so4) * log_gamZ(jA,mosaic.jh2so4) + + xmol(mosaic.jnh4hso4) * log_gamZ(jA,mosaic.jnh4hso4) + + xmol(mosaic.jlvcite) * log_gamZ(jA,mosaic.jlvcite) + + xmol(mosaic.jnh4so4) * log_gamZ(jA,mosaic.jnh4so4) + + xmol(mosaic.jnahso4) * log_gamZ(jA,mosaic.jnahso4) + + xmol(mosaic.jna3hso4) * log_gamZ(jA,mosaic.jna3hso4) + + xmol(mosaic.jna2so4) * log_gamZ(jA,mosaic.jna2so4) + + xmol(mosaic.jhno3) * log_gamZ(jA,mosaic.jhno3) + + xmol(mosaic.jhcl) * log_gamZ(jA,mosaic.jhcl); + gam(jA) = ats::pow(10.,log_gam(jA)); + + + // Na3H(SO4)2 + jA = mosaic.jna3hso4; + // log_gam(jA) = xmol(jh2so4) *log_gamZ(jA,jh2so4) + + // & xmol(jnh4hso4)*log_gamZ(jA,jnh4hso4)+ + // & xmol(jlvcite) *log_gamZ(jA,jlvcite) + + // & xmol(jnh4so4) *log_gamZ(jA,jnh4so4) + + // & xmol(jnahso4) *log_gamZ(jA,jnahso4) + + // & xmol(jna3hso4)*log_gamZ(jA,jna3hso4)+ + // & xmol(jna2so4) *log_gamZ(jA,jna2so4) + + // & xmol(jhno3) *log_gamZ(jA,jhno3) + + // & xmol(jhcl) *log_gamZ(jA,jhcl) + // gam(jA,ibin) = 10.**log_gam(jA) + gam(jA) = 1.0; + + + // Na2SO4 + jA = mosaic.jna2so4; + log_gam(jA) = xmol(mosaic.jh2so4) * log_gamZ(jA,mosaic.jh2so4) + + xmol(mosaic.jnh4hso4) * log_gamZ(jA,mosaic.jnh4hso4) + + xmol(mosaic.jlvcite) * log_gamZ(jA,mosaic.jlvcite) + + xmol(mosaic.jnh4so4) * log_gamZ(jA,mosaic.jnh4so4) + + xmol(mosaic.jnahso4) * log_gamZ(jA,mosaic.jnahso4) + + xmol(mosaic.jna3hso4) * log_gamZ(jA,mosaic.jna3hso4) + + xmol(mosaic.jna2so4) * log_gamZ(jA,mosaic.jna2so4) + + xmol(mosaic.jhno3) * log_gamZ(jA,mosaic.jhno3) + + xmol(mosaic.jhcl) * log_gamZ(jA,mosaic.jhcl); + gam(jA) = ats::pow(10.,log_gam(jA)); + + + // HNO3 + jA = mosaic.jhno3; + log_gam(jA) = xmol(mosaic.jh2so4) * log_gamZ(jA,mosiac.jh2so4) + + xmol(mosaic.jnh4hso4) * log_gamZ(jA,mosaic.jnh4hso4) + + xmol(mosaic.jlvcite) * log_gamZ(jA,mosaic.jlvcite) + + xmol(mosaic.jnh4so4) * log_gamZ(jA,mosaic.jnh4so4) + + xmol(mosaic.jnahso4) * log_gamZ(jA,mosaic.jnahso4) + + xmol(mosaic.jna3hso4) * log_gamZ(jA,mosaic.jna3hso4) + + xmol(mosaic.jna2so4) * log_gamZ(jA,mosaic.jna2so4) + + xmol(mosaic.jhno3) * log_gamZ(jA,mosaic.jhno3) + + xmol(mosaic.jhcl) * log_gamZ(jA,mosaic.jhcl); + gam(jA) = ats::pow(10.,log_gam(jA)); + + + // HCl + jA = mosaic.jhcl; + log_gam(jA) = xmol(mosaic.jh2so4) * log_gamZ(jA,mosaic.jh2so4) + + xmol(mosaic.jnh4hso4) * log_gamZ(jA,mosaic.jnh4hso4) + + xmol(mosaic.jlvcite) * log_gamZ(jA,mosaic.jlvcite) + + xmol(mosaic.jnh4so4) * log_gamZ(jA,mosaic.jnh4so4) + + xmol(moasic.jnahso4) * log_gamZ(jA,mosaic.jnahso4) + + xmol(mosaic.jna3hso4) * log_gamZ(jA,mosaic.jna3hso4) + + xmol(mosaic.jna2so4) * log_gamZ(jA,mosaic.jna2so4) + + xmol(mosaic.jhno3) * log_gamZ(jA,mosaic.jhno3) + + xmol(mosaic.jhcl) * log_gamZ(jA,mosaic.jhcl); + gam(jA) = ats::pow(10.,log_gam(jA)); + + // FIXME: duplicated code to avoid goto statement + gam(mosaic.jnh4no3) = 1.0; + gam(mosaic.jnh4cl) = 1.0; + gam(mosaic.jnano3) = 1.0; + gam(mosaic.jnacl) = 1.0; + gam(mosaic.jcano3) = 1.0; + gam(mosaic.jcacl2) = 1.0; + gam(mosaic.jnh4msa) = 1.0; + gam(mosaic.jnamsa) = 1.0; + gam(mosaic.jcamsa2) = 1.0; + + // compute equilibrium pH + // cation molalities (mol / kg water) + mc(mosaic.jc_ca) = 1.e-9 * aer_liquid(mosaic.ica_a) / water_a; + mc(mosaic.jc_nh4) = 1.e-9 * aer_liquid(mosaic.inh4_a) / water_a; + mc(mosaic.jc_na) = 1.e-9 * aer_liquid(mosaic.ina_a) / water_a; + + // anion molalities (mol / kg water) + mSULF = 1.e-9 * aer_liquid(mosaic.iso4_a) / water_a; + ma(mosaic.ja_hso4) = 0.0; + ma(mosaic.ja_so4) = 0.0; + ma(mosaic.ja_no3) = 1.e-9 * aer_liquid(mosaic.ino3_a) / water_a; + ma(mosaic.ja_cl) = 1.e-9 * aer_liquid(mosaic.icl_a) / water_a; + ma(mosaic.ja_msa) = 1.e-9 * aer_liquid(mosaic.imsa_a) / water_a; + + real_type dumK, c_bal, aq, bq, cq; + + gam_ratio = ats::pow(gam(mosaic.jnh4hso4),2.) / + ats::pow(gam(mosaic.jhhso4),2.); + fn_Keq(mosaic.Keq_ll_298(0), mosaic.Keq_a_ll(0), mosaic.Keq_b_ll(0), T_K, dumK); + dumK = dumK * ats::pow(gam(mosaic.jhhso4),2.) / + ats::pow(gam(mosaic.jh2so4),3.); + + + c_bal = mc(mosaic.jc_nh4) + mc(mosaic.jc_na) + 2. * mc(mosaic.jc_ca) - + ma(mosaic.ja_no3) - ma(mosaic.ja_cl) - mSULF - ma(mosaic.ja_msa); + + aq = 1.0; + bq = dumK + c_bal; + cq = dumK * (c_bal -mSULF); + + //--quadratic solution + if (bq != 0.0) { + real_type xq = 4. * (1./bq) * (cq/bq); + + } else { + real_type xq = 1.e6; + } + + if(atsabs(xq) < 1.e-6) { + real_type dum = xq * (0.5 + xq*(0.125 + xq*0.0625)); + real_type quad = (-0.5*bq/aq)*dum; + if (quad < 0.0) { + quad = -bq/aq - quad; + } + } else { + quad = 0.5 * (-bq + atssqrt(bq*bq - 4.*cq)); + } + + //--end of quadratic solution + mc(mosaic.jc_h) = atsmax(quad, 1.e-7); + ma(mosaic.ja_so4) = mSULF * dumK / (mc(mosaic.jc_h) + dumK); + ma(mosaic.ja_hso4) = mSULF - ma(mosaic.ja_so4); + + activity(mosaic.jcamsa2) = mc(mosaic.jc_ca) * + ats::pow(ma(mosaic.ja_msa),2.) * + ats::pow(gam(mosaic.jcamsa2),3.); + + activity(mosaic.jnh4so4) = ats::pow(mc(mosaic.jc_nh4),2.) * + ma(mosaic.ja_so4) * + ats::pow(gam(mosaic.jnh4so4),3.); + + activity(mosaic.jlvcite) = ats::pow(mc(mosaic.jc_nh4),3.) * + ma(mosaic.ja_hso4) * + ma(mosaic.ja_so4) * + ats::pow(gam(mosaic.lvcite),5.); + + activity(mosaic.jnh4hso4) = mc(mosaic.jc_nh4) * + ma(mosaic.ja_hso4) * + ats::pow(gam(mosaic.jnh4hso4),2.); + + activity(mosaic.jnh4msa) = mc(mosaic.jc_nh4) * + ma(mosaic.ja_msa) * + ats::pow(gam(mosaic.jnh4msa),2.); + + activity(mosaic.jna2so4) = ats::pow(mc(mosaic.jc_na),2.) * + ma(mosaic.ja_so4) * + ats::pow(gam(mosaic.jna2so4),3.); + + activity(mosaic.jnahso4) = mc(mosaic.jc_na) * + ma(mosaic.ja_hso4) * + ats::pow(gam(mosaic.jnahso4),2.); + + activity(mosaic.jnamsa) = mc(mosaic.jc_na) * + ma(mosaic.ja_msa) * + ats::pow(gam(mosaic.jnamsa),2.); + + // Note: these lines are also commented out in MOSAIC + // activity(jna3hso4,ibin)= mc(jc_na,ibin)**3 * ma(ja_hso4,ibin) * + // & ma(ja_so4,ibin) * gam(jna3hso4,ibin)**5 + + activity(mosaic.jna3hso4) = 0.0; + + activity(mosaic.jhno3) = mc(mosaic.jc_h) * + ma(mosaic.ja_no3) * + ats::pow(gam(mosaic.jhno3),2.); + + activity(mosaic.jhcl) = mc(mosaic.jc_h) * + ma(mosaic.ja_cl) * + ats::pow(gam(mosaic.jhcl),2.); + + activity(mosaic.jmsa) = mc(mosaic.jc_h) * + ma(mosaic.ja_msa) * + ats::pow(gam(mosaic.jmsa),2.); + + // sulfate-poor species + activity(mosaic.jnh4no3) = 0.0; + activity(mosaic.jnh4cl) = 0.0; + activity(mosaic.jnano3) = 0.0; + activity(mosaic.jnacl) = 0.0; + activity(mosaic.jcano3) = 0.0; + activity(mosaic.jcacl2) = 0.0; + } + } + + KOKKOS_INLINE_FUNCTION static void adjust_solid_aerosol(const MosaicModelData& mosaic, const real_type_1d_view_type& aer_solid, From 7022b40cf3e4b5bf4c63644498d98d3d8634c2b5 Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Thu, 21 Aug 2025 15:56:47 -0500 Subject: [PATCH 2/4] fixes to implementation of compute_activities --- src/core/impl/TChem_Impl_MOSAIC.hpp | 101 ++++++++++++++-------------- 1 file changed, 50 insertions(+), 51 deletions(-) diff --git a/src/core/impl/TChem_Impl_MOSAIC.hpp b/src/core/impl/TChem_Impl_MOSAIC.hpp index b5e764d7..9b2108e4 100644 --- a/src/core/impl/TChem_Impl_MOSAIC.hpp +++ b/src/core/impl/TChem_Impl_MOSAIC.hpp @@ -2423,9 +2423,12 @@ struct MOSAIC{ KOKKOS_INLINE_FUNCTION static void compute_activities(const MosaicModelData& mosaic, + const real_type_1d_view_type& molalities, + const real_type_1d_view_type& xmol, const real_type_1d_view_type& aer_liquid, const real_type_1d_view_type& ma, const real_type_1d_view_type& mc, + const real_type_1d_view_type& Keq_ll, const real_type_1d_view_type& electrolyte_solid, const real_type_1d_view_type& electrolyte_liquid, const real_type_1d_view_type& electrolyte_total, @@ -2433,17 +2436,17 @@ struct MOSAIC{ const real_type_2d_view_type& log_gamZ, const real_type_1d_view_type& gam, const real_type_1d_view_type& activity, - real_type& T_K) { + real_type& jaerosolstate, + real_type& jphase, + real_type& jhyst_leg, + real_type& aH2O_a) { // local variables - real_type_1d_view_type xmol("xmol", mosaic.nelectrolyte); - real_type a_c, Keq_ll, sum_elec, gam_ratio, water_a, mSULF, gam_ratio; + real_type a_c, sum_elec, gam_ratio, water_a, mSULF; ordinal_type jA; - // FIXME: these need to be set beforehand - ordinal_type jaerosolstate, jphase, jhyst_leg; // get aerosol water activity - aerosol_water(mosaic, electrolyte_liquid, water_a, jaerosolstate, jphase, jhyst_leg); + aerosol_water(mosaic,electrolyte_liquid,aH2O_a,molalities,jaerosolstate,jphase,jhyst_leg,water_a); if (water_a == 0.0) { return; @@ -2451,7 +2454,7 @@ struct MOSAIC{ // get sulfate ratio to determine regime real_type XT = 0.0; - calcuate_XT(aer_liquid, mosaic, XT); + calculate_XT(mosaic, aer_liquid, XT); if (XT > 2.0 || XT < 0.0) { // SULFATE POOR: fully dissociated electrolytes @@ -2466,7 +2469,7 @@ struct MOSAIC{ // cation molalities (mol / kg water) mc(mosaic.jc_ca) = 1.e-9 * aer_liquid(mosaic.ica_a) / water_a; mc(mosaic.jc_nh4) = 1.e-9 * aer_liquid(mosaic.inh4_a) / water_a; - mc(mosaic.jc_na) = 1.e-9 * aer_liquid(mosaic.na_a) / water_a; + mc(mosaic.jc_na) = 1.e-9 * aer_liquid(mosaic.ina_a) / water_a; a_c = ( (2. * ma(mosaic.ja_so4) + ma(mosaic.ja_no3) + @@ -2475,12 +2478,8 @@ struct MOSAIC{ (2. * mc(mosaic.jc_ca) + mc(mosaic.jc_nh4) + mc(mosaic.jc_na)) ); - // FIXME: consider adding update_thermodynamic_constants instead - // of computing equil. constants directly. - // Make sure temp. is handled - fn_Keq(mosaic.Keq_ll_298(2), mosaic.Keq_a_ll(2), mosiac.Keq_b_ll(2), T_K, Keq_ll); mc(mosaic.jc_h) = 0.5 * ( (a_c) + - (ats::sqrt(a_c*a_c + 4. * Keq_ll)) ); + (ats::sqrt(a_c*a_c + 4. * Keq_ll(2))) ); if (mc(mosaic.jc_h) == 0.0) { mc(mosaic.jc_h) = 1.e-10; @@ -2541,7 +2540,7 @@ struct MOSAIC{ xmol(mosaic.jhno3) = 2. * electrolyte_liquid(mosaic.jhno3) / sum_elec; xmol(mosaic.jhcl) = 2. * electrolyte_liquid(mosaic.jhcl) / sum_elec; - jA = mosaic.nh4so4; + jA = mosaic.jnh4so4; if (xmol(jA) > 0.0) { log_gam(jA) = xmol(mosaic.jnh4no3) * log_gamZ(jA,mosaic.jnh4no3) + xmol(mosaic.jnh4cl) * log_gamZ(jA,mosaic.jnh4cl) + @@ -2556,7 +2555,7 @@ struct MOSAIC{ gam(jA) = ats::pow(10., log_gam(jA)); activity(jA) = ats::pow(mc(mosaic.jc_nh4),2.) * ma(mosaic.ja_so4) * - ats::pow(gam(jA,3.)); + ats::pow(gam(jA),3.); } jA = mosaic.jnh4no3; @@ -2574,7 +2573,7 @@ struct MOSAIC{ gam(jA) = ats::pow(10., log_gam(jA)); activity(jA) = mc(mosaic.jc_nh4) * ma(mosaic.ja_no3) * - ats::pow(gam(jA,2.)); + ats::pow(gam(jA),2.); } jA = mosaic.jnh4cl; @@ -2592,7 +2591,7 @@ struct MOSAIC{ gam(jA) = ats::pow(10., log_gam(jA)); activity(jA) = mc(mosaic.jc_nh4) * ma(mosaic.ja_cl) * - ats::pow(gam(jA,2.)); + ats::pow(gam(jA),2.); } jA = mosaic.jna2so4; @@ -2610,7 +2609,7 @@ struct MOSAIC{ gam(jA) = ats::pow(10., log_gam(jA)); activity(jA) = ats::pow(mc(mosaic.jc_na),2.) * ma(mosaic.ja_so4) * - ats::pow(gam(jA,3.)); + ats::pow(gam(jA),3.); } jA = mosaic.jnano3; @@ -2628,7 +2627,7 @@ struct MOSAIC{ gam(jA) = ats::pow(10., log_gam(jA)); activity(jA) = mc(mosaic.jc_na) * ma(mosaic.ja_no3) * - ats::pow(gam(jA,2.)); + ats::pow(gam(jA),2.); } jA = mosaic.jnacl; @@ -2646,7 +2645,7 @@ struct MOSAIC{ gam(jA) = ats::pow(10., log_gam(jA)); activity(jA) = mc(mosaic.jc_na) * ma(mosaic.ja_cl) * - ats::pow(gam(jA,2.)); + ats::pow(gam(jA),2.); } // Note: these are commented out in MOSAIC also. @@ -2677,7 +2676,7 @@ struct MOSAIC{ gam(jA) = ats::pow(10., log_gam(jA)); activity(jA) = mc(mosaic.jc_ca) * ats::pow(ma(mosaic.ja_no3),2.) * - ats::pow(gam(jA,3.)); + ats::pow(gam(jA),3.); } jA = mosaic.jcacl2; @@ -2695,7 +2694,7 @@ struct MOSAIC{ gam(jA) = ats::pow(10., log_gam(jA)); activity(jA) = mc(mosaic.jc_ca) * ats::pow(ma(mosaic.ja_cl),2.) * - ats::pow(gam(jA,3.)); + ats::pow(gam(jA),3.); } jA = mosaic.jhno3; @@ -2713,7 +2712,7 @@ struct MOSAIC{ gam(jA) = ats::pow(10., log_gam(jA)); activity(jA) = mc(mosaic.jc_h) * ma(mosaic.ja_no3) * - ats::pow(gam(jA,2.)); + ats::pow(gam(jA),2.); } jA = mosaic.jhcl; @@ -2731,7 +2730,7 @@ struct MOSAIC{ gam(jA) = ats::pow(10., log_gam(jA)); activity(jA) = mc(mosaic.jc_h) * ma(mosaic.ja_cl) * - ats::pow(gam(jA,2.)); + ats::pow(gam(jA),2.); } // FIXME: duplicated code to avoid goto statement @@ -2774,7 +2773,7 @@ struct MOSAIC{ 2. * electrolyte_liquid(mosaic.jhcl); if (sum_elec == 0.0) { - for (ordinal_type jA = 0; jA < mosaic.nelectrolyte; j++) { + for (ordinal_type jA = 0; jA < mosaic.nelectrolyte; jA++) { gam(jA) = 1.0; } // FIXME: duplicated code to avoid goto statement @@ -2806,10 +2805,8 @@ struct MOSAIC{ gam_ratio = ats::pow(gam(mosaic.jnh4hso4),2.) / ats::pow(gam(mosaic.jhhso4),2.); - fn_Keq(mosaic.Keq_ll_298(0), mosaic.Keq_a_ll(0), mosaic.Keq_b_ll(0), T_K, dumK); - dumK = dumK * ats::pow(gam(mosaic.jhhso4),2.) / - ats::pow(gam(mosaic.jh2so4),3.); - + dumK = Keq_ll(0) * ats::pow(gam(mosaic.jhhso4),2.) / + ats::pow(gam(mosaic.jh2so4),3.); c_bal = mc(mosaic.jc_nh4) + mc(mosaic.jc_na) + 2. * mc(mosaic.jc_ca) - ma(mosaic.ja_no3) - ma(mosaic.ja_cl) - mSULF - ma(mosaic.ja_msa); @@ -2818,26 +2815,28 @@ struct MOSAIC{ bq = dumK + c_bal; cq = dumK * (c_bal -mSULF); + real_type xq; //--quadratic solution if (bq != 0.0) { - real_type xq = 4. * (1./bq) * (cq/bq); + xq = 4. * (1./bq) * (cq/bq); } else { - real_type xq = 1.e6; + xq = 1.e6; } - if(atsabs(xq) < 1.e-6) { + real_type quad = 0.0; + if(ats::abs(xq) < 1.e-6) { real_type dum = xq * (0.5 + xq*(0.125 + xq*0.0625)); - real_type quad = (-0.5*bq/aq)*dum; + quad = (-0.5*bq/aq)*dum; if (quad < 0.0) { quad = -bq/aq - quad; } } else { - quad = 0.5 * (-bq + atssqrt(bq*bq - 4.*cq)); + quad = 0.5 * (-bq + ats::sqrt(bq*bq - 4.*cq)); } //--end of quadratic solution - mc(mosaic.jc_h) = atsmax(quad, 1.e-7); + mc(mosaic.jc_h) = max(quad, 1.e-7); ma(mosaic.ja_so4) = mSULF * dumK / (mc(mosaic.jc_h) + dumK); ma(mosaic.ja_hso4) = mSULF - ma(mosaic.ja_so4); @@ -2852,7 +2851,7 @@ struct MOSAIC{ activity(mosaic.jlvcite) = ats::pow(mc(mosaic.jc_nh4),3.) * ma(mosaic.ja_hso4) * ma(mosaic.ja_so4) * - ats::pow(gam(mosaic.lvcite),5.); + ats::pow(gam(mosaic.jlvcite),5.); activity(mosaic.jnh4hso4) = mc(mosaic.jc_nh4) * ma(mosaic.ja_hso4) * @@ -2954,7 +2953,7 @@ struct MOSAIC{ // LETOVICITE jA = mosaic.jlvcite; - log_gam(jA) = xmol(mosaic.jh2so4) * log_gamZ(jA,mosiac.jh2so4) + + log_gam(jA) = xmol(mosaic.jh2so4) * log_gamZ(jA,mosaic.jh2so4) + xmol(mosaic.jnh4hso4) * log_gamZ(jA,mosaic.jnh4hso4) + xmol(mosaic.jlvcite) * log_gamZ(jA,mosaic.jlvcite) + xmol(mosaic.jnh4so4) * log_gamZ(jA,mosaic.jnh4so4) + @@ -3025,7 +3024,7 @@ struct MOSAIC{ // HNO3 jA = mosaic.jhno3; - log_gam(jA) = xmol(mosaic.jh2so4) * log_gamZ(jA,mosiac.jh2so4) + + log_gam(jA) = xmol(mosaic.jh2so4) * log_gamZ(jA,mosaic.jh2so4) + xmol(mosaic.jnh4hso4) * log_gamZ(jA,mosaic.jnh4hso4) + xmol(mosaic.jlvcite) * log_gamZ(jA,mosaic.jlvcite) + xmol(mosaic.jnh4so4) * log_gamZ(jA,mosaic.jnh4so4) + @@ -3043,7 +3042,7 @@ struct MOSAIC{ xmol(mosaic.jnh4hso4) * log_gamZ(jA,mosaic.jnh4hso4) + xmol(mosaic.jlvcite) * log_gamZ(jA,mosaic.jlvcite) + xmol(mosaic.jnh4so4) * log_gamZ(jA,mosaic.jnh4so4) + - xmol(moasic.jnahso4) * log_gamZ(jA,mosaic.jnahso4) + + xmol(mosaic.jnahso4) * log_gamZ(jA,mosaic.jnahso4) + xmol(mosaic.jna3hso4) * log_gamZ(jA,mosaic.jna3hso4) + xmol(mosaic.jna2so4) * log_gamZ(jA,mosaic.jna2so4) + xmol(mosaic.jhno3) * log_gamZ(jA,mosaic.jhno3) + @@ -3079,9 +3078,8 @@ struct MOSAIC{ gam_ratio = ats::pow(gam(mosaic.jnh4hso4),2.) / ats::pow(gam(mosaic.jhhso4),2.); - fn_Keq(mosaic.Keq_ll_298(0), mosaic.Keq_a_ll(0), mosaic.Keq_b_ll(0), T_K, dumK); - dumK = dumK * ats::pow(gam(mosaic.jhhso4),2.) / - ats::pow(gam(mosaic.jh2so4),3.); + dumK = Keq_ll(0) * ats::pow(gam(mosaic.jhhso4),2.) / + ats::pow(gam(mosaic.jh2so4),3.); c_bal = mc(mosaic.jc_nh4) + mc(mosaic.jc_na) + 2. * mc(mosaic.jc_ca) - @@ -3092,25 +3090,26 @@ struct MOSAIC{ cq = dumK * (c_bal -mSULF); //--quadratic solution + real_type xq; if (bq != 0.0) { - real_type xq = 4. * (1./bq) * (cq/bq); - + xq = 4. * (1./bq) * (cq/bq); } else { - real_type xq = 1.e6; + xq = 1.e6; } - if(atsabs(xq) < 1.e-6) { + real_type quad = 0.0; + if(ats::abs(xq) < 1.e-6) { real_type dum = xq * (0.5 + xq*(0.125 + xq*0.0625)); - real_type quad = (-0.5*bq/aq)*dum; + quad = (-0.5*bq/aq)*dum; if (quad < 0.0) { quad = -bq/aq - quad; } } else { - quad = 0.5 * (-bq + atssqrt(bq*bq - 4.*cq)); + quad = 0.5 * (-bq + ats::sqrt(bq*bq - 4.*cq)); } //--end of quadratic solution - mc(mosaic.jc_h) = atsmax(quad, 1.e-7); + mc(mosaic.jc_h) = max(quad, 1.e-7); ma(mosaic.ja_so4) = mSULF * dumK / (mc(mosaic.jc_h) + dumK); ma(mosaic.ja_hso4) = mSULF - ma(mosaic.ja_so4); @@ -3125,7 +3124,7 @@ struct MOSAIC{ activity(mosaic.jlvcite) = ats::pow(mc(mosaic.jc_nh4),3.) * ma(mosaic.ja_hso4) * ma(mosaic.ja_so4) * - ats::pow(gam(mosaic.lvcite),5.); + ats::pow(gam(mosaic.jlvcite),5.); activity(mosaic.jnh4hso4) = mc(mosaic.jc_nh4) * ma(mosaic.ja_hso4) * @@ -3173,7 +3172,7 @@ struct MOSAIC{ activity(mosaic.jcano3) = 0.0; activity(mosaic.jcacl2) = 0.0; } - } + } // compute_activities KOKKOS_INLINE_FUNCTION static From 9ecbc60d54662ac9a42a9c70b0e625e647b46761 Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Thu, 21 Aug 2025 15:57:13 -0500 Subject: [PATCH 3/4] add verification data --- .../mosaic/compute_activities_input_ts_0.yaml | 16 ++++++++++ .../mosaic/compute_activities_output_ts_0.py | 29 +++++++++++++++++++ 2 files changed, 45 insertions(+) create mode 100644 src/verification/data_sets/mosaic/compute_activities_input_ts_0.yaml create mode 100644 src/verification/data_sets/mosaic/compute_activities_output_ts_0.py diff --git a/src/verification/data_sets/mosaic/compute_activities_input_ts_0.yaml b/src/verification/data_sets/mosaic/compute_activities_input_ts_0.yaml new file mode 100644 index 00000000..bd43f518 --- /dev/null +++ b/src/verification/data_sets/mosaic/compute_activities_input_ts_0.yaml @@ -0,0 +1,16 @@ + +TChem-atm: + function: compute_activites +input: + fixed: + aer: [ 0.32545627929929049E-004, 0.13027173964710942E-004, 0.42984345673869127E-005, 0.82416864391955953E-004, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000] + ma: [ 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000] + mc: [ 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000] + Keq_ll: [ 0.12797357538504046E-001, 0.17313070459631448E-004, 0.59636669431836358E-014] + electrolyte: [ 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.32545627929929049E-004, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.13027173964710942E-004, 0.42984345673869127E-005, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.32545627929929049E-004, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.13027173964710942E-004, 0.42984345673869127E-005, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000] + log_gam: [ -0.92467696677458555E+000, -0.73558612830028736E+000, -0.44530213685759046E+000, 0.00000000000000000E+000, -0.48377416493643166E+000, -0.35781829712747104E+000, -0.92589087586115615E+000, 0.00000000000000000E+000, -0.43873079562974404E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.37269333516488135E+000, -0.18307245353063939E+000, 0.00000000000000000E+000, -0.12438617187313614E+001, 0.00000000000000000E+000, -0.70371630269205709E+000, -0.57720342022839766E+000, -0.68091938411426067E+000, 0.00000000000000000E+000, 0.00000000000000000E+000] + log_gamZ: [ -0.93189930255005704E+000, -0.73716119980468142E+000, -0.44577181799734644E+000, 0.00000000000000000E+000, -0.48516123661060950E+000, -0.35824982900254998E+000, -0.92838576344943036E+000, -0.71895656453287016E+000, -0.43924972001451246E+000, 0.00000000000000000E+000, -0.47002801416978679E+000, -0.34326084239873467E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.12485196995199486E+001, 0.00000000000000000E+000, -0.72387479422660039E+000, -0.59639681946609091E+000, -0.68315252942358407E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.92474688437480501E+000, -0.69926707746042083E+000, -0.36106819252785183E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.85566704443152819E+000, -0.63715711243095408E+000, -0.30931830533264870E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.98857066622077516E+000, 0.00000000000000000E+000, -0.49100585329946567E+000, -0.35049605916044158E+000, -0.40865287491596147E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.92429529547176159E+000, -0.65327411002576863E+000, -0.24672921352339339E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.81207287628608049E+000, -0.55232802614247523E+000, -0.16253576653315793E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.76116345224585191E+000, 0.00000000000000000E+000, -0.27132214435656321E+000, -0.13197884978368041E+000, -0.12419269655455178E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.84268251333158739E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.54072997317400118E+000, -0.38843884222877945E+000, -0.76909927637839148E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.48539306240114310E+000, -0.33330159506068568E+000, -0.62028144560086096E+000, -0.41735682707271637E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.31553841513790370E+000, -0.16329688481765947E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.74081719019098591E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.39996600807592153E+000, -0.25272268686704069E+000, -0.63394131280366162E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.31969834446098777E+000, -0.17223866159708801E+000, -0.43148961248034090E+000, -0.23495071802407430E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.21579872230026054E+000, -0.68573302435263161E-001, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.91386176601378111E+000, -0.69151260295722339E+000, -0.37414285555142657E+000, 0.00000000000000000E+000, -0.35750246132672459E+000, -0.24663825094759573E+000, -0.86715447959584879E+000, -0.65668126854752096E+000, -0.33965226709847846E+000, 0.00000000000000000E+000, -0.32291248585118193E+000, -0.21221513522311730E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.76212642211813897E+000, 0.00000000000000000E+000, -0.25206592389288485E+000, -0.14126199542059226E+000, -0.26879101160471208E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.88794598644565670E+000, -0.66004236981881526E+000, -0.31807210878087866E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.82460488646160601E+000, -0.60290255816614446E+000, -0.27053121266692093E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.65949517007424330E+000, 0.00000000000000000E+000, -0.17643545170660024E+000, -0.56732437173039330E-001, -0.14667486990129408E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.90375716410704676E+000, -0.63154831605807704E+000, -0.21856337793888311E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.79763442970371790E+000, -0.54101285577646596E+000, -0.14730118223821348E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.62159415809540219E+000, 0.00000000000000000E+000, -0.11966639210362451E+000, -0.39291344656318117E-002, -0.13311269959209948E-001, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.77829799606635319E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.48473747198382222E+000, -0.29750110495454352E+000, -0.68788383024268596E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.41595844139322358E+000, -0.23024203806043220E+000, -0.59489113169317809E+000, -0.34538240936901166E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.14000541414217343E+000, 0.47283865928515190E-001, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.69189741085576129E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.30627672640917325E+000, -0.18334337356266686E+000, -0.61184103977199111E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.24608403864154998E+000, -0.12281989101267698E+000, -0.48693706171468293E+000, -0.32553646848834550E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.36959699061424489E-001, 0.86191756701031075E-001, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.37643115195914645E+000, -0.21415579579471544E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.33847900247745488E+000, -0.17605266630977279E+000, -0.44909135338490402E+000, -0.23272026073336072E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.88409030365160235E-001, 0.73884502343463776E-001, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.20889666464046652E+000, -0.82604471626217313E-001, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.23673497803513835E+000, -0.11035688090914952E+000, -0.28957007910797139E+000, -0.12112804526524057E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.24576592120198582E-001, 0.15091502946842761E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.92199339280512693E+000, -0.57778783361446551E+000, -0.61413722635333706E-001, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.78262380664156628E+000, -0.45227977979490064E+000, 0.43202129297928549E-001, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.60793975837961933E+000, 0.00000000000000000E+000, -0.40556894604046878E-001, 0.43877659826288884E-001, 0.17412366449440908E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.97335984319178859E+000, -0.66884343044988237E+000, -0.21221446722070214E+000, 0.00000000000000000E+000, -0.31117078210445492E+000, -0.21564900814546784E+000, -0.84297665288144386E+000, -0.55159429161722595E+000, -0.11450399801479416E+000, 0.00000000000000000E+000, -0.21348555951083203E+000, -0.11785956687048982E+000, -0.37667313372424305E+000, -0.24915520902784838E+000, 0.00000000000000000E+000, -0.60502769633121067E+000, 0.00000000000000000E+000, -0.34987645133708556E-001, 0.60508989323065254E-001, 0.63992045030697220E-001, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.92903116364439087E+000, -0.55470528010339715E+000, 0.70222475748218383E-002, 0.00000000000000000E+000, -0.17300581312558116E+000, -0.94420335038213010E-001, -0.78490792416217114E+000, -0.42493027959645624E+000, 0.11513136751836117E+000, 0.00000000000000000E+000, -0.64950070557673545E-001, 0.13678990616324382E-001, -0.20787751606031502E+000, -0.10304066350517660E+000, 0.00000000000000000E+000, -0.65539316118879887E+000, 0.00000000000000000E+000, 0.32233577143845338E-001, 0.11087150116846933E+000, 0.21220896078684559E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000] + jaerosolstate: [ 2] + jphase: [ 2] + jhyst_leg: [ 1] + aH2O_a: [ 0.86918676121422112E+000] diff --git a/src/verification/data_sets/mosaic/compute_activities_output_ts_0.py b/src/verification/data_sets/mosaic/compute_activities_output_ts_0.py new file mode 100644 index 00000000..f435918b --- /dev/null +++ b/src/verification/data_sets/mosaic/compute_activities_output_ts_0.py @@ -0,0 +1,29 @@ + +# This file was generated by PartMC-mosaic. + +from math import nan as nan, inf as inf + +# Object is just a dynamic container that stores input/output data. +class Object(object): + pass +# Settings are stored here. +settings = Object() +# Input is stored here. +input = Object() +input.aer=[[ 0.32545627929929049E-004, 0.13027173964710942E-004, 0.42984345673869127E-005, 0.82416864391955953E-004, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000,],] +input.ma=[[ 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000,],] +input.mc=[[ 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000,],] +input.Keq_ll=[[ 0.12797357538504046E-001, 0.17313070459631448E-004, 0.59636669431836358E-014,],] +input.electrolyte=[[ 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.32545627929929049E-004, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.13027173964710942E-004, 0.42984345673869127E-005, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.32545627929929049E-004, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.13027173964710942E-004, 0.42984345673869127E-005, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000,],] +input.log_gam=[[ -0.92467696677458555E+000, -0.73558612830028736E+000, -0.44530213685759046E+000, 0.00000000000000000E+000, -0.48377416493643166E+000, -0.35781829712747104E+000, -0.92589087586115615E+000, 0.00000000000000000E+000, -0.43873079562974404E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.37269333516488135E+000, -0.18307245353063939E+000, 0.00000000000000000E+000, -0.12438617187313614E+001, 0.00000000000000000E+000, -0.70371630269205709E+000, -0.57720342022839766E+000, -0.68091938411426067E+000, 0.00000000000000000E+000, 0.00000000000000000E+000,],] +input.log_gamZ=[[ -0.93189930255005704E+000, -0.73716119980468142E+000, -0.44577181799734644E+000, 0.00000000000000000E+000, -0.48516123661060950E+000, -0.35824982900254998E+000, -0.92838576344943036E+000, -0.71895656453287016E+000, -0.43924972001451246E+000, 0.00000000000000000E+000, -0.47002801416978679E+000, -0.34326084239873467E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.12485196995199486E+001, 0.00000000000000000E+000, -0.72387479422660039E+000, -0.59639681946609091E+000, -0.68315252942358407E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.92474688437480501E+000, -0.69926707746042083E+000, -0.36106819252785183E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.85566704443152819E+000, -0.63715711243095408E+000, -0.30931830533264870E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.98857066622077516E+000, 0.00000000000000000E+000, -0.49100585329946567E+000, -0.35049605916044158E+000, -0.40865287491596147E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.92429529547176159E+000, -0.65327411002576863E+000, -0.24672921352339339E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.81207287628608049E+000, -0.55232802614247523E+000, -0.16253576653315793E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.76116345224585191E+000, 0.00000000000000000E+000, -0.27132214435656321E+000, -0.13197884978368041E+000, -0.12419269655455178E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.84268251333158739E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.54072997317400118E+000, -0.38843884222877945E+000, -0.76909927637839148E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.48539306240114310E+000, -0.33330159506068568E+000, -0.62028144560086096E+000, -0.41735682707271637E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.31553841513790370E+000, -0.16329688481765947E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.74081719019098591E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.39996600807592153E+000, -0.25272268686704069E+000, -0.63394131280366162E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.31969834446098777E+000, -0.17223866159708801E+000, -0.43148961248034090E+000, -0.23495071802407430E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.21579872230026054E+000, -0.68573302435263161E-001, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.91386176601378111E+000, -0.69151260295722339E+000, -0.37414285555142657E+000, 0.00000000000000000E+000, -0.35750246132672459E+000, -0.24663825094759573E+000, -0.86715447959584879E+000, -0.65668126854752096E+000, -0.33965226709847846E+000, 0.00000000000000000E+000, -0.32291248585118193E+000, -0.21221513522311730E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.76212642211813897E+000, 0.00000000000000000E+000, -0.25206592389288485E+000, -0.14126199542059226E+000, -0.26879101160471208E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.88794598644565670E+000, -0.66004236981881526E+000, -0.31807210878087866E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.82460488646160601E+000, -0.60290255816614446E+000, -0.27053121266692093E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.65949517007424330E+000, 0.00000000000000000E+000, -0.17643545170660024E+000, -0.56732437173039330E-001, -0.14667486990129408E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.90375716410704676E+000, -0.63154831605807704E+000, -0.21856337793888311E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.79763442970371790E+000, -0.54101285577646596E+000, -0.14730118223821348E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.62159415809540219E+000, 0.00000000000000000E+000, -0.11966639210362451E+000, -0.39291344656318117E-002, -0.13311269959209948E-001, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.77829799606635319E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.48473747198382222E+000, -0.29750110495454352E+000, -0.68788383024268596E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.41595844139322358E+000, -0.23024203806043220E+000, -0.59489113169317809E+000, -0.34538240936901166E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.14000541414217343E+000, 0.47283865928515190E-001, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.69189741085576129E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.30627672640917325E+000, -0.18334337356266686E+000, -0.61184103977199111E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.24608403864154998E+000, -0.12281989101267698E+000, -0.48693706171468293E+000, -0.32553646848834550E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.36959699061424489E-001, 0.86191756701031075E-001, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.37643115195914645E+000, -0.21415579579471544E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.33847900247745488E+000, -0.17605266630977279E+000, -0.44909135338490402E+000, -0.23272026073336072E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.88409030365160235E-001, 0.73884502343463776E-001, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.20889666464046652E+000, -0.82604471626217313E-001, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.23673497803513835E+000, -0.11035688090914952E+000, -0.28957007910797139E+000, -0.12112804526524057E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.24576592120198582E-001, 0.15091502946842761E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.92199339280512693E+000, -0.57778783361446551E+000, -0.61413722635333706E-001, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.78262380664156628E+000, -0.45227977979490064E+000, 0.43202129297928549E-001, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.60793975837961933E+000, 0.00000000000000000E+000, -0.40556894604046878E-001, 0.43877659826288884E-001, 0.17412366449440908E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.97335984319178859E+000, -0.66884343044988237E+000, -0.21221446722070214E+000, 0.00000000000000000E+000, -0.31117078210445492E+000, -0.21564900814546784E+000, -0.84297665288144386E+000, -0.55159429161722595E+000, -0.11450399801479416E+000, 0.00000000000000000E+000, -0.21348555951083203E+000, -0.11785956687048982E+000, -0.37667313372424305E+000, -0.24915520902784838E+000, 0.00000000000000000E+000, -0.60502769633121067E+000, 0.00000000000000000E+000, -0.34987645133708556E-001, 0.60508989323065254E-001, 0.63992045030697220E-001, 0.00000000000000000E+000, 0.00000000000000000E+000, -0.92903116364439087E+000, -0.55470528010339715E+000, 0.70222475748218383E-002, 0.00000000000000000E+000, -0.17300581312558116E+000, -0.94420335038213010E-001, -0.78490792416217114E+000, -0.42493027959645624E+000, 0.11513136751836117E+000, 0.00000000000000000E+000, -0.64950070557673545E-001, 0.13678990616324382E-001, -0.20787751606031502E+000, -0.10304066350517660E+000, 0.00000000000000000E+000, -0.65539316118879887E+000, 0.00000000000000000E+000, 0.32233577143845338E-001, 0.11087150116846933E+000, 0.21220896078684559E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000,],] +input.jaerosolstate=[[ 2],] +input.jphase=[[ 2],] +input.jhyst_leg=[[ 1],] +input.aH2O_a=[[ 0.86918676121422112E+000],] +# Output data is stored here. +output = Object() +output.jaerosolstate=[[ 2],] +output.jphase=[[ 2],] +output.jhyst_leg=[[ 1],] +output.activity=[[ 0.29345823239350449E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.86211053373726365E+000, 0.52539459929090893E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.53046328639500435E-008, 0.32388619957082867E-008, 0.00000000000000000E+000, 0.00000000000000000E+000, 0.00000000000000000E+000,],] From 4519812c9d56f30682c944431ef0f4a21c24b233 Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Thu, 21 Aug 2025 15:57:33 -0500 Subject: [PATCH 4/4] testing --- src/verification/mosaic/CMakeLists.txt | 4 + .../mosaic/compute_activities.cpp | 134 ++++++++++++++++++ src/verification/mosaic/mosaic_driver.cpp | 4 + 3 files changed, 142 insertions(+) create mode 100644 src/verification/mosaic/compute_activities.cpp diff --git a/src/verification/mosaic/CMakeLists.txt b/src/verification/mosaic/CMakeLists.txt index 8e9d7a8b..c6539004 100644 --- a/src/verification/mosaic/CMakeLists.txt +++ b/src/verification/mosaic/CMakeLists.txt @@ -5,6 +5,7 @@ add_executable(${DRIVER_NAME}.x ${DRIVER_NAME}.cpp check_aerosol_mass.cpp adjust_liquid_aerosol.cpp + compute_activities.cpp adjust_solid_aerosol.cpp do_full_deliquescence.cpp calculate_XT.cpp @@ -28,6 +29,7 @@ TARGET_LINK_LIBRARIES(${DRIVER_NAME}.x verification;${TCHEM_ATM_LINK_LIBRARIES}) SET(TEST_LIST check_aerosol_mass_input_ts_0 adjust_liquid_aerosol_input_ts_0 + compute_activities_input_ts_0 adjust_solid_aerosol_input_ts_0 do_full_deliquescence_input_ts_0 calculate_XT_input_ts_0 @@ -46,6 +48,7 @@ SET(TEST_LIST SET(TEST_BASELINE check_aerosol_mass_output_ts_0 adjust_liquid_aerosol_output_ts_0 + compute_activities_output_ts_0 adjust_solid_aerosol_output_ts_0 do_full_deliquescence_output_ts_0 calculate_XT_output_ts_0 @@ -65,6 +68,7 @@ SET(DEFAULT_TOL 1e-9) set(ERROR_THRESHOLDS ${DEFAULT_TOL} # check_aerosol_mass_output_ts_0 ${DEFAULT_TOL} # adjust_liquid_aerosol_output_ts_0 + ${DEFAULT_TOL} # compute_acivities_output_ts_0 ${DEFAULT_TOL} # adjust_solid_aerosol_output_ts_0 ${DEFAULT_TOL} # do_full_deliquescence_output_ts_0 ${DEFAULT_TOL} # calculate_XT_output_ts_0 diff --git a/src/verification/mosaic/compute_activities.cpp b/src/verification/mosaic/compute_activities.cpp new file mode 100644 index 00000000..c09ec7e1 --- /dev/null +++ b/src/verification/mosaic/compute_activities.cpp @@ -0,0 +1,134 @@ +#include "TChem.hpp" +#include "TChem_Impl_MOSAIC.hpp" +#include +#include "skywalker.hpp" + +using device_type = typename Tines::UseThisDevice::type; +using real_type_1d_view = TChem::real_type_1d_view; +using real_type_2d_view = TChem::real_type_2d_view; +using ordinal_type = TChem::ordinal_type; +using namespace skywalker; +using namespace TChem; + +void compute_activities(Ensemble *ensemble) { + ensemble->process([=](const Input &input, Output &output) { + + const auto aer_arr = input.get_array("aer"); + const auto ma_arr = input.get_array("ma"); + const auto mc_arr = input.get_array("mc"); + const auto Keq_ll_arr = input.get_array("Keq_ll"); + const auto electrolyte_arr = input.get_array("electrolyte"); + const auto log_gam_arr = input.get_array("log_gam"); + const auto log_gamZ_arr = input.get_array("log_gamZ"); + auto jaerosolstate_arr = input.get_array("jaerosolstate"); + auto jphase_arr = input.get_array("jphase"); + auto jhyst_leg_arr = input.get_array("jhyst_leg"); + const auto aH2O_a_arr = input.get_array("aH2O_a"); + + const auto mmd = TChem::Impl::MosaicModelData(); + + real_type_1d_view aer("aer", mmd.naer); + verification::convert_1d_vector_to_1d_view_device(aer_arr, aer); + + real_type_1d_view ma("ma", mmd.nelectrolyte); + verification::convert_1d_vector_to_1d_view_device(ma_arr, ma); + + real_type_1d_view mc("mc", mmd.nelectrolyte); + verification::convert_1d_vector_to_1d_view_device(mc_arr, mc); + + real_type_1d_view Keq_ll("Keq_ll", mmd.nrxn_aer_ll); + verification::convert_1d_vector_to_1d_view_device(Keq_ll_arr, Keq_ll); + + real_type_1d_view electrolyte_solid("electrolyte_solid", mmd.nelectrolyte); + real_type_1d_view electrolyte_liquid("electrolyte_liquid", mmd.nelectrolyte); + real_type_1d_view electrolyte_total("electrolyte_total", mmd.nelectrolyte); + + const ordinal_type n = 3; + const auto nsize_electrolyte = static_cast(electrolyte_arr.size())/n; + std::vector> electrolyte_arr_2d; + for (size_t i = 0; i < n*nsize_electrolyte; i += nsize_electrolyte) { + // Create a small vector from a segment of long_vector + std::vector electrolyte_temp(electrolyte_arr.begin() + i, electrolyte_arr.begin() + i + nsize_electrolyte); + electrolyte_arr_2d.push_back(electrolyte_temp); + } + + verification::convert_1d_vector_to_1d_view_device(electrolyte_arr_2d[0], electrolyte_solid); + verification::convert_1d_vector_to_1d_view_device(electrolyte_arr_2d[1], electrolyte_liquid); + verification::convert_1d_vector_to_1d_view_device(electrolyte_arr_2d[2], electrolyte_total); + + real_type_1d_view log_gam("log_gam", mmd.nelectrolyte); + verification::convert_1d_vector_to_1d_view_device(log_gam_arr, log_gam); + + real_type_2d_view log_gamZ("log_gamZ", mmd.nelectrolyte, mmd.nelectrolyte); + verification::convert_1d_vector_to_2d_view_device(log_gamZ_arr, log_gamZ); + + real_type_1d_view jaerosolstate("jaerosolstate", 1); + verification::convert_1d_vector_to_1d_view_device(jaerosolstate_arr, jaerosolstate); + + real_type_1d_view jphase("jphase", 1); + verification::convert_1d_vector_to_1d_view_device(jphase_arr, jphase); + + real_type_1d_view jhyst_leg("jhyst_leg", 1); + verification::convert_1d_vector_to_1d_view_device(jhyst_leg_arr, jhyst_leg); + + real_type_1d_view aH2O_a("aH2O_a", 1); + verification::convert_1d_vector_to_1d_view_device(aH2O_a_arr, aH2O_a); + + real_type_1d_view molalities("molalities", mmd.nelectrolyte); + real_type_1d_view xmol("xmol", mmd.nelectrolyte); + real_type_1d_view gam("gam", mmd.nelectrolyte); + real_type_1d_view activity("activity", mmd.nelectrolyte); + + std::string profile_name ="Verification_test_compute_activites"; + using policy_type = + typename TChem::UseThisTeamPolicy::type; + const auto exec_space_instance = TChem::exec_space(); + const auto host_exec_space = TChem::host_exec_space(); + policy_type policy(exec_space_instance, 1, Kokkos::AUTO()); + + // Check this routines on GPUs. + Kokkos::parallel_for( + profile_name, + policy, + KOKKOS_LAMBDA(const typename policy_type::member_type& member) { + + // Perform the adjustment calculation + TChem::Impl::MOSAIC::compute_activities( + mmd, + molalities, + xmol, + aer, + ma, + mc, + Keq_ll, + electrolyte_solid, + electrolyte_liquid, + electrolyte_total, + log_gam, + log_gamZ, + gam, + activity, + jaerosolstate(0), + jphase(0), + jhyst_leg(0), + aH2O_a(0)); + }); + + verification::convert_1d_view_device_to_1d_vector(jaerosolstate, jaerosolstate_arr); + output.set("jaerosolstate", jaerosolstate_arr); + + verification::convert_1d_view_device_to_1d_vector(jhyst_leg, jhyst_leg_arr); + output.set("jhyst_leg", jhyst_leg_arr); + + verification::convert_1d_view_device_to_1d_vector(jphase, jphase_arr); + output.set("jphase", jphase_arr); + + std::vector molalities_arr(mmd.nelectrolyte); + verification::convert_1d_view_device_to_1d_vector(molalities, molalities_arr); + output.set("molalities", molalities_arr); + + std::vector activity_arr(mmd.nelectrolyte); + verification::convert_1d_view_device_to_1d_vector(activity, activity_arr); + output.set("activity", activity_arr); + }); +} \ No newline at end of file diff --git a/src/verification/mosaic/mosaic_driver.cpp b/src/verification/mosaic/mosaic_driver.cpp index b42e38d2..eeaa21db 100644 --- a/src/verification/mosaic/mosaic_driver.cpp +++ b/src/verification/mosaic/mosaic_driver.cpp @@ -19,6 +19,8 @@ void check_aerosol_mass(Ensemble *ensemble); void adjust_liquid_aerosol(Ensemble *ensemble); +void compute_activities(Ensemble *ensemble); + void adjust_solid_aerosol(Ensemble *ensemble); void do_full_deliquescence(Ensemble *ensemble); @@ -72,6 +74,8 @@ int main(int argc, char **argv) { check_aerosol_mass(ensemble); } else if (func_name == "adjust_liquid_aerosol") { adjust_liquid_aerosol(ensemble); + } else if (func_name == "compute_activities") { + compute_activities(ensemble); } else if (func_name == "adjust_solid_aerosol") { adjust_solid_aerosol(ensemble); } else if (func_name == "do_full_deliquescence") {