Skip to content

API.md

Latest commit

 

History

History
1177 lines (895 loc) · 33.4 KB

API.md

File metadata and controls

1177 lines (895 loc) · 33.4 KB

libcloudph++ API Reference (generated by Claude 4.5; not checked for errors)

This document provides a comprehensive reference for the libcloudph++ programming interface, covering the three microphysics schemes and their usage patterns.

Table of Contents

  1. Overview
  2. Lagrangian Microphysics (lgrngn)
  3. Single-Moment Bulk Microphysics (blk_1m)
  4. Double-Moment Bulk Microphysics (blk_2m)
  5. Common Utilities

Overview

libcloudph++ provides three main microphysics schemes:

  • Lagrangian (lgrngn): Super-droplet method with detailed particle-scale processes
  • Single-moment bulk (blk_1m): Bulk scheme predicting mixing ratios only
  • Double-moment bulk (blk_2m): Bulk scheme predicting both mixing ratios and number concentrations

All schemes use SI units with Boost.Units for type safety. Template parameter real_t can be float or double.

Lagrangian Microphysics (lgrngn)

Namespace

namespace libcloudphxx::lgrngn

Backend Selection

enum backend_t { 
    undefined, 
    serial,      // Single-threaded CPU
    OpenMP,      // Multi-threaded CPU
    CUDA,        // Single GPU
    multi_CUDA   // Multiple GPUs
};

Factory Function

template <typename real_t>
particles_proto_t<real_t>* factory(
    const backend_t backend, 
    opts_init_t<real_t> opts_init
);

Description: Creates a particle system instance with the specified backend.

Parameters:

  • backend: Computational backend (serial, OpenMP, CUDA, multi_CUDA)
  • opts_init: Initialization options (see USER_OPTIONS.md)

Returns: Pointer to particles_proto_t<real_t> base class

Example:

#include <libcloudph++/lgrngn/factory.hpp>

// Configure initialization options
libcloudphxx::lgrngn::opts_init_t<double> opts_init;
opts_init.dt = 1.0;              // timestep [s]
opts_init.nx = 64;               // grid cells in x
opts_init.ny = 64;               // grid cells in y
opts_init.nz = 64;               // grid cells in z
opts_init.dx = 100.0;            // cell size x [m]
opts_init.dy = 100.0;            // cell size y [m]
opts_init.dz = 100.0;            // cell size z [m]
opts_init.sd_conc = 64;          // super-droplets per cell

// Create particle system
auto particles = libcloudphxx::lgrngn::factory<double>(
    libcloudphxx::lgrngn::CUDA,
    opts_init
);

Main Class: particles_t

template <typename real_t, backend_t backend>
struct particles_t : particles_proto_t<real_t>

Constructor

particles_t(opts_init_t<real_t> opts_init, int n_x_tot = 0);

Parameters:

  • opts_init: Initialization options
  • n_x_tot: Total number of cells in x-direction (for MPI, otherwise use opts_init.nx)

Initialization

void init(
    const arrinfo_t<real_t> th,          // potential temperature [K]
    const arrinfo_t<real_t> rv,          // water vapor mixing ratio [kg/kg]
    const arrinfo_t<real_t> rhod,        // dry air density [kg/m³]
    const arrinfo_t<real_t> p = arrinfo_t<real_t>(),            // pressure [Pa]
    const arrinfo_t<real_t> courant_x = arrinfo_t<real_t>(),   // Courant number in x
    const arrinfo_t<real_t> courant_y = arrinfo_t<real_t>(),   // Courant number in y
    const arrinfo_t<real_t> courant_z = arrinfo_t<real_t>(),   // Courant number in z
    const std::map<enum common::chem::chem_species_t, const arrinfo_t<real_t>> ambient_chem = {}
);

Description: Initialize the particle system with atmospheric state and place super-droplets in the domain.

Parameters:

  • th, rv, rhod: Thermodynamic state arrays (required)
  • p: Pressure field (required if not using anelastic approximation)
  • courant_x, courant_y, courant_z: Courant numbers for advection (optional)
  • ambient_chem: Ambient chemical species concentrations (for chemistry simulations)

Time-Stepping Methods

step_sync()
void step_sync(
    const opts_t<real_t> &opts,          // runtime options
    arrinfo_t<real_t> th,                // potential temperature [K]
    arrinfo_t<real_t> rv,                // water vapor mixing ratio [kg/kg]
    const arrinfo_t<real_t> rhod = arrinfo_t<real_t>(),         // dry air density [kg/m³]
    const arrinfo_t<real_t> courant_x = arrinfo_t<real_t>(),   // Courant number in x
    const arrinfo_t<real_t> courant_y = arrinfo_t<real_t>(),   // Courant number in y
    const arrinfo_t<real_t> courant_z = arrinfo_t<real_t>(),   // Courant number in z
    const arrinfo_t<real_t> diss_rate = arrinfo_t<real_t>(),   // TKE dissipation [m²/s³]
    std::map<enum common::chem::chem_species_t, arrinfo_t<real_t>> ambient_chem = {}
);

Description: Perform a synchronous time step including all microphysics processes. Updates th and rv arrays to reflect latent heat release and phase changes.

Parameters:

  • opts: Runtime options controlling which processes to execute (see USER_OPTIONS.md)
  • th, rv: Thermodynamic state (input/output)
  • rhod: Dry air density (required for anelastic models)
  • courant_x, courant_y, courant_z: Courant numbers for advection
  • diss_rate: TKE dissipation rate for turbulent collision enhancement
  • ambient_chem: Chemical species for aqueous chemistry

Processes (controlled by opts):

  • Advection (opts.adve)
  • Condensation/evaporation (opts.cond)
  • Coalescence (opts.coal)
  • Sedimentation (opts.sedi)
  • Subsidence (opts.subs)
  • Chemistry (opts.chem_switch)
  • Source/relaxation (opts.src, opts.rlx)
sync_in()
void sync_in(
    arrinfo_t<real_t> th,
    arrinfo_t<real_t> rv,
    const arrinfo_t<real_t> rhod = arrinfo_t<real_t>(),
    const arrinfo_t<real_t> courant_x = arrinfo_t<real_t>(),
    const arrinfo_t<real_t> courant_y = arrinfo_t<real_t>(),
    const arrinfo_t<real_t> courant_z = arrinfo_t<real_t>(),
    const arrinfo_t<real_t> diss_rate = arrinfo_t<real_t>(),
    std::map<enum common::chem::chem_species_t, arrinfo_t<real_t>> ambient_chem = {}
);

Description: Copy atmospheric state from host to device (for GPU backends) without performing microphysics.

step_cond()
void step_cond(
    const opts_t<real_t> &opts,
    arrinfo_t<real_t> th,
    arrinfo_t<real_t> rv,
    std::map<enum common::chem::chem_species_t, arrinfo_t<real_t>> ambient_chem = {}
);

Description: Perform only condensation/evaporation step. Used for operator splitting.

step_async()
void step_async(const opts_t<real_t> &opts);

Description: Perform asynchronous processes (advection, coalescence, sedimentation) that don't require immediate synchronization with Eulerian fields.

Diagnostic Methods

Super-Droplet Concentration
void diag_sd_conc();

Description: Calculate number of super-droplets per cell. Result available via outbuf().

Pressure and Temperature
void diag_pressure();     // Pressure [Pa]
void diag_temperature();  // Temperature [K]

Description: Diagnose pressure or temperature fields from current thermodynamic state.

Relative Humidity
void diag_RH();           // Relative humidity [%]

Description: Calculate relative humidity field.

Particle Size Range Selection
void diag_dry_rng(const real_t &r_min, const real_t &r_max);  // dry radius [m]
void diag_wet_rng(const real_t &r_min, const real_t &r_max);  // wet radius [m]
void diag_kappa_rng(const real_t &kappa_min, const real_t &kappa_max);  // hygroscopicity

Description: Select super-droplets within specified size or hygroscopicity range for subsequent moment calculations.

Parameters:

  • r_min, r_max: Radius range (dry or wet) in meters
  • kappa_min, kappa_max: Hygroscopicity parameter range
Consecutive Range Selection
void diag_dry_rng_cons(const real_t &r_min, const real_t &r_max);
void diag_wet_rng_cons(const real_t &r_min, const real_t &r_max);
void diag_kappa_rng_cons(const real_t &kappa_min, const real_t &kappa_max);

Description: Further filter previously selected super-droplets. Use immediately after diag_*_rng() calls.

Example:

// Select droplets with wet radius 0.5-1 μm AND kappa 0.1-0.2
particles->diag_wet_rng(0.5e-6, 1e-6);
particles->diag_kappa_rng_cons(0.1, 0.2);
particles->diag_wet_mom(3);  // 3rd moment of filtered droplets
Moment Calculations
void diag_dry_mom(const int &k);      // k-th moment of dry radius
void diag_wet_mom(const int &k);      // k-th moment of wet radius
void diag_kappa_mom(const int &k);    // k-th moment of hygroscopicity
void diag_incloud_time_mom(const int &k);  // k-th moment of in-cloud time

Description: Calculate k-th moment of particle size distribution (or other attribute). Result in outbuf().

Parameters:

  • k: Moment order (typically 0, 1, 2, 3 for concentration, mean radius, variance, mass)

Common moments:

  • k=0: Number concentration [#/m³]
  • k=1: Mean radius [m]
  • k=2: Variance-related
  • k=3: Mass/volume concentration [m³/m³] or [kg/kg]
Velocity Moments
void diag_up_mom(const int &k);  // k-th moment of u-velocity component
void diag_vp_mom(const int &k);  // k-th moment of v-velocity component
void diag_wp_mom(const int &k);  // k-th moment of w-velocity component
Mass Density
void diag_wet_mass_dens(const real_t &r_min, const real_t &r_max);

Description: Calculate total mass density of droplets in specified wet radius range.

Activation Diagnostics
void diag_rw_ge_rc();    // Select droplets with r_wet >= r_critical
void diag_RH_ge_Sc();    // Select droplets where RH >= critical supersaturation
void diag_all();         // Select all droplets
Precipitation
void diag_precip_rate();  // Precipitation rate [m/s]
void diag_max_rw();       // Maximum wet radius [m]
Velocity Divergence
void diag_vel_div();      // Velocity divergence [1/s]
Chemistry
void diag_chem(const enum common::chem::chem_species_t &species);

Description: Diagnose aqueous concentration of chemical species.

Surface Accumulation (Puddle)
std::map<libcloudphxx::common::output_t, real_t> diag_puddle();

Description: Get accumulated rainfall at surface.

Returns: Map with keys:

  • common::output_t::rl: Total liquid water at surface
  • common::output_t::rc: Cloud water at surface
  • common::output_t::rr: Rain water at surface
Direct Attribute Access
std::vector<real_t> get_attr(const std::string &attr_name);

Description: Get raw super-droplet attribute array.

Parameters:

  • attr_name: Attribute name (e.g., "rd", "rw", "kpa", "x", "y", "z", "vt")

Returns: Vector of attribute values for all super-droplets

Output Buffer Access
real_t* outbuf();

Description: Get pointer to output buffer containing most recent diagnostic result.

Returns: Pointer to array with dimensions matching grid (nx × ny × nz)

Usage Pattern:

particles->diag_wet_mom(3);
real_t *result = particles->outbuf();
// Copy result to your array
for (int i = 0; i < nx * ny * nz; ++i) {
    liquid_water_content[i] = result[i];
}

Data Structures

arrinfo_t

template <typename real_t>
struct arrinfo_t {
    real_t * const data;              // Pointer to array data
    const ptrdiff_t *strides;         // Array strides for multidimensional indexing
    const std::vector<ptrdiff_t> strvec;  // Optional: local storage for strides
    
    // Constructors
    arrinfo_t();  // Default: null array
    arrinfo_t(real_t *data, const ptrdiff_t *strides);
    arrinfo_t(real_t *data, const std::vector<ptrdiff_t> &strvec);
    
    // Methods
    bool is_null() const;
};

Description: Helper structure for passing multi-dimensional arrays to libcloudph++.

Parameters:

  • data: Pointer to contiguous or strided array
  • strides: Array of stride values for each dimension (in units of real_t)
  • strvec: Alternative: vector for local stride storage

Example:

// 1D array
std::vector<double> arr_1d(nx);
arrinfo_t<double> info_1d(arr_1d.data(), nullptr);

// 3D array with strides
double *arr_3d = new double[nx * ny * nz];
ptrdiff_t strides[3] = {1, nx, nx*ny};
arrinfo_t<double> info_3d(arr_3d, strides);

// Using vector for strides
std::vector<ptrdiff_t> str = {1, nx, nx*ny};
arrinfo_t<double> info_3d_v(arr_3d, str);

opts_init_t

See USER_OPTIONS.md for complete documentation of initialization options including:

  • Domain configuration (nx, ny, nz, dx, dy, dz)
  • Super-droplet configuration (sd_conc, sd_conc_mean)
  • Aerosol size distributions (dry_distros, dry_sizes)
  • Chemistry options (chem_switch, chem_rho)
  • GPU settings (dev_count, dev_id)

opts_t

See USER_OPTIONS.md for runtime options controlling:

  • Process switches (adve, sedi, cond, coal, chem_switch)
  • Timestep substepping (sstp_cond, sstp_coal, sstp_chem)
  • Boundary conditions (open_side_walls, periodic_topbot)

Single-Moment Bulk Microphysics (blk_1m)

Namespace

namespace libcloudphxx::blk_1m

Overview

The single-moment bulk scheme predicts mixing ratios for:

  • rc: Cloud water mixing ratio [kg/kg]
  • rr: Rain water mixing ratio [kg/kg]
  • ri: Ice mixing ratio [kg/kg]
  • rs: Snow mixing ratio [kg/kg]
  • rg: Graupel mixing ratio [kg/kg]

Saturation Adjustment

adj_cellwise_nwtrph()

template <typename real_t, class cont_t>
void adj_cellwise_nwtrph(
    const opts_t<real_t> &opts,
    const cont_t &rhod_cont,   // dry air density [kg/m³]
    const cont_t &p_cont,      // pressure [Pa]
    cont_t &th_cont,           // potential temperature [K] (input/output)
    cont_t &rv_cont,           // water vapor mixing ratio [kg/kg] (input/output)
    cont_t &rc_cont,           // cloud water mixing ratio [kg/kg] (input/output)
    const real_t &dt           // timestep [s]
);

Description: Saturation adjustment using Newton-Raphson method to solve for equilibrium between vapor and cloud water, accounting for latent heat release.

Parameters:

  • opts: Options controlling adjustment (see USER_OPTIONS.md)
  • rhod_cont: Dry air density (constant in anelastic models)
  • p_cont: Pressure field (required if opts.const_p == true)
  • th_cont: Potential temperature (updated with latent heat effects)
  • rv_cont: Water vapor mixing ratio (updated)
  • rc_cont: Cloud water mixing ratio (updated)
  • dt: Timestep duration

Notes:

  • Requires opts.const_p XOR opts.th_dry (exactly one must be true)
  • If opts.const_p == true: uses "standard" potential temperature
  • If opts.th_dry == true: uses dry potential temperature

adj_cellwise_constp()

template <typename real_t, class cont_t>
void adj_cellwise_constp(
    const opts_t<real_t> &opts,
    cont_t &th_cont,
    cont_t &rv_cont,
    cont_t &rc_cont,
    const cont_t &rhod_cont
);

Description: Saturation adjustment for constant pressure (simplified version).

Right-Hand Side Terms

rhs_cellwise()

template <typename real_t, class cont_t>
void rhs_cellwise(
    const opts_t<real_t> &opts,
    cont_t &dot_rc_cont,   // cloud water tendency [kg/kg/s] (output)
    cont_t &dot_rr_cont,   // rain water tendency [kg/kg/s] (output)
    const cont_t &rc_cont, // cloud water [kg/kg]
    const cont_t &rr_cont  // rain water [kg/kg]
);

Description: Calculate microphysical tendencies for cloud and rain water using Kessler (1969) parameterization.

Processes:

  • Autoconversion (if opts.conv == true): Cloud water → rain water

    • Formula: rate = k_acnv * max(rc - r_c0, 0)
    • Threshold: r_c0 (default 0.5 g/kg)
    • Time scale: 1/k_acnv (default 0.001 s)

  • Collection/Accretion (if opts.accr == true): Cloud water collected by rain

    • Formula: Kessler collection formula

Parameters:

  • opts: Options controlling processes
  • dot_rc_cont, dot_rr_cont: Output tendency arrays
  • rc_cont, rr_cont: Input mixing ratio arrays

Note: This version assumes rain evaporation is handled separately (requires opts.adj_nwtrph == false).

rhs_cellwise_revap()

template <typename real_t, class cont_t>
void rhs_cellwise_revap(
    const opts_t<real_t> &opts,
    cont_t &dot_th_cont,   // potential temperature tendency [K/s] (output)
    cont_t &dot_rv_cont,   // water vapor tendency [kg/kg/s] (output)
    cont_t &dot_rc_cont,   // cloud water tendency [kg/kg/s] (output)
    cont_t &dot_rr_cont,   // rain water tendency [kg/kg/s] (output)
    const cont_t &rhod_cont,
    const cont_t &p_cont,
    const cont_t &th_cont,
    const cont_t &rv_cont,
    const cont_t &rc_cont,
    const cont_t &rr_cont,
    const real_t &dt
);

Description: Calculate tendencies including rain evaporation (to be used with Newton-Raphson saturation adjustment).

Additional Process:

  • Rain evaporation: Parameterized evaporation of falling rain droplets

Parameters:

  • Same as rhs_cellwise() plus thermodynamic state for evaporation calculation

Note: Requires opts.adj_nwtrph == true.

Ice Processes

template <typename real_t, class cont_t>
void rhs_cellwise_all(
    const opts_t<real_t> &opts,
    cont_t &dot_th_cont,
    cont_t &dot_rv_cont,
    cont_t &dot_rc_cont,
    cont_t &dot_rr_cont,
    cont_t &dot_ri_cont,   // ice mixing ratio tendency
    cont_t &dot_rs_cont,   // snow mixing ratio tendency
    cont_t &dot_rg_cont,   // graupel mixing ratio tendency
    const cont_t &rhod_cont,
    const cont_t &p_cont,
    const cont_t &th_cont,
    const cont_t &rv_cont,
    const cont_t &rc_cont,
    const cont_t &rr_cont,
    const cont_t &ri_cont,
    const cont_t &rs_cont,
    const cont_t &rg_cont,
    const real_t &dt
);

Description: Full microphysics including ice processes (Grabowski 1999).

Processes:

  • Ice nucleation
  • Ice deposition growth
  • Riming
  • Melting
  • All warm rain processes

Options Structure

template <typename real_t>
struct opts_t {
    bool cond;      // condensation
    bool conv;      // autoconversion
    bool accr;      // accretion/collection
    bool revp;      // rain evaporation
    bool sedi;      // sedimentation
    
    // Additional ice options (for rhs_cellwise_all)
    bool nuc;       // ice nucleation
    bool dep;       // deposition
    bool rim;       // riming
    bool melt;      // melting
    
    // Thermodynamic options
    bool const_p;   // constant pressure (standard theta)
    bool th_dry;    // dry potential temperature
    
    // Parameters
    real_t r_c0;    // autoconversion threshold [kg/kg]
    real_t k_acnv;  // autoconversion rate [1/s]
    
    // Newton-Raphson options
    bool adj_nwtrph;  // use Newton-Raphson adjustment
    real_t eps;       // convergence tolerance
    int iters;        // max iterations
};

See USER_OPTIONS.md for defaults and details.

Double-Moment Bulk Microphysics (blk_2m)

Namespace

namespace libcloudphxx::blk_2m

Overview

Double-moment scheme predicts both mixing ratios and number concentrations:

  • rc, nc: Cloud water mass and number
  • rr, nr: Rain water mass and number

Right-Hand Side Calculation

rhs_cellwise()

template <typename real_t, class cont_t>
void rhs_cellwise(
    const opts_t<real_t> &opts,
    cont_t &dot_th_cont,    // potential temperature tendency [K/s]
    cont_t &dot_rv_cont,    // water vapor tendency [kg/kg/s]
    cont_t &dot_rc_cont,    // cloud water mass tendency [kg/kg/s]
    cont_t &dot_nc_cont,    // cloud droplet number tendency [1/kg/s]
    cont_t &dot_rr_cont,    // rain water mass tendency [kg/kg/s]
    cont_t &dot_nr_cont,    // rain drop number tendency [1/kg/s]
    const cont_t &rhod_cont,   // dry air density [kg/m³]
    const cont_t &th_cont,     // potential temperature [K]
    const cont_t &rv_cont,     // water vapor [kg/kg]
    const cont_t &rc_cont,     // cloud water [kg/kg]
    const cont_t &nc_cont,     // cloud droplet number [1/kg]
    const cont_t &rr_cont,     // rain water [kg/kg]
    const cont_t &nr_cont,     // rain drop number [1/kg]
    const real_t &dt,          // timestep [s]
    const cont_t &p_cont = cont_t()  // pressure [Pa] (if const_p)
);

Description: Calculate microphysical tendencies using double-moment parameterizations.

Processes (selected by options):

  1. Activation (opts.acti): Aerosol activation to form cloud droplets

    • Uses Twomey (1959) or other parameterizations
    • Requires aerosol specification via opts.dry_distros

  2. Autoconversion (opts.acnv): Cloud droplets → rain drops

    • Available formulations: Kessler, Beheng, Berry-Reinhardt
    • Affects both mass (rc→rr) and number (nc→nr)

  3. Accretion (opts.accr): Rain collecting cloud water

    • Separate parameterizations for mass and number

  4. Self-collection (opts.slfc): Droplet-droplet collisions within same category

    • Reduces number concentration while preserving mass

  5. Rain evaporation (opts.revp): Evaporation of rain below cloud

    • Ventilation effects included

  6. Sedimentation (opts.sedi): Gravitational settling

    • Mass-weighted and number-weighted fall speeds

Parameters:

  • opts: See USER_OPTIONS.md for complete options
  • Thermodynamic state: rhod, th, rv, p (if const_p)
  • Cloud state: rc, nc
  • Rain state: rr, nr
  • Output tendencies: dot_th, dot_rv, dot_rc, dot_nc, dot_rr, dot_nr

Notes:

  • Requires either opts.const_p == true OR opts.th_dry == true (not both)
  • If const_p == true: must provide p_cont and use standard potential temperature
  • If th_dry == true: uses dry potential temperature (for anelastic models)

Options Structure

template <typename real_t>
struct opts_t {
    // Process switches
    bool cond;      // condensation
    bool acti;      // activation
    bool acnv;      // autoconversion
    bool accr;      // accretion
    bool slfc;      // self-collection
    bool revp;      // rain evaporation
    bool sedi;      // sedimentation
    
    // Thermodynamics
    bool const_p;   // constant pressure
    bool th_dry;    // dry potential temperature
    
    // Aerosol specification (for activation)
    std::map<real_t, real_t> dry_distros;  // {mean_radius [m]: concentration [1/kg]}
    real_t kappa;                          // hygroscopicity parameter
    
    // Autoconversion parameterization
    enum acnv_t { Kessler, Beheng, Berry_Reinhardt } acnv_type;
    
    // Accretion parameterization  
    enum accr_t { Khairoutdinov_Kogan, Long } accr_type;
};

Activation Parameterizations

The double-moment scheme can use different aerosol activation schemes:

Available Formulations

  1. Twomey (1959): Simple power-law relationship

    N_c = C * S^k

  2. Abdul-Razzak and Ghan (2000): Physically-based with aerosol size distribution

  3. Khvorostyanov and Curry (2006): Detailed kinetic approach

Configuration:

opts.acti = true;
opts.dry_distros = {
    {0.02e-6, 80e6},   // mode 1: r=20 nm, N=80/cm³
    {0.075e-6, 40e6}   // mode 2: r=75 nm, N=40/cm³
};
opts.kappa = 0.61;      // ammonium sulfate

Autoconversion Formulations

Kessler Type

opts.acnv = true;
opts.acnv_type = opts_t<real_t>::Kessler;

Classic threshold-based autoconversion.

Beheng (1994)

opts.acnv_type = opts_t<real_t>::Beheng;

Stochastic collection equation solution accounting for droplet size distribution.

Berry-Reinhardt

opts.acnv_type = opts_t<real_t>::Berry_Reinhardt;

Considers turbulence effects on collection efficiency.

Common Utilities

Namespace

namespace libcloudphxx::common

Thermodynamic Functions

Available in #include <libcloudph++/common/...>

Moist Air Properties

#include <libcloudph++/common/moist_air.hpp>

// Saturation vapor pressure (Tetens formula)
quantity<si::pressure, real_t> p_vs(quantity<si::temperature, real_t> T);

// Saturation mixing ratio
quantity<si::dimensionless, real_t> r_vs(
    quantity<si::temperature, real_t> T,
    quantity<si::pressure, real_t> p
);

// Virtual temperature
quantity<si::temperature, real_t> T_v(
    quantity<si::temperature, real_t> T,
    quantity<si::dimensionless, real_t> rv
);

Potential Temperature Conversions

#include <libcloudph++/common/theta_std.hpp>
#include <libcloudph++/common/theta_dry.hpp>

// Standard potential temperature
namespace theta_std {
    quantity<si::temperature, real_t> T(
        quantity<si::temperature, real_t> theta,
        quantity<si::pressure, real_t> p
    );
    
    real_t exner(quantity<si::pressure, real_t> p);
}

// Dry potential temperature
namespace theta_dry {
    quantity<si::temperature, real_t> T(
        quantity<si::temperature, real_t> theta_dry,
        quantity<si::mass_density, real_t> rhod
    );
    
    quantity<si::pressure, real_t> p(
        quantity<si::mass_density, real_t> rhod,
        quantity<si::dimensionless, real_t> rv,
        quantity<si::temperature, real_t> T
    );
}

Terminal Velocity

#include <libcloudph++/common/vterm.hpp>

// Terminal fall speed for droplets/particles
quantity<divide_typeof_helper<si::length, si::time>::type, real_t> 
vterm_fall(
    quantity<si::length, real_t> r,          // particle radius
    quantity<si::mass_density, real_t> rho   // air density
);

Physical Constants

#include <libcloudph++/common/const_cp.hpp>

namespace const_cp {
    const quantity<...> c_pd;    // specific heat of dry air
    const quantity<...> c_pv;    // specific heat of water vapor
    const quantity<...> l_v;     // latent heat of vaporization
    const quantity<...> l_f;     // latent heat of fusion
    const quantity<...> R_d;     // gas constant for dry air
    const quantity<...> R_v;     // gas constant for water vapor
}

Output Enumerations

enum class output_t {
    // Mixing ratios
    rv,    // water vapor
    rc,    // cloud water
    rr,    // rain water
    ri,    // ice
    rs,    // snow
    rg,    // graupel
    
    // Number concentrations (for double-moment)
    nc,    // cloud droplet number
    nr,    // rain drop number
    
    // Diagnostics
    th,    // potential temperature
    rl,    // total liquid water
    precip_rate,  // precipitation rate
    
    // Chemistry
    chem_SO2,
    chem_O3,
    chem_H2O2,
    // ... other species
};

Usage Examples

Example 1: Basic Lagrangian Simulation

#include <libcloudph++/lgrngn/factory.hpp>
#include <vector>

int main() {
    // Grid setup
    const int nx = 64, ny = 64, nz = 64;
    const double dx = 100.0, dy = 100.0, dz = 100.0;
    const double dt = 1.0;
    
    // Initialize options
    libcloudphxx::lgrngn::opts_init_t<double> opts_init;
    opts_init.dt = dt;
    opts_init.nx = nx; opts_init.ny = ny; opts_init.nz = nz;
    opts_init.dx = dx; opts_init.dy = dy; opts_init.dz = dz;
    opts_init.sd_conc = 64;
    opts_init.dry_distros = {{0.04e-6, 60e6}};  // 40 nm mode
    opts_init.kappa = 0.61;
    
    // Create particle system
    auto particles = libcloudphxx::lgrngn::factory<double>(
        libcloudphxx::lgrngn::CUDA, opts_init
    );
    
    // Prepare atmospheric state
    std::vector<double> th(nx*ny*nz, 300.0);  // 300 K
    std::vector<double> rv(nx*ny*nz, 0.01);   // 10 g/kg
    std::vector<double> rhod(nx*ny*nz, 1.2);  // 1.2 kg/m³
    
    // Array info
    libcloudphxx::lgrngn::arrinfo_t<double> 
        th_info(th.data(), nullptr),
        rv_info(rv.data(), nullptr),
        rhod_info(rhod.data(), nullptr);
    
    // Initialize
    particles->init(th_info, rv_info, rhod_info);
    
    // Time loop
    libcloudphxx::lgrngn::opts_t<double> opts;
    opts.cond = true;
    opts.coal = true;
    
    for (int step = 0; step < 100; ++step) {
        particles->step_sync(opts, th_info, rv_info, rhod_info);
        
        // Diagnostics
        if (step % 10 == 0) {
            particles->diag_wet_mom(3);  // LWC
            double *lwc = particles->outbuf();
            // Process output...
        }
    }
    
    delete particles;
    return 0;
}

Example 2: Single-Moment Bulk Scheme

#include <libcloudph++/blk_1m/adj_cellwise.hpp>
#include <libcloudph++/blk_1m/rhs_cellwise.hpp>
#include <vector>

int main() {
    const int ncells = 1000;
    const double dt = 1.0;
    
    // Options
    libcloudphxx::blk_1m::opts_t<double> opts;
    opts.cond = true;
    opts.conv = true;
    opts.accr = true;
    opts.th_dry = true;
    opts.const_p = false;
    
    // State variables
    std::vector<double> th(ncells, 300.0);
    std::vector<double> rv(ncells, 0.01);
    std::vector<double> rc(ncells, 0.0);
    std::vector<double> rr(ncells, 0.0);
    std::vector<double> rhod(ncells, 1.2);
    std::vector<double> p(ncells, 100000.0);
    
    // Tendencies
    std::vector<double> dot_rc(ncells, 0.0);
    std::vector<double> dot_rr(ncells, 0.0);
    
    // Time step
    // 1. Calculate microphysics tendencies
    libcloudphxx::blk_1m::rhs_cellwise(opts, dot_rc, dot_rr, rc, rr);
    
    // 2. Apply tendencies
    for (int i = 0; i < ncells; ++i) {
        rc[i] += dot_rc[i] * dt;
        rr[i] += dot_rr[i] * dt;
    }
    
    // 3. Saturation adjustment
    libcloudphxx::blk_1m::adj_cellwise_nwtrph(
        opts, rhod, p, th, rv, rc, dt
    );
    
    return 0;
}

Example 3: Double-Moment Scheme with Activation

#include <libcloudph++/blk_2m/rhs_cellwise.hpp>
#include <vector>

int main() {
    const int ncells = 1000;
    const double dt = 1.0;
    
    // Options with aerosol specification
    libcloudphxx::blk_2m::opts_t<double> opts;
    opts.acti = true;
    opts.acnv = true;
    opts.accr = true;
    opts.const_p = true;
    opts.th_dry = false;
    opts.dry_distros = {{0.04e-6, 100e6}};
    opts.kappa = 0.61;
    
    // State variables
    std::vector<double> th(ncells, 300.0);
    std::vector<double> rv(ncells, 0.01);
    std::vector<double> rc(ncells, 0.0);
    std::vector<double> nc(ncells, 0.0);
    std::vector<double> rr(ncells, 0.0);
    std::vector<double> nr(ncells, 0.0);
    std::vector<double> rhod(ncells, 1.2);
    std::vector<double> p(ncells, 100000.0);
    
    // Tendencies
    std::vector<double> dot_th(ncells, 0.0);
    std::vector<double> dot_rv(ncells, 0.0);
    std::vector<double> dot_rc(ncells, 0.0);
    std::vector<double> dot_nc(ncells, 0.0);
    std::vector<double> dot_rr(ncells, 0.0);
    std::vector<double> dot_nr(ncells, 0.0);
    
    // Calculate tendencies
    libcloudphxx::blk_2m::rhs_cellwise(
        opts,
        dot_th, dot_rv, dot_rc, dot_nc, dot_rr, dot_nr,
        rhod, th, rv, rc, nc, rr, nr, dt, p
    );
    
    // Apply tendencies
    for (int i = 0; i < ncells; ++i) {
        th[i] += dot_th[i] * dt;
        rv[i] += dot_rv[i] * dt;
        rc[i] += dot_rc[i] * dt;
        nc[i] += dot_nc[i] * dt;
        rr[i] += dot_rr[i] * dt;
        nr[i] += dot_nr[i] * dt;
    }
    
    return 0;
}

Compilation and Linking

CMake Integration

find_package(libcloudphxx REQUIRED)

add_executable(myprogram main.cpp)
target_link_libraries(myprogram libcloudphxx::cloudphxx_lgrngn)

Compiler Requirements

  • C++11 or later
  • Boost libraries (Units, OdeInt)
  • Thrust (for GPU backends)
  • CUDA toolkit (for GPU backends)
  • OpenMP (for OpenMP backend)

Compiler Flags

# CPU only
g++ -std=c++11 myprogram.cpp -lcloudphxx_lgrngn -lboost_system

# With CUDA
nvcc -std=c++11 myprogram.cu -lcloudphxx_lgrngn

Thread Safety and Performance

CPU Backends (serial, OpenMP)

  • Thread-safe: No
  • Recommendation: Use separate particle systems for different threads

GPU Backends (CUDA, multi_CUDA)

  • Thread-safe: No
  • Recommendation: One particle system per GPU stream
  • Performance: Optimal with sd_conc ≥ 32 per cell

Best Practices

  1. Minimize host-device transfers: Use sync_in() and step_async() for better GPU utilization
  2. Batch diagnostics: Call multiple diag_* methods before copying results
  3. Choose appropriate sd_conc: Balance accuracy vs. computational cost (typically 32-128)
  4. Use operator splitting: Separate step_cond() from advection/sedimentation for better integration with host model

Error Handling

libcloudph++ uses assertions and exceptions for error handling:

try {
    particles->init(th_info, rv_info, rhod_info);
} catch (const std::exception &e) {
    std::cerr << "Initialization error: " << e.what() << std::endl;
}

Common errors:

  • Invalid array dimensions: Check that nx, ny, nz match array sizes
  • Negative mixing ratios: Ensure rv, rc, rr ≥ 0 before calling functions
  • Pressure consistency: Verify const_p XOR th_dry in options
  • GPU memory: Reduce sd_conc or domain size if out of memory

References

For algorithm details and physics, see:

  1. Shima et al. (2009): Super-droplet method - Q. J. Roy. Meteor. Soc. 135: 1307–1320
  2. Arabas et al. (2015): libcloudph++ 1.0 - Geosci. Model Dev. 8: 1677–1707
  3. Grabowski (1999): Ice microphysics - J. Atmos. Sci. 56: 2429–2454
  4. Kessler (1969): Autoconversion parameterization
  5. Khairoutdinov and Kogan (2000): Double-moment warm rain scheme - Mon. Wea. Rev. 128: 229–243

For user options details, see USER_OPTIONS.md.

Document Version: 1.0
Last Updated: 2024
Library Version: libcloudph++ 2.x