add magvel as a derived quantity

AMReX-Fluids · Apr 2, 2024 · 02b6c60 · 02b6c60
1 parent a732e84
commit 02b6c60
Show file tree

Hide file tree

Showing 5 changed files with 135 additions and 2 deletions.
diff --git a/src/derive/incflo_derive.cpp b/src/derive/incflo_derive.cpp
@@ -143,6 +143,65 @@ Real incflo::ComputeKineticEnergy ()
 }
 
 #if (AMREX_SPACEDIM == 2)
+void incflo::ComputeMagVel (int lev, Real /*time*/, MultiFab& magvel, MultiFab const& vel)
+{
+    BL_PROFILE("incflo::ComputeMagVel");
+
+#ifdef AMREX_USE_EB
+    const auto& fact = EBFactory(lev);
+    const auto& 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();
+            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            {
+                constexpr Real c0 = Real(-1.5);
+                constexpr Real c1 = Real( 2.0);
+                constexpr Real c2 = Real(-0.5);
+
+                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);
+                    magvel_fab(i,j,k) = std::sqrt(u*u + v*v);
+                }
+            });
+        }
+#else
+            amrex::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);
+                magvel_fab(i,j,k) = std::sqrt(u*u + v*v);
+            });
+#endif
+    } // mfi
+}
 void incflo::ComputeVorticity (int lev, Real /*time*/, MultiFab& vort, MultiFab const& vel)
 {
     BL_PROFILE("incflo::ComputeVorticity");
@@ -245,6 +304,63 @@ void incflo::ComputeVorticity (int lev, Real /*time*/, MultiFab& vort, MultiFab
 }
 
 #elif (AMREX_SPACEDIM == 3)
+void incflo::ComputeMagVel (int lev, Real /*time*/, MultiFab& magvel, MultiFab const& vel)
+{
+    BL_PROFILE("incflo::ComputeMagVel");
+
+#ifdef AMREX_USE_EB
+    const auto& fact = EBFactory(lev);
+    const auto& 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();
+            amrex::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);
+                    Real w = ccvel_fab(i,j,k,2);
+                    magvel_fab(i,j,k) = std::sqrt(u*u + v*v + w*w);
+                }
+            });
+        }
+#else
+            amrex::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);
+                Real w = ccvel_fab(i,j,k,2);
+                magvel_fab(i,j,k) = std::sqrt(u*u + v*v + w*w);
+            });
+#endif
+    } // mfi
+}
 void incflo::ComputeVorticity (int lev, Real /*time*/, MultiFab& vort, MultiFab const& vel)
 {
     BL_PROFILE("incflo::ComputeVorticity");

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
@@ -162,6 +162,8 @@ public:
     //
     ///////////////////////////////////////////////////////////////////////////
 
+    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);
@@ -531,6 +533,7 @@ private:
     int m_plt_macphi      = 0;
     int m_plt_eta         = 0;
     int m_plt_vort        = 1;
+    int m_plt_magvel      = 1;
     int m_plt_forcing     = 0;
     int m_plt_strainrate  = 0;
     int m_plt_divu        = 0;

diff --git a/src/setup/init.cpp b/src/setup/init.cpp
@@ -233,6 +233,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;
@@ -253,6 +254,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
@@ -381,6 +381,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;
 
@@ -582,6 +585,15 @@ void incflo::WritePlotFile()
         pltscaVarsName.push_back("vort");
         ++icomp;
     }
+    if (m_plt_magvel) {
+        for (int lev = 0; lev <= finest_level; ++lev) {
+            (m_leveldata[lev]->velocity).FillBoundary(geom[lev].periodicity());
+            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_forcing) {
         for (int lev = 0; lev <= finest_level; ++lev) {
             MultiFab forcing(mf[lev], amrex::make_alias, icomp, 3);