Fix calls to getTgivenRandRTh in moist path (#1677)

* In make_buoyancy, include qv in EOS calls under the moist path * Call diff dertemp function if moist * Update gold data -- moist temp output changed --------- Co-authored-by: Aaron M. Lattanzi <[email protected]>
erf-model · Jul 10, 2024 · 7878482 · 7878482
1 parent 59da3d2
commit 7878482
Show file tree

Hide file tree

Showing 7 changed files with 32,251 additions and 26,548 deletions.
diff --git a/Source/Derive.H b/Source/Derive.H
@@ -45,6 +45,17 @@ void erf_dertemp (
   const int* bcrec,
   const int level);
 
+void erf_dermoisttemp (
+  const amrex::Box& bx,
+  amrex::FArrayBox& derfab,
+  int dcomp,
+  int ncomp,
+  const amrex::FArrayBox& datfab,
+  const amrex::Geometry& geomdata,
+  amrex::Real time,
+  const int* bcrec,
+  const int level);
+
 void erf_dertheta (
   const amrex::Box& bx,
   amrex::FArrayBox& derfab,

diff --git a/Source/Derive.cpp b/Source/Derive.cpp
@@ -109,6 +109,29 @@ erf_dertemp (const Box& bx,
         tfab(i,j,k) = getTgivenRandRTh(rho,rhotheta);
     });
 }
+void
+erf_dermoisttemp (const Box& bx,
+             FArrayBox& derfab,
+             int /*dcomp*/,
+             int /*ncomp*/,
+             const FArrayBox& datfab,
+             const Geometry& /*geomdata*/,
+             Real /*time*/,
+             const int* /*bcrec*/,
+             const int /*level*/)
+{
+    auto const dat = datfab.array();
+    auto tfab      = derfab.array();
+
+    ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+    {
+        const Real rho = dat(i, j, k, Rho_comp);
+        const Real rhotheta = dat(i, j, k, RhoTheta_comp);
+        AMREX_ALWAYS_ASSERT(rhotheta > 0.);
+        const Real qv = dat(i, j, k, RhoQ1_comp) / rho;
+        tfab(i,j,k) = getTgivenRandRTh(rho,rhotheta,qv);
+    });
+}
 
 /**
  * Function to define the potential temperature by calling an EOS routine

diff --git a/Source/IO/ERF_WriteBndryPlanes.cpp b/Source/IO/ERF_WriteBndryPlanes.cpp
@@ -143,6 +143,9 @@ void WriteBndryPlanes::write_planes (const int t_step, const Real time,
     Box target_box_shifted(IntVect(0,0,0),new_hi);
     BoxArray ba_shifted(target_box_shifted);
 
+    int n_moist_var = NMOIST_max - (S.nComp() - NVAR_max);
+    bool ismoist = (n_moist_var >= 1);
+
     for (int i = 0; i < m_var_names.size(); i++)
     {
         std::string var_name = m_var_names[i];
@@ -169,7 +172,11 @@ void WriteBndryPlanes::write_planes (const int t_step, const Real time,
             for (MFIter mfi(Temp, TilingIfNotGPU()); mfi.isValid(); ++mfi)
             {
                 const Box& bx = mfi.tilebox();
-                derived::erf_dertemp(bx, Temp[mfi], 0, 1, S[mfi], m_geom[bndry_lev], time, nullptr, bndry_lev);
+                if (ismoist) {
+                    derived::erf_dermoisttemp(bx, Temp[mfi], 0, 1, S[mfi], m_geom[bndry_lev], time, nullptr, bndry_lev);
+                } else {
+                    derived::erf_dertemp(bx, Temp[mfi], 0, 1, S[mfi], m_geom[bndry_lev], time, nullptr, bndry_lev);
+                }
             }
             bndry.copyFrom(Temp, nghost, 0, 0, ncomp, m_geom[bndry_lev].periodicity());
         } else if (var_name == "scalar") {

diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp
@@ -300,9 +300,15 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
             }
         };
 
+        bool ismoist = (solverChoice.moisture_type != MoistureType::None);
+
         // Note: All derived variables must be computed in order of "derived_names" defined in ERF.H
         calculate_derived("soundspeed",  vars_new[lev][Vars::cons], derived::erf_dersoundspeed);
-        calculate_derived("temp",        vars_new[lev][Vars::cons], derived::erf_dertemp);
+        if (ismoist) {
+            calculate_derived("temp",        vars_new[lev][Vars::cons], derived::erf_dermoisttemp);
+        } else {
+            calculate_derived("temp",        vars_new[lev][Vars::cons], derived::erf_dertemp);
+        }
         calculate_derived("theta",       vars_new[lev][Vars::cons], derived::erf_dertheta);
         calculate_derived("KE",          vars_new[lev][Vars::cons], derived::erf_derKE);
         calculate_derived("QKE",         vars_new[lev][Vars::cons], derived::erf_derQKE);

diff --git a/Source/SourceTerms/ERF_make_buoyancy.cpp b/Source/SourceTerms/ERF_make_buoyancy.cpp
@@ -232,11 +232,15 @@ void make_buoyancy (Vector<MultiFab>& S_data,
                     // TODO: ice has not been dealt with (q1=qv, q2=qv, q3=qp)
                     ParallelFor(tbz, [=, buoyancy_type=solverChoice.buoyancy_type] AMREX_GPU_DEVICE (int i, int j, int k)
                     {
-                        Real tempp1d = getTgivenRandRTh(rho_d_ptr[k  ], rho_d_ptr[k  ]*theta_d_ptr[k  ]);
-                        Real tempm1d = getTgivenRandRTh(rho_d_ptr[k-1], rho_d_ptr[k-1]*theta_d_ptr[k-1]);
+                        Real tempp1d = getTgivenRandRTh(rho_d_ptr[k  ], rho_d_ptr[k  ]*theta_d_ptr[k  ], qv_d_ptr[k  ]);
+                        Real tempm1d = getTgivenRandRTh(rho_d_ptr[k-1], rho_d_ptr[k-1]*theta_d_ptr[k-1], qv_d_ptr[k-1]);
 
-                        Real tempp3d  = getTgivenRandRTh(cell_data(i,j,k  ,Rho_comp), cell_data(i,j,k  ,RhoTheta_comp));
-                        Real tempm3d  = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp), cell_data(i,j,k-1,RhoTheta_comp));
+                        Real tempp3d  = getTgivenRandRTh(cell_data(i,j,k  ,Rho_comp),
+                                                         cell_data(i,j,k  ,RhoTheta_comp),
+                                                         cell_data(i,j,k  ,RhoQ1_comp)/cell_data(i,j,k  ,Rho_comp));
+                        Real tempm3d  = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp),
+                                                         cell_data(i,j,k-1,RhoTheta_comp),
+                                                         cell_data(i,j,k-1,RhoQ1_comp)/cell_data(i,j,k-1,Rho_comp));
 
                         Real qplus, qminus;
 
@@ -300,11 +304,15 @@ void make_buoyancy (Vector<MultiFab>& S_data,
 
                     ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
                     {
-                        Real tempp1d = getTgivenRandRTh(rho_d_ptr[k  ], rho_d_ptr[k  ]*theta_d_ptr[k  ]);
-                        Real tempm1d = getTgivenRandRTh(rho_d_ptr[k-1], rho_d_ptr[k-1]*theta_d_ptr[k-1]);
-
-                        Real tempp3d  = getTgivenRandRTh(cell_data(i,j,k  ,Rho_comp), cell_data(i,j,k  ,RhoTheta_comp));
-                        Real tempm3d  = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp), cell_data(i,j,k-1,RhoTheta_comp));
+                        Real tempp1d = getTgivenRandRTh(rho_d_ptr[k  ], rho_d_ptr[k  ]*theta_d_ptr[k  ], qv_d_ptr[k  ]);
+                        Real tempm1d = getTgivenRandRTh(rho_d_ptr[k-1], rho_d_ptr[k-1]*theta_d_ptr[k-1], qv_d_ptr[k-1]);
+
+                        Real tempp3d  = getTgivenRandRTh(cell_data(i,j,k  ,Rho_comp),
+                                                         cell_data(i,j,k  ,RhoTheta_comp),
+                                                         cell_data(i,j,k  ,RhoQ1_comp)/cell_data(i,j,k  ,Rho_comp));
+                        Real tempm3d  = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp),
+                                                         cell_data(i,j,k-1,RhoTheta_comp),
+                                                         cell_data(i,j,k-1,RhoQ1_comp)/cell_data(i,j,k-1,Rho_comp));
 
                         Real qv_plus  = (n_moist_var >= 1) ? cell_prim(i,j,k  ,PrimQ1_comp) : 0.0;
                         Real qv_minus = (n_moist_var >= 1) ? cell_prim(i,j,k-1,PrimQ1_comp) : 0.0;