Add vorticity components (#1560)

* add components of vorticity to plotfile options

* add vorticity_z to the IsentropicVortexAdvecting CI test

Docs/sphinx_doc/Plotfiles.rst

@@ -201,6 +201,18 @@ Output Options
 |                             |                  |
 |                             |                  |
 +-----------------------------+------------------+
+| **vorticity_x*              | x-component of   |
+|                             | vorticity        |
+|                             |                  |
++-----------------------------+------------------+
+| **vorticity_y*              | y-component of   |
+|                             | vorticity        |
+|                             |                  |
++-----------------------------+------------------+
+| **vorticity_z*              | z-component of   |
+|                             | vorticity        |
+|                             |                  |
++-----------------------------+------------------+
 | **rhoadv_0**                | Conserved scalar |
 |                             |                  |
 |                             |                  |

Source/Derive.H

@@ -89,7 +89,29 @@ void erf_derQKE (
   const int* bcrec,
   const int level);
 
-void erf_dervort (
+void erf_dervortx (
+  const amrex::Box& bx,
+  amrex::FArrayBox& derfab,
+  int dcomp,
+  int ncomp,
+  const amrex::FArrayBox& datfab,
+  const amrex::Geometry& geomdata,
+  amrex::Real time,
+  const int* bcrec,
+  const int level);
+
+void erf_dervorty (
+  const amrex::Box& bx,
+  amrex::FArrayBox& derfab,
+  int dcomp,
+  int ncomp,
+  const amrex::FArrayBox& datfab,
+  const amrex::Geometry& geomdata,
+  amrex::Real time,
+  const int* bcrec,
+  const int level);
+
+void erf_dervortz (
   const amrex::Box& bx,
   amrex::FArrayBox& derfab,
   int dcomp,

Source/Derive.cpp

@@ -195,7 +195,7 @@ erf_derQKE (const Box& bx,
 }
 
 void
-erf_dervort(
+erf_dervortx (
   const amrex::Box& bx,
   amrex::FArrayBox& derfab,
   int dcomp,
@@ -206,13 +206,70 @@ erf_dervort(
   const int* /*bcrec*/,
   const int /*level*/)
 {
+    AMREX_ALWAYS_ASSERT(dcomp == 0);
     AMREX_ALWAYS_ASSERT(ncomp == 1);
 
     auto const dat = datfab.array(); // cell-centered velocity
-    auto tfab      = derfab.array(); // cell-centered vorticity
+    auto tfab      = derfab.array(); // cell-centered vorticity x-component
 
-    const Real dx = geomdata.CellSize(0);
     const Real dy = geomdata.CellSize(1);
+    const Real dz = geomdata.CellSize(2);
+
+    ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+    {
+        tfab(i,j,k,dcomp) = (dat(i,j+1,k,2) - dat(i,j-1,k,2)) / (2.0*dy)  // dw/dy
+                          - (dat(i,j,k+1,1) - dat(i,j,k-1,1)) / (2.0*dz); // dv/dz
+    });
+}
+
+void
+erf_dervorty (
+  const amrex::Box& bx,
+  amrex::FArrayBox& derfab,
+  int dcomp,
+  int ncomp,
+  const amrex::FArrayBox& datfab,
+  const amrex::Geometry& geomdata,
+  amrex::Real /*time*/,
+  const int* /*bcrec*/,
+  const int /*level*/)
+{
+    AMREX_ALWAYS_ASSERT(dcomp == 0);
+    AMREX_ALWAYS_ASSERT(ncomp == 1);
+
+    auto const dat = datfab.array(); // cell-centered velocity
+    auto tfab      = derfab.array(); // cell-centered vorticity y-component
+
+    const Real dx = geomdata.CellSize(0);
+    const Real dz = geomdata.CellSize(2);
+
+    ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+    {
+        tfab(i,j,k,dcomp) = (dat(i,j,k+1,0) - dat(i,j,k-1,0)) / (2.0*dz)  // du/dz
+                          - (dat(i+1,j,k,2) - dat(i-1,j,k,2)) / (2.0*dx); // dw/dx
+    });
+}
+
+void
+erf_dervortz (
+  const amrex::Box& bx,
+  amrex::FArrayBox& derfab,
+  int dcomp,
+  int ncomp,
+  const amrex::FArrayBox& datfab,
+  const amrex::Geometry& geomdata,
+  amrex::Real /*time*/,
+  const int* /*bcrec*/,
+  const int /*level*/)
+{
+    AMREX_ALWAYS_ASSERT(dcomp == 0);
+    AMREX_ALWAYS_ASSERT(ncomp == 1);
+
+    auto const dat = datfab.array(); // cell-centered velocity
+    auto tfab      = derfab.array(); // cell-centered vorticity z-component
+
+    const Real dx = geomdata.CellSize(0);
+    const Real dy = geomdata.CellSize(2);
 
     ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
     {
@@ -221,4 +278,5 @@ erf_dervort(
     });
 }
 
+
 } // namespace

Source/ERF.H

@@ -746,7 +746,8 @@ private:
                                                      "rhoQ4", "rhoQ5", "rhoQ6"};
 
     // Note that the order of variable names here must match the order in Derive.cpp
-    const amrex::Vector<std::string> derived_names {"soundspeed", "temp", "theta", "KE", "QKE", "scalar", "vorticity",
+    const amrex::Vector<std::string> derived_names {"soundspeed", "temp", "theta", "KE", "QKE", "scalar",
+                                                    "vorticity_x","vorticity_y","vorticity_z",
                                                     "pres_hse", "dens_hse", "pressure", "pert_pres", "pert_dens", "eq_pot_temp", "num_turb",
                                                     "dpdx", "dpdy", "pres_hse_x", "pres_hse_y",
                                                     "z_phys", "detJ" , "mapfac",

Source/IO/Plotfile.cpp

@@ -195,10 +195,12 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
     // Array of MultiFabs for cell-centered velocity
     Vector<MultiFab> mf_cc_vel(finest_level+1);
 
-    if (containerHasElement(plot_var_names, "x_velocity") ||
-        containerHasElement(plot_var_names, "y_velocity") ||
-        containerHasElement(plot_var_names, "z_velocity") ||
-        containerHasElement(plot_var_names, "vorticity") ) {
+    if (containerHasElement(plot_var_names, "x_velocity" ) ||
+        containerHasElement(plot_var_names, "y_velocity" ) ||
+        containerHasElement(plot_var_names, "z_velocity" ) ||
+        containerHasElement(plot_var_names, "vorticity_x") ||
+        containerHasElement(plot_var_names, "vorticity_y") ||
+        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));
@@ -209,7 +211,9 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
 
         // We need ghost cells if computing vorticity
         amrex::Interpolater* mapper = &cell_cons_interp;
-        if ( containerHasElement(plot_var_names, "vorticity") ) {
+        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]};
@@ -255,6 +259,7 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
         // Finally, check for any derived quantities and compute them, inserting
         // them into our output multifab
         auto calculate_derived = [&](const std::string& der_name,
+                                     MultiFab& src_mf,
                                      decltype(derived::erf_dernull)& der_function)
         {
             if (containerHasElement(plot_var_names, der_name)) {
@@ -266,7 +271,7 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
                 {
                     const Box& bx = mfi.tilebox();
                     auto& dfab = dmf[mfi];
-                    auto& sfab = vars_new[lev][Vars::cons][mfi];
+                    auto& sfab = src_mf[mfi];
                     der_function(bx, dfab, 0, 1, sfab, Geom(lev), t_new[0], nullptr, lev);
                 }
 
@@ -275,13 +280,15 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
         };
 
         // Note: All derived variables must be computed in order of "derived_names" defined in ERF.H
-        calculate_derived("soundspeed",  derived::erf_dersoundspeed);
-        calculate_derived("temp",        derived::erf_dertemp);
-        calculate_derived("theta",       derived::erf_dertheta);
-        calculate_derived("KE",          derived::erf_derKE);
-        calculate_derived("QKE",         derived::erf_derQKE);
-        calculate_derived("scalar",      derived::erf_derscalar);
-        calculate_derived("vorticity",   derived::erf_dervort);
+        calculate_derived("soundspeed",  vars_new[lev][Vars::cons], derived::erf_dersoundspeed);
+        calculate_derived("temp",        vars_new[lev][Vars::cons], derived::erf_dertemp);
+        calculate_derived("theta",       vars_new[lev][Vars::cons], derived::erf_dertheta);
+        calculate_derived("KE",          vars_new[lev][Vars::cons], derived::erf_derKE);
+        calculate_derived("QKE",         vars_new[lev][Vars::cons], derived::erf_derQKE);
+        calculate_derived("scalar",      vars_new[lev][Vars::cons], derived::erf_derscalar);
+        calculate_derived("vorticity_x", mf_cc_vel[lev]           , derived::erf_dervortx);
+        calculate_derived("vorticity_y", mf_cc_vel[lev]           , derived::erf_dervorty);
+        calculate_derived("vorticity_z", mf_cc_vel[lev]           , derived::erf_dervortz);
 
         MultiFab r_hse(base_state[lev], make_alias, 0, 1); // r_0 is first  component
         MultiFab p_hse(base_state[lev], make_alias, 1, 1); // p_0 is second component

Tests/ERFGoldFiles/IsentropicVortexAdvecting/Header

@@ -1,11 +1,12 @@
 HyperCLaw-V1.1
-7
+8
 density
 x_velocity
 y_velocity
 z_velocity
 temp
 theta
+vorticity_z
 pressure
 3
 0.005000000000000001
@@ -18,18 +19,12 @@ pressure
 0.5 0.5 0.5 
 0
 0
-0 4 0.005000000000000001
+0 2 0.005000000000000001
 10
--12 0
--12 0
--1 1
-0 12
+-12 12
 -12 0
 -1 1
--12 0
-0 12
--1 1
-0 12
+-12 12
 0 12
 -1 1
 Level_0/Cell

Tests/ERFGoldFiles/IsentropicVortexAdvecting/Level_0/Cell_H

@@ -1,28 +1,20 @@
 1
 1
-7
+8
 0
-(4 0
-((0,0,0) (23,23,3) (0,0,0))
-((24,0,0) (47,23,3) (0,0,0))
-((0,24,0) (23,47,3) (0,0,0))
-((24,24,0) (47,47,3) (0,0,0))
+(2 0
+((0,0,0) (47,23,3) (0,0,0))
+((0,24,0) (47,47,3) (0,0,0))
 )
-4
+2
 FabOnDisk: Cell_D_00000 0
 FabOnDisk: Cell_D_00001 0
-FabOnDisk: Cell_D_00002 0
-FabOnDisk: Cell_D_00003 0
 
-4,7
-1.1508058595959001e+00,2.9304936277992795e+02,2.5273011926225863e+02,0.0000000000000000e+00,2.9889822852087832e+02,2.9999999999999983e+02,9.8720490014937750e+04,
-1.0999792802791339e+00,2.9312820856682879e+02,2.3942911725482568e+02,0.0000000000000000e+00,2.9354610160293907e+02,2.9999999999999983e+02,9.2670758686574627e+04,
-1.0989693533104781e+00,2.4281456739867247e+02,1.3971339103521245e+02,0.0000000000000000e+00,2.9343826618210880e+02,2.9999999999999977e+02,9.2551662885714919e+04,
-6.7266824942273540e-01,9.8982769986740948e+01,8.2916282642571332e+01,0.0000000000000000e+00,2.4112535881735596e+02,2.9999999999999977e+02,4.6550646053037635e+04,
+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,
 
-4,7
-1.1658900337646418e+00,3.3570828211880183e+02,2.9360352434648462e+02,0.0000000000000000e+00,3.0045923092147848e+02,3.0000000000000017e+02,1.0053679536769103e+05,
-1.1654130664023867e+00,4.4005150729001542e+02,3.5146595531093556e+02,0.0000000000000000e+00,3.0041005755694704e+02,3.0000000000000023e+02,1.0047921842404548e+05,
-1.1656654423979138e+00,3.5574923061766560e+02,2.9371773867629133e+02,0.0000000000000000e+00,3.0043607798337888e+02,3.0000000000000023e+02,1.0050968272762457e+05,
-1.1631913051285259e+00,4.9220671551520189e+02,4.7708174379805695e+02,0.0000000000000000e+00,3.0018084387639641e+02,3.0000000000000017e+02,1.0021114355062916e+05,
+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,
 

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
+erf.plot_vars_1     = density x_velocity y_velocity z_velocity pressure theta temp vorticity_z
 
 # SOLVER CHOICE
 erf.alpha_T = 0.0