From b52b30cd811443844d47e5b530582bf442cc649a Mon Sep 17 00:00:00 2001
From: Thomas Howarth <t.howarth@itv.rwth-aachen.de>
Date: Fri, 15 Nov 2024 14:07:34 +0100
Subject: [PATCH] disable wbar

---
 Source/PeleLMeX_Diffusion.cpp | 21 +++++++++++----------
 Source/PeleLMeX_K.H           |  8 +++++---
 Source/PeleLMeX_Setup.cpp     |  8 +++++---
 3 files changed, 21 insertions(+), 16 deletions(-)

diff --git a/Source/PeleLMeX_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp
index 93e29134..5c43ced5 100644
--- a/Source/PeleLMeX_Diffusion.cpp
+++ b/Source/PeleLMeX_Diffusion.cpp
@@ -414,7 +414,7 @@ PeleLM::correctIsothermalBoundary(
       }
     }
   }
-  //TODO: wbar fluxes disabled for this case - boundary system becomes complex
+  // TODO: wbar fluxes disabled for this case - boundary system becomes complex
 }
 
 void
@@ -631,18 +631,18 @@ PeleLM::addWbarTerm(
 
       amrex::ParallelFor(
 
-        gbx,
-        [rho_arr, rhoY_arr, Wbar_arr, gradY_arr, Wbar_boundary_arr, domain,
-         have_boundary,
-         phys_bc = m_phys_bc,leosparm] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
-          getMwmixGivenRY(i, j, k, rho_arr, rhoY_arr, Wbar_arr,leosparm);
+        gbx, [rho_arr, rhoY_arr, Wbar_arr, gradY_arr, Wbar_boundary_arr, domain,
+              have_boundary, phys_bc = m_phys_bc,
+              leosparm] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+          getMwmixGivenRY(i, j, k, rho_arr, rhoY_arr, Wbar_arr, leosparm);
           if (have_boundary != 0) { // need to impose gradWbar on boundary for
             // computeGradient
             // for dirichlet boundaries, we'll overwrite inhomog neumann ones
-	    // NOTE: for now, this is skipped since wbar disabled for isothermal/soret
+            // NOTE: for now, this is skipped since wbar disabled for
+            // isothermal/soret
             Wbar_boundary_arr(i, j, k) = Wbar_arr(i, j, k);
             int idx[3] = {i, j, k};
-	    for (int idim = 0; idim < AMREX_SPACEDIM; idim++) {
+            for (int idim = 0; idim < AMREX_SPACEDIM; idim++) {
               const auto bc_lo = phys_bc.lo(idim);
               const auto bc_hi = phys_bc.hi(idim);
               bool on_lo = (bc_lo == BoundaryCondition::BCNoSlipWallIsotherm ||
@@ -651,9 +651,10 @@ PeleLM::addWbarTerm(
               bool on_hi = (bc_hi == BoundaryCondition::BCNoSlipWallIsotherm ||
                             bc_hi == BoundaryCondition::BCSlipWallIsotherm) &&
                            (idx[idim] > domain.bigEnd(idim));
-	      
+
               if (on_lo || on_hi) {
-		getGradMwmixGivengradYMwmix(i, j, k, gradY_arr, Wbar_arr, Wbar_boundary_arr,leosparm);
+                getGradMwmixGivengradYMwmix(
+                  i, j, k, gradY_arr, Wbar_arr, Wbar_boundary_arr, leosparm);
               }
             }
           }
diff --git a/Source/PeleLMeX_K.H b/Source/PeleLMeX_K.H
index b032bbb5..3ffbfa0a 100644
--- a/Source/PeleLMeX_K.H
+++ b/Source/PeleLMeX_K.H
@@ -571,16 +571,18 @@ getGradMwmixGivengradYMwmix(
   int k,
   amrex::Array4<const amrex::Real> const& gradY,
   amrex::Array4<const amrex::Real> const& Mwmix,
-  amrex::Array4<amrex::Real> const& gradMwmix) noexcept
+  amrex::Array4<amrex::Real> const& gradMwmix,
+  pele::physics::eos::EosParm<pele::physics::PhysicsType::eos_type> const*
+    eosparm) noexcept
 {
   using namespace amrex::literals;
 
-  auto eos = pele::physics::PhysicsType::eos();
+  auto eos = pele::physics::PhysicsType::eos(eosparm);
   amrex::Real imw[NUM_SPECIES] = {0.0_rt};
   eos.inv_molecular_weight(imw);
   gradMwmix(i, j, k) = 0.0;
   for (int n = 0; n < NUM_SPECIES; n++) {
-    // imw[n] *= 1000.0; // CGS -> MKS
+    imw[n] *= 1000.0; // CGS -> MKS
     gradMwmix(i, j, k) += gradY(i, j, k, n) * imw[n];
   }
   // gradWbar = -Wbar^2 * sum(grad Y_k / W_k)
diff --git a/Source/PeleLMeX_Setup.cpp b/Source/PeleLMeX_Setup.cpp
index 34ee0e6a..d4c99415 100644
--- a/Source/PeleLMeX_Setup.cpp
+++ b/Source/PeleLMeX_Setup.cpp
@@ -103,9 +103,10 @@ PeleLM::Setup()
             << "    Using mixture-averaged transport with Soret effects"
             << std::endl;
           if (m_soret_boundary_override != 0) {
-            amrex::Print() << "    Imposing inhomogeneous Neumann conditions "
-                              "for species on isothermal walls"
-                           << std::endl;
+            amrex::Print()
+              << "    Imposing inhomogeneous Neumann conditions "
+                 "for species on isothermal walls. WARNING: use_wbar disabled."
+              << std::endl;
           }
         }
       } else {
@@ -437,6 +438,7 @@ PeleLM::readParameters()
     }
     if (isothermal) {
       m_soret_boundary_override = 1;
+      m_use_wbar = 0;
 #if PELE_USE_EFIELD
       amrex::Abort("Isothermal walls with Soret incompatible with Efield");
 #endif