Crude balance with moisture and preliminary run.

Exec/RegTests/Bubble/inputs_kessler_bubble

@@ -14,11 +14,8 @@ zlo.type = "SlipWall"
 zhi.type = "SlipWall"
 
 # TIME STEP CONTROL
-#erf.cfl = 2.0
 erf.fixed_dt = 0.5
 erf.fixed_mri_dt_ratio = 4
-#erf.fixed_dt       = 3.0      # fixed time step [s] -- Straka et al 1993
-#erf.fixed_fast_dt  = 0.75     # fixed time step [s] -- Straka et al 1993
 
 # DIAGNOSTICS & VERBOSITY
 erf.sum_interval   = 1       # timesteps between computing mass
@@ -34,8 +31,8 @@ erf.check_int       = 100       # number of timesteps between checkpoints
 
 # PLOTFILES
 erf.plot_file_1     = plt        # prefix of plotfile name
-erf.plot_int_1      = 100          # number of timesteps between plotfiles
-erf.plot_vars_1     = density rhotheta rhoQt rhoQp x_velocity y_velocity z_velocity pressure theta temp pres_hse dens_hse pert_pres pert_dens qt qp qv qc qi 
+erf.plot_int_1      = 100        # number of timesteps between plotfiles
+erf.plot_vars_1     = density rhotheta rhoQt rhoQp rhoadv_0 x_velocity y_velocity z_velocity pressure theta scalar temp pres_hse dens_hse pert_pres pert_dens qt qp qv qc qi 
 
 # SOLVER CHOICES
 erf.use_gravity          = true
@@ -55,6 +52,7 @@ erf.moistscal_vert_adv_type  = "Upwind_3rd"
 erf.les_type        = "None"
 erf.pbl_type        = "None"
 erf.moisture_model  = "Kessler"
+erf.use_moist_background = true
 
 erf.molec_diff_type  = "ConstantAlpha"
 erf.rho0_trans       = 1.0 # [kg/m^3], used to convert input diffusivities
@@ -73,5 +71,5 @@ prob.z_c    =  2000.0
 prob.x_r    =  2000.0
 prob.z_r    =  2000.0
 prob.T_0    =   300.0
-prob.Qt_0   =       0.005
+prob.Qt_0   =       0.025
 prob.T_pert_is_airtemp = false # Perturb theta

Exec/RegTests/Bubble/prob.H

@@ -42,7 +42,7 @@ class Problem : public ProblemBase
 public:
     Problem();
 
-#include "Prob/init_density_hse_dry.H"
+#include "Prob/init_density_hse_moist.H"
 
     void init_custom_pert (
         const amrex::Box&  bx,
@@ -55,6 +55,7 @@ public:
         amrex::Array4<amrex::Real      > const& z_vel,
         amrex::Array4<amrex::Real      > const& r_hse,
         amrex::Array4<amrex::Real      > const& p_hse,
+        amrex::Array4<amrex::Real      > const& qv_arr,
         amrex::Array4<amrex::Real const> const& z_nd,
         amrex::Array4<amrex::Real const> const& z_cc,
 #if defined(ERF_USE_MOISTURE)
@@ -96,11 +97,15 @@ private:
  */
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 void
-perturb_rho_theta (const amrex::Real x,
+perturb_rho_theta (int i,
+                   int j,
+                   int k,
+                   const amrex::Real x,
                    const amrex::Real y,
                    const amrex::Real z,
                    const amrex::Real p_hse,
-                   const amrex::Real r_hse,
+                   const amrex::Real r_hse_dry,
+                   const amrex::Real qv,
                    const ProbParm pp,
                    const amrex::Real rdOcp,
                    amrex::Real& rho,
@@ -120,12 +125,11 @@ perturb_rho_theta (const amrex::Real x,
         dT = 0.0;
     }
     else {
-      //dT = pp.T_pert * (std::cos(PI*L) + 1.0)/2.0;
       dT = pp.T_pert * std::cos(PI*L/2.0)*std::cos(PI*L/2.0);
     }
 
     // Temperature that satisfies the EOS given the hydrostatically balanced (r,p)
-    const amrex::Real Tbar_hse = p_hse / (R_d * r_hse);
+    const amrex::Real Tbar_hse = p_hse / (r_hse_dry * (R_d + R_v*qv));
 
     // Note: theta_perturbed is theta PLUS perturbation in theta
     amrex::Real theta_perturbed;
@@ -143,13 +147,23 @@ perturb_rho_theta (const amrex::Real x,
         // - hydrostatic rebalance is needed (TODO: is this true?)
         // - this is the approach taken in the density current problem
         rhotheta = 0.0; // i.e., hydrostatically balanced pressure stays const
-        rho = getRhoThetagivenP(p_hse) / theta_perturbed - r_hse;
+        rho = getRhoThetagivenP(p_hse, qv) / theta_perturbed - r_hse_dry;
+
+        if (i==128 && j==1)
+            amrex::Print() << "IC PERT: " << amrex::IntVect(i,j,k) << ' '
+                           << Tbar_hse*std::pow(p_0/p_hse, rdOcp) << ' '
+                           << theta_perturbed << ' '
+                           << dT << ' '
+                           << getRhoThetagivenP(p_hse, qv) << ' '
+                           << r_hse_dry << ' '
+                           << p_hse << ' '
+                           << qv << "\n";
     }
     else
     {
         // this version perturbs rho*theta (i.e., p) but not rho
         rho = 0.0; // i.e., hydrostatically balanced density stays const
-        rhotheta = r_hse * theta_perturbed - getRhoThetagivenP(p_hse);
+        rhotheta = r_hse_dry * theta_perturbed - getRhoThetagivenP(p_hse, qv);
     }
 }
 

Exec/RegTests/Bubble/prob.cpp

@@ -12,6 +12,7 @@ amrex_probinit(
 
 Problem::Problem()
 {
+    amrex::Print() << "INIT PARMS\n";
   // Parse params
   ParmParse pp("prob");
   pp.query("T_0", parms.T_0);
@@ -31,7 +32,7 @@ Problem::Problem()
   pp.query("dampcoef", parms.dampcoef);
   pp.query("zdamp", parms.zdamp);
 
-  init_base_parms(parms.rho_0, parms.T_0);
+  init_base_parms(parms.rho_0, parms.T_0, parms.Qt_0);
 }
 
 void
@@ -46,6 +47,7 @@ Problem::init_custom_pert(
     Array4<Real      > const& z_vel,
     Array4<Real      > const& r_hse,
     Array4<Real      > const& p_hse,
+    Array4<Real      > const& qv_arr,
     Array4<Real const> const& /*z_nd*/,
     Array4<Real const> const& z_cc,
 #if defined(ERF_USE_MOISTURE)
@@ -72,15 +74,6 @@ Problem::init_custom_pert(
     const Real prob_lo_z = geomdata.ProbLo()[2];
     const Real rdOcp     = sc.rdOcp;
 
-#if 0
-    // These are at cell centers (unstaggered)
-    Vector<Real> h_r(khi+1);
-    Vector<Real> h_p(khi+1);
-
-    amrex::Gpu::DeviceVector<Real> d_r(khi+1);
-    amrex::Gpu::DeviceVector<Real> d_p(khi+1);
-#endif
-
     amrex::Print() << "Bubble delta T = " << parms.T_pert << " K" << std::endl;
     amrex::Print() << "  centered at ("
         << parms.x_c << " " << parms.y_c << " " << parms.z_c << ")" << std::endl;
@@ -96,30 +89,6 @@ Problem::init_custom_pert(
 
     if (z_cc) { // nonflat terrain
 
-#if 0
-        if (parms.T_0 > 0)
-        {
-            // Create a flat box with same horizontal extent but only one cell in vertical
-            Box b2d = surroundingNodes(bx); // Copy constructor
-            b2d.setRange(2,0);
-
-            ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int)
-            {
-                Array1D<Real,0,255> r;;
-                Array1D<Real,0,255> p;;
-
-                init_isentropic_hse_terrain(i,j,rho_sfc,thetabar,&(r(0)),&(p(0)),z_cc,khi);
-
-                for (int k = 0; k <= khi; k++) {
-                   r_hse(i,j,k) = r(k);
-                   p_hse(i,j,k) = p(k);
-                }
-                r_hse(i,j,   -1) = r_hse(i,j,0);
-                r_hse(i,j,khi+1) = r_hse(i,j,khi);
-            });
-        }
-#endif
-
         amrex::ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
         {
             // Geometry (note we must include these here to get the data on device)
@@ -130,7 +99,9 @@ Problem::init_custom_pert(
             const Real y = prob_lo[1] + (j + 0.5) * dx[1];
             const Real z = z_cc(i,j,k);
 
-            perturb_rho_theta(x, y, z, p_hse(i,j,k), r_hse(i,j,k),
+
+            perturb_rho_theta(i, j, k, x, y, z,
+                              p_hse(i,j,k), r_hse(i,j,k), 0.0,
                               parms, rdOcp,
                               state(i, j, k, Rho_comp),
                               state(i, j, k, RhoTheta_comp));
@@ -146,52 +117,44 @@ Problem::init_custom_pert(
 #endif
 
         });
-    } else {
 
-#if 0
-        if (parms.T_0 > 0)
-        {
-            init_isentropic_hse(rho_sfc,thetabar,h_r.data(),h_p.data(),dz,prob_lo_z,khi);
-
-            amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_r.begin(), h_r.end(), d_r.begin());
-            amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_p.begin(), h_p.end(), d_p.begin());
-
-            Real* r = d_r.data();
-            Real* p = d_p.data();
-
-            amrex::ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int) noexcept
-            {
-                for (int k = 0; k <= khi; k++) {
-                   r_hse(i,j,k) = r[k];
-                   p_hse(i,j,k) = p[k];
-                }
-                r_hse(i,j,   -1) = r_hse(i,j,0);
-                r_hse(i,j,khi+1) = r_hse(i,j,khi);
-            });
-        }
-#endif
+    } else {
 
         const Real Qt_0 = parms.Qt_0;
-
         amrex::ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
         {
             // Geometry (note we must include these here to get the data on device)
-            const auto prob_lo         = geomdata.ProbLo();
-            const auto dx              = geomdata.CellSize();
+            const auto prob_lo = geomdata.ProbLo();
+            const auto dx      = geomdata.CellSize();
 
             const Real x = prob_lo[0] + (i + 0.5) * dx[0];
             const Real y = prob_lo[1] + (j + 0.5) * dx[1];
             const Real z = prob_lo[2] + (k + 0.5) * dx[2];
 
-            perturb_rho_theta(x, y, z, p_hse(i,j,k), r_hse(i,j,k),
+            Real r_hse_dry = r_hse(i,j,k) / (1.0 + Qt_0);
+
+            perturb_rho_theta(i, j, k, x, y, z,
+                              p_hse(i,j,k), r_hse_dry, qv_arr(i,j,k),
                               parms, rdOcp,
                               state(i, j, k, Rho_comp),
                               state(i, j, k, RhoTheta_comp));
 
-            state(i, j, k, RhoScalar_comp) = 0.0;
+            amrex::Real L = 0.0;
+            if (parms.x_r > 0) L += std::pow((x - parms.x_c)/parms.x_r, 2);
+            if (parms.y_r > 0) L += std::pow((y - parms.y_c)/parms.y_r, 2);
+            if (parms.z_r > 0) L += std::pow((z - parms.z_c)/parms.z_r, 2);
+            L = std::sqrt(L);
+            amrex::Real dT;
+            if (L > 1.0) {
+                dT = 0.0;
+            }
+            else {
+                dT = 1.0* std::cos(PI*L/2.0)*std::cos(PI*L/2.0);
+            }
+            state(i, j, k, RhoScalar_comp) = dT; //0.0;
 
 #ifdef ERF_USE_MOISTURE
-            state(i, j, k, RhoQt_comp) = r_hse(i,j,k)*Qt_0;
+            state(i, j, k, RhoQt_comp) = r_hse_dry*Qt_0;
             state(i, j, k, RhoQp_comp) = 0.0;
 #endif
         });

Exec/SquallLine_2D/prob.H

@@ -44,6 +44,7 @@ public:
         amrex::Array4<amrex::Real      > const& z_vel,
         amrex::Array4<amrex::Real      > const& r_hse,
         amrex::Array4<amrex::Real      > const& p_hse,
+        amrex::Array4<amrex::Real      > const& /*qv_arr*/,
         amrex::Array4<amrex::Real const> const& z_nd,
         amrex::Array4<amrex::Real const> const& z_cc,
 #if defined(ERF_USE_MOISTURE)
@@ -62,6 +63,7 @@ public:
 
 
     void erf_init_dens_hse_moist (amrex::MultiFab& rho_hse,
+                                  amrex::MultiFab& /*qmoist*/,
                                   std::unique_ptr<amrex::MultiFab>& z_phys_nd,
                                   std::unique_ptr<amrex::MultiFab>& z_phys_cc,
                                   amrex::Geometry const& geom) override;

Exec/SquallLine_2D/prob.cpp

@@ -188,6 +188,7 @@ init_isentropic_hse_no_terrain(Real *theta, Real* r, Real* p, Real *q_v,
 
 void
 Problem::erf_init_dens_hse_moist (MultiFab& rho_hse,
+                                  amrex::MultiFab& /*qmoist*/,
                                   std::unique_ptr<MultiFab>& z_phys_nd,
                                   std::unique_ptr<MultiFab>& z_phys_cc,
                                   Geometry const& geom)
@@ -283,6 +284,7 @@ Problem::init_custom_pert (
     Array4<Real      > const& z_vel,
     Array4<Real      > const& /*r_hse*/,
     Array4<Real      > const& /*p_hse*/,
+    amrex::Array4<amrex::Real      > const& /*qv_arr*/,
     Array4<Real const> const& /*z_nd*/,
     Array4<Real const> const& /*z_cc*/,
 #if defined(ERF_USE_MOISTURE)

Exec/SuperCell/prob.H

@@ -38,6 +38,7 @@ public:
         amrex::Array4<amrex::Real      > const& z_vel,
         amrex::Array4<amrex::Real      > const& r_hse,
         amrex::Array4<amrex::Real      > const& p_hse,
+        amrex::Array4<amrex::Real      > const& /*qv_arr*/,
         amrex::Array4<amrex::Real const> const& z_nd,
         amrex::Array4<amrex::Real const> const& z_cc,
 #if defined(ERF_USE_MOISTURE)

Exec/SuperCell/prob.cpp

@@ -94,6 +94,7 @@ Problem::init_custom_pert(
     Array4<Real      > const& z_vel,
     Array4<Real      > const& r_hse,
     Array4<Real      > const& p_hse,
+    amrex::Array4<amrex::Real      > const& /*qv_arr*/,
     Array4<Real const> const& /*z_nd*/,
     Array4<Real const> const& /*z_cc*/,
 #if defined(ERF_USE_MOISTURE)

Source/ERF.H

@@ -366,6 +366,8 @@ private:
     // Compute a new MultiFab by copying from valid region and filling ghost cells -
     // works for single level and 2-level cases (fill fine grid ghost by interpolating from coarse)
     void FillPatchMoistVars (int lev, amrex::MultiFab& mf);
+
+    void init_from_hse_moist (int lev);
 #endif
 
     // compute new multifabs by copying in data from valid region and filling ghost cells

Source/ERF.cpp

@@ -856,6 +856,10 @@ ERF::init_only (int lev, Real time)
     lev_new[Vars::yvel].setVal(0.0); lev_old[Vars::yvel].setVal(0.0);
     lev_new[Vars::zvel].setVal(0.0); lev_old[Vars::zvel].setVal(0.0);
 
+#if defined(ERF_USE_MOISTURE)
+    qmoist[lev].setVal(0.);
+#endif
+
     // Initialize background flow (optional)
     if (init_type == "input_sounding") {
         // The base state is initialized by integrating vertically through the
@@ -886,12 +890,12 @@ ERF::init_only (int lev, Real time)
         // background field to be equal to the base state, calculated from the
         // problem-specific erf_init_dens_hse
         initHSE(lev); // need to call this first
-        init_from_hse(lev);
-    }
-
 #if defined(ERF_USE_MOISTURE)
-    qmoist[lev].setVal(0.);
+        init_from_hse_moist(lev);
+#else
+        init_from_hse(lev);
 #endif
+    }
 
     // Add problem-specific flow features
     //