Permalink
Browse files

sdm: adding the Kelvin term to the drop growth equation and to the eq…

…uilibrium radius calculations
  • Loading branch information...
1 parent d5133ef commit 23f7dcc4a4cbdcee39abf04a9f74a6796bca1210 Sylwester Arabas committed Nov 22, 2012
View
1 src/cmn/cmn_units.hpp
@@ -15,6 +15,7 @@
namespace si = boost::units::si;
using boost::units::quantity;
using boost::units::pow;
+using boost::units::root;
using boost::units::multiply_typeof_helper;
using boost::units::divide_typeof_helper;
using boost::units::power_typeof_helper;
View
104 src/phc/phc_kappa_koehler.hpp
@@ -7,40 +7,82 @@
*/
#pragma once
#include "phc.hpp"
+#include "phc_kelvin_term.hpp"
+#include <boost/math/tools/minima.hpp>
namespace phc
{
- /// @brief equilibrium wet radius to the third power for a given:
- /// @arg dry radius to the third power
- /// @arg the solubility parameter kappa
- /// @arg ratio of abmient vapour density/pressure to saturation vapour density/pressure for pure water
- ///
- /// the formula stems from applying the kappa-Koehler relation
- /// (eq. 6 in @copydetails Petters_and_Kreidenweis_2007) to a stationarity
- /// condition for a vapour diffusion equation, which translates to
- /// zero differece of vapour density: rho_ambient - rho_surface = 0
- ///
- /// since rho_surface = rho_surface_pure_water * a(r_w, r_d, kappa)
- /// one can derive r_w as a function of r_d, kappa and the ratio
- /// of abmient and surface vapour densities
- ///
- /// for the kappa-Koehler parameterisation rw3 is linear with rd3
- phc_declare_funct_macro quantity<si::volume, real_t> rw3_eq(
- quantity<si::volume, real_t> rd3,
- quantity<si::dimensionless, real_t> kappa,
- quantity<si::dimensionless, real_t> vap_ratio
- )
+ namespace kappa
{
- return rd3 * (1 - vap_ratio * (1 - kappa)) / (1 - vap_ratio);
- }
+ /// @brief equilibrium wet radius to the third power, with the Kelvin term discarded, for a given:
+ /// @arg dry radius to the third power
+ /// @arg the solubility parameter kappa
+ /// @arg ratio of abmient vapour density/pressure to saturation vapour density/pressure for pure water
+ ///
+ /// the formula stems from applying the kappa-Koehler relation
+ /// (eq. 6 in @copydetails Petters_and_Kreidenweis_2007) to a stationarity
+ /// condition for a vapour diffusion equation, which translates to
+ /// zero differece of vapour density: rho_ambient - rho_surface = 0
+ ///
+ /// since rho_surface = rho_surface_pure_water * a(r_w, r_d, kappa)
+ /// one can derive r_w as a function of r_d, kappa and the ratio
+ /// of abmient and surface vapour densities
+ ///
+ /// for the kappa-Koehler parameterisation rw3 is linear with rd3
+ phc_declare_funct_macro quantity<si::volume, real_t> rw3_eq_nokelvin(
+ quantity<si::volume, real_t> rd3,
+ quantity<si::dimensionless, real_t> kappa,
+ quantity<si::dimensionless, real_t> vap_ratio
+ )
+ {
+ return rd3 * (1 - vap_ratio * (1 - kappa)) / (1 - vap_ratio);
+ }
- /// @brief activity of water in solution (eqs. 1,6) in @copydetails Petters_and_Kreidenweis_2007
- phc_declare_funct_macro quantity<si::dimensionless, real_t> a_w(
- quantity<si::volume, real_t> rw3,
- quantity<si::volume, real_t> rd3,
- quantity<si::dimensionless, real_t> kappa
- )
- {
- return (rw3 - rd3) / (rw3 - rd3 * (real_t(1) - kappa));
- }
+
+ /// @brief activity of water in solution (eqs. 1,6) in @copydetails Petters_and_Kreidenweis_2007
+ phc_declare_funct_macro quantity<si::dimensionless, real_t> a_w(
+ quantity<si::volume, real_t> rw3,
+ quantity<si::volume, real_t> rd3,
+ quantity<si::dimensionless, real_t> kappa
+ )
+ {
+ return (rw3 - rd3) / (rw3 - rd3 * (real_t(1) - kappa));
+ }
+
+
+ // @brief equilibrium wet radius to the third power for a given:
+ /// @arg dry radius to the third power
+ /// @arg the solubility parameter kappa
+ /// @arg ratio of abmient vapour density/pressure to saturation vapour density/pressure for pure water
+ phc_declare_funct_macro quantity<si::volume, real_t> rw3_eq(
+ quantity<si::volume, real_t> rd3,
+ quantity<si::dimensionless, real_t> kappa,
+ quantity<si::dimensionless, real_t> vap_ratio,
+ quantity<si::temperature, real_t> T
+ )
+ {
+ // local functor to be passed to the minimisation func
+ struct f
+ {
+ quantity<si::dimensionless, real_t> vap_ratio;
+ quantity<si::volume, real_t> rd3;
+ quantity<si::dimensionless, real_t> kappa;
+ quantity<si::temperature, real_t> T;
+
+ real_t operator()(real_t rw3)
+ {
+ return vap_ratio
+ - a_w(rw3 * si::cubic_metres, rd3, kappa)
+ * phc::kelvin::klvntrm(pow(rw3, 1./3) * si::metres, T);
+ }
+ };
+
+ return real_t(boost::math::tools::brent_find_minima(
+ f({vap_ratio, rd3, kappa, T}), // the above-defined functor
+ real_t(rd3 / si::cubic_metres), // min
+ real_t(rw3_eq_nokelvin(rd3, kappa, vap_ratio) / si::cubic_metres), // max
+ sizeof(real_t) * 8 / 2 // the highest attainable precision with the algorithm according to Boost docs
+ ).first) * si::cubic_metres;
+ }
+ };
};
View
14 src/phc/phc_kelvin_term.hpp
@@ -17,11 +17,23 @@ namespace phc
// Kelvin term (for exapmle Kvorostyanov and Curry 2006)
+// TODO: document:
+// - why there's no exp()
+// - why there's no dependence on particle size!
+// - why the return type is length?
phc_declare_funct_macro quantity<si::length, real_t> A(
quantity<si::temperature, real_t> T
)
{
- return real_t(2.) * sg_surf<real_t>() / phc::R_v<real_t>() / T / phc::rho_w<real_t>();
+ return real_t(2.) * sg_surf<real_t>() / R_v<real_t>() / T / rho_w<real_t>();
+ }
+
+ phc_declare_funct_macro quantity<si::dimensionless, real_t> klvntrm(
+ quantity<si::length, real_t> r,
+ quantity<si::temperature, real_t> T
+ )
+ {
+ return exp(real_t(2.) * sg_surf<real_t>() / r / R_v<real_t>() / rho_w<real_t>() / T);
}
}
};
View
5 src/phc/phc_maxwell-mason.hpp
@@ -16,15 +16,16 @@ namespace phc
phc_declare_funct_macro quantity<divide_typeof_helper<si::area, si::time>::type, real_t> rdrdt(
const quantity<si::mass_density, real_t> rho_v, // ambient water vapour density
const quantity<si::temperature, real_t> T, // ambient temperature
- const quantity<si::dimensionless, real_t> a_w // water activity
+ const quantity<si::dimensionless, real_t> a_w, // water activity
+ const quantity<si::dimensionless, real_t> klvntrm // the Kelvin term
)
{
quantity<si::pressure, real_t> p_vs = phc::p_vs<real_t>(T);
quantity<divide_typeof_helper<si::energy, si::mass>::type, real_t> l_v = phc::l_v<real_t>(T);
quantity<divide_typeof_helper<si::energy, si::mass>::type, real_t> RvT = R_v<real_t>() * T;
return
( // vapour density difference (ambient - drop surface)
- rho_v - (p_vs / RvT * a_w) // TODO: the Kelvin term
+ rho_v - (p_vs / RvT * a_w) * klvntrm
)
/ phc::rho_w<real_t>() // water density
/ (
View
2 src/sdm/sdm.cpp
@@ -586,10 +586,8 @@ static void sd_diag(
for (auto moment : outmoments)
{
int k = moment.first, r = 0;
-std::cerr << "k=" << k << std::endl;
for (auto range : moment.second)
{
-std::cerr << "r in (" << range.first << " ... " << range.second << ")" << std::endl;
// zeroing the temporary var
thrust::fill(tmp_real.begin(), tmp_real.end(), 0); // TODO: is it needed
View
16 src/sdm/sdm_ode_cond.hpp
@@ -11,6 +11,7 @@
#include "../phc/phc_theta.hpp"
#include "../phc/phc_const_cp.hpp"
#include "../phc/phc_kappa_koehler.hpp"
+#include "../phc/phc_kelvin_term.hpp"
#include "../phc/phc_maxwell-mason.hpp"
namespace sdm
@@ -107,7 +108,7 @@ namespace sdm
void operator()(const thrust_size_t id)
{
thrust_size_t ij = stat.ij[id];
- real_t rw3_eq = phc::rw3_eq<real_t>( // TODO: allow choice among different Koehler curves
+ real_t rw3_eq = phc::kappa::rw3_eq<real_t>( // TODO: allow choice among different Koehler curves
stat.rd3[id] * si::cubic_metres,
real_t(stat.kpa[id]),
thrust::min(
@@ -116,7 +117,8 @@ namespace sdm
phc::R_v<real_t>() * (envi.T[ij] * si::kelvins) * (envi.rhod_rv[ij] * si::kilograms / si::cubic_metres)
/ (phc::p_vs<real_t>(envi.T[ij] * si::kelvins))
)
- )
+ ),
+ envi.T[ij] * si::kelvins
) / si::cubic_metres;
stat.xi[id] = this->xi_of_rw3(rw3_eq);
}
@@ -149,11 +151,12 @@ namespace sdm
real_t dxidt = phc::rdrdt<real_t>(
real_t(envi.rhod_rv[ij]) / si::cubic_metres * si::kilograms,
envi.T[ij] * si::kelvins,
- phc::a_w<real_t>(
+ phc::kappa::a_w<real_t>(
this->rw3_of_xi(stat.xi[id]) * si::cubic_metres, // rw3
stat.rd3[id] * si::cubic_metres, // rd3
real_t(stat.kpa[id]) // kappa
- )
+ ),
+ phc::kelvin::klvntrm<real_t>(this->rw_of_xi(stat.xi[id]) * si::metres, envi.T[ij] * si::kelvins) // the Kelvin term
)
* si::seconds / si::square_metres // to make it dimensionless
/ this->rw_of_xi(stat.xi[id]) // to get drdt
@@ -165,7 +168,7 @@ namespace sdm
// - Johnson 1980, JAS 37 2079-2085, disussion of inequality (5)
// - Heanel 1987, Beitr. Phys. Atmosph. 60 321-338, section 4.2
real_t xi_eq = this->xi_of_rw3(
- phc::rw3_eq<real_t>(
+ phc::kappa::rw3_eq<real_t>(
stat.rd3[id] * si::cubic_metres,
real_t(stat.kpa[id]),
real_t( // TODO: interpolation to drop positions?
@@ -174,7 +177,8 @@ namespace sdm
real_t(phc::R_v<real_t>() * (envi.T[ij] * si::kelvins) * (envi.rhod_rv[ij] * si::kilograms / si::cubic_metres)
/ (phc::p_vs<real_t>(envi.T[ij] * si::kelvins)))
)
- )
+ ),
+ envi.T[ij] * si::kelvins
) / si::cubic_metres);
// TODO: as an option
View
5 src/slv_parallel.hpp
@@ -95,6 +95,9 @@ class slv_parallel : public slv<real_t>
}
int n = 0;
+ // recording initial condition
+ record(sd, n, 0);
+
// adjustments for the initial condition
slvs[sd].apply_adjustments(n, setup.dt, true);
@@ -117,7 +120,7 @@ class slv_parallel : public slv<real_t>
bool fallback = (t == 0 && setup.advsch.time_levels() == 3);
n = fallback ? 0 : setup.advsch.time_levels() - 2;
- if (t % setup.nout == 0) record(sd, n, t / setup.nout);
+ if (t != 0 && t % setup.nout == 0) record(sd, n, t / setup.nout);
int last_group = -1;
for (int e = 0; e < setup.eqsys.n_vars(); ++e)

0 comments on commit 23f7dcc

Please sign in to comment.