Skip to content

Commit

Permalink
Clip the concentrations on fractional pows (#433)
Browse files Browse the repository at this point in the history
  • Loading branch information
@marchdf
marchdf authored Oct 25, 2023
1 parent 71a7780 commit b1b1892
Show file tree
Hide file tree
Showing 4 changed files with 131 additions and 51 deletions.
57 changes: 33 additions & 24 deletions Support/Mechanism/Models/JL4/mechanism.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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;
Expand All @@ -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;
Expand All @@ -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;
Expand Down Expand Up @@ -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
Expand All @@ -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]
Expand All @@ -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 =
Expand Down Expand Up @@ -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
Expand All @@ -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]
Expand Down Expand Up @@ -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
Expand All @@ -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]
Expand All @@ -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 =
Expand Down Expand Up @@ -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
Expand All @@ -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]
Expand Down
35 changes: 20 additions & 15 deletions Support/Mechanism/Models/chem-CH4-2step/mechanism.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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
Expand All @@ -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]
Expand Down Expand Up @@ -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
Expand All @@ -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]
Expand Down
24 changes: 20 additions & 4 deletions Support/ceptr/ceptr/jacobian.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand All @@ -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
Expand Down
Loading

0 comments on commit b1b1892

Please sign in to comment.