From ac193c4a4694438274eb138fac5a2c0e6e3a1b84 Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Mon, 30 Sep 2024 16:05:10 -0600 Subject: [PATCH 01/46] Initial commit. Needs debugging for source terms. --- CMake/BuildPeleExe.cmake | 1 + Exec/Make.PeleLMeX | 4 + .../RegTests/EB_DecayingODEQTY/CMakeLists.txt | 8 + Exec/RegTests/EB_DecayingODEQTY/GNUmakefile | 39 ++++ .../PeleLMeX_EBUserDefined.H | 58 +++++ .../PeleLMeX_ProblemSpecificFunctions.H | 19 ++ Exec/RegTests/EB_DecayingODEQTY/README.md | 2 + Exec/RegTests/EB_DecayingODEQTY/input.2d | 69 ++++++ .../EB_DecayingODEQTY/pelelmex_prob.H | 204 ++++++++++++++++++ .../EB_DecayingODEQTY/pelelmex_prob.cpp | 27 +++ .../EB_DecayingODEQTY/pelelmex_prob_parm.H | 21 ++ Source/Make.package | 1 + Source/PeleLMeX.H | 9 + Source/PeleLMeX_Advance.cpp | 3 + Source/PeleLMeX_Forces.cpp | 20 ++ Source/PeleLMeX_Index.H | 8 +- Source/PeleLMeX_Plot.cpp | 5 + Source/PeleLMeX_ProblemSpecificFunctions.H | 34 +++ Source/PeleLMeX_Setup.cpp | 16 ++ 19 files changed, 547 insertions(+), 1 deletion(-) create mode 100644 Exec/RegTests/EB_DecayingODEQTY/CMakeLists.txt create mode 100644 Exec/RegTests/EB_DecayingODEQTY/GNUmakefile create mode 100644 Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_EBUserDefined.H create mode 100644 Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.H create mode 100644 Exec/RegTests/EB_DecayingODEQTY/README.md create mode 100644 Exec/RegTests/EB_DecayingODEQTY/input.2d create mode 100644 Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob.H create mode 100644 Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob.cpp create mode 100644 Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob_parm.H create mode 100644 Source/PeleLMeX_ProblemSpecificFunctions.H diff --git a/CMake/BuildPeleExe.cmake b/CMake/BuildPeleExe.cmake index 4d6069e8b..1be05d024 100644 --- a/CMake/BuildPeleExe.cmake +++ b/CMake/BuildPeleExe.cmake @@ -50,6 +50,7 @@ function(build_pele_exe pele_exe_name pele_physics_lib_name) ${SRC_DIR}/PeleLMeX_PatchFlowVariables.cpp ${SRC_DIR}/PeleLMeX_Init.cpp ${SRC_DIR}/PeleLMeX_Plot.cpp + ${SRC_DIR}/PeleLMeX_ProblemSpecificFunctions.H ${SRC_DIR}/PeleLMeX_Projection.cpp ${SRC_DIR}/PeleLMeX_Reactions.cpp ${SRC_DIR}/PeleLMeX_Regrid.cpp diff --git a/Exec/Make.PeleLMeX b/Exec/Make.PeleLMeX index 3580e0fa3..4ef1c7734 100644 --- a/Exec/Make.PeleLMeX +++ b/Exec/Make.PeleLMeX @@ -57,6 +57,10 @@ ifeq ($(USE_RADIATION), TRUE) DEFINES += -DPELE_USE_RADIATION endif +ifeq ($(shell test 0$(NUM_ODE) -gt 0; echo $$?), 0) + DEFINES += -DNUM_ODE=$(NUM_ODE) +endif + Bpack += $(foreach dir, $(LMdirs), $(PELE_HOME)/$(dir)/Make.package) Blocs += $(foreach dir, $(LMdirs), $(PELE_HOME)/$(dir)) diff --git a/Exec/RegTests/EB_DecayingODEQTY/CMakeLists.txt b/Exec/RegTests/EB_DecayingODEQTY/CMakeLists.txt new file mode 100644 index 000000000..7169cbdba --- /dev/null +++ b/Exec/RegTests/EB_DecayingODEQTY/CMakeLists.txt @@ -0,0 +1,8 @@ +set(PELE_PHYSICS_EOS_MODEL Fuego) +set(PELE_PHYSICS_CHEMISTRY_MODEL air) +set(PELE_PHYSICS_TRANSPORT_MODEL Constant) +set(PELE_PHYSICS_ENABLE_SPRAY OFF) +set(PELE_PHYSICS_SPRAY_FUEL_NUM 0) +set(PELE_PHYSICS_ENABLE_SOOT OFF) +set(PELE_PHYSICS_ENABLE_RADIATION OFF) +include(BuildExeAndLib) diff --git a/Exec/RegTests/EB_DecayingODEQTY/GNUmakefile b/Exec/RegTests/EB_DecayingODEQTY/GNUmakefile new file mode 100644 index 000000000..9bce8ffc7 --- /dev/null +++ b/Exec/RegTests/EB_DecayingODEQTY/GNUmakefile @@ -0,0 +1,39 @@ +# AMReX +DIM = 2 +COMP = gnu +PRECISION = DOUBLE +USE_EB = TRUE +USE_HYPRE = FALSE + +# Profiling +PROFILE = FALSE +TINY_PROFILE = FALSE +COMM_PROFILE = FALSE +TRACE_PROFILE = FALSE +MEM_PROFILE = FALSE +USE_GPROF = FALSE + +# Performance +USE_MPI = TRUE +USE_OMP = FALSE +USE_CUDA = FALSE +USE_HIP = FALSE +USE_SYCL = FALSE + +# Debugging +DEBUG = FALSE +FSANITIZER = FALSE +THREAD_SANITIZER = FALSE + +# PeleLMeX +CEXE_headers += EBUserDefined.H +CEXE_headers += ProblemSpecificcFunctions.H +NUM_ODE = 1 + +# PelePhysics +Chemistry_Model = air +Eos_Model = Fuego +Transport_Model = Constant + +PELE_HOME ?= ../../.. +include $(PELE_HOME)/Exec/Make.PeleLMeX diff --git a/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_EBUserDefined.H b/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_EBUserDefined.H new file mode 100644 index 000000000..9de8406f9 --- /dev/null +++ b/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_EBUserDefined.H @@ -0,0 +1,58 @@ +#ifndef EBUSERDEFINED_H +#define EBUSERDEFINED_H + +using namespace amrex; + +#ifdef AMREX_USE_EB +#include +#include +void +EBUserDefined( + const Geometry& /*geom*/, + const int /*required_coarsening_level*/, + const int /*max_coarsening_level*/) +{ + // ParmParse your geometry parameters + + // Build geometry pieces using EB2::* methods + + // Build your geometry shop using EB2::makeShop + + // Build geom using EB2::Build + + // We shoulnd't be here, copy this file in your run folder + // and implement your geometry + Abort("In default EBUserDefined function! Shouldn't be here. Copy and edit " + "this file for your needs"); +} + +AMREX_GPU_DEVICE +AMREX_FORCE_INLINE +void +setEBState( + const amrex::Real xEBface[AMREX_SPACEDIM], + amrex::Real s_ext[NVAR], + const amrex::Real /*time*/, + amrex::GeometryData const& /*geomdata*/, + ProbParm const& /*prob_parm*/) +{ + if (xEBface[1] > 0.02) { + s_ext[TEMP] = 500.0; + } else { + s_ext[TEMP] = 300.0; + } +} + +AMREX_GPU_DEVICE +AMREX_FORCE_INLINE +void +setEBType( + const amrex::Real* /*xEBface[AMREX_SPACEDIM]*/, + amrex::Real& EBflagType, + amrex::GeometryData const& /*geomdata*/, + ProbParm const& /*prob_parm*/) +{ + EBflagType = 1.0; +} +#endif +#endif diff --git a/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.H b/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.H new file mode 100644 index 000000000..57a6b490c --- /dev/null +++ b/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.H @@ -0,0 +1,19 @@ +#ifndef PROBLEMSPECIFICFUNCTIONS_H +#define PROBLEMSPECIFICFUNCTIONS_H + +using namespace amrex; + +#include +#include + +void set_ode_names(amrex::Vector& a_ode_names) +{ + a_ode_names.resize(NUM_ODE); +#if NUM_ODE > 0 + for (int n = 0; n < NUM_ODE; n++) { + a_ode_names[n] = "MY_ODE_" + std::to_string(n); + } +#endif +} + +#endif \ No newline at end of file diff --git a/Exec/RegTests/EB_DecayingODEQTY/README.md b/Exec/RegTests/EB_DecayingODEQTY/README.md new file mode 100644 index 000000000..19fb1a405 --- /dev/null +++ b/Exec/RegTests/EB_DecayingODEQTY/README.md @@ -0,0 +1,2 @@ +## EB\_DecayingODEQty +A 2D (or 3D) inflow/outflow setup with an optional EB cylinder in the middle of the flow. Testing the "ODE" variables. diff --git a/Exec/RegTests/EB_DecayingODEQTY/input.2d b/Exec/RegTests/EB_DecayingODEQTY/input.2d new file mode 100644 index 000000000..de6ec7699 --- /dev/null +++ b/Exec/RegTests/EB_DecayingODEQTY/input.2d @@ -0,0 +1,69 @@ +#----------------------DOMAIN DEFINITION------------------------- +geometry.is_periodic = 0 1 # For each dir, 0: non-perio, 1: periodic +geometry.coord_sys = 0 # 0 => cart, 1 => RZ +geometry.prob_lo = -0.02 -0.02 # x_lo y_lo (z_lo) +geometry.prob_hi = 0.10 0.02 # x_hi y_hi (z_hi) + +#--------------------------BC FLAGS------------------------------ +# Interior, Inflow, Outflow, Symmetry, +# SlipWallAdiab, NoSlipWallAdiab, SlipWallIsotherm, NoSlipWallIsotherm +peleLM.lo_bc = Inflow Interior +peleLM.hi_bc = Outflow Interior + +#-------------------------AMR CONTROL---------------------------- +amr.n_cell = 192 64 # Level 0 number of cells in each direction +amr.v = 1 # AMR verbose +amr.max_level = 0 # maximum level number allowed +amr.ref_ratio = 2 2 2 2 # refinement ratio +amr.regrid_int = 5 # how often to regrid +amr.n_error_buf = 2 2 2 2 # number of buffer cells in error est +amr.grid_eff = 0.7 # what constitutes an efficient grid +amr.blocking_factor = 16 # block factor in grid generation (min box size) +amr.max_grid_size = 64 # max box size + +#--------------------------- Problem ---------------------------- +prob.T_mean = 300.0 +prob.P_mean = 101325.0 +prob.meanFlowMag = 4.255 +prob.meanFlowDir = 1 + +#-----------------INPUTS TO CONSTANT TRANSPORT ------------------ +transport.units = MKS +transport.const_viscosity = 0.0001 +transport.const_bulk_viscosity = 0.0 +transport.const_conductivity = 0.0 +transport.const_diffusivity = 0.0 + +#---------------------------TIME STEPPING------------------------ +amr.max_step = 1 +amr.stop_time = 0.025 +amr.dt_shrink = 0.1 +amr.dt_change_max = 1.1 +amr.cfl = 0.7 + +#-------------------------PELE CONTROLS---------------------- +peleLM.v = 1 + +#---------------------------IO CONTROL--------------------------- +#amr.restart = chk_07136 +amr.check_file = "chk_" +amr.check_per = -1 +amr.plot_file = "plt_" +amr.plot_per = 0.005 +amr.derive_plot_vars = avg_pressure mag_vort + +#------------------------- EB SETUP ----------------------------- +eb2.geom_type = sphere +eb2.sphere_radius = 0.005 +eb2.sphere_center = 0.03 0.0 +eb2.sphere_has_fluid_inside = 0 +eb2.small_volfrac = 1.0e-4 + +#--------------------REFINEMENT CONTROL-------------------------- +amr.refinement_indicators = VortL VortH +amr.VortL.max_level = 1 +amr.VortL.value_less = -1000 +amr.VortL.field_name = mag_vort +amr.VortH.max_level = 1 +amr.VortH.value_greater = 1000 +amr.VortH.field_name = mag_vort diff --git a/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob.H b/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob.H new file mode 100644 index 000000000..62fa67d1c --- /dev/null +++ b/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob.H @@ -0,0 +1,204 @@ +#ifndef PELELM_PROB_H +#define PELELM_PROB_H + +#include +#include +#include + +#include +#include +#include +#include + +AMREX_GPU_DEVICE +AMREX_FORCE_INLINE +void +pelelmex_initdata( + int i, + int j, + int k, + int is_incompressible, + amrex::Array4 const& state, + amrex::Array4 const& /*aux*/, + amrex::GeometryData const& geomdata, + ProbParm const& prob_parm, + pele::physics::PMF::PmfData::DataContainer const* /*pmf_data*/) +{ + auto eos = pele::physics::PhysicsType::eos(); + + amrex::Real massfrac[NUM_SPECIES] = {0.0}; + massfrac[O2_ID] = 0.233; + massfrac[N2_ID] = 0.767; + + switch (prob_parm.meanFlowDir) { + case 1: + AMREX_D_TERM(state(i, j, k, VELX) = prob_parm.meanFlowMag; + , state(i, j, k, VELY) = 0.0;, state(i, j, k, VELZ) = 0.0); + break; + case -1: + AMREX_D_TERM(state(i, j, k, VELX) = -prob_parm.meanFlowMag; + , state(i, j, k, VELY) = 0.0;, state(i, j, k, VELZ) = 0.0); + break; + case 2: + AMREX_D_TERM(state(i, j, k, VELX) = 0.0; + , state(i, j, k, VELY) = prob_parm.meanFlowMag; + , state(i, j, k, VELZ) = 0.0); + break; + case -2: + AMREX_D_TERM(state(i, j, k, VELX) = 0.0; + , state(i, j, k, VELY) = -prob_parm.meanFlowMag; + , state(i, j, k, VELZ) = 0.0); + break; + case 3: + AMREX_D_TERM(state(i, j, k, VELX) = prob_parm.meanFlowMag; + , state(i, j, k, VELY) = prob_parm.meanFlowMag; + , state(i, j, k, VELZ) = 0.0); + break; + case -3: + AMREX_D_TERM(state(i, j, k, VELX) = -prob_parm.meanFlowMag; + , state(i, j, k, VELY) = -prob_parm.meanFlowMag; + , state(i, j, k, VELZ) = 0.0); + break; + } + + if (!(is_incompressible == 0)) { + return; + } + + state(i, j, k, TEMP) = prob_parm.T_mean; + + amrex::Real P_cgs = prob_parm.P_mean * 10.0; + + // Density + amrex::Real rho_cgs = 0.0; + eos.PYT2R(P_cgs, massfrac, state(i, j, k, TEMP), rho_cgs); + state(i, j, k, DENSITY) = rho_cgs * 1.0e3; + + // Enthalpy + amrex::Real h_cgs = 0.0; + eos.TY2H(state(i, j, k, TEMP), massfrac, h_cgs); + state(i, j, k, RHOH) = h_cgs * 1.0e-4 * state(i, j, k, DENSITY); + + // Species mass + for (int n = 0; n < NUM_SPECIES; n++) { + state(i, j, k, FIRSTSPEC + n) = massfrac[n] * state(i, j, k, DENSITY); + } + + // Initial ODE quantity + { +#if NUM_ODE > 0 + for (int n = 0; n < NUM_ODE; n++) { + state(i, j, k, FIRSTODE + n) = 10.0; + } + /* + amrex::Real ode_qty[NUM_ODE] = {0.0}; + + // Initial state for ode quantities + amrex::Real ode_qty_0[NUM_ODE] = {0.0}; + for (int n = 0; n < NUM_ODE; n++) { + ode_qty_0[n] = (n + 1) * prob_parm.ode_IC; + } + + // Current x,y,z locations + const amrex::Real* prob_lo = geomdata.ProbLo(); + const amrex::Real* dx = geomdata.CellSize(); + const amrex::Real x[AMREX_SPACEDIM] = {AMREX_D_DECL( + prob_lo[0] + static_cast(i + 0.5) * dx[0], + prob_lo[1] + static_cast(j + 0.5) * dx[1], + prob_lo[2] + static_cast(k + 0.5) * dx[2])}; + + // Create NUM_ODE boxes from (x_strt,y_strt) to (x_end,y_end) + for (int n = 0; n < NUM_ODE; n++) { + amrex::Real x_strt = prob_parm.ode_xy_lo + n * 4. * prob_parm.ode_length; + amrex::Real x_end = x_strt + prob_parm.ode_length; + amrex::Real y_strt = prob_parm.ode_xy_lo; + amrex::Real y_end = y_strt + prob_parm.ode_height; + if ((x[0] > x_strt) && (x[1] > y_strt)) { + if ((x[0] < x_end) && (x[1] < y_end)) { + ode_qty[n] = ode_qty_0[n]; + } + } + } + // Add to state variable + for (int n = 0; n < NUM_ODE; n++) { + state(i, j, k, FIRSTODE + n) = ode_qty[n]; + } + */ +#endif + } +} + +AMREX_GPU_DEVICE +AMREX_FORCE_INLINE +void +bcnormal( + const amrex::Real* /*x[AMREX_SPACEDIM]*/, + const int /*m_node*/, + amrex::Real s_ext[NVAR], + const int /*idir*/, + const int /*sgn*/, + const amrex::Real /*time*/, + amrex::GeometryData const& /*geomdata*/, + ProbParm const& prob_parm, + pele::physics::PMF::PmfData::DataContainer const* /*pmf_data*/) +{ + // amrex::GpuArray pmf_vals = {0.0}; + amrex::Real massfrac[NUM_SPECIES] = {0.0}; + + auto eos = pele::physics::PhysicsType::eos(); + + switch (prob_parm.meanFlowDir) { + case 1: + AMREX_D_TERM(s_ext[VELX] = prob_parm.meanFlowMag;, s_ext[VELY] = 0.0; + , s_ext[VELZ] = 0.0); + break; + case 2: + AMREX_D_TERM(s_ext[VELY] = prob_parm.meanFlowMag;, s_ext[VELX] = 0.0; + , s_ext[VELZ] = 0.0); + break; + case 3: + AMREX_D_TERM(s_ext[VELZ] = prob_parm.meanFlowMag;, s_ext[VELX] = 0.0; + , s_ext[VELY] = 0.0); + break; + } + + massfrac[O2_ID] = 0.333; + massfrac[N2_ID] = 0.667; + + s_ext[TEMP] = prob_parm.T_mean; + + amrex::Real rho_cgs, P_cgs, RhoH_temp; + P_cgs = prob_parm.P_mean * 10.0; + + eos.PYT2R(P_cgs, massfrac, s_ext[TEMP], rho_cgs); + s_ext[DENSITY] = rho_cgs * 1.0e3; + + eos.TY2H(s_ext[TEMP], massfrac, RhoH_temp); + s_ext[RHOH] = RhoH_temp * 1.0e-4 * s_ext[DENSITY]; // CGS -> MKS conversion + + for (int n = 0; n < NUM_SPECIES; n++) { + s_ext[FIRSTSPEC + n] = massfrac[n] * s_ext[DENSITY]; + } +} + +AMREX_GPU_DEVICE +AMREX_FORCE_INLINE +void +zero_visc( + int i, + int j, + int k, + amrex::Array4 const& beta, + amrex::GeometryData const& geomdata, + amrex::Box const& domainBox, + const int dir, + const int beta_comp, + const int nComp) +{ + amrex::ignore_unused( + i, j, k, beta, geomdata, domainBox, dir, beta_comp, nComp); + // We treat species when beta_comp == 0 and nComp == NUM_SPECIES + // otherwise this routine could be called for other face diffusivity (Temp, + // velocity, ...) +} +#endif diff --git a/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob.cpp b/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob.cpp new file mode 100644 index 000000000..1d39f11bd --- /dev/null +++ b/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob.cpp @@ -0,0 +1,27 @@ +#include +#include + +void +PeleLM::readProbParm() // NOLINT(readability-make-member-function-const) +{ + amrex::ParmParse pp("prob"); + + pp.query("T_mean", prob_parm->T_mean); + pp.query("P_mean", prob_parm->P_mean); + pp.query("meanFlowDir", prob_parm->meanFlowDir); + pp.query("meanFlowMag", prob_parm->meanFlowMag); + +#if NUM_ODE > 0 + amrex::Print() << "\n!!!!!! NUM_ODE = " << NUM_ODE << " !!!!!!!!!" << "\n"; +#endif + + // if (!m_incompressible) { + // auto& trans_parm = PeleLM::trans_parms.host_parm(); + // amrex::ParmParse pptr("transport"); + // pp.query("const_viscosity", trans_parm.const_viscosity); + // pp.query("const_bulk_viscosity", trans_parm.const_bulk_viscosity); + // pp.query("const_conductivity", trans_parm.const_conductivity); + // pp.query("const_diffusivity", trans_parm.const_diffusivity); + // PeleLM::trans_parms.sync_to_device(); + // } +} diff --git a/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob_parm.H b/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob_parm.H new file mode 100644 index 000000000..34785b133 --- /dev/null +++ b/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob_parm.H @@ -0,0 +1,21 @@ +#ifndef PELELM_PROB_PARM_H +#define PELELM_PROB_PARM_H + +#include +#include + +using namespace amrex::literals; + +struct ProbParm : amrex::Gpu::Managed +{ + // Shared params + amrex::Real T_mean = 298.0_rt; + amrex::Real P_mean = 101325.0_rt; + amrex::Real meanFlowMag = 0.0; + amrex::Real ode_IC = 10.0; + amrex::Real ode_xy_lo = -0.015; + amrex::Real ode_length = 0.005; + amrex::Real ode_height = 0.010; + int meanFlowDir = 1; +}; +#endif diff --git a/Source/Make.package b/Source/Make.package index 1cd34c29f..eca6d42cd 100644 --- a/Source/Make.package +++ b/Source/Make.package @@ -10,6 +10,7 @@ CEXE_headers += PeleLMeX_EBUserDefined.H CEXE_headers += PeleLMeX_FlowControllerData.H CEXE_headers += PeleLMeX_BPatch.H CEXE_headers += PeleLMeX_PatchFlowVariables.H +CEXE_headers += PeleLMeX_ProblemSpecificFunctions.H ## Sources CEXE_sources += main.cpp diff --git a/Source/PeleLMeX.H b/Source/PeleLMeX.H index c0c13f673..598110218 100644 --- a/Source/PeleLMeX.H +++ b/Source/PeleLMeX.H @@ -874,6 +874,10 @@ public: void computeRadSource(const PeleLM::TimeStamp& a_timestamp); #endif +#if NUM_ODE > 0 + // Compute the ode source term + void computeODESource(const PeleLM::TimeStamp& a_timestamp); +#endif //----------------------------------------------------------------------------- // EOS @@ -1851,6 +1855,11 @@ public: amrex::Real m_divu_dtFactor = 0.5; amrex::Real m_divu_rhoMin = 0.1; + // ODE variable names + #if NUM_ODE > 0 + amrex::Vector m_ode_names; + #endif + // Aux data size int m_nAux = 0; diff --git a/Source/PeleLMeX_Advance.cpp b/Source/PeleLMeX_Advance.cpp index 905271bcf..110b86681 100644 --- a/Source/PeleLMeX_Advance.cpp +++ b/Source/PeleLMeX_Advance.cpp @@ -131,6 +131,9 @@ PeleLM::Advance(int is_initIter) BL_PROFILE_VAR_STOP(PLM_RAD); } #endif +#if NUM_ODE > 0 + //computeODESource(AmrOldTime); +#endif if (m_incompressible == 0) { floorSpecies(AmrOldTime); diff --git a/Source/PeleLMeX_Forces.cpp b/Source/PeleLMeX_Forces.cpp index 1c22d278c..a91854444 100644 --- a/Source/PeleLMeX_Forces.cpp +++ b/Source/PeleLMeX_Forces.cpp @@ -233,3 +233,23 @@ PeleLM::addSpark(const TimeStamp& a_timestamp) } } } + + +void +PeleLM::computeODESource(const TimeStamp& a_timestamp) +{ + for (int lev = 0; lev <= finest_level; lev++) { + for (int n = 0; n <= NUM_ODE; n++){ + //Real time = getTime(lev, a_timestamp); + //auto statema = getLevelDataPtr(lev, a_timestamp)->state.const_arrays(); + auto extma = m_extSource[lev]->arrays(); + + amrex::ParallelFor( + *m_extSource[lev], + [=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept { + extma[box_no](i, j, k, FIRSTODE + n) = 0.0; + }); + Gpu::streamSynchronize(); + } + } +} \ No newline at end of file diff --git a/Source/PeleLMeX_Index.H b/Source/PeleLMeX_Index.H index 39ca06f26..aa71ea5ae 100644 --- a/Source/PeleLMeX_Index.H +++ b/Source/PeleLMeX_Index.H @@ -22,7 +22,13 @@ #else #define NUMSOOTVAR 0 #endif -#define FIRSTAUX (RHORT + NUMEFIELD + NUMSOOTVAR + 1) + +#ifndef NUM_ODE +#define NUM_ODE 0 +#endif +#define FIRSTODE (RHORT + NUMEFIELD + NUMSOOTVAR + 1) + +#define FIRSTAUX (FIRSTODE + NUM_ODE + 1) #define NVAR (FIRSTAUX + 0) // TODO: handle NAUX diff --git a/Source/PeleLMeX_Plot.cpp b/Source/PeleLMeX_Plot.cpp index ebf9fe58f..8971f189b 100644 --- a/Source/PeleLMeX_Plot.cpp +++ b/Source/PeleLMeX_Plot.cpp @@ -201,6 +201,11 @@ PeleLM::WritePlotFile() plt_VarsName.push_back(sootname); } #endif +#if NUM_ODE > 0 + for (int n = 0; n < NUM_ODE; n++) { + plt_VarsName.push_back(m_ode_names[n]); + } +#endif #ifdef PELE_USE_RADIATION if (do_rad_solve) { plt_VarsName.push_back("rad.G"); diff --git a/Source/PeleLMeX_ProblemSpecificFunctions.H b/Source/PeleLMeX_ProblemSpecificFunctions.H new file mode 100644 index 000000000..21b72b8dc --- /dev/null +++ b/Source/PeleLMeX_ProblemSpecificFunctions.H @@ -0,0 +1,34 @@ +#ifndef PROBLEMSPECIFICFUNCTIONS_H +#define PROBLEMSPECIFICFUNCTIONS_H + +using namespace amrex; + +#include +#include + +void set_ode_names(amrex::Vector& a_ode_names) +{ + a_ode_names.resize(NUM_ODE); +#if NUM_ODE > 0 + for (int n = 0; n < NUM_ODE; n++) { + a_ode_names[n] = "ODE_" + std::to_string(n); + } +#endif +} + +void problem_modify_ext_sources( + amrex::Real /*time*/, + amrex::Real /*dt*/, + const amrex::MultiFab& /*state_old*/, + const amrex::MultiFab& /*state_new*/, + amrex::MultiFab& /*ext_src*/, + int /*ng*/, + amrex::GeometryData const& /*geomdata*/, + ProbParm const& /*prob_parm*/) +{ + /* Notes: ext_src contains sources from velocity forcing coming in + This function should add to rather than overwrite ext_src. + */ +} + +#endif \ No newline at end of file diff --git a/Source/PeleLMeX_Setup.cpp b/Source/PeleLMeX_Setup.cpp index b719ddfbe..0e37bdf7e 100644 --- a/Source/PeleLMeX_Setup.cpp +++ b/Source/PeleLMeX_Setup.cpp @@ -4,6 +4,8 @@ #include #include "PelePhysics.H" #include +#include + #ifdef PELE_USE_EFIELD #include "PeleLMeX_EOS_Extension.H" #endif @@ -15,6 +17,7 @@ #ifdef PELE_USE_SOOT #include "SootModel.H" #endif + using namespace amrex; static Box @@ -813,6 +816,13 @@ PeleLM::variablesSetup() stateComponents.emplace_back(FIRSTSOOT + mom, sootname); } setSootIndx(); +#endif +#if NUM_ODE > 0 + Print() << " First ODE: " << FIRSTODE << "\n"; + set_ode_names(m_ode_names); + for (int n = 0; n < NUM_ODE; n++) { + stateComponents.emplace_back(FIRSTODE + n, m_ode_names[n]); + } #endif } @@ -872,6 +882,12 @@ PeleLM::variablesSetup() m_AdvTypeState[FIRSTSOOT + mom] = 0; m_DiffTypeState[FIRSTSOOT + mom] = 0; } +#endif +#if NUM_ODE > 0 + for (int n = 0; n < NUM_ODE; n++) { + m_AdvTypeState[FIRSTODE + n] = 0; + m_DiffTypeState[FIRSTODE + n] = 0; + } #endif } From f790c8dc5f5546d4f0743559415ba66eda9677f4 Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Tue, 1 Oct 2024 08:48:04 -0600 Subject: [PATCH 02/46] Added ode_var to typical values. Still debugging initialization and source. --- Exec/RegTests/EB_DecayingODEQTY/input.2d | 1 + .../EB_DecayingODEQTY/pelelmex_prob.H | 68 +++++++++---------- .../EB_DecayingODEQTY/pelelmex_prob.cpp | 15 +--- .../EB_DecayingODEQTY/pelelmex_prob_parm.H | 2 +- Source/PeleLMeX_Advance.cpp | 2 +- Source/PeleLMeX_Setup.cpp | 6 -- Source/PeleLMeX_Utils.cpp | 14 ++++ 7 files changed, 51 insertions(+), 57 deletions(-) diff --git a/Exec/RegTests/EB_DecayingODEQTY/input.2d b/Exec/RegTests/EB_DecayingODEQTY/input.2d index de6ec7699..ca4f9a5ac 100644 --- a/Exec/RegTests/EB_DecayingODEQTY/input.2d +++ b/Exec/RegTests/EB_DecayingODEQTY/input.2d @@ -26,6 +26,7 @@ prob.T_mean = 300.0 prob.P_mean = 101325.0 prob.meanFlowMag = 4.255 prob.meanFlowDir = 1 +prob.ode_IC = 100.0 #-----------------INPUTS TO CONSTANT TRANSPORT ------------------ transport.units = MKS diff --git a/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob.H b/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob.H index 62fa67d1c..e1d2c4270 100644 --- a/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob.H +++ b/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob.H @@ -85,47 +85,45 @@ pelelmex_initdata( } // Initial ODE quantity - { #if NUM_ODE > 0 - for (int n = 0; n < NUM_ODE; n++) { - state(i, j, k, FIRSTODE + n) = 10.0; - } - /* - amrex::Real ode_qty[NUM_ODE] = {0.0}; + for (int n = 0; n < NUM_ODE; n++) { + state(i, j, k, FIRSTODE + n) = prob_parm.ode_IC; + } + /* + amrex::Real ode_qty[NUM_ODE] = {0.0}; - // Initial state for ode quantities - amrex::Real ode_qty_0[NUM_ODE] = {0.0}; - for (int n = 0; n < NUM_ODE; n++) { - ode_qty_0[n] = (n + 1) * prob_parm.ode_IC; - } + // Initial state for ode quantities + amrex::Real ode_qty_0[NUM_ODE] = {0.0}; + for (int n = 0; n < NUM_ODE; n++) { + ode_qty_0[n] = (n + 1) * prob_parm.ode_IC; + } - // Current x,y,z locations - const amrex::Real* prob_lo = geomdata.ProbLo(); - const amrex::Real* dx = geomdata.CellSize(); - const amrex::Real x[AMREX_SPACEDIM] = {AMREX_D_DECL( - prob_lo[0] + static_cast(i + 0.5) * dx[0], - prob_lo[1] + static_cast(j + 0.5) * dx[1], - prob_lo[2] + static_cast(k + 0.5) * dx[2])}; - - // Create NUM_ODE boxes from (x_strt,y_strt) to (x_end,y_end) - for (int n = 0; n < NUM_ODE; n++) { - amrex::Real x_strt = prob_parm.ode_xy_lo + n * 4. * prob_parm.ode_length; - amrex::Real x_end = x_strt + prob_parm.ode_length; - amrex::Real y_strt = prob_parm.ode_xy_lo; - amrex::Real y_end = y_strt + prob_parm.ode_height; - if ((x[0] > x_strt) && (x[1] > y_strt)) { - if ((x[0] < x_end) && (x[1] < y_end)) { - ode_qty[n] = ode_qty_0[n]; - } + // Current x,y,z locations + const amrex::Real* prob_lo = geomdata.ProbLo(); + const amrex::Real* dx = geomdata.CellSize(); + const amrex::Real x[AMREX_SPACEDIM] = {AMREX_D_DECL( + prob_lo[0] + static_cast(i + 0.5) * dx[0], + prob_lo[1] + static_cast(j + 0.5) * dx[1], + prob_lo[2] + static_cast(k + 0.5) * dx[2])}; + + // Create NUM_ODE boxes from (x_strt,y_strt) to (x_end,y_end) + for (int n = 0; n < NUM_ODE; n++) { + amrex::Real x_strt = prob_parm.ode_xy_lo + n * 4. * prob_parm.ode_length; + amrex::Real x_end = x_strt + prob_parm.ode_length; + amrex::Real y_strt = prob_parm.ode_xy_lo; + amrex::Real y_end = y_strt + prob_parm.ode_height; + if ((x[0] > x_strt) && (x[1] > y_strt)) { + if ((x[0] < x_end) && (x[1] < y_end)) { + ode_qty[n] = ode_qty_0[n]; } } - // Add to state variable - for (int n = 0; n < NUM_ODE; n++) { - state(i, j, k, FIRSTODE + n) = ode_qty[n]; - } - */ -#endif } + // Add to state variable + for (int n = 0; n < NUM_ODE; n++) { + state(i, j, k, FIRSTODE + n) = ode_qty[n]; + } + */ +#endif } AMREX_GPU_DEVICE diff --git a/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob.cpp b/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob.cpp index 1d39f11bd..f4d737f0f 100644 --- a/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob.cpp +++ b/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob.cpp @@ -10,18 +10,5 @@ PeleLM::readProbParm() // NOLINT(readability-make-member-function-const) pp.query("P_mean", prob_parm->P_mean); pp.query("meanFlowDir", prob_parm->meanFlowDir); pp.query("meanFlowMag", prob_parm->meanFlowMag); - -#if NUM_ODE > 0 - amrex::Print() << "\n!!!!!! NUM_ODE = " << NUM_ODE << " !!!!!!!!!" << "\n"; -#endif - - // if (!m_incompressible) { - // auto& trans_parm = PeleLM::trans_parms.host_parm(); - // amrex::ParmParse pptr("transport"); - // pp.query("const_viscosity", trans_parm.const_viscosity); - // pp.query("const_bulk_viscosity", trans_parm.const_bulk_viscosity); - // pp.query("const_conductivity", trans_parm.const_conductivity); - // pp.query("const_diffusivity", trans_parm.const_diffusivity); - // PeleLM::trans_parms.sync_to_device(); - // } + pp.query("ode_IC", prob_parm->ode_IC); } diff --git a/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob_parm.H b/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob_parm.H index 34785b133..7a8e080e9 100644 --- a/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob_parm.H +++ b/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob_parm.H @@ -12,7 +12,7 @@ struct ProbParm : amrex::Gpu::Managed amrex::Real T_mean = 298.0_rt; amrex::Real P_mean = 101325.0_rt; amrex::Real meanFlowMag = 0.0; - amrex::Real ode_IC = 10.0; + amrex::Real ode_IC = 10.0_rt; amrex::Real ode_xy_lo = -0.015; amrex::Real ode_length = 0.005; amrex::Real ode_height = 0.010; diff --git a/Source/PeleLMeX_Advance.cpp b/Source/PeleLMeX_Advance.cpp index 110b86681..ab3df0f2d 100644 --- a/Source/PeleLMeX_Advance.cpp +++ b/Source/PeleLMeX_Advance.cpp @@ -132,7 +132,7 @@ PeleLM::Advance(int is_initIter) } #endif #if NUM_ODE > 0 - //computeODESource(AmrOldTime); + computeODESource(AmrOldTime); #endif if (m_incompressible == 0) { diff --git a/Source/PeleLMeX_Setup.cpp b/Source/PeleLMeX_Setup.cpp index 619c2cff4..b8a53fbb0 100644 --- a/Source/PeleLMeX_Setup.cpp +++ b/Source/PeleLMeX_Setup.cpp @@ -890,12 +890,6 @@ PeleLM::variablesSetup() m_AdvTypeState[FIRSTSOOT + mom] = 0; m_DiffTypeState[FIRSTSOOT + mom] = 0; } -#endif -#if NUM_ODE > 0 - for (int n = 0; n < NUM_ODE; n++) { - m_AdvTypeState[FIRSTODE + n] = 0; - m_DiffTypeState[FIRSTODE + n] = 0; - } #endif } diff --git a/Source/PeleLMeX_Utils.cpp b/Source/PeleLMeX_Utils.cpp index 0be2abcce..fc3c8402e 100644 --- a/Source/PeleLMeX_Utils.cpp +++ b/Source/PeleLMeX_Utils.cpp @@ -1393,6 +1393,12 @@ PeleLM::setTypicalValues(const TimeStamp& a_time, int is_init) #ifdef PELE_USE_EFIELD typical_values[NE] = 0.5 * (stateMax[NE] + stateMin[NE]); #endif +#if NUM_ODE > 0 + for (int n = 0; n < NUM_ODE; n++) { + typical_values[FIRSTODE + n] = + 0.5 * (stateMax[FIRSTODE + n] + stateMin[FIRSTODE + n]); + } +#endif // Pass into chemsitry if requested updateTypicalValuesChem(); @@ -1421,6 +1427,14 @@ PeleLM::setTypicalValues(const TimeStamp& a_time, int is_init) } #ifdef PELE_USE_EFIELD Print() << "\tnE: " << typical_values[NE] << '\n'; +#endif +#if NUM_ODE > 0 + for (int n = 0; n < NUM_ODE; n++) { + Print() << "\t" << m_ode_names[n] + << std::setw( + std::max(0, static_cast(10 - m_ode_names[n].length()))) + << std::left << ":" << typical_values[FIRSTODE + n] << '\n'; + } #endif } Print() << PrettyLine; From 27eff1fede2004ae8b08d1f4b2887370a62dc67a Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Tue, 1 Oct 2024 09:46:10 -0600 Subject: [PATCH 03/46] Minor edit in set_ode_names function. --- .../EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.H | 2 +- Exec/RegTests/EB_DecayingODEQTY/input.2d | 2 ++ Source/PeleLMeX_ProblemSpecificFunctions.H | 2 +- 3 files changed, 4 insertions(+), 2 deletions(-) diff --git a/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.H b/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.H index 57a6b490c..d96b9a4bb 100644 --- a/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.H +++ b/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.H @@ -8,8 +8,8 @@ using namespace amrex; void set_ode_names(amrex::Vector& a_ode_names) { - a_ode_names.resize(NUM_ODE); #if NUM_ODE > 0 + a_ode_names.resize(NUM_ODE); for (int n = 0; n < NUM_ODE; n++) { a_ode_names[n] = "MY_ODE_" + std::to_string(n); } diff --git a/Exec/RegTests/EB_DecayingODEQTY/input.2d b/Exec/RegTests/EB_DecayingODEQTY/input.2d index ca4f9a5ac..22bcb9db1 100644 --- a/Exec/RegTests/EB_DecayingODEQTY/input.2d +++ b/Exec/RegTests/EB_DecayingODEQTY/input.2d @@ -44,6 +44,8 @@ amr.cfl = 0.7 #-------------------------PELE CONTROLS---------------------- peleLM.v = 1 +#peleLM.num_init_iter = 0 +#peleLM.do_init_proj = 0 #---------------------------IO CONTROL--------------------------- #amr.restart = chk_07136 diff --git a/Source/PeleLMeX_ProblemSpecificFunctions.H b/Source/PeleLMeX_ProblemSpecificFunctions.H index 21b72b8dc..22375bc62 100644 --- a/Source/PeleLMeX_ProblemSpecificFunctions.H +++ b/Source/PeleLMeX_ProblemSpecificFunctions.H @@ -8,8 +8,8 @@ using namespace amrex; void set_ode_names(amrex::Vector& a_ode_names) { - a_ode_names.resize(NUM_ODE); #if NUM_ODE > 0 + a_ode_names.resize(NUM_ODE); for (int n = 0; n < NUM_ODE; n++) { a_ode_names[n] = "ODE_" + std::to_string(n); } From 7badae19cbda4f51bd602f9b66fd99a5442753b5 Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Tue, 1 Oct 2024 14:33:58 -0600 Subject: [PATCH 04/46] Fixed misalignment in plotting. IC's appear to be correct now for multiple ode vars. Still working on source terms. --- Exec/Make.PeleLMeX | 4 +++ Exec/RegTests/EB_DecayingODEQTY/GNUmakefile | 2 +- Exec/RegTests/EB_DecayingODEQTY/input.2d | 5 +++- .../EB_DecayingODEQTY/pelelmex_prob.H | 10 +++---- .../EB_DecayingODEQTY/pelelmex_prob.cpp | 3 +++ .../EB_DecayingODEQTY/pelelmex_prob_parm.H | 6 ++--- Source/PeleLMeX_Forces.cpp | 5 ++-- Source/PeleLMeX_Plot.cpp | 26 +++++++++++++++---- Source/PeleLMeX_Utils.cpp | 4 +-- 9 files changed, 46 insertions(+), 19 deletions(-) diff --git a/Exec/Make.PeleLMeX b/Exec/Make.PeleLMeX index 4ef1c7734..9f2cc32b8 100644 --- a/Exec/Make.PeleLMeX +++ b/Exec/Make.PeleLMeX @@ -61,6 +61,10 @@ ifeq ($(shell test 0$(NUM_ODE) -gt 0; echo $$?), 0) DEFINES += -DNUM_ODE=$(NUM_ODE) endif +ifeq ($(shell test 0$(NUM_AUX) -gt 0; echo $$?), 0) + DEFINES += -DNUM_AUX=$(NUM_AUX) +endif + Bpack += $(foreach dir, $(LMdirs), $(PELE_HOME)/$(dir)/Make.package) Blocs += $(foreach dir, $(LMdirs), $(PELE_HOME)/$(dir)) diff --git a/Exec/RegTests/EB_DecayingODEQTY/GNUmakefile b/Exec/RegTests/EB_DecayingODEQTY/GNUmakefile index 9bce8ffc7..87d4ab32a 100644 --- a/Exec/RegTests/EB_DecayingODEQTY/GNUmakefile +++ b/Exec/RegTests/EB_DecayingODEQTY/GNUmakefile @@ -28,7 +28,7 @@ THREAD_SANITIZER = FALSE # PeleLMeX CEXE_headers += EBUserDefined.H CEXE_headers += ProblemSpecificcFunctions.H -NUM_ODE = 1 +NUM_ODE = 3 # PelePhysics Chemistry_Model = air diff --git a/Exec/RegTests/EB_DecayingODEQTY/input.2d b/Exec/RegTests/EB_DecayingODEQTY/input.2d index 22bcb9db1..c0c1bdfd9 100644 --- a/Exec/RegTests/EB_DecayingODEQTY/input.2d +++ b/Exec/RegTests/EB_DecayingODEQTY/input.2d @@ -26,7 +26,10 @@ prob.T_mean = 300.0 prob.P_mean = 101325.0 prob.meanFlowMag = 4.255 prob.meanFlowDir = 1 -prob.ode_IC = 100.0 +prob.ode_IC = 10.0 +prob.ode_xy_lo = -0.01 +prob.ode_length = 0.01 +prob.ode_height = 0.02 #-----------------INPUTS TO CONSTANT TRANSPORT ------------------ transport.units = MKS diff --git a/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob.H b/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob.H index e1d2c4270..9a3e42f85 100644 --- a/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob.H +++ b/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob.H @@ -86,10 +86,10 @@ pelelmex_initdata( // Initial ODE quantity #if NUM_ODE > 0 - for (int n = 0; n < NUM_ODE; n++) { - state(i, j, k, FIRSTODE + n) = prob_parm.ode_IC; - } - /* + //for (int n = 0; n < NUM_ODE; n++) { + // state(i, j, k, FIRSTODE + n) = prob_parm.ode_IC; + //} + amrex::Real ode_qty[NUM_ODE] = {0.0}; // Initial state for ode quantities @@ -122,7 +122,7 @@ pelelmex_initdata( for (int n = 0; n < NUM_ODE; n++) { state(i, j, k, FIRSTODE + n) = ode_qty[n]; } - */ + #endif } diff --git a/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob.cpp b/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob.cpp index f4d737f0f..41411c4f9 100644 --- a/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob.cpp +++ b/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob.cpp @@ -11,4 +11,7 @@ PeleLM::readProbParm() // NOLINT(readability-make-member-function-const) pp.query("meanFlowDir", prob_parm->meanFlowDir); pp.query("meanFlowMag", prob_parm->meanFlowMag); pp.query("ode_IC", prob_parm->ode_IC); + pp.query("ode_xy_lo", prob_parm->ode_xy_lo); + pp.query("ode_length", prob_parm->ode_length); + pp.query("ode_height", prob_parm->ode_height); } diff --git a/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob_parm.H b/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob_parm.H index 7a8e080e9..4e8e73760 100644 --- a/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob_parm.H +++ b/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob_parm.H @@ -13,9 +13,9 @@ struct ProbParm : amrex::Gpu::Managed amrex::Real P_mean = 101325.0_rt; amrex::Real meanFlowMag = 0.0; amrex::Real ode_IC = 10.0_rt; - amrex::Real ode_xy_lo = -0.015; - amrex::Real ode_length = 0.005; - amrex::Real ode_height = 0.010; + amrex::Real ode_xy_lo = -0.01; + amrex::Real ode_length = 0.01; + amrex::Real ode_height = 0.020; int meanFlowDir = 1; }; #endif diff --git a/Source/PeleLMeX_Forces.cpp b/Source/PeleLMeX_Forces.cpp index 62c8dd429..91237dac9 100644 --- a/Source/PeleLMeX_Forces.cpp +++ b/Source/PeleLMeX_Forces.cpp @@ -236,7 +236,7 @@ PeleLM::addSpark(const TimeStamp& a_timestamp) } } - +#if NUM_ODE > 0 void PeleLM::computeODESource(const TimeStamp& a_timestamp) { @@ -254,4 +254,5 @@ PeleLM::computeODESource(const TimeStamp& a_timestamp) Gpu::streamSynchronize(); } } -} \ No newline at end of file +} +#endif \ No newline at end of file diff --git a/Source/PeleLMeX_Plot.cpp b/Source/PeleLMeX_Plot.cpp index 1b52a135b..a8a11488d 100644 --- a/Source/PeleLMeX_Plot.cpp +++ b/Source/PeleLMeX_Plot.cpp @@ -201,11 +201,6 @@ PeleLM::WritePlotFile() plt_VarsName.push_back(sootname); } #endif -#if NUM_ODE > 0 - for (int n = 0; n < NUM_ODE; n++) { - plt_VarsName.push_back(m_ode_names[n]); - } -#endif #ifdef PELE_USE_RADIATION if (do_rad_solve) { plt_VarsName.push_back("rad.G"); @@ -286,6 +281,20 @@ PeleLM::WritePlotFile() plt_VarsName.push_back("viscturb"); } +#if NUM_ODE > 0 + for (int n = 0; n < NUM_ODE; n++) { + plt_VarsName.push_back(m_ode_names[n]); + } +#endif + +Print() << PrettyLine; +Print() << "Plotting debugging" << std::endl; +Print() << "ncomp = " << ncomp << std::endl; +Print() << "plt_VarsName.size() = " << plt_VarsName.size() << std::endl; +for(int n = 0; n < plt_VarsName.size(); n++){ + Print() << "plt_VarsName[" << n << "] = " << plt_VarsName[n] << std::endl; +} + //---------------------------------------------------------------- // Fill the plot MultiFabs for (int lev = 0; lev <= finest_level; ++lev) { @@ -406,6 +415,13 @@ PeleLM::WritePlotFile() cnt += m_ionsFluxes[lev]->nComp(); } #endif +#if NUM_ODE > 0 + MultiFab::Copy( + mf_plt[lev], m_leveldata_new[lev]->state, FIRSTODE, cnt, NUM_ODE, + 0); + cnt += NUM_ODE; + Print() << "cnt = " << cnt << std::endl; +#endif if (m_do_les && m_plot_les) { constexpr amrex::Real fact = 0.5 / AMREX_SPACEDIM; diff --git a/Source/PeleLMeX_Utils.cpp b/Source/PeleLMeX_Utils.cpp index fc3c8402e..0ec239cb6 100644 --- a/Source/PeleLMeX_Utils.cpp +++ b/Source/PeleLMeX_Utils.cpp @@ -1430,10 +1430,10 @@ PeleLM::setTypicalValues(const TimeStamp& a_time, int is_init) #endif #if NUM_ODE > 0 for (int n = 0; n < NUM_ODE; n++) { - Print() << "\t" << m_ode_names[n] + Print() << "\t" << m_ode_names[n] << std::setw( std::max(0, static_cast(10 - m_ode_names[n].length()))) - << std::left << ":" << typical_values[FIRSTODE + n] << '\n'; + << std::left << ":" << typical_values[FIRSTODE + n] <<'\n'; } #endif } From 9376aa911ecf71ed583516c7e44607e48a78157e Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Mon, 7 Oct 2024 15:33:08 -0600 Subject: [PATCH 05/46] Fix indexing. --- Source/PeleLMeX_Index.H | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/Source/PeleLMeX_Index.H b/Source/PeleLMeX_Index.H index aa71ea5ae..2d9ce9b95 100644 --- a/Source/PeleLMeX_Index.H +++ b/Source/PeleLMeX_Index.H @@ -22,14 +22,13 @@ #else #define NUMSOOTVAR 0 #endif - #ifndef NUM_ODE #define NUM_ODE 0 #endif #define FIRSTODE (RHORT + NUMEFIELD + NUMSOOTVAR + 1) - -#define FIRSTAUX (FIRSTODE + NUM_ODE + 1) -#define NVAR (FIRSTAUX + 0) +#define FIRSTAUX (RHORT + NUMEFIELD + NUMSOOTVAR + NUM_ODE + 1) +#define NUM_AUX 0 +#define NVAR (RHORT + NUMEFIELD + NUMSOOTVAR + NUM_ODE + NUM_AUX + 1) // TODO: handle NAUX #endif From 0751a38822187e0268188068e8b323c548c28349 Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Mon, 7 Oct 2024 22:00:29 -0600 Subject: [PATCH 06/46] Updated problem_modify_ext_sources. Works with local case definition. --- .../PeleLMeX_ProblemSpecificFunctions.H | 34 ++++++++++++-- Exec/RegTests/EB_DecayingODEQTY/input.2d | 2 +- Source/PeleLMeX.H | 4 -- Source/PeleLMeX_Advance.cpp | 28 ++++++++++-- Source/PeleLMeX_Forces.cpp | 23 +--------- Source/PeleLMeX_ProblemSpecificFunctions.H | 44 +++++++++++++------ 6 files changed, 88 insertions(+), 47 deletions(-) diff --git a/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.H b/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.H index d96b9a4bb..fba5c167c 100644 --- a/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.H +++ b/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.H @@ -1,12 +1,12 @@ #ifndef PROBLEMSPECIFICFUNCTIONS_H #define PROBLEMSPECIFICFUNCTIONS_H -using namespace amrex; - #include #include -void set_ode_names(amrex::Vector& a_ode_names) +using namespace amrex; + +static void set_ode_names(Vector& a_ode_names) { #if NUM_ODE > 0 a_ode_names.resize(NUM_ODE); @@ -16,4 +16,32 @@ void set_ode_names(amrex::Vector& a_ode_names) #endif } +static void problem_modify_ext_sources( + Real /*time*/, + Real /*dt*/, + int lev, + MultiArray4 const& state_old, + MultiArray4 const& state_new, + Vector>& a_extSource, + const GeometryData& /*geomdata*/, + ProbParm const& prob_parm) +{ + /* Notes: a_extSource contains sources from velocity forcing coming in + This function should add to rather than overwrite a_extSource. + + TODO: Add and test capabilities with species, velocity and enthalpy. + */ + + auto ext_source_arr = a_extSource[lev]->arrays(); + amrex::ParallelFor( + *a_extSource[lev], + [=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept { + for (int n = 0; n < NUM_ODE; n++){ + amrex::Real B_n = state_old[box_no](i, j, k, FIRSTODE + n); + ext_source_arr[box_no](i, j, k, FIRSTODE + n) += -1000.*B_n; + } + }); + amrex::Gpu::streamSynchronize(); +} + #endif \ No newline at end of file diff --git a/Exec/RegTests/EB_DecayingODEQTY/input.2d b/Exec/RegTests/EB_DecayingODEQTY/input.2d index c0c1bdfd9..1962974be 100644 --- a/Exec/RegTests/EB_DecayingODEQTY/input.2d +++ b/Exec/RegTests/EB_DecayingODEQTY/input.2d @@ -11,7 +11,7 @@ peleLM.lo_bc = Inflow Interior peleLM.hi_bc = Outflow Interior #-------------------------AMR CONTROL---------------------------- -amr.n_cell = 192 64 # Level 0 number of cells in each direction +amr.n_cell = 96 32 # Level 0 number of cells in each direction amr.v = 1 # AMR verbose amr.max_level = 0 # maximum level number allowed amr.ref_ratio = 2 2 2 2 # refinement ratio diff --git a/Source/PeleLMeX.H b/Source/PeleLMeX.H index 5a40e6d3f..13c839587 100644 --- a/Source/PeleLMeX.H +++ b/Source/PeleLMeX.H @@ -874,10 +874,6 @@ public: void computeRadSource(const PeleLM::TimeStamp& a_timestamp); #endif -#if NUM_ODE > 0 - // Compute the ode source term - void computeODESource(const PeleLM::TimeStamp& a_timestamp); -#endif //----------------------------------------------------------------------------- // EOS diff --git a/Source/PeleLMeX_Advance.cpp b/Source/PeleLMeX_Advance.cpp index ab3df0f2d..2bf9fadb2 100644 --- a/Source/PeleLMeX_Advance.cpp +++ b/Source/PeleLMeX_Advance.cpp @@ -1,5 +1,6 @@ #include #include +#include #include #include @@ -131,9 +132,30 @@ PeleLM::Advance(int is_initIter) BL_PROFILE_VAR_STOP(PLM_RAD); } #endif -#if NUM_ODE > 0 - computeODESource(AmrOldTime); -#endif + +// Additional user defined source terms +for (int lev = 0; lev <= finest_level; lev++) { + problem_modify_ext_sources(getTime(lev, AmrNewTime), m_dt, lev, + getLevelDataPtr(lev, AmrOldTime)->state.const_arrays(), + getLevelDataPtr(lev, AmrNewTime)->state.const_arrays(), + m_extSource, geom[lev].data(), *prob_parm_d); + + // !Test the call for modify ext_sources + /* + auto extma = m_extSource[lev]->arrays(); + amrex::ParallelFor( + *m_extSource[lev], + [=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept { + for(int n = 0; n < NUM_ODE; n++){ + if(extma[box_no](i, j, k, FIRSTODE + n) < 0){ + Print() << "extma[lev = "<ode_xy_lo); pp.query("ode_length", prob_parm->ode_length); pp.query("ode_height", prob_parm->ode_height); + pp.query("ode_srcstrength", prob_parm->ode_srcstrength); } diff --git a/Source/PeleLMeX_ProblemSpecificFunctions.H b/Source/PeleLMeX_ProblemSpecificFunctions.H index 8851c859f..804725484 100644 --- a/Source/PeleLMeX_ProblemSpecificFunctions.H +++ b/Source/PeleLMeX_ProblemSpecificFunctions.H @@ -40,8 +40,7 @@ static void problem_modify_ext_sources( 2) Requires compilation with USE_MODIFIED_SOURCES = TRUE TODO: - 1) Test geomdata, prob_parm - 2) Add and test capabilities with species, velocity and enthalpy. + 1) Add and test capabilities with species, velocity and enthalpy. // Example: adding a exponential decay source term to an ode_qty auto ext_source_arr = a_extSource[lev]->arrays(); From 26370253be8646c559b7e6e099a3f7bd0391773c Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Fri, 11 Oct 2024 09:23:17 -0600 Subject: [PATCH 10/46] Added ode_srcstrength to test case, fixed minor bugs. --- .../EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.H | 4 ++-- Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob_parm.H | 3 ++- Source/PeleLMeX_Advance.cpp | 7 +++---- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.H b/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.H index 7a18b097c..774b26dbc 100644 --- a/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.H +++ b/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.H @@ -31,7 +31,7 @@ static void problem_modify_ext_sources( MultiArray4 const& /*state_new*/, Vector>& a_extSource, const GeometryData& /*geomdata*/, - ProbParm const& /*prob_parm*/) + ProbParm const& prob_parm) { /* Notes: @@ -49,7 +49,7 @@ static void problem_modify_ext_sources( [=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept { for (int n = 0; n < NUM_ODE; n++){ amrex::Real B_n = state_old[box_no](i, j, k, FIRSTODE + n); - ext_source_arr[box_no](i, j, k, FIRSTODE + n) += prob_parm.ode_srcstrength.*B_n; + ext_source_arr[box_no](i, j, k, FIRSTODE + n) += prob_parm.ode_srcstrength*B_n; } }); amrex::Gpu::streamSynchronize(); diff --git a/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob_parm.H b/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob_parm.H index 4e8e73760..b74ae6ce2 100644 --- a/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob_parm.H +++ b/Exec/RegTests/EB_DecayingODEQTY/pelelmex_prob_parm.H @@ -12,10 +12,11 @@ struct ProbParm : amrex::Gpu::Managed amrex::Real T_mean = 298.0_rt; amrex::Real P_mean = 101325.0_rt; amrex::Real meanFlowMag = 0.0; - amrex::Real ode_IC = 10.0_rt; + amrex::Real ode_IC = 10.0; amrex::Real ode_xy_lo = -0.01; amrex::Real ode_length = 0.01; amrex::Real ode_height = 0.020; + amrex::Real ode_srcstrength = -100.0; int meanFlowDir = 1; }; #endif diff --git a/Source/PeleLMeX_Advance.cpp b/Source/PeleLMeX_Advance.cpp index 7cb8a1c28..9dd65ca25 100644 --- a/Source/PeleLMeX_Advance.cpp +++ b/Source/PeleLMeX_Advance.cpp @@ -133,7 +133,7 @@ PeleLM::Advance(int is_initIter) } #endif -// Additional user defined source terms +// Additional user defined source terms (only for ode qty currently) #ifdef PELE_USE_MODIFIED_SOURCES for (int lev = 0; lev <= finest_level; lev++) { problem_modify_ext_sources(getTime(lev, AmrNewTime), m_dt, lev, @@ -141,7 +141,7 @@ for (int lev = 0; lev <= finest_level; lev++) { getLevelDataPtr(lev, AmrNewTime)->state.const_arrays(), m_extSource, geom[lev].data(), *prob_parm_d); - // !Test the call for modify ext_sources + // Debugging: Test the call for modify ext_sources auto ext_src = m_extSource[lev]->arrays(); amrex::ParallelFor( *m_extSource[lev], @@ -152,8 +152,7 @@ for (int lev = 0; lev <= finest_level; lev++) { } } }); - - // !End of test + // Debugging: End of test } #endif From 59c987a9a6843538dc757c6ca2f6ce28d5cf30ee Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Fri, 11 Oct 2024 10:38:42 -0600 Subject: [PATCH 11/46] Changed flag for user defined external sources. Added future flag for including external sources in SDC loop. --- Exec/Make.PeleLMeX | 4 -- Exec/RegTests/EB_DecayingODEQTY/GNUmakefile | 1 - .../PeleLMeX_ProblemSpecificFunctions.H | 2 +- Exec/RegTests/EB_DecayingODEQTY/input.2d | 1 + Source/PeleLMeX.H | 6 +++ Source/PeleLMeX_Advance.cpp | 44 +++++++++---------- Source/PeleLMeX_ProblemSpecificFunctions.H | 2 +- Source/PeleLMeX_Setup.cpp | 7 +++ 8 files changed, 38 insertions(+), 29 deletions(-) diff --git a/Exec/Make.PeleLMeX b/Exec/Make.PeleLMeX index 17f1ae2eb..9f2cc32b8 100644 --- a/Exec/Make.PeleLMeX +++ b/Exec/Make.PeleLMeX @@ -65,10 +65,6 @@ ifeq ($(shell test 0$(NUM_AUX) -gt 0; echo $$?), 0) DEFINES += -DNUM_AUX=$(NUM_AUX) endif -ifeq ($(USE_MODIFIED_SOURCES), TRUE) - DEFINES += -DPELE_USE_MODIFIED_SOURCES -endif - Bpack += $(foreach dir, $(LMdirs), $(PELE_HOME)/$(dir)/Make.package) Blocs += $(foreach dir, $(LMdirs), $(PELE_HOME)/$(dir)) diff --git a/Exec/RegTests/EB_DecayingODEQTY/GNUmakefile b/Exec/RegTests/EB_DecayingODEQTY/GNUmakefile index 80bad1dfc..fb732dfe1 100644 --- a/Exec/RegTests/EB_DecayingODEQTY/GNUmakefile +++ b/Exec/RegTests/EB_DecayingODEQTY/GNUmakefile @@ -28,7 +28,6 @@ THREAD_SANITIZER = FALSE # PeleLMeX CEXE_headers += EBUserDefined.H CEXE_headers += ProblemSpecificFunctions.H -USE_MODIFIED_SOURCES = TRUE NUM_ODE = 3 # PelePhysics diff --git a/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.H b/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.H index 774b26dbc..3bc9a3960 100644 --- a/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.H +++ b/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.H @@ -37,7 +37,7 @@ static void problem_modify_ext_sources( Notes: 1) a_extSource contains sources from velocity forcing coming in. This function should add to rather than overwrite a_extSource. - 2) Requires compilation with USE_MODIFIED_SOURCES = TRUE + 2) Requires peleLM.user_defined_ext_sources = true in input file TODO: 1) Add and test capabilities with species, velocity and enthalpy. diff --git a/Exec/RegTests/EB_DecayingODEQTY/input.2d b/Exec/RegTests/EB_DecayingODEQTY/input.2d index 1962974be..223f100b4 100644 --- a/Exec/RegTests/EB_DecayingODEQTY/input.2d +++ b/Exec/RegTests/EB_DecayingODEQTY/input.2d @@ -47,6 +47,7 @@ amr.cfl = 0.7 #-------------------------PELE CONTROLS---------------------- peleLM.v = 1 +peleLM.user_defined_ext_sources = 0 #peleLM.num_init_iter = 0 #peleLM.do_init_proj = 0 diff --git a/Source/PeleLMeX.H b/Source/PeleLMeX.H index 13c839587..173c2b785 100644 --- a/Source/PeleLMeX.H +++ b/Source/PeleLMeX.H @@ -874,6 +874,12 @@ public: void computeRadSource(const PeleLM::TimeStamp& a_timestamp); #endif + //----------------------------------------------------------------------------- + // External Sources + bool m_user_defined_ext_sources; + bool m_ext_sources_SDC; + + //----------------------------------------------------------------------------- // EOS diff --git a/Source/PeleLMeX_Advance.cpp b/Source/PeleLMeX_Advance.cpp index 9dd65ca25..bc5427ff6 100644 --- a/Source/PeleLMeX_Advance.cpp +++ b/Source/PeleLMeX_Advance.cpp @@ -133,28 +133,28 @@ PeleLM::Advance(int is_initIter) } #endif -// Additional user defined source terms (only for ode qty currently) -#ifdef PELE_USE_MODIFIED_SOURCES -for (int lev = 0; lev <= finest_level; lev++) { - problem_modify_ext_sources(getTime(lev, AmrNewTime), m_dt, lev, - getLevelDataPtr(lev, AmrOldTime)->state.const_arrays(), - getLevelDataPtr(lev, AmrNewTime)->state.const_arrays(), - m_extSource, geom[lev].data(), *prob_parm_d); - - // Debugging: Test the call for modify ext_sources - auto ext_src = m_extSource[lev]->arrays(); - amrex::ParallelFor( - *m_extSource[lev], - [=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept { - for(int n = 0; n < NUM_ODE; n++){ - if(ext_src[box_no](i, j, k, FIRSTODE + n) < 0){ - Print() << "ext_src[lev = "< Date: Fri, 11 Oct 2024 10:39:49 -0600 Subject: [PATCH 12/46] Set peleLM.user_defined_ext_sources to true in input file. --- Exec/RegTests/EB_DecayingODEQTY/input.2d | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Exec/RegTests/EB_DecayingODEQTY/input.2d b/Exec/RegTests/EB_DecayingODEQTY/input.2d index 223f100b4..1a5093c9c 100644 --- a/Exec/RegTests/EB_DecayingODEQTY/input.2d +++ b/Exec/RegTests/EB_DecayingODEQTY/input.2d @@ -47,7 +47,7 @@ amr.cfl = 0.7 #-------------------------PELE CONTROLS---------------------- peleLM.v = 1 -peleLM.user_defined_ext_sources = 0 +peleLM.user_defined_ext_sources = 1 #peleLM.num_init_iter = 0 #peleLM.do_init_proj = 0 From 561e4a8fa1d044a466218744a1383495476e05b3 Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Fri, 11 Oct 2024 15:01:52 -0600 Subject: [PATCH 13/46] Added Euler step for advancing the ODE quantities in time. --- .../PeleLMeX_ProblemSpecificFunctions.H | 6 ++--- Source/Make.package | 4 +++ Source/PeleLMeX.H | 1 + Source/PeleLMeX_Advance.cpp | 15 ++++++++--- Source/PeleLMeX_ODEQty.cpp | 25 +++++++++++++++++++ Source/PeleLMeX_Plot.cpp | 13 +++++----- 6 files changed, 52 insertions(+), 12 deletions(-) create mode 100644 Source/PeleLMeX_ODEQty.cpp diff --git a/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.H b/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.H index 3bc9a3960..5df91a8f5 100644 --- a/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.H +++ b/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.H @@ -27,8 +27,8 @@ static void problem_modify_ext_sources( Real /*time*/, Real /*dt*/, int lev, - MultiArray4 const& state_old, - MultiArray4 const& /*state_new*/, + MultiArray4 const& state_old_arr, + MultiArray4 const& /*state_new_arr*/, Vector>& a_extSource, const GeometryData& /*geomdata*/, ProbParm const& prob_parm) @@ -48,7 +48,7 @@ static void problem_modify_ext_sources( *a_extSource[lev], [=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept { for (int n = 0; n < NUM_ODE; n++){ - amrex::Real B_n = state_old[box_no](i, j, k, FIRSTODE + n); + amrex::Real B_n = state_old_arr[box_no](i, j, k, FIRSTODE + n); ext_source_arr[box_no](i, j, k, FIRSTODE + n) += prob_parm.ode_srcstrength*B_n; } }); diff --git a/Source/Make.package b/Source/Make.package index eca6d42cd..de53387d3 100644 --- a/Source/Make.package +++ b/Source/Make.package @@ -51,6 +51,10 @@ ifeq ($(USE_SOOT), TRUE) CEXE_sources += PeleLMeX_Soot.cpp endif +ifeq ($(shell test 0$(NUM_ODE) -gt 0; echo $$?), 0) + CEXE_sources += PeleLMeX_ODEQty.cpp +endif + ifeq ($(USE_RADIATION), TRUE) CEXE_sources += PeleLMeX_Radiation.cpp endif diff --git a/Source/PeleLMeX.H b/Source/PeleLMeX.H index 173c2b785..4c90b2e4b 100644 --- a/Source/PeleLMeX.H +++ b/Source/PeleLMeX.H @@ -1865,6 +1865,7 @@ public: // ODE variable names #if NUM_ODE > 0 amrex::Vector m_ode_names; + void advanceODEQty(); #endif // Aux data size diff --git a/Source/PeleLMeX_Advance.cpp b/Source/PeleLMeX_Advance.cpp index bc5427ff6..8b1a4975d 100644 --- a/Source/PeleLMeX_Advance.cpp +++ b/Source/PeleLMeX_Advance.cpp @@ -147,9 +147,9 @@ PeleLM::Advance(int is_initIter) *m_extSource[lev], [=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept { for(int n = 0; n < NUM_ODE; n++){ - if(ext_src[box_no](i, j, k, FIRSTODE + n) < 0){ - Print() << "ext_src[lev = "< Date: Mon, 14 Oct 2024 11:36:43 -0600 Subject: [PATCH 16/46] Added .cpp file for ProblemSpecificFunctions to avoid need for static functions. --- CMake/BuildPeleExe.cmake | 1 + Exec/RegTests/EB_DecayingODEQTY/GNUmakefile | 2 +- ... => PeleLMeX_ProblemSpecificFunctions.cpp} | 18 +++--- Source/Make.package | 1 + Source/PeleLMeX_ProblemSpecificFunctions.H | 62 ++++--------------- Source/PeleLMeX_ProblemSpecificFunctions.cpp | 55 ++++++++++++++++ 6 files changed, 76 insertions(+), 63 deletions(-) rename Exec/RegTests/EB_DecayingODEQTY/{PeleLMeX_ProblemSpecificFunctions.H => PeleLMeX_ProblemSpecificFunctions.cpp} (80%) create mode 100644 Source/PeleLMeX_ProblemSpecificFunctions.cpp diff --git a/CMake/BuildPeleExe.cmake b/CMake/BuildPeleExe.cmake index 1be05d024..46e4114b3 100644 --- a/CMake/BuildPeleExe.cmake +++ b/CMake/BuildPeleExe.cmake @@ -51,6 +51,7 @@ function(build_pele_exe pele_exe_name pele_physics_lib_name) ${SRC_DIR}/PeleLMeX_Init.cpp ${SRC_DIR}/PeleLMeX_Plot.cpp ${SRC_DIR}/PeleLMeX_ProblemSpecificFunctions.H + ${SRC_DIR}/PeleLMeX_ProblemSpecificFunctions.cpp ${SRC_DIR}/PeleLMeX_Projection.cpp ${SRC_DIR}/PeleLMeX_Reactions.cpp ${SRC_DIR}/PeleLMeX_Regrid.cpp diff --git a/Exec/RegTests/EB_DecayingODEQTY/GNUmakefile b/Exec/RegTests/EB_DecayingODEQTY/GNUmakefile index fb732dfe1..e97efd845 100644 --- a/Exec/RegTests/EB_DecayingODEQTY/GNUmakefile +++ b/Exec/RegTests/EB_DecayingODEQTY/GNUmakefile @@ -27,7 +27,7 @@ THREAD_SANITIZER = FALSE # PeleLMeX CEXE_headers += EBUserDefined.H -CEXE_headers += ProblemSpecificFunctions.H +CEXE_sources += PeleLMeX_ProblemSpecificFunctions.cpp NUM_ODE = 3 # PelePhysics diff --git a/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.H b/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.cpp similarity index 80% rename from Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.H rename to Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.cpp index 5df91a8f5..6ed5497e4 100644 --- a/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.H +++ b/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.cpp @@ -1,19 +1,17 @@ -#ifndef PROBLEMSPECIFICFUNCTIONS_H -#define PROBLEMSPECIFICFUNCTIONS_H - -#include #include +#include +#include using namespace amrex; /* Problem specific functions: - This file must be copied locally to the case directory -- Add the following to GNUmakefile: CEXE_headers += ProblemSpecificFunctions.H -- Only copy the functions you want to use +- Add the following to GNUmakefile: CEXE_sources += PeleLMeX_ProblemSpecificFunctions.cpp */ -static void set_ode_names(Vector& a_ode_names) + +void set_ode_names(Vector& a_ode_names) { #if NUM_ODE > 0 a_ode_names.resize(NUM_ODE); @@ -23,7 +21,7 @@ static void set_ode_names(Vector& a_ode_names) #endif } -static void problem_modify_ext_sources( +void problem_modify_ext_sources( Real /*time*/, Real /*dt*/, int lev, @@ -53,6 +51,4 @@ static void problem_modify_ext_sources( } }); amrex::Gpu::streamSynchronize(); -} - -#endif \ No newline at end of file +} \ No newline at end of file diff --git a/Source/Make.package b/Source/Make.package index de53387d3..821172eff 100644 --- a/Source/Make.package +++ b/Source/Make.package @@ -46,6 +46,7 @@ CEXE_sources += PeleLMeX_FlowController.cpp CEXE_sources += PeleLMeX_DeriveUserDefined.cpp CEXE_sources += PeleLMeX_BPatch.cpp CEXE_sources += PeleLMeX_PatchFlowVariables.cpp +CEXE_sources += PeleLMeX_ProblemSpecificFunctions.cpp ifeq ($(USE_SOOT), TRUE) CEXE_sources += PeleLMeX_Soot.cpp diff --git a/Source/PeleLMeX_ProblemSpecificFunctions.H b/Source/PeleLMeX_ProblemSpecificFunctions.H index a020a5907..aeb92356c 100644 --- a/Source/PeleLMeX_ProblemSpecificFunctions.H +++ b/Source/PeleLMeX_ProblemSpecificFunctions.H @@ -1,59 +1,19 @@ #ifndef PROBLEMSPECIFICFUNCTIONS_H #define PROBLEMSPECIFICFUNCTIONS_H -#include -#include - using namespace amrex; -/* -Problem specific functions: -- This file must be copied locally to the case directory -- Add the following to GNUmakefile: CEXE_headers += ProblemSpecificFunctions.H -- Only copy the functions you want to use -*/ - -static void set_ode_names(Vector& a_ode_names) -{ -#if NUM_ODE > 0 - a_ode_names.resize(NUM_ODE); - for (int n = 0; n < NUM_ODE; n++) { - a_ode_names[n] = "ODE_" + std::to_string(n); - } -#endif -} - -static void problem_modify_ext_sources( - Real /*time*/, - Real /*dt*/, - int /*lev*/, - MultiArray4 const& /*state_old*/, - MultiArray4 const& /*state_new*/, - Vector>& /*a_extSource*/, - const GeometryData& /*geomdata*/, - ProbParm const& /*prob_parm*/) -{ - /* - Notes: - 1) a_extSource contains sources from velocity forcing coming in. - This function should add to rather than overwrite a_extSource. - 2) Requires peleLM.user_defined_ext_sources = true in input file - - TODO: - 1) Add and test capabilities with species, velocity and enthalpy. +class PeleLM; - // Example: adding a exponential decay source term to an ode_qty - auto ext_source_arr = a_extSource[lev]->arrays(); - amrex::ParallelFor( - *a_extSource[lev], - [=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept { - for (int n = 0; n < NUM_ODE; n++){ - amrex::Real B_n = state_old[box_no](i, j, k, FIRSTODE + n); - ext_source_arr[box_no](i, j, k, FIRSTODE + n) += -1000.*B_n; - } - }); - amrex::Gpu::streamSynchronize(); - */ -} +void set_ode_names(Vector& a_ode_names); +void problem_modify_ext_sources( + Real time, + Real dt, + int lev, + MultiArray4 const& state_old, + MultiArray4 const& state_new, + Vector>& a_extSource, + const GeometryData& geomdata, + ProbParm const& prob_parm); #endif \ No newline at end of file diff --git a/Source/PeleLMeX_ProblemSpecificFunctions.cpp b/Source/PeleLMeX_ProblemSpecificFunctions.cpp new file mode 100644 index 000000000..edd61cfb8 --- /dev/null +++ b/Source/PeleLMeX_ProblemSpecificFunctions.cpp @@ -0,0 +1,55 @@ +#include +#include +#include + +using namespace amrex; + +/* +Problem specific functions: +- This file must be copied in its entirety locally +- Add the following to GNUmakefile: CEXE_sources += PeleLMeX_ProblemSpecificFunctions.cpp +- Modify functions as needed +*/ + +void set_ode_names(Vector& a_ode_names) +{ +#if NUM_ODE > 0 + a_ode_names.resize(NUM_ODE); + for (int n = 0; n < NUM_ODE; n++) { + a_ode_names[n] = "ODE_" + std::to_string(n); + } +#endif +} + +static void problem_modify_ext_sources( + Real /*time*/, + Real /*dt*/, + int /*lev*/, + MultiArray4 const& /*state_old*/, + MultiArray4 const& /*state_new*/, + Vector>& /*a_extSource*/, + const GeometryData& /*geomdata*/, + ProbParm const& /*prob_parm*/) +{ + /* + Notes: + 1) a_extSource contains sources from velocity forcing coming in. + This function should add to rather than overwrite a_extSource. + 2) Requires peleLM.user_defined_ext_sources = true in input file + + TODO: + 1) Add and test capabilities with species, velocity and enthalpy. + + // Example: adding a exponential decay source term to an ode_qty + auto ext_source_arr = a_extSource[lev]->arrays(); + amrex::ParallelFor( + *a_extSource[lev], + [=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept { + for (int n = 0; n < NUM_ODE; n++){ + amrex::Real B_n = state_old[box_no](i, j, k, FIRSTODE + n); + ext_source_arr[box_no](i, j, k, FIRSTODE + n) += -10.*B_n; + } + }); + amrex::Gpu::streamSynchronize(); + */ +} \ No newline at end of file From 283b8b159e8d832e1ce8f2c2475537bc197bf01d Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Mon, 14 Oct 2024 21:02:32 -0600 Subject: [PATCH 17/46] Added convergence test for ODE Qty. --- .../PeleLMeX_ProblemSpecificFunctions.cpp | 1 + Exec/RegTests/EB_DecayingODEQTY/input.2d | 3 +- .../EB_DecayingODEQTY/input.convergence_0 | 79 +++++++++++ .../EB_DecayingODEQTY/run_convergence_test.py | 128 ++++++++++++++++++ Source/PeleLMeX.H | 1 + Source/PeleLMeX_ODEQty.cpp | 5 + Source/PeleLMeX_ProblemSpecificFunctions.cpp | 4 +- 7 files changed, 218 insertions(+), 3 deletions(-) create mode 100644 Exec/RegTests/EB_DecayingODEQTY/input.convergence_0 create mode 100644 Exec/RegTests/EB_DecayingODEQTY/run_convergence_test.py diff --git a/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.cpp b/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.cpp index 6ed5497e4..02a7f8fad 100644 --- a/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.cpp +++ b/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.cpp @@ -8,6 +8,7 @@ using namespace amrex; Problem specific functions: - This file must be copied locally to the case directory - Add the following to GNUmakefile: CEXE_sources += PeleLMeX_ProblemSpecificFunctions.cpp +- Modify as needed */ diff --git a/Exec/RegTests/EB_DecayingODEQTY/input.2d b/Exec/RegTests/EB_DecayingODEQTY/input.2d index 9faf145c6..ebf728a61 100644 --- a/Exec/RegTests/EB_DecayingODEQTY/input.2d +++ b/Exec/RegTests/EB_DecayingODEQTY/input.2d @@ -42,9 +42,10 @@ transport.const_diffusivity = 0.0 #---------------------------TIME STEPPING------------------------ #amr.max_step = 1 amr.stop_time = 0.06 +#amr.fixed_dt = 0.0001 amr.dt_shrink = 0.1 amr.dt_change_max = 1.1 -amr.cfl = 0.7 +amr.cfl = 0.9 #-------------------------PELE CONTROLS---------------------- peleLM.v = 1 diff --git a/Exec/RegTests/EB_DecayingODEQTY/input.convergence_0 b/Exec/RegTests/EB_DecayingODEQTY/input.convergence_0 new file mode 100644 index 000000000..d3d633b5c --- /dev/null +++ b/Exec/RegTests/EB_DecayingODEQTY/input.convergence_0 @@ -0,0 +1,79 @@ +#----------------------DOMAIN DEFINITION------------------------- +geometry.is_periodic = 0 1 # For each dir, 0: non-perio, 1: periodic +geometry.coord_sys = 0 # 0 => cart, 1 => RZ +geometry.prob_lo = -0.02 -0.02 # x_lo y_lo (z_lo) +geometry.prob_hi = 0.10 0.02 # x_hi y_hi (z_hi) + +#--------------------------BC FLAGS------------------------------ +# Interior, Inflow, Outflow, Symmetry, +# SlipWallAdiab, NoSlipWallAdiab, SlipWallIsotherm, NoSlipWallIsotherm +peleLM.lo_bc = Inflow Interior +peleLM.hi_bc = Outflow Interior + +#-------------------------AMR CONTROL---------------------------- +amr.n_cell = 96 32 # Level 0 number of cells in each direction +amr.v = 1 # AMR verbose +amr.max_level = 0 # maximum level number allowed +amr.ref_ratio = 2 2 2 2 # refinement ratio +amr.regrid_int = 5 # how often to regrid +amr.n_error_buf = 2 2 2 2 # number of buffer cells in error est +amr.grid_eff = 0.7 # what constitutes an efficient grid +amr.blocking_factor = 16 # block factor in grid generation (min box size) +amr.max_grid_size = 64 # max box size + +#--------------------------- Problem ---------------------------- +prob.T_mean = 300.0 +prob.P_mean = 101325.0 +prob.meanFlowMag = 4.255 +prob.meanFlowDir = 1 +prob.ode_IC = 10.0 +prob.ode_xy_lo = -0.01 +prob.ode_length = 0.01 +prob.ode_height = 0.02 +prob.ode_srcstrength = -100.0 + +#-----------------INPUTS TO CONSTANT TRANSPORT ------------------ +transport.units = MKS +transport.const_viscosity = 0.0001 +transport.const_bulk_viscosity = 0.0 +transport.const_conductivity = 0.0 +transport.const_diffusivity = 0.0 + +#---------------------------TIME STEPPING------------------------ +amr.stop_time = 0.05 +amr.fixed_dt = 0.001 + +#-------------------------PELE CONTROLS---------------------- +peleLM.v = 1 +peleLM.user_defined_ext_sources = 1 +#peleLM.num_init_iter = 0 +#peleLM.do_init_proj = 0 + +#---------------------------IO CONTROL--------------------------- +#amr.restart = chk_07136 +amr.check_file = "chk_" +amr.check_per = -1 +amr.plot_file = "plt_conv0_" +amr.plot_per = 0.001 +amr.derive_plot_vars = avg_pressure mag_vort + +#---------------------- Temporal CONTROL ------------------------- +peleLM.do_temporals = 1 # Turn temporals ON/OFF +peleLM.temporal_int = 1 +peleLM.do_extremas = 1 # Compute state extremas + +#------------------------- EB SETUP ----------------------------- +eb2.geom_type = sphere +eb2.sphere_radius = 0.005 +eb2.sphere_center = 0.03 0.0 +eb2.sphere_has_fluid_inside = 0 +eb2.small_volfrac = 1.0e-4 + +#--------------------REFINEMENT CONTROL-------------------------- +amr.refinement_indicators = VortL VortH +amr.VortL.max_level = 1 +amr.VortL.value_less = -1000 +amr.VortL.field_name = mag_vort +amr.VortH.max_level = 1 +amr.VortH.value_greater = 1000 +amr.VortH.field_name = mag_vort diff --git a/Exec/RegTests/EB_DecayingODEQTY/run_convergence_test.py b/Exec/RegTests/EB_DecayingODEQTY/run_convergence_test.py new file mode 100644 index 000000000..adc93560e --- /dev/null +++ b/Exec/RegTests/EB_DecayingODEQTY/run_convergence_test.py @@ -0,0 +1,128 @@ +import os +import pandas as pd +import numpy as np +import matplotlib.pyplot as plt + +def generate_input_file(final_time, dt, test, temporal_int, template_filename='input.convergence_0'): + new_input_filename = f'input.convergence_{test}' + + # Read the template input file + with open(template_filename, 'r') as file: + lines = file.readlines() + + # Modify amr.fixed_dt + for i, line in enumerate(lines): + line_start = line.strip().split() + if len(line_start) > 0: + line_start = line_start[0] + if line_start == 'amr.fixed_dt': + lines[i] = f'amr.fixed_dt = {dt}\n' + elif line_start == 'amr.plot_file': + lines[i] = f'amr.plot_file = "plt_conv{test}_"\n' + elif line_start == 'peleLM.temporal_int': + lines[i] = f'peleLM.temporal_int = {temporal_int}\n' + elif line_start == 'amr.stop_time': + lines[i] = f'amr.stop_time = {final_time}\n' + # Write the new input file + with open(new_input_filename, 'w') as file: + file.writelines(lines) + + return new_input_filename + +def run_PeleLMeX(input_file,num_proc): + if (num_proc > 1): + os.system(f'mpirun -np {num_proc} ./PeleLMeX2d.gnu.MPI.ex {input_file}') + else: + os.system(f'./PeleLMeX2d.gnu.MPI.ex {input_file}') + +def soln(t,IC): + k = -100.0 + return IC*np.exp(k*t) + +def dataAnalysis(vars,file_name): + data = pd.read_csv(file_name) + data_vars = data[vars].values + time = data.time.values + IC = data_vars[0,:] + + approx = data_vars[-1,:] + exact = soln(time[-1],IC) + abs_error = abs(exact - approx) + + return abs_error.reshape(1,len(vars)) + + +def removeAllPrevData(): + os.system('rm -r plt* temporals/*') + + +# Number of runs for convergence test +startWithNewData = 0 +nRuns = 4 +num_proc = 6 +dt_vec = [1e-3, 1e-4, 1e-5, 1e-6] +final_time = 0.05 +num_temporal_write = 60 +temporals_directory = 'temporals_simpleDecay_SDC0' +temporals_filename = temporals_directory+'/tempExtremas' + +NUM_ODE = 3 +vars = [f"max_MY_ODE_{i}" for i in range(0,NUM_ODE)] +errors = np.zeros((nRuns,NUM_ODE)) + +print("\n==============================================================================") +if startWithNewData: + removeAllPrevData() + +for n in range(0, nRuns): + dt = dt_vec[n] + Nt = round(final_time / dt) + temporal_int = round(Nt / num_temporal_write) + print(f"Running test #{n} with dt = {dt}, Nt = {Nt}, temporal_int = {temporal_int}") + + # Generate a new input file for the given dt + input_file = generate_input_file(final_time, dt, n, temporal_int) + + # Run the convergence test in PeleLMeX + if startWithNewData: + run_PeleLMeX(input_file,num_proc) + + # Rename the tempExtremas file and remove unused tempState file + new_temporals_filename = temporals_filename + f'_{n}' + os.system(f'mv {temporals_filename} {new_temporals_filename}') + os.system('rm temporals/tempState') + + print(f"Finished running test #{n}\n") + + # Read in extrema data from temporals + errors[n,:] = dataAnalysis(vars,new_temporals_filename) + +print("\n==============================================================================\n") + +plt.figure() +for j in range(NUM_ODE): + + + # Fit a line to the log-log data + log_dt_vec = np.log(dt_vec[:nRuns]) + log_errors = np.log(errors[:,j]) + slope, intercept = np.polyfit(log_dt_vec, log_errors, 1) + + # Plot scatter plot of errors and the line of best fit + best_fit_line = np.exp(intercept) * np.array(dt_vec[:nRuns])**slope + plt.scatter(dt_vec[:nRuns], errors[:,j], label=f'ODE {j} (rate: {slope:.2f})', marker='o') + plt.plot(dt_vec[:nRuns], best_fit_line, '--') + + # Reset the color cycle + #plt.gca().set_prop_cycle(None) + +# Set x and y scales to log-log +plt.xscale("log") +plt.yscale("log") + +# Add labels and legend +plt.xlabel('$\Delta t$') +plt.ylabel('Error') +plt.legend() + +plt.show() \ No newline at end of file diff --git a/Source/PeleLMeX.H b/Source/PeleLMeX.H index 880883273..b0af170fe 100644 --- a/Source/PeleLMeX.H +++ b/Source/PeleLMeX.H @@ -1866,6 +1866,7 @@ public: #if NUM_ODE > 0 amrex::Vector m_ode_names; void predictODEQty(); + void advanceODEQty(); #endif // Aux data size diff --git a/Source/PeleLMeX_ODEQty.cpp b/Source/PeleLMeX_ODEQty.cpp index b4759d39c..99eb9e1a8 100644 --- a/Source/PeleLMeX_ODEQty.cpp +++ b/Source/PeleLMeX_ODEQty.cpp @@ -23,4 +23,9 @@ void PeleLM::predictODEQty() }); Gpu::streamSynchronize(); } +} + +void PeleLM::advanceODEQty() +{ + // To do: If m_ext_sources_SDC = true, update with MISDC } \ No newline at end of file diff --git a/Source/PeleLMeX_ProblemSpecificFunctions.cpp b/Source/PeleLMeX_ProblemSpecificFunctions.cpp index edd61cfb8..35abc40e0 100644 --- a/Source/PeleLMeX_ProblemSpecificFunctions.cpp +++ b/Source/PeleLMeX_ProblemSpecificFunctions.cpp @@ -6,9 +6,9 @@ using namespace amrex; /* Problem specific functions: -- This file must be copied in its entirety locally +- This file must be copied locally to the case directory - Add the following to GNUmakefile: CEXE_sources += PeleLMeX_ProblemSpecificFunctions.cpp -- Modify functions as needed +- Modify as needed */ void set_ode_names(Vector& a_ode_names) From d433041904da9d7fa06b6efef16287a6241f6cb1 Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Tue, 15 Oct 2024 11:54:02 -0600 Subject: [PATCH 18/46] Updated the test case to include 3 odes with increasing stiffness. --- .../PeleLMeX_ProblemSpecificFunctions.cpp | 17 +++--- Exec/RegTests/EB_DecayingODEQTY/input.2d | 9 ++-- .../EB_DecayingODEQTY/input.convergence_0 | 19 ++++--- .../EB_DecayingODEQTY/run_convergence_test.py | 53 ++++++++++--------- Source/PeleLMeX_Advance.cpp | 13 ----- 5 files changed, 51 insertions(+), 60 deletions(-) diff --git a/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.cpp b/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.cpp index 02a7f8fad..065e1523e 100644 --- a/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.cpp +++ b/Exec/RegTests/EB_DecayingODEQTY/PeleLMeX_ProblemSpecificFunctions.cpp @@ -29,26 +29,31 @@ void problem_modify_ext_sources( MultiArray4 const& state_old_arr, MultiArray4 const& /*state_new_arr*/, Vector>& a_extSource, - const GeometryData& /*geomdata*/, + const GeometryData& geomdata, ProbParm const& prob_parm) { /* Notes: 1) a_extSource contains sources from velocity forcing coming in. This function should add to rather than overwrite a_extSource. - 2) Requires peleLM.user_defined_ext_sources = true in input file - - TODO: - 1) Add and test capabilities with species, velocity and enthalpy. + 2) Requires "peleLM.user_defined_ext_sources = true" in input file */ auto ext_source_arr = a_extSource[lev]->arrays(); + + const amrex::Real* prob_lo = geomdata.ProbLo(); + const amrex::Real* prob_hi = geomdata.ProbHi(); + const amrex::Real* dx = geomdata.CellSize(); + const amrex::Real Lx = prob_hi[0] - prob_lo[0]; + + amrex::ParallelFor( *a_extSource[lev], [=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept { for (int n = 0; n < NUM_ODE; n++){ amrex::Real B_n = state_old_arr[box_no](i, j, k, FIRSTODE + n); - ext_source_arr[box_no](i, j, k, FIRSTODE + n) += prob_parm.ode_srcstrength*B_n; + amrex::Real src = prob_parm.ode_srcstrength * pow(10.0,n+1) * B_n; + ext_source_arr[box_no](i, j, k, FIRSTODE + n) += src; } }); amrex::Gpu::streamSynchronize(); diff --git a/Exec/RegTests/EB_DecayingODEQTY/input.2d b/Exec/RegTests/EB_DecayingODEQTY/input.2d index ebf728a61..37a46d565 100644 --- a/Exec/RegTests/EB_DecayingODEQTY/input.2d +++ b/Exec/RegTests/EB_DecayingODEQTY/input.2d @@ -30,7 +30,7 @@ prob.ode_IC = 10.0 prob.ode_xy_lo = -0.01 prob.ode_length = 0.01 prob.ode_height = 0.02 -prob.ode_srcstrength = -100.0 +prob.ode_srcstrength = -10.0 #-----------------INPUTS TO CONSTANT TRANSPORT ------------------ transport.units = MKS @@ -41,17 +41,14 @@ transport.const_diffusivity = 0.0 #---------------------------TIME STEPPING------------------------ #amr.max_step = 1 -amr.stop_time = 0.06 -#amr.fixed_dt = 0.0001 +amr.stop_time = 0.02 amr.dt_shrink = 0.1 amr.dt_change_max = 1.1 -amr.cfl = 0.9 +amr.fixed_dt = 0.0001 #-------------------------PELE CONTROLS---------------------- peleLM.v = 1 peleLM.user_defined_ext_sources = 1 -#peleLM.num_init_iter = 0 -#peleLM.do_init_proj = 0 #---------------------------IO CONTROL--------------------------- #amr.restart = chk_07136 diff --git a/Exec/RegTests/EB_DecayingODEQTY/input.convergence_0 b/Exec/RegTests/EB_DecayingODEQTY/input.convergence_0 index d3d633b5c..595a6bdd9 100644 --- a/Exec/RegTests/EB_DecayingODEQTY/input.convergence_0 +++ b/Exec/RegTests/EB_DecayingODEQTY/input.convergence_0 @@ -11,26 +11,26 @@ peleLM.lo_bc = Inflow Interior peleLM.hi_bc = Outflow Interior #-------------------------AMR CONTROL---------------------------- -amr.n_cell = 96 32 # Level 0 number of cells in each direction +amr.n_cell = 48 16 # Level 0 number of cells in each direction amr.v = 1 # AMR verbose amr.max_level = 0 # maximum level number allowed amr.ref_ratio = 2 2 2 2 # refinement ratio amr.regrid_int = 5 # how often to regrid amr.n_error_buf = 2 2 2 2 # number of buffer cells in error est amr.grid_eff = 0.7 # what constitutes an efficient grid -amr.blocking_factor = 16 # block factor in grid generation (min box size) +amr.blocking_factor = 8 # block factor in grid generation (min box size) amr.max_grid_size = 64 # max box size #--------------------------- Problem ---------------------------- prob.T_mean = 300.0 prob.P_mean = 101325.0 -prob.meanFlowMag = 4.255 +prob.meanFlowMag = 1.0 prob.meanFlowDir = 1 prob.ode_IC = 10.0 prob.ode_xy_lo = -0.01 prob.ode_length = 0.01 prob.ode_height = 0.02 -prob.ode_srcstrength = -100.0 +prob.ode_srcstrength = -10.0 #-----------------INPUTS TO CONSTANT TRANSPORT ------------------ transport.units = MKS @@ -40,14 +40,12 @@ transport.const_conductivity = 0.0 transport.const_diffusivity = 0.0 #---------------------------TIME STEPPING------------------------ -amr.stop_time = 0.05 -amr.fixed_dt = 0.001 +amr.stop_time = 0.001 +amr.fixed_dt = 0.0001 #-------------------------PELE CONTROLS---------------------- peleLM.v = 1 peleLM.user_defined_ext_sources = 1 -#peleLM.num_init_iter = 0 -#peleLM.do_init_proj = 0 #---------------------------IO CONTROL--------------------------- #amr.restart = chk_07136 @@ -59,11 +57,12 @@ amr.derive_plot_vars = avg_pressure mag_vort #---------------------- Temporal CONTROL ------------------------- peleLM.do_temporals = 1 # Turn temporals ON/OFF -peleLM.temporal_int = 1 +peleLM.temporal_int = 5 peleLM.do_extremas = 1 # Compute state extremas #------------------------- EB SETUP ----------------------------- -eb2.geom_type = sphere +eb2.geom_type = "all_regular" +#eb2.geom_type = sphere eb2.sphere_radius = 0.005 eb2.sphere_center = 0.03 0.0 eb2.sphere_has_fluid_inside = 0 diff --git a/Exec/RegTests/EB_DecayingODEQTY/run_convergence_test.py b/Exec/RegTests/EB_DecayingODEQTY/run_convergence_test.py index adc93560e..5c7f2e65d 100644 --- a/Exec/RegTests/EB_DecayingODEQTY/run_convergence_test.py +++ b/Exec/RegTests/EB_DecayingODEQTY/run_convergence_test.py @@ -35,36 +35,43 @@ def run_PeleLMeX(input_file,num_proc): else: os.system(f'./PeleLMeX2d.gnu.MPI.ex {input_file}') -def soln(t,IC): - k = -100.0 - return IC*np.exp(k*t) - -def dataAnalysis(vars,file_name): +def dataAnalysis(vars, file_name): data = pd.read_csv(file_name) data_vars = data[vars].values time = data.time.values - IC = data_vars[0,:] + IC = data_vars[0, :] - approx = data_vars[-1,:] - exact = soln(time[-1],IC) - abs_error = abs(exact - approx) - - return abs_error.reshape(1,len(vars)) + # Numerical solution at final time + approx = data_vars[-1, :] + + # Exact solution + n = np.arange(len(IC)) + k = -10.0 * 10**(n + 1) + exact = IC * np.exp(k * time[-1]) + + # Absolute error + abs_error = np.abs(exact - approx) + + return abs_error.reshape(1, len(vars)) def removeAllPrevData(): os.system('rm -r plt* temporals/*') +# ----------------------------------------------------------------------------- +# Parameters +# ----------------------------------------------------------------------------- + +startWithNewData = 1 # erase all plt files and temporals if true +nRuns = 4 # number of tests +num_proc = 2 # number of processors to run tests +dt_vec = [1e-4, 1e-5, 1e-6, 1e-7] # dt for each test +final_time = 0.001 # final time for simulation -# Number of runs for convergence test -startWithNewData = 0 -nRuns = 4 -num_proc = 6 -dt_vec = [1e-3, 1e-4, 1e-5, 1e-6] -final_time = 0.05 -num_temporal_write = 60 -temporals_directory = 'temporals_simpleDecay_SDC0' +# Parameters for temporal files +temporals_directory = 'temporals' temporals_filename = temporals_directory+'/tempExtremas' +num_temporal_write = 2 NUM_ODE = 3 vars = [f"max_MY_ODE_{i}" for i in range(0,NUM_ODE)] @@ -101,9 +108,8 @@ def removeAllPrevData(): plt.figure() for j in range(NUM_ODE): - - # Fit a line to the log-log data + # Fit a log-log line to the data log_dt_vec = np.log(dt_vec[:nRuns]) log_errors = np.log(errors[:,j]) slope, intercept = np.polyfit(log_dt_vec, log_errors, 1) @@ -113,16 +119,13 @@ def removeAllPrevData(): plt.scatter(dt_vec[:nRuns], errors[:,j], label=f'ODE {j} (rate: {slope:.2f})', marker='o') plt.plot(dt_vec[:nRuns], best_fit_line, '--') - # Reset the color cycle - #plt.gca().set_prop_cycle(None) - # Set x and y scales to log-log plt.xscale("log") plt.yscale("log") # Add labels and legend plt.xlabel('$\Delta t$') -plt.ylabel('Error') +plt.ylabel('$\ell_\infty(approx - exact)$') plt.legend() plt.show() \ No newline at end of file diff --git a/Source/PeleLMeX_Advance.cpp b/Source/PeleLMeX_Advance.cpp index 749e2c2de..ec43fb700 100644 --- a/Source/PeleLMeX_Advance.cpp +++ b/Source/PeleLMeX_Advance.cpp @@ -140,19 +140,6 @@ PeleLM::Advance(int is_initIter) getLevelDataPtr(lev, AmrOldTime)->state.const_arrays(), getLevelDataPtr(lev, AmrNewTime)->state.const_arrays(), m_extSource, geom[lev].data(), *prob_parm_d); - - // Debugging: Test the call for modify ext_sources - auto ext_src = m_extSource[lev]->arrays(); - amrex::ParallelFor( - *m_extSource[lev], - [=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept { - for(int n = 0; n < NUM_ODE; n++){ - //if(ext_src[box_no](i, j, k, FIRSTODE + n) < 0){ - // Print() << "ext_src[lev = "<state.const_arrays(), - getLevelDataPtr(lev, AmrNewTime)->state.const_arrays(), - m_extSource, geom[lev].data(), *prob_parm_d); - } - } - + // External sources (soot, radiation, user defined, etc.) + getExternalSources(AmrOldTime,AmrNewTime); if (m_incompressible == 0) { floorSpecies(AmrOldTime); diff --git a/Source/PeleLMeX_Forces.cpp b/Source/PeleLMeX_Forces.cpp index 028427ae0..38f024bb0 100644 --- a/Source/PeleLMeX_Forces.cpp +++ b/Source/PeleLMeX_Forces.cpp @@ -1,5 +1,6 @@ #include #include +#include using namespace amrex; @@ -234,4 +235,47 @@ PeleLM::addSpark(const TimeStamp& a_timestamp) Gpu::streamSynchronize(); } } +} + +// Calculate additional external sources (soot, radiation, user defined, etc.) +void PeleLM::getExternalSources( + const PeleLM::TimeStamp& a_timestamp_old, + const PeleLM::TimeStamp& a_timestamp_new) +{ + if (m_n_sparks > 0) { + addSpark(a_timestamp_old); + } + +#ifdef PELE_USE_SPRAY + if (is_initIter == 0) { + SprayMKD(m_cur_time, m_dt); + } +#endif +#ifdef PELE_USE_SOOT + if (do_soot_solve) { + computeSootSource(a_timestamp_old, m_dt); + } +#endif +#ifdef PELE_USE_RADIATION + if (do_rad_solve) { + BL_PROFILE_VAR("PeleLM::advance::rad", PLM_RAD); + computeRadSource(a_timestamp_old); + BL_PROFILE_VAR_STOP(PLM_RAD); + } +#endif + +// User defined external sources + if(m_user_defined_ext_sources){ + for (int lev = 0; lev <= finest_level; lev++) { + problem_modify_ext_sources( + getTime(lev, a_timestamp_new), + m_dt, + lev, + getLevelDataPtr(lev, a_timestamp_old)->state.const_arrays(), + getLevelDataPtr(lev, a_timestamp_new)->state.const_arrays(), + m_extSource, + geom[lev].data(), + *prob_parm_d); + } + } } \ No newline at end of file diff --git a/Source/PeleLMeX_ODEQty.cpp b/Source/PeleLMeX_ODEQty.cpp index 99eb9e1a8..b4759d39c 100644 --- a/Source/PeleLMeX_ODEQty.cpp +++ b/Source/PeleLMeX_ODEQty.cpp @@ -23,9 +23,4 @@ void PeleLM::predictODEQty() }); Gpu::streamSynchronize(); } -} - -void PeleLM::advanceODEQty() -{ - // To do: If m_ext_sources_SDC = true, update with MISDC } \ No newline at end of file diff --git a/Source/PeleLMeX_ProblemSpecificFunctions.cpp b/Source/PeleLMeX_ProblemSpecificFunctions.cpp index 35abc40e0..3e5d8fcf9 100644 --- a/Source/PeleLMeX_ProblemSpecificFunctions.cpp +++ b/Source/PeleLMeX_ProblemSpecificFunctions.cpp @@ -7,7 +7,8 @@ using namespace amrex; /* Problem specific functions: - This file must be copied locally to the case directory -- Add the following to GNUmakefile: CEXE_sources += PeleLMeX_ProblemSpecificFunctions.cpp +- Add the following to GNUmakefile: + CEXE_sources += PeleLMeX_ProblemSpecificFunctions.cpp - Modify as needed */ @@ -25,8 +26,8 @@ static void problem_modify_ext_sources( Real /*time*/, Real /*dt*/, int /*lev*/, - MultiArray4 const& /*state_old*/, - MultiArray4 const& /*state_new*/, + MultiArray4 const& /*state_old_arr*/, + MultiArray4 const& /*state_new_arr*/, Vector>& /*a_extSource*/, const GeometryData& /*geomdata*/, ProbParm const& /*prob_parm*/) @@ -36,20 +37,18 @@ static void problem_modify_ext_sources( 1) a_extSource contains sources from velocity forcing coming in. This function should add to rather than overwrite a_extSource. 2) Requires peleLM.user_defined_ext_sources = true in input file - - TODO: - 1) Add and test capabilities with species, velocity and enthalpy. - // Example: adding a exponential decay source term to an ode_qty + // Example: Exponential decay ode quantity auto ext_source_arr = a_extSource[lev]->arrays(); - amrex::ParallelFor( + ParallelFor( *a_extSource[lev], [=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept { for (int n = 0; n < NUM_ODE; n++){ - amrex::Real B_n = state_old[box_no](i, j, k, FIRSTODE + n); - ext_source_arr[box_no](i, j, k, FIRSTODE + n) += -10.*B_n; + Real B_n = state_old_arr[box_no](i, j, k, FIRSTODE + n); + Real src_strength = -1.0 * pow(10.0,n+1); + ext_source_arr[box_no](i, j, k, FIRSTODE + n) += src_strength * B_n; } }); - amrex::Gpu::streamSynchronize(); + Gpu::streamSynchronize(); */ } \ No newline at end of file From 8f447ff50c1e29fbb383dd73adde1183da5e0a8e Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Tue, 15 Oct 2024 13:44:05 -0600 Subject: [PATCH 21/46] Update README for case. --- Exec/RegTests/EB_ODEQty/README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Exec/RegTests/EB_ODEQty/README.md b/Exec/RegTests/EB_ODEQty/README.md index 19fb1a405..a174227a6 100644 --- a/Exec/RegTests/EB_ODEQty/README.md +++ b/Exec/RegTests/EB_ODEQty/README.md @@ -1,2 +1,2 @@ -## EB\_DecayingODEQty -A 2D (or 3D) inflow/outflow setup with an optional EB cylinder in the middle of the flow. Testing the "ODE" variables. +## EB\_ODEQty +A 2D inflow/outflow setup with an optional EB cylinder in the middle of the flow. Demonstrates how to use ProblemSpecificFunctions and the ODE quantities. The ODE quantities experience simple exponential decay that gets stiffer for each quantity. Specifically, $\frac{\partial B_k}{\partial t} = k \cdot 10^{k+1} B_k$, for $k = 0, 1, ...,$ NUM_ODE. From e24675de0f7329d2f12cc5809791c4f241120ddf Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Tue, 15 Oct 2024 13:45:37 -0600 Subject: [PATCH 22/46] Update README for case. --- Exec/RegTests/EB_ODEQty/README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Exec/RegTests/EB_ODEQty/README.md b/Exec/RegTests/EB_ODEQty/README.md index a174227a6..a141c89cd 100644 --- a/Exec/RegTests/EB_ODEQty/README.md +++ b/Exec/RegTests/EB_ODEQty/README.md @@ -1,2 +1,2 @@ ## EB\_ODEQty -A 2D inflow/outflow setup with an optional EB cylinder in the middle of the flow. Demonstrates how to use ProblemSpecificFunctions and the ODE quantities. The ODE quantities experience simple exponential decay that gets stiffer for each quantity. Specifically, $\frac{\partial B_k}{\partial t} = k \cdot 10^{k+1} B_k$, for $k = 0, 1, ...,$ NUM_ODE. +A 2D inflow/outflow setup with an optional EB cylinder in the middle of the flow. Demonstrates how to use ProblemSpecificFunctions and the ODE quantities. The ODE quantities experience simple exponential decay that gets stiffer for each quantity. Specifically, $\frac{\partial B_k}{\partial t} = \gamma \cdot 10^{k+1} B_k$, for $k = 0, 1, ...,$ NUM_ODE and $\gamma < 0$. From bd21d5a58ce53f6b0c3634827415e79932939e57 Mon Sep 17 00:00:00 2001 From: Dave Montgomery Date: Tue, 15 Oct 2024 13:47:34 -0600 Subject: [PATCH 23/46] Update README.md --- Exec/RegTests/EB_ODEQty/README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Exec/RegTests/EB_ODEQty/README.md b/Exec/RegTests/EB_ODEQty/README.md index a141c89cd..c04720a5d 100644 --- a/Exec/RegTests/EB_ODEQty/README.md +++ b/Exec/RegTests/EB_ODEQty/README.md @@ -1,2 +1,2 @@ ## EB\_ODEQty -A 2D inflow/outflow setup with an optional EB cylinder in the middle of the flow. Demonstrates how to use ProblemSpecificFunctions and the ODE quantities. The ODE quantities experience simple exponential decay that gets stiffer for each quantity. Specifically, $\frac{\partial B_k}{\partial t} = \gamma \cdot 10^{k+1} B_k$, for $k = 0, 1, ...,$ NUM_ODE and $\gamma < 0$. +A 2D inflow/outflow setup with an optional EB cylinder in the middle of the flow. Demonstrates how to use ProblemSpecificFunctions and the ODE quantities. The ODE quantities experience simple exponential decay that gets stiffer for each quantity. Specifically, $\frac{\partial B_k}{\partial t} = \gamma \cdot 10^{k+1} B_k$, for $k = 0, 1,\dots,$ NUM_ODE and $\gamma < 0$. From 5700778ed14087257cc2401f973bcd4cccb61e3a Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Tue, 15 Oct 2024 13:53:22 -0600 Subject: [PATCH 24/46] Removed NUM_AUX for now. --- Exec/Make.PeleLMeX | 4 ---- Source/PeleLMeX_Index.H | 3 +-- 2 files changed, 1 insertion(+), 6 deletions(-) diff --git a/Exec/Make.PeleLMeX b/Exec/Make.PeleLMeX index 9f2cc32b8..4ef1c7734 100644 --- a/Exec/Make.PeleLMeX +++ b/Exec/Make.PeleLMeX @@ -61,10 +61,6 @@ ifeq ($(shell test 0$(NUM_ODE) -gt 0; echo $$?), 0) DEFINES += -DNUM_ODE=$(NUM_ODE) endif -ifeq ($(shell test 0$(NUM_AUX) -gt 0; echo $$?), 0) - DEFINES += -DNUM_AUX=$(NUM_AUX) -endif - Bpack += $(foreach dir, $(LMdirs), $(PELE_HOME)/$(dir)/Make.package) Blocs += $(foreach dir, $(LMdirs), $(PELE_HOME)/$(dir)) diff --git a/Source/PeleLMeX_Index.H b/Source/PeleLMeX_Index.H index 2d9ce9b95..09099b410 100644 --- a/Source/PeleLMeX_Index.H +++ b/Source/PeleLMeX_Index.H @@ -27,8 +27,7 @@ #endif #define FIRSTODE (RHORT + NUMEFIELD + NUMSOOTVAR + 1) #define FIRSTAUX (RHORT + NUMEFIELD + NUMSOOTVAR + NUM_ODE + 1) -#define NUM_AUX 0 -#define NVAR (RHORT + NUMEFIELD + NUMSOOTVAR + NUM_ODE + NUM_AUX + 1) +#define NVAR (RHORT + NUMEFIELD + NUMSOOTVAR + NUM_ODE + 1) // TODO: handle NAUX #endif From afb0aa7333a58ef9cec1d058b35ff12ca9c22ce2 Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Tue, 15 Oct 2024 14:01:48 -0600 Subject: [PATCH 25/46] Cleaning up --- .../EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp | 1 - Exec/RegTests/EB_ODEQty/pelelmex_prob.H | 1 - Exec/RegTests/EB_ODEQty/pelelmex_prob_parm.H | 2 +- Source/PeleLMeX_Plot.cpp | 7 ------- 4 files changed, 1 insertion(+), 10 deletions(-) diff --git a/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp b/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp index 482e493e0..0b61e698c 100644 --- a/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp +++ b/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp @@ -12,7 +12,6 @@ Problem specific functions: - Modify as needed */ - void set_ode_names(Vector& a_ode_names) { #if NUM_ODE > 0 diff --git a/Exec/RegTests/EB_ODEQty/pelelmex_prob.H b/Exec/RegTests/EB_ODEQty/pelelmex_prob.H index ab7dd004c..50d5182ad 100644 --- a/Exec/RegTests/EB_ODEQty/pelelmex_prob.H +++ b/Exec/RegTests/EB_ODEQty/pelelmex_prob.H @@ -136,7 +136,6 @@ bcnormal( ProbParm const& prob_parm, pele::physics::PMF::PmfData::DataContainer const* /*pmf_data*/) { - // amrex::GpuArray pmf_vals = {0.0}; amrex::Real massfrac[NUM_SPECIES] = {0.0}; auto eos = pele::physics::PhysicsType::eos(); diff --git a/Exec/RegTests/EB_ODEQty/pelelmex_prob_parm.H b/Exec/RegTests/EB_ODEQty/pelelmex_prob_parm.H index b74ae6ce2..ba51a72e5 100644 --- a/Exec/RegTests/EB_ODEQty/pelelmex_prob_parm.H +++ b/Exec/RegTests/EB_ODEQty/pelelmex_prob_parm.H @@ -16,7 +16,7 @@ struct ProbParm : amrex::Gpu::Managed amrex::Real ode_xy_lo = -0.01; amrex::Real ode_length = 0.01; amrex::Real ode_height = 0.020; - amrex::Real ode_srcstrength = -100.0; + amrex::Real ode_srcstrength = -10.0; int meanFlowDir = 1; }; #endif diff --git a/Source/PeleLMeX_Plot.cpp b/Source/PeleLMeX_Plot.cpp index cfc4378b2..3d8991acf 100644 --- a/Source/PeleLMeX_Plot.cpp +++ b/Source/PeleLMeX_Plot.cpp @@ -288,13 +288,6 @@ PeleLM::WritePlotFile() #endif Print() << PrettyLine; -// Debugging: -//Print() << "Plotting debugging" << std::endl; -//Print() << "ncomp = " << ncomp << std::endl; -//Print() << "plt_VarsName.size() = " << plt_VarsName.size() << std::endl; -//for(int n = 0; n < plt_VarsName.size(); n++){ -// Print() << "plt_VarsName[" << n << "] = " << plt_VarsName[n] << std::endl; -//} //---------------------------------------------------------------- // Fill the plot MultiFabs From 91d4b666b6c12befee8ee42a66ba078fc71c3320 Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Tue, 15 Oct 2024 15:03:52 -0600 Subject: [PATCH 26/46] Documentation for ode quantities. --- Docs/sphinx/manual/LMeXControls.rst | 2 ++ Docs/sphinx/manual/Model.rst | 6 ++++++ 2 files changed, 8 insertions(+) diff --git a/Docs/sphinx/manual/LMeXControls.rst b/Docs/sphinx/manual/LMeXControls.rst index 77cfeb965..bb3e757e6 100644 --- a/Docs/sphinx/manual/LMeXControls.rst +++ b/Docs/sphinx/manual/LMeXControls.rst @@ -269,6 +269,8 @@ PeleLMeX algorithm peleLM.spark1.radius = 1e-3 # [OPT] Radius of the spark [m] peleLM.spark1.duration = 1e-3 # [OPT] Duration of the spark [s] peleLM.spark1.time = 1e-2 # [OPT] Time when spark starts [s] + + peleLM.user_defined_ext_sources = 0 # [OPT, DEF=0] Enable user defined source terms. Requires local ProblemSpecificFunctions.cpp. Transport coefficients and LES diff --git a/Docs/sphinx/manual/Model.rst b/Docs/sphinx/manual/Model.rst index 62cd5b435..b83b55ce3 100644 --- a/Docs/sphinx/manual/Model.rst +++ b/Docs/sphinx/manual/Model.rst @@ -118,6 +118,12 @@ For the standard ideal gas EOS, the divergence constraint on velocity becomes: \nabla \cdot \boldsymbol{u} &= \frac{1}{\rho c_p T} \left(\nabla \cdot \lambda\nabla T - \sum_m \boldsymbol{\Gamma_m} \cdot \nabla h_m \right) \\ &- \frac{1}{\rho} \sum_m \frac{W}{W_m} \nabla \cdot \boldsymbol{\Gamma_m} + \frac{1}{\rho}\sum_m \left(\frac{W}{W_m} - \frac{h_m}{c_p T} \right) \dot \omega \equiv S . +In addition to the flow equations, `PeleLMeX` can also solve for a set of quantities that are neither advected nor diffused, satifying: + +.. math:: + + \frac{\partial B_k}{\partial t} = S_{\text{ext},B_k}. + Confined domain ambient pressure ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ From c1a94cc9a807fb0924863ee23082867b5702c939 Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Tue, 15 Oct 2024 15:06:26 -0600 Subject: [PATCH 27/46] Clang-format --- .../PeleLMeX_ProblemSpecificFunctions.cpp | 45 ++++++++++--------- Exec/RegTests/EB_ODEQty/pelelmex_prob.H | 2 +- Source/PeleLMeX.H | 11 +++-- Source/PeleLMeX_Advance.cpp | 11 +++-- Source/PeleLMeX_Forces.cpp | 17 +++---- Source/PeleLMeX_Index.H | 2 +- Source/PeleLMeX_ODEQty.cpp | 19 ++++---- Source/PeleLMeX_Plot.cpp | 17 ++++--- Source/PeleLMeX_ProblemSpecificFunctions.H | 16 +++---- Source/PeleLMeX_ProblemSpecificFunctions.cpp | 38 ++++++++-------- Source/PeleLMeX_Utils.cpp | 10 ++--- 11 files changed, 92 insertions(+), 96 deletions(-) diff --git a/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp b/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp index 0b61e698c..38139cef3 100644 --- a/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp +++ b/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp @@ -4,36 +4,38 @@ using namespace amrex; -/* +/* Problem specific functions: - This file must be copied locally to the case directory -- Add the following to GNUmakefile: +- Add the following to GNUmakefile: CEXE_sources += PeleLMeX_ProblemSpecificFunctions.cpp - Modify as needed */ -void set_ode_names(Vector& a_ode_names) +void +set_ode_names(Vector& a_ode_names) { #if NUM_ODE > 0 - a_ode_names.resize(NUM_ODE); - for (int n = 0; n < NUM_ODE; n++) { - a_ode_names[n] = "MY_ODE_" + std::to_string(n); - } + a_ode_names.resize(NUM_ODE); + for (int n = 0; n < NUM_ODE; n++) { + a_ode_names[n] = "MY_ODE_" + std::to_string(n); + } #endif } -void problem_modify_ext_sources( - Real /*time*/, - Real /*dt*/, - int lev, - MultiArray4 const& state_old_arr, - MultiArray4 const& /*state_new_arr*/, - Vector>& a_extSource, - const GeometryData& geomdata, - ProbParm const& prob_parm) +void +problem_modify_ext_sources( + Real /*time*/, + Real /*dt*/, + int lev, + MultiArray4 const& state_old_arr, + MultiArray4 const& /*state_new_arr*/, + Vector>& a_extSource, + const GeometryData& geomdata, + ProbParm const& prob_parm) { - /* - Notes: + /* + Notes: 1) a_extSource contains sources from velocity forcing coming in. This function should add to rather than overwrite a_extSource. 2) Requires "peleLM.user_defined_ext_sources = true" in input file @@ -42,14 +44,13 @@ void problem_modify_ext_sources( auto ext_source_arr = a_extSource[lev]->arrays(); ParallelFor( - *a_extSource[lev], + *a_extSource[lev], [=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept { - for (int n = 0; n < NUM_ODE; n++){ + for (int n = 0; n < NUM_ODE; n++) { Real B_n = state_old_arr[box_no](i, j, k, FIRSTODE + n); - Real src = prob_parm.ode_srcstrength * pow(10.0,n+1) * B_n; + Real src = prob_parm.ode_srcstrength * pow(10.0, n + 1) * B_n; ext_source_arr[box_no](i, j, k, FIRSTODE + n) += src; } }); Gpu::streamSynchronize(); - } \ No newline at end of file diff --git a/Exec/RegTests/EB_ODEQty/pelelmex_prob.H b/Exec/RegTests/EB_ODEQty/pelelmex_prob.H index 50d5182ad..79de0f3fd 100644 --- a/Exec/RegTests/EB_ODEQty/pelelmex_prob.H +++ b/Exec/RegTests/EB_ODEQty/pelelmex_prob.H @@ -118,7 +118,7 @@ pelelmex_initdata( for (int n = 0; n < NUM_ODE; n++) { state(i, j, k, FIRSTODE + n) = ode_qty[n]; } - + #endif } diff --git a/Source/PeleLMeX.H b/Source/PeleLMeX.H index ce82eb1d9..977f3864b 100644 --- a/Source/PeleLMeX.H +++ b/Source/PeleLMeX.H @@ -882,7 +882,6 @@ public: const PeleLM::TimeStamp& a_timestamp_old, const PeleLM::TimeStamp& a_timestamp_new); - //----------------------------------------------------------------------------- // EOS @@ -1865,11 +1864,11 @@ public: amrex::Real m_divu_dtFactor = 0.5; amrex::Real m_divu_rhoMin = 0.1; - // ODE variable names - #if NUM_ODE > 0 - amrex::Vector m_ode_names; - void predictODEQty(); - #endif +// ODE variable names +#if NUM_ODE > 0 + amrex::Vector m_ode_names; + void predictODEQty(); +#endif // Aux data size int m_nAux = 0; diff --git a/Source/PeleLMeX_Advance.cpp b/Source/PeleLMeX_Advance.cpp index b543921d8..d29ca27cd 100644 --- a/Source/PeleLMeX_Advance.cpp +++ b/Source/PeleLMeX_Advance.cpp @@ -111,7 +111,7 @@ PeleLM::Advance(int is_initIter) //---------------------------------------------------------------- // External sources (soot, radiation, user defined, etc.) - getExternalSources(AmrOldTime,AmrNewTime); + getExternalSources(AmrOldTime, AmrNewTime); if (m_incompressible == 0) { floorSpecies(AmrOldTime); @@ -140,10 +140,10 @@ PeleLM::Advance(int is_initIter) } #if NUM_ODE > 0 - // Euler step for predicting ode qty at tnp1 - if(m_user_defined_ext_sources){ - predictODEQty(); - } + // Euler step for predicting ode qty at tnp1 + if (m_user_defined_ext_sources) { + predictODEQty(); + } #endif BL_PROFILE_VAR_STOP(PLM_SETUP); //---------------------------------------------------------------- @@ -433,4 +433,3 @@ PeleLM::oneSDC( floorSpecies(AmrNewTime); setThermoPress(AmrNewTime); } - diff --git a/Source/PeleLMeX_Forces.cpp b/Source/PeleLMeX_Forces.cpp index 38f024bb0..f7de439d3 100644 --- a/Source/PeleLMeX_Forces.cpp +++ b/Source/PeleLMeX_Forces.cpp @@ -238,7 +238,8 @@ PeleLM::addSpark(const TimeStamp& a_timestamp) } // Calculate additional external sources (soot, radiation, user defined, etc.) -void PeleLM::getExternalSources( +void +PeleLM::getExternalSources( const PeleLM::TimeStamp& a_timestamp_old, const PeleLM::TimeStamp& a_timestamp_new) { @@ -264,18 +265,14 @@ void PeleLM::getExternalSources( } #endif -// User defined external sources - if(m_user_defined_ext_sources){ + // User defined external sources + if (m_user_defined_ext_sources) { for (int lev = 0; lev <= finest_level; lev++) { problem_modify_ext_sources( - getTime(lev, a_timestamp_new), - m_dt, - lev, + getTime(lev, a_timestamp_new), m_dt, lev, getLevelDataPtr(lev, a_timestamp_old)->state.const_arrays(), - getLevelDataPtr(lev, a_timestamp_new)->state.const_arrays(), - m_extSource, - geom[lev].data(), - *prob_parm_d); + getLevelDataPtr(lev, a_timestamp_new)->state.const_arrays(), + m_extSource, geom[lev].data(), *prob_parm_d); } } } \ No newline at end of file diff --git a/Source/PeleLMeX_Index.H b/Source/PeleLMeX_Index.H index 09099b410..8483b6202 100644 --- a/Source/PeleLMeX_Index.H +++ b/Source/PeleLMeX_Index.H @@ -22,7 +22,7 @@ #else #define NUMSOOTVAR 0 #endif -#ifndef NUM_ODE +#ifndef NUM_ODE #define NUM_ODE 0 #endif #define FIRSTODE (RHORT + NUMEFIELD + NUMSOOTVAR + 1) diff --git a/Source/PeleLMeX_ODEQty.cpp b/Source/PeleLMeX_ODEQty.cpp index b4759d39c..72e846acb 100644 --- a/Source/PeleLMeX_ODEQty.cpp +++ b/Source/PeleLMeX_ODEQty.cpp @@ -1,10 +1,10 @@ #include #include - using namespace amrex; -void PeleLM::predictODEQty() +void +PeleLM::predictODEQty() { // Uses forward Euler to predict values for ODE qty at tnp1 // If m_ext_sources_SDC = false, no SDC corrections used @@ -12,15 +12,14 @@ void PeleLM::predictODEQty() auto const& state_arrs = getLevelDataPtr(lev, AmrNewTime)->state.arrays(); auto const& ext_src_arrs = m_extSource[lev]->arrays(); ParallelFor( - *m_extSource[lev], - [state_arrs, ext_src_arrs, dt = m_dt] - AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept { - for (int n = 0; n < NUM_ODE; n++){ - Real const& B_n = state_arrs[box_no](i, j, k, FIRSTODE + n); - Real const& S_ext_n = ext_src_arrs[box_no](i, j, k, FIRSTODE + n); - state_arrs[box_no](i, j, k, FIRSTODE + n) = B_n + dt * S_ext_n; + *m_extSource[lev], [state_arrs, ext_src_arrs, dt = m_dt] AMREX_GPU_DEVICE( + int box_no, int i, int j, int k) noexcept { + for (int n = 0; n < NUM_ODE; n++) { + Real const& B_n = state_arrs[box_no](i, j, k, FIRSTODE + n); + Real const& S_ext_n = ext_src_arrs[box_no](i, j, k, FIRSTODE + n); + state_arrs[box_no](i, j, k, FIRSTODE + n) = B_n + dt * S_ext_n; } - }); + }); Gpu::streamSynchronize(); } } \ No newline at end of file diff --git a/Source/PeleLMeX_Plot.cpp b/Source/PeleLMeX_Plot.cpp index 3d8991acf..0d559427e 100644 --- a/Source/PeleLMeX_Plot.cpp +++ b/Source/PeleLMeX_Plot.cpp @@ -282,12 +282,12 @@ PeleLM::WritePlotFile() } #if NUM_ODE > 0 - for (int n = 0; n < NUM_ODE; n++) { - plt_VarsName.push_back(m_ode_names[n]); - } + for (int n = 0; n < NUM_ODE; n++) { + plt_VarsName.push_back(m_ode_names[n]); + } #endif -Print() << PrettyLine; + Print() << PrettyLine; //---------------------------------------------------------------- // Fill the plot MultiFabs @@ -410,11 +410,10 @@ Print() << PrettyLine; } #endif #if NUM_ODE > 0 - MultiFab::Copy( - mf_plt[lev], m_leveldata_new[lev]->state, FIRSTODE, cnt, NUM_ODE, - 0); - cnt += NUM_ODE; - Print() << "cnt = " << cnt << std::endl; + MultiFab::Copy( + mf_plt[lev], m_leveldata_new[lev]->state, FIRSTODE, cnt, NUM_ODE, 0); + cnt += NUM_ODE; + Print() << "cnt = " << cnt << std::endl; #endif if (m_do_les && m_plot_les) { diff --git a/Source/PeleLMeX_ProblemSpecificFunctions.H b/Source/PeleLMeX_ProblemSpecificFunctions.H index aeb92356c..697064b3b 100644 --- a/Source/PeleLMeX_ProblemSpecificFunctions.H +++ b/Source/PeleLMeX_ProblemSpecificFunctions.H @@ -8,12 +8,12 @@ class PeleLM; void set_ode_names(Vector& a_ode_names); void problem_modify_ext_sources( - Real time, - Real dt, - int lev, - MultiArray4 const& state_old, - MultiArray4 const& state_new, - Vector>& a_extSource, - const GeometryData& geomdata, - ProbParm const& prob_parm); + Real time, + Real dt, + int lev, + MultiArray4 const& state_old, + MultiArray4 const& state_new, + Vector>& a_extSource, + const GeometryData& geomdata, + ProbParm const& prob_parm); #endif \ No newline at end of file diff --git a/Source/PeleLMeX_ProblemSpecificFunctions.cpp b/Source/PeleLMeX_ProblemSpecificFunctions.cpp index 3e5d8fcf9..b0a5d95ce 100644 --- a/Source/PeleLMeX_ProblemSpecificFunctions.cpp +++ b/Source/PeleLMeX_ProblemSpecificFunctions.cpp @@ -4,36 +4,38 @@ using namespace amrex; -/* +/* Problem specific functions: - This file must be copied locally to the case directory -- Add the following to GNUmakefile: +- Add the following to GNUmakefile: CEXE_sources += PeleLMeX_ProblemSpecificFunctions.cpp - Modify as needed */ -void set_ode_names(Vector& a_ode_names) +void +set_ode_names(Vector& a_ode_names) { #if NUM_ODE > 0 - a_ode_names.resize(NUM_ODE); - for (int n = 0; n < NUM_ODE; n++) { - a_ode_names[n] = "ODE_" + std::to_string(n); - } + a_ode_names.resize(NUM_ODE); + for (int n = 0; n < NUM_ODE; n++) { + a_ode_names[n] = "ODE_" + std::to_string(n); + } #endif } -static void problem_modify_ext_sources( - Real /*time*/, - Real /*dt*/, - int /*lev*/, - MultiArray4 const& /*state_old_arr*/, - MultiArray4 const& /*state_new_arr*/, - Vector>& /*a_extSource*/, - const GeometryData& /*geomdata*/, - ProbParm const& /*prob_parm*/) +static void +problem_modify_ext_sources( + Real /*time*/, + Real /*dt*/, + int /*lev*/, + MultiArray4 const& /*state_old_arr*/, + MultiArray4 const& /*state_new_arr*/, + Vector>& /*a_extSource*/, + const GeometryData& /*geomdata*/, + ProbParm const& /*prob_parm*/) { /* - Notes: + Notes: 1) a_extSource contains sources from velocity forcing coming in. This function should add to rather than overwrite a_extSource. 2) Requires peleLM.user_defined_ext_sources = true in input file @@ -41,7 +43,7 @@ static void problem_modify_ext_sources( // Example: Exponential decay ode quantity auto ext_source_arr = a_extSource[lev]->arrays(); ParallelFor( - *a_extSource[lev], + *a_extSource[lev], [=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept { for (int n = 0; n < NUM_ODE; n++){ Real B_n = state_old_arr[box_no](i, j, k, FIRSTODE + n); diff --git a/Source/PeleLMeX_Utils.cpp b/Source/PeleLMeX_Utils.cpp index 0ec239cb6..707c00492 100644 --- a/Source/PeleLMeX_Utils.cpp +++ b/Source/PeleLMeX_Utils.cpp @@ -1430,11 +1430,11 @@ PeleLM::setTypicalValues(const TimeStamp& a_time, int is_init) #endif #if NUM_ODE > 0 for (int n = 0; n < NUM_ODE; n++) { - Print() << "\t" << m_ode_names[n] - << std::setw( - std::max(0, static_cast(10 - m_ode_names[n].length()))) - << std::left << ":" << typical_values[FIRSTODE + n] <<'\n'; - } + Print() << "\t" << m_ode_names[n] + << std::setw(std::max( + 0, static_cast(10 - m_ode_names[n].length()))) + << std::left << ":" << typical_values[FIRSTODE + n] << '\n'; + } #endif } Print() << PrettyLine; From 9a9de0ad3292efb54c05de422d210f7ad51974d2 Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Tue, 15 Oct 2024 15:29:15 -0600 Subject: [PATCH 28/46] Remove validation tests from case directory. --- Exec/RegTests/EB_ODEQty/input.convergence_0 | 78 ----------- .../EB_ODEQty/run_convergence_test.py | 131 ------------------ 2 files changed, 209 deletions(-) delete mode 100644 Exec/RegTests/EB_ODEQty/input.convergence_0 delete mode 100644 Exec/RegTests/EB_ODEQty/run_convergence_test.py diff --git a/Exec/RegTests/EB_ODEQty/input.convergence_0 b/Exec/RegTests/EB_ODEQty/input.convergence_0 deleted file mode 100644 index 595a6bdd9..000000000 --- a/Exec/RegTests/EB_ODEQty/input.convergence_0 +++ /dev/null @@ -1,78 +0,0 @@ -#----------------------DOMAIN DEFINITION------------------------- -geometry.is_periodic = 0 1 # For each dir, 0: non-perio, 1: periodic -geometry.coord_sys = 0 # 0 => cart, 1 => RZ -geometry.prob_lo = -0.02 -0.02 # x_lo y_lo (z_lo) -geometry.prob_hi = 0.10 0.02 # x_hi y_hi (z_hi) - -#--------------------------BC FLAGS------------------------------ -# Interior, Inflow, Outflow, Symmetry, -# SlipWallAdiab, NoSlipWallAdiab, SlipWallIsotherm, NoSlipWallIsotherm -peleLM.lo_bc = Inflow Interior -peleLM.hi_bc = Outflow Interior - -#-------------------------AMR CONTROL---------------------------- -amr.n_cell = 48 16 # Level 0 number of cells in each direction -amr.v = 1 # AMR verbose -amr.max_level = 0 # maximum level number allowed -amr.ref_ratio = 2 2 2 2 # refinement ratio -amr.regrid_int = 5 # how often to regrid -amr.n_error_buf = 2 2 2 2 # number of buffer cells in error est -amr.grid_eff = 0.7 # what constitutes an efficient grid -amr.blocking_factor = 8 # block factor in grid generation (min box size) -amr.max_grid_size = 64 # max box size - -#--------------------------- Problem ---------------------------- -prob.T_mean = 300.0 -prob.P_mean = 101325.0 -prob.meanFlowMag = 1.0 -prob.meanFlowDir = 1 -prob.ode_IC = 10.0 -prob.ode_xy_lo = -0.01 -prob.ode_length = 0.01 -prob.ode_height = 0.02 -prob.ode_srcstrength = -10.0 - -#-----------------INPUTS TO CONSTANT TRANSPORT ------------------ -transport.units = MKS -transport.const_viscosity = 0.0001 -transport.const_bulk_viscosity = 0.0 -transport.const_conductivity = 0.0 -transport.const_diffusivity = 0.0 - -#---------------------------TIME STEPPING------------------------ -amr.stop_time = 0.001 -amr.fixed_dt = 0.0001 - -#-------------------------PELE CONTROLS---------------------- -peleLM.v = 1 -peleLM.user_defined_ext_sources = 1 - -#---------------------------IO CONTROL--------------------------- -#amr.restart = chk_07136 -amr.check_file = "chk_" -amr.check_per = -1 -amr.plot_file = "plt_conv0_" -amr.plot_per = 0.001 -amr.derive_plot_vars = avg_pressure mag_vort - -#---------------------- Temporal CONTROL ------------------------- -peleLM.do_temporals = 1 # Turn temporals ON/OFF -peleLM.temporal_int = 5 -peleLM.do_extremas = 1 # Compute state extremas - -#------------------------- EB SETUP ----------------------------- -eb2.geom_type = "all_regular" -#eb2.geom_type = sphere -eb2.sphere_radius = 0.005 -eb2.sphere_center = 0.03 0.0 -eb2.sphere_has_fluid_inside = 0 -eb2.small_volfrac = 1.0e-4 - -#--------------------REFINEMENT CONTROL-------------------------- -amr.refinement_indicators = VortL VortH -amr.VortL.max_level = 1 -amr.VortL.value_less = -1000 -amr.VortL.field_name = mag_vort -amr.VortH.max_level = 1 -amr.VortH.value_greater = 1000 -amr.VortH.field_name = mag_vort diff --git a/Exec/RegTests/EB_ODEQty/run_convergence_test.py b/Exec/RegTests/EB_ODEQty/run_convergence_test.py deleted file mode 100644 index 5c7f2e65d..000000000 --- a/Exec/RegTests/EB_ODEQty/run_convergence_test.py +++ /dev/null @@ -1,131 +0,0 @@ -import os -import pandas as pd -import numpy as np -import matplotlib.pyplot as plt - -def generate_input_file(final_time, dt, test, temporal_int, template_filename='input.convergence_0'): - new_input_filename = f'input.convergence_{test}' - - # Read the template input file - with open(template_filename, 'r') as file: - lines = file.readlines() - - # Modify amr.fixed_dt - for i, line in enumerate(lines): - line_start = line.strip().split() - if len(line_start) > 0: - line_start = line_start[0] - if line_start == 'amr.fixed_dt': - lines[i] = f'amr.fixed_dt = {dt}\n' - elif line_start == 'amr.plot_file': - lines[i] = f'amr.plot_file = "plt_conv{test}_"\n' - elif line_start == 'peleLM.temporal_int': - lines[i] = f'peleLM.temporal_int = {temporal_int}\n' - elif line_start == 'amr.stop_time': - lines[i] = f'amr.stop_time = {final_time}\n' - # Write the new input file - with open(new_input_filename, 'w') as file: - file.writelines(lines) - - return new_input_filename - -def run_PeleLMeX(input_file,num_proc): - if (num_proc > 1): - os.system(f'mpirun -np {num_proc} ./PeleLMeX2d.gnu.MPI.ex {input_file}') - else: - os.system(f'./PeleLMeX2d.gnu.MPI.ex {input_file}') - -def dataAnalysis(vars, file_name): - data = pd.read_csv(file_name) - data_vars = data[vars].values - time = data.time.values - IC = data_vars[0, :] - - # Numerical solution at final time - approx = data_vars[-1, :] - - # Exact solution - n = np.arange(len(IC)) - k = -10.0 * 10**(n + 1) - exact = IC * np.exp(k * time[-1]) - - # Absolute error - abs_error = np.abs(exact - approx) - - return abs_error.reshape(1, len(vars)) - - -def removeAllPrevData(): - os.system('rm -r plt* temporals/*') - -# ----------------------------------------------------------------------------- -# Parameters -# ----------------------------------------------------------------------------- - -startWithNewData = 1 # erase all plt files and temporals if true -nRuns = 4 # number of tests -num_proc = 2 # number of processors to run tests -dt_vec = [1e-4, 1e-5, 1e-6, 1e-7] # dt for each test -final_time = 0.001 # final time for simulation - -# Parameters for temporal files -temporals_directory = 'temporals' -temporals_filename = temporals_directory+'/tempExtremas' -num_temporal_write = 2 - -NUM_ODE = 3 -vars = [f"max_MY_ODE_{i}" for i in range(0,NUM_ODE)] -errors = np.zeros((nRuns,NUM_ODE)) - -print("\n==============================================================================") -if startWithNewData: - removeAllPrevData() - -for n in range(0, nRuns): - dt = dt_vec[n] - Nt = round(final_time / dt) - temporal_int = round(Nt / num_temporal_write) - print(f"Running test #{n} with dt = {dt}, Nt = {Nt}, temporal_int = {temporal_int}") - - # Generate a new input file for the given dt - input_file = generate_input_file(final_time, dt, n, temporal_int) - - # Run the convergence test in PeleLMeX - if startWithNewData: - run_PeleLMeX(input_file,num_proc) - - # Rename the tempExtremas file and remove unused tempState file - new_temporals_filename = temporals_filename + f'_{n}' - os.system(f'mv {temporals_filename} {new_temporals_filename}') - os.system('rm temporals/tempState') - - print(f"Finished running test #{n}\n") - - # Read in extrema data from temporals - errors[n,:] = dataAnalysis(vars,new_temporals_filename) - -print("\n==============================================================================\n") - -plt.figure() -for j in range(NUM_ODE): - - # Fit a log-log line to the data - log_dt_vec = np.log(dt_vec[:nRuns]) - log_errors = np.log(errors[:,j]) - slope, intercept = np.polyfit(log_dt_vec, log_errors, 1) - - # Plot scatter plot of errors and the line of best fit - best_fit_line = np.exp(intercept) * np.array(dt_vec[:nRuns])**slope - plt.scatter(dt_vec[:nRuns], errors[:,j], label=f'ODE {j} (rate: {slope:.2f})', marker='o') - plt.plot(dt_vec[:nRuns], best_fit_line, '--') - -# Set x and y scales to log-log -plt.xscale("log") -plt.yscale("log") - -# Add labels and legend -plt.xlabel('$\Delta t$') -plt.ylabel('$\ell_\infty(approx - exact)$') -plt.legend() - -plt.show() \ No newline at end of file From 064562da5a77fd5abed2bf2d4066d983a1d79123 Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Tue, 15 Oct 2024 17:08:21 -0600 Subject: [PATCH 29/46] Minor fix for soot and spray --- Source/PeleLMeX.H | 1 + Source/PeleLMeX_Advance.cpp | 2 +- Source/PeleLMeX_Forces.cpp | 1 + Source/PeleLMeX_ProblemSpecificFunctions.cpp | 2 +- 4 files changed, 4 insertions(+), 2 deletions(-) diff --git a/Source/PeleLMeX.H b/Source/PeleLMeX.H index 977f3864b..d14cae2f7 100644 --- a/Source/PeleLMeX.H +++ b/Source/PeleLMeX.H @@ -879,6 +879,7 @@ public: bool m_user_defined_ext_sources; bool m_ext_sources_SDC; void getExternalSources( + int is_initIter, const PeleLM::TimeStamp& a_timestamp_old, const PeleLM::TimeStamp& a_timestamp_new); diff --git a/Source/PeleLMeX_Advance.cpp b/Source/PeleLMeX_Advance.cpp index d29ca27cd..b011dbbce 100644 --- a/Source/PeleLMeX_Advance.cpp +++ b/Source/PeleLMeX_Advance.cpp @@ -111,7 +111,7 @@ PeleLM::Advance(int is_initIter) //---------------------------------------------------------------- // External sources (soot, radiation, user defined, etc.) - getExternalSources(AmrOldTime, AmrNewTime); + getExternalSources(is_initIter, AmrOldTime, AmrNewTime); if (m_incompressible == 0) { floorSpecies(AmrOldTime); diff --git a/Source/PeleLMeX_Forces.cpp b/Source/PeleLMeX_Forces.cpp index f7de439d3..d3ed0cf46 100644 --- a/Source/PeleLMeX_Forces.cpp +++ b/Source/PeleLMeX_Forces.cpp @@ -240,6 +240,7 @@ PeleLM::addSpark(const TimeStamp& a_timestamp) // Calculate additional external sources (soot, radiation, user defined, etc.) void PeleLM::getExternalSources( + int is_initIter, const PeleLM::TimeStamp& a_timestamp_old, const PeleLM::TimeStamp& a_timestamp_new) { diff --git a/Source/PeleLMeX_ProblemSpecificFunctions.cpp b/Source/PeleLMeX_ProblemSpecificFunctions.cpp index b0a5d95ce..32f5387e5 100644 --- a/Source/PeleLMeX_ProblemSpecificFunctions.cpp +++ b/Source/PeleLMeX_ProblemSpecificFunctions.cpp @@ -23,7 +23,7 @@ set_ode_names(Vector& a_ode_names) #endif } -static void +void problem_modify_ext_sources( Real /*time*/, Real /*dt*/, From ab2f6e2aa0d665ec2d98ccbcb8ae95b54dae8305 Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Tue, 15 Oct 2024 17:57:39 -0600 Subject: [PATCH 30/46] Spelling and unused variable handling. --- Docs/sphinx/manual/Model.rst | 2 +- Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp | 4 ++-- Source/PeleLMeX_Forces.cpp | 3 +++ Source/PeleLMeX_ProblemSpecificFunctions.cpp | 4 ++-- 4 files changed, 8 insertions(+), 5 deletions(-) diff --git a/Docs/sphinx/manual/Model.rst b/Docs/sphinx/manual/Model.rst index b83b55ce3..dba5d68b3 100644 --- a/Docs/sphinx/manual/Model.rst +++ b/Docs/sphinx/manual/Model.rst @@ -118,7 +118,7 @@ For the standard ideal gas EOS, the divergence constraint on velocity becomes: \nabla \cdot \boldsymbol{u} &= \frac{1}{\rho c_p T} \left(\nabla \cdot \lambda\nabla T - \sum_m \boldsymbol{\Gamma_m} \cdot \nabla h_m \right) \\ &- \frac{1}{\rho} \sum_m \frac{W}{W_m} \nabla \cdot \boldsymbol{\Gamma_m} + \frac{1}{\rho}\sum_m \left(\frac{W}{W_m} - \frac{h_m}{c_p T} \right) \dot \omega \equiv S . -In addition to the flow equations, `PeleLMeX` can also solve for a set of quantities that are neither advected nor diffused, satifying: +In addition to the flow equations, `PeleLMeX` can also solve for a set of quantities that are neither advected nor diffused, satisfying: .. math:: diff --git a/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp b/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp index 38139cef3..21e6dffa2 100644 --- a/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp +++ b/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp @@ -12,16 +12,16 @@ Problem specific functions: - Modify as needed */ +#if NUM_ODE > 0 void set_ode_names(Vector& a_ode_names) { -#if NUM_ODE > 0 a_ode_names.resize(NUM_ODE); for (int n = 0; n < NUM_ODE; n++) { a_ode_names[n] = "MY_ODE_" + std::to_string(n); } -#endif } +#endif void problem_modify_ext_sources( diff --git a/Source/PeleLMeX_Forces.cpp b/Source/PeleLMeX_Forces.cpp index d3ed0cf46..55f05e150 100644 --- a/Source/PeleLMeX_Forces.cpp +++ b/Source/PeleLMeX_Forces.cpp @@ -244,6 +244,7 @@ PeleLM::getExternalSources( const PeleLM::TimeStamp& a_timestamp_old, const PeleLM::TimeStamp& a_timestamp_new) { + if (m_n_sparks > 0) { addSpark(a_timestamp_old); } @@ -252,6 +253,8 @@ PeleLM::getExternalSources( if (is_initIter == 0) { SprayMKD(m_cur_time, m_dt); } +#else + (void)is_initIter; #endif #ifdef PELE_USE_SOOT if (do_soot_solve) { diff --git a/Source/PeleLMeX_ProblemSpecificFunctions.cpp b/Source/PeleLMeX_ProblemSpecificFunctions.cpp index 32f5387e5..ff2d4f789 100644 --- a/Source/PeleLMeX_ProblemSpecificFunctions.cpp +++ b/Source/PeleLMeX_ProblemSpecificFunctions.cpp @@ -12,16 +12,16 @@ Problem specific functions: - Modify as needed */ +#if NUM_ODE > 0 void set_ode_names(Vector& a_ode_names) { -#if NUM_ODE > 0 a_ode_names.resize(NUM_ODE); for (int n = 0; n < NUM_ODE; n++) { a_ode_names[n] = "ODE_" + std::to_string(n); } -#endif } +#endif void problem_modify_ext_sources( From c4019653dd00318d10fb4b044b7416cc288509d8 Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Tue, 15 Oct 2024 18:06:05 -0600 Subject: [PATCH 31/46] Unused variable handling --- Source/PeleLMeX_Forces.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/Source/PeleLMeX_Forces.cpp b/Source/PeleLMeX_Forces.cpp index 55f05e150..3023a9c4b 100644 --- a/Source/PeleLMeX_Forces.cpp +++ b/Source/PeleLMeX_Forces.cpp @@ -244,7 +244,8 @@ PeleLM::getExternalSources( const PeleLM::TimeStamp& a_timestamp_old, const PeleLM::TimeStamp& a_timestamp_new) { - + amrex::ignore_unused(is_initIter); + if (m_n_sparks > 0) { addSpark(a_timestamp_old); } @@ -253,8 +254,6 @@ PeleLM::getExternalSources( if (is_initIter == 0) { SprayMKD(m_cur_time, m_dt); } -#else - (void)is_initIter; #endif #ifdef PELE_USE_SOOT if (do_soot_solve) { From d998f5e22ece2636b333925a3b1577b3d55aa89d Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Wed, 16 Oct 2024 08:40:05 -0600 Subject: [PATCH 32/46] Formatting --- Source/PeleLMeX_Forces.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/PeleLMeX_Forces.cpp b/Source/PeleLMeX_Forces.cpp index 3023a9c4b..0db931195 100644 --- a/Source/PeleLMeX_Forces.cpp +++ b/Source/PeleLMeX_Forces.cpp @@ -245,7 +245,7 @@ PeleLM::getExternalSources( const PeleLM::TimeStamp& a_timestamp_new) { amrex::ignore_unused(is_initIter); - + if (m_n_sparks > 0) { addSpark(a_timestamp_old); } From 44e6c4a59e7354ba4f6511f581620af6a2faf794 Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Wed, 16 Oct 2024 15:24:45 -0600 Subject: [PATCH 33/46] Removed unwanted print statements, renamed some files/vars for consistency. --- CMake/BuildPeleExe.cmake | 6 ++++++ Exec/Make.PeleLMeX | 4 ++-- Exec/RegTests/EB_ODEQty/GNUmakefile | 2 +- Exec/RegTests/EB_ODEQty/{input.2d => eb-odeqty-2d.inp} | 0 Source/Make.package | 2 +- Source/PeleLMeX_Plot.cpp | 3 --- 6 files changed, 10 insertions(+), 7 deletions(-) rename Exec/RegTests/EB_ODEQty/{input.2d => eb-odeqty-2d.inp} (100%) diff --git a/CMake/BuildPeleExe.cmake b/CMake/BuildPeleExe.cmake index 46e4114b3..6ad7eecac 100644 --- a/CMake/BuildPeleExe.cmake +++ b/CMake/BuildPeleExe.cmake @@ -90,6 +90,12 @@ function(build_pele_exe pele_exe_name pele_physics_lib_name) ${SRC_DIR}/PeleLMeX_Radiation.cpp) endif() + if(PELELM_NUM_ODE GREATER 0) + target_sources(${pele_exe_name} + PRIVATE + NUM_ODE=${PELELM_NUM_ODE}) + endif() + if(NOT "${pele_exe_name}" STREQUAL "${PROJECT_NAME}-UnitTests") target_sources(${pele_exe_name} PRIVATE diff --git a/Exec/Make.PeleLMeX b/Exec/Make.PeleLMeX index 4ef1c7734..6da5be0c5 100644 --- a/Exec/Make.PeleLMeX +++ b/Exec/Make.PeleLMeX @@ -57,8 +57,8 @@ ifeq ($(USE_RADIATION), TRUE) DEFINES += -DPELE_USE_RADIATION endif -ifeq ($(shell test 0$(NUM_ODE) -gt 0; echo $$?), 0) - DEFINES += -DNUM_ODE=$(NUM_ODE) +ifeq ($(shell test 0$(PELELM_NUM_ODE) -gt 0; echo $$?), 0) + DEFINES += -DNUM_ODE=$(PELELM_NUM_ODE) endif Bpack += $(foreach dir, $(LMdirs), $(PELE_HOME)/$(dir)/Make.package) diff --git a/Exec/RegTests/EB_ODEQty/GNUmakefile b/Exec/RegTests/EB_ODEQty/GNUmakefile index e97efd845..ea18ca64b 100644 --- a/Exec/RegTests/EB_ODEQty/GNUmakefile +++ b/Exec/RegTests/EB_ODEQty/GNUmakefile @@ -28,7 +28,7 @@ THREAD_SANITIZER = FALSE # PeleLMeX CEXE_headers += EBUserDefined.H CEXE_sources += PeleLMeX_ProblemSpecificFunctions.cpp -NUM_ODE = 3 +PELELM_NUM_ODE = 3 # PelePhysics Chemistry_Model = air diff --git a/Exec/RegTests/EB_ODEQty/input.2d b/Exec/RegTests/EB_ODEQty/eb-odeqty-2d.inp similarity index 100% rename from Exec/RegTests/EB_ODEQty/input.2d rename to Exec/RegTests/EB_ODEQty/eb-odeqty-2d.inp diff --git a/Source/Make.package b/Source/Make.package index 821172eff..313e9ee75 100644 --- a/Source/Make.package +++ b/Source/Make.package @@ -52,7 +52,7 @@ ifeq ($(USE_SOOT), TRUE) CEXE_sources += PeleLMeX_Soot.cpp endif -ifeq ($(shell test 0$(NUM_ODE) -gt 0; echo $$?), 0) +ifeq ($(shell test 0$(PELELM_NUM_ODE) -gt 0; echo $$?), 0) CEXE_sources += PeleLMeX_ODEQty.cpp endif diff --git a/Source/PeleLMeX_Plot.cpp b/Source/PeleLMeX_Plot.cpp index 0d559427e..a6634b4d6 100644 --- a/Source/PeleLMeX_Plot.cpp +++ b/Source/PeleLMeX_Plot.cpp @@ -287,8 +287,6 @@ PeleLM::WritePlotFile() } #endif - Print() << PrettyLine; - //---------------------------------------------------------------- // Fill the plot MultiFabs for (int lev = 0; lev <= finest_level; ++lev) { @@ -413,7 +411,6 @@ PeleLM::WritePlotFile() MultiFab::Copy( mf_plt[lev], m_leveldata_new[lev]->state, FIRSTODE, cnt, NUM_ODE, 0); cnt += NUM_ODE; - Print() << "cnt = " << cnt << std::endl; #endif if (m_do_les && m_plot_les) { From 838ca88583295f30541a0a428cd6b4c568ca8783 Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Wed, 16 Oct 2024 15:36:53 -0600 Subject: [PATCH 34/46] Added case to CI --- Exec/RegTests/EB_ODEQty/CMakeLists.txt | 1 + Tests/CMakeLists.txt | 8 ++++++++ 2 files changed, 9 insertions(+) diff --git a/Exec/RegTests/EB_ODEQty/CMakeLists.txt b/Exec/RegTests/EB_ODEQty/CMakeLists.txt index 7169cbdba..b1f0ab7b3 100644 --- a/Exec/RegTests/EB_ODEQty/CMakeLists.txt +++ b/Exec/RegTests/EB_ODEQty/CMakeLists.txt @@ -5,4 +5,5 @@ set(PELE_PHYSICS_ENABLE_SPRAY OFF) set(PELE_PHYSICS_SPRAY_FUEL_NUM 0) set(PELE_PHYSICS_ENABLE_SOOT OFF) set(PELE_PHYSICS_ENABLE_RADIATION OFF) +set(PELELM_NUM_ODE 3) include(BuildExeAndLib) diff --git a/Tests/CMakeLists.txt b/Tests/CMakeLists.txt index 58c53bc79..6d62f0f81 100644 --- a/Tests/CMakeLists.txt +++ b/Tests/CMakeLists.txt @@ -225,3 +225,11 @@ if(NOT PELE_ENABLE_EB) add_test_r(hit-les-${PELE_DIM}d HITDecay) endif() endif() + +if (PELE_ENABLE_EB) + if(PELE_DIM EQUAL 2) + if(PELELM_NUM_ODE GREATER 0) + add_test_r(eb-odeqty-${PELE_DIM}d EB_ODEQty) + endif() + endif() +endif() \ No newline at end of file From 88ca2d01a2e458697e432b0465e3104ad8b0243f Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Thu, 17 Oct 2024 08:03:41 -0600 Subject: [PATCH 35/46] Add USE_EB to EB_ODEQty/CmakeLists.txt --- Exec/RegTests/EB_ODEQty/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/Exec/RegTests/EB_ODEQty/CMakeLists.txt b/Exec/RegTests/EB_ODEQty/CMakeLists.txt index b1f0ab7b3..0283fda26 100644 --- a/Exec/RegTests/EB_ODEQty/CMakeLists.txt +++ b/Exec/RegTests/EB_ODEQty/CMakeLists.txt @@ -5,5 +5,6 @@ set(PELE_PHYSICS_ENABLE_SPRAY OFF) set(PELE_PHYSICS_SPRAY_FUEL_NUM 0) set(PELE_PHYSICS_ENABLE_SOOT OFF) set(PELE_PHYSICS_ENABLE_RADIATION OFF) +set(PELE_USE_EB ON) set(PELELM_NUM_ODE 3) include(BuildExeAndLib) From afbc707014d60764bd1869f2deeb64c3c46d6c65 Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Thu, 17 Oct 2024 08:12:15 -0600 Subject: [PATCH 36/46] Add EB_ODEQty to RegTests/CMakeLists.txt --- Exec/RegTests/CMakeLists.txt | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Exec/RegTests/CMakeLists.txt b/Exec/RegTests/CMakeLists.txt index d8ac60c5b..ec2826c50 100644 --- a/Exec/RegTests/CMakeLists.txt +++ b/Exec/RegTests/CMakeLists.txt @@ -6,6 +6,9 @@ if(PELE_ENABLE_EB) if(PELE_DIM EQUAL 3) add_subdirectory(EB_PipeFlow) endif() + if(PELE_DIM EQUAL 2) + add_subdirectory(EB_ODEQty) + endif() else() add_subdirectory(EnclosedFlame) add_subdirectory(EnclosedInjection) From fcdd8f31d42d068ca810d8c34885fb3ac405597d Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Thu, 17 Oct 2024 08:49:32 -0600 Subject: [PATCH 37/46] CI fix --- CMake/BuildPeleExe.cmake | 6 ------ Exec/RegTests/EB_ODEQty/CMakeLists.txt | 1 - 2 files changed, 7 deletions(-) diff --git a/CMake/BuildPeleExe.cmake b/CMake/BuildPeleExe.cmake index 6ad7eecac..46e4114b3 100644 --- a/CMake/BuildPeleExe.cmake +++ b/CMake/BuildPeleExe.cmake @@ -90,12 +90,6 @@ function(build_pele_exe pele_exe_name pele_physics_lib_name) ${SRC_DIR}/PeleLMeX_Radiation.cpp) endif() - if(PELELM_NUM_ODE GREATER 0) - target_sources(${pele_exe_name} - PRIVATE - NUM_ODE=${PELELM_NUM_ODE}) - endif() - if(NOT "${pele_exe_name}" STREQUAL "${PROJECT_NAME}-UnitTests") target_sources(${pele_exe_name} PRIVATE diff --git a/Exec/RegTests/EB_ODEQty/CMakeLists.txt b/Exec/RegTests/EB_ODEQty/CMakeLists.txt index 0283fda26..b1f0ab7b3 100644 --- a/Exec/RegTests/EB_ODEQty/CMakeLists.txt +++ b/Exec/RegTests/EB_ODEQty/CMakeLists.txt @@ -5,6 +5,5 @@ set(PELE_PHYSICS_ENABLE_SPRAY OFF) set(PELE_PHYSICS_SPRAY_FUEL_NUM 0) set(PELE_PHYSICS_ENABLE_SOOT OFF) set(PELE_PHYSICS_ENABLE_RADIATION OFF) -set(PELE_USE_EB ON) set(PELELM_NUM_ODE 3) include(BuildExeAndLib) From e070e3d5a487a961ed3a50bfa69c4bbbe0a8986a Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Thu, 17 Oct 2024 12:13:56 -0600 Subject: [PATCH 38/46] More CI --- CMake/BuildPeleExe.cmake | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/CMake/BuildPeleExe.cmake b/CMake/BuildPeleExe.cmake index 46e4114b3..ea99f1ec6 100644 --- a/CMake/BuildPeleExe.cmake +++ b/CMake/BuildPeleExe.cmake @@ -97,6 +97,10 @@ function(build_pele_exe pele_exe_name pele_physics_lib_name) ) endif() + if (PELELM_NUM_ODE GREATER 0) + target_compile_definitions(${pele_exe_name} PRIVATE NUM_ODE=${PELELM_NUM_ODE}) + endif() + if(PELE_ENABLE_CUDA) set(pctargets "${pele_exe_name};${pele_physics_lib_name}") foreach(tgt IN LISTS pctargets) From 9707476f95e5420adf4efec59f1823034cc97165 Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Thu, 17 Oct 2024 12:37:43 -0600 Subject: [PATCH 39/46] Add ODEQty.cpp to cmake --- CMake/BuildPeleExe.cmake | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/CMake/BuildPeleExe.cmake b/CMake/BuildPeleExe.cmake index ea99f1ec6..9e81f1912 100644 --- a/CMake/BuildPeleExe.cmake +++ b/CMake/BuildPeleExe.cmake @@ -46,9 +46,10 @@ function(build_pele_exe pele_exe_name pele_physics_lib_name) ${SRC_DIR}/PeleLMeX_Evolve.cpp ${SRC_DIR}/PeleLMeX_FlowController.cpp ${SRC_DIR}/PeleLMeX_Forces.cpp + ${SRC_DIR}/PeleLMeX_Init.cpp + ${SRC_DIR}/PeleLMeX_ODEQty.cpp ${SRC_DIR}/PeleLMeX_PatchFlowVariables.H ${SRC_DIR}/PeleLMeX_PatchFlowVariables.cpp - ${SRC_DIR}/PeleLMeX_Init.cpp ${SRC_DIR}/PeleLMeX_Plot.cpp ${SRC_DIR}/PeleLMeX_ProblemSpecificFunctions.H ${SRC_DIR}/PeleLMeX_ProblemSpecificFunctions.cpp From 9728a0f4e0b0fb30d5774741041344041812c10f Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Thu, 17 Oct 2024 13:21:19 -0600 Subject: [PATCH 40/46] More CI... --- CMake/BuildPeleExe.cmake | 8 ++++---- Exec/RegTests/CMakeLists.txt | 4 +++- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/CMake/BuildPeleExe.cmake b/CMake/BuildPeleExe.cmake index 9e81f1912..e04c85dbd 100644 --- a/CMake/BuildPeleExe.cmake +++ b/CMake/BuildPeleExe.cmake @@ -91,6 +91,10 @@ function(build_pele_exe pele_exe_name pele_physics_lib_name) ${SRC_DIR}/PeleLMeX_Radiation.cpp) endif() + if (PELELM_NUM_ODE GREATER 0) + target_compile_definitions(${pele_exe_name} PRIVATE NUM_ODE=${PELELM_NUM_ODE}) + endif() + if(NOT "${pele_exe_name}" STREQUAL "${PROJECT_NAME}-UnitTests") target_sources(${pele_exe_name} PRIVATE @@ -98,10 +102,6 @@ function(build_pele_exe pele_exe_name pele_physics_lib_name) ) endif() - if (PELELM_NUM_ODE GREATER 0) - target_compile_definitions(${pele_exe_name} PRIVATE NUM_ODE=${PELELM_NUM_ODE}) - endif() - if(PELE_ENABLE_CUDA) set(pctargets "${pele_exe_name};${pele_physics_lib_name}") foreach(tgt IN LISTS pctargets) diff --git a/Exec/RegTests/CMakeLists.txt b/Exec/RegTests/CMakeLists.txt index ec2826c50..4ea0d95c4 100644 --- a/Exec/RegTests/CMakeLists.txt +++ b/Exec/RegTests/CMakeLists.txt @@ -7,7 +7,9 @@ if(PELE_ENABLE_EB) add_subdirectory(EB_PipeFlow) endif() if(PELE_DIM EQUAL 2) - add_subdirectory(EB_ODEQty) + if(PELELM_NUM_ODE GREATER 0) + add_subdirectory(EB_ODEQty) + endif() endif() else() add_subdirectory(EnclosedFlame) From b0f9f350f888a228ce521260559df1a3ea772c07 Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Thu, 17 Oct 2024 16:37:04 -0600 Subject: [PATCH 41/46] Only define predictODEQty if NUM_ODE > 0 --- Source/PeleLMeX_ODEQty.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Source/PeleLMeX_ODEQty.cpp b/Source/PeleLMeX_ODEQty.cpp index 72e846acb..44bd75df2 100644 --- a/Source/PeleLMeX_ODEQty.cpp +++ b/Source/PeleLMeX_ODEQty.cpp @@ -3,6 +3,7 @@ using namespace amrex; +#if NUM_ODE > 0 void PeleLM::predictODEQty() { @@ -22,4 +23,5 @@ PeleLM::predictODEQty() }); Gpu::streamSynchronize(); } -} \ No newline at end of file +} +#endif \ No newline at end of file From 3e32ab3556efe51c4963cd77a8222ea6e53ba0fb Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Fri, 18 Oct 2024 15:43:36 -0600 Subject: [PATCH 42/46] Updated CMake to handle user defined EB's and ProblemSpecificFunctions --- CMake/BuildPeleExe.cmake | 22 +++++- Exec/RegTests/CMakeLists.txt | 6 +- .../EB_BackwardStepFlame/CMakeLists.txt | 1 + Exec/RegTests/EB_EnclosedFlame/CMakeLists.txt | 1 + .../RegTests/EB_EnclosedVortex/CMakeLists.txt | 1 + .../EB_FlowPastCylinder/CMakeLists.txt | 1 + Exec/RegTests/EB_ODEQty/CMakeLists.txt | 2 + Exec/RegTests/EB_ODEQty/GNUmakefile | 1 - .../EB_ODEQty/PeleLMeX_EBUserDefined.H | 58 -------------- .../PeleLMeX_ProblemSpecificFunctions.cpp | 2 +- Exec/RegTests/EB_ODEQty/eb-odeqty-2d.inp | 4 +- Exec/RegTests/EB_ODEQty/eb-odeqty-3d.inp | 76 +++++++++++++++++++ Exec/RegTests/EB_PipeFlow/CMakeLists.txt | 1 + Tests/CMakeLists.txt | 10 +-- 14 files changed, 109 insertions(+), 77 deletions(-) delete mode 100644 Exec/RegTests/EB_ODEQty/PeleLMeX_EBUserDefined.H create mode 100644 Exec/RegTests/EB_ODEQty/eb-odeqty-3d.inp diff --git a/CMake/BuildPeleExe.cmake b/CMake/BuildPeleExe.cmake index e04c85dbd..aaa164a56 100644 --- a/CMake/BuildPeleExe.cmake +++ b/CMake/BuildPeleExe.cmake @@ -25,7 +25,6 @@ function(build_pele_exe pele_exe_name pele_physics_lib_name) ${SRC_DIR}/PeleLMeX_DeriveUserDefined.cpp ${SRC_DIR}/PeleLMeX_DiffusionOp.H ${SRC_DIR}/PeleLMeX_DiffusionOp.cpp - ${SRC_DIR}/PeleLMeX_EBUserDefined.H ${SRC_DIR}/PeleLMeX_FlowControllerData.H ${SRC_DIR}/PeleLMeX.H ${SRC_DIR}/PeleLMeX.cpp @@ -52,7 +51,6 @@ function(build_pele_exe pele_exe_name pele_physics_lib_name) ${SRC_DIR}/PeleLMeX_PatchFlowVariables.cpp ${SRC_DIR}/PeleLMeX_Plot.cpp ${SRC_DIR}/PeleLMeX_ProblemSpecificFunctions.H - ${SRC_DIR}/PeleLMeX_ProblemSpecificFunctions.cpp ${SRC_DIR}/PeleLMeX_Projection.cpp ${SRC_DIR}/PeleLMeX_Reactions.cpp ${SRC_DIR}/PeleLMeX_Regrid.cpp @@ -72,6 +70,26 @@ function(build_pele_exe pele_exe_name pele_physics_lib_name) ${SRC_DIR}/main.cpp ) + if(PELE_EB_USER_DEFINED) + target_sources(${pele_exe_name} + PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/PeleLMeX_EBUserDefined.H) + else() + target_sources(${pele_exe_name} + PRIVATE + ${SRC_DIR}/PeleLMeX_EBUserDefined.H) + endif() + + if(PELELM_USER_DEFINED_EXT_SRC) + target_sources(${pele_exe_name} + PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/PeleLMeX_ProblemSpecificFunctions.cpp) + else() + target_sources(${pele_exe_name} + PRIVATE + ${SRC_DIR}/PeleLMeX_ProblemSpecificFunctions.cpp) + endif() + if(PELE_PHYSICS_ENABLE_SOOT) target_sources(${pele_exe_name} PRIVATE diff --git a/Exec/RegTests/CMakeLists.txt b/Exec/RegTests/CMakeLists.txt index 4ea0d95c4..88474eff5 100644 --- a/Exec/RegTests/CMakeLists.txt +++ b/Exec/RegTests/CMakeLists.txt @@ -3,14 +3,10 @@ if(PELE_ENABLE_EB) add_subdirectory(EB_EnclosedFlame) add_subdirectory(EB_EnclosedVortex) add_subdirectory(EB_FlowPastCylinder) + add_subdirectory(EB_ODEQty) if(PELE_DIM EQUAL 3) add_subdirectory(EB_PipeFlow) endif() - if(PELE_DIM EQUAL 2) - if(PELELM_NUM_ODE GREATER 0) - add_subdirectory(EB_ODEQty) - endif() - endif() else() add_subdirectory(EnclosedFlame) add_subdirectory(EnclosedInjection) diff --git a/Exec/RegTests/EB_BackwardStepFlame/CMakeLists.txt b/Exec/RegTests/EB_BackwardStepFlame/CMakeLists.txt index cd74295e8..a01666cf8 100644 --- a/Exec/RegTests/EB_BackwardStepFlame/CMakeLists.txt +++ b/Exec/RegTests/EB_BackwardStepFlame/CMakeLists.txt @@ -5,4 +5,5 @@ set(PELE_PHYSICS_ENABLE_SPRAY OFF) set(PELE_PHYSICS_SPRAY_FUEL_NUM 0) set(PELE_PHYSICS_ENABLE_SOOT OFF) set(PELE_PHYSICS_ENABLE_RADIATION OFF) +set(PELE_EB_USER_DEFINED ON) include(BuildExeAndLib) diff --git a/Exec/RegTests/EB_EnclosedFlame/CMakeLists.txt b/Exec/RegTests/EB_EnclosedFlame/CMakeLists.txt index cd74295e8..cf15fe304 100644 --- a/Exec/RegTests/EB_EnclosedFlame/CMakeLists.txt +++ b/Exec/RegTests/EB_EnclosedFlame/CMakeLists.txt @@ -5,4 +5,5 @@ set(PELE_PHYSICS_ENABLE_SPRAY OFF) set(PELE_PHYSICS_SPRAY_FUEL_NUM 0) set(PELE_PHYSICS_ENABLE_SOOT OFF) set(PELE_PHYSICS_ENABLE_RADIATION OFF) +set(PELE_EB_USER_DEFINED OFF) include(BuildExeAndLib) diff --git a/Exec/RegTests/EB_EnclosedVortex/CMakeLists.txt b/Exec/RegTests/EB_EnclosedVortex/CMakeLists.txt index 7169cbdba..2d2403318 100644 --- a/Exec/RegTests/EB_EnclosedVortex/CMakeLists.txt +++ b/Exec/RegTests/EB_EnclosedVortex/CMakeLists.txt @@ -5,4 +5,5 @@ set(PELE_PHYSICS_ENABLE_SPRAY OFF) set(PELE_PHYSICS_SPRAY_FUEL_NUM 0) set(PELE_PHYSICS_ENABLE_SOOT OFF) set(PELE_PHYSICS_ENABLE_RADIATION OFF) +set(PELE_EB_USER_DEFINED OFF) include(BuildExeAndLib) diff --git a/Exec/RegTests/EB_FlowPastCylinder/CMakeLists.txt b/Exec/RegTests/EB_FlowPastCylinder/CMakeLists.txt index 7169cbdba..6a10c2d9a 100644 --- a/Exec/RegTests/EB_FlowPastCylinder/CMakeLists.txt +++ b/Exec/RegTests/EB_FlowPastCylinder/CMakeLists.txt @@ -5,4 +5,5 @@ set(PELE_PHYSICS_ENABLE_SPRAY OFF) set(PELE_PHYSICS_SPRAY_FUEL_NUM 0) set(PELE_PHYSICS_ENABLE_SOOT OFF) set(PELE_PHYSICS_ENABLE_RADIATION OFF) +set(PELE_EB_USER_DEFINED ON) include(BuildExeAndLib) diff --git a/Exec/RegTests/EB_ODEQty/CMakeLists.txt b/Exec/RegTests/EB_ODEQty/CMakeLists.txt index b1f0ab7b3..06198d3ca 100644 --- a/Exec/RegTests/EB_ODEQty/CMakeLists.txt +++ b/Exec/RegTests/EB_ODEQty/CMakeLists.txt @@ -5,5 +5,7 @@ set(PELE_PHYSICS_ENABLE_SPRAY OFF) set(PELE_PHYSICS_SPRAY_FUEL_NUM 0) set(PELE_PHYSICS_ENABLE_SOOT OFF) set(PELE_PHYSICS_ENABLE_RADIATION OFF) +set(PELE_EB_USER_DEFINED OFF) +set(PELELM_USER_DEFINED_EXT_SRC ON) set(PELELM_NUM_ODE 3) include(BuildExeAndLib) diff --git a/Exec/RegTests/EB_ODEQty/GNUmakefile b/Exec/RegTests/EB_ODEQty/GNUmakefile index ea18ca64b..acda3985c 100644 --- a/Exec/RegTests/EB_ODEQty/GNUmakefile +++ b/Exec/RegTests/EB_ODEQty/GNUmakefile @@ -26,7 +26,6 @@ FSANITIZER = FALSE THREAD_SANITIZER = FALSE # PeleLMeX -CEXE_headers += EBUserDefined.H CEXE_sources += PeleLMeX_ProblemSpecificFunctions.cpp PELELM_NUM_ODE = 3 diff --git a/Exec/RegTests/EB_ODEQty/PeleLMeX_EBUserDefined.H b/Exec/RegTests/EB_ODEQty/PeleLMeX_EBUserDefined.H deleted file mode 100644 index 9de8406f9..000000000 --- a/Exec/RegTests/EB_ODEQty/PeleLMeX_EBUserDefined.H +++ /dev/null @@ -1,58 +0,0 @@ -#ifndef EBUSERDEFINED_H -#define EBUSERDEFINED_H - -using namespace amrex; - -#ifdef AMREX_USE_EB -#include -#include -void -EBUserDefined( - const Geometry& /*geom*/, - const int /*required_coarsening_level*/, - const int /*max_coarsening_level*/) -{ - // ParmParse your geometry parameters - - // Build geometry pieces using EB2::* methods - - // Build your geometry shop using EB2::makeShop - - // Build geom using EB2::Build - - // We shoulnd't be here, copy this file in your run folder - // and implement your geometry - Abort("In default EBUserDefined function! Shouldn't be here. Copy and edit " - "this file for your needs"); -} - -AMREX_GPU_DEVICE -AMREX_FORCE_INLINE -void -setEBState( - const amrex::Real xEBface[AMREX_SPACEDIM], - amrex::Real s_ext[NVAR], - const amrex::Real /*time*/, - amrex::GeometryData const& /*geomdata*/, - ProbParm const& /*prob_parm*/) -{ - if (xEBface[1] > 0.02) { - s_ext[TEMP] = 500.0; - } else { - s_ext[TEMP] = 300.0; - } -} - -AMREX_GPU_DEVICE -AMREX_FORCE_INLINE -void -setEBType( - const amrex::Real* /*xEBface[AMREX_SPACEDIM]*/, - amrex::Real& EBflagType, - amrex::GeometryData const& /*geomdata*/, - ProbParm const& /*prob_parm*/) -{ - EBflagType = 1.0; -} -#endif -#endif diff --git a/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp b/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp index 21e6dffa2..5c06c6cdc 100644 --- a/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp +++ b/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp @@ -31,7 +31,7 @@ problem_modify_ext_sources( MultiArray4 const& state_old_arr, MultiArray4 const& /*state_new_arr*/, Vector>& a_extSource, - const GeometryData& geomdata, + const GeometryData& /*geomdata*/, ProbParm const& prob_parm) { /* diff --git a/Exec/RegTests/EB_ODEQty/eb-odeqty-2d.inp b/Exec/RegTests/EB_ODEQty/eb-odeqty-2d.inp index 37a46d565..5ae43561b 100644 --- a/Exec/RegTests/EB_ODEQty/eb-odeqty-2d.inp +++ b/Exec/RegTests/EB_ODEQty/eb-odeqty-2d.inp @@ -14,8 +14,8 @@ peleLM.hi_bc = Outflow Interior amr.n_cell = 96 32 # Level 0 number of cells in each direction amr.v = 1 # AMR verbose amr.max_level = 0 # maximum level number allowed -amr.ref_ratio = 2 2 2 2 # refinement ratio -amr.regrid_int = 5 # how often to regrid +#amr.ref_ratio = 2 2 2 2 # refinement ratio +#amr.regrid_int = 5 # how often to regrid amr.n_error_buf = 2 2 2 2 # number of buffer cells in error est amr.grid_eff = 0.7 # what constitutes an efficient grid amr.blocking_factor = 16 # block factor in grid generation (min box size) diff --git a/Exec/RegTests/EB_ODEQty/eb-odeqty-3d.inp b/Exec/RegTests/EB_ODEQty/eb-odeqty-3d.inp new file mode 100644 index 000000000..d839ee738 --- /dev/null +++ b/Exec/RegTests/EB_ODEQty/eb-odeqty-3d.inp @@ -0,0 +1,76 @@ +#----------------------DOMAIN DEFINITION------------------------- +geometry.is_periodic = 0 1 1 # For each dir, 0: non-perio, 1: periodic +geometry.coord_sys = 0 # 0 => cart, 1 => RZ +geometry.prob_lo = -0.02 -0.02 -0.01 # x_lo y_lo z_lo +geometry.prob_hi = 0.10 0.02 0.01 # x_hi y_hi z_hi + +#--------------------------BC FLAGS------------------------------ +# Interior, Inflow, Outflow, Symmetry, +# SlipWallAdiab, NoSlipWallAdiab, SlipWallIsotherm, NoSlipWallIsotherm +peleLM.lo_bc = Inflow Interior Interior +peleLM.hi_bc = Outflow Interior Interior + +#-------------------------AMR CONTROL---------------------------- +amr.n_cell = 48 16 8 # Level 0 number of cells in each direction +amr.v = 1 # AMR verbose +amr.max_level = 1 # maximum level number allowed +#amr.ref_ratio = 2 2 2 2 # refinement ratio +#amr.regrid_int = 5 # how often to regrid +amr.n_error_buf = 2 2 2 2 # number of buffer cells in error est +amr.grid_eff = 0.7 # what constitutes an efficient grid +amr.blocking_factor = 8 # block factor in grid generation (min box size) +amr.max_grid_size = 64 # max box size + +#--------------------------- Problem ---------------------------- +prob.T_mean = 300.0 +prob.P_mean = 101325.0 +prob.meanFlowMag = 4.255 +prob.meanFlowDir = 1 +prob.ode_IC = 10.0 +prob.ode_xy_lo = -0.01 +prob.ode_length = 0.01 +prob.ode_height = 0.02 +prob.ode_srcstrength = -10.0 + +#-----------------INPUTS TO CONSTANT TRANSPORT ------------------ +transport.units = MKS +transport.const_viscosity = 0.0001 +transport.const_bulk_viscosity = 0.0 +transport.const_conductivity = 0.0 +transport.const_diffusivity = 0.0 + +#---------------------------TIME STEPPING------------------------ +#amr.max_step = 1 +amr.stop_time = 0.02 +#amr.dt_shrink = 0.1 +amr.dt_change_max = 1.1 +amr.fixed_dt = 0.0001 + +#-------------------------PELE CONTROLS---------------------- +peleLM.v = 1 +peleLM.user_defined_ext_sources = 1 + +#---------------------------IO CONTROL--------------------------- +#amr.restart = chk_07136 +amr.check_file = "chk_" +amr.check_per = -1 +amr.plot_file = "plt_" +amr.plot_per = 0.001 +amr.derive_plot_vars = avg_pressure mag_vort + +#------------------------- EB SETUP ----------------------------- +eb2.geom_type = cylinder +eb2.cylinder_radius = 0.005 +eb2.cylinder_direction = 2 +eb2.cylinder_center = 0.03 0.0 0.0 +eb2.cylinder_has_fluid_inside = 0 +eb2.small_volfrac = 1.0e-4 + +#--------------------REFINEMENT CONTROL-------------------------- +amr.refinement_indicators = VortL VortH +amr.VortL.max_level = 1 +amr.VortL.value_less = -1000 +amr.VortL.field_name = mag_vort +amr.VortH.max_level = 1 +amr.VortH.value_greater = 1000 +amr.VortH.field_name = mag_vort diff --git a/Exec/RegTests/EB_PipeFlow/CMakeLists.txt b/Exec/RegTests/EB_PipeFlow/CMakeLists.txt index 66d2fcdbb..44d49cf76 100644 --- a/Exec/RegTests/EB_PipeFlow/CMakeLists.txt +++ b/Exec/RegTests/EB_PipeFlow/CMakeLists.txt @@ -5,4 +5,5 @@ set(PELE_PHYSICS_ENABLE_SPRAY OFF) set(PELE_PHYSICS_SPRAY_FUEL_NUM 0) set(PELE_PHYSICS_ENABLE_SOOT OFF) set(PELE_PHYSICS_ENABLE_RADIATION OFF) +set(PELE_EB_USER_DEFINED ON) include(BuildExeAndLib) diff --git a/Tests/CMakeLists.txt b/Tests/CMakeLists.txt index 6d62f0f81..330b644c6 100644 --- a/Tests/CMakeLists.txt +++ b/Tests/CMakeLists.txt @@ -224,12 +224,6 @@ if(NOT PELE_ENABLE_EB) add_test_r(hit-${PELE_DIM}d HITDecay) add_test_r(hit-les-${PELE_DIM}d HITDecay) endif() -endif() - -if (PELE_ENABLE_EB) - if(PELE_DIM EQUAL 2) - if(PELELM_NUM_ODE GREATER 0) - add_test_r(eb-odeqty-${PELE_DIM}d EB_ODEQty) - endif() - endif() +else() + add_test_r(eb-odeqty-${PELE_DIM}d EB_ODEQty) endif() \ No newline at end of file From 19279b101bae1830e742f90489f2058903825269 Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Mon, 21 Oct 2024 16:19:12 -0600 Subject: [PATCH 43/46] Now passes in MultiFabs for state_old and state_new. This provides access to state_old.Factory() etc. --- .../EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp | 7 ++++--- Source/PeleLMeX_Forces.cpp | 8 ++++---- Source/PeleLMeX_ProblemSpecificFunctions.H | 6 +++--- Source/PeleLMeX_ProblemSpecificFunctions.cpp | 12 +++++++----- 4 files changed, 18 insertions(+), 15 deletions(-) diff --git a/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp b/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp index 5c06c6cdc..d7669544e 100644 --- a/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp +++ b/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp @@ -28,11 +28,11 @@ problem_modify_ext_sources( Real /*time*/, Real /*dt*/, int lev, - MultiArray4 const& state_old_arr, - MultiArray4 const& /*state_new_arr*/, + const MultiFab& state_old, + const MultiFab& /*state_new*/, Vector>& a_extSource, const GeometryData& /*geomdata*/, - ProbParm const& prob_parm) + const ProbParm& prob_parm) { /* Notes: @@ -42,6 +42,7 @@ problem_modify_ext_sources( */ auto ext_source_arr = a_extSource[lev]->arrays(); + auto const& state_old_arr = state_old.const_arrays(); ParallelFor( *a_extSource[lev], diff --git a/Source/PeleLMeX_Forces.cpp b/Source/PeleLMeX_Forces.cpp index 0db931195..164fac2ff 100644 --- a/Source/PeleLMeX_Forces.cpp +++ b/Source/PeleLMeX_Forces.cpp @@ -271,11 +271,11 @@ PeleLM::getExternalSources( // User defined external sources if (m_user_defined_ext_sources) { for (int lev = 0; lev <= finest_level; lev++) { + auto* ldata_p_old = getLevelDataPtr(lev, a_timestamp_old); + auto* ldata_p_new = getLevelDataPtr(lev, a_timestamp_new); problem_modify_ext_sources( - getTime(lev, a_timestamp_new), m_dt, lev, - getLevelDataPtr(lev, a_timestamp_old)->state.const_arrays(), - getLevelDataPtr(lev, a_timestamp_new)->state.const_arrays(), - m_extSource, geom[lev].data(), *prob_parm_d); + getTime(lev, a_timestamp_new), m_dt, lev, ldata_p_old->state, + ldata_p_new->state, m_extSource, geom[lev].data(), *prob_parm_d); } } } \ No newline at end of file diff --git a/Source/PeleLMeX_ProblemSpecificFunctions.H b/Source/PeleLMeX_ProblemSpecificFunctions.H index 697064b3b..534c19344 100644 --- a/Source/PeleLMeX_ProblemSpecificFunctions.H +++ b/Source/PeleLMeX_ProblemSpecificFunctions.H @@ -11,9 +11,9 @@ void problem_modify_ext_sources( Real time, Real dt, int lev, - MultiArray4 const& state_old, - MultiArray4 const& state_new, + const MultiFab& state_old, + const MultiFab& state_new, Vector>& a_extSource, const GeometryData& geomdata, - ProbParm const& prob_parm); + const ProbParm& prob_parm); #endif \ No newline at end of file diff --git a/Source/PeleLMeX_ProblemSpecificFunctions.cpp b/Source/PeleLMeX_ProblemSpecificFunctions.cpp index ff2d4f789..e228a088a 100644 --- a/Source/PeleLMeX_ProblemSpecificFunctions.cpp +++ b/Source/PeleLMeX_ProblemSpecificFunctions.cpp @@ -28,11 +28,11 @@ problem_modify_ext_sources( Real /*time*/, Real /*dt*/, int /*lev*/, - MultiArray4 const& /*state_old_arr*/, - MultiArray4 const& /*state_new_arr*/, + const MultiFab& /*state_old*/, + const MultiFab& /*state_new*/, Vector>& /*a_extSource*/, const GeometryData& /*geomdata*/, - ProbParm const& /*prob_parm*/) + const ProbParm& /*prob_parm*/) { /* Notes: @@ -42,13 +42,15 @@ problem_modify_ext_sources( // Example: Exponential decay ode quantity auto ext_source_arr = a_extSource[lev]->arrays(); + auto const& state_old_arr = state_old.const_arrays(); + ParallelFor( *a_extSource[lev], [=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept { - for (int n = 0; n < NUM_ODE; n++){ + for (int n = 0; n < NUM_ODE; n++) { Real B_n = state_old_arr[box_no](i, j, k, FIRSTODE + n); Real src_strength = -1.0 * pow(10.0,n+1); - ext_source_arr[box_no](i, j, k, FIRSTODE + n) += src_strength * B_n; + ext_source_arr[box_no](i, j, k, FIRSTODE + n) += src; } }); Gpu::streamSynchronize(); From 897bb41c5aa7e5839aa96e4a43552a0577381d94 Mon Sep 17 00:00:00 2001 From: Dave Montgomery Date: Mon, 21 Oct 2024 16:38:56 -0600 Subject: [PATCH 44/46] Update PeleLMeX_ProblemSpecificFunctions.cpp --- Source/PeleLMeX_ProblemSpecificFunctions.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/PeleLMeX_ProblemSpecificFunctions.cpp b/Source/PeleLMeX_ProblemSpecificFunctions.cpp index e228a088a..8c6cd905d 100644 --- a/Source/PeleLMeX_ProblemSpecificFunctions.cpp +++ b/Source/PeleLMeX_ProblemSpecificFunctions.cpp @@ -50,9 +50,9 @@ problem_modify_ext_sources( for (int n = 0; n < NUM_ODE; n++) { Real B_n = state_old_arr[box_no](i, j, k, FIRSTODE + n); Real src_strength = -1.0 * pow(10.0,n+1); - ext_source_arr[box_no](i, j, k, FIRSTODE + n) += src; + ext_source_arr[box_no](i, j, k, FIRSTODE + n) += src_strength * B_n; } }); Gpu::streamSynchronize(); */ -} \ No newline at end of file +} From 3452ff59f0c5e1647d0682e53e1da3eb8c0599a2 Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Tue, 22 Oct 2024 11:57:44 -0600 Subject: [PATCH 45/46] Removed lev argument in problem_modify_ext_sources --- .../PeleLMeX_ProblemSpecificFunctions.cpp | 14 ++++++-------- Source/PeleLMeX_Forces.cpp | 5 +++-- Source/PeleLMeX_ProblemSpecificFunctions.H | 3 +-- Source/PeleLMeX_ProblemSpecificFunctions.cpp | 11 +++++------ 4 files changed, 15 insertions(+), 18 deletions(-) diff --git a/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp b/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp index d7669544e..9e3dce025 100644 --- a/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp +++ b/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp @@ -27,30 +27,28 @@ void problem_modify_ext_sources( Real /*time*/, Real /*dt*/, - int lev, const MultiFab& state_old, const MultiFab& /*state_new*/, - Vector>& a_extSource, + std::unique_ptr& ext_src, const GeometryData& /*geomdata*/, const ProbParm& prob_parm) { /* Notes: - 1) a_extSource contains sources from velocity forcing coming in. - This function should add to rather than overwrite a_extSource. + 1) ext_src contains sources from velocity forcing coming in. + This function should add to rather than overwrite ext_src. 2) Requires "peleLM.user_defined_ext_sources = true" in input file */ - auto ext_source_arr = a_extSource[lev]->arrays(); + auto ext_src_arr = ext_src->arrays(); auto const& state_old_arr = state_old.const_arrays(); ParallelFor( - *a_extSource[lev], - [=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept { + *ext_src, [=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept { for (int n = 0; n < NUM_ODE; n++) { Real B_n = state_old_arr[box_no](i, j, k, FIRSTODE + n); Real src = prob_parm.ode_srcstrength * pow(10.0, n + 1) * B_n; - ext_source_arr[box_no](i, j, k, FIRSTODE + n) += src; + ext_src_arr[box_no](i, j, k, FIRSTODE + n) += src; } }); Gpu::streamSynchronize(); diff --git a/Source/PeleLMeX_Forces.cpp b/Source/PeleLMeX_Forces.cpp index 164fac2ff..d51708213 100644 --- a/Source/PeleLMeX_Forces.cpp +++ b/Source/PeleLMeX_Forces.cpp @@ -273,9 +273,10 @@ PeleLM::getExternalSources( for (int lev = 0; lev <= finest_level; lev++) { auto* ldata_p_old = getLevelDataPtr(lev, a_timestamp_old); auto* ldata_p_new = getLevelDataPtr(lev, a_timestamp_new); + auto& ext_src = m_extSource[lev]; problem_modify_ext_sources( - getTime(lev, a_timestamp_new), m_dt, lev, ldata_p_old->state, - ldata_p_new->state, m_extSource, geom[lev].data(), *prob_parm_d); + getTime(lev, a_timestamp_new), m_dt, ldata_p_old->state, + ldata_p_new->state, ext_src, geom[lev].data(), *prob_parm_d); } } } \ No newline at end of file diff --git a/Source/PeleLMeX_ProblemSpecificFunctions.H b/Source/PeleLMeX_ProblemSpecificFunctions.H index 534c19344..e504e23e7 100644 --- a/Source/PeleLMeX_ProblemSpecificFunctions.H +++ b/Source/PeleLMeX_ProblemSpecificFunctions.H @@ -10,10 +10,9 @@ void set_ode_names(Vector& a_ode_names); void problem_modify_ext_sources( Real time, Real dt, - int lev, const MultiFab& state_old, const MultiFab& state_new, - Vector>& a_extSource, + std::unique_ptr& ext_src, const GeometryData& geomdata, const ProbParm& prob_parm); #endif \ No newline at end of file diff --git a/Source/PeleLMeX_ProblemSpecificFunctions.cpp b/Source/PeleLMeX_ProblemSpecificFunctions.cpp index 8c6cd905d..82b3bf8f2 100644 --- a/Source/PeleLMeX_ProblemSpecificFunctions.cpp +++ b/Source/PeleLMeX_ProblemSpecificFunctions.cpp @@ -27,10 +27,9 @@ void problem_modify_ext_sources( Real /*time*/, Real /*dt*/, - int /*lev*/, const MultiFab& /*state_old*/, const MultiFab& /*state_new*/, - Vector>& /*a_extSource*/, + std::unique_ptr& /*ext_src*/, const GeometryData& /*geomdata*/, const ProbParm& /*prob_parm*/) { @@ -41,16 +40,16 @@ problem_modify_ext_sources( 2) Requires peleLM.user_defined_ext_sources = true in input file // Example: Exponential decay ode quantity - auto ext_source_arr = a_extSource[lev]->arrays(); + auto ext_src_arr = ext_src->arrays(); auto const& state_old_arr = state_old.const_arrays(); ParallelFor( - *a_extSource[lev], + *ext_src, [=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept { for (int n = 0; n < NUM_ODE; n++) { Real B_n = state_old_arr[box_no](i, j, k, FIRSTODE + n); - Real src_strength = -1.0 * pow(10.0,n+1); - ext_source_arr[box_no](i, j, k, FIRSTODE + n) += src_strength * B_n; + Real src = -1.0 * pow(10.0, n + 1) * B_n; + ext_src_arr[box_no](i, j, k, FIRSTODE + n) += src; } }); Gpu::streamSynchronize(); From 27c4bed466739322ebc4472d6211952cba9fe2c6 Mon Sep 17 00:00:00 2001 From: dmontgomeryNREL Date: Tue, 22 Oct 2024 13:08:24 -0600 Subject: [PATCH 46/46] Update comment --- Source/PeleLMeX_ProblemSpecificFunctions.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/PeleLMeX_ProblemSpecificFunctions.cpp b/Source/PeleLMeX_ProblemSpecificFunctions.cpp index 82b3bf8f2..03477e877 100644 --- a/Source/PeleLMeX_ProblemSpecificFunctions.cpp +++ b/Source/PeleLMeX_ProblemSpecificFunctions.cpp @@ -35,8 +35,8 @@ problem_modify_ext_sources( { /* Notes: - 1) a_extSource contains sources from velocity forcing coming in. - This function should add to rather than overwrite a_extSource. + 1) ext_src contains sources from velocity forcing coming in. + This function should add to rather than overwrite ext_src. 2) Requires peleLM.user_defined_ext_sources = true in input file // Example: Exponential decay ode quantity