add all components of vorticity as derived quantities (erf-model#1564)

* add all components of vorticity as derived quantities * fix for CUDA
ewquon · Apr 6, 2024 · dccea6f · dccea6f
1 parent d6ceb6e
commit dccea6f
Show file tree

Hide file tree

Showing 9 changed files with 131 additions and 31 deletions.
diff --git a/Source/BoundaryConditions/ERF_FillPatch.cpp b/Source/BoundaryConditions/ERF_FillPatch.cpp
@@ -532,3 +532,93 @@ ERF::FillCoarsePatch (int lev, Real time)
     // Since lev > 0 here we don't worry about m_r2d or wrfbdy data
     // ***************************************************************************
 }
+
+void
+ERF::FillBdyCCVels (Vector<MultiFab>& mf_cc_vel)
+{
+    // Impose bc's at domain boundaries
+    for (int lev = 0; lev <= finest_level; ++lev)
+    {
+        Box domain(Geom(lev).Domain());
+
+        int ihi = domain.bigEnd(0);
+        int jhi = domain.bigEnd(1);
+        int khi = domain.bigEnd(2);
+
+        // Impose periodicity first
+        mf_cc_vel[lev].FillBoundary(geom[lev].periodicity());
+
+        for (MFIter mfi(mf_cc_vel[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
+        {
+            const Box& bx = mfi.growntilebox(1);
+            const Array4<Real>& vel_arr = mf_cc_vel[lev].array(mfi);
+
+            if (!Geom(lev).isPeriodic(0)) {
+                // Low-x side
+                if (bx.smallEnd(0) <= domain.smallEnd(0)) {
+                    Real mult = (phys_bc_type[0] == ERF_BC::no_slip_wall) ? -1. : 1.;
+                    ParallelFor(makeSlab(bx,0,0), [=] AMREX_GPU_DEVICE(int , int j, int k) noexcept
+                    {
+                        vel_arr(-1,j,k,1) = mult*vel_arr(0,j,k,1); // v
+                        vel_arr(-1,j,k,2) = mult*vel_arr(0,j,k,2); // w
+                    });
+                }
+
+                // High-x side
+                if (bx.bigEnd(0) >= domain.bigEnd(0)) {
+                    Real mult = (phys_bc_type[3] == ERF_BC::no_slip_wall) ? -1. : 1.;
+                    ParallelFor(makeSlab(bx,0,0), [=] AMREX_GPU_DEVICE(int , int j, int k) noexcept
+                    {
+                        vel_arr(ihi+1,j,k,1) = mult*vel_arr(ihi,j,k,1); // v
+                        vel_arr(ihi+1,j,k,2) = mult*vel_arr(ihi,j,k,2); // w
+                    });
+                }
+            } // !periodic
+
+            if (!Geom(lev).isPeriodic(1)) {
+                // Low-y side
+                if (bx.smallEnd(1) <= domain.smallEnd(1)) {
+                    Real mult = (phys_bc_type[1] == ERF_BC::no_slip_wall) ? -1. : 1.;
+                    ParallelFor(makeSlab(bx,1,0), [=] AMREX_GPU_DEVICE(int i, int  , int k) noexcept
+                    {
+                        vel_arr(i,-1,k,0) = mult*vel_arr(i,0,k,0); // u
+                        vel_arr(i,-1,k,2) = mult*vel_arr(i,0,k,2); // w
+                    });
+                }
+
+                // High-y side
+                if (bx.bigEnd(1) >= domain.bigEnd(1)) {
+                    Real mult = (phys_bc_type[4] == ERF_BC::no_slip_wall) ? -1. : 1.;
+                    ParallelFor(makeSlab(bx,1,0), [=] AMREX_GPU_DEVICE(int i, int , int k) noexcept
+                    {
+                        vel_arr(i,jhi+1,k,0) = mult*vel_arr(i,jhi,k,0); // u
+                        vel_arr(i,jhi+1,k,2) = mult*-vel_arr(i,jhi,k,2); // w
+                    });
+                }
+            } // !periodic
+
+            if (!Geom(lev).isPeriodic(2)) {
+                // Low-z side
+                if (bx.smallEnd(2) <= domain.smallEnd(2)) {
+                    Real mult = (phys_bc_type[2] == ERF_BC::no_slip_wall) ? -1. : 1.;
+                    ParallelFor(makeSlab(bx,2,0), [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
+                    {
+                        vel_arr(i,j,-1,0) = mult*vel_arr(i,j,0,0); // u
+                        vel_arr(i,j,-1,1) = mult*vel_arr(i,j,0,1); // v
+                    });
+                }
+
+                // High-z side
+                if (bx.bigEnd(2) >= domain.bigEnd(2)) {
+                    Real mult = (phys_bc_type[5] == ERF_BC::no_slip_wall) ? -1. : 1.;
+                    ParallelFor(makeSlab(bx,2,0), [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
+                    {
+                        vel_arr(i,j,khi+1,0) = mult*vel_arr(i,j,khi,0); // u
+                        vel_arr(i,j,khi+1,1) = mult*vel_arr(i,j,khi,1); // v
+                    });
+                }
+            } // !periodic
+        } // MFIter
+
+    } // lev
+}
diff --git a/Source/ERF.H b/Source/ERF.H
@@ -147,6 +147,9 @@ public:
     void write_1D_profiles (amrex::Real time);
     void write_1D_profiles_stag (amrex::Real time);
 
+    // Fill the physical boundary conditions for cell-centered velocity (diagnostic only)
+    void FillBdyCCVels (amrex::Vector<amrex::MultiFab>& mf_cc_vel);
+
     void sample_points (int lev, amrex::Real time, amrex::IntVect cell, amrex::MultiFab& mf);
     void sample_lines (int lev, amrex::Real time, amrex::IntVect cell, amrex::MultiFab& mf);
 

diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp
@@ -203,31 +203,36 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
         containerHasElement(plot_var_names, "vorticity_z") ) {
 
         for (int lev = 0; lev <= finest_level; ++lev) {
-            mf_cc_vel[lev].define(grids[lev], dmap[lev], AMREX_SPACEDIM, IntVect(1,1,0));
+            mf_cc_vel[lev].define(grids[lev], dmap[lev], AMREX_SPACEDIM, IntVect(1,1,1));
             average_face_to_cellcenter(mf_cc_vel[lev],0,
-                                       Array<const MultiFab*,3>{&vars_new[lev][Vars::xvel],&vars_new[lev][Vars::yvel],&vars_new[lev][Vars::zvel]});
-            mf_cc_vel[lev].FillBoundary(geom[lev].periodicity());
+                                       Array<const MultiFab*,3>{&vars_new[lev][Vars::xvel],
+                                                                &vars_new[lev][Vars::yvel],
+                                                                &vars_new[lev][Vars::zvel]});
         } // lev
+    } // if (vel or vort)
 
-        // We need ghost cells if computing vorticity
+    // We need ghost cells if computing vorticity
+    if ( containerHasElement(plot_var_names, "vorticity_x")||
+         containerHasElement(plot_var_names, "vorticity_y") ||
+         containerHasElement(plot_var_names, "vorticity_z") )
+    {
         amrex::Interpolater* mapper = &cell_cons_interp;
-        if ( containerHasElement(plot_var_names, "vorticity_x")||
-             containerHasElement(plot_var_names, "vorticity_y") ||
-             containerHasElement(plot_var_names, "vorticity_z") ) {
-            for (int lev = 1; lev <= finest_level; ++lev) {
-                Vector<MultiFab*> fmf = {&(mf_cc_vel[lev]), &(mf_cc_vel[lev])};
-                Vector<Real> ftime    = {t_new[lev], t_new[lev]};
-                Vector<MultiFab*> cmf = {&mf_cc_vel[lev-1], &mf_cc_vel[lev-1]};
-                Vector<Real> ctime    = {t_new[lev], t_new[lev]};
-
-                MultiFab mf_to_fill;
-                amrex::FillPatchTwoLevels(mf_cc_vel[lev], t_new[lev], cmf, ctime, fmf, ftime,
-                                          0, 0, AMREX_SPACEDIM, geom[lev-1], geom[lev],
-                                          null_bc_for_fill, 0, null_bc_for_fill, 0, refRatio(lev-1),
-                                          mapper, domain_bcs_type, 0);
-            } // lev
-        } // if
-    } // if
+        for (int lev = 1; lev <= finest_level; ++lev)
+        {
+            Vector<MultiFab*> fmf = {&(mf_cc_vel[lev]), &(mf_cc_vel[lev])};
+            Vector<Real> ftime    = {t_new[lev], t_new[lev]};
+            Vector<MultiFab*> cmf = {&mf_cc_vel[lev-1], &mf_cc_vel[lev-1]};
+            Vector<Real> ctime    = {t_new[lev], t_new[lev]};
+
+            amrex::FillPatchTwoLevels(mf_cc_vel[lev], t_new[lev], cmf, ctime, fmf, ftime,
+                                      0, 0, AMREX_SPACEDIM, geom[lev-1], geom[lev],
+                                      null_bc_for_fill, 0, null_bc_for_fill, 0, refRatio(lev-1),
+                                      mapper, domain_bcs_type, 0);
+        } // lev
+
+        // Impose bc's at domain boundaries at all levels
+        FillBdyCCVels(mf_cc_vel);
+    } // if (vort)
 
     for (int lev = 0; lev <= finest_level; ++lev)
     {

diff --git a/Tests/CTestList.cmake b/Tests/CTestList.cmake
@@ -39,7 +39,7 @@ function(add_test_r TEST_NAME TEST_EXE PLTFILE)
 
     set(TEST_EXE ${CMAKE_BINARY_DIR}/Exec/${TEST_EXE})
     set(FCOMPARE_TOLERANCE "-r 2e-10 --abs_tol 2.0e-10")
-    set(FCOMPARE_FLAGS "-a ${FCOMPARE_TOLERANCE}")
+    set(FCOMPARE_FLAGS "--abort_if_not_all_found -a ${FCOMPARE_TOLERANCE}")
     set(test_command sh -c "${MPI_COMMANDS} ${TEST_EXE} ${CURRENT_TEST_BINARY_DIR}/${TEST_NAME}.i ${RUNTIME_OPTIONS} > ${TEST_NAME}.log && ${MPI_FCOMP_COMMANDS} ${FCOMPARE_EXE} ${FCOMPARE_FLAGS} ${PLOT_GOLD} ${CURRENT_TEST_BINARY_DIR}/${PLTFILE}")
 
     add_test(${TEST_NAME} ${test_command})

diff --git a/Tests/ERFGoldFiles/IsentropicVortexAdvecting/Header b/Tests/ERFGoldFiles/IsentropicVortexAdvecting/Header
@@ -1,11 +1,13 @@
 HyperCLaw-V1.1
-8
+10
 density
 x_velocity
 y_velocity
 z_velocity
 temp
 theta
+vorticity_x
+vorticity_y
 vorticity_z
 pressure
 3

diff --git a/Tests/ERFGoldFiles/IsentropicVortexAdvecting/Level_0/Cell_D_00000 b/Tests/ERFGoldFiles/IsentropicVortexAdvecting/Level_0/Cell_D_00000
diff --git a/Tests/ERFGoldFiles/IsentropicVortexAdvecting/Level_0/Cell_D_00001 b/Tests/ERFGoldFiles/IsentropicVortexAdvecting/Level_0/Cell_D_00001
diff --git a/Tests/ERFGoldFiles/IsentropicVortexAdvecting/Level_0/Cell_H b/Tests/ERFGoldFiles/IsentropicVortexAdvecting/Level_0/Cell_H
@@ -1,6 +1,6 @@
 1
 1
-8
+10
 0
 (2 0
 ((0,0,0) (47,23,3) (0,0,0))
@@ -10,11 +10,11 @@
 FabOnDisk: Cell_D_00000 0
 FabOnDisk: Cell_D_00001 0
 
-2,8
-1.0999792802791342e+00,2.9304936277992795e+02,2.3942911725482566e+02,0.0000000000000000e+00,2.9354610160293907e+02,2.9999999999999983e+02,-9.4494008092048944e+01,9.2670758686574656e+04,
-6.7266824942273507e-01,9.8982769986740919e+01,8.2916282642571346e+01,0.0000000000000000e+00,2.4112535881735599e+02,2.9999999999999977e+02,-9.4845419180726424e+01,4.6550646053037621e+04,
+2,10
+1.0999792802791342e+00,2.9304936277992795e+02,2.3942911725482566e+02,0.0000000000000000e+00,2.9354610160293907e+02,2.9999999999999983e+02,0.0000000000000000e+00,0.0000000000000000e+00,-9.4494008092048944e+01,9.2670758686574656e+04,
+6.7266824942273507e-01,9.8982769986740919e+01,8.2916282642571346e+01,0.0000000000000000e+00,2.4112535881735599e+02,2.9999999999999977e+02,0.0000000000000000e+00,0.0000000000000000e+00,-9.4845419180726424e+01,4.6550646053037621e+04,
 
-2,8
-1.1658900337646416e+00,4.4005150729001542e+02,3.5146595531093567e+02,0.0000000000000000e+00,3.0045923092147854e+02,3.0000000000000023e+02,2.3952311355344591e-01,1.0053679536769103e+05,
-1.1656654423979136e+00,4.9220671551520172e+02,4.7708174379805683e+02,0.0000000000000000e+00,3.0043607798337882e+02,3.0000000000000023e+02,5.6424931117062329e+02,1.0050968272762455e+05,
+2,10
+1.1658900337646416e+00,4.4005150729001542e+02,3.5146595531093567e+02,0.0000000000000000e+00,3.0045923092147854e+02,3.0000000000000023e+02,0.0000000000000000e+00,0.0000000000000000e+00,2.3952311355344591e-01,1.0053679536769103e+05,
+1.1656654423979136e+00,4.9220671551520172e+02,4.7708174379805683e+02,0.0000000000000000e+00,3.0043607798337882e+02,3.0000000000000023e+02,0.0000000000000000e+00,0.0000000000000000e+00,5.6424931117062329e+02,1.0050968272762455e+05,
 
diff --git a/Tests/test_files/IsentropicVortexAdvecting/IsentropicVortexAdvecting.i b/Tests/test_files/IsentropicVortexAdvecting/IsentropicVortexAdvecting.i
@@ -34,7 +34,7 @@ erf.check_int       = 100        # number of timesteps between checkpoints
 # PLOTFILES
 erf.plot_file_1     = plt        # number of timesteps between plotfiles
 erf.plot_int_1      = 10         # number of timesteps between plotfiles
-erf.plot_vars_1     = density x_velocity y_velocity z_velocity pressure theta temp vorticity_z
+erf.plot_vars_1     = density x_velocity y_velocity z_velocity pressure theta temp vorticity_x vorticity_y vorticity_z
 
 # SOLVER CHOICE
 erf.alpha_T = 0.0