From ec4446e5be229791cb88cb2740ac44bee596c863 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Wed, 28 May 2025 18:15:34 -0600 Subject: [PATCH 1/2] vectorize the function using pybind11/numpy.h, achieving handling numpy array, float and dtype --- src/calcmod.cpp | 48 ++++--------------------------------- src/virtual_temperature.cpp | 17 ++----------- src/virtual_temperature.hpp | 8 ++----- 3 files changed, 9 insertions(+), 64 deletions(-) diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 06d4d6098a..94300b2d3b 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,7 @@ 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("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/virtual_temperature.cpp b/src/virtual_temperature.cpp index 592365cce0..e80586887a 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -1,19 +1,6 @@ #include "virtual_temperature.hpp" //#include -std::vector VirtualTemperature( - const std::vector& temperature, - const std::vector& mixing_ratio, - double epsilon -) { - - 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)); - } - - 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..3c627418ec 100644 --- a/src/virtual_temperature.hpp +++ b/src/virtual_temperature.hpp @@ -1,12 +1,8 @@ #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 VirtualTemperature(double temperature, double mixing_ratio, double epsilon); #endif // VIRTUAL_TEMPERATURE_HPP From 504e0ef0c685f93645dad4a90b2997e026732b87 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Fri, 30 May 2025 16:41:59 -0600 Subject: [PATCH 2/2] add dewpoint, add constants.hpp, adding other functions WIP --- src/calcmod.cpp | 4 ++++ src/constants.hpp | 37 +++++++++++++++++++++++++++++++++++++ src/virtual_temperature.cpp | 23 +++++++++++++++++++++++ src/virtual_temperature.hpp | 2 ++ 4 files changed, 66 insertions(+) create mode 100644 src/constants.hpp diff --git a/src/calcmod.cpp b/src/calcmod.cpp index 94300b2d3b..b32aeb41be 100644 --- a/src/calcmod.cpp +++ b/src/calcmod.cpp @@ -15,6 +15,10 @@ PYBIND11_MODULE(_calc_mod, m) { m.def("add", &add, "Add two numbers"); // Unified binding with default epsilon + 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 e80586887a..0204a07066 100644 --- a/src/virtual_temperature.cpp +++ b/src/virtual_temperature.cpp @@ -1,6 +1,29 @@ +#include +#include "constants.hpp" #include "virtual_temperature.hpp" //#include +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; + + // 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); +} + 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 3c627418ec..4b871bebd0 100644 --- a/src/virtual_temperature.hpp +++ b/src/virtual_temperature.hpp @@ -3,6 +3,8 @@ //#include +double DewPoint(double vapor_pressure); + double VirtualTemperature(double temperature, double mixing_ratio, double epsilon); #endif // VIRTUAL_TEMPERATURE_HPP