Skip to content

Commit

Permalink
Merge branch 'development' into Probe_Implementation
Browse files Browse the repository at this point in the history 
including probe implementation
  • Loading branch information
@SreejithNREL
SreejithNREL committed Oct 10, 2024
2 parents a79812c + 45ae8ba commit 3b718ff
Show file tree
Hide file tree
Showing 8 changed files with 67 additions and 88 deletions.
2 changes: 2 additions & 0 deletions Docs/sphinx/Utility.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,6 +77,7 @@ provided below. ::
turbinflow.low.turb_conv_vel = 5. # Velocity to move through the 3rd dimension to simulate time evolution
turbinflow.low.turb_nplane = 32 # Number of planes to read and store at a time
turbinflow.low.time_offset = 0.0 # Offset in time for reading through the 3rd dimension
turbinflow.low.verbose = 0 # verbosity level

turbinflow.high.turb_file = TurbFileHIT/TurbTEST # All same as above, but for second injection patch
turbinflow.high.dir = 1
Expand All @@ -87,6 +88,7 @@ provided below. ::
turbinflow.high.turb_conv_vel = 5.
turbinflow.high.turb_nplane = 32
turbinflow.high.time_offset = 0.0006
turbinflow.high.verbose = 2

Plt File Management
===================
Expand Down
8 changes: 4 additions & 4 deletions Source/Eos/GammaLaw.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -478,10 +478,10 @@ struct GammaLaw
AMREX_GPU_HOST_DEVICE
GammaLaw(const EosParm<GammaLaw>* eparm)
{
AMREX_ASSERT_WITH_MESSAGE(
eparm != nullptr,
"GammaLaw Eos initialized with eosparm, but eosparm was not initialized");
gamma = eparm->gamma;
// This should eventually error out if eparm is nullptr
// But for now we require that eparm->gamma=gamma anyway,
// and this supports some potential legacy use cases for now
gamma = (eparm != nullptr) ? eparm->gamma : gamma;
}

template <class... Args>
Expand Down
55 changes: 20 additions & 35 deletions Source/Reactions/ReactorArkode.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -199,31 +199,31 @@ ReactorArkode::react(
if (use_erkstep == 0) {
arkode_mem = ARKStepCreate(
cF_RHS, nullptr, time, y, *amrex::sundials::The_Sundials_Context());
ARKStepSetUserData(arkode_mem, static_cast<void*>(user_data));
ARKodeSetUserData(arkode_mem, static_cast<void*>(user_data));
utils::set_sundials_solver_tols<Ordering>(
*amrex::sundials::The_Sundials_Context(), arkode_mem, user_data->ncells,
relTol, absTol, m_typ_vals, "arkstep", verbose);
ARKStepSetTableNum(
arkode_mem, ARKODE_DIRK_NONE, static_cast<ARKODE_ERKTableID>(rk_method));
int flag = ARKStepSetAdaptController(arkode_mem, sun_controller);
utils::check_flag(&flag, "ARKStepSetAdaptController", 1);
int flag = ARKodeSetAdaptController(arkode_mem, sun_controller);
utils::check_flag(&flag, "ARKodeSetAdaptController", 1);
BL_PROFILE_VAR(
"Pele::ReactorArkode::react():ARKStepEvolve", AroundARKEvolve);
ARKStepEvolve(arkode_mem, time_out, y, &time_init, ARK_NORMAL);
ARKodeEvolve(arkode_mem, time_out, y, &time_init, ARK_NORMAL);
BL_PROFILE_VAR_STOP(AroundARKEvolve);
} else {
arkode_mem =
ERKStepCreate(cF_RHS, time, y, *amrex::sundials::The_Sundials_Context());
ERKStepSetUserData(arkode_mem, static_cast<void*>(user_data));
ARKodeSetUserData(arkode_mem, static_cast<void*>(user_data));
utils::set_sundials_solver_tols<Ordering>(
*amrex::sundials::The_Sundials_Context(), arkode_mem, user_data->ncells,
relTol, absTol, m_typ_vals, "erkstep", verbose);
ERKStepSetTableNum(arkode_mem, static_cast<ARKODE_ERKTableID>(rk_method));
int flag = ERKStepSetAdaptController(arkode_mem, sun_controller);
utils::check_flag(&flag, "ERKStepSetAdaptController", 1);
int flag = ARKodeSetAdaptController(arkode_mem, sun_controller);
utils::check_flag(&flag, "ARKodeSetAdaptController", 1);
BL_PROFILE_VAR(
"Pele::ReactorArkode::react():ERKStepEvolve", AroundERKEvolve);
ERKStepEvolve(arkode_mem, time_out, y, &time_init, ARK_NORMAL);
ARKodeEvolve(arkode_mem, time_out, y, &time_init, ARK_NORMAL);
BL_PROFILE_VAR_STOP(AroundERKEvolve);
}

Expand All @@ -250,11 +250,7 @@ ReactorArkode::react(
user_data->rhoe_init, d_nfe, dt_react);

N_VDestroy(y);
if (use_erkstep == 0) {
ARKStepFree(&arkode_mem);
} else {
ERKStepFree(&arkode_mem);
}
ARKodeFree(&arkode_mem);

delete user_data;

Expand Down Expand Up @@ -339,24 +335,24 @@ ReactorArkode::react(
if (use_erkstep == 0) {
arkode_mem = ARKStepCreate(
cF_RHS, nullptr, time, y, *amrex::sundials::The_Sundials_Context());
ARKStepSetUserData(arkode_mem, static_cast<void*>(user_data));
ARKodeSetUserData(arkode_mem, static_cast<void*>(user_data));
utils::set_sundials_solver_tols<Ordering>(
*amrex::sundials::The_Sundials_Context(), arkode_mem, user_data->ncells,
relTol, absTol, m_typ_vals, "arkstep", verbose);
BL_PROFILE_VAR(
"Pele::ReactorArkode::react():ARKStepEvolve", AroundARKEvolve);
ARKStepEvolve(arkode_mem, time_out, y, &time_init, ARK_NORMAL);
ARKodeEvolve(arkode_mem, time_out, y, &time_init, ARK_NORMAL);
BL_PROFILE_VAR_STOP(AroundARKEvolve);
} else {
arkode_mem =
ERKStepCreate(cF_RHS, time, y, *amrex::sundials::The_Sundials_Context());
ERKStepSetUserData(arkode_mem, static_cast<void*>(user_data));
ARKodeSetUserData(arkode_mem, static_cast<void*>(user_data));
utils::set_sundials_solver_tols<Ordering>(
*amrex::sundials::The_Sundials_Context(), arkode_mem, user_data->ncells,
relTol, absTol, m_typ_vals, "erkstep", verbose);
BL_PROFILE_VAR(
"Pele::ReactorArkode::react():ERKStepEvolve", AroundERKEvolve);
ERKStepEvolve(arkode_mem, time_out, y, &time_init, ARK_NORMAL);
ARKodeEvolve(arkode_mem, time_out, y, &time_init, ARK_NORMAL);
BL_PROFILE_VAR_STOP(AroundERKEvolve);
}
#ifdef MOD_REACTOR
Expand Down Expand Up @@ -386,11 +382,7 @@ ReactorArkode::react(
}

N_VDestroy(y);
if (use_erkstep == 0) {
ARKStepFree(&arkode_mem);
} else {
ERKStepFree(&arkode_mem);
}
ARKodeFree(&arkode_mem);

delete user_data;

Expand Down Expand Up @@ -436,23 +428,16 @@ ReactorArkode::print_final_stats(void* arkode_mem)
long int nst, nst_a, netf, nfe, nfi;
int flag;

flag = ARKodeGetNumSteps(arkode_mem, &nst);
utils::check_flag(&flag, "ARKodeGetNumSteps", 1);
flag = ARKodeGetNumStepAttempts(arkode_mem, &nst_a);
utils::check_flag(&flag, "ARKodeGetNumStepAttempts", 1);
flag = ARKodeGetNumErrTestFails(arkode_mem, &netf);
utils::check_flag(&flag, "ARKodeGetNumErrTestFails", 1);
if (use_erkstep != 0) {
flag = ERKStepGetNumSteps(arkode_mem, &nst);
utils::check_flag(&flag, "ERKStepGetNumSteps", 1);
flag = ERKStepGetNumStepAttempts(arkode_mem, &nst_a);
utils::check_flag(&flag, "ERKStepGetNumStepAttempts", 1);
flag = ERKStepGetNumErrTestFails(arkode_mem, &netf);
utils::check_flag(&flag, "ERKStepGetNumErrTestFails", 1);
flag = ERKStepGetNumRhsEvals(arkode_mem, &nfe);
utils::check_flag(&flag, "ERKStepGetNumRhsEvals", 1);

} else {
flag = ARKStepGetNumSteps(arkode_mem, &nst);
utils::check_flag(&flag, "ARKStepGetNumSteps", 1);
flag = ARKStepGetNumStepAttempts(arkode_mem, &nst_a);
utils::check_flag(&flag, "ARKStepGetNumStepAttempts", 1);
flag = ARKStepGetNumErrTestFails(arkode_mem, &netf);
utils::check_flag(&flag, "ARKStepGetNumErrTestFails", 1);
flag = ARKStepGetNumRhsEvals(arkode_mem, &nfe, &nfi);
utils::check_flag(&flag, "ARKStepGetNumRhsEvals", 1);
}
Expand Down
6 changes: 2 additions & 4 deletions Source/Reactions/ReactorUtils.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -146,10 +146,8 @@ set_sundials_solver_tols(
int flag;
if (solvername == "cvode") {
flag = CVodeSVtolerances(sundials_mem, relTol, atol);
} else if (solvername == "arkstep") {
flag = ARKStepSVtolerances(sundials_mem, relTol, atol);
} else if (solvername == "erkstep") {
flag = ERKStepSVtolerances(sundials_mem, relTol, atol);
} else if ((solvername == "arkstep") || (solvername == "erkstep")) {
flag = ARKodeSVtolerances(sundials_mem, relTol, atol);
} else {
amrex::Abort("setSundialsSolverTols not implemented for this solver type");
}
Expand Down
10 changes: 7 additions & 3 deletions Source/Utility/Diagnostics/DiagFramePlane.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,10 @@ DiagFramePlane::init(const std::string& a_prefix, std::string_view a_diagName)
{
DiagBase::init(a_prefix, a_diagName);

AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
AMREX_SPACEDIM == 3, "DiagFramePlane::init(): DiagFramePlane is only "
"supported for 3 dimensional simulations");

if (m_filters.empty()) {
amrex::Print() << " Filters are not available on DiagFramePlane and will "
"be discarded \n";
Expand Down Expand Up @@ -117,9 +121,9 @@ DiagFramePlane::prepare(
cdim += 1;
}
}
initDomain.setRange(2, 0, 1);
initRealBox.setLo(2, 0.0);
initRealBox.setHi(2, dxlcl[2]);
initDomain.setRange(AMREX_SPACEDIM - 1, 0, 1);
initRealBox.setLo(AMREX_SPACEDIM - 1, 0.0);
initRealBox.setHi(AMREX_SPACEDIM - 1, dxlcl[2]);
m_geomLev0.define(
initDomain, initRealBox, a_geoms[0].Coord(),
amrex::Array<int, AMREX_SPACEDIM>({AMREX_D_DECL(0, 0, 0)}));
Expand Down
13 changes: 8 additions & 5 deletions Source/Utility/TurbInflow/turbinflow.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,12 +13,15 @@ struct TurbParm
// Turbulence data file
std::string m_turb_file;

// Verbosity
int verbose = 1;

// Geometry of the turbulent data
amrex::GpuArray<int, AMREX_SPACEDIM> npboxcells = {{0}};
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> pboxlo = {{0.0}};
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> dx = {{0.0}};
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> dxinv = {{0.0}};
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> pboxsize = {{0.0}};
amrex::GpuArray<int, 3> npboxcells = {{0}};
amrex::GpuArray<amrex::Real, 3> pboxlo = {{0.0}};
amrex::GpuArray<amrex::Real, 3> dx = {{0.0}};
amrex::GpuArray<amrex::Real, 3> dxinv = {{0.0}};
amrex::GpuArray<amrex::Real, 3> pboxsize = {{0.0}};
int dir; // Direction where to use this TurbParm
amrex::Orientation::Side side; // Side on which to use this TurbParm
amrex::Real time_shift =
Expand Down
59 changes: 23 additions & 36 deletions Source/Utility/TurbInflow/turbinflow.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,9 @@ TurbInflow::init(amrex::Geometry const& /*geom*/)
n_tp = ppr.countval("turbinflows");
amrex::Vector<std::string> tp_list;
if (n_tp > 0) {
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
AMREX_SPACEDIM == 3, "TurbInflow::init(): TurbInflows are only supported "
"for 3 dimensional simulations for now");
tp.resize(n_tp);
tp_list.resize(n_tp);
for (int n = 0; n < n_tp; n++) {
Expand Down Expand Up @@ -41,12 +44,15 @@ TurbInflow::init(amrex::Geometry const& /*geom*/)
pp.query("time_offset", tp[n].time_shift);
pp.query("turb_scale_loc", tp[n].turb_scale_loc);
pp.query("turb_scale_vel", tp[n].turb_scale_vel);
amrex::Print() << "Initializing turbInflow " << tp_list[n]
<< " with file " << tp[n].m_turb_file
<< " (location coordinates in will be scaled by "
<< tp[n].turb_scale_loc
<< " and velocity out to be scaled by "
<< tp[n].turb_scale_vel << ") \n";
pp.query("verbose", tp[n].verbose);
if (tp[n].verbose > 0) {
amrex::Print() << "Initializing turbInflow " << tp_list[n]
<< " with file " << tp[n].m_turb_file
<< " (location coordinates in will be scaled by "
<< tp[n].turb_scale_loc
<< " and velocity out to be scaled by "
<< tp[n].turb_scale_vel << ") \n";
}

// Get the turbcenter on the injection face
amrex::Vector<amrex::Real> turb_center(AMREX_SPACEDIM - 1, 0);
Expand Down Expand Up @@ -189,16 +195,6 @@ TurbInflow::add_turb(

// Moving it into data
set_turb(dir, tdir1, tdir2, v, data, dcomp);

#if 0
std::string junk = "TurbV_AftTP"+std::to_string(n)+"_D";
std::ofstream os;
os.precision(15);
os.open(junk.c_str());
data.writeOn(os);
os.close();
amrex::Abort();
#endif
}

void
Expand Down Expand Up @@ -246,8 +242,8 @@ TurbInflow::read_one_turb_plane(TurbParm& a_tp, int iplane, int k)
}

amrex::Box dstBox = a_tp.sdata->box();
dstBox.setSmall(2, iplane);
dstBox.setBig(2, iplane);
dstBox.setSmall(AMREX_SPACEDIM - 1, iplane);
dstBox.setBig(AMREX_SPACEDIM - 1, iplane);

for (int n = 0; n < AMREX_SPACEDIM; ++n) {

Expand Down Expand Up @@ -280,25 +276,17 @@ TurbInflow::read_turb_planes(TurbParm& a_tp, amrex::Real z)
a_tp.szlo = static_cast<amrex::Real>(izlo) * a_tp.dx[2];
a_tp.szhi = static_cast<amrex::Real>(izhi) * a_tp.dx[2];

#if 0
amrex::AllPrint() << "read_turb_planes filling " << izlo << " to " << izhi
<< " covering " << a_tp.szlo + 0.5 * a_tp.dx[2]
<< " to " << a_tp.szhi - 0.5 * a_tp.dx[2] << " for z = " << z << std::endl;
#endif
if (a_tp.verbose > 1) {
amrex::Print() << "read_turb_planes filling " << izlo << " to " << izhi
<< " covering " << a_tp.szlo + 0.5 * a_tp.dx[2] << " to "
<< a_tp.szhi - 0.5 * a_tp.dx[2] << " for z = " << z
<< std::endl;
}

for (int iplane = 1; iplane <= a_tp.nplane; ++iplane) {
int k = (izlo + iplane - 1) % (a_tp.npboxcells[2] - 2);
read_one_turb_plane(a_tp, iplane, k);
}
#if 0
int myproc = amrex::ParallelDescriptor::MyProc();
std::string junk = "TurbData_proc"+std::to_string(myproc)+"_D";
std::ofstream os;
os.precision(15);
os.open(junk.c_str());
a_tp.sdata->writeOn(os);
os.close();
#endif
}

void
Expand All @@ -311,12 +299,11 @@ TurbInflow::fill_turb_plane(
{
if (
(z < a_tp.szlo + 0.5 * a_tp.dx[2]) || (z > a_tp.szhi - 0.5 * a_tp.dx[2])) {
#if 0
{
amrex::AllPrint() << "Reading new data because z " << z << " is outside " << a_tp.szlo + 0.5 * a_tp.dx[2] << " and "
if (a_tp.verbose > 1) {
amrex::Print() << "Reading new data because z " << z << " is outside "
<< a_tp.szlo + 0.5 * a_tp.dx[2] << " and "
<< a_tp.szhi - 0.5 * a_tp.dx[2] << std::endl;
}
#endif
read_turb_planes(a_tp, z);
}

Expand Down
2 changes: 1 addition & 1 deletion Submodules/sundials
Submodule sundials updated 928 files

0 comments on commit 3b718ff

Please sign in to comment.