Implicit Euler Method for Condensation¶

Header file: <libs/superdrops/impliciteuler.hpp> [source]

struct ImplicitIterations¶

Struct for performing iterations of the Implicit Euler Method.

This struct defines parameters and methods for performing an iterations of the implicit method.

_Note: abbreviation NR = Newton Raphson (Method)

Public Functions

inline ImplicitIterations(const size_t maxniters, const double rtol, const double atol)¶

Constructor for ImplicitIterations class.

Parameters:

maxniters – Maximum no. iterations of Newton Raphson Method.
rtol – Relative tolerance for implicit Euler method.
atol – Absolute tolerance for implicit Euler method.

double integrate_condensation_ode(const ODEConstants &odeconsts, const double subdelt, const double rprev, double ziter) const¶

Integrates the condensation / evaporation ODE for radius^2 from t -> t+ subdelt.

Employs the Implicit Euler method (with potential sub-timestepping based on uniqueness criteria of Matsushima et. al), 2023 to forward timestep previous radius ‘rprev’ by subdelt according to the condensation/evaporation ODE. Implict timestepping equation defined in section 5.1.2 of Shima et al. 2009 and is root of polynomial g(z) = 0, where z = [R_i(t+delt)]^squared.

Uses at least 2 iterations of the Newton Raphson method and then checks if convergence criteria has been met (if a root of the g(Z) polynomial has been converged upon), else performs upto maxniters number of further iterations, checking for convergence after each one.

Uses at least ‘niters’ iterations of the Newton Raphson method and then checks if convergence criteria has been met (if a root of the g(Z) polynomial has been converged upon), else performs upto maxniters number of further iterations, checking for convergence after each one.

Parameters:

odeconsts – Constants of ODE during integration
subdelt – Time over which to integrate ODE
rprev – Radius of droplet at previous timestep.
ziter – Initial value for ziter.
odeconsts – Constants of ODE during integration
subdelt – Time over which to integrate ODE
rprev – Radius of droplet at previous timestep.

double initialguess(const ODEConstants &odeconsts, const double rprev) const¶

Returns appropriate initial guess (ie. a reasonable guess) for the Newton-Raphson method.

This method returns an initial guess for the Newton-Raphson method based on the given radius from the previous timestep and the current supersaturation ratio.

Guess is supposed to be a reasonable value for initial ‘ziter’ to use as first iteration of NR method in rootfinding algorithm for timestepping condensation/evaporation ODE. Here the guess criteria are as in SCALE-SDM for making initial guess for given droplet much greater than its (activation radius)^2 if the supersaturation > its activation supersaturation.

Parameters:

odeconsts – Constants of ODE during integration
rprev – Radius of droplet at previous timestep.

Returns:

Initial guess for ziter.

Private Functions

Kokkos::pair<double, bool> newtonraphson_niterations(const ODEConstants &odeconsts, const double subdelt, const double rprev, double ziter, const size_t niters) const¶

Performs niters number of Newton-Raphson iterations.

Function integrates (timesteps) condensation ODE by delt given initial guess for ziter, (which is usually radius^squared from previous timestep). Uses Newton Raphson iterative method with ‘niters’ number of iterations then returns updated ziter and boolean which is true if rootfinding has passed convergence test.

Parameters:

odeconsts – Constants of ODE during integration
subdelt – Time over which to integrate ODE
rprev – Radius at previous timestep
ziter – The current guess for ziter.
niters – Number of iterations of NR method to perform

Returns:

The updated value of ziter.

double newtonraphson_untilconverged(const ODEConstants &odeconsts, const size_t niterslimit, const double subdelt, const double rprev, double ziter) const¶

Performs Newton-Raphson iterations until convergence or maximum number of iterations is reached.

After every iteration, convergence criteria is tested and error is raised if method does not converge within ‘niterslimit’ iterations. Otherwise once convergence test is passed, function returns the new value for the ziter (which is the radius^2 at timestep ‘t+delt’). Refer to section 5.1.2 Shima et al. 2009 and section 3.3.3 of Matsushima et al. 2023 for more details.

After every iteration, convergence criteria is tested and error is raised if method does not converge within ‘maxniters’ iterations. Otherwise once convergence test is passed, function returns the new value for the ziter (which is the radius^2 at timestep ‘t+delt’). Refer to section 5.1.2 Shima et al. 2009 and section 3.3.3 of Matsushima et al. 2023 for more details.

Parameters:

odeconsts – Constants of ODE during integration
subdelt – Time over which to integrate ODE
niterslimit – The maxiumum number of iterations to attempt.
rprev – Radius at the previous timestep.
ziter – The current guess for ziter.
odeconsts – Constants of ODE during integration
subdelt – Time over which to integrate ODE
niterslimit – The maxiumum number of iterations to attempt.
rprev – Radius at the previous timestep.
ziter – The current guess for ziter.

Returns:

The updated value of ziter.

Returns:

The updated value of ziter.

Kokkos::pair<double, bool> iterate_rootfinding_algorithm(const ODEConstants &odeconsts, const double subdelt, const double rprev, double ziter) const¶

Perform one iteration of the Newton-Raphson rootfinding algorithm.

This function performs one iteration of the Newton-Raphson rootfinding algorithm, i.e. ziter^(m) -> ziter^(m+1) for iteration m+1 starting at m=1. Returns the updated value of ziter alongside a boolean which is true if the new value of ziter passes the convergence test.

Note: ziter is limited to >= 1e-8 so it’s always > 0.0

Parameters:

odeconsts – Constants of ODE during integration
subdelt – Time over which to integrate ODE
rprev – Radius at the previous timestep.
ziter – The current guess for ziter.

Returns:

A pair of the updated ziter and a boolean which is true if a root is converged upon.

double ode_gfunc(const ODEConstants &odeconsts, const double subdelt, const double rprev, const double rsqrd) const¶

Returns the value of g(z) / z * subdelt for the ODE.

This method computes the value of g(z) / z * subdelt used in the root-finding Newton-Raphson method for the dr/dt condensation / evaporation ODE.

ODE is for radial growth/shrink of each superdroplet due to condensation and diffusion of water vapour according to equations from “An Introduction To Clouds….” (see note at top of file).

Note: z = ziter = radius^2.

Parameters:

odeconsts – Constants of ODE during integration
subdelt – Time over which to integrate ODE
rprev – Radius at the previous timestep.
rsqrd – Current radius squared.

Returns:

RHS of g(z) / z * subdelt evaluted at rqrd.

double ode_gfuncderivative(const ODEConstants &odeconsts, const double subdelt, const double rsqrd) const¶

Returns the value of the derivative of g(z) with respect to z.

This method computes the value of dg(z)/dz * subdelt, where dg(z)/dz is the derivative of g(z) with respect to z=rsqr. g(z) is polynomial to find root of using Newton Raphson Method consistent as in ode_gfunc(…).

Parameters:

odeconsts – Constants of ODE during integration
subdelt – Time over which to integrate ODE
rsqrd – Current radius squared.

Returns:

RHS of dg(z)/dz * subdelt evaluted at rqrd.

inline bool check_for_convergence(const double gfunciter, const double gfuncprev) const¶

Returns true if the Newton-Raphson iterations have converged.

This method checks if the Newton-Raphson iterations have converged based on a standard local error test: |iteration - previous iteration| < RTOL * |iteration| + ATOL.

Parameters:

gfunciter – Value proportional to g(z) for current iteration
gfuncprev – Value proportional to g(z) for previous iteration

Returns:

Boolean=true if converged, false otherwise

Private Members

size_t maxniters¶: Maximum no. iterations of Newton Raphson Method

double rtol¶: Relative tolerance for convergence of NR method.

double atol¶: Absolute tolerance for convergence of NR method.

struct ODEConstants¶

Struct for constants of ODE during integration.

Public Members

double s_ratio¶: Supersaturation ratio.

double akoh¶: Kelvin factor in Kohler theory “a”.

double bkoh¶: Raoult factor in Kohler theory “b”.

double ffactor¶: Sum of heat and vapor diffusion factors.

class ImplicitEuler¶

Class for Implicit Euler (IE) integration of superdroplet condensational growth / evaporational shrinking ODE.

This class performs Implicit Euler integration of the superdroplet condensation / evaporation ODE using a Newton-Raphson root finding method to solve the implicit timestep equation of a stiff ODE.

Public Functions

inline ImplicitEuler(const double delt, const size_t maxniters, const double rtol, const double atol, const double minsubdelt)¶

Constructor for ImplicitEuler class.

Parameters:

delt – Time over which to integrate ODE using implcit Euler method.
maxniters – Maximum no. iterations of Newton Raphson Method.
rtol – Relative tolerance for implicit Euler method.
atol – Absolute tolerance for implicit Euler method.
minsubdelt – Minimum subtimestep in cases of substepping implicit Euler method.

double solve_condensation(const double s_ratio, const Kokkos::pair<double, double> kohler_ab, const double ffactor, const double rprev) const¶

Integrates the condensation / evaporation ODE employing the Implicit Euler method similarly to Matsushima et. al, 2023.

Forward timestep previous radius ‘rprev’ by delt using an Implicit Euler method (possibly with sub-timestepping) to integrate the condensation/evaporation ODE using fixed thermodynamics from the start of the timestep. Sub-timestepping employed when unique solution to g(Z) within required radius range is not guarenteed. Criteria as in appendix C of Matsushima et. al, 2023 except minimum sub-timestep is limited by minsubdelt.

Parameters:

s_ratio – The saturation ratio.
kohler_ab – A pair containing ‘a’ and ‘b’ factors for Kohler curve in that order.
ffactor – The sum of the diffusion factors.
rprev – Previous radius at time = t
s_ratio – The saturation ratio.
kohler_ab – A pair containing ‘a’ and ‘b’ factors for Kohler curve in that order.
ffactor – The sum of the diffusion factors.
rprev – Previous radius at time = t

Returns:

Updated radius for time = t + delt

Returns:

Updated radius for time = t + delt

Private Functions

bool first_unique_criteria(const ImplicitIterations::ODEConstants &odeconsts, const double rprev, const double ziter) const¶

Test of uniqueness criteria for un-activated droplets in environment with supersaturation less than its activation supersaturation.

Returns true if solution to g(Z) is guarenteed to be unique because it meets the uniquenes criteria of Case 2 from Matsushima et al. 2023 (see appendix C), namely that there is only one real root to g(Z) in the range 0 < Z < critical_R^2, where critical_R is the critical i.e. activation radius of the droplet. Here we use the less stringent constrain that S <= S_crit rather than S <= 1, and we ensure the current value for ziter is also less than the critical_R as it must be to guarentee solution in range 0 < R < critical_R is converged upon.

Parameters:

odeconsts – Constants of ODE during integration
rprev – Radius at previous timestep
ziter – Current guess for ziter.

Returns:

Boolean = true if solution is guarenteed to be unique.

inline double critial_timestep(const ImplicitIterations::ODEConstants &odeconsts) const¶

Calculates largest timestep which guarentees uniqueness of solution to g(Z) polynomial.

Returns the largest possible timestep that can be undertaken in which g(Z) has only one real root to g(Z) in the range 0 < Z < infinity. See Case 1 from Matsushima et al. 2023 (and derivation in appendix C).

Parameters:: odeconsts – Constants of ODE during integration
Returns:: Critical time step for unique solution.

inline bool second_unique_criteria(const ImplicitIterations::ODEConstants &odeconsts, const double subdelt) const¶

Test of uniqueness criteria for small enough timestep.

Returns true if solution to g(Z) is guarenteed to be unique because it meets the uniquenes criteria of Case 1 from Matsushima et al. 2023 (see appendix C), namely that the timestep is small enough to guarentee there is only one real root to g(Z) in the range 0 < Z < infinity.

Parameters:

odeconsts – Constants of ODE during integration
subdelt – Time over which to integrate ODE over.

Returns:

Boolean = true if solution is guarenteed to be unique.

double solve_with_adaptive_subtimestepping(const ImplicitIterations::ODEConstants &odeconsts, const double delt, double rprev, double ziter) const¶

Integrates the condensation / evaporation ODE employing the Implicit Euler method similarly to Matsushima et. al, 2023 with an adaptive timestepping subroutine.

Forward timestep previous radius ‘rprev’ by delt using an Implicit Euler method with sub-timestepping to integrate the condensation/evaporation ODE using fixed thermodynamics from the start of the timestep. Sub-timestepping employed to try to ensure uique solution to g(Z) as Matsushima et. al, 2023 except minimum sub-timestep is limited by minsubdelt. If critdelt < minsubdelt the uniqueness is not guarenteed. Reducing minsubdelt therefore increases the likelyhood of having a unique solution to g(Z), i.e. the accuracy of the solver is increased.

Parameters:

odeconsts – Constants of ODE during integration
delt – Time over which to integrate ODE over.
rprev – Previous radius at time = t
ziter – Initial guess for ziter.

Returns:

Updated radius^2 for time = t + delt

Private Members

double delt¶: Timestep of ODE solver (at each step implicit method is called).

double minsubdelt¶: Minimum subtimestep in cases of substepping

ImplicitIterations implit¶: Performs Newton Raphson Iterations of Implicit Method