Merge pull request #137 from E3SM-Project/celdred/horiz-vert-diffusio…

…n-split

Split dissipation into horizontal + vertical parts, also split velocity dissipation into vorticity + divergence

dynamics/spam/CMakeLists.txt

@@ -14,7 +14,6 @@ target_include_directories(dycore INTERFACE src/operators)
 target_include_directories(dycore INTERFACE src/timesteppers)
 target_include_directories(dycore INTERFACE src/hamiltonians)
 target_include_directories(dycore INTERFACE src/grids)
-target_include_directories(dycore INTERFACE src/common)
 target_include_directories(dycore INTERFACE src/io)
 target_include_directories(dycore INTERFACE src/core)
 
@@ -84,5 +83,3 @@ endif()
 if ("${PAMC_IO}" STREQUAL "none")
   target_compile_definitions(dycore INTERFACE -DPAMC_NOIO)
 endif()
-
-

dynamics/spam/src/extrudedmodel-common.h

@@ -101,9 +101,14 @@ class ModelParameters : public Parameters {
   bool dycore_tracer_pos[ntracers_dycore + GPU_PAD];
   bool acoustic_balance;
   bool uniform_vertical;
-  real scalar_diffusion_coeff;
-  real scalar_diffusion_subtract_refstate;
-  real velocity_diffusion_coeff;
+  real scalar_horiz_diffusion_coeff;
+  real scalar_vert_diffusion_coeff;
+  real velocity_vort_horiz_diffusion_coeff;
+  real velocity_vort_vert_diffusion_coeff;
+  real velocity_div_horiz_diffusion_coeff;
+  real velocity_div_vert_diffusion_coeff;
+  bool scalar_diffusion_subtract_refstate;
+
   // forces reference state to be in perfect hydrostatic balance by subtracting
   // the hydrostatic balance equation evaluated at the reference state in
   // the velocity tendency

dynamics/spam/src/models/extrudedmodel.h

@@ -349,9 +349,13 @@ struct TotalDensityFunctor {
 // *******   Tendencies   ***********//
 
 class ModelTendencies : public ExtrudedTendencies {
-  real scalar_diffusion_coeff;
+  real scalar_horiz_diffusion_coeff;
+  real scalar_vert_diffusion_coeff;
   real scalar_diffusion_subtract_refstate;
-  real velocity_diffusion_coeff;
+  real velocity_vort_horiz_diffusion_coeff;
+  real velocity_vort_vert_diffusion_coeff;
+  real velocity_div_horiz_diffusion_coeff;
+  real velocity_div_vert_diffusion_coeff;
   bool force_refstate_hydrostatic_balance;
   bool check_anelastic_constraint;
 #if defined PAMC_AN || defined PAMC_MAN
@@ -366,12 +370,15 @@ class ModelTendencies : public ExtrudedTendencies {
                   const Geometry<Twisted> &dual_geom) {
 
     ExtrudedTendencies::initialize(params, equations, primal_geom, dual_geom);
-    scalar_diffusion_coeff = params.scalar_diffusion_coeff;
-    scalar_diffusion_subtract_refstate =
-        params.scalar_diffusion_subtract_refstate;
-    velocity_diffusion_coeff = params.velocity_diffusion_coeff;
-    force_refstate_hydrostatic_balance =
-        params.force_refstate_hydrostatic_balance;
+    scalar_horiz_diffusion_coeff = params.scalar_horiz_diffusion_coeff;
+    scalar_vert_diffusion_coeff = params.scalar_vert_diffusion_coeff;
+    velocity_vort_horiz_diffusion_coeff = params.velocity_vort_horiz_diffusion_coeff;
+    velocity_vort_vert_diffusion_coeff = params.velocity_vort_vert_diffusion_coeff;
+    velocity_div_horiz_diffusion_coeff = params.velocity_div_horiz_diffusion_coeff;
+    velocity_div_vert_diffusion_coeff = params.velocity_div_vert_diffusion_coeff;
+
+    scalar_diffusion_subtract_refstate = params.scalar_diffusion_subtract_refstate;
+    force_refstate_hydrostatic_balance = params.force_refstate_hydrostatic_balance;
     check_anelastic_constraint = params.check_anelastic_constraint;
 
 #if defined PAMC_AN || defined PAMC_MAN
@@ -1272,7 +1279,7 @@ class ModelTendencies : public ExtrudedTendencies {
     yakl::timer_stop("compute_recons");
   }
 
-  void add_scalar_diffusion(real scalar_diffusion_coeff, real5d denstendvar,
+  void add_scalar_diffusion(real scalar_horiz_diffusion_coeff, real scalar_vert_diffusion_coeff, real5d denstendvar,
                             const real5d densvar, const real5d dens0var,
                             const real5d Fdiffvar, const real5d FWdiffvar,
                             FieldSet<nauxiliary> &auxiliary_vars) {
@@ -1373,8 +1380,7 @@ class ModelTendencies : public ExtrudedTendencies {
               primal_geometry, dual_geometry, pis, pjs, pks, i, j, k, n);
 
           for (int d = 0; d < VS::ndensity_diffused; ++d) {
-            const real diff_tend = -scalar_diffusion_coeff * rho * Hn1bar_diag *
-                                   (hdiv(d) + vdiv(d));
+            const real diff_tend = -scalar_horiz_diffusion_coeff * rho * Hn1bar_diag * hdiv(d) + -scalar_vert_diffusion_coeff * rho * Hn1bar_diag * vdiv(d);
             int dens_id = varset.diffused_dens_ids(d);
             denstendvar(dens_id, k + pks, j + pjs, i + pis, n) += diff_tend;
           }
@@ -1384,7 +1390,8 @@ class ModelTendencies : public ExtrudedTendencies {
   }
 
   void add_velocity_diffusion_2d(
-      real velocity_coeff, real5d Vtendvar, real5d Wtendvar, const real5d Vvar,
+      real velocity_vort_horiz_diffusion_coeff, real velocity_vort_vert_diffusion_coeff, real velocity_div_horiz_diffusion_coeff, real velocity_div_vert_diffusion_coeff,
+      real5d Vtendvar, real5d Wtendvar, const real5d Vvar,
       const real5d Wvar, const real5d qhzedgereconvar, const real5d qhzvar,
       const real5d dens0var, const real5d Kvar, const real5d Fvar,
       const real5d FWvar, FieldSet<nauxiliary> &auxiliary_vars) {
@@ -1465,21 +1472,21 @@ class ModelTendencies : public ExtrudedTendencies {
         SimpleBounds<4>(primal_topology.ni, primal_topology.n_cells_y,
                         primal_topology.n_cells_x, primal_topology.nens),
         YAKL_LAMBDA(int k, int j, int i, int n) {
-          // *d(*d)
+          // *d(*d) = vort
           const real Dvert =
               compute_D0bar_vert<1>(qhzvar, dis, djs, dks, i, j, k, n);
           SArray<real, 1, ndims> Hnm11bar_diag;
           Hnm11bar_diagonal(Hnm11bar_diag, primal_geometry, dual_geometry, dis,
                             djs, dks, i, j, k, n);
           Vtendvar(0, k + pks, j + pjs, i + pis, n) -=
-              velocity_coeff * Dvert * Hnm11bar_diag(0);
+              velocity_vort_horiz_diffusion_coeff * Dvert * Hnm11bar_diag(0);
 
-          // d(*d*)
+          // d(*d*) = div
           SArray<real, 2, 1, ndims> vdiff;
           compute_D0<1>(vdiff, dens0var, pis, pjs, pks, i, j, k, n);
           for (int d = 0; d < ndims; ++d) {
             Vtendvar(d, pks + k, pjs + j, pis + i, n) -=
-                velocity_coeff * vdiff(0, d);
+                velocity_div_horiz_diffusion_coeff * vdiff(0, d);
           }
         });
 
@@ -1488,26 +1495,27 @@ class ModelTendencies : public ExtrudedTendencies {
         SimpleBounds<4>(primal_topology.nl, primal_topology.n_cells_y,
                         primal_topology.n_cells_x, primal_topology.nens),
         YAKL_LAMBDA(int k, int j, int i, int n) {
-          // *d(*d)
+          // *d(*d) = vort
           SArray<real, 1, ndims> Dhorz;
           compute_D0bar_ext<1>(Dhorz, qhzvar, dis, djs, dks, i, j, k + 1, n);
           const real Hn0bar_diag = Hn0bar_diagonal(
               primal_geometry, dual_geometry, dis, djs, dks, i, j, k, n);
           Wtendvar(0, k + pks, j + pjs, i + pis, n) -=
-              velocity_coeff * Dhorz(0) * Hn0bar_diag;
+              velocity_vort_vert_diffusion_coeff * Dhorz(0) * Hn0bar_diag;
 
-          // d(*d*)
+          // d(*d*) = div
           SArray<real, 1, 1> wdiff;
           compute_D0_vert<1>(wdiff, dens0var, pis, pjs, pks, i, j, k, n);
           Wtendvar(0, pks + k, pjs + j, pis + i, n) -=
-              velocity_coeff * wdiff(0);
+              velocity_div_vert_diffusion_coeff * wdiff(0);
         });
 
     yakl::timer_stop("add_velocity_diffusion");
   }
 
   void add_velocity_diffusion_3d(
-      real velocity_coeff, real5d Vtendvar, real5d Wtendvar, const real5d Vvar,
+      real velocity_vort_horiz_diffusion_coeff, real velocity_vort_vert_diffusion_coeff, real velocity_div_horiz_diffusion_coeff, real velocity_div_vert_diffusion_coeff,
+      real5d Vtendvar, real5d Wtendvar, const real5d Vvar,
       const real5d Wvar, const real5d qhzedgereconvar, const real5d qhzvar,
       const real5d dens0var, const real5d Kvar, const real5d Fvar,
       const real5d FWvar, const real5d qxyedgereconvar, const real5d qxyvar,
@@ -1529,7 +1537,7 @@ class ModelTendencies : public ExtrudedTendencies {
     YAKL_SCOPE(primal_geometry, this->primal_geometry);
     YAKL_SCOPE(dual_geometry, this->dual_geometry);
 
-    // *d*d
+    // *d*d = vort
     parallel_for(
         "Velocity diffusion 1",
         SimpleBounds<4>(primal_topology.nl, primal_topology.n_cells_y,
@@ -1599,7 +1607,7 @@ class ModelTendencies : public ExtrudedTendencies {
                                                   i, j, k, n);
           for (int d = 0; d < ndims; ++d) {
             Vtendvar(d, k + pks, j + pjs, i + pis, n) +=
-                velocity_coeff * vdiff(d);
+                velocity_vort_horiz_diffusion_coeff * vdiff(d);
           }
         });
 
@@ -1613,10 +1621,10 @@ class ModelTendencies : public ExtrudedTendencies {
               wdiff, FWvar, primal_geometry, dual_geometry, pis, pjs, pks, i, j,
               k, n);
           Wtendvar(0, k + pks, j + pjs, i + pis, n) +=
-              velocity_coeff * wdiff(0);
+              velocity_vort_vert_diffusion_coeff * wdiff(0);
         });
 
-    // d*d*
+    // d*d* = div
     parallel_for(
         "Velocity diffusion 7",
         SimpleBounds<4>(dual_topology.nl, dual_topology.n_cells_y,
@@ -1671,7 +1679,7 @@ class ModelTendencies : public ExtrudedTendencies {
           compute_D0<1>(vdiff, dens0var, pis, pjs, pks, i, j, k, n);
           for (int d = 0; d < ndims; ++d) {
             Vtendvar(d, pks + k, pjs + j, pis + i, n) -=
-                velocity_coeff * vdiff(0, d);
+                velocity_div_horiz_diffusion_coeff * vdiff(0, d);
           }
         });
 
@@ -1683,7 +1691,7 @@ class ModelTendencies : public ExtrudedTendencies {
           SArray<real, 1, 1, 1> wdiff;
           compute_D0_vert<1>(wdiff, dens0var, pis, pjs, pks, i, j, k, n);
           Wtendvar(0, pks + k, pjs + j, pis + i, n) -=
-              velocity_coeff * wdiff(0);
+              velocity_div_vert_diffusion_coeff * wdiff(0);
         });
 
     yakl::timer_stop("add_velocity_diffusion");
@@ -2483,17 +2491,17 @@ class ModelTendencies : public ExtrudedTendencies {
               : std::nullopt);
     }
 
-    if (scalar_diffusion_coeff > 0) {
+    if (scalar_horiz_diffusion_coeff > 0 or scalar_vert_diffusion_coeff > 0) {
       add_scalar_diffusion(
-          scalar_diffusion_coeff, xtend.fields_arr[DENSVAR].data,
+          scalar_horiz_diffusion_coeff, scalar_vert_diffusion_coeff, xtend.fields_arr[DENSVAR].data,
           x.fields_arr[DENSVAR].data, auxiliary_vars.fields_arr[DENS0VAR].data,
           auxiliary_vars.fields_arr[FDIFFVAR].data,
           auxiliary_vars.fields_arr[FWDIFFVAR].data, auxiliary_vars);
     }
-    if (velocity_diffusion_coeff > 0) {
+    if (velocity_vort_horiz_diffusion_coeff > 0 or velocity_vort_vert_diffusion_coeff > 0 or velocity_div_horiz_diffusion_coeff > 0 or velocity_div_vert_diffusion_coeff > 0) {
       if (ndims == 1) {
         add_velocity_diffusion_2d(
-            velocity_diffusion_coeff, xtend.fields_arr[VVAR].data,
+            velocity_vort_horiz_diffusion_coeff, velocity_vort_vert_diffusion_coeff, velocity_div_horiz_diffusion_coeff, velocity_div_vert_diffusion_coeff, xtend.fields_arr[VVAR].data,
             xtend.fields_arr[WVAR].data, x.fields_arr[VVAR].data,
             x.fields_arr[WVAR].data,
             auxiliary_vars.fields_arr[QHZEDGERECONVAR].data,
@@ -2504,7 +2512,7 @@ class ModelTendencies : public ExtrudedTendencies {
             auxiliary_vars.fields_arr[FWVAR].data, auxiliary_vars);
       } else {
         add_velocity_diffusion_3d(
-            velocity_diffusion_coeff, xtend.fields_arr[VVAR].data,
+            velocity_vort_horiz_diffusion_coeff, velocity_vort_vert_diffusion_coeff, velocity_div_horiz_diffusion_coeff, velocity_div_vert_diffusion_coeff, xtend.fields_arr[VVAR].data,
             xtend.fields_arr[WVAR].data, x.fields_arr[VVAR].data,
             x.fields_arr[WVAR].data,
             auxiliary_vars.fields_arr[QHZEDGERECONVAR].data,
@@ -3785,17 +3793,18 @@ void read_model_params_file(std::string inFile, ModelParameters &params,
   params.acoustic_balance = config["balance_initial_density"].as<bool>(false);
   params.uniform_vertical = (config["vcoords"].as<std::string>() == "uniform");
   // Read diffusion coefficients
-  params.scalar_diffusion_coeff = config["scalar_diffusion_coeff"].as<real>(0);
-  params.scalar_diffusion_subtract_refstate =
-      config["scalar_diffusion_subtract_refstate"].as<bool>(true);
-  params.velocity_diffusion_coeff =
-      config["velocity_diffusion_coeff"].as<real>(0);
+  params.scalar_horiz_diffusion_coeff = config["scalar_horiz_diffusion_coeff"].as<real>(0);
+  params.scalar_vert_diffusion_coeff = config["scalar_vert_diffusion_coeff"].as<real>(0);
+  params.velocity_vort_horiz_diffusion_coeff = config["velocity_vort_horiz_diffusion_coeff"].as<real>(0);
+  params.velocity_vort_vert_diffusion_coeff = config["velocity_vort_vert_diffusion_coeff"].as<real>(0);
+  params.velocity_div_horiz_diffusion_coeff = config["velocity_div_horiz_diffusion_coeff"].as<real>(0);
+  params.velocity_div_vert_diffusion_coeff = config["velocity_div_vert_diffusion_coeff"].as<real>(0);
+
+  params.scalar_diffusion_subtract_refstate = config["scalar_diffusion_subtract_refstate"].as<bool>(true);
   // Read the data initialization options
   params.init_data = config["init_data"].as<std::string>();
-  params.force_refstate_hydrostatic_balance =
-      config["force_refstate_hydrostatic_balance"].as<bool>(false);
-  params.check_anelastic_constraint =
-      config["check_anelastic_constraint"].as<bool>(false);
+  params.force_refstate_hydrostatic_balance = config["force_refstate_hydrostatic_balance"].as<bool>(false);
+  params.check_anelastic_constraint = config["check_anelastic_constraint"].as<bool>(false);
 
   for (int i = 0; i < ntracers_dycore; i++) {
     params.init_dycore_tracer[i] =
@@ -3825,9 +3834,13 @@ void read_model_params_coupler(ModelParameters &params, Parallel &par,
 
   params.acoustic_balance = false;
   params.uniform_vertical = false;
-  params.scalar_diffusion_coeff = 0;
+  params.scalar_horiz_diffusion_coeff = 0;
+  params.scalar_vert_diffusion_coeff = 0;
+  params.velocity_vort_horiz_diffusion_coeff = 0;
+  params.velocity_vort_vert_diffusion_coeff = 0;
+  params.velocity_div_horiz_diffusion_coeff = 0;
+  params.velocity_div_vert_diffusion_coeff = 0;
   params.scalar_diffusion_subtract_refstate = true;
-  params.velocity_diffusion_coeff = 0;
   params.init_data = "coupler";
   params.force_refstate_hydrostatic_balance = true;
   params.check_anelastic_constraint = false;
@@ -3853,15 +3866,13 @@ void check_and_print_model_parameters(const ModelParameters &params,
     serial_print("acoustically balanced: " +
                      std::to_string(params.acoustic_balance),
                  par.masterproc);
-    serial_print("scalar_diffusion_coeff: " +
-                     std::to_string(params.scalar_diffusion_coeff),
-                 par.masterproc);
-    serial_print("scalar_diffusion_subtract_refstate: " +
-                     std::to_string(params.scalar_diffusion_subtract_refstate),
-                 par.masterproc);
-    serial_print("velocity_diffusion_coeff: " +
-                     std::to_string(params.velocity_diffusion_coeff),
-                 par.masterproc);
+    serial_print("scalar_horiz_diffusion_coeff: " + std::to_string(params.scalar_horiz_diffusion_coeff), par.masterproc);
+    serial_print("scalar_vert_diffusion_coeff: " + std::to_string(params.scalar_vert_diffusion_coeff), par.masterproc);
+    serial_print("velocity_vort_horiz_diffusion_coeff: " + std::to_string(params.velocity_vort_horiz_diffusion_coeff), par.masterproc);
+    serial_print("velocity_vort_vert_diffusion_coeff: " + std::to_string(params.velocity_vort_vert_diffusion_coeff), par.masterproc);
+    serial_print("velocity_div_horiz_diffusion_coeff: " + std::to_string(params.velocity_div_horiz_diffusion_coeff), par.masterproc);
+    serial_print("velocity_div_vert_diffusion_coeff: " + std::to_string(params.velocity_div_vert_diffusion_coeff), par.masterproc);
+    serial_print("scalar_diffusion_subtract_refstate: " + std::to_string(params.scalar_diffusion_subtract_refstate), par.masterproc);
 
     for (int i = 0; i < ntracers_dycore; i++) {
       serial_print("Dycore Tracer" + std::to_string(i) +

standalone/mmf_simplified/inputs/pamc_idealized/pamc_input_extruded_densitycurrent.yaml

@@ -4,8 +4,12 @@ verbose : true
 # Simulation time in seconds
 sim_time : 900
 
-scalar_diffusion_coeff : 75
-velocity_diffusion_coeff : 75
+scalar_horiz_diffusion_coeff : 75
+scalar_vert_diffusion_coeff : 75
+velocity_vort_horiz_diffusion_coeff : 75
+velocity_vort_vert_diffusion_coeff : 75
+velocity_div_horiz_diffusion_coeff : 75
+velocity_div_vert_diffusion_coeff : 75
 
 # Number of cells to use
 crm_nx : 256

standalone/mmf_simplified/pam-c/run_extruded_script.sh

@@ -1,9 +1,34 @@
+#clean up any existing files
 rm *.png *.nc
 cd ../build
 rm driver
-./../pam-c/make_script.sh $3
-mpirun.mpich -n $3 ./driver ../inputs/pamc_input_extruded_$2.yaml
+
+#build model
+source ../../machines/linux_laptop_gnu_mpi_cpu.env
+./cmakescript_pamc.sh PAM_SGS=none PAM_MICRO=none PAMC_MODEL=extrudedmodel PAMC_HAMIL=an PAMC_THERMO=idealgaspottemp PAMC_IO=serial
+make -j 4
+
+#linux_laptop_gnu_mpi_cpu_debug linux_laptop_gnu_mpi_cpu
+#p3 none
+#shoc none
+#ce mce_rho an man
+#idealgaspottemp constkappavirpottemp
+
+
+#run model
+mpirun.mpich -n $1 ./driver ../inputs/pamc_idealized/pamc_input_extruded_densitycurrent.yaml
 cd ../pam-c
 mv ../build/*.nc .
-python3 plot_extrudedmodel2D.py $1
-#python3 plot_extrudedmodel2D_parallel.py $1 ../inputs/pamc_input_extruded_$2.yaml $3
+
+#pamc_input_extruded_densitycurrent
+#pamc_input_extruded_gravitywave
+#pamc_input_extruded_largerisingbubble
+#pamc_input_extruded_moistrisingbubble
+#pamc_input_extruded_risingbubble
+#pamc_input_extruded_twobubbles
+
+
+#plot model
+python3 plot_extrudedmodel2D.py an
+
+#an ce man mce