Initialize non-zero base state for wrfinput made from ideal.exe (#3730)

Author: asalmgren

Source/Initialization/ERF_InitFromWRFInput.cpp

@@ -163,58 +163,58 @@ read_base_state_params_from_wrfinput (const std::string& fname,
     if (ParallelDescriptor::IOProcessor()) {
         auto ncf = ncutils::NCFile::open(fname, NC_CLOBBER | NC_NETCDF4);
 
+        // Remember what was passed in so we can fall back to it below
+        const Real T00_def = T00;
+        const Real P00_def = P00;
+        const Real TLP_def = TLP;
+        const Real TISO_def = TISO;
+        const Real TLP_STRAT_def = TLP_STRAT;
+        const Real P_STRAT_def   = P_STRAT;
+
         std::vector<size_t> shape;
         std::vector<size_t> start;
-        int success = ncf.has_var("T00");
-        if (success) {
-            Print() << "Reading T00 from wrfinput\n";
-            shape = ncf.var("T00").shape();
-            start.resize(shape.size(), 0);
-            ncf.var("T00").get(&T00 ,start, shape);
-        }
-
-        success = ncf.has_var("P00");
-        if (success) {
-            Print() << "Reading P00 from wrfinput\n";
-            shape = ncf.var("P00").shape();
-            start.resize(shape.size(), 0);
-            ncf.var("P00").get(&P00 ,start, shape);
-        }
-
-        success = ncf.has_var("TLP");
-        if (success) {
-            Print() << "Reading TLP from wrfinput\n";
-            shape = ncf.var("TLP").shape();
+        auto read_scalar = [&] (const std::string& name, Real& val)
+        {
+            if (!ncf.has_var(name)) return;
+            Print() << "Reading " << name << " from wrfinput\n";
+            shape = ncf.var(name).shape();
+            start.clear();
             start.resize(shape.size(), 0);
-            ncf.var("TLP").get(&TLP ,start, shape);
-        }
+            ncf.var(name).get(&val, start, shape);
+        };
 
-        success = ncf.has_var("TISO");
-        if (success) {
-            Print() << "Reading TISO from wrfinput\n";
-            shape = ncf.var("TISO").shape();
-            start.resize(shape.size(), 0);
-            ncf.var("TISO").get(&TISO ,start, shape);
-        }
+        read_scalar("T00"      , T00);
+        read_scalar("P00"      , P00);
+        read_scalar("TLP"      , TLP);
+        read_scalar("TISO"     , TISO);
+        read_scalar("TLP_STRAT", TLP_STRAT);
+        read_scalar("P_STRAT"  , P_STRAT);
 
-        success = ncf.has_var("TLP_STRAT");
-        if (success) {
-            Print() << "Reading TLP_STRAT from wrfinput\n";
-            shape = ncf.var("TLP_STRAT").shape();
-            start.resize(shape.size(), 0);
-            ncf.var("TLP_STRAT").get(&TLP_STRAT ,start, shape);
-        }
+        ncf.close();
 
-        success = ncf.has_var("P_STRAT");
-        if (success) {
-            Print() << "Reading P_STRAT from wrfinput\n";
-            shape = ncf.var("P_STRAT").shape();
-            start.resize(shape.size(), 0);
-            ncf.var("P_STRAT").get(&P_STRAT ,start, shape);
+        // Idealized WRF cases (and some hand-built wrfinput files) declare these
+        // variables but leave them zero-filled.  T00, P00 and TISO are absolute
+        // temperatures/pressures, so a non-positive value is never meaningful; if
+        // any of them is bad we discard the whole group and keep the defaults.
+        // (TLP, TLP_STRAT and P_STRAT may legitimately be zero, so they are only
+        // reset alongside a bad T00/P00/TISO.)
+        const bool params_ok = std::isfinite(T00)  && (T00  > Real(0)) &&
+                               std::isfinite(P00)  && (P00  > Real(0)) &&
+                               std::isfinite(TISO) && (TISO > Real(0)) &&
+                               std::isfinite(TLP)  && std::isfinite(TLP_STRAT) &&
+                               std::isfinite(P_STRAT);
+
+        if (!params_ok) {
+            Print() << "WARNING: WRF base state parameters read from " << fname
+                    << " are invalid: (T00, P00, TLP, TISO, TLP_STRAT, P_STRAT) = ("
+                    << T00 << ", " << P00 << ", " << TLP << ", " << TISO << ", "
+                    << TLP_STRAT << ", " << P_STRAT << ")\n";
+            Print() << "         T00, P00 and TISO must all be positive and finite; "
+                       "reverting to ERF defaults.\n";
+            T00 = T00_def;   P00 = P00_def;   TLP = TLP_def;
+            TISO = TISO_def; TLP_STRAT = TLP_STRAT_def; P_STRAT = P_STRAT_def;
         }
 
-        ncf.close();
-
         Print() << "WRF base state parameters (T00, P00, TLP, TISO, TLP_STRAT, P_STRAT) are: ("
                 << T00 << ", " << P00 << ", " << TLP << ", " << TISO << ", " << TLP_STRAT << ", "
                 << P_STRAT << ") \n";
@@ -1629,6 +1629,12 @@ init_base_state_from_wrfinput (const Box& subdomain,
     int k_dom_lo = dom_lo.z;
     int k_dom_hi = dom_hi.z;
 
+    // The vertical integration below is seeded with the surface values (P00,T00),
+    // so these must be valid regardless of whether ALB was present in the file
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(std::isfinite(P00) && (P00 > Real(0)) &&
+                                     std::isfinite(T00) && (T00 > Real(0)),
+                                     "Cannot rebalance the WRF base state: P00 and T00 must be positive");
+
 #ifdef AMREX_USE_FLOAT
     Real tol  = Real(1.0e-6);
 #else

Source/Utils/ERF_HSEUtils.H

@@ -55,9 +55,22 @@ namespace HSEutils
                         Real& F,
                         const bool& maintain_Th)
     {
+        // A non-finite residual (or a non-positive initial guess) would make the
+        // convergence test below evaluate to false on the very first pass, so we
+        // would silently return the initial guess and leave, e.g., an all-zero
+        // base state.  Trap that here rather than propagating it downstream.
+        if (!(P > amrex::Real(0)) || !(rd > amrex::Real(0)) || !std::isfinite(F)) {
+            AMREX_DEVICE_PRINTF("ERROR: HSE Newton started from an invalid state: P = %e, rd = %e, F = %e\n",
+                                double(P), double(rd), double(F));
+            AMREX_DEVICE_PRINTF("       Check the surface pressure/temperature used to seed the integration%s\n", "");
+        }
+        AMREX_ALWAYS_ASSERT(P > amrex::Real(0) && rd > amrex::Real(0) && std::isfinite(F));
+
         int iter=0;
         int max_iter=20;
-        while (std::abs(F)>m_tol && iter<max_iter) {
+        // Written as !(|F| <= tol) so that a NaN residual enters the loop instead
+        // of being mistaken for convergence
+        while (!(std::abs(F)<=m_tol) && iter<max_iter) {
             // Compute change in pressure
             Real dRdP = (maintain_Th) ? iGamma * rd / P : rd / P;
             Real dFdp = one + myhalf*dRdP*g*dz;