Update the divergence routine we call when writing a plotfile, and up…

…date the moving geometry substepping to account for time-varying r0 in the buoyancy term

CMake/BuildERFExe.cmake

@@ -186,10 +186,10 @@ function(build_erf_lib erf_lib_name)
        ${SRC_DIR}/TimeIntegration/ERF_FastRhs_MT.cpp
        ${SRC_DIR}/Utils/ERF_AverageDown.cpp
        ${SRC_DIR}/Utils/ERF_ChopGrids.cpp
+       ${SRC_DIR}/Utils/ERF_ComputeDivergence.cpp 
        ${SRC_DIR}/Utils/ERF_MomentumToVelocity.cpp
        ${SRC_DIR}/LinearSolvers/ERF_PoissonSolve.cpp
        ${SRC_DIR}/LinearSolvers/ERF_PoissonSolve_tb.cpp
-       ${SRC_DIR}/LinearSolvers/ERF_ComputeDivergence.cpp 
        ${SRC_DIR}/LinearSolvers/ERF_SolveWithGMRES.cpp
        ${SRC_DIR}/LinearSolvers/ERF_SolveWithMLMG.cpp
        ${SRC_DIR}/LinearSolvers/ERF_TerrainPoisson.cpp

Source/ERF.H

@@ -164,8 +164,8 @@ public:
     void WriteMyEBSurface ();
 #endif
 
-    // Compute the divergence of the momenta
-    void compute_divergence (int lev, amrex::MultiFab& rhs, amrex::Vector<amrex::MultiFab>& mom_mf,
+    // Compute the divergence -- whether EB, no-terrain, flat terrain or general terrain
+    void compute_divergence (int lev, amrex::MultiFab& rhs, amrex::Array<amrex::MultiFab const*,AMREX_SPACEDIM> rho0_u_const,
                              amrex::Geometry const& geom_at_lev);
 
     // Project the velocities to be divergence-free -- this is only relevant if anelastic == 1
@@ -847,6 +847,7 @@ private:
 
     amrex::Vector<amrex::MultiFab> base_state;
     amrex::Vector<amrex::MultiFab> base_state_new;
+    amrex::Vector<amrex::MultiFab> base_state_src;
 
     // Wave coupling data
     amrex::Vector<std::unique_ptr<amrex::MultiFab>> Hwave;

Source/ERF.cpp

@@ -288,6 +288,7 @@ ERF::ERF_shared ()
     // Base state
     base_state.resize(nlevs_max);
     base_state_new.resize(nlevs_max);
+    base_state_src.resize(nlevs_max);
 
     // Wave coupling data
     Hwave.resize(nlevs_max);
@@ -575,7 +576,8 @@ ERF::post_timestep (int nstep, Real time, Real dt_lev0)
         // Copy z_phs_nd and detJ_cc at end of timestep
         MultiFab::Copy(*z_phys_nd[lev], *z_phys_nd_new[lev], 0, 0, 1, z_phys_nd[lev]->nGrowVect());
         MultiFab::Copy(  *detJ_cc[lev],   *detJ_cc_new[lev], 0, 0, 1,   detJ_cc[lev]->nGrowVect());
-        MultiFab::Copy(base_state[lev],base_state_new[lev],0,0,BaseState::num_comps,base_state[lev].nGrowVect());
+        MultiFab::Copy(base_state[lev]    ,base_state_new[lev],0,0,BaseState::num_comps,base_state[lev].nGrowVect());
+        MultiFab::Copy(base_state_src[lev],base_state_new[lev],0,0,BaseState::num_comps,base_state[lev].nGrowVect());
 
         make_zcc(geom[lev],*z_phys_nd[lev],*z_phys_cc[lev]);
       }
@@ -968,9 +970,10 @@ ERF::InitData_post ()
         // For moving terrain only
         if (solverChoice.terrain_type == TerrainType::Moving) {
             MultiFab::Copy(base_state_new[lev],base_state[lev],0,0,BaseState::num_comps,base_state[lev].nGrowVect());
+            MultiFab::Copy(base_state_src[lev],base_state[lev],0,0,BaseState::num_comps,base_state[lev].nGrowVect());
             base_state_new[lev].FillBoundary(geom[lev].periodicity());
+            base_state_src[lev].FillBoundary(geom[lev].periodicity());
         }
-
     }
 
     // Allow idealized cases over water, used to set lmask

Source/ERF_MakeNewArrays.cpp

@@ -37,6 +37,8 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm,
     if (solverChoice.terrain_type == TerrainType::Moving) {
         base_state_new[lev].define(ba,dm,BaseState::num_comps,base_state[lev].nGrowVect());
         base_state_new[lev].setVal(0.);
+        base_state_src[lev].define(ba,dm,BaseState::num_comps,base_state[lev].nGrowVect());
+        base_state_src[lev].setVal(0.);
     }
 
     // ********************************************************************************************

Source/IO/ERF_Plotfile.cpp

@@ -1,7 +1,6 @@
 #include <ERF_EOS.H>
 #include <ERF.H>
 #include "AMReX_Interp_3D_C.H"
-#include "AMReX_PlotFileUtil.H"
 #include "ERF_TerrainMetrics.H"
 #include "ERF_Constants.H"
 
@@ -368,7 +367,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector<std::string> p
             u[0] = &(vars_new[lev][Vars::xvel]);
             u[1] = &(vars_new[lev][Vars::yvel]);
             u[2] = &(vars_new[lev][Vars::zvel]);
-            computeDivergence(dmf, u, geom[lev]);
+            compute_divergence (lev, dmf, u, geom[lev]);
             mf_comp += 1;
         }
 

Source/LinearSolvers/ERF_PoissonSolve.cpp

@@ -72,7 +72,12 @@ void ERF::project_velocities (int lev, Real l_dt, Vector<MultiFab>& mom_mf, Mult
     // Compute divergence which will form RHS
     // Note that we replace "rho0w" with the contravariant momentum, Omega
     // ****************************************************************************
-    compute_divergence(lev, rhs[0], mom_mf, geom_tmp[0]);
+    Array<MultiFab const*, AMREX_SPACEDIM> rho0_u_const;
+    rho0_u_const[0] = &mom_mf[IntVars::xmom];
+    rho0_u_const[1] = &mom_mf[IntVars::ymom];
+    rho0_u_const[2] = &mom_mf[IntVars::zmom];
+
+    compute_divergence(lev, rhs[0], rho0_u_const, geom_tmp[0]);
 
     Real rhsnorm = rhs[0].norm0();
 
@@ -212,7 +217,7 @@ void ERF::project_velocities (int lev, Real l_dt, Vector<MultiFab>& mom_mf, Mult
     //
     if (mg_verbose > 0)
     {
-        compute_divergence(lev, rhs[0], mom_mf, geom_tmp[0]);
+        compute_divergence(lev, rhs[0], rho0_u_const, geom_tmp[0]);
 
         Print() << "Max/L2 norm of divergence  after solve at level " << lev << " : " << rhs[0].norm0() << " " << rhs[0].norm2() << std::endl;
 

Source/LinearSolvers/Make.package

@@ -1,10 +1,9 @@
 CEXE_headers += ERF_TerrainPoisson_3D_K.H
 CEXE_headers += ERF_TerrainPoisson.H
-CEXE_sources += ERF_TerrainPoisson.cpp
 
-CEXE_sources += ERF_ComputeDivergence.cpp
 CEXE_sources += ERF_PoissonSolve.cpp
 CEXE_sources += ERF_PoissonSolve_tb.cpp
+CEXE_sources += ERF_TerrainPoisson.cpp
 CEXE_sources += ERF_SolveWithGMRES.cpp
 CEXE_sources += ERF_SolveWithMLMG.cpp
 

Source/SourceTerms/ERF_MakeBuoyancy.cpp

@@ -40,9 +40,6 @@ void make_buoyancy (Vector<MultiFab>& S_data,
 {
     BL_PROFILE("make_buoyancy()");
 
-    const    Array<Real,AMREX_SPACEDIM> grav{0.0, 0.0, -solverChoice.gravity};
-    const GpuArray<Real,AMREX_SPACEDIM> grav_gpu{grav[0], grav[1], grav[2]};
-
     const int klo = geom.Domain().smallEnd()[2];
     const int khi = geom.Domain().bigEnd()[2] + 1;
 
@@ -53,6 +50,8 @@ void make_buoyancy (Vector<MultiFab>& S_data,
     MultiFab p0 (base_state, make_alias, BaseState::p0_comp , 1);
     MultiFab th0(base_state, make_alias, BaseState::th0_comp, 1);
 
+    Real grav_gpu = -solverChoice.gravity;
+
 #ifdef _OPENMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
@@ -83,7 +82,7 @@ void make_buoyancy (Vector<MultiFab>& S_data,
                 //
                 // Return -rho0 g (thetaprime / theta0)
                 //
-                buoyancy_fab(i, j, k) = buoyancy_dry_anelastic(i,j,k,grav_gpu[2],
+                buoyancy_fab(i, j, k) = buoyancy_dry_anelastic(i,j,k,grav_gpu,
                                                                r0_arr,th0_arr,cell_data);
             });
         }
@@ -97,7 +96,7 @@ void make_buoyancy (Vector<MultiFab>& S_data,
                 //
                 // Return -rho0 g (thetaprime / theta0)
                 //
-                buoyancy_fab(i, j, k) = buoyancy_moist_anelastic(i,j,k,grav_gpu[2],rv_over_rd,
+                buoyancy_fab(i, j, k) = buoyancy_moist_anelastic(i,j,k,grav_gpu,rv_over_rd,
                                                                  r0_arr,th0_arr,cell_data);
             });
         }
@@ -114,7 +113,7 @@ void make_buoyancy (Vector<MultiFab>& S_data,
                     //
                     // Return -rho0 g (thetaprime / theta0)
                     //
-                    buoyancy_fab(i, j, k) = buoyancy_rhopert(i,j,k,n_q_dry,grav_gpu[2],
+                    buoyancy_fab(i, j, k) = buoyancy_rhopert(i,j,k,n_q_dry,grav_gpu,
                                                              r0_arr,cell_data);
                 });
             }
@@ -125,7 +124,7 @@ void make_buoyancy (Vector<MultiFab>& S_data,
                     //
                     // Return -rho0 g (Tprime / T0)
                     //
-                    buoyancy_fab(i, j, k) = buoyancy_dry_Tpert(i,j,k,grav_gpu[2],rd_over_cp,
+                    buoyancy_fab(i, j, k) = buoyancy_dry_Tpert(i,j,k,grav_gpu,rd_over_cp,
                                                                r0_arr,p0_arr,th0_arr,cell_data);
                 });
             }
@@ -136,7 +135,7 @@ void make_buoyancy (Vector<MultiFab>& S_data,
                     //
                     // Return -rho0 g (Theta_prime / Theta_0)
                     //
-                    buoyancy_fab(i, j, k) = buoyancy_dry_Thpert(i,j,k,grav_gpu[2],
+                    buoyancy_fab(i, j, k) = buoyancy_dry_Thpert(i,j,k,grav_gpu,
                                                                 r0_arr,th0_arr,cell_data);
                 });
             } // buoyancy_type for dry compressible
@@ -158,7 +157,7 @@ void make_buoyancy (Vector<MultiFab>& S_data,
             {
                 ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
                 {
-                    buoyancy_fab(i, j, k) = buoyancy_rhopert(i,j,k,n_qstate,grav_gpu[2],
+                    buoyancy_fab(i, j, k) = buoyancy_rhopert(i,j,k,n_qstate,grav_gpu,
                                                              r0_arr,cell_data);
                 });
             }
@@ -167,7 +166,7 @@ void make_buoyancy (Vector<MultiFab>& S_data,
 
                 ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
                 {
-                    buoyancy_fab(i, j, k) = buoyancy_moist_Tpert(i,j,k,n_qstate,grav_gpu[2],rd_over_cp,
+                    buoyancy_fab(i, j, k) = buoyancy_moist_Tpert(i,j,k,n_qstate,grav_gpu,rd_over_cp,
                                                                  r0_arr,th0_arr,p0_arr,
                                                                  cell_prim,cell_data);
                 });
@@ -176,7 +175,7 @@ void make_buoyancy (Vector<MultiFab>& S_data,
             {
                 ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
                 {
-                    buoyancy_fab(i, j, k) = buoyancy_moist_Thpert(i,j,k,n_qstate,grav_gpu[2],
+                    buoyancy_fab(i, j, k) = buoyancy_moist_Thpert(i,j,k,n_qstate,grav_gpu,
                                                                   r0_arr,th0_arr,cell_prim);
                     });
             }

Source/TimeIntegration/ERF_FastRhs_MT.cpp

@@ -29,7 +29,8 @@ using namespace amrex;
  * @param[in   ] detJ_cc_old Jacobian of the metric transformation at old time
  * @param[in   ] detJ_cc_new Jacobian of the metric transformation at new time
  * @param[in   ] detJ_cc_stg Jacobian of the metric transformation at previous stage
- * @param[in   ] dtau fast time step
+ * @param[in   ] delta_r0
+ * @param[in   ] dtau        fast time step
  * @param[in   ] beta_s  Coefficient which determines how implicit vs explicit the solve is
  * @param[in   ] facinv inverse factor for time-averaging the momenta
  * @param[in   ] mapfac_m map factor at cell centers
@@ -63,6 +64,7 @@ void erf_fast_rhs_MT (int step, int nrk,
                       std::unique_ptr<MultiFab>& detJ_cc_old,        // at previous substep time (tau)
                       std::unique_ptr<MultiFab>& detJ_cc_new,        // at      new substep time (tau + delta tau)
                       std::unique_ptr<MultiFab>& detJ_cc_stg,        // at last RK stg
+                      const MultiFab& delta_r0,
                       const Real dtau, const Real beta_s,
                       const Real facinv,
                       std::unique_ptr<MultiFab>& mapfac_m,
@@ -174,6 +176,8 @@ void erf_fast_rhs_MT (int step, int nrk,
 
         const Array4<Real>& theta_extrap = extrap.array(mfi);
 
+        const Array4<Real const>& delta_r0_arr = delta_r0.const_array(mfi);
+
         // Map factors
         const Array4<const Real>& mf_m = mapfac_m->const_array(mfi);
         const Array4<const Real>& mf_u = mapfac_u->const_array(mfi);
@@ -444,7 +448,10 @@ void erf_fast_rhs_MT (int step, int nrk,
 
             // line 1
             RHS_a(i,j,k) = dJ_old_kface * prev_zmom(i,j,k) - dJ_stg_kface * stg_zmom(i,j,k)
-                            + dtau *(slow_rhs_rho_w(i,j,k) + R0_tmp + dtau*beta_2*R1_tmp );
+                            + dtau *(slow_rhs_rho_w(i,j,k) + R0_tmp + dtau*beta_2*R1_tmp);
+
+            // Subtract the RHS with the time-varying base state adjustment
+            RHS_a(i,j,k) -= dtau * halfg * (delta_r0_arr(i,j,k) + delta_r0_arr(i,j,k-1));
 
             // We cannot use omega_arr here since that was built with old_rho_u and old_rho_v ...
             Real UppVpp = dJ_new_kface * OmegaFromW(i,j,k,0.,cur_xmom,cur_ymom,z_nd_new,dxInv)

Source/TimeIntegration/ERF_MRI.H

@@ -22,8 +22,8 @@ private:
     std::function<void(T&, T&, T&, T&, const amrex::Real, const amrex::Real, const amrex::Real, const int)> slow_rhs_pre;
     std::function<void(T&, T&, T&, T&, const amrex::Real, const amrex::Real, const amrex::Real, const int)> slow_rhs_inc;
     std::function<void(T&, T&, T&, T&, T&, const amrex::Real, const amrex::Real, const amrex::Real, const int )> slow_rhs_post;
-    std::function<void(int, int, int, T&, const T&, T&, T&, T&, const amrex::Real, const amrex::Real,
-                                                                const amrex::Real, const amrex::Real)> fast_rhs;
+    std::function<void(int, int, int, T&, const T&, T&, T&, T&, const amrex::Real, const amrex::Real, const amrex::Real,
+                                                                const amrex::Real, const amrex::Real, const amrex::Real)> fast_rhs;
 
    /**
     * \brief Integrator timestep size (Real)
@@ -145,6 +145,7 @@ public:
     }
 
     void set_fast_rhs (std::function<void(int, int, int, T&, const T&, T&, T&, T&,
+                                          const amrex::Real, const amrex::Real,
                                           const amrex::Real, const amrex::Real,
                                           const amrex::Real, const amrex::Real)> F)
     {
@@ -289,8 +290,10 @@ public:
                 // *******************************************************************************
                 for (int ks = 0; ks < nsubsteps; ++ks)
                 {
+                    // We pass time and time_stage so that -- in the case of moving terrain --
+                    //    we can interpolate the base state between old_step_time and time_stage
                     fast_rhs(ks, nsubsteps, nrk, *F_slow, S_old, S_new, *S_sum, *S_scratch, dtau, inv_fac,
-                             time + ks*dtau, time + (ks+1) * dtau);
+                             time, time_stage, time + ks*dtau, time + (ks+1) * dtau);
 
                 } // ks
 

Source/TimeIntegration/ERF_TI_fast_headers.H

@@ -91,6 +91,7 @@ void erf_fast_rhs_MT (int step, int nrk, int level, int finest_level,
                       std::unique_ptr<amrex::MultiFab>& detJ_cc_old,
                       std::unique_ptr<amrex::MultiFab>& detJ_cc_new,
                       std::unique_ptr<amrex::MultiFab>& detJ_cc_stg,
+                      const amrex::MultiFab& delta_r0_at_fast_rhs_time,
                       const amrex::Real dtau, const amrex::Real beta_s,
                       const amrex::Real facinv,
                       std::unique_ptr<amrex::MultiFab>& mapfac_m,

Source/TimeIntegration/ERF_TI_slow_rhs_fun.H

@@ -157,6 +157,7 @@
                                   xflux_imask[level], yflux_imask[level], zflux_imask[level],
                                   thin_xforce[level], thin_yforce[level], thin_zforce[level]);
 
+
             // We define and evolve (rho theta)_0 in order to re-create p_0 in a way that is consistent
             //    with our update of (rho theta) but does NOT maintain dp_0 / dz = -rho_0 g.  This is why
             //    we no longer discretize the vertical pressure gradient in perturbational form.
@@ -175,6 +176,10 @@
 
             const Real l_rdOcp   = solverChoice.rdOcp;
 
+            // First copy base_state_new into base_state_src, so that base_state_src always holds the
+            //  base state at the time where we most recently computed the slow RHS
+            MultiFab::Copy(base_state_src[level], base_state_new[level], 0, 0, BaseState::num_comps,base_state_new[level].nGrowVect());
+
 #ifdef _OPENMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif