run isothermal Soret channel as a test

Docs/sphinx/manual/LMeXControls.rst

@@ -271,8 +271,8 @@ PeleLMeX algorithm
     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
 ------------------------------
 
@@ -291,6 +291,13 @@ Transport coefficients and LES
     peleLM.les_cs_sigma = 1.35             # [OPT, DEF=1.35] If using Sigma LES model, provides model coefficient
     peleLM.les_v = 0                       # [OPT, DEF=0] Verbosity level for LES model
     peleLM.plot_les = 0                    # [OPT, DEF=0] If doing LES, whether to plot the turbulent viscosity
+    transport.use_soret = 0                # [OPT, DEF=0] Compute diffusion including the Soret effect (note, this option is inherited from PelePhysics)
+
+.. note::
+   When using the Soret effect, boundary condition corrections are needed at isothermal boundaries,
+   which are not fully supported. Currently a correction for all terms except the wbar term
+   is applied at isothermal domain boundaries (this is likely sufficient), while no corrections are applied
+   at isothermal embedded boundaries (so use caution for isothermal EBs with Soret diffusion active).
 
 Chemistry integrator
 --------------------

Exec/RegTests/CMakeLists.txt

@@ -21,6 +21,7 @@ else()
   endif()
   if(PELE_DIM EQUAL 2)
     add_subdirectory(TripleFlame)
+    add_subdirectory(IsothermalSoretChannel)
   endif()
   if(PELE_DIM EQUAL 3)
     add_subdirectory(TurbInflow)

Exec/RegTests/IsothermalSoretChannel/input.inp → Exec/RegTests/IsothermalSoretChannel/isothermal-soret.inp

@@ -25,7 +25,7 @@ prob.P_mean = 101325.0
 prob.Vin = 10.0
 prob.Tlow = 400.0
 prob.Thigh = 800.0
-prob.phi = 0.5
+prob.phi = 1.0
 #---------------------- Transport --------------------------------
 transport.use_soret = 1
 
@@ -46,7 +46,7 @@ peleLM.do_species_balance = 1             # Compute species balance
 amr.max_step = 2000                       # Maximum number of time steps
 amr.stop_time = 4.00                      # final simulation physical time [s]
 amr.max_wall_time = 1                     # Maximum simulation run time [h]
-amr.cfl = 0.5                             # CFL number for hyperbolic system 
+amr.cfl = 0.5                             # CFL number for hyperbolic system
 amr.dt_shrink = 0.01                     # Scale back initial timestep
 amr.dt_change_max = 1.1                   # Maximum dt increase btw successive steps
 
@@ -86,6 +86,6 @@ amr.gradT.adjacent_difference_greater = 200
 amr.gradT.field_name = temp
 
 #---------------------- Debug/HPC CONTROL-------------------------
-amrex.fpe_trap_invalid = 1
-amrex.fpe_trap_zero = 1
-amrex.fpe_trap_overflow = 1
+#amrex.fpe_trap_invalid = 1
+#amrex.fpe_trap_zero = 1
+#amrex.fpe_trap_overflow = 1

Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.H

@@ -59,11 +59,8 @@ pelelmex_initdata(
   const amrex::Real* prob_hi = geomdata.ProbHi();
   const amrex::Real* dx = geomdata.CellSize();
 
-  const amrex::Real x = prob_lo[0] + (i + 0.5) * dx[0];
   const amrex::Real y = prob_lo[1] + (j + 0.5) * dx[1];
   const amrex::Real y_low = prob_lo[1];
-
-  const amrex::Real Lx = prob_hi[0] - prob_lo[0];
   const amrex::Real Ly = prob_hi[1] - prob_lo[1];
 
   auto eos = pele::physics::PhysicsType::eos();
@@ -98,19 +95,16 @@ AMREX_GPU_HOST_DEVICE
 AMREX_FORCE_INLINE
 void
 bcnormal(
-  const amrex::Real x[AMREX_SPACEDIM],
+  const amrex::Real* /*x[AMREX_SPACEDIM]*/,
   const int /*m_nAux*/,
   amrex::Real s_ext[NVAR],
   const int /*idir*/,
   const int sgn,
-  const amrex::Real time,
+  const amrex::Real /*time*/,
   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};
-
   if (sgn == 1) {
     s_ext[TEMP] = prob_parm.Tlow;
   } else if (sgn == -1) {

Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.cpp

@@ -10,6 +10,7 @@ PeleLM::readProbParm()
   pp.query("Vin", PeleLM::prob_parm->Vin);
   pp.query("Thigh", PeleLM::prob_parm->Thigh);
   pp.query("Tlow", PeleLM::prob_parm->Tlow);
+  pp.query("phi", PeleLM::prob_parm->phi);
 
   PeleLM::prob_parm->bathID = N2_ID;
   PeleLM::prob_parm->oxidID = O2_ID;

Exec/RegTests/IsothermalSoretChannel/test.py

@@ -0,0 +1,31 @@
+import os
+import numpy.testing as npt
+import pandas as pd
+import unittest
+
+class SpeciesBalTestCase(unittest.TestCase):
+    """Test species balance with isothermal walls and soret"""
+
+    def test_composition(self):
+        """Verify species conservation"""
+
+        # Load the data
+        fdir = os.path.abspath(".")
+        fname = os.path.join(fdir, "temporals/tempSpecies")
+        df = pd.read_csv(fname)
+        print(df)
+        for col in df.columns:
+            if col.startswith('rhoYnew'):
+                init_value = df[col][0]
+                if init_value != 0.0:
+                    print('testing (relative)', col)
+                    npt.assert_allclose(df[col], init_value, rtol=1e-13)
+                else:
+                    print('testing (absolute)', col)
+                    npt.assert_allclose(df[col], init_value, atol=1e-13)
+            if col.startswith('netFlux'):
+                print('testing (absolute)', col)
+                npt.assert_allclose(df[col], 0.0, atol=1e-13)
+
+if __name__ == "__main__":
+    unittest.main()

Tests/CMakeLists.txt

@@ -236,6 +236,7 @@ if(NOT PELE_ENABLE_EB)
   endif()
   if(PELE_DIM EQUAL 2)
     add_test_r(tripleflame-${PELE_DIM}d TripleFlame)
+    add_test_rt(isothermal-soret IsothermalSoretChannel)
   endif()
   if(PELE_DIM EQUAL 3)
     add_test_r(hit-${PELE_DIM}d HITDecay)