diff --git a/components/eamxx/docs/technical/physics/p3/autoconversion.md b/components/eamxx/docs/technical/physics/p3/autoconversion.md new file mode 100644 index 000000000000..3cf9d2b74228 --- /dev/null +++ b/components/eamxx/docs/technical/physics/p3/autoconversion.md @@ -0,0 +1,298 @@ +# P3 Autoconversion + +This document details the implementation and verification of the cloud water to +rain autoconversion process in the P3 microphysics scheme. + +## Physical Formulation + +The autoconversion parameterization follows the formulation of **Khairoutdinov +and Kogan (2000)**. It represents the collision and coalescence of cloud +droplets to form rain drops. + +### Autoconversion Rate + +The time tendency of rain mass mixing ratio due to autoconversion, +$\left(\frac{\partial q_r}{\partial t}\right)_{\text{auto}}$, is calculated as a +power law function of the cloud water mixing ratio ($q_c$) and the cloud droplet +number concentration ($N_c$). + +$$ +\left(\frac{\partial q_r}{\partial t}\right)_{\text{auto}} = A \, q_c^{\alpha} \, (N_{c, \text{vol}})^{-\beta} +$$ + +Where: + +* $A$ is the autoconversion prefactor ( `autoconversion_prefactor` ). +* $\alpha$ is the cloud water exponent ( `autoconversion_qc_exponent` ). +* $\beta$ is the cloud number exponent ( `autoconversion_nc_exponent` ). + * **Note**: $\beta$ is stored as a positive value (1.79 in default + configuration) following KK2000 convention. In the implementation this + appears as + `pow(nc_incld * sp(1.e-6) * rho, -autoconversion_nc_exponent)`, which is + the code equivalent of the analytic form + $\mathrm{pow}(N_{c, \text{vol}}, -\beta)$. +* $N_{c, \text{vol}}$ is the cloud droplet number concentration in units of + $\text{cm}^{-3}$. + +In the implementation, $N_{c, \text{vol}}$ is derived from the in-cloud number +mixing ratio $N_c$ ($\#/\text{kg}$) and air density $\rho$ ($\text{kg}/\text{m}^3$), +where $N_c$ corresponds to the variable `nc_incld`: + +$$ +N_{c,\text{vol}} = N_c \cdot \rho \cdot 10^{-6} +\;\equiv\; \texttt{nc\_incld} \cdot \texttt{sp(1.e-6)} \cdot \rho +$$ + +\equiv \texttt{nc_incld} \cdot \texttt{sp(1.e-6)} \cdot \rho +### Number Tendencies + +Consistent with the 2-moment approach, the parameterization calculates the +tendencies for both mass and number. + +#### 1. Cloud-to-Rain Number Transfer + +The rate at which cloud droplet number converts to rain droplet number is +determined by preserving the specific number concentration: + +$$ +\left|\frac{\partial N_c}{\partial t}\right|_{\text{auto}} = \left(\frac{\partial q_r}{\partial t}\right)_{\text{auto}} \frac{N_c}{q_c} +$$ + +**Implementation**: The variable `nc2nr_autoconv_tend` stores the positive +magnitude of this transfer rate. The actual tendency of cloud droplet number is +the negative of this value: + +$$ +\frac{\partial N_c}{\partial t} += +-\texttt{nc2nr\_autoconv\_tend} +$$ + +#### 2. Rain Droplet Number Source + +The source of rain droplet number is determined by the characteristic mass of +newly formed rain embryos: + +$$ +\left(\frac{\partial N_r}{\partial t}\right)_{\text{auto}} = \frac{\left(\frac{\partial q_r}{\partial t}\right)_{\text{auto}}}{m_{\text{drop}}(r_{\text{auto}})} +$$ + +where $m_{\text{drop}}(r) = \frac{4}{3}\pi\rho_w r^3$ and $r_{\text{auto}}$ is +the characteristic autoconversion radius (default: 25 $\mu\text{m}$). + +**Implementation**: The variable `ncautr` stores this rain number source. It is +computed as: + +```cpp +ncautr = qc2qr_autoconv_tend * CONS3 +``` + +where `CONS3` is computed locally as: + +$$ +\text{CONS3} = \frac{1}{\frac{4\pi}{3}\rho_w r_{\text{auto}}^3} = \frac{1}{m_{\text{drop}}(r_{\text{auto}})} +$$ + +**Important**: Note that `ncautr` (rain number source) is **not equal** to +`nc2nr_autoconv_tend` (cloud number transfer magnitude). The relationship is: +$$ +\texttt{ncautr} += +\texttt{nc2nr\_autoconv\_tend} +\times \frac{q_c}{N_c} +\times \texttt{CONS3} +$$ + +## Implementation Details + +* **Reference**: This logic is implemented in `cloud_water_autoconversion()` + within [p3_autoconversion_impl.hpp](https://github.com/E3SM-Project/E3SM/blob/master/components/eamxx/src/physics/p3/impl/p3_autoconversion_impl.hpp). +* **Thresholds**: The process is active only when $q_c > 10^{-8} \, \text{kg}/\text{kg}$. +* **Consistency Enforcement**: The implementation ensures that if either the + mass or number tendency is zero, both are set to zero, preventing numerical + artifacts in extreme regimes: + + ```cpp + nc2nr_autoconv_tend.set(qc2qr_autoconv_tend == 0 && context, 0); + qc2qr_autoconv_tend.set(nc2nr_autoconv_tend == 0 && context, 0); + ``` + +* **Subgrid Variance**: The code includes a placeholder for subgrid variance + scaling (Jensen's inequality effect), but it is currently disabled: + + ```cpp + // sgs_var_coef = subgrid_variance_scaling(inv_qc_relvar, sp(2.47) ); + sgs_var_coef = 1; + ``` + + **Note**: When enabled, the scaling factor uses the same exponent + ($\alpha = 2.47$) as the $q_c$ dependency, consistent with Jensen's + inequality for power-law functions: + + $$ + \langle q_c^{\alpha} \rangle > \langle q_c \rangle^{\alpha} + $$ + + when variance exists. + +## Property Tests + +The unit testing suite (`p3_autoconversion_unit_tests.cpp`) employs a "Physics +Property Test" strategy. We validate that the implementation adheres to +fundamental physical principles across a wide parameter space. + +### Tolerance Philosophy + +The test suite uses two categories of tolerances based on what is being +validated: + +#### 1. Identity Tolerances (Precision-Dependent) + +These validate that the implementation follows its documented mathematical +formulas exactly. Since floating-point precision differs between single and +double precision builds, these tolerances are conditionally compiled: + +* **Double precision**: `1e-14` (~100× machine epsilon) +* **Single precision**: `1e-7` (~1000× machine epsilon) + +**Examples**: + +* Verifying `nc2nr_autoconv_tend = qc2qr_autoconv_tend × (Nc/qc)` +* Conservation laws with exact arithmetic + +**Rationale**: These are code correctness checks. Failure indicates +implementation bugs, not physics issues. The tolerance accounts for ~10 +floating-point operations while remaining strict enough to catch formula errors. + +#### 2. Physics Tolerances (Precision-Independent) + +These validate physical approximations and configuration parameters. They do not +depend on floating-point precision: + +* **Standard**: `1%` (0.01) + +**Examples**: + +* Embryo radius consistency (configuration parameter) +* KK2000 parameterization accuracy +* Monotonicity detection thresholds + +**Rationale**: Physics uncertainties exceed numerical precision. A 1% +parameterization error is 1% regardless of whether calculations use float or +double. + +#### 3. Regime Thresholds (Physical Relevance) + +These define when values are physically negligible rather than mathematically +zero: + +* **Small droplet regime floor**: `1e-15` $\text{kg}/\text{kg}/\text{s}$ (physically irrelevant + rate for $r < 1 \, \mu\text{m}$) +* **Absolute floor**: `1e-30` (proxy for numerical zero) + +**Rationale**: Power-law formulas produce non-zero values even in physically +irrelevant regimes. The small droplet threshold distinguishes when mean droplet +sizes are too small for efficient collision-coalescence ($r < 1 \, \mu\text{m}$, +representing incipient cloud or sub-threshold conditions) from regimes where +autoconversion is physically active. + +### Test Strategy + +We perform a dense sampling of the phase space: + +* **Cloud Water ($q_c$)**: Logarithmic sweep from $5 \times 10^{-9}$ + (sub-threshold regime) to $10^{-2} \, \text{kg}/\text{kg}$ (heavy precipitation). + * Below $10^{-8}$: autoconversion inactive (incipient cloud) + * Above $10^{-8}$: active cloud-to-rain transition +* **Cloud Number ($N_c$)**: Logarithmic sweep from $10^{6}$ to $10^{9} \, \#/\text{m}^3$. + +### Parameter Validation + +Before running property checks, the test suite verifies that the runtime +configuration parameters match the expected Khairoutdinov and Kogan (2000) +constants (e.g., $A=1350$, $\alpha=2.47$, $r_{\text{auto}}=25 \, \mu\text{m}$). + +### 1. Threshold Behavior + +For $q_c < 10^{-8} \, \text{kg}/\text{kg}$, the test verifies that the autoconversion +rate is strictly zero. + +### 2. Monotonicity & Sensitivity + +We verify the sign of the partial derivatives with respect to the input state +variables: + +* **Colloidal Stability (Inverse $N_c$ Dependency)**: + Increasing droplet number concentration while holding water content fixed must + strictly decrease the rate. The test enforces this with **no relative + tolerance**, as the physics dictates a strict inverse relationship. + + $$ + \frac{\partial}{\partial N_c} \left(\frac{\partial q_r}{\partial t}\right)_{\text{auto}} < 0 + $$ + +* **Water Content Sensitivity (Positive $q_c$ Dependency)**: + Increasing water content while holding number concentration fixed should + increase the rate. + + $$ + \frac{\partial}{\partial q_c} \left(\frac{\partial q_r}{\partial t}\right)_{\text{auto}} > 0 + $$ + +### 3. Consistency Constraints + +We check that the derived number tendencies match their physical definitions. + +* **Specific Loss Conservation** (Mathematical Identity): + The code checks that the number loss rate satisfies the conservation + relationship: $(dNc/dt)_{\text{auto}} = (dqr/dt)_{\text{auto}} \cdot (N_c / q_c)$. + + **Tolerance**: Precision-dependent (`1e-14` for double, `1e-7` for float). + This is an exact mathematical identity in the implementation, so failures + indicate code bugs rather than physics issues. + +* **Rain Embryo Mass** (Configuration Consistency): + The implicit mass of the newly formed rain drops is recovered from the ratio + of mass tendency to number tendency: + + $$ + m_{\text{effective}} = \frac{\left(\partial q_r / \partial t\right)}{\left(\partial N_r / \partial t\right)} + $$ + + The test verifies that this mass corresponds to the characteristic radius + $r_{\text{auto}}$ (25 $\mu\text{m}$) within **1% tolerance** (physics-based, + precision-independent). This confirms the `autoconversion_radius` + configuration parameter is correctly preserved through the calculation. + +### 4. Physical Limits + +* **Small Droplet Regime Limit** (Physical Relevance Check): + When the mean droplet radius is very small ($r < 1 \, \mu\text{m}$, characteristic of + incipient cloud or sub-threshold conditions), autoconversion must be + physically negligible due to low collision efficiency. **The test uses a + threshold of $10^{-15} \, \text{kg}/\text{kg}/\text{s}$** to distinguish "physically irrelevant" + from "mathematically zero." This accounts for the power-law formula producing + non-zero values in regimes where collision-coalescence has no physical + significance. This is a regime threshold, not a precision-dependent tolerance. + +### 5. Variance Scaling (Disabled) + +The test suite checks for subgrid variance scaling. + +* **Principle**: Due to the nonlinearity of the autoconversion rate, subgrid + variability in $q_c$ should enhance the grid-mean rate. +* **Status**: Since the feature is currently disabled in the implementation + (`sgs_var_coef = 1`), the test detects this state and reports it as + "DISABLED" without failing. If the rates differ significantly but do not show + enhancement, it flags a failure. + +### Tolerance Summary + +| Test Category | Tolerance Type | Value (Double/Single) | Rationale | +| ------------- | -------------- | --------------------- | --------- | +| **Specific Loss Conservation** | Identity | `1e-14` / `1e-7` | Mathematical identity (2 FLOPs) | +| **Rain Embryo Mass** | Physics | `1%` (both) | Configuration parameter accuracy | +| **Monotonicity (Nc)** | Strict | No tolerance | Physics requires strict inequality | +| **Monotonicity (qc)** | Physics | `0.1%` | Physical sensitivity detection | +| **Small Droplet Regime** | Physical | `1e-15` $\text{kg}/\text{kg}/\text{s}$ | Relevance threshold ($r < 1 \, \mu\text{m}$) | +| **Variance Detection** | Feature | `0.1%` | Distinguish on/off state | diff --git a/components/eamxx/mkdocs.yml b/components/eamxx/mkdocs.yml index d8edc87f3987..68e38db5e956 100644 --- a/components/eamxx/mkdocs.yml +++ b/components/eamxx/mkdocs.yml @@ -54,6 +54,9 @@ nav: - 'Wrap Torch emulators': 'developer/torch_support.md' - 'Technical Guide': - 'Overview': 'technical/index.md' + - 'Physics': + - 'P3': + - 'Autoconversion': 'technical/physics/p3/autoconversion.md' - 'AeroCom cloud top': 'technical/aerocom_cldtop.md' edit_uri: "" diff --git a/components/eamxx/src/physics/p3/tests/p3_autoconversion_unit_tests.cpp b/components/eamxx/src/physics/p3/tests/p3_autoconversion_unit_tests.cpp index ffa0869ab7d7..53fdd780a4dd 100644 --- a/components/eamxx/src/physics/p3/tests/p3_autoconversion_unit_tests.cpp +++ b/components/eamxx/src/physics/p3/tests/p3_autoconversion_unit_tests.cpp @@ -10,6 +10,7 @@ #include #include #include +#include #include // std::setprecision namespace scream { @@ -123,33 +124,409 @@ void cloud_water_autoconversion_unit_bfb_tests() { cloud_water_autoconversion_unit_bfb_tests(); } - KOKKOS_FUNCTION static void autoconversion_is_positive(const Int &i, Int &errors){ - - const Spack rho(1.0), inv_qc_relvar(1.0); - Spack qc_incld, nc_incld(1e7), qc2qr_autoconv_tend(0.0), nc2nr_autoconv_tend(0.0), ncautr(0.0); - for(int si=0; si::P3Runtime(); + + const Real expected_prefactor = 1350.0; // s^-1 + const Real expected_qc_exp = 2.47; // dimensionless + const Real expected_nc_exp = 1.79; // dimensionless + const Real expected_radius = 25.0e-6; // m (25 um) + + const Real tolerance = 0.01; // 1% tolerance + + // Validate prefactor + Real prefactor_error = std::abs(runtime.autoconversion_prefactor - expected_prefactor) + / expected_prefactor; + if (prefactor_error > tolerance) { + std::cout << "WARNING: Prefactor mismatch\n"; + std::cout << " Expected: " << expected_prefactor << " s^-1\n"; + std::cout << " Actual: " << runtime.autoconversion_prefactor << " s^-1\n"; + } else { + std::cout << " Prefactor: " << runtime.autoconversion_prefactor << " s^-1\n"; } - Functions::cloud_water_autoconversion( - rho, qc_incld, nc_incld, inv_qc_relvar, qc2qr_autoconv_tend, - nc2nr_autoconv_tend, ncautr, - p3::Functions::P3Runtime()); - if((qc2qr_autoconv_tend < 0.0).any()) { - errors++; + + // Validate qc exponent + Real qc_exp_error = std::abs(runtime.autoconversion_qc_exponent - expected_qc_exp) + / expected_qc_exp; + if (qc_exp_error > tolerance) { + std::cout << "WARNING: qc exponent mismatch\n"; + std::cout << " Expected: " << expected_qc_exp << "\n"; + std::cout << " Actual: " << runtime.autoconversion_qc_exponent << "\n"; + } else { + std::cout << " qc exponent: " << runtime.autoconversion_qc_exponent << "\n"; + } + + // Validate Nc exponent + Real nc_exp_error = std::abs(runtime.autoconversion_nc_exponent - expected_nc_exp) + / expected_nc_exp; + if (nc_exp_error > tolerance) { + std::cout << "WARNING: Nc exponent mismatch\n"; + std::cout << " Expected: " << expected_nc_exp << "\n"; + std::cout << " Actual: " << runtime.autoconversion_nc_exponent << "\n"; + } else { + std::cout << " Nc exponent: " << runtime.autoconversion_nc_exponent << "\n"; + } + + // Validate characteristic radius + Real radius_error = std::abs(runtime.autoconversion_radius - expected_radius) + / expected_radius; + if (radius_error > tolerance) { + std::cout << "INFO: Characteristic radius differs from documented value\n"; + std::cout << " Documentation: " << expected_radius*1e6 << " um\n"; + std::cout << " Implementation: " << runtime.autoconversion_radius*1e6 << " um\n"; + } else { + std::cout << " Characteristic radius: " << runtime.autoconversion_radius*1e6 << " um\n"; } + + std::cout << "===========================\n\n"; } - void run_physics(){ + void run_physics() { + validate_autoconversion_parameters(); + + // ========================================================================= + // Test Tolerance Definitions + // ========================================================================= + // Identity tolerances: precision-dependent (mathematical exactness) + // These verify that the implementation follows its documented formulas. +#ifdef SCREAM_DOUBLE_PRECISION + const Real identity_tol = 1e-14; // ~100x machine epsilon (double) +#else + const Real identity_tol = 1e-7; // ~1000x machine epsilon (float) +#endif + + // Physics tolerances: precision-independent (empirical/configuration) + // These verify physical approximations and configuration parameters. + const Real physics_tol = 0.01; // 1% (matches KK2000 accuracy, config params) + + // Regime thresholds: physical relevance cutoffs (DO NOT MODIFY) + const Real small_droplet_regime_floor = 1e-15; // Physically negligible rate for r < 1μm + const Real absolute_floor = 1e-30; // Numerical zero proxy + const Real detect_threshold = 1e-3; // Feature detection (0.1%) + + // ========================================================================= + // P3 Autoconversion Physics Property Testing Plan + // ========================================================================= + // This suite validates that the `cloud_water_autoconversion` implementation + // adheres to the fundamental physical principles of the Khairoutdinov and + // Kogan (2000) parameterization. + // + // Strategy: + // Uses a deterministic grid-sampling strategy that sweeps the parameter + // space logarithmically to cover regimes ranging from sub-threshold + // (qc < 1e-8) through incipient cloud to heavy precipitation, and + // pristine to polluted conditions. + // ========================================================================= + + // Grid Setup: qc from 5e-9 to 1e-2 kg/kg + // Lower bound: sub-threshold regime (qc < 1e-8, autoconversion inactive) + // Upper bound: heavy precipitation (qc ~ 1e-2, vigorous autoconversion) + const int n_qc = 40; + const int n_nc = 40; + const int num_cases = n_qc * n_nc; + + const Real log_qc_min = std::log10(5.0e-9); + const Real log_qc_max = std::log10(1.0e-2); + const Real log_nc_min = std::log10(1.0e6); + const Real log_nc_max = std::log10(1.0e9); + + view_1d qc_dev("qc_dev", num_cases); + view_1d nc_dev("nc_dev", num_cases); + view_1d rho_dev("rho_dev", num_cases); + view_1d inv_qc_relvar_dev("inv_qc_relvar_dev", num_cases); + + view_1d qc2qr_dev("qc2qr_dev", num_cases); + view_1d nc2nr_dev("nc2nr_dev", num_cases); + view_1d ncautr_dev("ncautr_dev", num_cases); + + auto qc_host = Kokkos::create_mirror_view(qc_dev); + auto nc_host = Kokkos::create_mirror_view(nc_dev); + auto rho_host = Kokkos::create_mirror_view(rho_dev); + auto inv_qc_relvar_host = Kokkos::create_mirror_view(inv_qc_relvar_dev); + + for (int j = 0; j < n_nc; ++j) { + for (int i = 0; i < n_qc; ++i) { + int idx = j * n_qc + i; + Real alpha_qc = (Real)i / (Real)(n_qc - 1); + Real alpha_nc = (Real)j / (Real)(n_nc - 1); + + Real qc = std::pow(10.0, log_qc_min + alpha_qc * (log_qc_max - log_qc_min)); + Real nc = std::pow(10.0, log_nc_min + alpha_nc * (log_nc_max - log_nc_min)); + + qc_host(idx) = qc; + nc_host(idx) = nc; + rho_host(idx) = 1.0; + inv_qc_relvar_host(idx) = 1.0; + } + } - int nerr = 0; + Kokkos::deep_copy(qc_dev, qc_host); + Kokkos::deep_copy(nc_dev, nc_host); + Kokkos::deep_copy(rho_dev, rho_host); + Kokkos::deep_copy(inv_qc_relvar_dev, inv_qc_relvar_host); + + const int pack_size = Spack::n; + const int num_packs = (num_cases + pack_size - 1) / pack_size; + + Kokkos::parallel_for("P3_Property_Test_Kernel", num_packs, KOKKOS_LAMBDA(const Int& ip) { + const int offset = ip * pack_size; + + Spack rho, qc, nc, inv_qc_relvar; + Spack qc2qr, nc2nr, ncautr; + + bool any_valid = false; + for (int s = 0; s < pack_size; ++s) { + int idx = offset + s; + if (idx < num_cases) { + rho[s] = rho_dev(idx); + qc[s] = qc_dev(idx); + nc[s] = nc_dev(idx); + inv_qc_relvar[s] = inv_qc_relvar_dev(idx); + any_valid = true; + } else { + rho[s] = 1.0; + qc[s] = 0.0; + nc[s] = 1.0e6; + inv_qc_relvar[s] = 1.0; + } + } + + if (any_valid) { + Functions::cloud_water_autoconversion( + rho, qc, nc, inv_qc_relvar, + qc2qr, nc2nr, ncautr, + p3::Functions::P3Runtime() + ); + } + + for (int s = 0; s < pack_size; ++s) { + int idx = offset + s; + if (idx < num_cases) { + qc2qr_dev(idx) = qc2qr[s]; + nc2nr_dev(idx) = nc2nr[s]; + ncautr_dev(idx) = ncautr[s]; + } + } + }); - Kokkos::parallel_reduce("TestAutoConversionPositive", 1000, KOKKOS_LAMBDA(const Int& i, Int& errors) { - autoconversion_is_positive(i, errors); - }, nerr); + Kokkos::fence(); + + auto qc2qr_host = Kokkos::create_mirror_view(qc2qr_dev); + auto nc2nr_host = Kokkos::create_mirror_view(nc2nr_dev); + auto ncautr_host = Kokkos::create_mirror_view(ncautr_dev); + + Kokkos::deep_copy(qc2qr_host, qc2qr_dev); + Kokkos::deep_copy(nc2nr_host, nc2nr_dev); + Kokkos::deep_copy(ncautr_host, ncautr_dev); + + int failures = 0; + + for (int j = 0; j < n_nc; ++j) { + for (int i = 0; i < n_qc; ++i) { + int idx = j * n_qc + i; + Real qc = qc_host(idx); + Real nc = nc_host(idx); + Real R = qc2qr_host(idx); + Real N_loss = nc2nr_host(idx); + Real R_ncautr = ncautr_host(idx); + + // Check for non-zero below threshold + if (qc < 1e-8) { + if (R > 1e-30) { + std::cout << "FAIL: Non-zero rate below threshold (1e-8). qc=" << qc << " R=" << R << "\n"; + failures++; + } + continue; + } + + if (R < 1e-30) R = 0.0; + if (N_loss < 1e-30) N_loss = 0.0; + if (R_ncautr < 1e-30) R_ncautr = 0.0; + + // ===================================================================== + // 1. Physical Sensitivity Tests (Monotonicity) + // ===================================================================== + // Uses strict inequality for Nc (physics requirement) and 0.1% tolerance for qc + // (sensitivity detection). Both reference the global absolute_floor (1e-30). + + const Real relative_tolerance = 1e-3; // 0.1% sensitivity threshold for qc monotonicity + + if (j < n_nc - 1) { + int idx_next = (j + 1) * n_qc + i; + Real R_next = qc2qr_host(idx_next); + if (R_next < 1e-30) R_next = 0.0; + + // A. Colloidal Stability (Inverse Nc Dependency): dR/dNc < 0 + // R_next (high Nc) must be strictly less than R (low Nc). + if (R > absolute_floor) { + // FAIL if R_next is greater than or equal to R. + // We do not use the loose 'relative_tolerance' here because the physics + // requires strict monotonicity. + if (R_next >= R) { + std::cout << "Monotonicity Nc Fail: R increased or stayed same with Nc. qc=" << qc + << " R=" << R << " R_next=" << R_next << "\n"; + failures++; + } + } + } + + if (i < n_qc - 1) { + int idx_next = j * n_qc + (i + 1); + Real R_next = qc2qr_host(idx_next); + if (R_next < 1e-30) R_next = 0.0; + + // B. Water Content Sensitivity (Positive qc Dependency): dR/dqc > 0 + // R_next (high qc) should be greater than R (low qc). + if (R_next > absolute_floor) { + // If R_next < R, it's a failure. + Real min_allowed_R_next = R * (1.0 - relative_tolerance); + if (R_next < min_allowed_R_next) { + std::cout << "Monotonicity qc Fail: R decreased with qc. qc=" << qc + << " R=" << R << " R_next=" << R_next << "\n"; + failures++; + } + } + } + + // ===================================================================== + // 2. Consistency Tests (Constraint Satisfaction) + // ===================================================================== + + if (R > absolute_floor) { + // A. Specific Loss Conservation (Mathematical Identity) + // Verifies: nc2nr_autoconv_tend = qc2qr_autoconv_tend * (nc/qc) + Real expected_N_loss = R * nc / qc; + Real rel_error = std::abs(N_loss - expected_N_loss) / std::max(absolute_floor, expected_N_loss); + + if (rel_error > identity_tol) { + std::cout << "Specific Loss Conservation Fail (should be near machine precision):\n"; + std::cout << " Expected: " << expected_N_loss << "\n"; + std::cout << " Actual: " << N_loss << "\n"; + std::cout << " Rel Error: " << rel_error << " (tolerance: " << identity_tol << ")\n"; + failures++; + } + } + + // B. Rain Embryo Size Consistency + if (R > 1e-20 && R_ncautr > 1e-20) { + // B. Rain Embryo Size Consistency (Configuration Parameter) + // Verifies: autoconversion_radius parameter is preserved through calculation + Real mass_drop = R / R_ncautr; + + using C = scream::physics::Constants; + const Real r_expected = 25.0e-6; // Must match autoconversion_radius + const Real expected_mass = (4.0/3.0) * C::Pi * 1000.0 * std::pow(r_expected, 3); + + if (std::abs(mass_drop - expected_mass) / expected_mass > physics_tol) { + std::cout << "Embryo Size Fail (configuration mismatch):\n"; + std::cout << " Computed mass: " << mass_drop << " kg\n"; + std::cout << " Expected mass: " << expected_mass << " kg (r=25um)\n"; + std::cout << " Tolerance: " << physics_tol*100 << "%\n"; + failures++; + } + } + + // ===================================================================== + // 3. Limit & Singularity Tests + // ===================================================================== + + Real mean_mass = qc / nc; + // Correct formula: r = (3*m / (4*pi*rho))^(1/3) + // rho_w = 1000. + using C = scream::physics::Constants; + Real mean_rad = std::pow((3.0 * mean_mass) / (4.0 * C::Pi * 1000.0), 1.0/3.0); + + // A. Small Droplet Regime Limit + if (mean_rad < 1e-6) { + // A. Small Droplet Regime Limit (Physical Relevance Check) + // For very small activated droplets (r < 1 μm), autoconversion rates + // should be physically negligible due to low collision efficiency. + // This regime represents incipient cloud or sub-threshold conditions. + // Note: Implementation threshold is qc < 1e-8, not radius-based. + if (R > small_droplet_regime_floor) { + std::cout << "Small Droplet Regime Fail: r=" << mean_rad << " m, R=" << R << " kg/kg/s\n"; + std::cout << " (Rate exceeds physical relevance threshold: " << small_droplet_regime_floor << ")\n"; + failures++; + } + } + } + } + // ========================================================================= + // CHECK 4: Subgrid Variance Scaling + // ========================================================================= + + view_1d var_rate_1("rate_1", n_qc); + view_1d var_rate_2("rate_2", n_qc); + + Real fixed_nc = 100.0e6; + + Kokkos::parallel_for("Variance_Check", n_qc, KOKKOS_LAMBDA(const int& i) { + Real alpha = (Real)i / (Real)(n_qc - 1); + Real qc = std::pow(10.0, -6.0 + alpha * 4.0); + + Spack s_rho(1.0), s_qc(qc), s_nc(fixed_nc); + Spack s_rate_1, s_rate_2, dummy; + + // Run 1: Homogeneous + Functions::cloud_water_autoconversion( + s_rho, s_qc, s_nc, Spack(1.0), + s_rate_1, dummy, dummy, + p3::Functions::P3Runtime() + ); + + // Run 2: High Variance + Functions::cloud_water_autoconversion( + s_rho, s_qc, s_nc, Spack(0.1), + s_rate_2, dummy, dummy, + p3::Functions::P3Runtime() + ); + + var_rate_1(i) = s_rate_1[0]; + var_rate_2(i) = s_rate_2[0]; + }); + Kokkos::fence(); - REQUIRE(nerr == 0); + auto h_rate_1 = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), var_rate_1); + auto h_rate_2 = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), var_rate_2); + + bool is_variance_scaling_active = false; + int variance_suppression_failures = 0; + const Real variance_detect_eps = detect_threshold; // 0.1% change indicates active scaling + const Real variance_validate_eps = 1e-12; // Numerical noise floor for suppression check + + for(int i=0; i 1e-30) { + Real ratio = h_rate_2(i) / h_rate_1(i); + + // Detection: Do standard and variance rates differ? + if (std::abs(ratio - 1.0) > variance_detect_eps) { + is_variance_scaling_active = true; + + // Validation: If active, is the rate suppressed? (Physical violation) + // Jensen's inequality for convex functions implies enhancement (ratio > 1). + // We fail if ratio < 1.0 (minus tolerance). + if (ratio < 1.0 - variance_validate_eps) { + variance_suppression_failures++; + } + } + } + } + + if (is_variance_scaling_active) { + if (variance_suppression_failures > 0) { + std::cout << "Variance Scaling FAIL: " << variance_suppression_failures << " cases suppressed rate (physical violation).\n"; + failures += variance_suppression_failures; + } else { + std::cout << "Variance Scaling: ACTIVE and working correctly (rates enhanced).\n"; + } + } else { + std::cout << "Variance Scaling: DISABLED (sgs_var_coef approx 1)\n"; + } + REQUIRE(failures == 0); } }; // TestP3CloudWaterAutoconversion