add magvel as an optional derived quantity (AMReX-Fluids#109)

siddanib · Jun 20, 2024 · a101f71 · a101f71
1 parent a9a7df4
commit a101f71
Show file tree

Hide file tree

Showing 5 changed files with 95 additions and 5 deletions.
diff --git a/src/derive/incflo_derive.cpp b/src/derive/incflo_derive.cpp
@@ -142,6 +142,78 @@ Real incflo::ComputeKineticEnergy ()
     return 0;
 }
 
+#ifdef AMREX_USE_EB
+void incflo::ComputeMagVel (int lev,
+#else
+void incflo::ComputeMagVel (int /*lev*/,
+#endif
+                            Real /*time*/,
+                            MultiFab& magvel, MultiFab const& vel)
+{
+    BL_PROFILE("incflo::ComputeMagVel");
+
+#ifdef AMREX_USE_EB
+    auto const& fact = EBFactory(lev);
+    auto const& flags_mf = fact.getMultiEBCellFlagFab();
+#endif
+
+#ifdef _OPENMP
+#pragma omp parallel if (Gpu::notInLaunchRegion())
+#endif
+    for(MFIter mfi(vel, TilingIfNotGPU()); mfi.isValid(); ++mfi)
+    {
+        Box bx = mfi.tilebox();
+        Array4<Real const> const& ccvel_fab = vel.const_array(mfi);
+        Array4<Real> const& magvel_fab = magvel.array(mfi);
+
+#ifdef AMREX_USE_EB
+        const EBCellFlagFab& flags = flags_mf[mfi];
+        auto typ = flags.getType(bx);
+        if (typ == FabType::covered)
+        {
+            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            {
+                magvel_fab(i,j,k) = Real(0.0);
+            });
+        }
+        else
+        {
+            const auto& flag_fab = flags.const_array();
+            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            {
+                if (flag_fab(i,j,k).isCovered())
+                {
+                    magvel_fab(i,j,k) = Real(0.0);
+                }
+                else
+                {
+                    Real u = ccvel_fab(i,j,k,0);
+                    Real v = ccvel_fab(i,j,k,1);
+#if (AMREX_SPACEDIM == 2)
+                    magvel_fab(i,j,k) = std::sqrt(u*u + v*v);
+#elif  (AMREX_SPACEDIM == 3)
+                    Real w = ccvel_fab(i,j,k,2);
+                    magvel_fab(i,j,k) = std::sqrt(u*u + v*v + w*w);
+#endif
+                }
+            });
+        }
+#else       // No EB
+            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            {
+                Real u = ccvel_fab(i,j,k,0);
+                Real v = ccvel_fab(i,j,k,1);
+#if (AMREX_SPACEDIM == 2)
+                magvel_fab(i,j,k) = std::sqrt(u*u + v*v);
+#elif (AMREX_SPACEDIM == 3)
+                Real w = ccvel_fab(i,j,k,2);
+                magvel_fab(i,j,k) = std::sqrt(u*u + v*v + w*w);
+#endif
+            });
+#endif
+    } // mfi
+}
+
 #if (AMREX_SPACEDIM == 2)
 void incflo::ComputeVorticity (int lev, Real /*time*/, MultiFab& vort, MultiFab const& vel)
 {

diff --git a/src/derive/incflo_derive_K.H b/src/derive/incflo_derive_K.H
@@ -1,5 +1,5 @@
-#ifndef DERIVE_K_H_
-#define DERIVE_K_H_
+#ifndef INCFLO_DERIVE_K_H_
+#define INCFLO_DERIVE_K_H_
 
 #include <AMReX_FArrayBox.H>
 #include <cmath>

diff --git a/src/incflo.H b/src/incflo.H
@@ -174,11 +174,14 @@ public:
     //
     ///////////////////////////////////////////////////////////////////////////
 
-    void ComputeVorticity (int lev, amrex::Real time, amrex::MultiFab& vort,
-                           amrex::MultiFab const& vel);
+#ifndef AMREX_USE_EB
+    static
+#endif
+    void ComputeMagVel    (int lev, amrex::Real time, amrex::MultiFab& magvel, amrex::MultiFab const& vel);
+    void ComputeVorticity (int lev, amrex::Real time, amrex::MultiFab&   vort, amrex::MultiFab const& vel);
     void ComputeDivU (amrex::Real time);
     void ComputeDrag ();
-    [[nodiscard]] static amrex::Real ComputeKineticEnergy () ;
+    [[nodiscard]] static amrex::Real ComputeKineticEnergy ();
 
     void DiffFromExact (int lev, amrex::Geometry& lev_geom, amrex::Real time, amrex::Real dt,
                         amrex::MultiFab& error, int soln_comp, int err_comp) const;
@@ -544,6 +547,7 @@ private:
     int m_plt_p_nd        = 0;
     int m_plt_macphi      = 0;
     int m_plt_eta         = 0;
+    int m_plt_magvel      = 1;
     int m_plt_vort        = 1;
     int m_plt_forcing     = 0;
     int m_plt_strainrate  = 0;

diff --git a/src/setup/init.cpp b/src/setup/init.cpp
@@ -247,6 +247,7 @@ void incflo::ReadIOParameters()
         m_plt_macphi     = 0;
         m_plt_eta        = 0;
         m_plt_vort       = 0;
+        m_plt_magvel     = 0;
         m_plt_strainrate = 0;
         m_plt_divu       = 0;
         m_plt_vfrac      = 0;
@@ -270,6 +271,7 @@ void incflo::ReadIOParameters()
     pp.query("plt_p_nd",       m_plt_p_nd  );
     pp.query("plt_macphi",     m_plt_macphi);
     pp.query("plt_eta",        m_plt_eta   );
+    pp.query("plt_magvel",     m_plt_magvel);
     pp.query("plt_vort",       m_plt_vort  );
     pp.query("plt_strainrate", m_plt_strainrate);
     pp.query("plt_divu",       m_plt_divu  );

diff --git a/src/utilities/io.cpp b/src/utilities/io.cpp
@@ -393,6 +393,9 @@ void incflo::WritePlotFile()
     // Apparent viscosity
     if(m_plt_eta) ++ncomp;
 
+    // Magnitude of velocity
+    if(m_plt_magvel) ++ncomp;
+
     // Vorticity
     if(m_plt_vort) ++ncomp;
 
@@ -590,6 +593,15 @@ void incflo::WritePlotFile()
         pltscaVarsName.push_back("eta");
         ++icomp;
     }
+    if (m_plt_magvel) {
+        for (int lev = 0; lev <= finest_level; ++lev) {
+            MultiFab magvel(mf[lev], amrex::make_alias, icomp, 1);
+            ComputeMagVel(lev, m_cur_time, magvel, m_leveldata[lev]->velocity);
+        }
+        pltscaVarsName.push_back("magvel");
+        ++icomp;
+    }
+
     if (m_plt_vort) {
         for (int lev = 0; lev <= finest_level; ++lev) {
             (m_leveldata[lev]->velocity).FillBoundary(geom[lev].periodicity());