SHARPlib v1.4.0

CMakeLists.txt

@@ -31,6 +31,10 @@ file(GLOB_RECURSE SOURCES CONFIGURE_DEPENDS
 
 add_library(SHARPlib ${SOURCES})
 
+target_compile_definitions(SHARPlib PUBLIC 
+    $<BUILD_INTERFACE:FMT_HEADER_ONLY>
+)
+
 target_include_directories(SHARPlib PUBLIC 
     $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
     $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/external/fmt/include>

docs/reference/cpp_api/parcel.rst

@@ -9,6 +9,12 @@ Parcel Lifters
 .. doxygenstruct:: sharp::lifter_cm1
    :members:
 
+.. doxygenstruct:: sharp::lut_data
+   :members:
+
+.. doxygenstruct:: sharp::lifter_lut
+   :members:
+
 Parcels
 -------
 .. doxygenenum:: sharp::LPL

docs/reference/cpp_api/thermo.rst

@@ -32,8 +32,8 @@ Moist Variables
 .. doxygenfunction:: sharp::buoyancy(const float, const float)
 .. doxygenfunction:: sharp::buoyancy(const float[], const float[], float[], std::ptrdiff_t)
 .. doxygenfunction:: sharp::thetae
-.. doxygenfunction:: sharp::wetbulb(Lft, float, float, float)
-.. doxygenfunction:: sharp::theta_wetbulb(Lft, float, float, float)
+.. doxygenfunction:: sharp::wetbulb(Lft&, float, float, float)
+.. doxygenfunction:: sharp::theta_wetbulb(Lft&, float, float, float)
 
 Moist Ascent
 ------------

docs/reference/python_api/parcel.rst

@@ -9,17 +9,36 @@ SHARPlib supports configurable parcel ascent via parcel lifter 'functors'. These
 
 .. autoclass:: nwsspc.sharp.calc.parcel.lifter_wobus
 
+   .. autoproperty:: nwsspc.sharp.calc.parcel.lifter_wobus.ma_type
+   .. autoproperty:: nwsspc.sharp.calc.parcel.lifter_wobus.lift_from_lcl
    .. autofunction:: nwsspc.sharp.calc.parcel.lifter_wobus.__call__
    .. autofunction:: nwsspc.sharp.calc.parcel.lifter_wobus.parcel_virtual_temperature
 
 .. autoclass:: nwsspc.sharp.calc.parcel.lifter_cm1
 
    .. autoproperty:: nwsspc.sharp.calc.parcel.lifter_cm1.converge
    .. autoproperty:: nwsspc.sharp.calc.parcel.lifter_cm1.ma_type
+   .. autoproperty:: nwsspc.sharp.calc.parcel.lifter_cm1.lift_from_lcl
    .. autoproperty:: nwsspc.sharp.calc.parcel.lifter_cm1.pressure_incr
    .. autofunction:: nwsspc.sharp.calc.parcel.lifter_cm1.__call__
    .. autofunction:: nwsspc.sharp.calc.parcel.lifter_cm1.parcel_virtual_temperature
 
+.. autofunction:: nwsspc.sharp.calc.parcel.lut_data
+.. autoclass:: nwsspc.sharp.calc.parcel.lut_data_wobus
+.. autoclass:: nwsspc.sharp.calc.parcel.lut_data_cm1
+.. autofunction:: nwsspc.sharp.calc.parcel.lifter_lut
+.. autoclass:: nwsspc.sharp.calc.parcel.lifter_lut_wobus
+
+   .. autoproperty:: nwsspc.sharp.calc.parcel.lifter_lut_wobus.lift_from_lcl
+   .. autofunction:: nwsspc.sharp.calc.parcel.lifter_lut_wobus.__call__
+   .. autofunction:: nwsspc.sharp.calc.parcel.lifter_lut_wobus.parcel_virtual_temperature
+
+.. autoclass:: nwsspc.sharp.calc.parcel.lifter_lut_cm1
+
+   .. autoproperty:: nwsspc.sharp.calc.parcel.lifter_lut_cm1.lift_from_lcl
+   .. autofunction:: nwsspc.sharp.calc.parcel.lifter_lut_cm1.__call__
+   .. autofunction:: nwsspc.sharp.calc.parcel.lifter_lut_cm1.parcel_virtual_temperature
+
 Parcels
 -------
 Parcel creation, parcel ascent, and integrated CAPE/CINH.

examples/Python/hrrr-example.ipynb

@@ -54,6 +54,7 @@
     "\n",
     "client = Client(\n",
     "    n_workers=8, \n",
+    "    threads_per_worker=1,\n",
     "    memory_limit='8GB', \n",
     "    silence_logs=logging.ERROR\n",
     ")\n",
@@ -195,6 +196,7 @@
     "    pres, hght, \n",
     "    tmpk, spfh, \n",
     "    uwin, vwin, \n",
+    "    lut_data=None\n",
     "):\n",
     "    def qc(val):\n",
     "        if (val == constants.MISSING): return np.nan\n",
@@ -220,12 +222,17 @@
     "        theta,\n",
     "        mixr\n",
     "    )\n",
-    "\n",
-    "    lifter = parcel.lifter_cm1()\n",
-    "    lifter.ma_type = thermo.adiabat.pseudo_liq\n",
-    "    lifter.converge = 0.15 # K\n",
-    "    lifter.pressure_incr = 500.0 # 5 hPa\n",
-    "    \n",
+    "    ## If a lookup table is not provided,\n",
+    "    ## construct a parcel lifter for use\n",
+    "    lifter = None\n",
+    "    if lut_data is None:\n",
+    "        lifter = parcel.lifter_cm1()\n",
+    "        lifter.ma_type = thermo.adiabat.pseudo_liq\n",
+    "        lifter.converge = 0.15 # K\n",
+    "        lifter.pressure_incr = 500.0 # 5 hPa\n",
+    "    else:\n",
+    "        lifter = parcel.lifter_lut(lut_data)\n",
+    "        \n",
     "    # lift the parcel and get CAPE/CINH\n",
     "    pcl_vtmp = mlpcl.lift_parcel(lifter, pres)\n",
     "    pcl_buoy = thermo.buoyancy(pcl_vtmp, vtmp)\n",
@@ -273,28 +280,16 @@
     "        vwin\n",
     "    )\n",
     "\n",
-    "    if mupcl.eql_pressure != constants.MISSING:\n",
-    "        # Get the effective bulk wind difference\n",
-    "        eql_hght = interp.interp_pressure(\n",
-    "            mupcl.eql_pressure,\n",
-    "            pres,\n",
-    "            hght\n",
-    "        )\n",
-    "        depth = (eql_hght - eil_hght.bottom)*0.5\n",
-    "        ebwd_lyr = layer.HeightLayer(\n",
-    "            eil_hght.bottom, \n",
-    "            float(eil_hght.bottom + depth)\n",
-    "        )\n",
-    "        ebwd_cmp = winds.wind_shear(\n",
-    "            ebwd_lyr,\n",
-    "            hght,\n",
-    "            uwin,\n",
-    "            vwin\n",
-    "        )\n",
-    "        ebwd = winds.vector_magnitude(ebwd_cmp.u, ebwd_cmp.v)\n",
-    "    else:\n",
-    "        ebwd_cmp = winds.WindComponents(constants.MISSING, constants.MISSING)\n",
-    "        ebwd = constants.MISSING\n",
+    "    ebwd_cmp = params.effective_bulk_wind_difference(\n",
+    "        pres, \n",
+    "        hght,\n",
+    "        uwin,\n",
+    "        vwin,\n",
+    "        eil,\n",
+    "        mupcl.eql_pressure\n",
+    "    )\n",
+    "\n",
+    "    ebwd = winds.vector_magnitude(ebwd_cmp.u, ebwd_cmp.v)\n",
     "    \n",
     "    stp = params.significant_tornado_parameter(\n",
     "        mlpcl,\n",
@@ -396,6 +391,90 @@
     "print(result_ds)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "cc3224ba-5fc3-42f2-a94e-5e9aa52306e7",
+   "metadata": {},
+   "source": [
+    "## **NEW in SHARPlib v1.4.0 - Lookup Table (LUT) Based Parcel Ascent!**\n",
+    "It is a fairly common performance optimization (especially when processing large gridded datasets) to precompute a lookup table for parcel ascent. Instead of solving the ODEs directly for every grid point and every parcel lift operation, you create a sufficiently dense lookup table that you can perform interpolation operations on. This can dramatically improve performance, and as of version 1.4.0, SHARPlib supports generating lookup table data with any of the configurable parcel lifters. The catch: this only works for ***pseudoadiabatic*** ascent, and does not work for reversible (adiabatic) ascent. These next cells demonstrate the construction of the LUT data and the usage of ```lifter_lut```.\n",
+    "\n",
+    "The LUT data stored in the ```lut_data``` struct can be shared across multiple threads without making copies (it will make copies if using multiple nodes) and is thread safe. Each thread should have its own ```lifter_lut```, however."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e2ec4b2b-6935-4b18-be29-2f4a59d0d842",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%time\n",
+    "lifter_cm1 = parcel.lifter_cm1()\n",
+    "lifter_cm1.ma_type = thermo.adiabat.pseudo_liq\n",
+    "lifter_cm1.converge = 0.001 # K\n",
+    "lifter_cm1.pressure_incr = 50.0 # every 0.5 hPa\n",
+    "\n",
+    "## This will construct the LUT with the defaults. The LUT \n",
+    "## bounds and discretization can be specified. These \n",
+    "## defaults are in the documentation, but provided for \n",
+    "## clarity below:\n",
+    "## min_pressure: 5000.0 Pa\n",
+    "## max_pressure: 110000.0 Pa\n",
+    "## min_thetae: 210.0 K\n",
+    "## max_thetae: 430.0f\n",
+    "## num_logp: 201\n",
+    "## num_thetae: 221\n",
+    "lut_cm1 = parcel.lut_data(lifter_cm1)\n",
+    "print(help(parcel.lut_data))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "69e01d81-babe-40d1-b79d-dcba6db54b55",
+   "metadata": {},
+   "source": [
+    "Now we can re-run our computations, this time passing the lookup table data to the computation function and compare the execution time! A note on the comparison, however: when no LUT is provided, the function configures the ```lifter_cm1``` function with an integration step of `5 hPa` and a convergence accuracy of `0.15 K`. The lookup table, however, was constructed with an integration step of `0.5 hPa` and a convergence accuracy of `0.001 K`. So, the lookup table should not only be faster than directly calling `lifter_cm1`, but also more accurate!"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9b41dd91-93a9-4118-b7b6-b1505c91da1a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%time\n",
+    "result_func = xr.apply_ufunc(\n",
+    "    compute_everything, \n",
+    "    ds_merged[\"pres\"], ds_merged[\"gh\"],\n",
+    "    ds_merged[\"t\"], ds_merged[\"q\"],\n",
+    "    ds_merged[\"u\"], ds_merged[\"v\"],\n",
+    "    input_core_dims=input_core_dims,\n",
+    "    output_core_dims=output_core_dims,\n",
+    "    vectorize=True,\n",
+    "    dask=\"parallelized\",\n",
+    "    output_dtypes=output_dtypes,\n",
+    "    ## pass our LUT as a keyword \n",
+    "    ## argument to the compute \n",
+    "    ## function\n",
+    "    kwargs={\"lut_data\": lut_cm1},\n",
+    ")\n",
+    "result_dict = {name: data for name, data in zip(output_names, result_func)}\n",
+    "result_ds = xr.Dataset(result_dict)\n",
+    "result_ds = result_ds.compute()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fe39f48f-8909-47ad-a16c-adc81ac661c4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(result_ds)"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "6252d8f4-c279-444b-a310-75c57f415e16",

include/SHARPlib/layer.h

@@ -479,20 +479,20 @@ template <typename L>
                                           const std::ptrdiff_t N,
                                           const int integ_sign = 0,
                                           const bool weighted = false) {
-    float var_lyr_bottom;
-    float coord_lyr_bottom = layer.bottom;
-
-    float var_lyr_top;
-    float coord_lyr_top = layer.top;
-
     float integrated = 0.0;
     float weights = 0.0;
 
     const bool isign = std::signbit(integ_sign);
 
-    // using constexpr means that this if statement optimizes
-    // away at compile time since the layer coordinate is known
     LayerIndex idx = get_layer_index(layer, coord_array, N);
+
+    float var_lyr_bottom;
+    float var_lyr_top;
+
+    auto include = [integ_sign, isign](float layer_avg) -> float {
+        return (float)((!integ_sign) | (isign == std::signbit(layer_avg)));
+    };
+
     if constexpr (layer.coord == LayerCoordinate::height) {
         var_lyr_bottom = interp_height(layer.bottom, coord_array, var_array, N);
         var_lyr_top = interp_height(layer.top, coord_array, var_array, N);
@@ -502,64 +502,72 @@ template <typename L>
         var_lyr_top = interp_pressure(layer.top, coord_array, var_array, N);
     }
 
+    // When the layer falls entirely between two native levels
+    // with no interior points, handle as a single trapezoid.
+    // NOTE: weights is accumulated by _integ_trapz through its
+    // reference parameter whenever weighted=true. The same weights
+    // variable must be passed to every _integ_trapz call in this
+    // function for the final division to be correct.
     if (idx.ktop < idx.kbot) {
         if ((var_lyr_bottom == MISSING) || (var_lyr_top == MISSING)) {
             return MISSING;
         }
         float layer_avg = _integ_trapz(var_lyr_top, var_lyr_bottom, layer.top,
                                        layer.bottom, weights, weighted);
 
-        float cond =
-            (float)((!integ_sign) | (isign == std::signbit(layer_avg)));
-        integrated += cond * layer_avg;
+        integrated += include(layer_avg) * layer_avg;
     } else {
-        for (std::ptrdiff_t k = idx.kbot; k < idx.ktop; ++k) {
+        // Start from the interpolated layer bottom and walk
+        // through native levels to the interpolated top.
+        // When a level is MISSING, the next valid level forms a
+        // wider trapezoid from the last valid point, bridging
+        // the gap via linear interpolation.
+        float coord_prev = layer.bottom;
+        float var_prev = var_lyr_bottom;
+        bool have_prev = (var_lyr_bottom != MISSING);
+        bool any_segment = false;
+
+        for (std::ptrdiff_t k = idx.kbot; k <= idx.ktop; ++k) {
 #ifndef NO_QC
             if (var_array[k] == MISSING) {
                 continue;
             }
 #endif
+            if (have_prev) {
+                float layer_avg =
+                    _integ_trapz(var_array[k], var_prev, coord_array[k],
+                                 coord_prev, weights, weighted);
 
-            float coord_bottom = coord_array[k];
-            float var_bottom = var_array[k];
-
-            float coord_top = coord_array[k + 1];
-            float var_top = var_array[k + 1];
-
-            float layer_avg = _integ_trapz(var_top, var_bottom, coord_top,
-                                           coord_bottom, weights, weighted);
+                integrated += include(layer_avg) * layer_avg;
+                any_segment = true;
+            }
 
-            float cond =
-                (float)((!integ_sign) | (isign == std::signbit(layer_avg)));
-            integrated += cond * layer_avg;
+            coord_prev = coord_array[k];
+            var_prev = var_array[k];
+            have_prev = true;
         }
 
-        // interpolated bottom of layer
+        // Final segment: last valid native level to the
+        // interpolated layer top
+        if (have_prev && (var_lyr_top != MISSING)) {
+            float layer_avg = _integ_trapz(var_lyr_top, var_prev, layer.top,
+                                           coord_prev, weights, weighted);
 
-        if (var_lyr_bottom != sharp::MISSING) {
-            float layer_avg = _integ_trapz(var_array[idx.kbot], var_lyr_bottom,
-                                           coord_array[idx.kbot],
-                                           coord_lyr_bottom, weights, weighted);
-            float cond = ((!integ_sign) | (isign == std::signbit(layer_avg)));
-            integrated += cond * layer_avg;
-        }
-
-        // interpolated top of layer
-        if (var_lyr_top != sharp::MISSING) {
-            float layer_avg =
-                _integ_trapz(var_lyr_top, var_array[idx.ktop], coord_lyr_top,
-                             coord_array[idx.ktop], weights, weighted);
-            float cond = ((!integ_sign) | (isign == std::signbit(layer_avg)));
-            integrated += cond * layer_avg;
+            integrated += include(layer_avg) * layer_avg;
+            any_segment = true;
         }
+        if (!any_segment) return MISSING;
     }
 
     if constexpr (layer.coord == LayerCoordinate::pressure) {
         integrated *= -1.0f;
         weights *= -1.0f;
     }
 
-    if (weighted) integrated /= weights;
+    if (weighted) {
+        if (weights == 0.0f) return MISSING;
+        integrated /= weights;
+    }
     return integrated;
 }