Remove advection from auxiliary quantities B_k (#839)

Co-authored-by: Bruce Perry <[email protected]>

CMake/BuildPeleExe.cmake

@@ -22,10 +22,10 @@ function(build_pele_exe pele_exe_name pele_physics_lib_name)
   target_include_directories(${pele_exe_name} PRIVATE ${CMAKE_SOURCE_DIR}/Source/Params/param_includes)
   #Adv and Aux Variables
   if (PELEC_NUM_ADV GREATER 0)
-    target_compile_definitions(${pelec_exe_name} PRIVATE NUM_ADV=${PELEC_NUM_ADV})
+    target_compile_definitions(${pele_exe_name} PRIVATE NUM_ADV=${PELEC_NUM_ADV})
   endif()
   if (PELEC_NUM_AUX GREATER 0)
-    target_compile_definitions(${pelec_exe_name} PRIVATE NUM_AUX=${PELEC_NUM_AUX})
+    target_compile_definitions(${pele_exe_name} PRIVATE NUM_AUX=${PELEC_NUM_AUX})
   endif()
 
   target_sources(${pele_exe_name}

Docs/sphinx/Equations.rst

@@ -13,7 +13,7 @@ Equations
 Conservative system
 -------------------
 
-PeleC advances the following set of fully compressible equations for the conserved state vector: :math:`\mathbf{U} = (\rho, \rho \mathbf{u}, \rho E, \rho Y_k, \rho A_k, \rho B_k):`
+PeleC advances the following set of fully compressible equations for the conserved state vector: :math:`\mathbf{U} = (\rho, \rho \mathbf{u}, \rho E, \rho Y_k, \rho A_k, B_k):`
 
 .. math::
  
@@ -28,7 +28,7 @@ PeleC advances the following set of fully compressible equations for the conserv
 
   \frac{\partial (\rho A_k)}{\partial t} &=& - \nabla \cdot (\rho \mathbf{u} A_k) + S_{{\rm ext},\rho A_k},
 
-  \frac{\partial (\rho B_k)}{\partial t} &=& - \nabla \cdot (\rho \mathbf{u} B_k) + S_{{\rm ext},\rho B_k}.
+  \frac{\partial B_k}{\partial t} &=&  S_{{\rm ext}, B_k}.
 
 
 Here :math:`\rho, \mathbf{u}, T`, and :math:`p` are the density, velocity,
@@ -38,7 +38,7 @@ and chemical energy (species heats of formation) and is conserved across chemica
 :math:`Y_k` is the mass fraction of the :math:`k^{\rm th}` species,
 with associated production rate, :math:`\dot\omega_k`.  Here :math:`\mathbf{g}` is the gravitational vector, and
 :math:`S_{{\rm ext},\rho}, \mathbf{S}_{{\rm ext},\rho\mathbf{u}}`, etc., are user-specified
-source terms.  :math:`A_k` is an advected quantity, i.e., a tracer.  Also
+source terms.  :math:`A_k` is an advected quantity, i.e., a tracer. :math:`B_k` is spatially stationary quantity.  Also
 :math:`\boldsymbol{\mathcal{F}}_{m}, \mathbf{\Pi}`, and :math:`\boldsymbol{\mathcal{Q}}` are
 the diffusive transport fluxes for species, momentum and heat.  Note that the internal
 energy for species :math:`k` includes its heat of formation (and can therefore take on negative and
@@ -95,7 +95,7 @@ The inviscid equations for primitive variables namely density, velocity, and pre
   && \quad\qquad\qquad\qquad+\ S_{{\rm ext},\rho E} - \mathbf{u}\cdot\left(\mathbf{S}_{{\rm ext},\rho\mathbf{u}} - \frac{\mathbf{u}}{2}S_{{\rm ext},\rho}\right)\Biggr] 
   
 
-The advected quantities appear as:
+The advected and auxiliary quantities appear as:
 
 .. math::
   
@@ -105,8 +105,7 @@ The advected quantities appear as:
   \frac{\partial A_k}{\partial t} &=& -\mathbf{u}\cdot\nabla A_k + \frac{1}{\rho}
                                      ( S_{{\rm ext},\rho A_k} - A_k S_{{\rm ext},\rho} ),
 
-  \frac{\partial B_k}{\partial t} &=& -\mathbf{u}\cdot\nabla B_k + \frac{1}{\rho} 
-                                     ( S_{{\rm ext},\rho B_k}  - B_k S_{{\rm ext},\rho} ).
+  \frac{\partial B_k}{\partial t} &=&   S_{{\rm ext}, B_k}.
   
 
 
@@ -140,7 +139,7 @@ accounted for in the characteristic integration in the PPM algorithm.  The sourc
     \nabla\cdot k_{\rm th} \nabla T + S_{{\rm ext},\rho E} \\
     \frac{1}{\rho}S_{{\rm ext},\rho Y_k} \\
     \frac{1}{\rho}S_{{\rm ext},\rho A_k} \\
-    \frac{1}{\rho}S_{{\rm ext},\rho B_k}
+    S_{{\rm ext},B_k}
     \end{array}\right)^n,
     
 
@@ -153,7 +152,7 @@ accounted for in the characteristic integration in the PPM algorithm.  The sourc
     S_{\rho E} \\
     S_{\rho Y_k} \\
     S_{\rho A_k} \\
-    S_{\rho B_k}
+    S_{ B_k}
     \end{array}\right)^n
     =
     \left(\begin{array}{c}
@@ -162,5 +161,5 @@ accounted for in the characteristic integration in the PPM algorithm.  The sourc
     \rho \mathbf{u} \cdot \mathbf{g} + \nabla\cdot k_{\rm th} \nabla T + S_{{\rm ext},\rho E} \\
     S_{{\rm ext},\rho Y_k} \\
     S_{{\rm ext},\rho A_k} \\
-    S_{{\rm ext},\rho B_k}
+    S_{{\rm ext}, B_k}
     \end{array}\right)^n.

Docs/sphinx/ebverification/FlowPastCylinder/README.rst

@@ -16,7 +16,7 @@ where :math:`\rho` is the density and :math:`\mu` is the dynamic viscosity.
 
 Domain definition
 ##################################
-The cylinder with radius :math:`r=5` is centered at (0,0) on the upstream side of the domain as seen below. 
+The cylinder with radius :math:`r=0.5` [cm] is centered at (0,0) on the upstream side of the domain as seen below. 
 
 .. image:: /ebverification/FlowPastCylinder/domain.png
    :width: 450pt

Exec/RegTests/AuxQuantities/CMakeLists.txt

@@ -0,0 +1,9 @@
+set(PELE_PHYSICS_EOS_MODEL GammaLaw)
+set(PELE_PHYSICS_CHEMISTRY_MODEL Null)
+set(PELE_PHYSICS_TRANSPORT_MODEL Constant)
+set(PELE_PHYSICS_ENABLE_SOOT OFF)
+set(PELE_PHYSICS_ENABLE_SPRAY OFF)
+set(PELE_PHYSICS_SPRAY_FUEL_NUM 0)
+set(PELEC_NUM_ADV 2)
+set(PELEC_NUM_AUX 2)
+include(BuildExeAndLib)

Exec/RegTests/AuxQuantities/GNUmakefile

@@ -0,0 +1,40 @@
+# AMReX
+DIM = 2
+COMP = gnu
+PRECISION = DOUBLE
+
+# 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
+
+# PeleC
+PELE_CVODE_FORCE_YCORDER = FALSE
+PELE_USE_MAGMA = FALSE
+PELE_COMPILE_AJACOBIAN = FALSE
+Eos_Model := GammaLaw
+Chemistry_Model := Null
+Transport_Model := Constant
+PELEC_NUM_AUX = 2
+PELEC_NUM_ADV = 2
+
+# GNU Make
+Bpack := ./Make.package
+Blocs := .
+PELE_HOME := ../../..
+include $(PELE_HOME)/Exec/Make.PeleC

Exec/RegTests/AuxQuantities/Make.package

@@ -0,0 +1,3 @@
+CEXE_headers += prob.H
+CEXE_headers += prob_parm.H
+CEXE_sources += prob.cpp

Exec/RegTests/AuxQuantities/README.md

@@ -0,0 +1,14 @@
+# Advected and auxiliary quantities
+
+This demonstrates the use of the advected and auxiliary quantities described in the PeleC model [equations](https://amrex-combustion.github.io/PeleC/Equations.html). The case introduces two advected and two auxiliary quantities into the domain.  The auxiliary quantities experience simple exponential decay with a source term of $S_{ext,B_k} = -30 B_k$.
+
+
+## Short case description
+
+|                    | description                                         |
+|:-------------------|:----------------------------------------------------|
+| Problem definition | advected and auxiliaru quantities                   |
+| EB geometry        | embedded cylinder (optional)                        |
+| EOS                | GammaLaw                                            |
+| Multi-level        | yes                                                 |
+| Metric             | numerical solutions for $B_k$ match theory          |

Exec/RegTests/AuxQuantities/auxquantities.inp

@@ -0,0 +1,83 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+#max_step = 100
+stop_time = 0.1
+
+# PROBLEM SIZE & GEOMETRY
+geometry.is_periodic =  0  0  0
+geometry.coord_sys   =  0       # 0 => cart
+geometry.prob_lo     =  -2.   -2.  -2.
+geometry.prob_hi     =  10.  2.   2.
+amr.n_cell           =  96 32 8
+
+# >>>>>>>>>>>>>  BC KEYWORDS <<<<<<<<<<<<<<<<<<<<<<
+# Interior, UserBC, Symmetry, SlipWall, NoSlipWall
+# >>>>>>>>>>>>>  BC KEYWORDS <<<<<<<<<<<<<<<<<<<<<<
+
+pelec.lo_bc       =  "Hard" "SlipWall" "FOExtrap"
+pelec.hi_bc       =  "Hard" "SlipWall" "FOExtrap"
+
+# Problem setup
+prob.p = 1013250.
+prob.T = 100.
+prob.Re = 100.
+prob.Pr = 0.7
+prob.U = 100. 
+prob.aux_xy_lo = -1.
+prob.aux_length = 1.
+prob.aux_height = 2.
+prob.aux_srcstrength = -30.
+
+# Add user specified source term(s)
+pelec.add_ext_src = 1
+
+# WHICH PHYSICS
+pelec.do_hydro = 1
+pelec.do_mol = 0
+pelec.do_react = 0
+pelec.allow_negative_energy = 0
+pelec.diffuse_temp = 1
+pelec.diffuse_vel  = 1
+pelec.diffuse_spec = 0
+pelec.diffuse_enth = 1
+
+# TIME STEP CONTROL
+pelec.dt_cutoff      = 5.e-20  # level 0 timestep below which we halt
+pelec.cfl            = 0.7     # cfl number for hyperbolic system
+pelec.init_shrink    = 0.1     # scale back initial timestep
+pelec.change_max     = 1.05    # maximum increase in dt over successive steps
+
+# DIAGNOSTICS & VERBOSITY
+pelec.sum_interval   = 1       # timesteps between computing mass
+pelec.v              = 1       # verbosity in PeleC cpp files
+amr.v                = 1       # verbosity in Amr.cpp
+
+# REFINEMENT / REGRIDDING
+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.blocking_factor = 8        # block factor in grid generation
+amr.max_grid_size   = 64
+
+# TAGGING
+tagging.refinement_indicators = adv_0 adv_1
+tagging.adv_0.field_name = rho_A_0
+tagging.adv_0.adjacent_difference_greater = 0.01
+tagging.adv_0.max_level = 1
+tagging.adv_1.field_name = rho_A_1
+tagging.adv_1.adjacent_difference_greater = 0.02
+tagging.adv_1.max_level = 1
+
+# CHECKPOINT FILES
+amr.checkpoint_files_output = 0
+amr.check_file      = chk       # root name of checkpoint file
+amr.check_int       = 10000     # number of timesteps between checkpoints
+
+# PLOTFILES
+amr.plot_files_output = 1
+amr.plot_file       = plt
+amr.plot_int        = 1000
+amr.derive_plot_vars= x_velocity y_velocity z_velocity
+
+# EB
+ebd.boundary_grad_stencil_type = 0
+eb2.geom_type = "all_regular"

Exec/RegTests/AuxQuantities/auxquantities_eb.inp

@@ -0,0 +1,91 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+#max_step = 1
+stop_time = 0.1
+
+# PROBLEM SIZE & GEOMETRY
+geometry.is_periodic =  0  0  0
+geometry.coord_sys   =  0       # 0 => cart
+geometry.prob_lo     =  -2.   -2.  -2.
+geometry.prob_hi     =  10.  2.   2.
+amr.n_cell           =  96 32 8
+
+# >>>>>>>>>>>>>  BC KEYWORDS <<<<<<<<<<<<<<<<<<<<<<
+# Interior, UserBC, Symmetry, SlipWall, NoSlipWall
+# >>>>>>>>>>>>>  BC KEYWORDS <<<<<<<<<<<<<<<<<<<<<<
+
+pelec.lo_bc       =  "Hard" "SlipWall" "FOExtrap"
+pelec.hi_bc       =  "Hard" "SlipWall" "FOExtrap"
+
+# Problem setup
+prob.p = 1013250.
+prob.T = 100.
+prob.Re = 100.
+prob.Pr = 0.7
+prob.U = 100. 
+prob.aux_xy_lo = -1.
+prob.aux_length = 1.
+prob.aux_height = 2.
+prob.aux_srcstrength = -30.
+
+# Add user specified source term(s)
+pelec.add_ext_src = 1
+
+# WHICH PHYSICS
+pelec.do_hydro = 1
+pelec.do_mol = 0
+pelec.do_react = 0
+pelec.allow_negative_energy = 0
+pelec.diffuse_temp = 1
+pelec.diffuse_vel  = 1
+pelec.diffuse_spec = 0
+pelec.diffuse_enth = 1
+
+# TIME STEP CONTROL
+pelec.dt_cutoff      = 5.e-20  # level 0 timestep below which we halt
+pelec.cfl            = 0.7     # cfl number for hyperbolic system
+pelec.init_shrink    = 0.1     # scale back initial timestep
+pelec.change_max     = 1.05    # maximum increase in dt over successive steps
+
+# DIAGNOSTICS & VERBOSITY
+pelec.sum_interval   = 1       # timesteps between computing mass
+pelec.v              = 1       # verbosity in PeleC cpp files
+amr.v                = 1       # verbosity in Amr.cpp
+
+# REFINEMENT / REGRIDDING
+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.blocking_factor = 8       # block factor in grid generation
+amr.max_grid_size   = 64
+
+# TAGGING
+tagging.refinement_indicators = adv_0 adv_1
+tagging.adv_0.field_name = rho_A_0
+tagging.adv_0.adjacent_difference_greater = 0.01
+tagging.adv_0.max_level = 1
+tagging.adv_1.field_name = rho_A_1
+tagging.adv_1.adjacent_difference_greater = 0.02
+tagging.adv_1.max_level = 1
+
+# CHECKPOINT FILES
+amr.checkpoint_files_output = 0
+amr.check_file      = chk      # root name of checkpoint file
+amr.check_int       = 10000    # number of timesteps between checkpoints
+
+# PLOTFILES
+amr.plot_files_output = 1
+amr.plot_file       = plt
+amr.plot_int        = 1000
+amr.derive_plot_vars= x_velocity y_velocity z_velocity vfrac
+
+# EB
+ebd.boundary_grad_stencil_type = 0
+
+# for 2D (need a sphere):
+eb2.geom_type = "sphere"
+eb2.sphere_radius = 0.5
+eb2.sphere_center = 3 0 0
+eb2.sphere_has_fluid_inside = 0
+eb2.small_volfrac = 1.0e-4
+pelec.eb_boundary_T = 100.
+pelec.eb_isothermal = 1