Soret isothermal wall fix

Exec/RegTests/IsothermalSoretChannel/CMakeLists.txt

@@ -0,0 +1,8 @@
+set(PELE_PHYSICS_EOS_MODEL Fuego)
+set(PELE_PHYSICS_CHEMISTRY_MODEL LiDryer)
+set(PELE_PHYSICS_TRANSPORT_MODEL Simple)
+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)

Exec/RegTests/IsothermalSoretChannel/GNUmakefile

@@ -0,0 +1,37 @@
+# AMReX
+DIM = 2
+COMP = gnu
+PRECISION = DOUBLE
+USE_EB = FALSE
+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
+USE_EFIELD = FALSE
+
+# PelePhysics
+Chemistry_Model = LiDryer
+Eos_Model = Fuego
+Transport_Model = Simple
+
+PELE_HOME ?= ../../..
+include $(PELE_HOME)/Exec/Make.PeleLMeX

Exec/RegTests/IsothermalSoretChannel/README.md

@@ -0,0 +1,3 @@
+## Soret Isothermal wall test
+
+2D case of periodic channel with isothermal walls and Soret effects enabled. Test to ensure zero species flux at walls.

Exec/RegTests/IsothermalSoretChannel/input.inp

@@ -0,0 +1,91 @@
+#---------------------- DOMAIN DEFINITION ------------------------
+geometry.is_periodic = 1 0                # For each dir, 0: non-perio, 1: periodic
+geometry.coord_sys   = 0                  # 0 => cart, 1 => RZ
+geometry.prob_lo     = 0.0 0.0 0.0        # x_lo y_lo (z_lo)
+geometry.prob_hi     = 0.1 0.0125 0.0      # x_hi y_hi (z_hi)
+
+#---------------------- BC FLAGS ---------------------------------
+# Interior, Inflow, Outflow, Symmetry,
+# SlipWallAdiab, NoSlipWallAdiab, SlipWallIsotherm, NoSlipWallIsotherm
+peleLM.lo_bc = Interior NoSlipWallIsotherm   # bc in x_lo y_lo (z_lo)
+peleLM.hi_bc = Interior NoSlipWallIsotherm   # bc in x_hi y_hi (z_hi)
+
+#---------------------- AMR CONTROL ------------------------------
+amr.n_cell          = 256 32                # Level 0 number of cells in each direction
+amr.max_level       = 1                   # maximum level number allowed
+amr.ref_ratio       = 2 2 2 2             # refinement ratio
+amr.regrid_int      = 2                   # how often to regrid
+amr.n_error_buf     = 1 1 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   = 32                  # max box size
+
+#---------------------- Problem ----------------------------------
+prob.P_mean = 101325.0
+prob.Vin = 10.0
+prob.Tlow = 400.0
+prob.Thigh = 800.0
+prob.phi = 0.5
+#---------------------- Transport --------------------------------
+transport.use_soret = 1
+
+#---------------------- PeleLM CONTROL ---------------------------
+peleLM.v = 1                              # PeleLMeX verbose
+peleLM.use_wbar = 1                       # Include Wbar term in species diffusion fluxes
+peleLM.sdc_iterMax = 2                    # Number of SDC iterations
+peleLM.num_init_iter = 3                  # Number of initial iterations
+
+#---------------------- Temporal CONTROL -------------------------
+peleLM.do_temporals = 1                   # Turn temporals ON/OFF
+peleLM.temporal_int = 10                  # Frequency of temporals
+peleLM.do_extremas = 1                    # Compute state extremas
+peleLM.do_mass_balance = 1                # Compute mass balance
+peleLM.do_species_balance = 1             # Compute species balance
+
+#---------------------- Time Stepping CONTROL --------------------
+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.dt_shrink = 0.01                     # Scale back initial timestep
+amr.dt_change_max = 1.1                   # Maximum dt increase btw successive steps
+
+#---------------------- IO CONTROL -------------------------------
+#amr.restart = chk01000                   # Restart checkpoint file
+amr.check_int = 2000                      # Frequency of checkpoint output
+amr.plot_int = 20                         # Frequency of pltfile output
+amr.derive_plot_vars = avg_pressure mag_vort mass_fractions mixture_fraction
+
+#---------------------- Derived CONTROLS -------------------------
+peleLM.fuel_name = H2
+peleLM.mixtureFraction.format = Cantera
+peleLM.mixtureFraction.type   = mass
+peleLM.mixtureFraction.oxidTank = O2:0.233 N2:0.767
+peleLM.mixtureFraction.fuelTank = H2:1.0
+
+#---------------------- Reactor CONTROL --------------------------
+peleLM.chem_integrator = "ReactorNull"
+#peleLM.use_typ_vals_chem = 1              # Use species/temp typical values in CVODE
+#ode.rtol = 1.0e-6                         # Relative tolerance of the chemical solve
+#ode.atol = 1.0e-5                         # Absolute tolerance factor applied on typical values
+#cvode.solve_type = denseAJ_direct         # CVODE Linear solve type (for Newton direction)
+#cvode.max_order  = 4                      # CVODE max BDF order.
+
+#---------------------- Linear solver CONTROL --------------------
+mac_proj.verbose = 0
+nodal_proj.verbose = 0
+
+#---------------------- Refinement CONTROL------------------------
+amr.refinement_indicators = highT gradT
+amr.highT.max_level     = 1
+amr.highT.value_greater = 600
+amr.highT.field_name    = temp
+
+amr.gradT.max_level = 2
+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

Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.H

@@ -0,0 +1,141 @@
+#ifndef PELELM_PROB_H
+#define PELELM_PROB_H
+
+#include "PeleLMeX_Index.H"
+#include "pelelmex_prob_parm.H"
+#include "PMFData.H"
+#include "PelePhysics.H"
+
+#include <AMReX_Geometry.H>
+#include <AMReX_FArrayBox.H>
+#include <AMReX_SPACE.H>
+
+AMREX_GPU_HOST_DEVICE
+AMREX_FORCE_INLINE
+void
+set_Y_from_Phi(ProbParm const& prob_parm, amrex::Real Y[])
+{
+  amrex::Real phi = prob_parm.phi;
+  auto eos = pele::physics::PhysicsType::eos();
+  amrex::Real Xt[NUM_SPECIES] = {0.0};
+  amrex::Real a = 0.0;
+  switch (prob_parm.fuelID) {
+#if defined(H2_ID)
+  case H2_ID:
+    a = 0.5;
+    break;
+#endif
+#if defined(CH4_ID)
+  case CH4_ID:
+    a = 2.0;
+    break;
+#endif
+  default:
+    amrex::Abort("fuelID not defined");
+  }
+
+  Xt[prob_parm.oxidID] = 1.0 / (1.0 + phi / a + 0.79 / 0.21);
+  Xt[prob_parm.fuelID] = phi * Xt[prob_parm.oxidID] / a;
+  Xt[prob_parm.bathID] = 1.0 - Xt[prob_parm.oxidID] - Xt[prob_parm.fuelID];
+
+  eos.X2Y(Xt, Y);
+}
+
+AMREX_GPU_DEVICE
+AMREX_FORCE_INLINE
+void
+pelelmex_initdata(
+  int i,
+  int j,
+  int k,
+  int /*is_incompressible*/,
+  amrex::Array4<amrex::Real> const& state,
+  amrex::Array4<amrex::Real> const& /*aux*/,
+  amrex::GeometryData const& geomdata,
+  ProbParm const& prob_parm,
+  pele::physics::PMF::PmfData::DataContainer const* /*pmf_data*/)
+{
+  const amrex::Real* prob_lo = geomdata.ProbLo();
+  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();
+  amrex::Real massfrac[NUM_SPECIES] = {0.0};
+
+  set_Y_from_Phi(prob_parm, massfrac);
+
+  const amrex::Real T_low = prob_parm.Tlow;
+  const amrex::Real T_high = prob_parm.Thigh;
+
+  state(i, j, k, TEMP) = (y - y_low) / (Ly) * (T_high - T_low) + T_low;
+
+  state(i, j, k, VELX) = prob_parm.Vin;
+  state(i, j, k, VELY) = 0.0;
+
+  amrex::Real rho_cgs, P_cgs;
+  P_cgs = prob_parm.P_mean * 10.0;
+
+  eos.PYT2R(P_cgs, massfrac, state(i, j, k, TEMP), rho_cgs);
+  state(i, j, k, DENSITY) = rho_cgs * 1.0e3; // CGS -> MKS conversion
+
+  eos.TY2H(state(i, j, k, TEMP), massfrac, state(i, j, k, RHOH));
+  state(i, j, k, RHOH) *=
+    1.0e-4 * state(i, j, k, DENSITY); // CGS -> MKS conversion
+
+  for (int n = 0; n < NUM_SPECIES; n++) {
+    state(i, j, k, FIRSTSPEC + n) = massfrac[n] * state(i, j, k, DENSITY);
+  }
+}
+
+AMREX_GPU_HOST_DEVICE
+AMREX_FORCE_INLINE
+void
+bcnormal(
+  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,
+  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) {
+    s_ext[TEMP] = prob_parm.Thigh;
+  }
+}
+
+AMREX_GPU_DEVICE
+AMREX_FORCE_INLINE
+void
+zero_visc(
+  int i,
+  int j,
+  int k,
+  amrex::Array4<amrex::Real> 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

Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.cpp

@@ -0,0 +1,33 @@
+#include <PeleLMeX.H>
+#include <AMReX_ParmParse.H>
+
+void
+PeleLM::readProbParm()
+{
+  amrex::ParmParse pp("prob");
+
+  pp.query("P_mean", PeleLM::prob_parm->P_mean);
+  pp.query("Vin", PeleLM::prob_parm->Vin);
+  pp.query("Thigh", PeleLM::prob_parm->Thigh);
+  pp.query("Tlow", PeleLM::prob_parm->Tlow);
+
+  PeleLM::prob_parm->bathID = N2_ID;
+  PeleLM::prob_parm->oxidID = O2_ID;
+  // get the fuel name from the peleLM.fuel_name declaration
+  amrex::ParmParse pp_pele("peleLM");
+  std::string fuelName = "";
+  pp_pele.get("fuel_name", fuelName);
+#if defined(H2_ID)
+  if (fuelName == "H2") {
+    PeleLM::prob_parm->fuelID = H2_ID;
+  } else
+#endif
+#if defined(CH4_ID)
+    if (fuelName == "CH4") {
+    PeleLM::prob_parm->fuelID = CH4_ID;
+  } else
+#endif
+  {
+    amrex::Abort("fuel_name not recognised!");
+  }
+}

Exec/RegTests/IsothermalSoretChannel/pelelmex_prob_parm.H

@@ -0,0 +1,21 @@
+#ifndef PELELM_PROB_PARM_H
+#define PELELM_PROB_PARM_H
+
+#include <AMReX_REAL.H>
+#include <PeleLMeX_FlowControllerData.H>
+
+using namespace amrex::literals;
+
+struct ProbParm
+{
+  amrex::Real P_mean = 101325.0_rt;
+  amrex::Real Vin = 10.0;
+  amrex::Real Thigh = 700.0;
+  amrex::Real Tlow = 300.0;
+  amrex::Real phi = 1.0;
+
+  int bathID{-1};
+  int fuelID{-1};
+  int oxidID{-1};
+};
+#endif

Source/Efield/PeleLMeX_EFIonDrift.cpp

@@ -41,10 +41,10 @@ PeleLM::ionDriftVelocity(std::unique_ptr<AdvanceAdvData>& advData)
   auto bcRecPhiV = fetchBCRecArray(PHIV, 1);
   getDiffusionOp()->computeGradient(
     GetVecOfArrOfPtrs(gphiVOld), {}, // don't need the laplacian out
-    GetVecOfConstPtrs(getPhiVVect(AmrOldTime)), bcRecPhiV[0], do_avgDown);
+    GetVecOfConstPtrs(getPhiVVect(AmrOldTime)), {}, bcRecPhiV[0], do_avgDown);
   getDiffusionOp()->computeGradient(
     GetVecOfArrOfPtrs(gphiVNew), {}, // don't need the laplacian out
-    GetVecOfConstPtrs(getPhiVVect(AmrNewTime)), bcRecPhiV[0], do_avgDown);
+    GetVecOfConstPtrs(getPhiVVect(AmrNewTime)), {}, bcRecPhiV[0], do_avgDown);
 
   //----------------------------------------------------------------
   // TODO : this assumes that all the ions are grouped together at th end ...

Source/Efield/PeleLMeX_EFNLSolve.cpp

@@ -53,7 +53,7 @@ PeleLM::implicitNonLinearSolve(
   auto bcRecPhiV = fetchBCRecArray(PHIV, 1);
   getDiffusionOp()->computeGradient(
     getNLgradPhiVVect(), {}, // don't need the laplacian out
-    GetVecOfConstPtrs(getPhiVVect(AmrOldTime)), bcRecPhiV[0], do_avgDown);
+    GetVecOfConstPtrs(getPhiVVect(AmrOldTime)), {}, bcRecPhiV[0], do_avgDown);
 
   // Stash away a copy of umac
   for (int lev = 0; lev <= finest_level; ++lev) {
@@ -438,7 +438,7 @@ PeleLM::nonLinearResidual(
   auto bcRecPhiV = fetchBCRecArray(PHIV, 1);
   getDiffusionOp()->computeGradient(
     GetVecOfArrOfPtrs(gradPhiVCur), GetVecOfPtrs(laplacian),
-    GetVecOfConstPtrs(phiV), bcRecPhiV[0], do_avgDown);
+    GetVecOfConstPtrs(phiV), {}, bcRecPhiV[0], do_avgDown);
 
   // Get nE diffusion term
   Vector<MultiFab> diffnE(finest_level + 1);

Source/Efield/PeleLMeX_EFPoisson.cpp

@@ -47,5 +47,5 @@ PeleLM::poissonSolveEF(const TimeStamp& a_time)
   // Solve for PhiV
   getDiffusionOp()->diffuse_scalar(
     GetVecOfPtrs(getPhiVVect(a_time)), 0, GetVecOfConstPtrs(rhsPoisson), 0, {},
-    0, {}, {}, {}, 0, bcRecPhiV, 1, 1, -eps0 * epsr);
+    0, {}, {}, {}, 0, bcRecPhiV, 1, 1, -eps0 * epsr, {});
 }