Added function to sum Kinetic Energy (Exawind#34)

* added function that computes kinetic energy and an input to control frequency. Added 3D Taylor-Green Vortex initial condition and input file * moved a FAB based volmask to an iMultiFab level_mask. Level mask does not have to perform box intersection everytime kinetic energy is called but instead only when regrid is called. * cleaned up a bit and now compiles with latest AMReX development branch * updated ComputeKineticEnergy to use the three MultiFab version of ReduceSum and divided by freestream density ro_0. Use latest version of AMReX development to compile. * Removed level_mask setup function since it already exists in AMReX as makeFineMask, thanks Weiqun. * moved the ComputeKineticEnergy to public for CUDA lambdas to work, and flipped the masking number (0,1)->(1,0), thanks Weiqun
ewquon · Oct 31, 2019 · c9b991d · c9b991d
1 parent 2396738
commit c9b991d
Show file tree

Hide file tree

Showing 7 changed files with 177 additions and 3 deletions.
diff --git a/exec/inputs.taylor_green_vortices3d b/exec/inputs.taylor_green_vortices3d
@@ -0,0 +1,59 @@
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#            SIMULATION STOP            #
+#.......................................#
+stop_time               =   3.183098861837907   # Max (simulated) time to evolve
+#max_step                =   1          # Max number of time steps
+steady_state            =   0           # Steady-state solver? 
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#         TIME STEP COMPUTATION         #
+#.......................................#
+incflo.fixed_dt         =   .005         # Use this constant dt if > 0
+incflo.cfl              =   0.7         # CFL factor
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#            INPUT AND OUTPUT           #
+#.......................................#
+amr.KE_int = 1
+amr.plot_int            =   10          # Steps between plot files
+#amr.plot_per            =   0.1         # Steps between plot files
+amr.check_int           =   -1          # Steps between checkpoint files
+amr.restart             =   ""          # Checkpoint to restart from 
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#               PHYSICS                 #
+#.......................................#
+incflo.gravity          =   0.  0.  0.  # Gravitational force (3D)
+incflo.ro_0             =   1.          # Reference density 
+# Re = rho*u*L/mu = 1*1*1*L/mu = 1600, L = 1/2/pi, mu = 1/1600/2/pi
+incflo.mu               =  0.00009947183943  # Dynamic viscosity coefficient
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#        ADAPTIVE MESH REFINEMENT       #
+#.......................................#
+amr.n_cell              =   64 64 64       # Grid cells at coarsest AMRlevel
+amr.max_level           =   0           # Max AMR level in hierarchy 
+##amr.blocking_factor = 16
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#              GEOMETRY                 #
+#.......................................#
+geometry.prob_lo        =   0.  0.  0.  # Lo corner coordinates
+geometry.prob_hi        =   1.  1.  1.  # Hi corner coordinates
+geometry.is_periodic    =   1   1   1   # Periodicity x y z (0/1)
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#          INITIAL CONDITIONS           #
+#.......................................#
+incflo.probtype         =   3
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#          NUMERICAL PARAMETERS         #
+#.......................................#
+incflo.steady_state_tol =   1.e-5       # Tolerance for steady-state
+amrex.fpe_trap_invalid  =   1           # Trap NaNs
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#              VERBOSITY                #
+#.......................................#
+incflo.verbose          =   1           # incflo_level
+mac.verbose             =   0           # MacProjector
diff --git a/src/derive/derive.cpp b/src/derive/derive.cpp
@@ -219,6 +219,49 @@ void incflo::ComputeStrainrate(Real time_in)
     }
 }
 
+
+Real incflo::ComputeKineticEnergy()
+{
+    BL_PROFILE("incflo::ComputeKineticEnergy");
+
+    // integrated total Kinetic energy
+    Real KE = 0.0;
+
+    for(int lev = 0; lev <= finest_level; lev++)
+    {
+        Real cell_vol = geom[lev].CellSize()[0]*geom[lev].CellSize()[1]*geom[lev].CellSize()[2];
+
+        KE += amrex::ReduceSum(*density[lev],*vel[lev],*level_mask[lev],0,
+        [=] AMREX_GPU_HOST_DEVICE (Box const& bx,
+                                   Array4<Real const> const& den_arr,
+                                   Array4<Real const> const& vel_arr,
+                                   Array4<int const>  const& mask_arr) -> Real
+        {
+            Real KE_Fab = 0.0;
+
+            amrex::Loop(bx, [=,&KE_Fab] (int i, int j, int k) noexcept
+            {
+                KE_Fab += cell_vol*mask_arr(i,j,k)*den_arr(i,j,k)*( vel_arr(i,j,k,0)*vel_arr(i,j,k,0)
+                                                                   +vel_arr(i,j,k,1)*vel_arr(i,j,k,1)
+                                                                   +vel_arr(i,j,k,2)*vel_arr(i,j,k,2));
+
+            });
+            return KE_Fab;
+
+        });
+    }
+
+    // total volume of grid on level 0
+    Real total_vol = geom[0].ProbDomain().volume();
+
+    KE *= 0.5/total_vol/ro_0;
+
+    ParallelDescriptor::ReduceRealSum(KE);
+
+    return KE;
+
+}
+
 void incflo::ComputeVorticity(Real time_in)
 {
 	BL_PROFILE("incflo::ComputeVorticity");

diff --git a/src/incflo.H b/src/incflo.H
@@ -159,6 +159,8 @@ public:
     //////////////////////////////////////////////////////////////////////////////////////////////
 
     void ComputeVorticity(Real time_in);
+    double ComputeKineticEnergy();
+
 
     //////////////////////////////////////////////////////////////////////////////////////////////
     //
@@ -264,7 +266,7 @@ private:
     // Resize arrays to fit (up to) max_level + 1 AMR levels
 	void ResizeArrays();
 
-    // Post-initialization: set BCs, apply ICs, initial velocity projection, pressure iterations
+    // Post-initialization: set BCs, apply ICs, initial velocity projection, pressure iterations, level_mask
 	void PostInit(int restart_flag);
 	void SetBCTypes();
     void InitFluid();
@@ -370,7 +372,7 @@ private:
     //////////////////////////////////////////////////////////////////////////////////////////////
 
     void ComputeDrag();
-
+    
     //////////////////////////////////////////////////////////////////////////////////////////////
     //
     // Embedded Boundaries
@@ -484,6 +486,7 @@ private:
 
     int check_int = -1;
     int last_chk = -1;
+    int KE_int = -1;
     std::string check_file{"chk"};
     std::string restart_file{""};
 
@@ -588,6 +591,10 @@ private:
     Vector<std::unique_ptr<MultiFab>> m_u_mac;
     Vector<std::unique_ptr<MultiFab>> m_v_mac;
     Vector<std::unique_ptr<MultiFab>> m_w_mac;
+
+    // integer MultiFab specifying if a finer level exists above
+    // value is 1 if a finer level is not above and 0 otherwise
+    Vector<std::unique_ptr<iMultiFab>> level_mask;
 
     //////////////////////////////////////////////////////////////////////////////////////////////
     //

diff --git a/src/incflo.cpp b/src/incflo.cpp
@@ -97,6 +97,10 @@ void incflo::InitData()
         WritePlotFile();
         last_plt = 0;
     }
+    if(KE_int > 0 && !restart_flag)
+    {
+        amrex::Print() << "Time, Kinetic Energy: " << cur_time << ", " << ComputeKineticEnergy() << std::endl;
+    }
 
     ParmParse pp("incflo");
     bool write_eb_surface = 0;
@@ -200,6 +204,7 @@ void incflo::Evolve()
                     incflo_setup_solvers();
                 }
             }
+         
         }*/
 
         // Advance to time t + dt
@@ -219,6 +224,11 @@ void incflo::Evolve()
             WriteCheckPointFile();
             last_chk = nstep;
         }
+
+        if(KE_int > 0 && (nstep % KE_int == 0))
+        {
+            amrex::Print() << "Time, Kinetic Energy: " << cur_time << ", " << ComputeKineticEnergy() << std::endl;
+        }
 
         // Mechanism to terminate incflo normally.
         do_not_evolve = (steady_state && SteadyStateReached()) ||

diff --git a/src/setup/incflo_arrays.cpp b/src/setup/incflo_arrays.cpp
@@ -8,6 +8,13 @@ void incflo::AllocateArrays(int lev)
     // Cell-based arrays
     // ********************************************************************************
 
+    if(lev < finest_level){
+        level_mask[lev].reset(new iMultiFab(makeFineMask(grids[lev],dmap[lev], grids[lev+1], IntVect(2), 1, 0)));
+    } else {
+        level_mask[lev].reset(new iMultiFab(grids[lev], dmap[lev], 1, 0, MFInfo() /*, default factory*/));
+        level_mask[lev]->setVal(1);
+    }
+
     // Current Density
     density[lev].reset(new MultiFab(grids[lev], dmap[lev], 1, nghost, MFInfo(), *ebfactory[lev]));
     density[lev]->setVal(0.);
@@ -172,6 +179,13 @@ void incflo::RegridArrays(int lev)
     // FillBoundary().
     //
 
+    if(lev < finest_level){
+        level_mask[lev].reset(new iMultiFab(makeFineMask(grids[lev],dmap[lev], grids[lev+1], IntVect(2), 1, 0)));
+    } else {
+        level_mask[lev].reset(new iMultiFab(grids[lev], dmap[lev], 1, 0, MFInfo() /*, default factory*/));
+        level_mask[lev]->setVal(1);
+    }
+
    // Density
    std::unique_ptr<MultiFab> density_new(new MultiFab(grids[lev], dmap[lev], 1, nghost, 
                                                       MFInfo(), *ebfactory[lev]));
@@ -506,6 +520,8 @@ void incflo::ResizeArrays()
 
     // EB factory
     ebfactory.resize(max_level + 1);
+
+    level_mask.resize(max_level + 1);
 }
 
 void incflo::MakeBCArrays()

diff --git a/src/setup/init.cpp b/src/setup/init.cpp
@@ -33,6 +33,8 @@ void incflo::ReadParameters()
 	pp.query("plot_int", plot_int);
 	pp.query("plot_per", plot_per);
 
+    pp.query("KE_int", KE_int);
+
         // Which variables to write to plotfile
         pltVarCount = 0;
 
@@ -247,7 +249,7 @@ void incflo::PostInit(int restart_flag)
     {
         InitFluid();
     }
-
+    
     // Set the background pressure and gradients in "DELP" cases
     SetBackgroundPressure();
 

diff --git a/src/setup/init_fluid.f90 b/src/setup/init_fluid.f90
@@ -45,6 +45,7 @@ subroutine init_fluid(slo, shi, lo, hi, &
 
       if (probtype ==  1) call taylor_green(lo, hi, vel, slo, shi, dx, dy, dz, domlo)
       if (probtype ==  2) call double_shear_layer(lo, hi, vel, slo, shi, dx, dy, dz, domlo)
+      if (probtype ==  3) call taylor_green3d(lo, hi, vel, slo, shi, dx, dy, dz, domlo)
       if (probtype == 31) call plane_poiseuille(lo, hi, vel, tracer, slo, shi, dx, dy, dz, domlo, domhi, probtype)
       if (probtype == 32) call plane_poiseuille(lo, hi, vel, tracer, slo, shi, dx, dy, dz, domlo, domhi, probtype)
       if (probtype == 33) call plane_poiseuille(lo, hi, vel, tracer, slo, shi, dx, dy, dz, domlo, domhi, probtype)
@@ -90,6 +91,42 @@ subroutine taylor_green(lo, hi, vel, slo, shi, dx, dy, dz, domlo)
 
    end subroutine taylor_green
 
+    subroutine taylor_green3d(lo, hi, vel, slo, shi, dx, dy, dz, domlo)
+
+       use constant, only: zero, half, one
+
+       implicit none
+
+       ! Array bounds
+       integer(c_int), intent(in   ) :: lo(3), hi(3)
+       integer(c_int), intent(in   ) :: domlo(3)
+       integer(c_int), intent(in   ) :: slo(3), shi(3)
+
+       ! Arrays
+       real(rt),       intent(inout) :: vel(slo(1):shi(1),slo(2):shi(2),slo(3):shi(3),3)
+
+       real(rt),       intent(in   ) :: dx, dy, dz
+
+       ! Local variables
+       integer(c_int)                  :: i, j, k
+       real(rt)                        :: x, y, z
+       real(rt)                        :: twopi = 8.0_rt * atan(one)
+
+       do j = lo(2), hi(2)
+          y =  ( real(j,rt) + half ) * dy
+          do i = lo(1), hi(1)
+             x =  ( real(i,rt) + half ) * dx
+             do k = lo(3), hi(3)
+                z =  ( real(k,rt) + half ) * dz
+                vel(i,j,k,1) =  sin(twopi * x) * cos(twopi * y) * cos(twopi * z)
+                vel(i,j,k,2) = -cos(twopi * x) * sin(twopi * y) * cos(twopi * z)
+                vel(i,j,k,3) = zero
+             end do
+          end do
+       end do
+
+    end subroutine taylor_green3d
+
    subroutine double_shear_layer(lo, hi, vel, slo, shi, dx, dy, dz, domlo)
 
       use constant, only: zero, half, one