diff --git a/python/csrc/pythermo.cpp b/python/csrc/pythermo.cpp index d525d5f..6eb74a8 100644 --- a/python/csrc/pythermo.cpp +++ b/python/csrc/pythermo.cpp @@ -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")); diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 90d5880..0b47286 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -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} diff --git a/src/func_table.cpp b/src/func_table.cpp index 41bf4de..c5bab1e 100644 --- a/src/func_table.cpp +++ b/src/func_table.cpp @@ -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 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 ///// diff --git a/src/func_table.cu b/src/func_table.cu index c7a49a7..c5593fe 100644 --- a/src/func_table.cu +++ b/src/func_table.cu @@ -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; diff --git a/src/kinetics/kinetics.hpp b/src/kinetics/kinetics.hpp index 224a998..f0fa3fa 100644 --- a/src/kinetics/kinetics.hpp +++ b/src/kinetics/kinetics.hpp @@ -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" diff --git a/src/thermo/extrapolate_ad.cpp b/src/thermo/extrapolate_ad.cpp index bd6f5a0..98936a0 100644 --- a/src/thermo/extrapolate_ad.cpp +++ b/src/thermo/extrapolate_ad.cpp @@ -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}); @@ -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 temp_list; - // std::vector pres_list; - // std::vector 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; @@ -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() << ", " - << temp.max().item() << "] K" << std::endl; - std::cout << " pres = [" << pres.min().item() << ", " - << pres.max().item() << "] 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() << ", " - << pres.max().item() << "] Pa" << std::endl; - std::cout << " temp = [" << temp.min().item() << ", " - << temp.max().item() << "] K" << std::endl; - std::cout << " dlnp = [" << dlnp.min().item() << ", " - << dlnp.max().item() << "]" << 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() << ", " + << pres.max().item() << "] Pa" << std::endl; + std::cout << " temp = [" << temp.min().item() << ", " + << temp.max().item() << "] K" << std::endl; + std::cout << " dlnp = [" << dlnp.min().item() << ", " + << dlnp.max().item() << "]" << std::endl; } - //} } if (iter >= options->max_iter()) { diff --git a/src/thermo/thermo.hpp b/src/thermo/thermo.hpp index 7afd91a..7edf731 100644 --- a/src/thermo/thermo.hpp +++ b/src/thermo/thermo.hpp @@ -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" @@ -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 diff --git a/src/utils/utils_dispatch.cpp b/src/utils/utils_dispatch.cpp index 7182a22..cc8c17a 100644 --- a/src/utils/utils_dispatch.cpp +++ b/src/utils/utils_dispatch.cpp @@ -8,8 +8,7 @@ #include // kintera -#include - +#include "user_funcs.hpp" #include "utils_dispatch.hpp" namespace kintera { @@ -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); diff --git a/src/utils/utils_dispatch.cu b/src/utils/utils_dispatch.cu index 37ae424..3edac15 100644 --- a/src/utils/utils_dispatch.cu +++ b/src/utils/utils_dispatch.cu @@ -5,8 +5,8 @@ #include // kintera -#include #include +#include "user_funcs.hpp" #include "utils_dispatch.hpp" namespace kintera { @@ -94,9 +94,9 @@ void call_func3_cuda(at::TensorIterator &iter, std::vector 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); diff --git a/src/vapors/vapor_functions.h b/src/vapors/vapor_functions.h index a509246..974b8e3 100644 --- a/src/vapors/vapor_functions.h +++ b/src/vapors/vapor_functions.h @@ -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); } + +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.); +} + +DISPATCH_MACRO +inline double na_h2s_visscher_ddT(double T) { + double b = 13889.; + return 2. * b * log(10.) / (T * T); +} + } // namespace kintera