Also saving mu(I)-stress ratio when m_fluid_model_second is FluidMode…

…l::DataDrivenMPMD
siddanib · Nov 1, 2024 · 663f410 · 663f410
1 parent 819a575
commit 663f410
Show file tree

Hide file tree

Showing 4 changed files with 49 additions and 9 deletions.
diff --git a/src/incflo.H b/src/incflo.H
@@ -637,6 +637,9 @@ private:
     int m_plt_error_p      = 0;
     int m_plt_error_mac_p  = 0;
     int m_plt_inertial_num = 0;
+#ifdef USE_AMREX_MPMD
+    int m_plt_mu_I         = 0;
+#endif
 
 #ifdef INCFLO_USE_PARTICLES
     int m_plt_particle_count = 0;

diff --git a/src/rheology/incflo_rheology.cpp b/src/rheology/incflo_rheology.cpp
@@ -53,7 +53,7 @@ void incflo::compute_viscosity (Vector<MultiFab*> const& vel_eta,
                                 Real time, int nghost)
 {
 #ifdef USE_AMREX_MPMD
-    // Call to indicate this is not final data-transfer
+    // Call to indicate this is from incflo::compute_viscosity
     if (ParallelDescriptor::MyProc() == 0) {
         Vector<int> last_call;
         last_call.push_back(0);

diff --git a/src/setup/init.cpp b/src/setup/init.cpp
@@ -300,6 +300,9 @@ void incflo::ReadIOParameters()
     pp.query("plt_error_p",    m_plt_error_p );
     pp.query("plt_error_mac_p",m_plt_error_mac_p );
     pp.query("plt_inertial_num", m_plt_inertial_num);
+#ifdef USE_AMREX_MPMD
+    pp.query("plt_mu_I", m_plt_mu_I);
+#endif
 
 #ifdef INCFLO_USE_PARTICLES
     pp.query("plt_particle_count", m_plt_particle_count );

diff --git a/src/utilities/io.cpp b/src/utilities/io.cpp
@@ -420,11 +420,14 @@ void incflo::WritePlotFile()
 
     // Nodal Inertial Number in mu(I)
     if (m_nodal_vel_eta && m_plt_inertial_num
-        && (m_fluid_model == FluidModel::DataDrivenMPMD
-            || m_fluid_model == FluidModel::Rauter
-            || m_fluid_model_second == FluidModel::DataDrivenMPMD
+        && (m_fluid_model_second == FluidModel::DataDrivenMPMD
             || m_fluid_model_second == FluidModel::Rauter)) ++ncomp;
 
+#ifdef USE_AMREX_MPMD
+    if (m_nodal_vel_eta && m_plt_mu_I
+        && m_fluid_model_second == FluidModel::DataDrivenMPMD) ++ncomp;
+#endif
+
     Vector<MultiFab> mf(finest_level + 1);
     for (int lev = 0; lev <= finest_level; ++lev) {
         mf[lev].define(grids[lev], dmap[lev], ncomp, 0, MFInfo(), Factory(lev));
@@ -592,10 +595,10 @@ void incflo::WritePlotFile()
 
     if (m_plt_eta) {
 #ifdef USE_AMREX_MPMD
-    // Call to indicate this is not final data-transfer
+    // Call to indicate this is from incflo::WritePlotFile
     if (ParallelDescriptor::MyProc() == 0) {
         Vector<int> last_call;
-        last_call.push_back(0);
+        last_call.push_back(-1);
         MPI_Send(last_call.data(), last_call.size(), MPI_INT,m_mpmd_other_root,94,MPI_COMM_WORLD);
     }
 #endif
@@ -714,9 +717,7 @@ void incflo::WritePlotFile()
 #endif
 
     if (m_nodal_vel_eta && m_plt_inertial_num
-        && (m_fluid_model == FluidModel::DataDrivenMPMD
-            || m_fluid_model == FluidModel::Rauter
-            || m_fluid_model_second == FluidModel::DataDrivenMPMD
+        && (m_fluid_model_second == FluidModel::DataDrivenMPMD
             || m_fluid_model_second == FluidModel::Rauter)) {
         for (int lev = 0; lev <= finest_level; ++lev) {
             MultiFab inertial_num(amrex::convert(mf[lev].boxArray(),
@@ -737,6 +738,39 @@ void incflo::WritePlotFile()
         ++icomp;
     }
 
+#ifdef USE_AMREX_MPMD
+    if (m_nodal_vel_eta && m_plt_mu_I
+        && m_fluid_model_second == FluidModel::DataDrivenMPMD) {
+        // Call to indicate this is from incflo::WritePlotFile
+        if (ParallelDescriptor::MyProc() == 0) {
+            Vector<int> last_call;
+            last_call.push_back(-1);
+            MPI_Send(last_call.data(), last_call.size(), MPI_INT,m_mpmd_other_root,94,MPI_COMM_WORLD);
+        }
+        for (int lev = 0; lev <= finest_level; ++lev) {
+            MultiFab mu_I(amrex::convert(mf[lev].boxArray(),
+                                IndexType::TheNodeType().ixType()),
+                                mf[lev].DistributionMap(),1,0);
+            compute_nodal_inertial_num_at_level(lev,
+                                           &mu_I,
+                                           &m_leveldata[lev]->velocity,
+                                           &m_leveldata[lev]->p_nd,
+                                           m_mu_p_eps_second,
+                                           m_ro_grain_second,
+                                           m_diam_second,
+                                           Geom(lev),
+                                           m_cur_time, 0);
+            // Copier send inertial_num
+            mpmd_copiers_send_lev(mu_I,0,1,lev);
+            // Receive mu(I) = stress ratio
+            mpmd_copiers_recv_lev(mu_I,0,1,lev);
+            amrex::average_node_to_cellcenter(mf[lev],icomp,mu_I,0,1);
+        }
+        pltscaVarsName.push_back("mu_I");
+        ++icomp;
+    }
+#endif
+
 #ifdef AMREX_USE_EB
     for (int lev = 0; lev <= finest_level; ++lev) {
         EB_set_covered(mf[lev], 0.0);