AMReX-Combustion · baperry2 · Oct 22, 2024 · Sep 30, 2024 · Sep 30, 2024 · Oct 1, 2024
diff --git 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
@@ -46,10 +45,12 @@ 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_Projection.cpp
        ${SRC_DIR}/PeleLMeX_Reactions.cpp
        ${SRC_DIR}/PeleLMeX_Regrid.cpp
@@ -69,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
@@ -88,6 +109,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

diff --git 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
@@ -143,6 +143,12 @@ Thus, the terms containing :math:`S_{\text{ext},\rho}` cancel and the divergence
     \nabla \cdot \boldsymbol{u} = \frac{1}{\rho c_p T} \Big(-\nabla \cdot \boldsymbol{Q} + S_{\text{ext},\rho h} \Big) 
     +  \sum_m \bigg( \frac{W}{\rho W_m} -  \frac{h_m}{\rho c_p T}\bigg)\bigg( - \nabla \cdot \boldsymbol{\mathcal{F}}_m + \rho \dot \omega_m + S_{\text{ext},\rho Y_m} \bigg) \equiv S .
 
+In addition to the flow equations, `PeleLMeX` can also solve for a set of quantities that are neither advected nor diffused, satisfying:
+
+.. math::
+
+    \frac{\partial B_k}{\partial t} = S_{\text{ext},B_k}.
+
 Confined domain ambient pressure
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 

diff --git 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$(PELELM_NUM_ODE) -gt 0; echo $$?), 0)
+   DEFINES += -DNUM_ODE=$(PELELM_NUM_ODE)
+endif
+
 Bpack += $(foreach dir, $(LMdirs), $(PELE_HOME)/$(dir)/Make.package)
 Blocs += $(foreach dir, $(LMdirs), $(PELE_HOME)/$(dir))
 

diff --git a/Exec/RegTests/CMakeLists.txt b/Exec/RegTests/CMakeLists.txt
@@ -3,6 +3,7 @@ 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()

diff --git 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
@@ -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
@@ -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
@@ -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
@@ -0,0 +1,11 @@
+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)
+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
@@ -0,0 +1,38 @@
+# 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_sources += PeleLMeX_ProblemSpecificFunctions.cpp
+PELELM_NUM_ODE = 3
+
+# PelePhysics
+Chemistry_Model = air
+Eos_Model = Fuego
+Transport_Model = Constant
+
+PELE_HOME ?= ../../..
+include $(PELE_HOME)/Exec/Make.PeleLMeX
diff --git a/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp b/Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp
@@ -0,0 +1,55 @@
+#include <PeleLMeX.H>
+#include <PeleLMeX_K.H>
+#include <PeleLMeX_ProblemSpecificFunctions.H>
+
+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
+*/
+
+#if NUM_ODE > 0
+void
+set_ode_names(Vector<std::string>& a_ode_names)
+{
+  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*/,
+  const MultiFab& state_old,
+  const MultiFab& /*state_new*/,
+  std::unique_ptr<MultiFab>& ext_src,
+  const GeometryData& /*geomdata*/,
+  const ProbParm& prob_parm)
+{
+  /*
+  Notes:
+    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_src_arr = ext_src->arrays();
+  auto const& state_old_arr = state_old.const_arrays();
+
+  ParallelFor(
+    *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_src_arr[box_no](i, j, k, FIRSTODE + n) += src;
+      }
+    });
+  Gpu::streamSynchronize();
+}
diff --git a/Exec/RegTests/EB_ODEQty/README.md b/Exec/RegTests/EB_ODEQty/README.md
@@ -0,0 +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,\dots,$ NUM_ODE and $\gamma < 0$.
diff --git a/Exec/RegTests/EB_ODEQty/eb-odeqty-2d.inp b/Exec/RegTests/EB_ODEQty/eb-odeqty-2d.inp
@@ -0,0 +1,75 @@
+#----------------------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 = -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 = 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/eb-odeqty-3d.inp 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
Original file line number	Diff line number	Diff line change
		@@ -0,0 +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,\dots,$ NUM_ODE and $\gamma < 0$.