Merge pull request #147 from E3SM-Project/celdred/cfl-w-limit

adding W clipping and CRM dt adjustment based on vertical cfl

dynamics/spam/Dycore.h

@@ -245,13 +245,27 @@ class Dycore {
         coupler, prognostic_vars, auxiliary_vars, constant_vars,
         params.couple_wind, params.couple_wind_exact_inverse);
 
-    // Time stepping loop
+    // adjust crm_per_phys and dtcrm if needed
+
+    if (params.adjust_crm_per_phys_using_vert_cfl) {
+      tendencies.adjust_crm_per_phys_using_vert_cfl(params, prognostic_vars);
+    }
+    // tendencies.a
+    //  Time stepping loop
     debug_print("start time stepping loop", par.masterproc);
     for (uint nstep = 0; nstep < params.crm_per_phys; nstep++) {
       yakl::fence();
 
       time_integrator->step_forward(params.dtcrm);
 
+      if (params.clip_vertical_velocities) {
+        tendencies.clip_vertical_velocities(params, prognostic_vars);
+      }
+
+      if (params.clip_negative_densities) {
+        tendencies.clip_negative_densities(prognostic_vars);
+      }
+
 #if defined PAMC_AN || defined PAMC_MAN
       if (params.check_anelastic_constraint) {
         real max_div = tendencies.compute_max_anelastic_constraint(
@@ -287,10 +301,6 @@ class Dycore {
     }
     prevstep += params.crm_per_phys;
 
-    if (params.clip_negative_densities) {
-      tendencies.clip_negative_densities(prognostic_vars);
-    }
-
     // convert dynamics state to Coupler state
     tendencies.convert_dynamics_to_coupler_state(
         coupler, prognostic_vars, constant_vars, params.couple_wind,

dynamics/spam/src/core/params.h

@@ -42,6 +42,11 @@ class Parameters {
   // during conversion to coupler state when coupling winds
   bool couple_wind_exact_inverse = false;
   bool clip_negative_densities = true;
+
+  bool clip_vertical_velocities = false;
+  bool adjust_crm_per_phys_using_vert_cfl = false;
+  real target_cfl = 0.7;
+  real max_w = 50.0;
 };
 
 void readParamsFile(std::string inFile, Parameters &params, Parallel &par,
@@ -82,6 +87,15 @@ void readParamsFile(std::string inFile, Parameters &params, Parallel &par,
   params.clip_negative_densities =
       config["clip_negative_densities"].as<bool>(true);
 
+  params.clip_vertical_velocities =
+      config["clip_vertical_velocities"].as<bool>(false);
+
+  params.adjust_crm_per_phys_using_vert_cfl =
+      config["adjust_crm_per_phys_using_vert_cfl"].as<bool>(false);
+
+  params.target_cfl = config["max_w"].as<real>(0.7);
+  params.max_w = config["max_w"].as<real>(50.0);
+
   // ADD A CHECK HERE THAT TOTAL TIME IS EXACTLY DIVISIBLE BY STAT_FREQ
   if (params.stat_freq >= 0.) {
     params.statSize = params.sim_time / params.stat_freq;

dynamics/spam/src/models/extrudedmodel.h

@@ -244,6 +244,71 @@ class ModelTendencies : public ExtrudedTendencies {
     check_anelastic_constraint = params.check_anelastic_constraint;
   }
 
+  void clip_vertical_velocities(ModelParameters &params,
+                                FieldSet<nprognostic> &x) {
+    real max_w = params.max_w;
+
+    const auto &primal_topology = primal_geometry.topology;
+
+    const auto wvar = x.fields_arr[WVAR].data;
+    const int pis = primal_topology.is;
+    const int pjs = primal_topology.js;
+    const int pks = primal_topology.ks;
+
+    parallel_for(
+        "Clip vertical velocity",
+        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) {
+          real wedgelength =
+              primal_geometry.get_area_01entity(k + pks, j + pjs, i + pis, n);
+          wvar(0, k + pks, j + pjs, i + pis, n) =
+              std::max(-max_w * wedgelength,
+                       std::min(max_w * wedgelength,
+                                wvar(0, k + pks, j + pjs, i + pis, n)));
+        });
+  }
+
+  void adjust_crm_per_phys_using_vert_cfl(ModelParameters &params,
+                                          FieldSet<nprognostic> &x) {
+    // cfl = c * dt / dx
+    // dt = cfl * dx / c
+
+    const auto &primal_topology = primal_geometry.topology;
+    const auto wvar = x.fields_arr[WVAR].data;
+    const int pis = primal_topology.is;
+    const int pjs = primal_topology.js;
+    const int pks = primal_topology.ks;
+
+    real4d min_dt_4d("min_dt_4d", primal_topology.nl, primal_topology.n_cells_y,
+                     primal_topology.n_cells_x, primal_topology.nens);
+
+    real target_cfl = params.target_cfl;
+    real eps = 1e-8;
+    parallel_for(
+        "compute adjusted dt",
+        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) {
+          real wedgelength =
+              primal_geometry.get_area_01entity(k + pks, j + pjs, i + pis, n);
+          real wval =
+              std::abs(wvar(0, k + pks, j + pjs, i + pis, n)) / wedgelength;
+          min_dt_4d(k, j, i, n) = target_cfl * wedgelength / (wval + eps);
+        });
+    real min_dt = yakl::intrinsics::minval(min_dt_4d);
+    real adj_min_dt = std::min(params.dtcrm, min_dt);
+
+    params.crm_per_phys = (int)std::ceil(params.dt_crm_phys / adj_min_dt);
+    params.dtcrm = params.dt_crm_phys / params.crm_per_phys;
+
+    //  std::cout << "min_dt:         " << min_dt << "\n";
+    //  std::cout << "adj_min_dt:         " << adj_min_dt << "\n";
+    //  std::cout << "dtcrm:         " << params.dtcrm << "\n";
+    //  std::cout << "dt_crm_phys:         " << params.dt_crm_phys << "\n";
+    //  std::cout << "crm per phys:     " << params.crm_per_phys << "\n";
+  }
+
   void convert_dynamics_to_coupler_state(PamCoupler &coupler,
                                          const FieldSet<nprognostic> &prog_vars,
                                          const FieldSet<nconstant> &const_vars,

dynamics/spam/src/models/layermodel.h

@@ -131,6 +131,12 @@ class ModelTendencies : public Tendencies {
     }
   }
 
+  void clip_vertical_velocities(ModelParameters &params,
+                                FieldSet<nprognostic> &x) {}
+
+  void adjust_crm_per_phys_using_vert_cfl(ModelParameters &params,
+                                          FieldSet<nprognostic> &x) {}
+
   void compute_constants(FieldSet<nconstant> &const_vars,
                          FieldSet<nprognostic> &x) override {}