Skip to content

Add kcl na2s #68

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
chengcli merged 5 commits into from
Jan 2, 2026
Merged

Add kcl na2s #68

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
eed1556
add kcl
chengcli Dec 17, 2025
68185a1
wip
chengcli Dec 21, 2025
65f0be7
wip
chengcli Jan 1, 2026
9d2cb6a
wip
chengcli Jan 2, 2026
8ba9a02
wip
chengcli Jan 2, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 10 additions & 6 deletions python/csrc/pythermo.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -75,25 +75,29 @@ void bind_thermo(py::module &m) {
.def(
"extrapolate_dlnp",
[](kintera::ThermoXImpl &self, torch::Tensor temp, torch::Tensor pres,
torch::Tensor xfrac, double dlnp, double ds_dlnp, bool verbose) {
torch::Tensor xfrac, double dlnp, double ds_dlnp, bool rainout,
bool verbose) {
kintera::ExtrapOptions opts;
opts.dlnp(dlnp).ds_dlnp(ds_dlnp).verbose(verbose);
opts.dlnp(dlnp).ds_dlnp(ds_dlnp).rainout(rainout).verbose(verbose);
return self.extrapolate_dlnp(temp, pres, xfrac, opts);
},
py::arg("temp"), py::arg("pres"), py::arg("xfrac"), py::arg("dlnp"),
py::arg("ds_dlnp") = 0., py::arg("verbose") = false)
py::arg("ds_dlnp") = 0., py::arg("rainout") = false,
py::arg("verbose") = false)

.def(
"extrapolate_dz",
[](kintera::ThermoXImpl &self, torch::Tensor temp, torch::Tensor pres,
torch::Tensor xfrac, double dz, double grav, double ds_dz,
bool verbose) {
bool rainout, bool verbose) {
kintera::ExtrapOptions opts;
opts.dz(dz).grav(grav).ds_dz(ds_dz).verbose(verbose);
opts.dz(dz).grav(grav).ds_dz(ds_dz).rainout(rainout).verbose(
verbose);
return self.extrapolate_dz(temp, pres, xfrac, opts);
},
py::arg("temp"), py::arg("pres"), py::arg("xfrac"), py::arg("dz"),
py::arg("grav"), py::arg("ds_dz") = 0., py::arg("verbose") = false);
py::arg("grav"), py::arg("ds_dz") = 0., py::arg("rainout") = false,
py::arg("verbose") = false);

m.def("relative_humidity", &kintera::relative_humidity, py::arg("temp"),
py::arg("conc"), py::arg("stoich"), py::arg("op"));
Expand Down
2 changes: 1 addition & 1 deletion src/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,8 +65,8 @@ add_library(kintera::kintera ALIAS ${namel}_${buildl})
if (CUDAToolkit_FOUND)
file(GLOB cu_src_files
*.cu
thermo/*.cu
utils/*.cu
thermo/*.cu
)

add_library(${namel}_cuda_${buildl}
Expand Down
20 changes: 12 additions & 8 deletions src/func_table.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,16 +8,20 @@ namespace kintera {
///// func1 registry /////

user_func1 func1_table_cpu[] = {
h2o_ideal, h2o_ideal_ddT, nh3_ideal, nh3_ideal_ddT,
nh3_h2s_lewis, nh3_h2s_lewis_ddT, h2s_ideal, h2s_ideal_ddT,
h2s_antoine, h2s_antoine_ddT, ch4_ideal, ch4_ideal_ddT,
so2_antoine, so2_antoine_ddT, co2_antoine, co2_antoine_ddT};
h2o_ideal, h2o_ideal_ddT, nh3_ideal, nh3_ideal_ddT,
nh3_h2s_lewis, nh3_h2s_lewis_ddT, h2s_ideal, h2s_ideal_ddT,
h2s_antoine, h2s_antoine_ddT, ch4_ideal, ch4_ideal_ddT,
so2_antoine, so2_antoine_ddT, co2_antoine, co2_antoine_ddT,
kcl_lodders, kcl_lodders_ddT, na_h2s_visscher, na_h2s_visscher_ddT};

std::vector<std::string> func1_names = {
"h2o_ideal", "h2o_ideal_ddT", "nh3_ideal", "nh3_ideal_ddT",
"nh3_h2s_lewis", "nh3_h2s_lewis_ddT", "h2s_ideal", "h2s_ideal_ddT",
"h2s_antoine", "h2s_antoine_ddT", "ch4_ideal", "ch4_ideal_ddT",
"so2_antoine", "so2_antoine_ddT", "co2_antoine", "co2_antoine_ddT"};
"h2o_ideal", "h2o_ideal_ddT", "nh3_ideal",
"nh3_ideal_ddT", "nh3_h2s_lewis", "nh3_h2s_lewis_ddT",
"h2s_ideal", "h2s_ideal_ddT", "h2s_antoine",
"h2s_antoine_ddT", "ch4_ideal", "ch4_ideal_ddT",
"so2_antoine", "so2_antoine_ddT", "co2_antoine",
"co2_antoine_ddT", "kcl_lodders", "kcl_lodders_ddT",
"na_h2s_visscher", "na_h2s_visscher_ddT"};

///// func2 registry /////

Expand Down
6 changes: 5 additions & 1 deletion src/func_table.cu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,11 @@ __device__ user_func1 func1_table_cuda[] = {
so2_antoine,
so2_antoine_ddT,
co2_antoine,
co2_antoine_ddT
co2_antoine_ddT,
kcl_lodders,
kcl_lodders_ddT,
na_h2s_visscher,
na_h2s_visscher_ddT
};

__device__ __constant__ user_func1* func1_table_device_ptr = func1_table_cuda;
Expand Down
1 change: 1 addition & 0 deletions src/kinetics/kinetics.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@ struct KineticsOptionsImpl final : public SpeciesThermoImpl {
YAML::Node const& config, bool verbose = false);

void report(std::ostream& os) const {
os << "-- kinetics options --\n";
os << "* Tref = " << Tref() << " K\n"
<< "* Pref = " << Pref() << " Pa\n"
<< "* offset_zero = " << (offset_zero() ? "true" : "false") << "\n"
Expand Down
83 changes: 33 additions & 50 deletions src/thermo/extrapolate_ad.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -79,11 +79,18 @@ void ThermoXImpl::extrapolate_dlnp(torch::Tensor temp, torch::Tensor pres,

int iter = 0;
pres *= exp(dlnp);
auto xfrac0 = xfrac.clone();

while (iter++ < options->max_iter()) {
xfrac.copy_(xfrac0);
auto gain = forward(temp, pres, xfrac);
if (opts.rainout()) {
xfrac.narrow(-1, options->vapor_ids().size(),
options->cloud_ids().size()) = 0.;
xfrac /= xfrac.sum(-1, true);
}

conc = compute("TPX->V", {temp, pres, xfrac});

auto cp_mole = effective_cp(temp, pres, xfrac, gain, conc);

entropy_vol = compute("TPV->S", {temp, pres, conc});
Expand Down Expand Up @@ -158,22 +165,24 @@ void ThermoXImpl::extrapolate_dz(torch::Tensor temp, torch::Tensor pres,

auto mmw = (mu * xfrac).sum(-1);
auto pres0 = pres.clone();
auto xfrac0 = xfrac.clone();

// using isothermal as an initial guess
pres.set_(pres0 * exp(-grav * mmw * dz / (constants::Rgas * temp)));

int iter = 0;
// std::vector<torch::Tensor> temp_list;
// std::vector<torch::Tensor> pres_list;
// std::vector<torch::Tensor> entropy_list;

while (iter++ < options->max_iter()) {
// temp_list.push_back(temp.clone());
// pres_list.push_back(pres.clone());

xfrac.copy_(xfrac0);
auto gain = forward(temp, pres, xfrac);
if (opts.rainout()) {
xfrac.narrow(-1, options->vapor_ids().size(),
options->cloud_ids().size()) = 0.;
xfrac /= xfrac.sum(-1, true);
}

conc = compute("TPX->V", {temp, pres, xfrac});
auto cp_mole = effective_cp(temp, pres, xfrac, gain, conc);
auto entropy_mole = compute("TPV->S", {temp, pres, conc}) / conc.sum(-1);
// entropy_list.push_back(entropy_mole.clone());

if (verbose) {
std::cout << "Iter " << iter << std::endl;
Expand All @@ -199,52 +208,26 @@ void ThermoXImpl::extrapolate_dz(torch::Tensor temp, torch::Tensor pres,
break;
}

/*if (iter % 2 == 0) {
// weight by distance from previous iterate
auto w1 = entropy_list[0] - entropy_target;
auto w2 = entropy_list[1] - entropy_target;
w1 = w1 * w1 / (w1 * w1 + w2 * w2 + 1e-10);
w2 = 1. - w1;

temp.set_(w2 * temp_list[0] + w1 * temp_list[1]);
pres.set_(w2 * pres_list[0] + w1 * pres_list[1]);

if (verbose) {
std::cout << " Weighted averaging over last " << temp_list.size()
<< " iterates." << std::endl;
std::cout << " temp = [" << temp.min().item<double>() << ", "
<< temp.max().item<double>() << "] K" << std::endl;
std::cout << " pres = [" << pres.min().item<double>() << ", "
<< pres.max().item<double>() << "] Pa" << std::endl;
}

temp_list.clear();
pres_list.clear();
entropy_list.clear();
} else {*/
int sub_iter = 3;
auto pres1 = pres.clone();
auto temp1 = temp.clone();

// total gas mole fractions
auto xg = xfrac.narrow(-1, 0, options->vapor_ids().size()).sum(-1);
while (sub_iter-- > 0) {
auto rho = compute("V->D", {conc});
pres.set_(pres0 - 0.5 * (rho + rho0) * grav * dz);
auto dlnp = pres.log() - pres1.log();
temp.set_(temp1 * (1. + (entropy_target - entropy_mole +
xg * constants::Rgas * dlnp) /
cp_mole));
conc = compute("TPX->V", {temp, pres, xfrac});
if (verbose) {
std::cout << " Sub-iter: pres = [" << pres.min().item<double>() << ", "
<< pres.max().item<double>() << "] Pa" << std::endl;
std::cout << " temp = [" << temp.min().item<double>() << ", "
<< temp.max().item<double>() << "] K" << std::endl;
std::cout << " dlnp = [" << dlnp.min().item<double>() << ", "
<< dlnp.max().item<double>() << "]" << std::endl;
}
auto rho = compute("V->D", {conc});
pres.set_(pres0 - 0.5 * (rho + rho0) * grav * dz);
auto dlnp = pres.log() - pres1.log();
temp.set_(temp1 * (1. + (entropy_target - entropy_mole +
xg * constants::Rgas * dlnp) /
cp_mole));
conc = compute("TPX->V", {temp, pres, xfrac});
if (verbose) {
std::cout << " Sub-iter: pres = [" << pres.min().item<double>() << ", "
<< pres.max().item<double>() << "] Pa" << std::endl;
std::cout << " temp = [" << temp.min().item<double>() << ", "
<< temp.max().item<double>() << "] K" << std::endl;
std::cout << " dlnp = [" << dlnp.min().item<double>() << ", "
<< dlnp.max().item<double>() << "]" << std::endl;
}
//}
}

if (iter >= options->max_iter()) {
Expand Down
7 changes: 7 additions & 0 deletions src/thermo/thermo.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,7 @@ struct ThermoOptionsImpl final : public SpeciesThermoImpl {
bool verbose = false);

void report(std::ostream& os) const {
os << "-- thermodynamics options --\n";
os << "* Tref = " << Tref() << " K\n"
<< "* Pref = " << Pref() << " Pa\n"
<< "* max_iter = " << max_iter() << "\n"
Expand Down Expand Up @@ -84,13 +85,19 @@ struct ExtrapOptions {
// Logarithmic change in pressure (dlnp = ln(p_new / p_old))
ADD_ARG(double, dlnp) = 0.;
ADD_ARG(double, dz) = 0.;

// Gravitational acceleration (m/s^2), positive
ADD_ARG(double, grav) = 0;
ADD_ARG(double, ds_dlnp) = 0;

// Height change (m)
ADD_ARG(double, ds_dz) = 0;

// If true, print debug information
ADD_ARG(bool, verbose) = false;

// Drop condensates
ADD_ARG(bool, rainout) = false;
};

//! Mass Thermodynamics
Expand Down
7 changes: 3 additions & 4 deletions src/utils/utils_dispatch.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,8 +8,7 @@
#include <configure.h>

// kintera
#include <kintera/utils/user_funcs.hpp>

#include "user_funcs.hpp"
#include "utils_dispatch.hpp"

namespace kintera {
Expand Down Expand Up @@ -99,13 +98,13 @@ void call_func3_cpu(at::TensorIterator &iter,

namespace at::native {

#ifdef DISABLE_CUDA
// #ifdef DISABLE_CUDA

DEFINE_DISPATCH(call_func1);
DEFINE_DISPATCH(call_func2);
DEFINE_DISPATCH(call_func3);

#endif
// #endif

REGISTER_ALL_CPU_DISPATCH(call_func1, &kintera::call_func1_cpu);
REGISTER_ALL_CPU_DISPATCH(call_func2, &kintera::call_func2_cpu);
Expand Down
8 changes: 4 additions & 4 deletions src/utils/utils_dispatch.cu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,8 +5,8 @@
#include <c10/cuda/CUDAGuard.h>

// kintera
#include <kintera/utils/user_funcs.hpp>
#include <kintera/loops.cuh>
#include "user_funcs.hpp"
#include "utils_dispatch.hpp"

namespace kintera {
Expand Down Expand Up @@ -94,9 +94,9 @@ void call_func3_cuda(at::TensorIterator &iter, std::vector<std::string> const& f

namespace at::native {

DEFINE_DISPATCH(call_func1);
DEFINE_DISPATCH(call_func2);
DEFINE_DISPATCH(call_func3);
//DEFINE_DISPATCH(call_func1);
//DEFINE_DISPATCH(call_func2);
//DEFINE_DISPATCH(call_func3);

REGISTER_CUDA_DISPATCH(call_func1, &kintera::call_func1_cuda);
REGISTER_CUDA_DISPATCH(call_func2, &kintera::call_func2_cuda);
Expand Down
24 changes: 24 additions & 0 deletions src/vapors/vapor_functions.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -161,4 +161,28 @@ inline double co2_antoine_ddT(double T) {
return logsvp_antoine_ddT(T, B, C);
}

DISPATCH_MACRO
inline double kcl_lodders(double T) {
double logp = 7.611 - 11382. / T;
return log(1.E5) + logp;
}

DISPATCH_MACRO
inline double kcl_lodders_ddT(double T) { return 11382. / (T * T); }
Comment on lines +166 to +171
Copy link

Copilot AI Jan 2, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The function lacks documentation explaining the source of the constants (7.611, 11382.) and the units expected for T and the returned value. Consider adding a brief comment describing the Lodders formulation and parameter units.

Suggested change
double logp = 7.611 - 11382. / T;
return log(1.E5) + logp;
}
DISPATCH_MACRO
inline double kcl_lodders_ddT(double T) { return 11382. / (T * T); }
// Lodders-type parameterization of KCl saturation vapor pressure.
// T: temperature in Kelvin.
// The constants 7.611 (dimensionless intercept) and 11382. (in Kelvin)
// are empirical coefficients from the Lodders formulation for log(p_sat).
// This function returns the natural logarithm of the saturation vapor
// pressure p_sat expressed in Pascals.
double logp = 7.611 - 11382. / T;
return log(1.E5) + logp;
}
DISPATCH_MACRO
inline double kcl_lodders_ddT(double T) {
// Temperature derivative of the Lodders KCl log(p_sat) expression.
// T is in Kelvin; the returned value is d/dT [log(p_sat [Pa])].
return 11382. / (T * T);
}

Copilot uses AI. Check for mistakes.

DISPATCH_MACRO
inline double na_h2s_visscher(double T) {
// double log10p = 8.55 - 13889. / T - 0.5 * log10(pH2S / 1E5);
// return 1.E5 * pow(10., log10p);
double a = 8.55;
double b = 13889.;
return (15. + 2. * a - 2. * b / T) * log(10.);
Comment on lines +174 to +179
Copy link

Copilot AI Jan 2, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The commented-out code and the magic number 15. in the implementation need explanation. Add documentation clarifying the relationship between the commented formula and the actual implementation, and explain the derivation of the constant 15.

Copilot uses AI. Check for mistakes.
}

DISPATCH_MACRO
inline double na_h2s_visscher_ddT(double T) {
double b = 13889.;
return 2. * b * log(10.) / (T * T);
}

} // namespace kintera
Loading