From b1b18924e8090b8fcf3861b8d51060953ec0ec53 Mon Sep 17 00:00:00 2001 From: "Marc T. Henry de Frahan" Date: Wed, 25 Oct 2023 08:39:18 -0600 Subject: [PATCH] Clip the concentrations on fractional pows (#433) --- Support/Mechanism/Models/JL4/mechanism.H | 57 +++++++++------- .../Models/chem-CH4-2step/mechanism.H | 35 +++++----- Support/ceptr/ceptr/jacobian.py | 24 +++++-- Support/ceptr/ceptr/utilities.py | 66 ++++++++++++++++--- 4 files changed, 131 insertions(+), 51 deletions(-) diff --git a/Support/Mechanism/Models/JL4/mechanism.H b/Support/Mechanism/Models/JL4/mechanism.H index 698d23b6d..5bc2fee8e 100644 --- a/Support/Mechanism/Models/JL4/mechanism.H +++ b/Support/Mechanism/Models/JL4/mechanism.H @@ -1913,15 +1913,16 @@ comp_qfqr( { // reaction 0: 2 CH4 + O2 => 2 CO + 4 H2 - qf[0] = pow(sc[0], 0.500000) * pow(sc[1], 1.250000); + qf[0] = std::sqrt(std::max(sc[0], 0.0)) * pow(std::max(sc[1], 0.0), 1.250000); qr[0] = 0.0; // reaction 1: CH4 + H2O <=> CO + 3 H2 qf[1] = sc[0] * sc[2]; - qr[1] = sc[4] * pow(sc[6], 3.000000); + qr[1] = sc[4] * sc[6] * sc[6] * sc[6]; // reaction 2: 2 H2 + O2 => 2 H2O - qf[2] = pow(sc[1], 1.500000) * pow(sc[6], 0.250000); + qf[2] = + pow(std::max(sc[1], 0.0), 1.500000) * pow(std::max(sc[6], 0.0), 0.250000); qr[2] = 0.0; // reaction 3: CO + H2O <=> CO2 + H2 @@ -1997,7 +1998,8 @@ productionRate(amrex::Real* wdot, const amrex::Real* sc, const amrex::Real T) { // reaction 0: 2 CH4 + O2 => 2 CO + 4 H2 const amrex::Real k_f = 1236450565.12584 * exp(-(15096.4999741416) * invT); - const amrex::Real qf = k_f * (pow(sc[0], 0.500000) * pow(sc[1], 1.250000)); + const amrex::Real qf = k_f * (std::sqrt(std::max(sc[0], 0.0)) * + pow(std::max(sc[1], 0.0), 1.250000)); const amrex::Real qr = 0.0; const amrex::Real qdot = qf - qr; wdot[0] -= 2.000000 * qdot; @@ -2012,7 +2014,7 @@ productionRate(amrex::Real* wdot, const amrex::Real* sc, const amrex::Real T) const amrex::Real qf = k_f * (sc[0] * sc[2]); const amrex::Real qr = k_f * exp(-(g_RT[0] + g_RT[2] - g_RT[4] - 3.000000 * g_RT[6])) * - ((refCinv * refCinv)) * (sc[4] * pow(sc[6], 3.000000)); + ((refCinv * refCinv)) * (sc[4] * sc[6] * sc[6] * sc[6]); const amrex::Real qdot = qf - qr; wdot[0] -= qdot; wdot[2] -= qdot; @@ -2024,7 +2026,8 @@ productionRate(amrex::Real* wdot, const amrex::Real* sc, const amrex::Real T) // reaction 2: 2 H2 + O2 => 2 H2O const amrex::Real k_f = 191159684557179 * exp((-1) * tc[0] - (20128.6666321888) * invT); - const amrex::Real qf = k_f * (pow(sc[1], 1.500000) * pow(sc[6], 0.250000)); + const amrex::Real qf = k_f * (pow(std::max(sc[1], 0.0), 1.500000) * + pow(std::max(sc[6], 0.0), 0.250000)); const amrex::Real qr = 0.0; const amrex::Real qdot = qf - qr; wdot[1] -= qdot; @@ -2324,7 +2327,7 @@ aJacobian_precond( // reaction 0: 2 CH4 + O2 => 2 CO + 4 H2 // a non-third-body and non-pressure-fall-off reaction // forward - phi_f = pow(sc[0], 0.500000) * pow(sc[1], 1.250000); + phi_f = std::sqrt(std::max(sc[0], 0.0)) * pow(std::max(sc[1], 0.0), 1.250000); k_f = 1236450565.12584 * exp(-(15096.4999741416) * invT); dlnkfdT = (15096.4999741416) * invT2; // rate of progress @@ -2336,14 +2339,15 @@ aJacobian_precond( wdot[4] += 2 * q; // CO wdot[6] += 4 * q; // H2 // d()/d[CH4] - dqdci = +k_f * 0.500000 * pow(std::max(sc[0], 1e-12), -0.500000) * - pow(sc[1], 1.250000); + dqdci = +k_f * 0.500000 * pow(std::max(sc[0], 1e-16), -0.500000) * + pow(std::max(sc[1], 0.0), 1.250000); J[0] += -2 * dqdci; // dwdot[CH4]/d[CH4] J[1] -= dqdci; // dwdot[O2]/d[CH4] J[4] += 2 * dqdci; // dwdot[CO]/d[CH4] J[6] += 4 * dqdci; // dwdot[H2]/d[CH4] // d()/d[O2] - dqdci = +k_f * pow(sc[0], 0.500000) * 1.250000 * pow(sc[1], 0.250000); + dqdci = +k_f * std::sqrt(std::max(sc[0], 0.0)) * 1.250000 * + pow(std::max(sc[1], 0.0), 0.250000); J[8] += -2 * dqdci; // dwdot[CH4]/d[O2] J[9] -= dqdci; // dwdot[O2]/d[O2] J[12] += 2 * dqdci; // dwdot[CO]/d[O2] @@ -2361,7 +2365,7 @@ aJacobian_precond( k_f = 300000 * exp(-(15096.4999741416) * invT); dlnkfdT = (15096.4999741416) * invT2; // reverse - phi_r = sc[4] * pow(sc[6], 3.000000); + phi_r = sc[4] * sc[6] * sc[6] * sc[6]; Kc = (refC * refC) * exp(g_RT[0] + g_RT[2] - g_RT[4] - 3.000000 * g_RT[6]); k_r = k_f / Kc; dlnKcdT = @@ -2408,7 +2412,8 @@ aJacobian_precond( // reaction 2: 2 H2 + O2 => 2 H2O // a non-third-body and non-pressure-fall-off reaction // forward - phi_f = pow(sc[1], 1.500000) * pow(sc[6], 0.250000); + phi_f = + pow(std::max(sc[1], 0.0), 1.500000) * pow(std::max(sc[6], 0.0), 0.250000); k_f = 191159684557179 * exp(-1 * tc[0] - (20128.6666321888) * invT); dlnkfdT = -1 * invT + (20128.6666321888) * invT2; // rate of progress @@ -2419,13 +2424,14 @@ aJacobian_precond( wdot[2] += 2 * q; // H2O wdot[6] -= 2 * q; // H2 // d()/d[O2] - dqdci = +k_f * 1.500000 * pow(sc[1], 0.500000) * pow(sc[6], 0.250000); + dqdci = +k_f * 1.500000 * std::sqrt(std::max(sc[1], 0.0)) * + pow(std::max(sc[6], 0.0), 0.250000); J[9] -= dqdci; // dwdot[O2]/d[O2] J[10] += 2 * dqdci; // dwdot[H2O]/d[O2] J[14] += -2 * dqdci; // dwdot[H2]/d[O2] // d()/d[H2] - dqdci = +k_f * pow(sc[1], 1.500000) * 0.250000 * - pow(std::max(sc[6], 1e-12), -0.750000); + dqdci = +k_f * pow(std::max(sc[1], 0.0), 1.500000) * 0.250000 * + pow(std::max(sc[6], 1e-16), -0.750000); J[49] -= dqdci; // dwdot[O2]/d[H2] J[50] += 2 * dqdci; // dwdot[H2O]/d[H2] J[54] += -2 * dqdci; // dwdot[H2]/d[H2] @@ -2598,7 +2604,7 @@ aJacobian( // reaction 0: 2 CH4 + O2 => 2 CO + 4 H2 // a non-third-body and non-pressure-fall-off reaction // forward - phi_f = pow(sc[0], 0.500000) * pow(sc[1], 1.250000); + phi_f = std::sqrt(std::max(sc[0], 0.0)) * pow(std::max(sc[1], 0.0), 1.250000); k_f = 1236450565.12584 * exp(-(15096.4999741416) * invT); dlnkfdT = (15096.4999741416) * invT2; // rate of progress @@ -2610,14 +2616,15 @@ aJacobian( wdot[4] += 2 * q; // CO wdot[6] += 4 * q; // H2 // d()/d[CH4] - dqdci = +k_f * 0.500000 * pow(std::max(sc[0], 1e-12), -0.500000) * - pow(sc[1], 1.250000); + dqdci = +k_f * 0.500000 * pow(std::max(sc[0], 1e-16), -0.500000) * + pow(std::max(sc[1], 0.0), 1.250000); J[0] += -2 * dqdci; // dwdot[CH4]/d[CH4] J[1] -= dqdci; // dwdot[O2]/d[CH4] J[4] += 2 * dqdci; // dwdot[CO]/d[CH4] J[6] += 4 * dqdci; // dwdot[H2]/d[CH4] // d()/d[O2] - dqdci = +k_f * pow(sc[0], 0.500000) * 1.250000 * pow(sc[1], 0.250000); + dqdci = +k_f * std::sqrt(std::max(sc[0], 0.0)) * 1.250000 * + pow(std::max(sc[1], 0.0), 0.250000); J[8] += -2 * dqdci; // dwdot[CH4]/d[O2] J[9] -= dqdci; // dwdot[O2]/d[O2] J[12] += 2 * dqdci; // dwdot[CO]/d[O2] @@ -2635,7 +2642,7 @@ aJacobian( k_f = 300000 * exp(-(15096.4999741416) * invT); dlnkfdT = (15096.4999741416) * invT2; // reverse - phi_r = sc[4] * pow(sc[6], 3.000000); + phi_r = sc[4] * sc[6] * sc[6] * sc[6]; Kc = (refC * refC) * exp(g_RT[0] + g_RT[2] - g_RT[4] - 3.000000 * g_RT[6]); k_r = k_f / Kc; dlnKcdT = @@ -2682,7 +2689,8 @@ aJacobian( // reaction 2: 2 H2 + O2 => 2 H2O // a non-third-body and non-pressure-fall-off reaction // forward - phi_f = pow(sc[1], 1.500000) * pow(sc[6], 0.250000); + phi_f = + pow(std::max(sc[1], 0.0), 1.500000) * pow(std::max(sc[6], 0.0), 0.250000); k_f = 191159684557179 * exp(-1 * tc[0] - (20128.6666321888) * invT); dlnkfdT = -1 * invT + (20128.6666321888) * invT2; // rate of progress @@ -2693,13 +2701,14 @@ aJacobian( wdot[2] += 2 * q; // H2O wdot[6] -= 2 * q; // H2 // d()/d[O2] - dqdci = +k_f * 1.500000 * pow(sc[1], 0.500000) * pow(sc[6], 0.250000); + dqdci = +k_f * 1.500000 * std::sqrt(std::max(sc[1], 0.0)) * + pow(std::max(sc[6], 0.0), 0.250000); J[9] -= dqdci; // dwdot[O2]/d[O2] J[10] += 2 * dqdci; // dwdot[H2O]/d[O2] J[14] += -2 * dqdci; // dwdot[H2]/d[O2] // d()/d[H2] - dqdci = +k_f * pow(sc[1], 1.500000) * 0.250000 * - pow(std::max(sc[6], 1e-12), -0.750000); + dqdci = +k_f * pow(std::max(sc[1], 0.0), 1.500000) * 0.250000 * + pow(std::max(sc[6], 1e-16), -0.750000); J[49] -= dqdci; // dwdot[O2]/d[H2] J[50] += 2 * dqdci; // dwdot[H2O]/d[H2] J[54] += -2 * dqdci; // dwdot[H2]/d[H2] diff --git a/Support/Mechanism/Models/chem-CH4-2step/mechanism.H b/Support/Mechanism/Models/chem-CH4-2step/mechanism.H index a337258c4..8f9171b96 100644 --- a/Support/Mechanism/Models/chem-CH4-2step/mechanism.H +++ b/Support/Mechanism/Models/chem-CH4-2step/mechanism.H @@ -1787,7 +1787,8 @@ comp_qfqr( qr[0] = 0.0; // reaction 1: 2 CO + O2 => 2 CO2 - qf[1] = pow(sc[0], 0.250000) * pow(sc[1], 0.500000) * sc[3]; + qf[1] = pow(std::max(sc[0], 0.0), 0.250000) * + std::sqrt(std::max(sc[1], 0.0)) * sc[3]; qr[1] = 0.0; // reaction 2: 2 CO2 => 2 CO + O2 @@ -1872,8 +1873,8 @@ productionRate(amrex::Real* wdot, const amrex::Real* sc, const amrex::Real T) { // reaction 1: 2 CO + O2 => 2 CO2 const amrex::Real k_f = 6296094821.39524 * exp(-(20128.6666321888) * invT); - const amrex::Real qf = - k_f * (pow(sc[0], 0.250000) * pow(sc[1], 0.500000) * sc[3]); + const amrex::Real qf = k_f * (pow(std::max(sc[0], 0.0), 0.250000) * + std::sqrt(std::max(sc[1], 0.0)) * sc[3]); const amrex::Real qr = 0.0; const amrex::Real qdot = qf - qr; wdot[0] -= qdot; @@ -2195,7 +2196,8 @@ aJacobian_precond( // reaction 1: 2 CO + O2 => 2 CO2 // a non-third-body and non-pressure-fall-off reaction // forward - phi_f = pow(sc[0], 0.250000) * pow(sc[1], 0.500000) * sc[3]; + phi_f = pow(std::max(sc[0], 0.0), 0.250000) * + std::sqrt(std::max(sc[1], 0.0)) * sc[3]; k_f = 6296094821.39524 * exp(-(20128.6666321888) * invT); dlnkfdT = (20128.6666321888) * invT2; // rate of progress @@ -2206,19 +2208,20 @@ aJacobian_precond( wdot[3] -= 2 * q; // CO wdot[4] += 2 * q; // CO2 // d()/d[O2] - dqdci = +k_f * 0.250000 * pow(std::max(sc[0], 1e-12), -0.750000) * - pow(sc[1], 0.500000) * sc[3]; + dqdci = +k_f * 0.250000 * pow(std::max(sc[0], 1e-16), -0.750000) * + std::sqrt(std::max(sc[1], 0.0)) * sc[3]; J[0] -= dqdci; // dwdot[O2]/d[O2] J[3] += -2 * dqdci; // dwdot[CO]/d[O2] J[4] += 2 * dqdci; // dwdot[CO2]/d[O2] // d()/d[H2O] - dqdci = +k_f * pow(sc[0], 0.250000) * 0.500000 * - pow(std::max(sc[1], 1e-12), -0.500000) * sc[3]; + dqdci = +k_f * pow(std::max(sc[0], 0.0), 0.250000) * 0.500000 * + pow(std::max(sc[1], 1e-16), -0.500000) * sc[3]; J[7] -= dqdci; // dwdot[O2]/d[H2O] J[10] += -2 * dqdci; // dwdot[CO]/d[H2O] J[11] += 2 * dqdci; // dwdot[CO2]/d[H2O] // d()/d[CO] - dqdci = +k_f * pow(sc[0], 0.250000) * pow(sc[1], 0.500000); + dqdci = +k_f * pow(std::max(sc[0], 0.0), 0.250000) * + std::sqrt(std::max(sc[1], 0.0)); J[21] -= dqdci; // dwdot[O2]/d[CO] J[24] += -2 * dqdci; // dwdot[CO]/d[CO] J[25] += 2 * dqdci; // dwdot[CO2]/d[CO] @@ -2396,7 +2399,8 @@ aJacobian( // reaction 1: 2 CO + O2 => 2 CO2 // a non-third-body and non-pressure-fall-off reaction // forward - phi_f = pow(sc[0], 0.250000) * pow(sc[1], 0.500000) * sc[3]; + phi_f = pow(std::max(sc[0], 0.0), 0.250000) * + std::sqrt(std::max(sc[1], 0.0)) * sc[3]; k_f = 6296094821.39524 * exp(-(20128.6666321888) * invT); dlnkfdT = (20128.6666321888) * invT2; // rate of progress @@ -2407,19 +2411,20 @@ aJacobian( wdot[3] -= 2 * q; // CO wdot[4] += 2 * q; // CO2 // d()/d[O2] - dqdci = +k_f * 0.250000 * pow(std::max(sc[0], 1e-12), -0.750000) * - pow(sc[1], 0.500000) * sc[3]; + dqdci = +k_f * 0.250000 * pow(std::max(sc[0], 1e-16), -0.750000) * + std::sqrt(std::max(sc[1], 0.0)) * sc[3]; J[0] -= dqdci; // dwdot[O2]/d[O2] J[3] += -2 * dqdci; // dwdot[CO]/d[O2] J[4] += 2 * dqdci; // dwdot[CO2]/d[O2] // d()/d[H2O] - dqdci = +k_f * pow(sc[0], 0.250000) * 0.500000 * - pow(std::max(sc[1], 1e-12), -0.500000) * sc[3]; + dqdci = +k_f * pow(std::max(sc[0], 0.0), 0.250000) * 0.500000 * + pow(std::max(sc[1], 1e-16), -0.500000) * sc[3]; J[7] -= dqdci; // dwdot[O2]/d[H2O] J[10] += -2 * dqdci; // dwdot[CO]/d[H2O] J[11] += 2 * dqdci; // dwdot[CO2]/d[H2O] // d()/d[CO] - dqdci = +k_f * pow(sc[0], 0.250000) * pow(sc[1], 0.500000); + dqdci = +k_f * pow(std::max(sc[0], 0.0), 0.250000) * + std::sqrt(std::max(sc[1], 0.0)); J[21] -= dqdci; // dwdot[O2]/d[CO] J[24] += -2 * dqdci; // dwdot[CO]/d[CO] J[25] += 2 * dqdci; // dwdot[CO2]/d[CO] diff --git a/Support/ceptr/ceptr/jacobian.py b/Support/ceptr/ceptr/jacobian.py index a007435ba..b0c1847fa 100644 --- a/Support/ceptr/ceptr/jacobian.py +++ b/Support/ceptr/ceptr/jacobian.py @@ -1452,15 +1452,26 @@ def dphase_space(mechanism, species_info, reagents, r, reaction_orders, syms): [f"sc[{species_info.ordered_idx_map[symbol]}]"] * int(exponent) ) + elif exponent == 0.5: + idx = species_info.ordered_idx_map[symbol] + conc = ( + f"std::sqrt(std::max(sc[{idx}], {cu.sc_cutoff(0.5)}))" + ) else: idx = species_info.ordered_idx_map[symbol] - conc = f"pow(sc[{idx}],{exponent:f})" + conc = ( + f"pow(std::max(sc[{idx}]," + f" {cu.sc_cutoff(exponent)}),{exponent:f})" + ) phi += [conc] elif order < 1: phi += [f"{order:f}"] exponent = order - 1.0 idx = species_info.ordered_idx_map[symbol] - conc = f"pow(std::max(sc[{idx}], 1e-12),{exponent:f})" + conc = ( + f"pow(std::max(sc[{idx}]," + f" {cu.sc_cutoff(exponent)}),{exponent:f})" + ) phi += [conc] else: if order == 1.0: @@ -1470,10 +1481,15 @@ def dphase_space(mechanism, species_info, reagents, r, reaction_orders, syms): conc = "*".join( [f"sc[{species_info.ordered_idx_map[symbol]}]"] * int(order) ) + elif order == 0.5: + conc = ( + f"std::sqrt(std::max(sc[{species_info.ordered_idx_map[symbol]}]," + f" {cu.sc_cutoff(0.5)}))" + ) else: conc = ( - f"pow(sc[{species_info.ordered_idx_map[symbol]}]," - f" {order:f})" + f"pow(std::max(sc[{species_info.ordered_idx_map[symbol]}]," + f" {cu.sc_cutoff(order)}), {order:f})" ) phi += [conc] # Symbol is in qssa_species_list diff --git a/Support/ceptr/ceptr/utilities.py b/Support/ceptr/ceptr/utilities.py index df5e129a7..5b219e2bc 100644 --- a/Support/ceptr/ceptr/utilities.py +++ b/Support/ceptr/ceptr/utilities.py @@ -13,6 +13,14 @@ def intersection(lst1, lst2): return list(set(lst1).intersection(lst2)) +def sc_cutoff(exponent): + """Return cutoff for sc when using a fractional exponent.""" + if exponent < 0: + return "1e-16" + else: + return "0.0" + + def qss_sorted_phase_space(mechanism, species_info, reaction, reagents, syms=None): """Get string of phase space.""" record_symbolic_operations = True @@ -68,10 +76,17 @@ def qss_sorted_phase_space(mechanism, species_info, reaction, reagents, syms=Non [f"sc_qss[{species_info.ordered_idx_map[symbol] - n_species}]"] * int(order) ) + elif float(order) == 0.5: + conc = ( + "std::sqrt(std::max(sc_qss" + f"[{species_info.ordered_idx_map[symbol] - n_species}]," + f" {sc_cutoff(0.5)}))" + ) else: conc = ( - f"pow(sc_qss[{species_info.ordered_idx_map[symbol] - n_species}]," - f" {float(order):f})" + "pow(sc_qss[std::max(" + f"{species_info.ordered_idx_map[symbol] - n_species}]," + f" {sc_cutoff(order)}), {float(order):f})" ) if record_symbolic_operations: conc_smp = pow( @@ -99,10 +114,28 @@ def qss_sorted_phase_space(mechanism, species_info, reaction, reagents, syms=Non syms.sc_smp[species_info.ordered_idx_map[symbol]] * syms.sc_smp[species_info.ordered_idx_map[symbol]] ) + elif order.is_integer(): + conc = "*".join( + [f"sc[{species_info.ordered_idx_map[symbol]}]"] * int(order) + ) + if record_symbolic_operations: + conc_smp = syms.sc_smp[ + species_info.ordered_idx_map[symbol] + ] ** int(order) + elif float(order) == 0.5: + conc = ( + f"std::sqrt(std::max(sc[{species_info.ordered_idx_map[symbol]}]," + f" {sc_cutoff(0.5)}))" + ) + if record_symbolic_operations: + conc_smp = pow( + syms.sc_smp[species_info.ordered_idx_map[symbol]], + float(order), + ) else: conc = ( - f"pow(sc[{species_info.ordered_idx_map[symbol]}]," - f" {float(order):f})" + f"pow(std::max(sc[{species_info.ordered_idx_map[symbol]}]," + f" {sc_cutoff(order)}), {float(order):f})" ) if record_symbolic_operations: conc_smp = pow( @@ -183,7 +216,9 @@ def fkc_conv_inv(self, mechanism, reaction, syms=None): if record_symbolic_operations: conversion_smp *= syms.refCinv_smp * syms.refCinv_smp else: - conversion = "*".join([f"pow(refCinv, {dim:f})"]) + conversion = "*".join( + [f"pow(std::max(refCinv, {sc_cutoff(dim)}), {dim:f})"] + ) if record_symbolic_operations: conversion_smp *= syms.refCinv_smp**dim else: @@ -192,7 +227,12 @@ def fkc_conv_inv(self, mechanism, reaction, syms=None): if record_symbolic_operations: conversion_smp *= syms.refC_smp else: - conversion = "*".join([f"pow(refC, {abs(dim):f})"]) + if dim.is_integer(): + conversion = "*".join(["refC"] * int(dim)) + else: + conversion = "*".join( + [f"pow(std::max(refC, {sc_cutoff(abs(dim))}), {abs(dim):f})"] + ) if record_symbolic_operations: conversion_smp *= syms.refC_smp**dim @@ -221,12 +261,22 @@ def kc_conv(mechanism, reaction): if dim == 2.0: conversion = "*".join(["(refC * refC)"]) else: - conversion = "*".join([f"pow(refC,{dim:f})"]) + if dim.is_integer(): + conversion = "*".join(["refC"] * int(dim)) + else: + conversion = "*".join( + [f"pow(std::max(refC, {sc_cutoff(dim)}),{dim:f})"] + ) else: if dim == -1.0: conversion = "*".join(["refCinv"]) else: - conversion = "*".join([f"pow(refCinv,{abs(dim):f})"]) + if dim.is_integer(): + conversion = "*".join(["refCinv"] * int(dim)) + else: + conversion = "*".join( + [f"pow(std::max(refCinv, {sc_cutoff(abs(dim))}),{abs(dim):f})"] + ) return conversion