From 9d58cd77d408a47f7dd2ca591b20ddc01ab9e9ec Mon Sep 17 00:00:00 2001
From: Bhargav Sriram Siddani <siddanibhargav@gmail.com>
Date: Mon, 15 Jul 2024 22:12:03 -0700
Subject: [PATCH] Two-fluid effective viscosity: arithmetic & harmonic options

---
 src/incflo.H                     |  1 +
 src/rheology/incflo_rheology.cpp | 11 +++++++++--
 src/setup/init.cpp               |  1 +
 3 files changed, 11 insertions(+), 2 deletions(-)

diff --git a/src/incflo.H b/src/incflo.H
index cc9b3e9b..ff941330 100644
--- a/src/incflo.H
+++ b/src/incflo.H
@@ -549,6 +549,7 @@ private:
     amrex::Real m_eta_max = amrex::Real(1.0e+20);
     // Adding properties for a second fluid
     bool m_two_fluid = false;
+    bool m_two_fluid_eta_harmonic = false;
     amrex::Real m_ro_0_second = 1.0;
     FluidModel m_fluid_model_second;
     amrex::Real m_mu_second = 1.0;
diff --git a/src/rheology/incflo_rheology.cpp b/src/rheology/incflo_rheology.cpp
index 350b1ebd..369bde73 100644
--- a/src/rheology/incflo_rheology.cpp
+++ b/src/rheology/incflo_rheology.cpp
@@ -506,13 +506,20 @@ void incflo::compute_nodal_viscosity_at_level (int /*lev*/,
            Array4<Real const> const& eta_arr_second = vel_eta_second.const_array(mfi);
            Array4<Real> const& eta_arr = vel_eta->array(mfi);
            const Real min_conc_scnd = m_min_conc_second;
+           const bool eta_harmonic = m_two_fluid_eta_harmonic;
            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
            {
               if (conc_second_arr(i,j,k) > min_conc_scnd) {
                 // Using weighted harmonic mean for vel_eta
-                eta_arr(i,j,k) = ((Real(1.0)-conc_second_arr(i,j,k))/eta_arr(i,j,k))
+                if (eta_harmonic) {
+                    eta_arr(i,j,k) = ((Real(1.0)-conc_second_arr(i,j,k))/eta_arr(i,j,k))
                                   + (conc_second_arr(i,j,k)/eta_arr_second(i,j,k));
-                eta_arr(i,j,k) = Real(1.0)/eta_arr(i,j,k);
+                    eta_arr(i,j,k) = Real(1.0)/eta_arr(i,j,k);
+                }
+                else {
+                    eta_arr(i,j,k) = (Real(1.0)-conc_second_arr(i,j,k))*eta_arr(i,j,k) +
+                                     conc_second_arr(i,j,k)*eta_arr_second(i,j,k);
+                }
               }
            });
        }
diff --git a/src/setup/init.cpp b/src/setup/init.cpp
index 6324e276..513385dd 100644
--- a/src/setup/init.cpp
+++ b/src/setup/init.cpp
@@ -142,6 +142,7 @@ void incflo::ReadParameters ()
 
         // Check if there is a second fluid
         pp.query("two_fluid", m_two_fluid);
+        pp.query("two_fluid_eta_harmonic",m_two_fluid_eta_harmonic);
 
         // Density (if constant)
         pp.query("ro_0", m_ro_0);