From 5fd80015e9bca523ac52f22ca5b4fb662f76c2fb Mon Sep 17 00:00:00 2001 From: Ryan May Date: Fri, 6 Sep 2024 17:18:29 -0600 Subject: [PATCH 01/41] Update to build using scikit-build-core --- pyproject.toml | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 33bf3be94c..04ef962d24 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,9 +1,9 @@ [build-system] -requires = ["setuptools>=61", "wheel", "setuptools_scm[toml]>=3.4"] -build-backend = "setuptools.build_meta" +requires = ["setuptools>=61", "setuptools-scm[toml]>=3.4", "scikit-build-core", "pybind11"] +build-backend = "scikit_build_core.build" [project] -name = "MetPy" +name = "metpy" description = "Collection of tools for reading, visualizing and performing calculations with weather data." readme = "README.md" dynamic = ["version"] @@ -174,5 +174,9 @@ max-complexity = 61 [tool.ruff.lint.pydocstyle] convention = "numpy" +[tool.scikit-build] +metadata.version.provider = "scikit_build_core.metadata.setuptools_scm" +cmake.source-dir = "src" + [tool.setuptools_scm] version_scheme = "post-release" From 19399ff0d99b627c22aae3cbfc02a95a804150a9 Mon Sep 17 00:00:00 2001 From: Ryan May Date: Fri, 6 Sep 2024 17:18:57 -0600 Subject: [PATCH 02/41] WIP: Add stub module that uses pybind11 --- src/CMakeLists.txt | 15 +++++++++++++++ src/calcmod.cpp | 10 ++++++++++ 2 files changed, 25 insertions(+) create mode 100644 src/CMakeLists.txt create mode 100644 src/calcmod.cpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt new file mode 100644 index 0000000000..1080051ca0 --- /dev/null +++ b/src/CMakeLists.txt @@ -0,0 +1,15 @@ +cmake_minimum_required(VERSION 3.15...3.26) +project( + "${SKBUILD_PROJECT_NAME}" + LANGUAGES CXX + VERSION "${SKBUILD_PROJECT_VERSION}") + +find_package(Python COMPONENTS Interpreter Development.Module REQUIRED) +find_package(pybind11 CONFIG REQUIRED) + +pybind11_add_module(_calc_mod MODULE calcmod.cpp) + +target_compile_definitions(_calc_mod + PRIVATE VERSION_INFO=${PROJECT_VERSION}) + +install(TARGETS _calc_mod DESTINATION metpy) \ No newline at end of file diff --git a/src/calcmod.cpp b/src/calcmod.cpp new file mode 100644 index 0000000000..0d7837d01e --- /dev/null +++ b/src/calcmod.cpp @@ -0,0 +1,10 @@ +#include + +int add(int i, int j) { + return i + j; +} + +PYBIND11_MODULE(_calc_mod, m) { + m.doc() = "accelerator module docstring"; + m.def("add", &add, "Add two numbers"); +} \ No newline at end of file From e76004c4dd61246f8bccc2f7bd3fe03755448394 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Mon, 26 May 2025 14:44:51 -0600 Subject: [PATCH 03/41] add virtual_temperature.cpp and link it using pybind11 --- src/CMakeLists.txt | 11 ++++++++-- src/calcmod.cpp | 25 ++++++++++++++++++++-- src/virtual_temperature.cpp | 42 +++++++++++++++++++++++++++++++++++++ src/virtual_temperature.hpp | 24 +++++++++++++++++++++ 4 files changed, 98 insertions(+), 4 deletions(-) create mode 100644 src/virtual_temperature.cpp create mode 100644 src/virtual_temperature.hpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 1080051ca0..59aeba7608 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -4,12 +4,19 @@ project( LANGUAGES CXX VERSION "${SKBUILD_PROJECT_VERSION}") +# Enforce C++14 or higher +set(CMAKE_CXX_STANDARD 14) +set(CMAKE_CXX_STANDARD_REQUIRED ON) + find_package(Python COMPONENTS Interpreter Development.Module REQUIRED) find_package(pybind11 CONFIG REQUIRED) -pybind11_add_module(_calc_mod MODULE calcmod.cpp) +pybind11_add_module(_calc_mod MODULE + calcmod.cpp + virtual_temperature.cpp +) target_compile_definitions(_calc_mod PRIVATE VERSION_INFO=${PROJECT_VERSION}) -install(TARGETS _calc_mod DESTINATION metpy) \ No newline at end of file +install(TARGETS _calc_mod DESTINATION metpy) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 0d7837d01e..3c48600549 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -1,10 +1,31 @@ #include +#include +#include "virtual_temperature.hpp" + +namespace py = pybind11; int add(int i, int j) { - return i + j; + return i - j; } + PYBIND11_MODULE(_calc_mod, m) { m.doc() = "accelerator module docstring"; + m.def("add", &add, "Add two numbers"); -} \ No newline at end of file + + m.def("virtual_temperature", + py::overload_cast&, const std::vector&, double>(&VirtualTemperature), + py::arg("temperature"), py::arg("mixing_ratio"), py::arg("epsilon"), + "Compute virtual temperature from temperature and mixing ratio arrays"); + + m.def("virtual_temperature", + py::overload_cast&, double>(&VirtualTemperature), + py::arg("temperature"), py::arg("mixing_ratio"), py::arg("epsilon"), + "Compute virtual temperature from scalar temperature and mixing ratio array"); + + m.def("virtual_temperature", + py::overload_cast&, double, double>(&VirtualTemperature), + py::arg("temperature"), py::arg("mixing_ratio"), py::arg("epsilon"), + "Compute virtual temperature from temperature array and scalar mixing ratio"); +} diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp new file mode 100644 index 0000000000..b07236c872 --- /dev/null +++ b/src/virtual_temperature.cpp @@ -0,0 +1,42 @@ +#include "virtual_temperature.hpp" +#include + +std::vector VirtualTemperature( + const std::vector& temperature, + const std::vector& mixing_ratio, + double epsilon +) { + if (temperature.size() != mixing_ratio.size()) { + throw std::invalid_argument("Temperature and mixing ratio vectors must be the same size."); + } + + std::vector result; + result.reserve(temperature.size()); + + for (size_t i = 0; i < temperature.size(); ++i) { + double T = temperature[i]; + double w = mixing_ratio[i]; + result.push_back(T * (w + epsilon) / (epsilon * (1 + w))); + } + + return result; +} + +std::vector VirtualTemperature( + double temperature, + const std::vector& mixing_ratio, + double epsilon +) { + std::vector temperature_vec(mixing_ratio.size(), temperature); + return VirtualTemperature(temperature_vec, mixing_ratio, epsilon); +} + +std::vector VirtualTemperature( + const std::vector& temperature, + double mixing_ratio, + double epsilon +) { + std::vector mixing_ratio_vec(temperature.size(), mixing_ratio); + return VirtualTemperature(temperature, mixing_ratio_vec, epsilon); +} + diff --git a/src/virtual_temperature.hpp b/src/virtual_temperature.hpp new file mode 100644 index 0000000000..aaceac9055 --- /dev/null +++ b/src/virtual_temperature.hpp @@ -0,0 +1,24 @@ +#ifndef VIRTUAL_TEMPERATURE_HPP // if not defined +#define VIRTUAL_TEMPERATURE_HPP // define the header file + +#include + +// Compute virtual temperature: vector + vector +std::vector VirtualTemperature( + const std::vector& temperature, + const std::vector& mixing_ratio, + double epsilon = 0.622); + +// Overload: scalar + vector +std::vector VirtualTemperature( + double temperature, + const std::vector& mixing_ratio, + double epsilon = 0.622); + +// Overload: vector + scalar +std::vector VirtualTemperature( + const std::vector& temperature, + double mixing_ratio, + double epsilon = 0.622); + +#endif // VIRTUAL_TEMPERATURE_HPP From 358161ada6f6483ca23bffa8719ebe8edc16793a Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Tue, 27 May 2025 11:50:53 -0600 Subject: [PATCH 04/41] debug 1.input data format compatibility 2.set default epsilon 0.622 --- src/calcmod.cpp | 60 +++++++++++++++++++++++++++---------- src/virtual_temperature.cpp | 35 ++++------------------ src/virtual_temperature.hpp | 14 +-------- 3 files changed, 51 insertions(+), 58 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 3c48600549..06d4d6098a 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -1,31 +1,59 @@ +#include "virtual_temperature.hpp" #include #include -#include "virtual_temperature.hpp" namespace py = pybind11; +// Flexible dispatcher that supports scalar/list combinations +std::vector virtual_temperature_py(py::object temperature, py::object mixing_ratio, double epsilon) { + std::vector temperature_vec; + std::vector mixing_ratio_vec; + + // Case 1: both are lists + if (py::isinstance(temperature) && py::isinstance(mixing_ratio)) { + temperature_vec = temperature.cast>(); + mixing_ratio_vec = mixing_ratio.cast>(); + if (temperature_vec.size() != mixing_ratio_vec.size()) { + throw std::invalid_argument("Temperature and mixing ratio lists must be the same length."); + } + } + // Case 2: temperature is float, mixing_ratio is list + else if (py::isinstance(temperature) && py::isinstance(mixing_ratio)) { + mixing_ratio_vec = mixing_ratio.cast>(); + temperature_vec = std::vector(mixing_ratio_vec.size(), temperature.cast()); + } + // Case 3: temperature is list, mixing_ratio is float + else if (py::isinstance(temperature) && py::isinstance(mixing_ratio)) { + temperature_vec = temperature.cast>(); + mixing_ratio_vec = std::vector(temperature_vec.size(), mixing_ratio.cast()); + } + // Case 4: both are floats + else if (py::isinstance(temperature) && py::isinstance(mixing_ratio)) { + temperature_vec = {temperature.cast()}; + mixing_ratio_vec = {mixing_ratio.cast()}; + } + else { + throw std::invalid_argument("Inputs must be float or list."); + } + + return VirtualTemperature(temperature_vec, mixing_ratio_vec, epsilon); +} + int add(int i, int j) { return i - j; } - PYBIND11_MODULE(_calc_mod, m) { m.doc() = "accelerator module docstring"; m.def("add", &add, "Add two numbers"); - m.def("virtual_temperature", - py::overload_cast&, const std::vector&, double>(&VirtualTemperature), - py::arg("temperature"), py::arg("mixing_ratio"), py::arg("epsilon"), - "Compute virtual temperature from temperature and mixing ratio arrays"); - - m.def("virtual_temperature", - py::overload_cast&, double>(&VirtualTemperature), - py::arg("temperature"), py::arg("mixing_ratio"), py::arg("epsilon"), - "Compute virtual temperature from scalar temperature and mixing ratio array"); - - m.def("virtual_temperature", - py::overload_cast&, double, double>(&VirtualTemperature), - py::arg("temperature"), py::arg("mixing_ratio"), py::arg("epsilon"), - "Compute virtual temperature from temperature array and scalar mixing ratio"); + // Unified binding with default epsilon + m.def("virtual_temperature", &virtual_temperature_py, + py::arg("temperature"), py::arg("mixing_ratio"), py::arg("epsilon") = 0.622, + "Compute virtual temperature.\n" + "Accepts:\n" + " - two lists of equal length\n" + " - one scalar and one list\n" + "Defaults to epsilon = 0.622"); } diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index b07236c872..592365cce0 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -1,42 +1,19 @@ #include "virtual_temperature.hpp" -#include +//#include std::vector VirtualTemperature( const std::vector& temperature, const std::vector& mixing_ratio, double epsilon ) { - if (temperature.size() != mixing_ratio.size()) { - throw std::invalid_argument("Temperature and mixing ratio vectors must be the same size."); - } - - std::vector result; - result.reserve(temperature.size()); + std::vector result(temperature.size()); + double T, w; for (size_t i = 0; i < temperature.size(); ++i) { - double T = temperature[i]; - double w = mixing_ratio[i]; - result.push_back(T * (w + epsilon) / (epsilon * (1 + w))); + T = temperature[i]; + w = mixing_ratio[i]; + result[i] = T * (w + epsilon) / (epsilon * (1 + w)); } return result; } - -std::vector VirtualTemperature( - double temperature, - const std::vector& mixing_ratio, - double epsilon -) { - std::vector temperature_vec(mixing_ratio.size(), temperature); - return VirtualTemperature(temperature_vec, mixing_ratio, epsilon); -} - -std::vector VirtualTemperature( - const std::vector& temperature, - double mixing_ratio, - double epsilon -) { - std::vector mixing_ratio_vec(temperature.size(), mixing_ratio); - return VirtualTemperature(temperature, mixing_ratio_vec, epsilon); -} - diff --git a/src/virtual_temperature.hpp b/src/virtual_temperature.hpp index aaceac9055..2fc91564ef 100644 --- a/src/virtual_temperature.hpp +++ b/src/virtual_temperature.hpp @@ -7,18 +7,6 @@ std::vector VirtualTemperature( const std::vector& temperature, const std::vector& mixing_ratio, - double epsilon = 0.622); - -// Overload: scalar + vector -std::vector VirtualTemperature( - double temperature, - const std::vector& mixing_ratio, - double epsilon = 0.622); - -// Overload: vector + scalar -std::vector VirtualTemperature( - const std::vector& temperature, - double mixing_ratio, - double epsilon = 0.622); + double epsilon); #endif // VIRTUAL_TEMPERATURE_HPP From f18d93bd782766c8c11af6d089b8fa25f4903436 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Wed, 28 May 2025 15:55:49 -0600 Subject: [PATCH 05/41] remove setting c++14 in CMakeLists.txt --- src/CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 59aeba7608..968e79ddc7 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -5,8 +5,8 @@ project( VERSION "${SKBUILD_PROJECT_VERSION}") # Enforce C++14 or higher -set(CMAKE_CXX_STANDARD 14) -set(CMAKE_CXX_STANDARD_REQUIRED ON) +#set(CMAKE_CXX_STANDARD 14) +#set(CMAKE_CXX_STANDARD_REQUIRED ON) find_package(Python COMPONENTS Interpreter Development.Module REQUIRED) find_package(pybind11 CONFIG REQUIRED) From d12c77d9f3493b0cf1cd1627ef56bc713db27376 Mon Sep 17 00:00:00 2001 From: Linfeng Li <99068959+gitlinffff@users.noreply.github.com> Date: Tue, 3 Jun 2025 18:16:22 -0400 Subject: [PATCH 06/41] Pybind11_vectorize (#1) * vectorize the function using pybind11/numpy.h, achieving handling numpy array, float and dtype * add dewpoint, add constants.hpp, adding other functions WIP --- src/calcmod.cpp | 52 +++++++------------------------------ src/constants.hpp | 37 ++++++++++++++++++++++++++ src/virtual_temperature.cpp | 36 +++++++++++++++---------- src/virtual_temperature.hpp | 10 +++---- 4 files changed, 73 insertions(+), 62 deletions(-) create mode 100644 src/constants.hpp diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 06d4d6098a..b32aeb41be 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -1,46 +1,12 @@ #include "virtual_temperature.hpp" #include +#include #include namespace py = pybind11; -// Flexible dispatcher that supports scalar/list combinations -std::vector virtual_temperature_py(py::object temperature, py::object mixing_ratio, double epsilon) { - std::vector temperature_vec; - std::vector mixing_ratio_vec; - - // Case 1: both are lists - if (py::isinstance(temperature) && py::isinstance(mixing_ratio)) { - temperature_vec = temperature.cast>(); - mixing_ratio_vec = mixing_ratio.cast>(); - if (temperature_vec.size() != mixing_ratio_vec.size()) { - throw std::invalid_argument("Temperature and mixing ratio lists must be the same length."); - } - } - // Case 2: temperature is float, mixing_ratio is list - else if (py::isinstance(temperature) && py::isinstance(mixing_ratio)) { - mixing_ratio_vec = mixing_ratio.cast>(); - temperature_vec = std::vector(mixing_ratio_vec.size(), temperature.cast()); - } - // Case 3: temperature is list, mixing_ratio is float - else if (py::isinstance(temperature) && py::isinstance(mixing_ratio)) { - temperature_vec = temperature.cast>(); - mixing_ratio_vec = std::vector(temperature_vec.size(), mixing_ratio.cast()); - } - // Case 4: both are floats - else if (py::isinstance(temperature) && py::isinstance(mixing_ratio)) { - temperature_vec = {temperature.cast()}; - mixing_ratio_vec = {mixing_ratio.cast()}; - } - else { - throw std::invalid_argument("Inputs must be float or list."); - } - - return VirtualTemperature(temperature_vec, mixing_ratio_vec, epsilon); -} - int add(int i, int j) { - return i - j; + return i + j; } PYBIND11_MODULE(_calc_mod, m) { @@ -49,11 +15,11 @@ PYBIND11_MODULE(_calc_mod, m) { m.def("add", &add, "Add two numbers"); // Unified binding with default epsilon - m.def("virtual_temperature", &virtual_temperature_py, - py::arg("temperature"), py::arg("mixing_ratio"), py::arg("epsilon") = 0.622, - "Compute virtual temperature.\n" - "Accepts:\n" - " - two lists of equal length\n" - " - one scalar and one list\n" - "Defaults to epsilon = 0.622"); + m.def("dewpoint", py::vectorize(DewPoint), + "Calculate dew point from water vapor partial pressure.", + py::arg("vapor_pressure")); + + m.def("virtual_temperature", py::vectorize(VirtualTemperature), + "Calculate virtual temperature from temperature and mixing ratio.", + py::arg("temperature"), py::arg("mixing_ratio"), py::arg("epsilon") = 0.622); } diff --git a/src/constants.hpp b/src/constants.hpp new file mode 100644 index 0000000000..8e9d596a61 --- /dev/null +++ b/src/constants.hpp @@ -0,0 +1,37 @@ +#ifndef CONSTANTS_HPP +#define CONSTANTS_HPP + +namespace metpy_constants { +// Gas constants (J / kg / K) +constexpr double R = 8.314462618; // Universal gas constant (J / mol / K) + +// Dry air +constexpr double Md = 28.96546e-3; // Molar mass of dry air (kg / mol) +constexpr double Rd = R / Md; // Dry air (J / kg / K) +constexpr double dry_air_spec_heat_ratio = 1.4; +constexpr double Cp_d = Rd * dry_air_spec_heat_ratio / (dry_air_spec_heat_ratio - 1.0); // (J / kg / K) +constexpr double Cv_d = Cp_d / dry_air_spec_heat_ratio; // (J / kg / K) + +// Water +constexpr double Mw = 18.015268e-3; // Molar mass of water (kg / mol) +constexpr double Rv = R / Mw; // Water vapor (J / kg / K) +constexpr double wv_spec_heat_ratio = 1.33; +constexpr double Cp_v = Rv * wv_spec_heat_ratio / (wv_spec_heat_ratio - 1.0); // (J / kg / K) +constexpr double Cv_v = Cp_v / wv_spec_heat_ratio; // (J / kg / K) +constexpr double Cp_l = 4.2194e3; // Specific heat capacity of liquid water (J / kg / K) +constexpr double Lv = 2.50084e6; // Latent heat of vaporization of water (J / kg) +constexpr double T0 = 273.16; // Triple point of water (K) + +// General Meteorological constants +constexpr double epsilon = Mw / Md; // ≈ 0.622 + + +// Standard gravity +constexpr double g = 9.80665; // m / s^2 + +// Reference pressure +constexpr double P0 = 100000.0; // Pa (hPa = 1000) + +} + +#endif // CONSTANTS_HPP diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index 592365cce0..0204a07066 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -1,19 +1,29 @@ +#include +#include "constants.hpp" #include "virtual_temperature.hpp" //#include -std::vector VirtualTemperature( - const std::vector& temperature, - const std::vector& mixing_ratio, - double epsilon -) { +double water_latent_heat_vaporization(double temperature) { + // Calculate the latent heat of vaporization of water in J/kg at a given temperature. + using namespace metpy_constants; + return Lv - (Cp_l - Cp_v) * (temperature - T0); +} + +double _saturation_vapor_pressure(double temperature) { + // Calculate saturation (equilibrium) water vapor (partial) pressure over liquid water. + // Constants for the Magnus-Tetens approximation + //const double a = 17.67; + //const double b = 243.5; - std::vector result(temperature.size()); - double T, w; - for (size_t i = 0; i < temperature.size(); ++i) { - T = temperature[i]; - w = mixing_ratio[i]; - result[i] = T * (w + epsilon) / (epsilon * (1 + w)); - } + // Calculate saturation vapor pressure using the Magnus-Tetens formula + //return 6.112 * exp((a * temperature) / (b + temperature)); +} + +double DewPoint(double vapor_pressure) { + double val = log(vapor_pressure / 6.112); + return 243.5 * val / (17.67 - val); +} - return result; +double VirtualTemperature(double temperature, double mixing_ratio, double epsilon) { + return temperature * (mixing_ratio + epsilon) / (epsilon * (1. + mixing_ratio)); } diff --git a/src/virtual_temperature.hpp b/src/virtual_temperature.hpp index 2fc91564ef..4b871bebd0 100644 --- a/src/virtual_temperature.hpp +++ b/src/virtual_temperature.hpp @@ -1,12 +1,10 @@ #ifndef VIRTUAL_TEMPERATURE_HPP // if not defined #define VIRTUAL_TEMPERATURE_HPP // define the header file -#include +//#include -// Compute virtual temperature: vector + vector -std::vector VirtualTemperature( - const std::vector& temperature, - const std::vector& mixing_ratio, - double epsilon); +double DewPoint(double vapor_pressure); + +double VirtualTemperature(double temperature, double mixing_ratio, double epsilon); #endif // VIRTUAL_TEMPERATURE_HPP From 31dbf18284b5422ff0cc339764c023300492eb95 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Tue, 3 Jun 2025 18:16:47 -0600 Subject: [PATCH 07/41] wip --- src/calcmod.cpp | 3 +++ src/constants.cpp | 16 ++++++++++++++++ src/constants.hpp | 34 ++++------------------------------ src/constants_backup | 37 +++++++++++++++++++++++++++++++++++++ src/virtual_temperature.cpp | 14 +++++++++++--- 5 files changed, 71 insertions(+), 33 deletions(-) create mode 100644 src/constants.cpp create mode 100644 src/constants_backup diff --git a/src/calcmod.cpp b/src/calcmod.cpp index b32aeb41be..baf21ddc28 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -1,3 +1,4 @@ +#include "constants.hpp" #include "virtual_temperature.hpp" #include #include @@ -12,6 +13,8 @@ int add(int i, int j) { PYBIND11_MODULE(_calc_mod, m) { m.doc() = "accelerator module docstring"; + metpy_constants::load_constants_from_python(); + m.def("add", &add, "Add two numbers"); // Unified binding with default epsilon diff --git a/src/constants.cpp b/src/constants.cpp new file mode 100644 index 0000000000..7f54c29d35 --- /dev/null +++ b/src/constants.cpp @@ -0,0 +1,16 @@ +// constants.cpp +#include "constants.hpp" +#include + +namespace py = pybind11; +using namespace metpy_constants; + +// Definitions (required in C++11) +double sat_pressure_0c; + +void load_constants_from_python() { + py::object mod = py::module_::import("metpy.constants.default"); + + sat_pressure_0c = mod.attr("sat_pressure_0c").attr("magnitude").cast(); +} + diff --git a/src/constants.hpp b/src/constants.hpp index 8e9d596a61..fef19add59 100644 --- a/src/constants.hpp +++ b/src/constants.hpp @@ -2,36 +2,10 @@ #define CONSTANTS_HPP namespace metpy_constants { -// Gas constants (J / kg / K) -constexpr double R = 8.314462618; // Universal gas constant (J / mol / K) - -// Dry air -constexpr double Md = 28.96546e-3; // Molar mass of dry air (kg / mol) -constexpr double Rd = R / Md; // Dry air (J / kg / K) -constexpr double dry_air_spec_heat_ratio = 1.4; -constexpr double Cp_d = Rd * dry_air_spec_heat_ratio / (dry_air_spec_heat_ratio - 1.0); // (J / kg / K) -constexpr double Cv_d = Cp_d / dry_air_spec_heat_ratio; // (J / kg / K) - -// Water -constexpr double Mw = 18.015268e-3; // Molar mass of water (kg / mol) -constexpr double Rv = R / Mw; // Water vapor (J / kg / K) -constexpr double wv_spec_heat_ratio = 1.33; -constexpr double Cp_v = Rv * wv_spec_heat_ratio / (wv_spec_heat_ratio - 1.0); // (J / kg / K) -constexpr double Cv_v = Cp_v / wv_spec_heat_ratio; // (J / kg / K) -constexpr double Cp_l = 4.2194e3; // Specific heat capacity of liquid water (J / kg / K) -constexpr double Lv = 2.50084e6; // Latent heat of vaporization of water (J / kg) -constexpr double T0 = 273.16; // Triple point of water (K) - -// General Meteorological constants -constexpr double epsilon = Mw / Md; // ≈ 0.622 - - -// Standard gravity -constexpr double g = 9.80665; // m / s^2 - -// Reference pressure -constexpr double P0 = 100000.0; // Pa (hPa = 1000) + extern double sat_pressure_0c; + void load_constants_from_python(); // call once in your PYBIND11_MODULE } -#endif // CONSTANTS_HPP +#endif + diff --git a/src/constants_backup b/src/constants_backup new file mode 100644 index 0000000000..8e9d596a61 --- /dev/null +++ b/src/constants_backup @@ -0,0 +1,37 @@ +#ifndef CONSTANTS_HPP +#define CONSTANTS_HPP + +namespace metpy_constants { +// Gas constants (J / kg / K) +constexpr double R = 8.314462618; // Universal gas constant (J / mol / K) + +// Dry air +constexpr double Md = 28.96546e-3; // Molar mass of dry air (kg / mol) +constexpr double Rd = R / Md; // Dry air (J / kg / K) +constexpr double dry_air_spec_heat_ratio = 1.4; +constexpr double Cp_d = Rd * dry_air_spec_heat_ratio / (dry_air_spec_heat_ratio - 1.0); // (J / kg / K) +constexpr double Cv_d = Cp_d / dry_air_spec_heat_ratio; // (J / kg / K) + +// Water +constexpr double Mw = 18.015268e-3; // Molar mass of water (kg / mol) +constexpr double Rv = R / Mw; // Water vapor (J / kg / K) +constexpr double wv_spec_heat_ratio = 1.33; +constexpr double Cp_v = Rv * wv_spec_heat_ratio / (wv_spec_heat_ratio - 1.0); // (J / kg / K) +constexpr double Cv_v = Cp_v / wv_spec_heat_ratio; // (J / kg / K) +constexpr double Cp_l = 4.2194e3; // Specific heat capacity of liquid water (J / kg / K) +constexpr double Lv = 2.50084e6; // Latent heat of vaporization of water (J / kg) +constexpr double T0 = 273.16; // Triple point of water (K) + +// General Meteorological constants +constexpr double epsilon = Mw / Md; // ≈ 0.622 + + +// Standard gravity +constexpr double g = 9.80665; // m / s^2 + +// Reference pressure +constexpr double P0 = 100000.0; // Pa (hPa = 1000) + +} + +#endif // CONSTANTS_HPP diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index 0204a07066..2cc5cccbed 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -1,12 +1,15 @@ #include +#include #include "constants.hpp" #include "virtual_temperature.hpp" //#include +namespace py = pybind11; + double water_latent_heat_vaporization(double temperature) { // Calculate the latent heat of vaporization of water in J/kg at a given temperature. - using namespace metpy_constants; - return Lv - (Cp_l - Cp_v) * (temperature - T0); +// using namespace metpy_constants; +// return Lv - (Cp_l - Cp_v) * (temperature - T0); } double _saturation_vapor_pressure(double temperature) { @@ -20,7 +23,12 @@ double _saturation_vapor_pressure(double temperature) { } double DewPoint(double vapor_pressure) { - double val = log(vapor_pressure / 6.112); + // fetch constants from python module + //py::object default_mod = py::module_::import("metpy.constants.default"); + // unit issue ignored for now + //double sat_pressure_0c = default_mod.attr("sat_pressure_0c").attr("magnitude").cast(); + + double val = log(vapor_pressure / sat_pressure_0c); return 243.5 * val / (17.67 - val); } From 8e1ecccb44c8f247882def973111a03856f10fd8 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Wed, 4 Jun 2025 17:55:31 -0600 Subject: [PATCH 08/41] create constants.cpp and .hpp to fetch all constants from python; converting saturation vapor pressure WIP --- src/CMakeLists.txt | 1 + src/calcmod.cpp | 8 ++++++++ src/constants.cpp | 25 ++++++++++++++++++------- src/constants.hpp | 31 +++++++++++++++++++++++++++++++ src/virtual_temperature.cpp | 25 +++++++++---------------- src/virtual_temperature.hpp | 2 ++ 6 files changed, 69 insertions(+), 23 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 968e79ddc7..c6d61218a2 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -13,6 +13,7 @@ find_package(pybind11 CONFIG REQUIRED) pybind11_add_module(_calc_mod MODULE calcmod.cpp + constants.cpp virtual_temperature.cpp ) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index baf21ddc28..573a1ad35d 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -17,6 +17,14 @@ PYBIND11_MODULE(_calc_mod, m) { m.def("add", &add, "Add two numbers"); + m.def("water_latent_heat_vaporization", py::vectorize(WaterLatentHeatVaporization), + "Calculate water latent heat vaporization from temperature.", + py::arg("temperature")); + + m.def("_saturation_vapor_pressure_liquid", py::vectorize(_SaturationVaporPressureLiquid), + "Calculate saturation vapor pressure from temperature.", + py::arg("temperature")); + // Unified binding with default epsilon m.def("dewpoint", py::vectorize(DewPoint), "Calculate dew point from water vapor partial pressure.", diff --git a/src/constants.cpp b/src/constants.cpp index 7f54c29d35..d99bfd3bba 100644 --- a/src/constants.cpp +++ b/src/constants.cpp @@ -3,14 +3,25 @@ #include namespace py = pybind11; -using namespace metpy_constants; -// Definitions (required in C++11) -double sat_pressure_0c; +namespace metpy_constants { + double Mw; + double Rv; + double Cp_v; + double Cp_l; + double Lv; + double sat_pressure_0c; + double T0; -void load_constants_from_python() { - py::object mod = py::module_::import("metpy.constants.default"); + void load_constants_from_python() { + py::object mod = py::module_::import("metpy.constants.default"); - sat_pressure_0c = mod.attr("sat_pressure_0c").attr("magnitude").cast(); + Mw = mod.attr("Mw").attr("to")("kg / mol").attr("magnitude").cast(); + Rv = mod.attr("Rv").attr("magnitude").cast(); + Cp_v = mod.attr("Cp_v").attr("magnitude").cast(); + Cp_l = mod.attr("Cp_l").attr("magnitude").cast(); + Lv = mod.attr("Lv").attr("magnitude").cast(); + sat_pressure_0c = mod.attr("sat_pressure_0c").attr("magnitude").cast(); + T0 = mod.attr("T0").attr("magnitude").cast(); + } } - diff --git a/src/constants.hpp b/src/constants.hpp index fef19add59..95ad4cb856 100644 --- a/src/constants.hpp +++ b/src/constants.hpp @@ -2,7 +2,38 @@ #define CONSTANTS_HPP namespace metpy_constants { + + // Gas constants (J / kg / K) +// extern double R; + + // Dry air +// extern double Md; +// extern double Rd; +// extern double dry_air_spec_heat_ratio; +// extern double Cp_d = Rd * dry_air_spec_heat_ratio / (dry_air_spec_heat_ratio - 1.0); // (J / kg / K) +// extern double Cv_d = Cp_d / dry_air_spec_heat_ratio; // (J / kg / K) + + // Water + extern double Mw; + extern double Rv; +// extern double wv_spec_heat_ratio = 1.33; + extern double Cp_v; +// extern double Cv_v = Cp_v / wv_spec_heat_ratio; // (J / kg / K) + extern double Cp_l; + extern double Lv; extern double sat_pressure_0c; + extern double T0; + + // General Meteorological constants +// extern double epsilon = Mw / Md; // ≈ 0.622 + + + // Standard gravity +// extern double g = 9.80665; // m / s^2 + + // Reference pressure +// extern double P0 = 100000.0; // Pa (hPa = 1000) + void load_constants_from_python(); // call once in your PYBIND11_MODULE } diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index 2cc5cccbed..b296043144 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -5,30 +5,23 @@ //#include namespace py = pybind11; +namespace mc = metpy_constants; -double water_latent_heat_vaporization(double temperature) { - // Calculate the latent heat of vaporization of water in J/kg at a given temperature. -// using namespace metpy_constants; -// return Lv - (Cp_l - Cp_v) * (temperature - T0); +double WaterLatentHeatVaporization(double temperature) { + return mc::Lv - (mc::Cp_l - mc::Cp_v) * (temperature - mc::T0); } -double _saturation_vapor_pressure(double temperature) { - // Calculate saturation (equilibrium) water vapor (partial) pressure over liquid water. - // Constants for the Magnus-Tetens approximation - //const double a = 17.67; - //const double b = 243.5; +double _SaturationVaporPressureLiquid(double temperature) { + double latent_heat = WaterLatentHeatVaporization(temperature); + double heat_power = (mc::Cp_l - mc::Cp_v) / mc::Rv; + double exp_term = (mc::Lv / mc::T0 - latent_heat / temperature) / mc::Rv; - // Calculate saturation vapor pressure using the Magnus-Tetens formula - //return 6.112 * exp((a * temperature) / (b + temperature)); + return mc::sat_pressure_0c * exp(exp_term) * pow(mc::T0 / temperature, heat_power); } double DewPoint(double vapor_pressure) { - // fetch constants from python module - //py::object default_mod = py::module_::import("metpy.constants.default"); - // unit issue ignored for now - //double sat_pressure_0c = default_mod.attr("sat_pressure_0c").attr("magnitude").cast(); - double val = log(vapor_pressure / sat_pressure_0c); + double val = log(vapor_pressure / mc::sat_pressure_0c); return 243.5 * val / (17.67 - val); } diff --git a/src/virtual_temperature.hpp b/src/virtual_temperature.hpp index 4b871bebd0..230b3737f6 100644 --- a/src/virtual_temperature.hpp +++ b/src/virtual_temperature.hpp @@ -2,7 +2,9 @@ #define VIRTUAL_TEMPERATURE_HPP // define the header file //#include +double WaterLatentHeatVaporization(double temperature); +double _SaturationVaporPressureLiquid(double temperature); double DewPoint(double vapor_pressure); double VirtualTemperature(double temperature, double mixing_ratio, double epsilon); From 6b81f0e89468e6dc6db079238898ff01e35b0fac Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Thu, 5 Jun 2025 16:17:53 -0600 Subject: [PATCH 09/41] c++ fetch nounit constants from nounit.py --- src/constants.cpp | 25 +++++++++++++++++-------- src/metpy/constants/nounit.py | 1 + src/virtual_temperature.cpp | 2 +- 3 files changed, 19 insertions(+), 9 deletions(-) diff --git a/src/constants.cpp b/src/constants.cpp index d99bfd3bba..55ebbbe656 100644 --- a/src/constants.cpp +++ b/src/constants.cpp @@ -14,14 +14,23 @@ namespace metpy_constants { double T0; void load_constants_from_python() { - py::object mod = py::module_::import("metpy.constants.default"); + py::object mod = py::module_::import("metpy.constants.nounit"); - Mw = mod.attr("Mw").attr("to")("kg / mol").attr("magnitude").cast(); - Rv = mod.attr("Rv").attr("magnitude").cast(); - Cp_v = mod.attr("Cp_v").attr("magnitude").cast(); - Cp_l = mod.attr("Cp_l").attr("magnitude").cast(); - Lv = mod.attr("Lv").attr("magnitude").cast(); - sat_pressure_0c = mod.attr("sat_pressure_0c").attr("magnitude").cast(); - T0 = mod.attr("T0").attr("magnitude").cast(); +// Mw = mod.attr("Mw").attr("magnitude").cast(); +// Rv = mod.attr("Rv").attr("magnitude").cast(); +// Cp_v = mod.attr("Cp_v").attr("magnitude").cast(); +// Cp_l = mod.attr("Cp_l").attr("magnitude").cast(); +// Lv = mod.attr("Lv").attr("magnitude").cast(); +// sat_pressure_0c = mod.attr("sat_pressure_0c").cast(); +// T0 = mod.attr("T0").attr("magnitude").cast(); + + + Mw = mod.attr("Mw").cast(); + Rv = mod.attr("Rv").cast(); + Cp_v = mod.attr("Cp_v").cast(); + Cp_l = mod.attr("Cp_l").cast(); + Lv = mod.attr("Lv").cast(); + sat_pressure_0c = mod.attr("sat_pressure_0c").cast(); + T0 = mod.attr("T0").cast(); } } diff --git a/src/metpy/constants/nounit.py b/src/metpy/constants/nounit.py index 2fd6e8a78f..d1b5e83730 100644 --- a/src/metpy/constants/nounit.py +++ b/src/metpy/constants/nounit.py @@ -6,6 +6,7 @@ from . import default from ..units import units +Mw = default.Mw.m_as('kg / mol') Rd = default.Rd.m_as('m**2 / K / s**2') Rv = default.Rv.m_as('m**2 / K / s**2') Lv = default.Lv.m_as('m**2 / s**2') diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index b296043144..273e1858cd 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -22,7 +22,7 @@ double _SaturationVaporPressureLiquid(double temperature) { double DewPoint(double vapor_pressure) { double val = log(vapor_pressure / mc::sat_pressure_0c); - return 243.5 * val / (17.67 - val); + return 243.5 * val / (17.67 - val); // use SI units instead } double VirtualTemperature(double temperature, double mixing_ratio, double epsilon) { From 0cb07ba3462b829aa97730407320475d0565c260 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Fri, 6 Jun 2025 11:15:17 -0600 Subject: [PATCH 10/41] finish converting saturation vapor pressure --- src/calcmod.cpp | 10 +++++++++- src/constants.cpp | 4 ++++ src/constants.hpp | 2 ++ src/virtual_temperature.cpp | 32 +++++++++++++++++++++++++++++++- src/virtual_temperature.hpp | 6 ++++++ 5 files changed, 52 insertions(+), 2 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 573a1ad35d..2732d1b1cc 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -20,11 +20,19 @@ PYBIND11_MODULE(_calc_mod, m) { m.def("water_latent_heat_vaporization", py::vectorize(WaterLatentHeatVaporization), "Calculate water latent heat vaporization from temperature.", py::arg("temperature")); - + + m.def("water_latent_heat_sublimation", py::vectorize(WaterLatentHeatSublimation), + "Calculate water latent heat sublimation from temperature.", + py::arg("temperature")); + m.def("_saturation_vapor_pressure_liquid", py::vectorize(_SaturationVaporPressureLiquid), "Calculate saturation vapor pressure from temperature.", py::arg("temperature")); + m.def("saturation_vapor_pressure", py::vectorize(SaturationVaporPressure), + "Calculate saturation vapor pressure from temperature.", + py::arg("temperature"), py::arg("phase") = "liquid"); + // Unified binding with default epsilon m.def("dewpoint", py::vectorize(DewPoint), "Calculate dew point from water vapor partial pressure.", diff --git a/src/constants.cpp b/src/constants.cpp index 55ebbbe656..978eb93f9c 100644 --- a/src/constants.cpp +++ b/src/constants.cpp @@ -12,6 +12,8 @@ namespace metpy_constants { double Lv; double sat_pressure_0c; double T0; + double Ls; + double Cp_i; void load_constants_from_python() { py::object mod = py::module_::import("metpy.constants.nounit"); @@ -32,5 +34,7 @@ namespace metpy_constants { Lv = mod.attr("Lv").cast(); sat_pressure_0c = mod.attr("sat_pressure_0c").cast(); T0 = mod.attr("T0").cast(); + Ls = mod.attr("Ls").cast(); + Cp_i = mod.attr("Cp_i").cast(); } } diff --git a/src/constants.hpp b/src/constants.hpp index 95ad4cb856..08d76e07e8 100644 --- a/src/constants.hpp +++ b/src/constants.hpp @@ -23,6 +23,8 @@ namespace metpy_constants { extern double Lv; extern double sat_pressure_0c; extern double T0; + extern double Ls; + extern double Cp_i; // General Meteorological constants // extern double epsilon = Mw / Md; // ≈ 0.622 diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index 273e1858cd..1c20b8670c 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -1,8 +1,9 @@ #include +#include #include #include "constants.hpp" #include "virtual_temperature.hpp" -//#include +#include namespace py = pybind11; namespace mc = metpy_constants; @@ -11,6 +12,10 @@ double WaterLatentHeatVaporization(double temperature) { return mc::Lv - (mc::Cp_l - mc::Cp_v) * (temperature - mc::T0); } +double WaterLatentHeatSublimation(double temperature) { + return mc::Ls - (mc::Cp_i - mc::Cp_v) * (temperature - mc::T0); +} + double _SaturationVaporPressureLiquid(double temperature) { double latent_heat = WaterLatentHeatVaporization(temperature); double heat_power = (mc::Cp_l - mc::Cp_v) / mc::Rv; @@ -19,6 +24,31 @@ double _SaturationVaporPressureLiquid(double temperature) { return mc::sat_pressure_0c * exp(exp_term) * pow(mc::T0 / temperature, heat_power); } +double _SaturationVaporPressureSolid(double temperature) { + double latent_heat = WaterLatentHeatSublimation(temperature); + double heat_power = (mc::Cp_i - mc::Cp_v) / mc::Rv; + double exp_term = (mc::Ls / mc::T0 - latent_heat / temperature) / mc::Rv; + + return mc::sat_pressure_0c * exp(exp_term) * pow(mc::T0 / temperature, heat_power); +} + +double SaturationVaporPressure(double temperature, std::string phase) { + if (phase == "liquid") { + return _SaturationVaporPressureLiquid(temperature); + } else if (phase == "solid") { + return _SaturationVaporPressureSolid(temperature); + } else if (phase == "auto") { + if (temperature > mc::T0) { + return _SaturationVaporPressureLiquid(temperature); + } else { + return _SaturationVaporPressureSolid(temperature); + } + } else { + throw std::invalid_argument("'" + phase + "' is not a valid option for phase. " + "Valid options are 'liquid', 'solid', or 'auto'."); + } +} + double DewPoint(double vapor_pressure) { double val = log(vapor_pressure / mc::sat_pressure_0c); diff --git a/src/virtual_temperature.hpp b/src/virtual_temperature.hpp index 230b3737f6..d8182c2d0c 100644 --- a/src/virtual_temperature.hpp +++ b/src/virtual_temperature.hpp @@ -2,9 +2,15 @@ #define VIRTUAL_TEMPERATURE_HPP // define the header file //#include +#include + double WaterLatentHeatVaporization(double temperature); +double WaterLatentHeatSublimation(double temperature); double _SaturationVaporPressureLiquid(double temperature); +double _SaturationVaporPressureSolid(double temperature); +double SaturationVaporPressure(double temperature, std::string phase); + double DewPoint(double vapor_pressure); double VirtualTemperature(double temperature, double mixing_ratio, double epsilon); From 5514401f542347feb2779abb1c9e9d5c4980bdd8 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Fri, 6 Jun 2025 15:08:06 -0600 Subject: [PATCH 11/41] calc.thermo calls c++ function WIP --- src/calcmod.cpp | 2 +- src/metpy/calc/thermo.py | 12 +++++++++--- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 2732d1b1cc..c3a6ad866d 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -40,5 +40,5 @@ PYBIND11_MODULE(_calc_mod, m) { m.def("virtual_temperature", py::vectorize(VirtualTemperature), "Calculate virtual temperature from temperature and mixing ratio.", - py::arg("temperature"), py::arg("mixing_ratio"), py::arg("epsilon") = 0.622); + py::arg("temperature"), py::arg("mixing_ratio"), py::arg("epsilon")); } diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 40274192b2..d1ee282828 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -16,6 +16,7 @@ import scipy.optimize as so from scipy.special import lambertw import xarray as xr +import metpy._calc_mod as _calc_mod from .exceptions import InvalidSoundingError from .tools import (_greater_or_close, _less_or_close, _remove_nans, find_bounding_indices, @@ -2123,9 +2124,14 @@ def virtual_temperature( Renamed ``mixing`` parameter to ``mixing_ratio`` """ - return temperature * ((mixing_ratio + molecular_weight_ratio) - / (molecular_weight_ratio * (1 + mixing_ratio))) - + print(temperature) + print(type(temperature)) + exit() + T_mag = temperature.to('K').magnitude + w_mag = mixing_ratio.to('').magnitude + + return _calc_mod.virtual_temperature( + T_mag, w_mag, molecular_weight_ratio) * temperature.units @exporter.export @preprocess_and_wrap(wrap_like='temperature', broadcast=('pressure', From 697afd6eccb8ea1f87d9e8f1f95f14e398d154a9 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Fri, 6 Jun 2025 16:23:55 -0600 Subject: [PATCH 12/41] apply c++ functions in thermo.py --- src/calcmod.cpp | 4 ++-- src/metpy/calc/thermo.py | 47 +++++++++------------------------------- 2 files changed, 12 insertions(+), 39 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index c3a6ad866d..d773b9fb85 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -29,9 +29,9 @@ PYBIND11_MODULE(_calc_mod, m) { "Calculate saturation vapor pressure from temperature.", py::arg("temperature")); - m.def("saturation_vapor_pressure", py::vectorize(SaturationVaporPressure), + m.def("_saturation_vapor_pressure_solid", py::vectorize(_SaturationVaporPressureSolid), "Calculate saturation vapor pressure from temperature.", - py::arg("temperature"), py::arg("phase") = "liquid"); + py::arg("temperature")); // Unified binding with default epsilon m.def("dewpoint", py::vectorize(DewPoint), diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index d1ee282828..566b7d6629 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -186,10 +186,8 @@ def water_latent_heat_vaporization(temperature): Eq 15, [Ambaum2020]_, using MetPy-defined constants in place of cited values. """ - return (mpconsts.nounit.Lv - - (mpconsts.nounit.Cp_l - mpconsts.nounit.Cp_v) - * (temperature - mpconsts.nounit.T0)) - + # Calling c++ calculation module + return _calc_mod.water_latent_heat_vaporization(temperature) @exporter.export @preprocess_and_wrap(wrap_like='temperature') @@ -228,10 +226,8 @@ def water_latent_heat_sublimation(temperature): Eq 18, [Ambaum2020]_, using MetPy-defined constants in place of cited values. """ - return (mpconsts.nounit.Ls - - (mpconsts.nounit.Cp_i - mpconsts.nounit.Cp_v) - * (temperature - mpconsts.nounit.T0)) - + # Calling c++ calculation module + return _calc_mod.water_latent_heat_sublimation(temperature) @exporter.export @preprocess_and_wrap(wrap_like='temperature') @@ -1624,17 +1620,8 @@ def _saturation_vapor_pressure_liquid(temperature): .. math:: e = e_{s0} \frac{T_0}{T}^{(c_{pl} - c_{pv}) / R_v} \exp{ \frac{L_0}{R_v T_0} - \frac{L}{R_v T}} """ - latent_heat = water_latent_heat_vaporization._nounit(temperature) - heat_power = (mpconsts.nounit.Cp_l - mpconsts.nounit.Cp_v) / mpconsts.nounit.Rv - exp_term = ((mpconsts.nounit.Lv / mpconsts.nounit.T0 - latent_heat / temperature) - / mpconsts.nounit.Rv) - - return ( - mpconsts.nounit.sat_pressure_0c - * (mpconsts.nounit.T0 / temperature) ** heat_power - * np.exp(exp_term) - ) - + # Calling c++ calculation module + return _calc_mod._saturation_vapor_pressure_liquid(temperature) @preprocess_and_wrap(wrap_like='temperature') @process_units({'temperature': '[temperature]'}, '[pressure]') @@ -1661,17 +1648,8 @@ def _saturation_vapor_pressure_solid(temperature): .. math:: e_i = e_{i0} \frac{T_0}{T}^{(c_{pi} - c_{pv}) / R_v} \exp{ \frac{L_{s0}}{R_v T_0} - \frac{L_s}{R_v T}} """ - latent_heat = water_latent_heat_sublimation._nounit(temperature) - heat_power = (mpconsts.nounit.Cp_i - mpconsts.nounit.Cp_v) / mpconsts.nounit.Rv - exp_term = ((mpconsts.nounit.Ls / mpconsts.nounit.T0 - latent_heat / temperature) - / mpconsts.nounit.Rv) - - return ( - mpconsts.nounit.sat_pressure_0c - * (mpconsts.nounit.T0 / temperature) ** heat_power - * np.exp(exp_term) - ) - + # Calling c++ calculation module + return _calc_mod._saturation_vapor_pressure_solid(temperature) @exporter.export @preprocess_and_wrap(wrap_like='temperature', broadcast=('temperature', 'relative_humidity')) @@ -2124,14 +2102,9 @@ def virtual_temperature( Renamed ``mixing`` parameter to ``mixing_ratio`` """ - print(temperature) - print(type(temperature)) - exit() - T_mag = temperature.to('K').magnitude - w_mag = mixing_ratio.to('').magnitude - + # Calling c++ calculation module return _calc_mod.virtual_temperature( - T_mag, w_mag, molecular_weight_ratio) * temperature.units + temperature, mixing_ratio, molecular_weight_ratio) @exporter.export @preprocess_and_wrap(wrap_like='temperature', broadcast=('pressure', From 5072b95bb1e30a3472a29d61b43ec3480c56d1e3 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Fri, 6 Jun 2025 16:28:55 -0600 Subject: [PATCH 13/41] wip --- src/calcmod.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index d773b9fb85..7d4f51e94c 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -32,6 +32,10 @@ PYBIND11_MODULE(_calc_mod, m) { m.def("_saturation_vapor_pressure_solid", py::vectorize(_SaturationVaporPressureSolid), "Calculate saturation vapor pressure from temperature.", py::arg("temperature")); + + m.def("saturation_vapor_pressure", py::vectorize(SaturationVaporPressure), + "Calculate saturation vapor pressure from temperature.", + py::arg("temperature"), py::arg("phase") = "liquid"); // Unified binding with default epsilon m.def("dewpoint", py::vectorize(DewPoint), From 40c94190be70681f336dcd3f1a78ce971677ae8c Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Mon, 9 Jun 2025 09:34:55 -0600 Subject: [PATCH 14/41] wip --- src/constants.cpp | 2 ++ src/constants.hpp | 1 + src/metpy/calc/thermo.py | 3 +-- src/virtual_temperature.cpp | 2 +- 4 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/constants.cpp b/src/constants.cpp index 978eb93f9c..1a64786a30 100644 --- a/src/constants.cpp +++ b/src/constants.cpp @@ -14,6 +14,7 @@ namespace metpy_constants { double T0; double Ls; double Cp_i; + double zero_degc; void load_constants_from_python() { py::object mod = py::module_::import("metpy.constants.nounit"); @@ -36,5 +37,6 @@ namespace metpy_constants { T0 = mod.attr("T0").cast(); Ls = mod.attr("Ls").cast(); Cp_i = mod.attr("Cp_i").cast(); + zero_degc = mod.attr("zero_degc").cast(); } } diff --git a/src/constants.hpp b/src/constants.hpp index 08d76e07e8..4d2f0bc7a4 100644 --- a/src/constants.hpp +++ b/src/constants.hpp @@ -25,6 +25,7 @@ namespace metpy_constants { extern double T0; extern double Ls; extern double Cp_i; + extern double zero_degc; // General Meteorological constants // extern double epsilon = Mw / Md; // ≈ 0.622 diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 566b7d6629..0b0dc786f8 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -1729,9 +1729,8 @@ def dewpoint(vapor_pressure): Renamed ``e`` parameter to ``vapor_pressure`` """ - val = np.log(vapor_pressure / mpconsts.nounit.sat_pressure_0c) - return mpconsts.nounit.zero_degc + 243.5 * val / (17.67 - val) + return _calc_mod.dewpoint(vapor_pressure) @exporter.export @preprocess_and_wrap(wrap_like='partial_press', broadcast=('partial_press', 'total_press')) diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index 1c20b8670c..29a94a8a3e 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -52,7 +52,7 @@ double SaturationVaporPressure(double temperature, std::string phase) { double DewPoint(double vapor_pressure) { double val = log(vapor_pressure / mc::sat_pressure_0c); - return 243.5 * val / (17.67 - val); // use SI units instead + return mc::zero_degc + 243.5 * val / (17.67 - val); } double VirtualTemperature(double temperature, double mixing_ratio, double epsilon) { From f913ea7ab587ae5830399817e763df688f934ef8 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Mon, 9 Jun 2025 15:47:40 -0600 Subject: [PATCH 15/41] wip --- src/calcmod.cpp | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 7d4f51e94c..7f3491da79 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -25,6 +25,10 @@ PYBIND11_MODULE(_calc_mod, m) { "Calculate water latent heat sublimation from temperature.", py::arg("temperature")); + m.def("saturation_vapor_pressure", py::vectorize(SaturationVaporPressure), + "Calculate saturation vapor pressure from temperature.", + py::arg("temperature"), py::arg("phase") = "liquid"); + m.def("_saturation_vapor_pressure_liquid", py::vectorize(_SaturationVaporPressureLiquid), "Calculate saturation vapor pressure from temperature.", py::arg("temperature")); @@ -33,11 +37,6 @@ PYBIND11_MODULE(_calc_mod, m) { "Calculate saturation vapor pressure from temperature.", py::arg("temperature")); - m.def("saturation_vapor_pressure", py::vectorize(SaturationVaporPressure), - "Calculate saturation vapor pressure from temperature.", - py::arg("temperature"), py::arg("phase") = "liquid"); - - // Unified binding with default epsilon m.def("dewpoint", py::vectorize(DewPoint), "Calculate dew point from water vapor partial pressure.", py::arg("vapor_pressure")); From e7d870ef9f09c9622a1aff2c9a4f7cdcf2f438fb Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Wed, 11 Jun 2025 16:27:44 -0600 Subject: [PATCH 16/41] add c++ functions: mixing_ratio, saturation_mixing_ratio, specific_humidity_from_mixing_ratio, specific_humidity_from_dewpoint --- src/calcmod.cpp | 17 +++++++++++++++++ src/constants.cpp | 2 ++ src/constants.hpp | 1 + src/metpy/calc/thermo.py | 5 +++-- src/virtual_temperature.cpp | 25 +++++++++++++++++++++++++ src/virtual_temperature.hpp | 7 +++++++ 6 files changed, 55 insertions(+), 2 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 7f3491da79..7fe49686b7 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -41,6 +41,23 @@ PYBIND11_MODULE(_calc_mod, m) { "Calculate dew point from water vapor partial pressure.", py::arg("vapor_pressure")); + m.def("mixing_ratio", py::vectorize(MixingRatio), + "Calculate the mixing ratio of a gas.", + py::arg("partial_press"), py::arg("total_press"), py::arg("epsilon")); + + m.def("saturation_mixing_ratio", py::vectorize(SaturationMixingRatio), + "Calculate the saturation mixing ratio of water vapor given total atmospheric pressure and temperature.", + py::arg("total_press"), py::arg("temperature"), py::arg("phase")); + + m.def("specific_humidity_from_mixing_ratio", py::vectorize(SpecificHumidityFromMixingRatio), + "Calculate the specific humidity from the mixing ratio.", + py::arg("mixing_ratio")); + + m.def("specific_humidity_from_dewpoint", py::vectorize(SpecificHumidityFromDewPoint), + "Calculate the specific humidity from the dewpoint temperature and pressure.", + py::arg("pressure"), py::arg("dewpoint"), py::arg("phase")); + + m.def("virtual_temperature", py::vectorize(VirtualTemperature), "Calculate virtual temperature from temperature and mixing ratio.", py::arg("temperature"), py::arg("mixing_ratio"), py::arg("epsilon")); diff --git a/src/constants.cpp b/src/constants.cpp index 1a64786a30..a83a7ed737 100644 --- a/src/constants.cpp +++ b/src/constants.cpp @@ -15,6 +15,7 @@ namespace metpy_constants { double Ls; double Cp_i; double zero_degc; + double epsilon; void load_constants_from_python() { py::object mod = py::module_::import("metpy.constants.nounit"); @@ -38,5 +39,6 @@ namespace metpy_constants { Ls = mod.attr("Ls").cast(); Cp_i = mod.attr("Cp_i").cast(); zero_degc = mod.attr("zero_degc").cast(); + epsilon = mod.attr("epsilon").cast(); } } diff --git a/src/constants.hpp b/src/constants.hpp index 4d2f0bc7a4..4707c7a046 100644 --- a/src/constants.hpp +++ b/src/constants.hpp @@ -26,6 +26,7 @@ namespace metpy_constants { extern double Ls; extern double Cp_i; extern double zero_degc; + extern double epsilon; // General Meteorological constants // extern double epsilon = Mw / Md; // ≈ 0.622 diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 0b0dc786f8..e36796cc50 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -1729,7 +1729,7 @@ def dewpoint(vapor_pressure): Renamed ``e`` parameter to ``vapor_pressure`` """ - + # Calling c++ calculation module return _calc_mod.dewpoint(vapor_pressure) @exporter.export @@ -1790,7 +1790,8 @@ def mixing_ratio(partial_press, total_press, molecular_weight_ratio=mpconsts.nou Renamed ``part_press``, ``tot_press`` parameters to ``partial_press``, ``total_press`` """ - return molecular_weight_ratio * partial_press / (total_press - partial_press) + # Calling c++ calculation module + return _calc_mod.mixing_ratio(partial_press, total_press, molecular_weight_ratio) @exporter.export diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index 29a94a8a3e..6255455a45 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -4,6 +4,8 @@ #include "constants.hpp" #include "virtual_temperature.hpp" #include +#include // for std::cerr +#include // for std::numeric_limits namespace py = pybind11; namespace mc = metpy_constants; @@ -55,6 +57,29 @@ double DewPoint(double vapor_pressure) { return mc::zero_degc + 243.5 * val / (17.67 - val); } +double MixingRatio(double partial_press, double total_press, double epsilon) { + return epsilon * partial_press / (total_press - partial_press); +} + +double SaturationMixingRatio(double total_press, double temperature, std::string phase) { + double e_s = SaturationVaporPressure(temperature, phase); + if (e_s >= total_press) { + std::cerr << "Total pressure must be greater than the saturation vapor pressure " + << "for liquid water to be in equilibrium.\n"; + return std::numeric_limits::quiet_NaN(); + } + return MixingRatio(e_s, total_press); +} + +double SpecificHumidityFromMixingRatio(double mixing_ratio) { + return mixing_ratio / (mixing_ratio + 1.0); +} + +double SpecificHumidityFromDewPoint(double pressure, double dew_point, std::string phase) { + double mixing_ratio = SaturationMixingRatio(pressure, dew_point, phase); + return SpecificHumidityFromMixingRatio(mixing_ratio); +} + double VirtualTemperature(double temperature, double mixing_ratio, double epsilon) { return temperature * (mixing_ratio + epsilon) / (epsilon * (1. + mixing_ratio)); } diff --git a/src/virtual_temperature.hpp b/src/virtual_temperature.hpp index d8182c2d0c..32c6de0915 100644 --- a/src/virtual_temperature.hpp +++ b/src/virtual_temperature.hpp @@ -3,6 +3,9 @@ //#include #include +#include "constants.hpp" + +namespace mc = metpy_constants; double WaterLatentHeatVaporization(double temperature); double WaterLatentHeatSublimation(double temperature); @@ -12,6 +15,10 @@ double _SaturationVaporPressureSolid(double temperature); double SaturationVaporPressure(double temperature, std::string phase); double DewPoint(double vapor_pressure); +double MixingRatio(double partial_press, double total_press, double epsilon=mc::epsilon); +double SaturationMixingRatio(double total_press, double temperature, std::string phase); +double SpecificHumidityFromMixingRatio(double mixing_ratio); +double SpecificHumidityFromDewPoint(double pressure, double dew_point, std::string phase); double VirtualTemperature(double temperature, double mixing_ratio, double epsilon); From ad18ce35fdc67a08b46c065a338a07b52f759bdf Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Wed, 11 Jun 2025 17:50:46 -0600 Subject: [PATCH 17/41] add c++ functions: moist_air_gas_constant, moist_air_specific_heat_pressure; WIP: LCL --- src/calcmod.cpp | 12 ++++++++++++ src/constants.cpp | 4 ++++ src/constants.hpp | 2 ++ src/virtual_temperature.cpp | 28 +++++++++++++++++++++++++--- src/virtual_temperature.hpp | 5 +++++ 5 files changed, 48 insertions(+), 3 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 7fe49686b7..16b4b13f7c 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -16,6 +16,14 @@ PYBIND11_MODULE(_calc_mod, m) { metpy_constants::load_constants_from_python(); m.def("add", &add, "Add two numbers"); + + m.def("moist_air_gas_constant", py::vectorize(MoistAirGasConstant), + "Calculate R_m, the gas constant for moist air.", + py::arg("specific_humidity")); + + m.def("moist_air_specific_heat_pressure", py::vectorize(MoistAirSpecificHeatPressure), + "Calculate C_pm, the specific heat of moist air at constant pressure.", + py::arg("specific_humidity")); m.def("water_latent_heat_vaporization", py::vectorize(WaterLatentHeatVaporization), "Calculate water latent heat vaporization from temperature.", @@ -24,6 +32,10 @@ PYBIND11_MODULE(_calc_mod, m) { m.def("water_latent_heat_sublimation", py::vectorize(WaterLatentHeatSublimation), "Calculate water latent heat sublimation from temperature.", py::arg("temperature")); + + m.def("lcl", py::vectorize(LCL), + "Calculate the lifting condensation level (LCL) from pressure, temperature and dewpoint.", + py::arg("pressure"), py::arg("temperature"), py::arg("dewpoint")); m.def("saturation_vapor_pressure", py::vectorize(SaturationVaporPressure), "Calculate saturation vapor pressure from temperature.", diff --git a/src/constants.cpp b/src/constants.cpp index a83a7ed737..92c6327bea 100644 --- a/src/constants.cpp +++ b/src/constants.cpp @@ -6,7 +6,9 @@ namespace py = pybind11; namespace metpy_constants { double Mw; + double Rd; double Rv; + double Cp_d; double Cp_v; double Cp_l; double Lv; @@ -30,7 +32,9 @@ namespace metpy_constants { Mw = mod.attr("Mw").cast(); + Rd = mod.attr("Rd").cast(); Rv = mod.attr("Rv").cast(); + Cp_d = mod.attr("Cp_d").cast(); Cp_v = mod.attr("Cp_v").cast(); Cp_l = mod.attr("Cp_l").cast(); Lv = mod.attr("Lv").cast(); diff --git a/src/constants.hpp b/src/constants.hpp index 4707c7a046..a1c182623a 100644 --- a/src/constants.hpp +++ b/src/constants.hpp @@ -15,8 +15,10 @@ namespace metpy_constants { // Water extern double Mw; + extern double Rd; extern double Rv; // extern double wv_spec_heat_ratio = 1.33; + extern double Cp_d; extern double Cp_v; // extern double Cv_v = Cp_v / wv_spec_heat_ratio; // (J / kg / K) extern double Cp_l; diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index 6255455a45..63e188e2af 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -10,6 +10,14 @@ namespace py = pybind11; namespace mc = metpy_constants; +double MoistAirGasConstant(double specific_humidity) { + return (1.0 - specific_humidity) * mc::Rd + specific_humidity * mc::Rv; +} + +double MoistAirSpecificHeatPressure(double specific_humidity) { + return (1.0 - specific_humidity) * mc::Cp_d + specific_humidity * mc::Cp_v; +} + double WaterLatentHeatVaporization(double temperature) { return mc::Lv - (mc::Cp_l - mc::Cp_v) * (temperature - mc::T0); } @@ -18,6 +26,21 @@ double WaterLatentHeatSublimation(double temperature) { return mc::Ls - (mc::Cp_i - mc::Cp_v) * (temperature - mc::T0); } +double LCL(double pressure, double temperature, double dewpoint) { + if (temperature <= dewpoint) { + std::cerr << "Temperature must be greater than dew point for LCL calculation.\n"; + return std::numeric_limits::quiet_NaN(); + } + + double q = SpecificHumidityFromDewPoint(pressure, dewpoint, "liquid"); + double moist_heat_ratio = MoistAirSpecificHeatPressure(q) / MoistAirGasConstant(q); + double spec_heat_diff = mc::Cp_l - mc::Cp_v; + + double a = moist_heat_ratio + spec_heat_diff / mc::Rv; + + return ; +} + double _SaturationVaporPressureLiquid(double temperature) { double latent_heat = WaterLatentHeatVaporization(temperature); double heat_power = (mc::Cp_l - mc::Cp_v) / mc::Rv; @@ -52,7 +75,6 @@ double SaturationVaporPressure(double temperature, std::string phase) { } double DewPoint(double vapor_pressure) { - double val = log(vapor_pressure / mc::sat_pressure_0c); return mc::zero_degc + 243.5 * val / (17.67 - val); } @@ -75,8 +97,8 @@ double SpecificHumidityFromMixingRatio(double mixing_ratio) { return mixing_ratio / (mixing_ratio + 1.0); } -double SpecificHumidityFromDewPoint(double pressure, double dew_point, std::string phase) { - double mixing_ratio = SaturationMixingRatio(pressure, dew_point, phase); +double SpecificHumidityFromDewPoint(double pressure, double dewpoint, std::string phase) { + double mixing_ratio = SaturationMixingRatio(pressure, dewpoint, phase); return SpecificHumidityFromMixingRatio(mixing_ratio); } diff --git a/src/virtual_temperature.hpp b/src/virtual_temperature.hpp index 32c6de0915..3352012e1c 100644 --- a/src/virtual_temperature.hpp +++ b/src/virtual_temperature.hpp @@ -7,9 +7,14 @@ namespace mc = metpy_constants; +double MoistAirGasConstant(double specific_humidity); +double MoistAirSpecificHeatPressure(double specific_humidity); + double WaterLatentHeatVaporization(double temperature); double WaterLatentHeatSublimation(double temperature); +double LCL(double pressure, double temperature, double dewpoint); + double _SaturationVaporPressureLiquid(double temperature); double _SaturationVaporPressureSolid(double temperature); double SaturationVaporPressure(double temperature, std::string phase); From 89a17abd57124b9cee2e5d680840c12ca9c99604 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Wed, 11 Jun 2025 17:52:35 -0600 Subject: [PATCH 18/41] wip --- src/calcmod.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 16b4b13f7c..87c79c954e 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -50,7 +50,7 @@ PYBIND11_MODULE(_calc_mod, m) { py::arg("temperature")); m.def("dewpoint", py::vectorize(DewPoint), - "Calculate dew point from water vapor partial pressure.", + "Calculate dewpoint from water vapor partial pressure.", py::arg("vapor_pressure")); m.def("mixing_ratio", py::vectorize(MixingRatio), From 0af355b0657ca587008a335116522da8c8a3c028 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Thu, 12 Jun 2025 11:15:30 -0600 Subject: [PATCH 19/41] add function lambert_wm1: Lambert W function -1 section; LCL c++ WIP --- src/CMakeLists.txt | 1 + src/calcmod.cpp | 6 +++++- src/constants_backup | 37 ------------------------------------- src/math.cpp | 27 +++++++++++++++++++++++++++ src/math.hpp | 7 +++++++ src/virtual_temperature.cpp | 16 +++++++++++++++- src/virtual_temperature.hpp | 2 ++ 7 files changed, 57 insertions(+), 39 deletions(-) delete mode 100644 src/constants_backup create mode 100644 src/math.cpp create mode 100644 src/math.hpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index c6d61218a2..aa59e45f65 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -14,6 +14,7 @@ find_package(pybind11 CONFIG REQUIRED) pybind11_add_module(_calc_mod MODULE calcmod.cpp constants.cpp + math.cpp virtual_temperature.cpp ) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 87c79c954e..3de2918ca9 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -32,7 +32,11 @@ PYBIND11_MODULE(_calc_mod, m) { m.def("water_latent_heat_sublimation", py::vectorize(WaterLatentHeatSublimation), "Calculate water latent heat sublimation from temperature.", py::arg("temperature")); - + + m.def("relative_humidity_from_dewpoint", py::vectorize(RelativeHumidityFromDewPoint), + "Calculate relative humidity from temperature and dewpoint.", + py::arg("temperature"), py::arg("dewpoint"), py::arg("phase")); + m.def("lcl", py::vectorize(LCL), "Calculate the lifting condensation level (LCL) from pressure, temperature and dewpoint.", py::arg("pressure"), py::arg("temperature"), py::arg("dewpoint")); diff --git a/src/constants_backup b/src/constants_backup deleted file mode 100644 index 8e9d596a61..0000000000 --- a/src/constants_backup +++ /dev/null @@ -1,37 +0,0 @@ -#ifndef CONSTANTS_HPP -#define CONSTANTS_HPP - -namespace metpy_constants { -// Gas constants (J / kg / K) -constexpr double R = 8.314462618; // Universal gas constant (J / mol / K) - -// Dry air -constexpr double Md = 28.96546e-3; // Molar mass of dry air (kg / mol) -constexpr double Rd = R / Md; // Dry air (J / kg / K) -constexpr double dry_air_spec_heat_ratio = 1.4; -constexpr double Cp_d = Rd * dry_air_spec_heat_ratio / (dry_air_spec_heat_ratio - 1.0); // (J / kg / K) -constexpr double Cv_d = Cp_d / dry_air_spec_heat_ratio; // (J / kg / K) - -// Water -constexpr double Mw = 18.015268e-3; // Molar mass of water (kg / mol) -constexpr double Rv = R / Mw; // Water vapor (J / kg / K) -constexpr double wv_spec_heat_ratio = 1.33; -constexpr double Cp_v = Rv * wv_spec_heat_ratio / (wv_spec_heat_ratio - 1.0); // (J / kg / K) -constexpr double Cv_v = Cp_v / wv_spec_heat_ratio; // (J / kg / K) -constexpr double Cp_l = 4.2194e3; // Specific heat capacity of liquid water (J / kg / K) -constexpr double Lv = 2.50084e6; // Latent heat of vaporization of water (J / kg) -constexpr double T0 = 273.16; // Triple point of water (K) - -// General Meteorological constants -constexpr double epsilon = Mw / Md; // ≈ 0.622 - - -// Standard gravity -constexpr double g = 9.80665; // m / s^2 - -// Reference pressure -constexpr double P0 = 100000.0; // Pa (hPa = 1000) - -} - -#endif // CONSTANTS_HPP diff --git a/src/math.cpp b/src/math.cpp new file mode 100644 index 0000000000..816cca8f47 --- /dev/null +++ b/src/math.cpp @@ -0,0 +1,27 @@ +#include +#include +#include "math.hpp" + +double lambert_wm1(double x, double tol, int max_iter) { + if (x >= 0 || x < -1.0 / std::exp(1.0)) { + throw std::domain_error("W₋₁(x) is only defined for -1/e < x < 0"); + } + + // Initial guess for W₋₁, from Corless et al. (1996) // need to check if this is + // correct + double L1 = std::log(-x); + double L2 = std::log(-L1); + double w = L1 - L2 + (L2 / L1); + + for (int i = 0; i < max_iter; ++i) { + double ew = std::exp(w); + double wew = w * ew; + double diff = (wew - x) / (ew * (w + 1) - ((w + 2) * (wew - x)) / (2 * w + 2)); + w -= diff; + if (std::abs(diff) < tol) { + return w; + } + } + + throw std::runtime_error("lambert_wm1 did not converge"); +} diff --git a/src/math.hpp b/src/math.hpp new file mode 100644 index 0000000000..c4c99ea255 --- /dev/null +++ b/src/math.hpp @@ -0,0 +1,7 @@ +#ifndef MATH_HPP +#define MATH_HPP + +double lambert_wm1(double x, double tol = 1e-10, int max_iter = 100); + +#endif // MATH_HPP + diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index 63e188e2af..5fbbab35c9 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -1,6 +1,7 @@ #include #include #include +#include "math.hpp" #include "constants.hpp" #include "virtual_temperature.hpp" #include @@ -26,6 +27,12 @@ double WaterLatentHeatSublimation(double temperature) { return mc::Ls - (mc::Cp_i - mc::Cp_v) * (temperature - mc::T0); } +double RelativeHumidityFromDewPoint(double temperature, double dewpoint, std::string phase) { + double e_s = SaturationVaporPressure(temperature, phase); + double e = SaturationVaporPressure(dewpoint, phase); + return e / e_s; +} + double LCL(double pressure, double temperature, double dewpoint) { if (temperature <= dewpoint) { std::cerr << "Temperature must be greater than dew point for LCL calculation.\n"; @@ -37,8 +44,15 @@ double LCL(double pressure, double temperature, double dewpoint) { double spec_heat_diff = mc::Cp_l - mc::Cp_v; double a = moist_heat_ratio + spec_heat_diff / mc::Rv; + double b = -(mc::Lv + spec_heat_diff * mc::T0) / (mc::Rv * temperature); + double c = b / a; + + double rh = RelativeHumidityFromDewPoint(temperature, dewpoint, "liquid"); + double w_minus1 = lambert_wm1(pow(rh, 1.0 / a) * c * exp(c)); + double t_lcl = c / w_minus1 * temperature; + double p_lcl = pressure * pow(t_lcl / temperature, moist_heat_ratio); - return ; + return t_lcl; // returning t_lcl and p_lcl together is needed } double _SaturationVaporPressureLiquid(double temperature) { diff --git a/src/virtual_temperature.hpp b/src/virtual_temperature.hpp index 3352012e1c..c360006c5b 100644 --- a/src/virtual_temperature.hpp +++ b/src/virtual_temperature.hpp @@ -13,6 +13,8 @@ double MoistAirSpecificHeatPressure(double specific_humidity); double WaterLatentHeatVaporization(double temperature); double WaterLatentHeatSublimation(double temperature); +double RelativeHumidityFromDewPoint(double temperature, double dewpoint, std::string phase="liquid"); + double LCL(double pressure, double temperature, double dewpoint); double _SaturationVaporPressureLiquid(double temperature); From 94803bd748c754e9ade0af275afb7fde8b4c9284 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Sun, 15 Jun 2025 15:38:37 -0600 Subject: [PATCH 20/41] achieved manually vectorization of function LCL, support np.array of multiple dimensions --- src/calcmod.cpp | 64 +++++++++++++++++++++++++++++++++++-- src/virtual_temperature.cpp | 7 ++-- src/virtual_temperature.hpp | 4 +-- 3 files changed, 67 insertions(+), 8 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 3de2918ca9..a91f8b2c3d 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -3,6 +3,8 @@ #include #include #include +#include // For std::pair +#include // For std::make_tuple namespace py = pybind11; @@ -37,9 +39,65 @@ PYBIND11_MODULE(_calc_mod, m) { "Calculate relative humidity from temperature and dewpoint.", py::arg("temperature"), py::arg("dewpoint"), py::arg("phase")); - m.def("lcl", py::vectorize(LCL), - "Calculate the lifting condensation level (LCL) from pressure, temperature and dewpoint.", - py::arg("pressure"), py::arg("temperature"), py::arg("dewpoint")); +// m.def("lcl", py::vectorize(LCL), +// "Calculate the lifting condensation level (LCL) from pressure, temperature and dewpoint.", +// py::arg("pressure"), py::arg("temperature"), py::arg("dewpoint")); + + m.def("lcl", [](py::array_t pressure, + py::array_t temperature, + py::array_t dewpoint) { + + // This helper ensures the arrays are in C-style contiguous memory. + // If an input array is already contiguous, it's a zero-cost operation. + // If it's a slice or has a different memory layout, it creates a copy. + // This makes the subsequent looping simple and safe. + auto p_contig = py::array::ensure(pressure, py::array::c_style); + auto t_contig = py::array::ensure(temperature, py::array::c_style); + auto d_contig = py::array::ensure(dewpoint, py::array::c_style); + + // --- Step 1: Check that all input arrays have the same shape --- + if (p_contig.ndim() != t_contig.ndim() || p_contig.ndim() != d_contig.ndim()) { + throw std::runtime_error("Input arrays must have the same number of dimensions."); + } + for (int i = 0; i < p_contig.ndim(); ++i) { + if (p_contig.shape(i) != t_contig.shape(i) || p_contig.shape(i) != d_contig.shape(i)) { + throw std::runtime_error("Input arrays must have the same shape."); + } + } + + // --- Step 2: Create output arrays with the exact same N-D shape as the inputs --- + auto p_lcl = py::array_t(p_contig.request().shape); + auto t_lcl = py::array_t(p_contig.request().shape); + + // --- Step 3: Get the total number of elements to loop over --- + ssize_t size = p_contig.size(); + + // --- Step 4: Get direct pointers to the (now contiguous) data buffers --- + const double* p_ptr = static_cast(p_contig.request().ptr); + const double* t_ptr = static_cast(t_contig.request().ptr); + const double* d_ptr = static_cast(d_contig.request().ptr); + double* p_lcl_ptr = p_lcl.mutable_data(); + double* t_lcl_ptr = t_lcl.mutable_data(); + + // --- Step 5: Loop through the data as if it were a single flat 1D array --- + for (ssize_t i = 0; i < size; i++) { + // Call the scalar c++ function for each element + std::pair result = LCL(p_ptr[i], t_ptr[i], d_ptr[i]); + + p_lcl_ptr[i] = result.first; + t_lcl_ptr[i] = result.second; + } + + // --- Step 6: Return a tuple of the two new, N-dimensional arrays --- + return std::make_tuple(p_lcl, t_lcl); + + }, "Calculate the lifting condensation level (LCL) from pressure, temperature and dewpoint.", + py::arg("pressure"), py::arg("temperature"), py::arg("dewpoint")); + + + + + m.def("saturation_vapor_pressure", py::vectorize(SaturationVaporPressure), "Calculate saturation vapor pressure from temperature.", diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index 5fbbab35c9..d523cbbb79 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -1,5 +1,6 @@ #include #include +#include // For std::pair #include #include "math.hpp" #include "constants.hpp" @@ -33,10 +34,10 @@ double RelativeHumidityFromDewPoint(double temperature, double dewpoint, std::st return e / e_s; } -double LCL(double pressure, double temperature, double dewpoint) { +std::pair LCL(double pressure, double temperature, double dewpoint) { if (temperature <= dewpoint) { std::cerr << "Temperature must be greater than dew point for LCL calculation.\n"; - return std::numeric_limits::quiet_NaN(); + return {std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()}; } double q = SpecificHumidityFromDewPoint(pressure, dewpoint, "liquid"); @@ -52,7 +53,7 @@ double LCL(double pressure, double temperature, double dewpoint) { double t_lcl = c / w_minus1 * temperature; double p_lcl = pressure * pow(t_lcl / temperature, moist_heat_ratio); - return t_lcl; // returning t_lcl and p_lcl together is needed + return {p_lcl, t_lcl}; // returning t_lcl and p_lcl together is needed } double _SaturationVaporPressureLiquid(double temperature) { diff --git a/src/virtual_temperature.hpp b/src/virtual_temperature.hpp index c360006c5b..51488eeabb 100644 --- a/src/virtual_temperature.hpp +++ b/src/virtual_temperature.hpp @@ -1,8 +1,8 @@ #ifndef VIRTUAL_TEMPERATURE_HPP // if not defined #define VIRTUAL_TEMPERATURE_HPP // define the header file -//#include #include +#include // For std::pair #include "constants.hpp" namespace mc = metpy_constants; @@ -15,7 +15,7 @@ double WaterLatentHeatSublimation(double temperature); double RelativeHumidityFromDewPoint(double temperature, double dewpoint, std::string phase="liquid"); -double LCL(double pressure, double temperature, double dewpoint); +std::pair LCL(double pressure, double temperature, double dewpoint); double _SaturationVaporPressureLiquid(double temperature); double _SaturationVaporPressureSolid(double temperature); From 8f50f44a5605aa1d38bbd06533b4a07e12568a97 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Mon, 16 Jun 2025 09:41:54 -0600 Subject: [PATCH 21/41] add reference to Lambert W function calculation --- src/math.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/math.cpp b/src/math.cpp index 816cca8f47..5b89a5f4cc 100644 --- a/src/math.cpp +++ b/src/math.cpp @@ -3,12 +3,13 @@ #include "math.hpp" double lambert_wm1(double x, double tol, int max_iter) { + // Corless, R.M., Gonnet, G.H., Hare, D.E.G. et al. On the LambertW function. + // Adv Comput Math 5, 329–359 (1996). https://doi.org/10.1007/BF02124750 + if (x >= 0 || x < -1.0 / std::exp(1.0)) { throw std::domain_error("W₋₁(x) is only defined for -1/e < x < 0"); } - // Initial guess for W₋₁, from Corless et al. (1996) // need to check if this is - // correct double L1 = std::log(-x); double L2 = std::log(-L1); double w = L1 - L2 + (L2 / L1); From 76ce4c2952c3aba6b2f971f53e6f14de61a5e011 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Mon, 16 Jun 2025 16:52:51 -0600 Subject: [PATCH 22/41] set warning messages and return NAN value --- src/math.cpp | 10 +++++++--- src/virtual_temperature.cpp | 8 +++++--- 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/src/math.cpp b/src/math.cpp index 5b89a5f4cc..2c0dc4f39d 100644 --- a/src/math.cpp +++ b/src/math.cpp @@ -1,5 +1,5 @@ #include -#include +#include // for std::cerr #include "math.hpp" double lambert_wm1(double x, double tol, int max_iter) { @@ -7,7 +7,9 @@ double lambert_wm1(double x, double tol, int max_iter) { // Adv Comput Math 5, 329–359 (1996). https://doi.org/10.1007/BF02124750 if (x >= 0 || x < -1.0 / std::exp(1.0)) { - throw std::domain_error("W₋₁(x) is only defined for -1/e < x < 0"); + std::cerr << "Warning in function '" << __func__ + << "': lambert_wm1 is only defined for -1/e < x < 0\n"; + return std::numeric_limits::quiet_NaN(); } double L1 = std::log(-x); @@ -24,5 +26,7 @@ double lambert_wm1(double x, double tol, int max_iter) { } } - throw std::runtime_error("lambert_wm1 did not converge"); + std::cerr << "Warning in function '" << __func__ + << "': lambert_wm1 did not converge.\n"; + return std::numeric_limits::quiet_NaN(); } diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index d523cbbb79..5222b14bf0 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -36,7 +36,8 @@ double RelativeHumidityFromDewPoint(double temperature, double dewpoint, std::st std::pair LCL(double pressure, double temperature, double dewpoint) { if (temperature <= dewpoint) { - std::cerr << "Temperature must be greater than dew point for LCL calculation.\n"; + std::cerr << "Warning in function '" << __func__ + << "': Temperature must be greater than dew point for LCL calculation.\n"; return {std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()}; } @@ -101,8 +102,9 @@ double MixingRatio(double partial_press, double total_press, double epsilon) { double SaturationMixingRatio(double total_press, double temperature, std::string phase) { double e_s = SaturationVaporPressure(temperature, phase); if (e_s >= total_press) { - std::cerr << "Total pressure must be greater than the saturation vapor pressure " - << "for liquid water to be in equilibrium.\n"; + std::cerr << "Warning in function '" << __func__ + << "': Total pressure must be greater than the saturation vapor pressure " + << "for liquid water to be in equilibrium.\n"; return std::numeric_limits::quiet_NaN(); } return MixingRatio(e_s, total_press); From 097844add4b2fd36f7f5b00bb3051589a8a8636f Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Wed, 18 Jun 2025 17:21:15 -0600 Subject: [PATCH 23/41] add c++ dry_lapse --- src/calcmod.cpp | 6 +++--- src/constants.cpp | 2 ++ src/constants.hpp | 1 + src/metpy/calc/thermo.py | 24 ++++++++++++++++++++++++ src/virtual_temperature.cpp | 35 ++++++++++++++++++++++++++++++++--- src/virtual_temperature.hpp | 4 ++++ 6 files changed, 66 insertions(+), 6 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index a91f8b2c3d..24e007f239 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -39,9 +39,9 @@ PYBIND11_MODULE(_calc_mod, m) { "Calculate relative humidity from temperature and dewpoint.", py::arg("temperature"), py::arg("dewpoint"), py::arg("phase")); -// m.def("lcl", py::vectorize(LCL), -// "Calculate the lifting condensation level (LCL) from pressure, temperature and dewpoint.", -// py::arg("pressure"), py::arg("temperature"), py::arg("dewpoint")); + m.def("dry_lapse", py::vectorize(DryLapse), + "Calculate the temperature at a level assuming only dry adiabatic process.", + py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure")); m.def("lcl", [](py::array_t pressure, py::array_t temperature, diff --git a/src/constants.cpp b/src/constants.cpp index 92c6327bea..6cc7bed626 100644 --- a/src/constants.cpp +++ b/src/constants.cpp @@ -18,6 +18,7 @@ namespace metpy_constants { double Cp_i; double zero_degc; double epsilon; + double kappa; void load_constants_from_python() { py::object mod = py::module_::import("metpy.constants.nounit"); @@ -44,5 +45,6 @@ namespace metpy_constants { Cp_i = mod.attr("Cp_i").cast(); zero_degc = mod.attr("zero_degc").cast(); epsilon = mod.attr("epsilon").cast(); + kappa = mod.attr("kappa").cast(); } } diff --git a/src/constants.hpp b/src/constants.hpp index a1c182623a..02271e28f6 100644 --- a/src/constants.hpp +++ b/src/constants.hpp @@ -29,6 +29,7 @@ namespace metpy_constants { extern double Cp_i; extern double zero_degc; extern double epsilon; + extern double kappa; // General Meteorological constants // extern double epsilon = Mw / Md; // ≈ 0.622 diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index e36796cc50..a1bc96913e 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -513,6 +513,30 @@ def dry_lapse(pressure, temperature, reference_pressure=None, vertical_dim=0): return temperature * (pressure / reference_pressure)**mpconsts.kappa + +@exporter.export +@preprocess_and_wrap( + wrap_like='temperature', + broadcast=('pressure', 'temperature', 'reference_pressure') +) +@process_units( + { + 'pressure': '[pressure]', + 'temperature': '[temperature]', + 'reference_pressure': '[pressure]' + }, + '[temperature]' +) +def dry_lapse_linfel(pressure, temperature, reference_pressure=None, vertical_dim=0): + """ + Linfeng's version of 'dry_lapse'. Added on Jun18 2025 + """ + if reference_pressure is None: + reference_pressure = pressure[0] + return _calc_mod.dry_lapse(pressure, temperature, reference_pressure) + + + @exporter.export @preprocess_and_wrap( wrap_like='temperature', diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index 5222b14bf0..12d835019f 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -1,13 +1,14 @@ #include #include +#include #include // For std::pair #include -#include "math.hpp" -#include "constants.hpp" -#include "virtual_temperature.hpp" #include #include // for std::cerr #include // for std::numeric_limits +#include "math.hpp" +#include "constants.hpp" +#include "virtual_temperature.hpp" namespace py = pybind11; namespace mc = metpy_constants; @@ -34,6 +35,12 @@ double RelativeHumidityFromDewPoint(double temperature, double dewpoint, std::st return e / e_s; } +double DryLapse(double pressure, double ref_temperature, double ref_pressure) { + // calculate temperature at pressure given reference temperature and pressure + // assuming dry adiabatic process + return ref_temperature * pow(pressure / ref_pressure, mc::kappa); +} + std::pair LCL(double pressure, double temperature, double dewpoint) { if (temperature <= dewpoint) { std::cerr << "Warning in function '" << __func__ @@ -57,6 +64,28 @@ std::pair LCL(double pressure, double temperature, double dewpoi return {p_lcl, t_lcl}; // returning t_lcl and p_lcl together is needed } +bool _CheckPressure(const std::vector& pressure) { + for (size_t i = 0; i + 1 < pressure.size(); ++i) { + if (pressure[i] < pressure[i + 1]) { + return false; // pressure increases (invalid) + } + } + return true; // strictly non-increasing +} + +void _ParcelProfileHelper(const std::vector& pressure, double temperature, double dewpoint) { + // Check that pressure does not increase. + if (!_CheckPressure(pressure)) { + throw std::runtime_error( + "Pressure increases between at least two points in your sounding. " + "Using a smoothing filter (e.g., scipy.signal.medfilt) may fix this."); + } + + // Find the LCL + //press_lcl, temp_lcl = LCL(pressure[0], temperature, dewpoint); + +} + double _SaturationVaporPressureLiquid(double temperature) { double latent_heat = WaterLatentHeatVaporization(temperature); double heat_power = (mc::Cp_l - mc::Cp_v) / mc::Rv; diff --git a/src/virtual_temperature.hpp b/src/virtual_temperature.hpp index 51488eeabb..8620dee62a 100644 --- a/src/virtual_temperature.hpp +++ b/src/virtual_temperature.hpp @@ -2,6 +2,7 @@ #define VIRTUAL_TEMPERATURE_HPP // define the header file #include +#include #include // For std::pair #include "constants.hpp" @@ -14,9 +15,12 @@ double WaterLatentHeatVaporization(double temperature); double WaterLatentHeatSublimation(double temperature); double RelativeHumidityFromDewPoint(double temperature, double dewpoint, std::string phase="liquid"); +double DryLapse(double pressure, double ref_temperature, double ref_pressure); std::pair LCL(double pressure, double temperature, double dewpoint); +bool _CheckPressure(const std::vector& pressure); + double _SaturationVaporPressureLiquid(double temperature); double _SaturationVaporPressureSolid(double temperature); double SaturationVaporPressure(double temperature, std::string phase); From 94fde62b07404f355bc487ba4b92293dae77996e Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Wed, 18 Jun 2025 18:10:11 -0600 Subject: [PATCH 24/41] add DryLapseProfile for c++ internal vectorization of DryLapse; _ParcelProfileHelper wip --- src/virtual_temperature.cpp | 25 +++++++++++++++++++++++-- src/virtual_temperature.hpp | 4 ++++ 2 files changed, 27 insertions(+), 2 deletions(-) diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index 12d835019f..0d846413bb 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -41,6 +41,22 @@ double DryLapse(double pressure, double ref_temperature, double ref_pressure) { return ref_temperature * pow(pressure / ref_pressure, mc::kappa); } +std::vector DryLapseProfile(const std::vector& pressure_profile, + double ref_temperature, + double ref_pressure) { + // Vectorized version of DryLapse for a pressure profile. C++ internally use. + + // calculate temperature profile of an air parcel lifting dry adiabatically + // through the given pressure profile. + std::vector temperature_profile; + temperature_profile.reserve(pressure_profile.size()); + + for (double p : pressure_profile) { + temperature_profile.push_back(DryLapse(p, ref_temperature, ref_pressure)); + } + return temperature_profile; +} + std::pair LCL(double pressure, double temperature, double dewpoint) { if (temperature <= dewpoint) { std::cerr << "Warning in function '" << __func__ @@ -61,7 +77,7 @@ std::pair LCL(double pressure, double temperature, double dewpoi double t_lcl = c / w_minus1 * temperature; double p_lcl = pressure * pow(t_lcl / temperature, moist_heat_ratio); - return {p_lcl, t_lcl}; // returning t_lcl and p_lcl together is needed + return {p_lcl, t_lcl}; } bool _CheckPressure(const std::vector& pressure) { @@ -82,7 +98,12 @@ void _ParcelProfileHelper(const std::vector& pressure, double temperatur } // Find the LCL - //press_lcl, temp_lcl = LCL(pressure[0], temperature, dewpoint); + std::pair result = LCL(pressure[0], temperature, dewpoint); + double press_lcl = result.first; + double temp_lcl = result.second; + + std::vector press_lower; + std::vector temp_lower = DryLapseProfile(press_lower, temperature, press_lower[0]); } diff --git a/src/virtual_temperature.hpp b/src/virtual_temperature.hpp index 8620dee62a..2c28fcf283 100644 --- a/src/virtual_temperature.hpp +++ b/src/virtual_temperature.hpp @@ -15,7 +15,11 @@ double WaterLatentHeatVaporization(double temperature); double WaterLatentHeatSublimation(double temperature); double RelativeHumidityFromDewPoint(double temperature, double dewpoint, std::string phase="liquid"); + double DryLapse(double pressure, double ref_temperature, double ref_pressure); +std::vector DryLapseProfile(const std::vector& pressure_profile, + double ref_temperature, + double ref_pressure); std::pair LCL(double pressure, double temperature, double dewpoint); From e9feef389308af479cd8f2fd4d05c0a360df9ccb Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Fri, 20 Jun 2025 22:28:59 -0600 Subject: [PATCH 25/41] add moist lapse, use rk4 to integrate dlnT_dlnP --- src/calcmod.cpp | 6 ++- src/virtual_temperature.cpp | 83 +++++++++++++++++++++++++++++++++++-- src/virtual_temperature.hpp | 7 ++++ 3 files changed, 92 insertions(+), 4 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 24e007f239..c430380f39 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -40,9 +40,13 @@ PYBIND11_MODULE(_calc_mod, m) { py::arg("temperature"), py::arg("dewpoint"), py::arg("phase")); m.def("dry_lapse", py::vectorize(DryLapse), - "Calculate the temperature at a level assuming only dry adiabatic process.", + "Calculate the temperature at pressure levels assuming dry adiabatic process.", py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure")); + m.def("moist_lapse", py::vectorize(MoistLapse), + "Calculate the temperature at pressure levels assuming saturated adiabatic process.", + py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure"), py::arg("rk_nstep")); + m.def("lcl", [](py::array_t pressure, py::array_t temperature, py::array_t dewpoint) { diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index 0d846413bb..a967c66926 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -57,6 +57,61 @@ std::vector DryLapseProfile(const std::vector& pressure_profile, return temperature_profile; } +double CaldlnTdlnP(double temperature, double pressure) { + // Calculate dlnT/dlnP for a moist (saturated) adiabatic process. + double rs = SaturationMixingRatio(pressure, temperature, "liquid"); + double dlnT_dlnP_linfel = (mc::Rd + rs * mc::Rv) / (mc::Cp_d + rs * mc::Cp_v + + (mc::Lv * mc::Lv * rs * mc::epsilon) / (mc::Rd * temperature * temperature)); + + double dlnT_dlnP_Bakhshaii2013 = (mc::Rd + mc::Lv * rs / temperature) / (mc::Cp_d + + (mc::Lv * mc::Lv * rs * mc::epsilon) / (mc::Rd * temperature * temperature)); + + return dlnT_dlnP_Bakhshaii2013; +} + +double MoistLapse(double pressure, double ref_temperature, double ref_pressure, int nstep) { + // calculate temperature at pressure given reference temperature and pressure + // assuming moist adiabatic expansion (vapor condenses and removed from the air + // parcel) + + double dlnP = log(pressure / ref_pressure) / (double)nstep; + double T1 = ref_temperature; + double P1 = ref_pressure; + double k[4]; + + for (int i = 0; i < nstep; ++i) { + k[0] = CaldlnTdlnP(T1, P1); + k[1] = CaldlnTdlnP(T1 * exp(k[0] * dlnP/2.), P1 * exp(dlnP/2.)); + k[2] = CaldlnTdlnP(T1 * exp(k[1] * dlnP/2.), P1 * exp(dlnP/2.)); + k[3] = CaldlnTdlnP(T1 * exp(k[2] * dlnP), P1 * exp(dlnP)); + + T1 = T1 * exp((k[0] + 2.0 * k[1] + 2.0 * k[2] + k[3]) * dlnP / 6.0); + P1 = P1 * exp(dlnP); + } + + return T1; // check final T1 P1 +} + +std::vector MoistLapseProfile(const std::vector& pressure_profile, + double ref_temperature, + double ref_pressure) { + // Vectorized version of MoistLapse for a pressure profile. C++ internally use. + + // calculate temperature profile of an air parcel lifting saturated adiabatically + // through the given pressure profile. + std::vector temperature_profile; + temperature_profile.reserve(pressure_profile.size()); + + //double T1 = ref_temperature; + //double P1 = ref_pressure; + double T; + for (size_t i = 0; i < pressure_profile.size(); ++i) { + T = MoistLapse(pressure_profile[i], ref_temperature, ref_pressure, 50); + temperature_profile.push_back(T); + } + return temperature_profile; +} + std::pair LCL(double pressure, double temperature, double dewpoint) { if (temperature <= dewpoint) { std::cerr << "Warning in function '" << __func__ @@ -98,13 +153,35 @@ void _ParcelProfileHelper(const std::vector& pressure, double temperatur } // Find the LCL - std::pair result = LCL(pressure[0], temperature, dewpoint); - double press_lcl = result.first; - double temp_lcl = result.second; + auto [press_lcl, temp_lcl] = LCL(pressure[0], temperature, dewpoint); + // Establish profile below LCL std::vector press_lower; + for (double p : pressure) { + if (p >= press_lcl) { + press_lower.push_back(p); + } + } + press_lower.push_back(press_lcl); std::vector temp_lower = DryLapseProfile(press_lower, temperature, press_lower[0]); + // Early return if profile ends before reaching above LCL + if (pressure.back() >= press_lcl) { + press_lower.pop_back(); + temp_lower.pop_back(); +// return {press_lower, {}, press_lcl, temp_lower, {}, temp_lcl}; + } + + // Establish profile above LCL + std::vector press_upper; + press_upper.push_back(press_lcl); + for (double p : pressure) { + if (p < press_lcl) { + press_upper.push_back(p); + } + } + + } double _SaturationVaporPressureLiquid(double temperature) { diff --git a/src/virtual_temperature.hpp b/src/virtual_temperature.hpp index 2c28fcf283..e362246d7a 100644 --- a/src/virtual_temperature.hpp +++ b/src/virtual_temperature.hpp @@ -21,6 +21,13 @@ std::vector DryLapseProfile(const std::vector& pressure_profile, double ref_temperature, double ref_pressure); +double CaldlnTdlnP(double temperature, double pressure); +double MoistLapse(double pressure, double ref_temperature, double ref_pressure, int nstep); +std::vector MoistLapseProfile(const std::vector& pressure_profile, + double ref_temperature, + double ref_pressure); + + std::pair LCL(double pressure, double temperature, double dewpoint); bool _CheckPressure(const std::vector& pressure); From 569c92f2352b2a83eea3379bddf1709e21f0dfac Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Mon, 23 Jun 2025 13:38:53 -0600 Subject: [PATCH 26/41] use Bakhshaii2013 for c++ moist_lapse --- src/metpy/calc/thermo.py | 24 ++++++++++++++++++++++++ src/virtual_temperature.cpp | 6 +++--- 2 files changed, 27 insertions(+), 3 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index a1bc96913e..3859fb1b05 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -675,6 +675,30 @@ def dt(p, t): return ret.squeeze() +@exporter.export +@preprocess_and_wrap( + wrap_like='temperature', + broadcast=('pressure', 'temperature', 'reference_pressure') +) +@process_units( + { + 'pressure': '[pressure]', + 'temperature': '[temperature]', + 'reference_pressure': '[pressure]' + }, + '[temperature]' +) +def moist_lapse_linfel(pressure, temperature, reference_pressure=None): + """ + Linfeng's version of 'moist_lapse'. Added on Jun23 2025 + """ + if reference_pressure is None: + reference_pressure = pressure[0] + # nstep for RK4 is set to 30 + return _calc_mod.moist_lapse(pressure, temperature, reference_pressure, 30) + + + @exporter.export @preprocess_and_wrap() @process_units( diff --git a/src/virtual_temperature.cpp b/src/virtual_temperature.cpp index a967c66926..cfcaea4178 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -60,9 +60,9 @@ std::vector DryLapseProfile(const std::vector& pressure_profile, double CaldlnTdlnP(double temperature, double pressure) { // Calculate dlnT/dlnP for a moist (saturated) adiabatic process. double rs = SaturationMixingRatio(pressure, temperature, "liquid"); - double dlnT_dlnP_linfel = (mc::Rd + rs * mc::Rv) / (mc::Cp_d + rs * mc::Cp_v + - (mc::Lv * mc::Lv * rs * mc::epsilon) / (mc::Rd * temperature * temperature)); - + + //double dlnT_dlnP_linfel = (mc::Rd + rs * mc::Rv) / (mc::Cp_d + rs * mc::Cp_v + + // (mc::Lv * mc::Lv * rs * mc::epsilon) / (mc::Rd * temperature * temperature)); double dlnT_dlnP_Bakhshaii2013 = (mc::Rd + mc::Lv * rs / temperature) / (mc::Cp_d + (mc::Lv * mc::Lv * rs * mc::epsilon) / (mc::Rd * temperature * temperature)); From ca222aa26e0a02a4318a416d91e845f724c2dbd7 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Mon, 23 Jun 2025 16:02:24 -0600 Subject: [PATCH 27/41] import c++ lcl into metpy thermo.py --- src/metpy/calc/thermo.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 3859fb1b05..f9fafbc431 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -784,6 +784,21 @@ def lcl(pressure, temperature, dewpoint, max_iters=None, eps=None): return p_lcl, t_lcl +@exporter.export +@preprocess_and_wrap() +@process_units( + {'pressure': '[pressure]', 'temperature': '[temperature]', 'dewpoint': '[temperature]'}, + ('[pressure]', '[temperature]') +) +def lcl_linfel(pressure, temperature, dewpoint, max_iters=None, eps=None): + """ + Linfeng's version of 'lcl'. Added on Jun23 2025 + """ + p_lcl, t_lcl = _calc_mod.lcl(pressure, temperature, dewpoint) + return p_lcl, t_lcl + + + @exporter.export @preprocess_and_wrap() @check_units('[pressure]', '[temperature]', '[temperature]') From 7986a8acf931ca929cd0e5b13da642e3b8ce75d4 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Mon, 23 Jun 2025 16:36:53 -0600 Subject: [PATCH 28/41] change file name 'virtual_temperature.cpp' -> 'thermo.cpp' --- src/CMakeLists.txt | 2 +- src/calcmod.cpp | 2 +- src/{virtual_temperature.cpp => thermo.cpp} | 2 +- src/{virtual_temperature.hpp => thermo.hpp} | 6 +++--- 4 files changed, 6 insertions(+), 6 deletions(-) rename src/{virtual_temperature.cpp => thermo.cpp} (99%) rename src/{virtual_temperature.hpp => thermo.hpp} (93%) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index aa59e45f65..d9be3d2855 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -15,7 +15,7 @@ pybind11_add_module(_calc_mod MODULE calcmod.cpp constants.cpp math.cpp - virtual_temperature.cpp + thermo.cpp ) target_compile_definitions(_calc_mod diff --git a/src/calcmod.cpp b/src/calcmod.cpp index c430380f39..5ffd42afaa 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -1,5 +1,5 @@ #include "constants.hpp" -#include "virtual_temperature.hpp" +#include "thermo.hpp" #include #include #include diff --git a/src/virtual_temperature.cpp b/src/thermo.cpp similarity index 99% rename from src/virtual_temperature.cpp rename to src/thermo.cpp index cfcaea4178..dbcb70be13 100644 --- a/src/virtual_temperature.cpp +++ b/src/thermo.cpp @@ -8,7 +8,7 @@ #include // for std::numeric_limits #include "math.hpp" #include "constants.hpp" -#include "virtual_temperature.hpp" +#include "thermo.hpp" namespace py = pybind11; namespace mc = metpy_constants; diff --git a/src/virtual_temperature.hpp b/src/thermo.hpp similarity index 93% rename from src/virtual_temperature.hpp rename to src/thermo.hpp index e362246d7a..415ce469a1 100644 --- a/src/virtual_temperature.hpp +++ b/src/thermo.hpp @@ -1,5 +1,5 @@ -#ifndef VIRTUAL_TEMPERATURE_HPP // if not defined -#define VIRTUAL_TEMPERATURE_HPP // define the header file +#ifndef THERMO_HPP // if not defined +#define THERMO_HPP // define the header file #include #include @@ -44,4 +44,4 @@ double SpecificHumidityFromDewPoint(double pressure, double dew_point, std::stri double VirtualTemperature(double temperature, double mixing_ratio, double epsilon); -#endif // VIRTUAL_TEMPERATURE_HPP +#endif // THERMO_HPP From 1f127493c54c1fc358e891cd6bb89d887d19cb9e Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Mon, 23 Jun 2025 17:13:09 -0600 Subject: [PATCH 29/41] parcel profile wip --- src/thermo.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/thermo.cpp b/src/thermo.cpp index dbcb70be13..8ad53e52f2 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -180,6 +180,9 @@ void _ParcelProfileHelper(const std::vector& pressure, double temperatur press_upper.push_back(p); } } + std::vector temp_lower = MoistLapseProfile(press_upper, temp_lcl, press_lcl); + + return //return type } From 8eb21f0f2de1f19e7d77d46700f2d5d508113d63 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Mon, 23 Jun 2025 18:08:46 -0600 Subject: [PATCH 30/41] add a return struct for _ParcelProfileHelper wip --- src/thermo.cpp | 14 +++++++++----- src/thermo.hpp | 10 ++++++++++ 2 files changed, 19 insertions(+), 5 deletions(-) diff --git a/src/thermo.cpp b/src/thermo.cpp index 8ad53e52f2..4451132fd9 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -144,7 +144,8 @@ bool _CheckPressure(const std::vector& pressure) { return true; // strictly non-increasing } -void _ParcelProfileHelper(const std::vector& pressure, double temperature, double dewpoint) { + +ParcelProfile _ParcelProfileHelper(const std::vector& pressure, double temperature, double dewpoint) { // Check that pressure does not increase. if (!_CheckPressure(pressure)) { throw std::runtime_error( @@ -169,7 +170,7 @@ void _ParcelProfileHelper(const std::vector& pressure, double temperatur if (pressure.back() >= press_lcl) { press_lower.pop_back(); temp_lower.pop_back(); -// return {press_lower, {}, press_lcl, temp_lower, {}, temp_lcl}; + return {press_lower, {}, press_lcl, temp_lower, {}, temp_lcl}; } // Establish profile above LCL @@ -180,11 +181,14 @@ void _ParcelProfileHelper(const std::vector& pressure, double temperatur press_upper.push_back(p); } } - std::vector temp_lower = MoistLapseProfile(press_upper, temp_lcl, press_lcl); - - return //return type + std::vector temp_upper = MoistLapseProfile(press_upper, temp_lower.back(), press_lcl); + press_lower.pop_back(); + temp_lower.pop_back(); + press_upper.erase(press_upper.begin()); + temp_upper.erase(temp_upper.begin()); + return {press_lower, press_upper, press_lcl, temp_lower, temp_upper, temp_lcl}; } double _SaturationVaporPressureLiquid(double temperature) { diff --git a/src/thermo.hpp b/src/thermo.hpp index 415ce469a1..f100468e52 100644 --- a/src/thermo.hpp +++ b/src/thermo.hpp @@ -32,6 +32,16 @@ std::pair LCL(double pressure, double temperature, double dewpoi bool _CheckPressure(const std::vector& pressure); +// Return struct for _ParcelProfileHelper +struct ParcelProfile { + std::vector press_lower, press_upper; + double press_lcl; + std::vector temp_lower, temp_upper; + double temp_lcl; +}; + +ParcelProfile _ParcelProfileHelper(const std::vector& pressure, double temperature, double dewpoint); + double _SaturationVaporPressureLiquid(double temperature); double _SaturationVaporPressureSolid(double temperature); double SaturationVaporPressure(double temperature, std::string phase); From e33ff90c11753587f4879597d3f1300352a3f5a3 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Tue, 24 Jun 2025 16:24:31 -0600 Subject: [PATCH 31/41] c++ ParcelProfile finished, assembled dry and moist part of the profile --- src/calcmod.cpp | 11 +++++++++-- src/metpy/calc/thermo.py | 17 +++++++++++++++++ src/thermo.cpp | 19 ++++++++++++++++++- src/thermo.hpp | 8 ++++++-- 4 files changed, 50 insertions(+), 5 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 5ffd42afaa..95398e765d 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -99,8 +99,15 @@ PYBIND11_MODULE(_calc_mod, m) { py::arg("pressure"), py::arg("temperature"), py::arg("dewpoint")); - - + m.def("parcel_profile", + [](py::array_t pressure, double temperature, double dewpoint) { + // pressure.data() gives the beginning pointer of the data buffer + std::vector pressure_vec(pressure.data(), pressure.data() + pressure.size()); + std::vector temp_prof = ParcelProfile(pressure_vec, temperature, dewpoint); + return py::array_t(temp_prof.size(), temp_prof.data()); + }, + "Compute the parcel temperature profile as it rises from a given pressure and temperature.", + py::arg("pressure"), py::arg("temperature"), py::arg("dewpoint")); m.def("saturation_vapor_pressure", py::vectorize(SaturationVaporPressure), diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index f9fafbc431..a19d40f4fe 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -1312,6 +1312,23 @@ def parcel_profile(pressure, temperature, dewpoint): return concatenate((t_l, t_u)) +@exporter.export +@preprocess_and_wrap(wrap_like='pressure') +@process_units( + { + 'pressure': '[pressure]', + 'temperature': '[temperature]', + 'dewpoint': '[temperature]' + }, + '[temperature]' +) # process units because no unit should be passed to c++ function +def parcel_profile_linfel(pressure, temperature, dewpoint): + """ + Linfeng's version of 'parcel_profile'. Added on Jun 24 2025 + """ + return _calc_mod.parcel_profile(pressure, temperature, dewpoint) + + @exporter.export @preprocess_and_wrap() @check_units('[pressure]', '[temperature]', '[temperature]') diff --git a/src/thermo.cpp b/src/thermo.cpp index 4451132fd9..f79825f27d 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -145,7 +145,24 @@ bool _CheckPressure(const std::vector& pressure) { } -ParcelProfile _ParcelProfileHelper(const std::vector& pressure, double temperature, double dewpoint) { +std::vector ParcelProfile(const std::vector& pressure, + double temperature, + double dewpoint) { + // Returns a vector of temperatures corresponding to the input pressure levels. + + ParProStruct profile = _ParcelProfileHelper(pressure, temperature, dewpoint); + + // Combine lower and upper temperature profiles + std::vector combined_temp; + combined_temp.reserve(pressure.size()); + + combined_temp.insert(combined_temp.end(), profile.temp_lower.begin(), profile.temp_lower.end()); + combined_temp.insert(combined_temp.end(), profile.temp_upper.begin(), profile.temp_upper.end()); + + return combined_temp; +} + +ParProStruct _ParcelProfileHelper(const std::vector& pressure, double temperature, double dewpoint) { // Check that pressure does not increase. if (!_CheckPressure(pressure)) { throw std::runtime_error( diff --git a/src/thermo.hpp b/src/thermo.hpp index f100468e52..f4f3c0379e 100644 --- a/src/thermo.hpp +++ b/src/thermo.hpp @@ -33,14 +33,18 @@ std::pair LCL(double pressure, double temperature, double dewpoi bool _CheckPressure(const std::vector& pressure); // Return struct for _ParcelProfileHelper -struct ParcelProfile { +struct ParProStruct { std::vector press_lower, press_upper; double press_lcl; std::vector temp_lower, temp_upper; double temp_lcl; }; -ParcelProfile _ParcelProfileHelper(const std::vector& pressure, double temperature, double dewpoint); +std::vector ParcelProfile(const std::vector& pressure, + double temperature, + double dewpoint); + +ParProStruct _ParcelProfileHelper(const std::vector& pressure, double temperature, double dewpoint); double _SaturationVaporPressureLiquid(double temperature); double _SaturationVaporPressureSolid(double temperature); From 44f1dddf364a398559229b1c1f1b1fecf585b452 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Wed, 25 Jun 2025 18:07:10 -0600 Subject: [PATCH 32/41] c++ version moist_lapse can handle 2D multiple starting temperature now --- src/calcmod.cpp | 49 +++++++++++++++++++++++++++++++++++++--- src/metpy/calc/thermo.py | 5 +++- src/thermo.cpp | 23 ++++++++----------- src/thermo.hpp | 6 ++--- 4 files changed, 63 insertions(+), 20 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 95398e765d..70aa270230 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -5,6 +5,7 @@ #include #include // For std::pair #include // For std::make_tuple +#include namespace py = pybind11; @@ -43,9 +44,51 @@ PYBIND11_MODULE(_calc_mod, m) { "Calculate the temperature at pressure levels assuming dry adiabatic process.", py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure")); - m.def("moist_lapse", py::vectorize(MoistLapse), - "Calculate the temperature at pressure levels assuming saturated adiabatic process.", - py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure"), py::arg("rk_nstep")); + m.def("moist_lapse", [](py::array_t pressure, + py::array_t ref_temperature, + double ref_pressure, + int rk_nstep) { + // This function calculates the moist adiabatic profile for multiple starting + // temperatures (2D surface) and a single communal starting pressure, along a + // 1D pressure profile. + + // --- Step 1: Prepare the C++ vector for pressure levels --- + if (pressure.ndim() != 1) { + throw std::runtime_error("Input 'pressure' array must be 1D."); + } + std::vector pressure_vec(pressure.data(), pressure.data() + pressure.size()); + + // --- Step 2: Ensure the reference temperature array is contiguous --- + auto ref_temp_contig = py::array::ensure(ref_temperature, py::array::c_style); + + // --- Step 3: Define the shape of the output array: (N+1) dimension--- + // Create a vector to hold the shape dimensions. + std::vector out_shape; + for(int i = 0; i < ref_temp_contig.ndim(); ++i) { + out_shape.push_back(ref_temp_contig.shape(i)); + } + ssize_t profile_len = pressure_vec.size(); + out_shape.push_back(profile_len); + + // Create the final output array with the correct N+1 dimensional shape. + auto out_array = py::array_t(out_shape); + + // --- Step 4: Get direct pointers to data buffers for fast access --- + const double* ref_temp_ptr = static_cast(ref_temp_contig.request().ptr); + double* out_array_ptr = out_array.mutable_data(); + ssize_t num_profiles = ref_temp_contig.size(); // Total number of profiles to calculate + + // --- Step 5: Loop through each reference temperature --- + for (ssize_t i = 0; i < num_profiles; ++i) { + for (ssize_t j = 0; j < profile_len; ++j) { + out_array_ptr[i * profile_len + j] = MoistLapse(pressure_vec[j], ref_temp_ptr[i], ref_pressure, rk_nstep); + } + } + + return out_array; + }, "Calculate the temperature along a pressure profile assuming saturated adiabatic process.", + py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure"), py::arg("rk_nstep")); + m.def("lcl", [](py::array_t pressure, py::array_t temperature, diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index a19d40f4fe..774c15240a 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -690,7 +690,10 @@ def dt(p, t): ) def moist_lapse_linfel(pressure, temperature, reference_pressure=None): """ - Linfeng's version of 'moist_lapse'. Added on Jun23 2025 + Linfeng's version of 'moist_lapse'. Added on Jun 25 2025 + This function calculates the moist adiabatic profile for multiple starting + temperatures (2D surface) and a single communal starting pressure, along a + 1D pressure profile. """ if reference_pressure is None: reference_pressure = pressure[0] diff --git a/src/thermo.cpp b/src/thermo.cpp index f79825f27d..ba1c7bd0f7 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -92,24 +92,21 @@ double MoistLapse(double pressure, double ref_temperature, double ref_pressure, return T1; // check final T1 P1 } -std::vector MoistLapseProfile(const std::vector& pressure_profile, +std::vector MoistLapseProfile(const std::vector& press_profile, double ref_temperature, - double ref_pressure) { - // Vectorized version of MoistLapse for a pressure profile. C++ internally use. + double ref_pressure, + int nstep) { + // MoistLapse for one full pressure profile given one ref_temperature. C++ internally use. // calculate temperature profile of an air parcel lifting saturated adiabatically // through the given pressure profile. - std::vector temperature_profile; - temperature_profile.reserve(pressure_profile.size()); + std::vector temp_profile; + temp_profile.reserve(press_profile.size()); - //double T1 = ref_temperature; - //double P1 = ref_pressure; - double T; - for (size_t i = 0; i < pressure_profile.size(); ++i) { - T = MoistLapse(pressure_profile[i], ref_temperature, ref_pressure, 50); - temperature_profile.push_back(T); + for (double p : press_profile) { + temp_profile.push_back(MoistLapse(p, ref_temperature, ref_pressure, nstep)); } - return temperature_profile; + return temp_profile; } std::pair LCL(double pressure, double temperature, double dewpoint) { @@ -198,7 +195,7 @@ ParProStruct _ParcelProfileHelper(const std::vector& pressure, double te press_upper.push_back(p); } } - std::vector temp_upper = MoistLapseProfile(press_upper, temp_lower.back(), press_lcl); + std::vector temp_upper = MoistLapseProfile(press_upper, temp_lower.back(), press_lcl, 30); press_lower.pop_back(); temp_lower.pop_back(); diff --git a/src/thermo.hpp b/src/thermo.hpp index f4f3c0379e..2623bc366b 100644 --- a/src/thermo.hpp +++ b/src/thermo.hpp @@ -23,10 +23,10 @@ std::vector DryLapseProfile(const std::vector& pressure_profile, double CaldlnTdlnP(double temperature, double pressure); double MoistLapse(double pressure, double ref_temperature, double ref_pressure, int nstep); -std::vector MoistLapseProfile(const std::vector& pressure_profile, +std::vector MoistLapseProfile(const std::vector& press_profile, double ref_temperature, - double ref_pressure); - + double ref_pressure, + int nstep); std::pair LCL(double pressure, double temperature, double dewpoint); From e0d1f1e8904bd5c294c2b9f975b3f2a25925d62f Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Thu, 26 Jun 2025 08:24:00 -0600 Subject: [PATCH 33/41] put MoistLapseVectorized in thermo.cpp, no need to vectorize in pybind11 module --- src/calcmod.cpp | 90 ++++++++++++++++++++++++++----------------------- src/thermo.cpp | 53 ++++++++++++++++++++++++++--- src/thermo.hpp | 11 ++++-- 3 files changed, 104 insertions(+), 50 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 70aa270230..5e8013be14 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -44,51 +44,55 @@ PYBIND11_MODULE(_calc_mod, m) { "Calculate the temperature at pressure levels assuming dry adiabatic process.", py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure")); - m.def("moist_lapse", [](py::array_t pressure, - py::array_t ref_temperature, - double ref_pressure, - int rk_nstep) { - // This function calculates the moist adiabatic profile for multiple starting - // temperatures (2D surface) and a single communal starting pressure, along a - // 1D pressure profile. - - // --- Step 1: Prepare the C++ vector for pressure levels --- - if (pressure.ndim() != 1) { - throw std::runtime_error("Input 'pressure' array must be 1D."); - } - std::vector pressure_vec(pressure.data(), pressure.data() + pressure.size()); - - // --- Step 2: Ensure the reference temperature array is contiguous --- - auto ref_temp_contig = py::array::ensure(ref_temperature, py::array::c_style); - - // --- Step 3: Define the shape of the output array: (N+1) dimension--- - // Create a vector to hold the shape dimensions. - std::vector out_shape; - for(int i = 0; i < ref_temp_contig.ndim(); ++i) { - out_shape.push_back(ref_temp_contig.shape(i)); - } - ssize_t profile_len = pressure_vec.size(); - out_shape.push_back(profile_len); - - // Create the final output array with the correct N+1 dimensional shape. - auto out_array = py::array_t(out_shape); - - // --- Step 4: Get direct pointers to data buffers for fast access --- - const double* ref_temp_ptr = static_cast(ref_temp_contig.request().ptr); - double* out_array_ptr = out_array.mutable_data(); - ssize_t num_profiles = ref_temp_contig.size(); // Total number of profiles to calculate - - // --- Step 5: Loop through each reference temperature --- - for (ssize_t i = 0; i < num_profiles; ++i) { - for (ssize_t j = 0; j < profile_len; ++j) { - out_array_ptr[i * profile_len + j] = MoistLapse(pressure_vec[j], ref_temp_ptr[i], ref_pressure, rk_nstep); - } - } - - return out_array; - }, "Calculate the temperature along a pressure profile assuming saturated adiabatic process.", + m.def("moist_lapse", &MoistLapseVectorized, + "Calculate the temperature along a pressure profile assuming saturated adiabatic process.", py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure"), py::arg("rk_nstep")); +// m.def("moist_lapse", [](py::array_t pressure, +// py::array_t ref_temperature, +// double ref_pressure, +// int rk_nstep) { +// // This function calculates the moist adiabatic profile for multiple starting +// // temperatures (2D surface) and a single communal starting pressure, along a +// // 1D pressure profile. +// +// // --- Step 1: Prepare the C++ vector for pressure levels --- +// if (pressure.ndim() > 1) { +// throw std::runtime_error("Input 'pressure' must be 1D array or a single value."); +// } +// std::vector pressure_vec(pressure.data(), pressure.data() + pressure.size()); +// +// // --- Step 2: Ensure the reference temperature array is contiguous --- +// auto ref_temp_contig = py::array::ensure(ref_temperature, py::array::c_style); +// +// // --- Step 3: Define the shape of the output array: (N+1) dimension--- +// // Create a vector to hold the shape dimensions. +// std::vector out_shape; +// for(int i = 0; i < ref_temp_contig.ndim(); ++i) { +// out_shape.push_back(ref_temp_contig.shape(i)); +// } +// ssize_t profile_len = pressure_vec.size(); +// out_shape.push_back(profile_len); +// +// // Create the final output array with the correct N+1 dimensional shape. +// auto out_array = py::array_t(out_shape); +// +// // --- Step 4: Get direct pointers to data buffers for fast access --- +// const double* ref_temp_ptr = static_cast(ref_temp_contig.request().ptr); +// double* out_array_ptr = out_array.mutable_data(); +// ssize_t num_profiles = ref_temp_contig.size(); // Total number of profiles to calculate +// +// // --- Step 5: Loop through each reference temperature --- +// for (ssize_t i = 0; i < num_profiles; ++i) { +// for (ssize_t j = 0; j < profile_len; ++j) { +// out_array_ptr[i * profile_len + j] = MoistLapse(pressure_vec[j], ref_temp_ptr[i], ref_pressure, rk_nstep); +// } +// } +// +// return out_array; +// }, "Calculate the temperature along a pressure profile assuming saturated adiabatic process.", +// py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure"), py::arg("rk_nstep")); + m.def("lcl", [](py::array_t pressure, py::array_t temperature, diff --git a/src/thermo.cpp b/src/thermo.cpp index ba1c7bd0f7..0ef694251f 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -3,6 +3,7 @@ #include #include // For std::pair #include +#include #include #include // for std::cerr #include // for std::numeric_limits @@ -69,17 +70,17 @@ double CaldlnTdlnP(double temperature, double pressure) { return dlnT_dlnP_Bakhshaii2013; } -double MoistLapse(double pressure, double ref_temperature, double ref_pressure, int nstep) { +double MoistLapse(double pressure, double ref_temperature, double ref_pressure, int rk_nstep) { // calculate temperature at pressure given reference temperature and pressure // assuming moist adiabatic expansion (vapor condenses and removed from the air // parcel) - double dlnP = log(pressure / ref_pressure) / (double)nstep; + double dlnP = log(pressure / ref_pressure) / (double)rk_nstep; double T1 = ref_temperature; double P1 = ref_pressure; double k[4]; - for (int i = 0; i < nstep; ++i) { + for (int i = 0; i < rk_nstep; ++i) { k[0] = CaldlnTdlnP(T1, P1); k[1] = CaldlnTdlnP(T1 * exp(k[0] * dlnP/2.), P1 * exp(dlnP/2.)); k[2] = CaldlnTdlnP(T1 * exp(k[1] * dlnP/2.), P1 * exp(dlnP/2.)); @@ -95,7 +96,7 @@ double MoistLapse(double pressure, double ref_temperature, double ref_pressure, std::vector MoistLapseProfile(const std::vector& press_profile, double ref_temperature, double ref_pressure, - int nstep) { + int rk_nstep) { // MoistLapse for one full pressure profile given one ref_temperature. C++ internally use. // calculate temperature profile of an air parcel lifting saturated adiabatically @@ -104,11 +105,53 @@ std::vector MoistLapseProfile(const std::vector& press_profile, temp_profile.reserve(press_profile.size()); for (double p : press_profile) { - temp_profile.push_back(MoistLapse(p, ref_temperature, ref_pressure, nstep)); + temp_profile.push_back(MoistLapse(p, ref_temperature, ref_pressure, rk_nstep)); } return temp_profile; } +py::array_t MoistLapseVectorized(py::array_t pressure, + py::array_t ref_temperature, + double ref_pressure, + int rk_nstep) { + // This function calculates the moist adiabatic profile for multiple starting + // temperatures (2D surface) and a single communal starting pressure, along a + // 1D pressure profile. + // --- Step 1: Prepare the C++ vector for pressure levels --- + if (pressure.ndim() > 1) { + throw std::runtime_error("Input 'pressure' must be 1D array or a single value."); + } + std::vector pressure_vec(pressure.data(), pressure.data() + pressure.size()); + + // --- Step 2: Ensure the reference temperature array is contiguous --- + auto ref_temp_contig = py::array::ensure(ref_temperature, py::array::c_style); + + // --- Step 3: Define the shape of the output array: (N+1) dimension--- + std::vector out_shape; + for(int i = 0; i < ref_temp_contig.ndim(); ++i) { + out_shape.push_back(ref_temp_contig.shape(i)); + } + ssize_t profile_len = pressure_vec.size(); + out_shape.push_back(profile_len); + + auto out_array = py::array_t(out_shape); + + // --- Step 4: Get direct pointers to data buffers for fast access --- + const double* ref_temp_ptr = static_cast(ref_temp_contig.request().ptr); + double* out_array_ptr = out_array.mutable_data(); + ssize_t num_profiles = ref_temp_contig.size(); + + // --- Step 5: Loop through each reference temperature --- + for (ssize_t i = 0; i < num_profiles; ++i) { + for (ssize_t j = 0; j < profile_len; ++j) { + out_array_ptr[i * profile_len + j] = MoistLapse(pressure_vec[j], ref_temp_ptr[i], ref_pressure, rk_nstep); + } + } + + return out_array; +} + + std::pair LCL(double pressure, double temperature, double dewpoint) { if (temperature <= dewpoint) { std::cerr << "Warning in function '" << __func__ diff --git a/src/thermo.hpp b/src/thermo.hpp index 2623bc366b..b1632333de 100644 --- a/src/thermo.hpp +++ b/src/thermo.hpp @@ -4,8 +4,11 @@ #include #include #include // For std::pair +#include +#include #include "constants.hpp" +namespace py = pybind11; namespace mc = metpy_constants; double MoistAirGasConstant(double specific_humidity); @@ -22,11 +25,15 @@ std::vector DryLapseProfile(const std::vector& pressure_profile, double ref_pressure); double CaldlnTdlnP(double temperature, double pressure); -double MoistLapse(double pressure, double ref_temperature, double ref_pressure, int nstep); +double MoistLapse(double pressure, double ref_temperature, double ref_pressure, int rk_nstep); std::vector MoistLapseProfile(const std::vector& press_profile, double ref_temperature, double ref_pressure, - int nstep); + int rk_nstep); +py::array_t MoistLapseVectorized(py::array_t pressure, + py::array_t ref_temperature, + double ref_pressure, + int rk_nstep); std::pair LCL(double pressure, double temperature, double dewpoint); From 853a213634afdd9b811b6b2984254890a0d150b8 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Thu, 26 Jun 2025 08:58:11 -0600 Subject: [PATCH 34/41] add DryLapseVectorized in thermo.cpp, no need to vectorize in pybind11 module --- src/calcmod.cpp | 49 ++----------------------------------------------- src/thermo.cpp | 41 +++++++++++++++++++++++++++++++++++++++++ src/thermo.hpp | 3 +++ 3 files changed, 46 insertions(+), 47 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 5e8013be14..f6f1c1aca0 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -40,59 +40,14 @@ PYBIND11_MODULE(_calc_mod, m) { "Calculate relative humidity from temperature and dewpoint.", py::arg("temperature"), py::arg("dewpoint"), py::arg("phase")); - m.def("dry_lapse", py::vectorize(DryLapse), - "Calculate the temperature at pressure levels assuming dry adiabatic process.", + m.def("dry_lapse", &DryLapseVectorized, + "Calculate the temperature along a pressure profile assuming dry adiabatic process.", py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure")); m.def("moist_lapse", &MoistLapseVectorized, "Calculate the temperature along a pressure profile assuming saturated adiabatic process.", py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure"), py::arg("rk_nstep")); -// m.def("moist_lapse", [](py::array_t pressure, -// py::array_t ref_temperature, -// double ref_pressure, -// int rk_nstep) { -// // This function calculates the moist adiabatic profile for multiple starting -// // temperatures (2D surface) and a single communal starting pressure, along a -// // 1D pressure profile. -// -// // --- Step 1: Prepare the C++ vector for pressure levels --- -// if (pressure.ndim() > 1) { -// throw std::runtime_error("Input 'pressure' must be 1D array or a single value."); -// } -// std::vector pressure_vec(pressure.data(), pressure.data() + pressure.size()); -// -// // --- Step 2: Ensure the reference temperature array is contiguous --- -// auto ref_temp_contig = py::array::ensure(ref_temperature, py::array::c_style); -// -// // --- Step 3: Define the shape of the output array: (N+1) dimension--- -// // Create a vector to hold the shape dimensions. -// std::vector out_shape; -// for(int i = 0; i < ref_temp_contig.ndim(); ++i) { -// out_shape.push_back(ref_temp_contig.shape(i)); -// } -// ssize_t profile_len = pressure_vec.size(); -// out_shape.push_back(profile_len); -// -// // Create the final output array with the correct N+1 dimensional shape. -// auto out_array = py::array_t(out_shape); -// -// // --- Step 4: Get direct pointers to data buffers for fast access --- -// const double* ref_temp_ptr = static_cast(ref_temp_contig.request().ptr); -// double* out_array_ptr = out_array.mutable_data(); -// ssize_t num_profiles = ref_temp_contig.size(); // Total number of profiles to calculate -// -// // --- Step 5: Loop through each reference temperature --- -// for (ssize_t i = 0; i < num_profiles; ++i) { -// for (ssize_t j = 0; j < profile_len; ++j) { -// out_array_ptr[i * profile_len + j] = MoistLapse(pressure_vec[j], ref_temp_ptr[i], ref_pressure, rk_nstep); -// } -// } -// -// return out_array; -// }, "Calculate the temperature along a pressure profile assuming saturated adiabatic process.", -// py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure"), py::arg("rk_nstep")); - m.def("lcl", [](py::array_t pressure, py::array_t temperature, diff --git a/src/thermo.cpp b/src/thermo.cpp index 0ef694251f..122c48a935 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -58,6 +58,47 @@ std::vector DryLapseProfile(const std::vector& pressure_profile, return temperature_profile; } + +py::array_t DryLapseVectorized(py::array_t pressure, + py::array_t ref_temperature, + double ref_pressure) { + // This function calculates the dry adiabatic profile for multiple starting + // temperatures (2D surface) and a single communal starting pressure, along a + // 1D pressure profile. + // --- Step 1: Prepare the C++ vector for pressure levels --- + if (pressure.ndim() > 1) { + throw std::runtime_error("Input 'pressure' must be 1D array or a single value."); + } + std::vector pressure_vec(pressure.data(), pressure.data() + pressure.size()); + + // --- Step 2: Ensure the reference temperature array is contiguous --- + auto ref_temp_contig = py::array::ensure(ref_temperature, py::array::c_style); + + // --- Step 3: Define the shape of the output array: (N+1) dimension--- + std::vector out_shape; + for(int i = 0; i < ref_temp_contig.ndim(); ++i) { + out_shape.push_back(ref_temp_contig.shape(i)); + } + ssize_t profile_len = pressure_vec.size(); + out_shape.push_back(profile_len); + + auto out_array = py::array_t(out_shape); + + // --- Step 4: Get direct pointers to data buffers for fast access --- + const double* ref_temp_ptr = static_cast(ref_temp_contig.request().ptr); + double* out_array_ptr = out_array.mutable_data(); + ssize_t num_profiles = ref_temp_contig.size(); + + // --- Step 5: Loop through each reference temperature --- + for (ssize_t i = 0; i < num_profiles; ++i) { + for (ssize_t j = 0; j < profile_len; ++j) { + out_array_ptr[i * profile_len + j] = DryLapse(pressure_vec[j], ref_temp_ptr[i], ref_pressure); + } + } + + return out_array; +} + double CaldlnTdlnP(double temperature, double pressure) { // Calculate dlnT/dlnP for a moist (saturated) adiabatic process. double rs = SaturationMixingRatio(pressure, temperature, "liquid"); diff --git a/src/thermo.hpp b/src/thermo.hpp index b1632333de..ea48f2ae0f 100644 --- a/src/thermo.hpp +++ b/src/thermo.hpp @@ -23,6 +23,9 @@ double DryLapse(double pressure, double ref_temperature, double ref_pressure); std::vector DryLapseProfile(const std::vector& pressure_profile, double ref_temperature, double ref_pressure); +py::array_t DryLapseVectorized(py::array_t pressure, + py::array_t ref_temperature, + double ref_pressure); double CaldlnTdlnP(double temperature, double pressure); double MoistLapse(double pressure, double ref_temperature, double ref_pressure, int rk_nstep); From e7679607451694f546f1cac38011212bccbaed49 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Fri, 27 Jun 2025 06:53:42 -0600 Subject: [PATCH 35/41] modify conditions for LCL to execute --- src/thermo.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/thermo.cpp b/src/thermo.cpp index 122c48a935..6aa15c91ac 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -194,7 +194,7 @@ py::array_t MoistLapseVectorized(py::array_t pressure, std::pair LCL(double pressure, double temperature, double dewpoint) { - if (temperature <= dewpoint) { + if (temperature < dewpoint) { std::cerr << "Warning in function '" << __func__ << "': Temperature must be greater than dew point for LCL calculation.\n"; return {std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()}; From 0afe0cc22533280243a9f0e98b8ed269aa66074e Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Mon, 30 Jun 2025 09:35:17 -0600 Subject: [PATCH 36/41] add c++ virtual_temperature_from_dewpoint --- src/calcmod.cpp | 6 ++++- src/metpy/calc/thermo.py | 53 +++++++++++++++++++++++++++++++++++++++- src/thermo.cpp | 7 ++++++ src/thermo.hpp | 4 ++- 4 files changed, 67 insertions(+), 3 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index f6f1c1aca0..fd99061c15 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -143,9 +143,13 @@ PYBIND11_MODULE(_calc_mod, m) { m.def("specific_humidity_from_dewpoint", py::vectorize(SpecificHumidityFromDewPoint), "Calculate the specific humidity from the dewpoint temperature and pressure.", py::arg("pressure"), py::arg("dewpoint"), py::arg("phase")); - m.def("virtual_temperature", py::vectorize(VirtualTemperature), "Calculate virtual temperature from temperature and mixing ratio.", py::arg("temperature"), py::arg("mixing_ratio"), py::arg("epsilon")); + + m.def("virtual_temperature_from_dewpoint", py::vectorize(VirtualTemperatureFromDewpoint), + "Calculate virtual temperature from dewpoint.", + py::arg("pressure"), py::arg("temperature"), py::arg("dewpoint"), py::arg("epsilon"), py::arg("phase")); + } diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 774c15240a..ed832f9dbb 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -2185,10 +2185,29 @@ def virtual_temperature( Renamed ``mixing`` parameter to ``mixing_ratio`` """ - # Calling c++ calculation module + return temperature * ((mixing_ratio + molecular_weight_ratio) + / (molecular_weight_ratio * (1 + mixing_ratio))) + +@exporter.export +@preprocess_and_wrap(wrap_like='temperature', broadcast=('temperature', 'mixing_ratio')) +@process_units( + { + 'temperature': '[temperature]', + 'mixing_ratio': '[dimensionless]', + 'molecular_weight_ratio': '[dimensionless]' + }, + '[temperature]', + ignore_inputs_for_output=('molecular_weight_ratio',) +) +def virtual_temperature_linfel( + temperature, mixing_ratio, molecular_weight_ratio=mpconsts.nounit.epsilon): + """ + Linfeng's version of 'virtual_temperature'. Added on Jun 30 2025 + """ return _calc_mod.virtual_temperature( temperature, mixing_ratio, molecular_weight_ratio) + @exporter.export @preprocess_and_wrap(wrap_like='temperature', broadcast=('pressure', 'temperature', @@ -2257,6 +2276,38 @@ def virtual_temperature_from_dewpoint( return virtual_temperature(temperature, mixing_ratio, molecular_weight_ratio) +@exporter.export +@preprocess_and_wrap(wrap_like='temperature', broadcast=('pressure', + 'temperature', + 'dewpoint')) +@process_units( + { + 'pressure': '[pressure]', + 'temperature': '[temperature]', + 'dewpoint': '[temperature]', + 'molecular_weight_ratio': '[dimensionless]' + }, + '[temperature]', + ignore_inputs_for_output=('molecular_weight_ratio',) +) +def virtual_temperature_from_dewpoint_linfel( + pressure, + temperature, + dewpoint, + molecular_weight_ratio=mpconsts.nounit.epsilon, + *, + phase='liquid' +): + """ + Linfeng's version of 'virtual_temperature_from_dewpoint'. Added on Jun 30 2025 + """ + return _calc_mod.virtual_temperature_from_dewpoint(pressure, + temperature, + dewpoint, + molecular_weight_ratio, + phase) + + @exporter.export @preprocess_and_wrap( wrap_like='temperature', diff --git a/src/thermo.cpp b/src/thermo.cpp index 6aa15c91ac..c96528b526 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -354,3 +354,10 @@ double SpecificHumidityFromDewPoint(double pressure, double dewpoint, std::strin double VirtualTemperature(double temperature, double mixing_ratio, double epsilon) { return temperature * (mixing_ratio + epsilon) / (epsilon * (1. + mixing_ratio)); } + +double VirtualTemperatureFromDewpoint(double pressure, double temperature, + double dewpoint, double epsilon, + std::string phase) { + double mixing_ratio = SaturationMixingRatio(pressure, dewpoint, phase); + return VirtualTemperature(temperature, mixing_ratio, epsilon); +} diff --git a/src/thermo.hpp b/src/thermo.hpp index ea48f2ae0f..1333c3a028 100644 --- a/src/thermo.hpp +++ b/src/thermo.hpp @@ -67,5 +67,7 @@ double SpecificHumidityFromMixingRatio(double mixing_ratio); double SpecificHumidityFromDewPoint(double pressure, double dew_point, std::string phase); double VirtualTemperature(double temperature, double mixing_ratio, double epsilon); - +double VirtualTemperatureFromDewpoint(double pressure, double temperature, + double dewpoint, double epsilon, + std::string phase); #endif // THERMO_HPP From 47845b53f317ae5b727a50e4cca2022aa195e661 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Mon, 7 Jul 2025 09:50:54 -0400 Subject: [PATCH 37/41] change ssize_t to size_t, ssize_t is not supported on windows --- src/calcmod.cpp | 4 ++-- src/thermo.cpp | 20 ++++++++++---------- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index fd99061c15..f742dfa974 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -76,7 +76,7 @@ PYBIND11_MODULE(_calc_mod, m) { auto t_lcl = py::array_t(p_contig.request().shape); // --- Step 3: Get the total number of elements to loop over --- - ssize_t size = p_contig.size(); + size_t size = p_contig.size(); // --- Step 4: Get direct pointers to the (now contiguous) data buffers --- const double* p_ptr = static_cast(p_contig.request().ptr); @@ -86,7 +86,7 @@ PYBIND11_MODULE(_calc_mod, m) { double* t_lcl_ptr = t_lcl.mutable_data(); // --- Step 5: Loop through the data as if it were a single flat 1D array --- - for (ssize_t i = 0; i < size; i++) { + for (size_t i = 0; i < size; i++) { // Call the scalar c++ function for each element std::pair result = LCL(p_ptr[i], t_ptr[i], d_ptr[i]); diff --git a/src/thermo.cpp b/src/thermo.cpp index c96528b526..578f1b6e5b 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -75,11 +75,11 @@ py::array_t DryLapseVectorized(py::array_t pressure, auto ref_temp_contig = py::array::ensure(ref_temperature, py::array::c_style); // --- Step 3: Define the shape of the output array: (N+1) dimension--- - std::vector out_shape; + std::vector out_shape; for(int i = 0; i < ref_temp_contig.ndim(); ++i) { out_shape.push_back(ref_temp_contig.shape(i)); } - ssize_t profile_len = pressure_vec.size(); + size_t profile_len = pressure_vec.size(); out_shape.push_back(profile_len); auto out_array = py::array_t(out_shape); @@ -87,11 +87,11 @@ py::array_t DryLapseVectorized(py::array_t pressure, // --- Step 4: Get direct pointers to data buffers for fast access --- const double* ref_temp_ptr = static_cast(ref_temp_contig.request().ptr); double* out_array_ptr = out_array.mutable_data(); - ssize_t num_profiles = ref_temp_contig.size(); + size_t num_profiles = ref_temp_contig.size(); // --- Step 5: Loop through each reference temperature --- - for (ssize_t i = 0; i < num_profiles; ++i) { - for (ssize_t j = 0; j < profile_len; ++j) { + for (size_t i = 0; i < num_profiles; ++i) { + for (size_t j = 0; j < profile_len; ++j) { out_array_ptr[i * profile_len + j] = DryLapse(pressure_vec[j], ref_temp_ptr[i], ref_pressure); } } @@ -168,11 +168,11 @@ py::array_t MoistLapseVectorized(py::array_t pressure, auto ref_temp_contig = py::array::ensure(ref_temperature, py::array::c_style); // --- Step 3: Define the shape of the output array: (N+1) dimension--- - std::vector out_shape; + std::vector out_shape; for(int i = 0; i < ref_temp_contig.ndim(); ++i) { out_shape.push_back(ref_temp_contig.shape(i)); } - ssize_t profile_len = pressure_vec.size(); + size_t profile_len = pressure_vec.size(); out_shape.push_back(profile_len); auto out_array = py::array_t(out_shape); @@ -180,11 +180,11 @@ py::array_t MoistLapseVectorized(py::array_t pressure, // --- Step 4: Get direct pointers to data buffers for fast access --- const double* ref_temp_ptr = static_cast(ref_temp_contig.request().ptr); double* out_array_ptr = out_array.mutable_data(); - ssize_t num_profiles = ref_temp_contig.size(); + size_t num_profiles = ref_temp_contig.size(); // --- Step 5: Loop through each reference temperature --- - for (ssize_t i = 0; i < num_profiles; ++i) { - for (ssize_t j = 0; j < profile_len; ++j) { + for (size_t i = 0; i < num_profiles; ++i) { + for (size_t j = 0; j < profile_len; ++j) { out_array_ptr[i * profile_len + j] = MoistLapse(pressure_vec[j], ref_temp_ptr[i], ref_pressure, rk_nstep); } } From 87b31eada8d4762b1b065dc840639d385564f1f6 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Mon, 7 Jul 2025 18:41:21 -0400 Subject: [PATCH 38/41] move lcl python binding from calcmod.cpp to thermo.cpp; debug syntax 'structured binding' --- src/calcmod.cpp | 54 +++-------------------------------------------- src/thermo.cpp | 56 +++++++++++++++++++++++++++++++++++++++++++++++-- src/thermo.hpp | 4 ++++ 3 files changed, 61 insertions(+), 53 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index f742dfa974..58ec3c7243 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -49,57 +49,9 @@ PYBIND11_MODULE(_calc_mod, m) { py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure"), py::arg("rk_nstep")); - m.def("lcl", [](py::array_t pressure, - py::array_t temperature, - py::array_t dewpoint) { - - // This helper ensures the arrays are in C-style contiguous memory. - // If an input array is already contiguous, it's a zero-cost operation. - // If it's a slice or has a different memory layout, it creates a copy. - // This makes the subsequent looping simple and safe. - auto p_contig = py::array::ensure(pressure, py::array::c_style); - auto t_contig = py::array::ensure(temperature, py::array::c_style); - auto d_contig = py::array::ensure(dewpoint, py::array::c_style); - - // --- Step 1: Check that all input arrays have the same shape --- - if (p_contig.ndim() != t_contig.ndim() || p_contig.ndim() != d_contig.ndim()) { - throw std::runtime_error("Input arrays must have the same number of dimensions."); - } - for (int i = 0; i < p_contig.ndim(); ++i) { - if (p_contig.shape(i) != t_contig.shape(i) || p_contig.shape(i) != d_contig.shape(i)) { - throw std::runtime_error("Input arrays must have the same shape."); - } - } - - // --- Step 2: Create output arrays with the exact same N-D shape as the inputs --- - auto p_lcl = py::array_t(p_contig.request().shape); - auto t_lcl = py::array_t(p_contig.request().shape); - - // --- Step 3: Get the total number of elements to loop over --- - size_t size = p_contig.size(); - - // --- Step 4: Get direct pointers to the (now contiguous) data buffers --- - const double* p_ptr = static_cast(p_contig.request().ptr); - const double* t_ptr = static_cast(t_contig.request().ptr); - const double* d_ptr = static_cast(d_contig.request().ptr); - double* p_lcl_ptr = p_lcl.mutable_data(); - double* t_lcl_ptr = t_lcl.mutable_data(); - - // --- Step 5: Loop through the data as if it were a single flat 1D array --- - for (size_t i = 0; i < size; i++) { - // Call the scalar c++ function for each element - std::pair result = LCL(p_ptr[i], t_ptr[i], d_ptr[i]); - - p_lcl_ptr[i] = result.first; - t_lcl_ptr[i] = result.second; - } - - // --- Step 6: Return a tuple of the two new, N-dimensional arrays --- - return std::make_tuple(p_lcl, t_lcl); - - }, "Calculate the lifting condensation level (LCL) from pressure, temperature and dewpoint.", - py::arg("pressure"), py::arg("temperature"), py::arg("dewpoint")); - + m.def("lcl", &LCLVectorized, + "Calculate the lifting condensation level (LCL) from pressure, temperature and dewpoint.", + py::arg("pressure"), py::arg("temperature"), py::arg("dewpoint")); m.def("parcel_profile", [](py::array_t pressure, double temperature, double dewpoint) { diff --git a/src/thermo.cpp b/src/thermo.cpp index 578f1b6e5b..4bf70e6932 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -1,6 +1,7 @@ #include #include #include +#include // For std::make_tuple #include // For std::pair #include #include @@ -212,10 +213,59 @@ std::pair LCL(double pressure, double temperature, double dewpoi double w_minus1 = lambert_wm1(pow(rh, 1.0 / a) * c * exp(c)); double t_lcl = c / w_minus1 * temperature; double p_lcl = pressure * pow(t_lcl / temperature, moist_heat_ratio); - + return {p_lcl, t_lcl}; } + +std::tuple, py::array_t> LCLVectorized(py::array_t pressure, + py::array_t temperature, + py::array_t dewpoint) { + + // This helper ensures the arrays are in C-style contiguous memory. + // If an input array is already contiguous, it's a zero-cost operation. + // If it's a slice or has a different memory layout, it creates a copy. + // This makes the subsequent looping simple and safe. + auto p_contig = py::array::ensure(pressure, py::array::c_style); + auto t_contig = py::array::ensure(temperature, py::array::c_style); + auto d_contig = py::array::ensure(dewpoint, py::array::c_style); + + // --- Step 1: Check that all input arrays have the same shape --- + if (p_contig.ndim() != t_contig.ndim() || p_contig.ndim() != d_contig.ndim()) { + throw std::runtime_error("Input arrays must have the same number of dimensions."); + } + for (int i = 0; i < p_contig.ndim(); ++i) { + if (p_contig.shape(i) != t_contig.shape(i) || p_contig.shape(i) != d_contig.shape(i)) { + throw std::runtime_error("Input arrays must have the same shape."); + } + } + + // --- Step 2: Create output arrays with the exact same N-D shape as the inputs --- + auto p_lcl = py::array_t(p_contig.request().shape); + auto t_lcl = py::array_t(p_contig.request().shape); + + // --- Step 3: Get the total number of elements to loop over --- + size_t size = p_contig.size(); + + // --- Step 4: Get direct pointers to the (now contiguous) data buffers --- + const double* p_ptr = static_cast(p_contig.request().ptr); + const double* t_ptr = static_cast(t_contig.request().ptr); + const double* d_ptr = static_cast(d_contig.request().ptr); + double* p_lcl_ptr = p_lcl.mutable_data(); + double* t_lcl_ptr = t_lcl.mutable_data(); + + // --- Step 5: Loop through the data as if it were a single flat 1D array --- + for (size_t i = 0; i < size; i++) { + // Call the scalar c++ function for each element + std::pair result = LCL(p_ptr[i], t_ptr[i], d_ptr[i]); + p_lcl_ptr[i] = result.first; + t_lcl_ptr[i] = result.second; + } + + // --- Step 6: Return a tuple of the two new, N-dimensional arrays --- + return std::make_tuple(p_lcl, t_lcl); +} + bool _CheckPressure(const std::vector& pressure) { for (size_t i = 0; i + 1 < pressure.size(); ++i) { if (pressure[i] < pressure[i + 1]) { @@ -252,7 +302,9 @@ ParProStruct _ParcelProfileHelper(const std::vector& pressure, double te } // Find the LCL - auto [press_lcl, temp_lcl] = LCL(pressure[0], temperature, dewpoint); + std::pair lcl_result = LCL(pressure[0], temperature, dewpoint); + double press_lcl = lcl_result.first; + double temp_lcl = lcl_result.second; // Establish profile below LCL std::vector press_lower; diff --git a/src/thermo.hpp b/src/thermo.hpp index 1333c3a028..7af7f2c45d 100644 --- a/src/thermo.hpp +++ b/src/thermo.hpp @@ -4,6 +4,7 @@ #include #include #include // For std::pair +#include // For std::tuple #include #include #include "constants.hpp" @@ -39,6 +40,9 @@ py::array_t MoistLapseVectorized(py::array_t pressure, int rk_nstep); std::pair LCL(double pressure, double temperature, double dewpoint); +std::tuple, py::array_t> LCLVectorized(py::array_t pressure, + py::array_t temperature, + py::array_t dewpoint); bool _CheckPressure(const std::vector& pressure); From c93e6aa47014ac31da42b8517a7ef7e2f9b72748 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Fri, 11 Jul 2025 10:40:18 -0600 Subject: [PATCH 39/41] add DryLapseVectorized_3D --- src/calcmod.cpp | 4 ++++ src/thermo.cpp | 49 +++++++++++++++++++++++++++++++++++++++++++++++++ src/thermo.hpp | 3 +++ 3 files changed, 56 insertions(+) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 58ec3c7243..50929f8b59 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -44,6 +44,10 @@ PYBIND11_MODULE(_calc_mod, m) { "Calculate the temperature along a pressure profile assuming dry adiabatic process.", py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure")); + m.def("dry_lapse_3d", &DryLapseVectorized_3D, + "Calculate the temperature along multiple pressure profiles assuming dry adiabatic process.", + py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure")); + m.def("moist_lapse", &MoistLapseVectorized, "Calculate the temperature along a pressure profile assuming saturated adiabatic process.", py::arg("pressure"), py::arg("ref_temperature"), py::arg("ref_pressure"), py::arg("rk_nstep")); diff --git a/src/thermo.cpp b/src/thermo.cpp index 4bf70e6932..6b6e1ed43b 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -100,6 +100,55 @@ py::array_t DryLapseVectorized(py::array_t pressure, return out_array; } +py::array_t DryLapseVectorized_3D(py::array_t pressure, + py::array_t ref_temperature, + py::array_t ref_pressure) { + // --- Step 1: Ensure all input arrays are using a contiguous memory layout --- + auto p_contig = py::array::ensure(pressure, py::array::c_style); + auto ref_temp_contig = py::array::ensure(ref_temperature, py::array::c_style); + auto ref_press_contig = py::array::ensure(ref_pressure, py::array::c_style); + + // --- Step 2: Perform comprehensive shape validation --- + if (ref_temp_contig.ndim() != ref_press_contig.ndim()) { + throw std::runtime_error("Input 'ref_temperature' and 'ref_pressure' must have the same number of dimensions."); + } + if (p_contig.ndim() != ref_temp_contig.ndim() + 1) { + throw std::runtime_error("Input 'pressure' must have one more dimension than 'ref_temperature'."); + } + for (int i = 0; i < ref_temp_contig.ndim(); ++i) { + if (ref_temp_contig.shape(i) != ref_press_contig.shape(i) || + p_contig.shape(i+1) != ref_temp_contig.shape(i)) { + throw std::runtime_error("The horizontal dimensions of all input arrays must match."); + } + } + + // --- Step 3: Define the shape of the output array --- + // The output shape will be identical to the input pressure array's shape. + auto out_array = py::array_t(p_contig.request().shape); + + // --- Step 4: Get direct pointers to data buffers for fast access --- + const double* pressure_ptr = static_cast(p_contig.request().ptr); + const double* ref_temp_ptr = static_cast(ref_temp_contig.request().ptr); + const double* ref_press_ptr = static_cast(ref_press_contig.request().ptr); + double* out_array_ptr = out_array.mutable_data(); + + // --- Step 5: Define loop boundaries --- + size_t num_profiles = ref_temp_contig.size(); // Total number of horizontal points + size_t profile_len = p_contig.shape(0); // Length of the vertical dimension + + // --- Step 6: Loop through each horizontal point and its vertical profile --- + for (int j = 0; j < profile_len; ++j) { + for (int i = 0; i < num_profiles; ++i) { + // Calculate the index for the current point in the flattened arrays. + out_array_ptr[i+j*num_profiles] = DryLapse(pressure_ptr[i+j*num_profiles], + ref_temp_ptr[i], + ref_press_ptr[i]); + } + } + + return out_array; +} + double CaldlnTdlnP(double temperature, double pressure) { // Calculate dlnT/dlnP for a moist (saturated) adiabatic process. double rs = SaturationMixingRatio(pressure, temperature, "liquid"); diff --git a/src/thermo.hpp b/src/thermo.hpp index 7af7f2c45d..656a6b37c6 100644 --- a/src/thermo.hpp +++ b/src/thermo.hpp @@ -27,6 +27,9 @@ std::vector DryLapseProfile(const std::vector& pressure_profile, py::array_t DryLapseVectorized(py::array_t pressure, py::array_t ref_temperature, double ref_pressure); +py::array_t DryLapseVectorized_3D(py::array_t pressure, + py::array_t ref_temperature, + py::array_t ref_pressure); double CaldlnTdlnP(double temperature, double pressure); double MoistLapse(double pressure, double ref_temperature, double ref_pressure, int rk_nstep); From 57995ebabc25f7c2a11403d3947af2afd06f4323 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Wed, 16 Jul 2025 09:36:31 -0600 Subject: [PATCH 40/41] debug: number of dimension inconsistency when lcl has inputs of a scaler and an array with one element --- src/metpy/calc/thermo.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index ed832f9dbb..e387bf0a43 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -797,6 +797,7 @@ def lcl_linfel(pressure, temperature, dewpoint, max_iters=None, eps=None): """ Linfeng's version of 'lcl'. Added on Jun23 2025 """ + pressure, temperature, dewpoint = np.atleast_1d(pressure, temperature, dewpoint) p_lcl, t_lcl = _calc_mod.lcl(pressure, temperature, dewpoint) return p_lcl, t_lcl From 9dbc34871cd423e82073858cc048f835b260f95c Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Mon, 21 Jul 2025 09:53:33 -0600 Subject: [PATCH 41/41] code convention check --- src/calcmod.cpp | 4 ++-- src/math.cpp | 4 ++-- src/thermo.cpp | 4 ++-- src/thermo.hpp | 4 ++-- 4 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 50929f8b59..8c2cdbe64d 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -1,8 +1,8 @@ -#include "constants.hpp" -#include "thermo.hpp" #include #include #include +#include "constants.hpp" +#include "thermo.hpp" #include // For std::pair #include // For std::make_tuple #include diff --git a/src/math.cpp b/src/math.cpp index 2c0dc4f39d..12c89469a0 100644 --- a/src/math.cpp +++ b/src/math.cpp @@ -18,8 +18,8 @@ double lambert_wm1(double x, double tol, int max_iter) { for (int i = 0; i < max_iter; ++i) { double ew = std::exp(w); - double wew = w * ew; - double diff = (wew - x) / (ew * (w + 1) - ((w + 2) * (wew - x)) / (2 * w + 2)); + double w_ew = w * ew; + double diff = (w_ew - x) / (ew * (w + 1) - ((w + 2) * (w_ew - x)) / (2 * w + 2)); w -= diff; if (std::abs(diff) < tol) { return w; diff --git a/src/thermo.cpp b/src/thermo.cpp index 6b6e1ed43b..9f04ab82c8 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -1,10 +1,10 @@ +#include +#include #include #include #include #include // For std::make_tuple #include // For std::pair -#include -#include #include #include // for std::cerr #include // for std::numeric_limits diff --git a/src/thermo.hpp b/src/thermo.hpp index 656a6b37c6..555bbc437a 100644 --- a/src/thermo.hpp +++ b/src/thermo.hpp @@ -1,12 +1,12 @@ #ifndef THERMO_HPP // if not defined #define THERMO_HPP // define the header file +#include +#include #include #include #include // For std::pair #include // For std::tuple -#include -#include #include "constants.hpp" namespace py = pybind11;