Skip to content

Chemical mask implementation #504

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
SreejithNREL wants to merge 8 commits into
base: development
Choose a base branch
from
Open

Chemical mask implementation #504

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
d1f5961
Modified RK64 to include chem reaction mask
SreejithNREL Apr 20, 2024
ea68b74
Added chem mask for Arkode
SreejithNREL May 24, 2024
db216d1
Corrected Cvode
SreejithNREL May 24, 2024
e742d2f
Clang tidy ied
SreejithNREL May 27, 2024
fa66c36
Update ReactorCvodeUtils.H
SreejithNREL May 27, 2024
c748573
Update ReactorCvodeUtils.H
SreejithNREL May 27, 2024
377232a
Update ReactorArkode.cpp
SreejithNREL May 27, 2024
2c9f963
Update ReactorCvode.cpp
SreejithNREL May 27, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
19 changes: 19 additions & 0 deletions Source/Reactions/ReactorArkode.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -128,6 +128,25 @@ public:
rhoE, frcEExt, FC_in, y_vect, vect_energy, FCunt, dt);
}

void unflatten(
const amrex::Box& box,
const int ncells,
amrex::Array4<amrex::Real> const& rhoY,
amrex::Array4<amrex::Real> const& temperature,
amrex::Array4<amrex::Real> const& rhoE,
amrex::Array4<amrex::Real> const& frcEExt,
amrex::Array4<amrex::Real> const& FC_in,
amrex::Array4<int> const& mask,
Copy link
Contributor

@marchdf marchdf May 28, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

the mask should probably come in as amrex::Array4<const int> const&

amrex::Real* y_vect,
amrex::Real* vect_energy,
long int* FCunt,
amrex::Real dt)
{
flatten_ops.unflatten(
box, ncells, m_reactor_type, m_clean_init_massfrac, rhoY, temperature,
rhoE, frcEExt, FC_in, mask, y_vect, vect_energy, FCunt, dt);
}

private:
amrex::Real relTol{1e-6};
amrex::Real absTol{1e-10};
Expand Down
4 changes: 2 additions & 2 deletions Source/Reactions/ReactorArkode.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -135,7 +135,7 @@ ReactorArkode::react(
amrex::Array4<amrex::Real> const& rEner_in,
amrex::Array4<amrex::Real> const& rEner_src_in,
amrex::Array4<amrex::Real> const& FC_in,
amrex::Array4<int> const& /*mask*/,
amrex::Array4<int> const& mask,
Copy link
Contributor

@marchdf marchdf May 28, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same here.

amrex::Real& dt_react,
amrex::Real& time
#ifdef AMREX_USE_GPU
Expand Down Expand Up @@ -239,7 +239,7 @@ ReactorArkode::react(
amrex::Gpu::DeviceVector<long int> v_nfe(ncells, nfe);
long int* d_nfe = v_nfe.data();
unflatten(
box, ncells, rY_in, T_in, rEner_in, rEner_src_in, FC_in, yvec_d,
box, ncells, rY_in, T_in, rEner_in, rEner_src_in, FC_in, mask, yvec_d,
user_data->rhoe_init, d_nfe, dt_react);

N_VDestroy(y);
Expand Down
217 changes: 110 additions & 107 deletions Source/Reactions/ReactorBDF.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -319,7 +319,7 @@ ReactorBDF::react(
amrex::Array4<amrex::Real> const& rEner_in,
amrex::Array4<amrex::Real> const& rEner_src_in,
amrex::Array4<amrex::Real> const& FC_in,
amrex::Array4<int> const& /*mask*/,
amrex::Array4<int> const& mask,
amrex::Real& dt_react,
amrex::Real& time
#ifdef AMREX_USE_GPU
Expand Down Expand Up @@ -355,127 +355,130 @@ ReactorBDF::react(
int* d_nsteps = v_nsteps.data();

amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
amrex::Real soln[NUM_SPECIES + 1] = {0.0};
amrex::Real soln_n[NUM_SPECIES + 1] = {0.0}; // at time level n
amrex::Real soln_nm1[NUM_SPECIES + 1] = {0.0}; // at time level n-1
amrex::Real soln_nm2[NUM_SPECIES + 1] = {0.0}; // at time level n-2
amrex::Real dsoln[NUM_SPECIES + 1] = {
0.0}; // newton_soln_k+1 - newton_soln_k
amrex::Real dsoln0[NUM_SPECIES + 1] = {
0.0}; // initial newton_soln_k+1 -newton_soln_k
amrex::Real rYsrc_ext[NUM_SPECIES] = {0.0};
const int neq = (NUM_SPECIES + 1);

// initialization of variables before timestepping
amrex::Real current_time = time_init;
amrex::Real dt = dt_react / amrex::Real(captured_nsubsteps);
amrex::Real rho = 0.0;
amrex::Real massfrac[NUM_SPECIES] = {0.0};
for (int sp = 0; sp < NUM_SPECIES; sp++) {
soln_n[sp] = rY_in(i, j, k, sp);
}
get_rho_and_massfracs(soln_n, rho, massfrac);

amrex::Real temp = T_in(i, j, k, 0);
amrex::Real Enrg_loc = rEner_in(i, j, k, 0) / rho;
auto eos = pele::physics::PhysicsType::eos();
if (captured_reactor_type == ReactorTypes::e_reactor_type) {
eos.REY2T(rho, Enrg_loc, massfrac, temp);
} else if (captured_reactor_type == ReactorTypes::h_reactor_type) {
eos.RHY2T(rho, Enrg_loc, massfrac, temp);
} else {
amrex::Abort("Wrong reactor type. Choose between 1 (e) or 2 (h).");
}
soln_n[NUM_SPECIES] = temp;
amrex::Real rhoe_init[] = {rEner_in(i, j, k, 0)};
amrex::Real rhoesrc_ext[] = {rEner_src_in(i, j, k, 0)};
for (int sp = 0; sp < NUM_SPECIES; sp++) {
rYsrc_ext[sp] = rYsrc_in(i, j, k, sp);
}
for (int ii = 0; ii < neq; ii++) {
soln[ii] = soln_n[ii];
soln_nm1[ii] = soln_n[ii];
soln_nm2[ii] = soln_n[ii];
}

// begin timestepping
int printflag = 0;
// if BDF2 use trapz or BDF1 for first step
int first_tstepscheme =
(captured_tstepscheme >= BDF2SCHEME) ? TRPZSCHEME : captured_tstepscheme;
int schemechangestep = captured_tstepscheme - 2; // 0 for BDF2 and 1 for
// BDF3
amrex::Real rhs[(NUM_SPECIES + 1)] = {0.0};
amrex::Real Jmat2d[NUM_SPECIES + 1][NUM_SPECIES + 1] = {{0.0}};
for (int nsteps = 0; nsteps < captured_nsubsteps; nsteps++) {
// shift to BDF2 after first step
int tstepscheme =
(nsteps > schemechangestep) ? captured_tstepscheme : first_tstepscheme;

if (mask(i, j, k) != -1) {
amrex::Real soln[NUM_SPECIES + 1] = {0.0};
amrex::Real soln_n[NUM_SPECIES + 1] = {0.0}; // at time level n
amrex::Real soln_nm1[NUM_SPECIES + 1] = {0.0}; // at time level n-1
amrex::Real soln_nm2[NUM_SPECIES + 1] = {0.0}; // at time level n-2
amrex::Real dsoln[NUM_SPECIES + 1] = {
0.0}; // newton_soln_k+1 - newton_soln_k
amrex::Real dsoln0[NUM_SPECIES + 1] = {
0.0}; // initial newton_soln_k+1 -newton_soln_k
amrex::Real rYsrc_ext[NUM_SPECIES] = {0.0};
const int neq = (NUM_SPECIES + 1);

// initialization of variables before timestepping
amrex::Real current_time = time_init;
amrex::Real dt = dt_react / amrex::Real(captured_nsubsteps);
amrex::Real rho = 0.0;
amrex::Real massfrac[NUM_SPECIES] = {0.0};
for (int sp = 0; sp < NUM_SPECIES; sp++) {
soln_n[sp] = rY_in(i, j, k, sp);
}
get_rho_and_massfracs(soln_n, rho, massfrac);

amrex::Real temp = T_in(i, j, k, 0);
amrex::Real Enrg_loc = rEner_in(i, j, k, 0) / rho;
auto eos = pele::physics::PhysicsType::eos();
if (captured_reactor_type == ReactorTypes::e_reactor_type) {
eos.REY2T(rho, Enrg_loc, massfrac, temp);
} else if (captured_reactor_type == ReactorTypes::h_reactor_type) {
eos.RHY2T(rho, Enrg_loc, massfrac, temp);
} else {
amrex::Abort("Wrong reactor type. Choose between 1 (e) or 2 (h).");
}
soln_n[NUM_SPECIES] = temp;
amrex::Real rhoe_init[] = {rEner_in(i, j, k, 0)};
amrex::Real rhoesrc_ext[] = {rEner_src_in(i, j, k, 0)};
for (int sp = 0; sp < NUM_SPECIES; sp++) {
rYsrc_ext[sp] = rYsrc_in(i, j, k, sp);
}
for (int ii = 0; ii < neq; ii++) {
dsoln0[ii] = 0.0;
soln[ii] = soln_n[ii];
soln_nm1[ii] = soln_n[ii];
soln_nm2[ii] = soln_n[ii];
}
// non-linear iterations for each timestep
for (int nlit = 0; nlit < captured_nonlinear_iters; nlit++) {
for (int ii = 0; ii < neq; ii++) {
dsoln0[ii] = dsoln[ii];
}
get_bdf_matrix_and_rhs(
soln, soln_n, soln_nm1, soln_nm2, captured_reactor_type, tstepscheme,
dt, rhoe_init, rhoesrc_ext, rYsrc_ext, current_time, time_init,
Jmat2d, rhs);

performgmres(
Jmat2d, rhs, dsoln0, dsoln, captured_gmres_precond,
captured_gmres_restarts, captured_gmres_kspiters, captured_gmres_tol,
printflag);
// begin timestepping
int printflag = 0;
// if BDF2 use trapz or BDF1 for first step
int first_tstepscheme = (captured_tstepscheme >= BDF2SCHEME)
? TRPZSCHEME
: captured_tstepscheme;
int schemechangestep = captured_tstepscheme - 2; // 0 for BDF2 and 1 for
// BDF3
amrex::Real rhs[(NUM_SPECIES + 1)] = {0.0};
amrex::Real Jmat2d[NUM_SPECIES + 1][NUM_SPECIES + 1] = {{0.0}};
for (int nsteps = 0; nsteps < captured_nsubsteps; nsteps++) {
// shift to BDF2 after first step
int tstepscheme = (nsteps > schemechangestep) ? captured_tstepscheme
: first_tstepscheme;

for (int ii = 0; ii < neq; ii++) {
soln[ii] += dsoln[ii];
dsoln0[ii] = 0.0;
}
// non-linear iterations for each timestep
for (int nlit = 0; nlit < captured_nonlinear_iters; nlit++) {
for (int ii = 0; ii < neq; ii++) {
dsoln0[ii] = dsoln[ii];
}
get_bdf_matrix_and_rhs(
soln, soln_n, soln_nm1, soln_nm2, captured_reactor_type,
tstepscheme, dt, rhoe_init, rhoesrc_ext, rYsrc_ext, current_time,
time_init, Jmat2d, rhs);

performgmres(
Jmat2d, rhs, dsoln0, dsoln, captured_gmres_precond,
captured_gmres_restarts, captured_gmres_kspiters,
captured_gmres_tol, printflag);

for (int ii = 0; ii < neq; ii++) {
soln[ii] += dsoln[ii];
}

amrex::Real norm = 0.0;
for (int ii = 0; ii < neq; ii++) {
norm += rhs[ii] * rhs[ii];
}
norm = std::sqrt(norm);
if (norm <= captured_nonlin_tol) {
break;
}
}

amrex::Real norm = 0.0;
// copy non-linear solution onto soln_n,
// soln_n to soln_nm1
for (int ii = 0; ii < neq; ii++) {
norm += rhs[ii] * rhs[ii];
}
norm = std::sqrt(norm);
if (norm <= captured_nonlin_tol) {
break;
soln_nm2[ii] = soln_nm1[ii];
soln_nm1[ii] = soln_n[ii];
soln_n[ii] = soln[ii];
}
current_time += dt;
}

// copy non-linear solution onto soln_n,
// soln_n to soln_nm1
for (int ii = 0; ii < neq; ii++) {
soln_nm2[ii] = soln_nm1[ii];
soln_nm1[ii] = soln_n[ii];
soln_n[ii] = soln[ii];
}
current_time += dt;
}

// copy data back
int icell = (k - lo.z) * len.x * len.y + (j - lo.y) * len.x + (i - lo.x);
d_nsteps[icell] = captured_nsubsteps;
// copy data back
int icell = (k - lo.z) * len.x * len.y + (j - lo.y) * len.x + (i - lo.x);
d_nsteps[icell] = captured_nsubsteps;

get_rho_and_massfracs(soln_n, rho, massfrac);
for (int sp = 0; sp < NUM_SPECIES; sp++) {
rY_in(i, j, k, sp) = soln_n[sp];
}
get_rho_and_massfracs(soln_n, rho, massfrac);
for (int sp = 0; sp < NUM_SPECIES; sp++) {
rY_in(i, j, k, sp) = soln_n[sp];
}

temp = soln_n[NUM_SPECIES];
rEner_in(i, j, k, 0) = rhoe_init[0] + dt_react * rhoesrc_ext[0];
Enrg_loc = rEner_in(i, j, k, 0) / rho;
temp = soln_n[NUM_SPECIES];
rEner_in(i, j, k, 0) = rhoe_init[0] + dt_react * rhoesrc_ext[0];
Enrg_loc = rEner_in(i, j, k, 0) / rho;

if (captured_reactor_type == ReactorTypes::e_reactor_type) {
eos.REY2T(rho, Enrg_loc, massfrac, temp);
} else if (captured_reactor_type == ReactorTypes::h_reactor_type) {
eos.RHY2T(rho, Enrg_loc, massfrac, temp);
} else {
amrex::Abort("Wrong reactor type. Choose between 1 (e) or 2 (h).");
if (captured_reactor_type == ReactorTypes::e_reactor_type) {
eos.REY2T(rho, Enrg_loc, massfrac, temp);
} else if (captured_reactor_type == ReactorTypes::h_reactor_type) {
eos.RHY2T(rho, Enrg_loc, massfrac, temp);
} else {
amrex::Abort("Wrong reactor type. Choose between 1 (e) or 2 (h).");
}
T_in(i, j, k, 0) = temp;
FC_in(i, j, k, 0) = captured_nsubsteps;
}
T_in(i, j, k, 0) = temp;
FC_in(i, j, k, 0) = captured_nsubsteps;
});

#ifdef MOD_REACTOR
Expand Down
19 changes: 19 additions & 0 deletions Source/Reactions/ReactorCvode.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -146,6 +146,25 @@ public:
rhoE, frcEExt, FC_in, y_vect, vect_energy, FCunt, dt);
}

void unflatten(
const amrex::Box& box,
const int ncells,
amrex::Array4<amrex::Real> const& rhoY,
amrex::Array4<amrex::Real> const& temperature,
amrex::Array4<amrex::Real> const& rhoE,
amrex::Array4<amrex::Real> const& frcEExt,
amrex::Array4<amrex::Real> const& FC_in,
amrex::Array4<int> const& mask,
amrex::Real* y_vect,
amrex::Real* vect_energy,
long int* FCunt,
amrex::Real dt)
{
flatten_ops.unflatten(
box, ncells, m_reactor_type, m_clean_init_massfrac, rhoY, temperature,
rhoE, frcEExt, FC_in, mask, y_vect, vect_energy, FCunt, dt);
}

private:
void checkCvodeOptions(
const std::string& a_solve_type_str,
Expand Down
2 changes: 1 addition & 1 deletion Source/Reactions/ReactorCvode.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1373,7 +1373,7 @@ ReactorCvode::react(
amrex::Gpu::DeviceVector<long int> v_nfe(ncells, nfe);
long int* d_nfe = v_nfe.data();
unflatten(
box, ncells, rY_in, T_in, rEner_in, rEner_src_in, FC_in, yvec_d,
box, ncells, rY_in, T_in, rEner_in, rEner_src_in, FC_in, mask, yvec_d,
udata->rhoe_init, d_nfe, dt_react);

if (udata->verbose > 1) {
Expand Down
Loading