add heating source term to erf (erf-model#1395)

wiersema1 · Jan 24, 2024 · 6a0ba3b · 6a0ba3b
1 parent 9ea9813
commit 6a0ba3b
Show file tree

Hide file tree

Showing 5 changed files with 38 additions and 2 deletions.
diff --git a/Source/ERF.H b/Source/ERF.H
@@ -529,6 +529,7 @@ private:
 
 #if defined(ERF_USE_RRTMGP)
     Radiation rad;
+    amrex::Vector<amrex::MultiFab> qheating_rates;  // radiation heating rate source terms
 #endif
 
     // Fillpatcher classes for coarse-fine boundaries

diff --git a/Source/ERF.cpp b/Source/ERF.cpp
@@ -118,6 +118,9 @@ ERF::ERF ()
     // NOTE: size micro before readparams (chooses the model at all levels)
     micro.ReSize(nlevs_max);
     qmoist.resize(nlevs_max);
+#if defined(ERF_USE_RRTMGP)
+    qheating_rates.resize(nlevs_max);
+#endif
 
     ReadParameters();
     const std::string& pv1 = "plot_vars_1"; setPlotVariables(pv1,plot_var_names_1);
@@ -1488,6 +1491,9 @@ ERF::ERF (const amrex::RealBox& rb, int max_level_in,
     // NOTE: size micro before readparams (chooses the model at all levels)
     micro.ReSize(nlevs_max);
     qmoist.resize(nlevs_max);
+#if defined(ERF_USE_RRTMGP)
+    qheating_rates.resize(nlevs_max);
+#endif
 
     ReadParameters();
     const std::string& pv1 = "plot_vars_1"; setPlotVariables(pv1,plot_var_names_1);

diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp
@@ -84,6 +84,9 @@ void erf_slow_rhs_pre (int level, int finest_level,
                        MultiFab& Omega,
                        const MultiFab& source,
                        const MultiFab& buoyancy,
+#ifdef ERF_USE_RRTMGP
+                       const MultiFab& qheating_rates,
+#endif
                        MultiFab* Tau11, MultiFab* Tau22, MultiFab* Tau33,
                        MultiFab* Tau12, MultiFab* Tau13, MultiFab* Tau21,
                        MultiFab* Tau23, MultiFab* Tau31, MultiFab* Tau32,
@@ -699,6 +702,21 @@ void erf_slow_rhs_pre (int level, int finest_level,
                     cell_rhs(i,j,k,RhoTheta_comp) += src_arr(i,j,k,RhoTheta_comp);
                 });
             }
+
+#ifdef ERF_USE_RRTMGP
+            auto const& qheating_arr = qheating_rates.const_array(mfi);
+            if (l_use_terrain && l_moving_terrain) {
+                amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                {
+                    cell_rhs(i,j,k,RhoTheta_comp) += (qheating_arr(i,j,k,0) + qheating_arr(i,j,k,1)) / detJ_arr(i,j,k);
+                });
+            } else {
+                amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                {
+                    cell_rhs(i,j,k,RhoTheta_comp) += qheating_arr(i,j,k,0) + qheating_arr(i,j,k,1);
+                });
+            }
+#endif
         }
 
         // Add Rayleigh damping

diff --git a/Source/TimeIntegration/TI_headers.H b/Source/TimeIntegration/TI_headers.H
@@ -26,6 +26,9 @@ void erf_slow_rhs_pre (int level, int finest_level, int nrk,
                             amrex::MultiFab& Omega,
                       const amrex::MultiFab& source,
                       const amrex::MultiFab& buoyancy,
+#if defined(ERF_USE_RRTMGP)
+                      const amrex::MultiFab& qheating_rates,
+#endif
                             amrex::MultiFab* Tau11,
                             amrex::MultiFab* Tau22,
                             amrex::MultiFab* Tau33,

diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H
@@ -94,7 +94,11 @@
 
             erf_slow_rhs_pre(level, finest_level, nrk, slow_dt, S_rhs, S_data, S_prim, S_scratch,
                              xvel_new, yvel_new, zvel_new,
-                             z_t_rk[level], Omega, source, buoyancy, Tau11, Tau22, Tau33, Tau12,
+                             z_t_rk[level], Omega, source, buoyancy,
+#if defined(ERF_USE_RRTMGP)
+                             qheating_rates[level],
+#endif
+                             Tau11, Tau22, Tau33, Tau12,
                              Tau13, Tau21,  Tau23, Tau31, Tau32, SmnSmn, eddyDiffs,
                              Hfx3, Diss,
                              fine_geom, solverChoice, m_most, domain_bcs_type_d, domain_bcs_type,
@@ -185,7 +189,11 @@
 
             erf_slow_rhs_pre(level, finest_level, nrk, slow_dt, S_rhs, S_data, S_prim, S_scratch,
                              xvel_new, yvel_new, zvel_new,
-                             z_t_rk[level], Omega, source, buoyancy, Tau11, Tau22, Tau33, Tau12,
+                             z_t_rk[level], Omega, source, buoyancy,
+#if defined(ERF_USE_RRTMGP)
+                             qheating_rates[level],
+#endif
+                             Tau11, Tau22, Tau33, Tau12,
                              Tau13, Tau21,  Tau23, Tau31, Tau32, SmnSmn, eddyDiffs,
                              Hfx3, Diss,
                              fine_geom, solverChoice, m_most, domain_bcs_type_d, domain_bcs_type,