From c116a74b952bca83611374046bf523a09e9923e2 Mon Sep 17 00:00:00 2001 From: "Sebastian D. Eastham" Date: Tue, 18 Aug 2026 19:35:23 -0400 Subject: [PATCH 01/14] Implemented semi-lagrangian advection solver in place of diffusive multi-step Eulerian --- .../include/FVM_ANDS/AdvDiffSystem.hpp | 4 + Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp | 156 ++++++++++++++++++ Code.v05-00/src/FVM_ANDS/FVM_Solver.cpp | 48 +----- Code.v05-00/tests/test_adv_diff_solver.cpp | 67 ++++++++ 4 files changed, 233 insertions(+), 42 deletions(-) diff --git a/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp b/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp index 136ff7335..6c79dd3b9 100644 --- a/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp +++ b/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp @@ -11,6 +11,9 @@ namespace FVM_ANDS{ // Separate the SOR solver for testing without having to build an AdvDiffSystem object void sor_solve(const Eigen::SparseMatrix &A, const Eigen::VectorXd &rhs, Eigen::VectorXd &phi, double omega = 1.0, double threshold = 1e-3, int n_iters = 3); + // 1D Semi-Lagrangian advection helper + void semiLagrangianAdvection1D(std::vector& slice, double velocity, double dt, double ds, double bc_left, double bc_right); + struct AdvDiffParams { AdvDiffParams(double u, double v, double shear, double Dh, double Dv, double dt){ this->u = u; @@ -35,6 +38,7 @@ namespace FVM_ANDS{ const Eigen::VectorXd& calcRHS(); void applyBoundaryCondition(); void updateBoundaryCondition(const BoundaryConditions& bc); + void semiLagrangianAdvection(double dt); Eigen::VectorXd forwardEulerAdvection(bool operatorSplit = false, bool parallelAdvection = false) const noexcept; // Breakup the implementation of sor_solve to allow for easy testing by inputing an arbitrary linear system to solve: // Implementation is moved outside of the class, and make class method to be used in code diff --git a/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp b/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp index d4b6c823d..7e5cae952 100644 --- a/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp +++ b/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp @@ -541,6 +541,162 @@ namespace FVM_ANDS{ applyBoundaryCondition(); //need this to calculate minmod function at some timestep. } + void semiLagrangianAdvection1D( + std::vector& slice, + double velocity, + double dt, + double ds, + double bc_left, + double bc_right) + { + const int N = static_cast(slice.size()); + if (N == 0) return; + + double disp = velocity * dt; + if (std::abs(disp) < 1.0e-14 || dt <= 0.0) return; + + if (velocity > 0.0) { + int k = static_cast(disp / ds); + double rem_disp = disp - k * ds; + double rem_dt = rem_disp / velocity; + + // Step 1: Discrete shift to the right (positive direction) + if (k > 0) { + if (k >= N) { + std::fill(slice.begin(), slice.end(), bc_left); + } else { + for (int m = N - 1; m >= k; --m) { + slice[m] = slice[m - k]; + } + for (int m = 0; m < k; ++m) { + slice[m] = bc_left; + } + } + } + + // Step 2: Fractional Forward Euler step + if (rem_disp > 1.0e-12 && rem_dt > 0.0) { + std::vector face_flux(N + 1, 0.0); + face_flux[0] = bc_left; + + auto minmod = [](double a, double b) -> double { + if (a * b <= 0.0) return 0.0; + return (a > 0.0) ? std::min(a, b) : std::max(a, b); + }; + + for (int m = 0; m < N - 1; ++m) { + double diff_up = (m == 0) ? (2.0 * (slice[0] - bc_left)) : (slice[m] - slice[m - 1]); + double diff_down = slice[m + 1] - slice[m]; + double slope = minmod(diff_up, diff_down); + face_flux[m + 1] = slice[m] + 0.5 * slope; + } + + // Outflow face at m = N + double diff_up_last = (N >= 2) ? (slice[N - 1] - slice[N - 2]) : (2.0 * (slice[0] - bc_left)); + double diff_down_last = 2.0 * (bc_right - slice[N - 1]); + double slope_last = minmod(diff_up_last, diff_down_last); + face_flux[N] = slice[N - 1] + 0.5 * slope_last; + + double cfl_frac = velocity * rem_dt / ds; + for (int m = 0; m < N; ++m) { + slice[m] -= cfl_frac * (face_flux[m + 1] - face_flux[m]); + } + } + } else { // velocity < 0.0 + double abs_vel = -velocity; + double abs_disp = -disp; + int k = static_cast(abs_disp / ds); + double rem_disp = abs_disp - k * ds; + double rem_dt = rem_disp / abs_vel; + + // Step 1: Discrete shift to the left (negative direction) + if (k > 0) { + if (k >= N) { + std::fill(slice.begin(), slice.end(), bc_right); + } else { + for (int m = 0; m < N - k; ++m) { + slice[m] = slice[m + k]; + } + for (int m = N - k; m < N; ++m) { + slice[m] = bc_right; + } + } + } + + // Step 2: Fractional Forward Euler step + if (rem_disp > 1.0e-12 && rem_dt > 0.0) { + std::vector face_flux(N + 1, 0.0); + face_flux[N] = bc_right; + + auto minmod = [](double a, double b) -> double { + if (a * b <= 0.0) return 0.0; + return (a > 0.0) ? std::min(a, b) : std::max(a, b); + }; + + for (int m = 0; m < N - 1; ++m) { + double diff_up = (m + 1 == N - 1) ? (2.0 * (bc_right - slice[N - 1])) : (slice[m + 2] - slice[m + 1]); + double diff_down = slice[m + 1] - slice[m]; + double slope = minmod(diff_down, diff_up); + face_flux[m + 1] = slice[m + 1] - 0.5 * slope; + } + + // Outflow face at m = 0 + double diff_up_0 = 2.0 * (slice[0] - bc_left); + double diff_down_0 = (N >= 2) ? (slice[1] - slice[0]) : diff_up_0; + double slope_0 = minmod(diff_down_0, diff_up_0); + face_flux[0] = slice[0] - 0.5 * slope_0; + + double cfl_frac = abs_vel * rem_dt / ds; + for (int m = 0; m < N; ++m) { + slice[m] -= cfl_frac * (face_flux[m] - face_flux[m + 1]); + } + } + } + } + + void AdvDiffSystem::semiLagrangianAdvection(double dt) { + // 1. Horizontal Advection along X (row by row) + #pragma omp parallel for default(shared) schedule(static) + for (int j = 0; j < ny_; ++j) { + double u_j = u_double_ - yCoord_[j] * shear_; + if (std::abs(u_j) > 1.0e-14) { + std::vector row(nx_); + for (int i = 0; i < nx_; ++i) { + int idx = twoDIdx_to_vecIdx(i, j, nx_, ny_, format_); + row[i] = phi_[idx]; + } + double bc_left = bcVals_left_.empty() ? 0.0 : bcVals_left_[j]; + double bc_right = bcVals_right_.empty() ? 0.0 : bcVals_right_[j]; + semiLagrangianAdvection1D(row, u_j, dt, dx_, bc_left, bc_right); + for (int i = 0; i < nx_; ++i) { + int idx = twoDIdx_to_vecIdx(i, j, nx_, ny_, format_); + phi_[idx] = row[i]; + } + } + } + + // 2. Vertical Advection along Y (column by column) + if (std::abs(v_double_) > 1.0e-14) { + #pragma omp parallel for default(shared) schedule(static) + for (int i = 0; i < nx_; ++i) { + std::vector col(ny_); + for (int j = 0; j < ny_; ++j) { + int idx = twoDIdx_to_vecIdx(i, j, nx_, ny_, format_); + col[j] = phi_[idx]; + } + double bc_bot = bcVals_bot_.empty() ? 0.0 : bcVals_bot_[i]; + double bc_top = bcVals_top_.empty() ? 0.0 : bcVals_top_[i]; + semiLagrangianAdvection1D(col, v_double_, dt, dy_, bc_bot, bc_top); + for (int j = 0; j < ny_; ++j) { + int idx = twoDIdx_to_vecIdx(i, j, nx_, ny_, format_); + phi_[idx] = col[j]; + } + } + } + + applyBoundaryCondition(); + } + Eigen::VectorXd AdvDiffSystem::forwardEulerAdvection(bool operatorSplit, bool parallelAdvection) const noexcept{ Eigen::VectorXd soln(nTotalPoints_); // double avgBackgroundCalcTime = 0; diff --git a/Code.v05-00/src/FVM_ANDS/FVM_Solver.cpp b/Code.v05-00/src/FVM_ANDS/FVM_Solver.cpp index 719e4d1fb..dd906b060 100644 --- a/Code.v05-00/src/FVM_ANDS/FVM_Solver.cpp +++ b/Code.v05-00/src/FVM_ANDS/FVM_Solver.cpp @@ -66,27 +66,15 @@ namespace FVM_ANDS{ #endif //Strang Splitting - //Step 1: Calculate explicit advection timestep based on CFL condition set - bool operatorSplit = true; - double courant = advDiffSys_.courant(); double dt_max = advDiffSys_.timestep(); - double dt_adv = dt_max * (courant_max / courant); - - int n_timesteps_advection_half = std::ceil((0.5 * dt_max) / dt_adv); - dt_adv = (0.5 * dt_max) / n_timesteps_advection_half; #ifdef ENABLE_TIMING - std::cout << " N Advection timesteps = 2 * " << n_timesteps_advection_half << std::endl; auto start = std::chrono::high_resolution_clock::now(); #endif - //Step 2: Solve Advection for half timestep - advDiffSys_.updateTimestep(dt_adv); - for(int i = 0; i < n_timesteps_advection_half; i++){ - advDiffSys_.updatePhi(advDiffSys_.forwardEulerAdvection(operatorSplit, parallelAdvection)); - advDiffSys_.applyBoundaryCondition(); - } + //Step 1: Solve Advection for half timestep via Semi-Lagrangian advection + advDiffSys_.semiLagrangianAdvection(0.5 * dt_max); #ifdef ENABLE_TIMING auto stop = std::chrono::high_resolution_clock::now(); @@ -96,7 +84,7 @@ namespace FVM_ANDS{ start = std::chrono::high_resolution_clock::now(); #endif - //Step 3: Implicitly solve diffusion (first to help smoothen out potential steep gradients) + //Step 2: Implicitly solve diffusion (first to help smoothen out potential steep gradients) advDiffSys_.updateTimestep(dt_max); // Should never happen given using operatorSplit is hard coded into runTransport and above, but serves as @@ -128,12 +116,6 @@ namespace FVM_ANDS{ start = std::chrono::high_resolution_clock::now(); #endif - // auto mat = advDiffSys_.getCoefMatrix(); - // auto b = advDiffSys_.getRHS(); - // solver_.compute(mat); - // Eigen::VectorXd solution = solver_.solveWithGuess(b, advDiffSys_.phi()); - // advDiffSys_.updatePhi(std::move(solution)); - advDiffSys_.sor_solve(); #ifdef ENABLE_TIMING @@ -144,15 +126,8 @@ namespace FVM_ANDS{ start = std::chrono::high_resolution_clock::now(); #endif - //Step 4: Explicitly solve advection to full timestep - - advDiffSys_.updateTimestep(dt_adv); - for(int i = 0; i < n_timesteps_advection_half; i++){ - advDiffSys_.updatePhi(advDiffSys_.forwardEulerAdvection(operatorSplit)); - advDiffSys_.applyBoundaryCondition(); - } - - advDiffSys_.updateTimestep(dt_max); + //Step 3: Solve advection for second half timestep via Semi-Lagrangian advection + advDiffSys_.semiLagrangianAdvection(0.5 * dt_max); #ifdef ENABLE_TIMING stop = std::chrono::high_resolution_clock::now(); @@ -185,19 +160,8 @@ namespace FVM_ANDS{ advDiffSys_.updatePhi(vec_Eigen); advDiffSys_.updateBoundaryCondition(bc); - bool operatorSplit = true; - double courant = advDiffSys_.courant(); double dt_max = advDiffSys_.timestep(); - double dt_adv = dt_max * (courant_max / courant); - - int n_timesteps_advection_half = std::ceil((0.5 * dt_max) / dt_adv); - dt_adv = (0.5 * dt_max) / n_timesteps_advection_half; - - advDiffSys_.updateTimestep(dt_adv); - for(int i = 0; i < n_timesteps_advection_half; i++){ - advDiffSys_.updatePhi(advDiffSys_.forwardEulerAdvection(operatorSplit)); - advDiffSys_.applyBoundaryCondition(); - } + advDiffSys_.semiLagrangianAdvection(0.5 * dt_max); vec = eigenVec_to_std2dVec(advDiffSys_.phi(), vec[0].size(), vec.size()); } diff --git a/Code.v05-00/tests/test_adv_diff_solver.cpp b/Code.v05-00/tests/test_adv_diff_solver.cpp index cde4c8313..313e82ecf 100644 --- a/Code.v05-00/tests/test_adv_diff_solver.cpp +++ b/Code.v05-00/tests/test_adv_diff_solver.cpp @@ -505,4 +505,71 @@ namespace FVM_ANDS{ REQUIRE(std::abs(maxy-0.381) < 0.01); } + + TEST_CASE("Semi-Lagrangian 1D Advection", "[advection]"){ + SECTION("Pure integer shift right (positive velocity)"){ + std::vector slice = {1.0, 2.0, 3.0, 4.0, 5.0}; + double velocity = 3.0; // dx = 1.0, dt = 1.0 -> disp = 3.0 -> shift by 3 + double dt = 1.0; + double ds = 1.0; + double bc_left = 0.0; + double bc_right = 0.0; + semiLagrangianAdvection1D(slice, velocity, dt, ds, bc_left, bc_right); + REQUIRE(slice[0] == Catch::Approx(0.0)); + REQUIRE(slice[1] == Catch::Approx(0.0)); + REQUIRE(slice[2] == Catch::Approx(0.0)); + REQUIRE(slice[3] == Catch::Approx(1.0)); + REQUIRE(slice[4] == Catch::Approx(2.0)); + } + + SECTION("Pure integer shift left (negative velocity)"){ + std::vector slice = {1.0, 2.0, 3.0, 4.0, 5.0}; + double velocity = -2.0; // dx = 1.0, dt = 1.0 -> disp = -2.0 -> shift left by 2 + double dt = 1.0; + double ds = 1.0; + double bc_left = 0.0; + double bc_right = 0.0; + semiLagrangianAdvection1D(slice, velocity, dt, ds, bc_left, bc_right); + REQUIRE(slice[0] == Catch::Approx(3.0)); + REQUIRE(slice[1] == Catch::Approx(4.0)); + REQUIRE(slice[2] == Catch::Approx(5.0)); + REQUIRE(slice[3] == Catch::Approx(0.0)); + REQUIRE(slice[4] == Catch::Approx(0.0)); + } + + SECTION("Inflow boundary condition padding"){ + std::vector slice = {1.0, 2.0, 3.0, 4.0, 5.0}; + double velocity = 2.0; + double dt = 1.0; + double ds = 1.0; + double bc_left = 9.9; + double bc_right = 0.0; + semiLagrangianAdvection1D(slice, velocity, dt, ds, bc_left, bc_right); + REQUIRE(slice[0] == Catch::Approx(9.9)); + REQUIRE(slice[1] == Catch::Approx(9.9)); + REQUIRE(slice[2] == Catch::Approx(1.0)); + REQUIRE(slice[3] == Catch::Approx(2.0)); + REQUIRE(slice[4] == Catch::Approx(3.0)); + } + + SECTION("Mass conservation on interior pulse"){ + int N = 50; + std::vector slice(N, 0.0); + slice[20] = 1.0; + slice[21] = 2.0; + slice[22] = 1.0; + double initial_mass = 4.0; + + double velocity = 15.5; // moves ~15.5 cells + double dt = 1.0; + double ds = 1.0; + semiLagrangianAdvection1D(slice, velocity, dt, ds, 0.0, 0.0); + + double final_mass = 0.0; + for (double val : slice) { + final_mass += val; + } + REQUIRE(final_mass == Catch::Approx(initial_mass).margin(1e-10)); + } + } } \ No newline at end of file From 68a143c28c18377fca01d7d92a5f90acc4cf79b4 Mon Sep 17 00:00:00 2001 From: "Sebastian D. Eastham" Date: Wed, 19 Aug 2026 17:53:33 -0400 Subject: [PATCH 02/14] Correction of error in S-L scheme and addition of docstring --- .../include/FVM_ANDS/AdvDiffSystem.hpp | 70 ++++++++++++++++++- Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp | 20 ++++-- 2 files changed, 83 insertions(+), 7 deletions(-) diff --git a/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp b/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp index 6c79dd3b9..772dbcc7b 100644 --- a/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp +++ b/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp @@ -11,7 +11,75 @@ namespace FVM_ANDS{ // Separate the SOR solver for testing without having to build an AdvDiffSystem object void sor_solve(const Eigen::SparseMatrix &A, const Eigen::VectorXd &rhs, Eigen::VectorXd &phi, double omega = 1.0, double threshold = 1e-3, int n_iters = 3); - // 1D Semi-Lagrangian advection helper + /** + * @brief 1D Flux-Form Semi-Lagrangian (FFSL) advection with Lax-Wendroff TVD subgrid reconstruction. + * + * @details + * Solves the 1D linear advection equation for a scalar field \f$\phi(s, t)\f$: + * \f[ + * \frac{\partial \phi}{\partial t} + v \frac{\partial \phi}{\partial s} = 0 + * \f] + * across a uniform grid with spacing \f$\Delta s\f$ and arbitrary timestep \f$\Delta t\f$ without CFL restrictions. + * + * ### Method Overview + * + * The method combines two key concepts from atmospheric transport and hyperbolic conservation laws: + * + * 1. **Integer-Shift Trajectory Decomposition (Ritchie, 1986)**: + * The advective displacement \f$\Delta s_{\text{total}} = v \Delta t\f$ is decomposed into an integer + * grid-cell translation \f$k\f$ and a subgrid residual displacement \f$\delta s\f$: + * \f[ + * k = \left\lfloor \frac{|v| \Delta t}{\Delta s} \right\rfloor \in \mathbb{Z}_{\ge 0}, \quad + * \delta s = |v| \Delta t - k \Delta s \in [0, \Delta s) + * \f] + * The fractional Courant number is \f$c_{\text{frac}} = \frac{\delta s}{\Delta s} \in [0, 1)\f$, + * and the residual timestep is \f$\Delta t_{\text{rem}} = \frac{\delta s}{|v|}\f$. + * - For \f$v > 0\f$: cells are shifted rightward by \f$k\f$ positions (\f$\phi_m \leftarrow \phi_{m-k}\f$), + * with inflow \f$m \in [0, k-1]\f$ padded by \f$\text{bc\_left}\f$. + * - For \f$v < 0\f$: cells are shifted leftward by \f$k\f$ positions (\f$\phi_m \leftarrow \phi_{m+k}\f$), + * with inflow \f$m \in [N-k, N-1]\f$ padded by \f$\text{bc\_right}\f$. + * + * 2. **Flux-Form Subgrid Advection with TVD Limiter (Lin & Rood, 1996; LeVeque, 2002)**: + * The remaining subgrid displacement is evolved via a single-step conservative finite-volume update: + * \f[ + * \phi_m^{n+1} = \phi_m^n - c_{\text{frac}} \left( F_{m+1/2} - F_{m-1/2} \right) + * \f] + * To maintain second-order accuracy in time and total-variation-diminishing (TVD) monotonicity, + * numerical interface fluxes \f$F_{m+1/2}\f$ represent the time-average of the characteristic departure + * interval \f$[x_{m+1/2} - v \Delta t_{\text{rem}}, \, x_{m+1/2}]\f$. As derived in LeVeque (2002, Ch. 6), + * evaluating the linear reconstruction at the centroid of this interval introduces the Lax-Wendroff factor + * \f$(1 - c_{\text{frac}})\f$: + * \f[ + * F_{m+1/2} = + * \begin{cases} + * \phi_m + \frac{1}{2} (1 - c_{\text{frac}}) \, \text{minmod}(\Delta \phi_{m-1/2}, \Delta \phi_{m+1/2}) & \text{if } v > 0 \\ + * \phi_{m+1} - \frac{1}{2} (1 - c_{\text{frac}}) \, \text{minmod}(\Delta \phi_{m+1/2}, \Delta \phi_{m+3/2}) & \text{if } v < 0 + * \end{cases} + * \f] + * where \f$\text{minmod}(a, b) = \text{sgn}(a) \max\left(0, \min(|a|, b \cdot \text{sgn}(a))\right)\f$. + * + * ### Key Properties + * - **Strict Mass Conservation**: Guaranteed by the conservative flux-differencing formulation (Lin & Rood, 1996). + * - **Monotonicity (TVD)**: The minmod limiter with the \f$(1 - c_{\text{frac}})\f$ Lax-Wendroff correction + * prevents spurious numerical oscillations for all Courant numbers. + * - **Unconditional Stability**: The integer translation ensures the residual Eulerian step always satisfies + * \f$c_{\text{frac}} < 1\f$. + * + * ### References + * - Ritchie, H. (1986). Eliminating the interpolation associated with the semi-Lagrangian scheme. + * *Monthly Weather Review*, 114(1), 135–146. + * - Lin, S.-J., & Rood, R. B. (1996). Multidimensional flux-form semi-Lagrangian transport schemes. + * *Monthly Weather Review*, 124(9), 2046–2070. + * - LeVeque, R. J. (2002). *Finite Volume Methods for Hyperbolic Problems*. Cambridge University Press, + * Chapters 6 (TVD Limiters) & 9 (Variable-Coefficient and Large Time-Step Methods). + * + * @param[in,out] slice 1D vector of cell-centered scalar values across the slice (modified in place). + * @param[in] velocity Advecting velocity (\f$v\f$) along the coordinate direction [m/s]. + * @param[in] dt Advective timestep (\f$\Delta t\f$) [s]. + * @param[in] ds Grid cell spacing (\f$\Delta s\f$) [m]. + * @param[in] bc_left Dirichlet boundary value at the left (inflow/outflow) face. + * @param[in] bc_right Dirichlet boundary value at the right (inflow/outflow) face. + */ void semiLagrangianAdvection1D(std::vector& slice, double velocity, double dt, double ds, double bc_left, double bc_right); struct AdvDiffParams { diff --git a/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp b/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp index 7e5cae952..abce979c5 100644 --- a/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp +++ b/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp @@ -541,6 +541,10 @@ namespace FVM_ANDS{ applyBoundaryCondition(); //need this to calculate minmod function at some timestep. } + /** + * @brief 1D Flux-Form Semi-Lagrangian (FFSL) advection with Lax-Wendroff TVD subgrid reconstruction. + * @see AdvDiffSystem.hpp for detailed algorithmic documentation and literature references. + */ void semiLagrangianAdvection1D( std::vector& slice, double velocity, @@ -576,6 +580,9 @@ namespace FVM_ANDS{ // Step 2: Fractional Forward Euler step if (rem_disp > 1.0e-12 && rem_dt > 0.0) { + const double cfl_frac = velocity * rem_dt / ds; + const double slope_weight = 0.5 * (1.0 - cfl_frac); + std::vector face_flux(N + 1, 0.0); face_flux[0] = bc_left; @@ -588,16 +595,15 @@ namespace FVM_ANDS{ double diff_up = (m == 0) ? (2.0 * (slice[0] - bc_left)) : (slice[m] - slice[m - 1]); double diff_down = slice[m + 1] - slice[m]; double slope = minmod(diff_up, diff_down); - face_flux[m + 1] = slice[m] + 0.5 * slope; + face_flux[m + 1] = slice[m] + slope_weight * slope; } // Outflow face at m = N double diff_up_last = (N >= 2) ? (slice[N - 1] - slice[N - 2]) : (2.0 * (slice[0] - bc_left)); double diff_down_last = 2.0 * (bc_right - slice[N - 1]); double slope_last = minmod(diff_up_last, diff_down_last); - face_flux[N] = slice[N - 1] + 0.5 * slope_last; + face_flux[N] = slice[N - 1] + slope_weight * slope_last; - double cfl_frac = velocity * rem_dt / ds; for (int m = 0; m < N; ++m) { slice[m] -= cfl_frac * (face_flux[m + 1] - face_flux[m]); } @@ -625,6 +631,9 @@ namespace FVM_ANDS{ // Step 2: Fractional Forward Euler step if (rem_disp > 1.0e-12 && rem_dt > 0.0) { + const double cfl_frac = abs_vel * rem_dt / ds; + const double slope_weight = 0.5 * (1.0 - cfl_frac); + std::vector face_flux(N + 1, 0.0); face_flux[N] = bc_right; @@ -637,16 +646,15 @@ namespace FVM_ANDS{ double diff_up = (m + 1 == N - 1) ? (2.0 * (bc_right - slice[N - 1])) : (slice[m + 2] - slice[m + 1]); double diff_down = slice[m + 1] - slice[m]; double slope = minmod(diff_down, diff_up); - face_flux[m + 1] = slice[m + 1] - 0.5 * slope; + face_flux[m + 1] = slice[m + 1] - slope_weight * slope; } // Outflow face at m = 0 double diff_up_0 = 2.0 * (slice[0] - bc_left); double diff_down_0 = (N >= 2) ? (slice[1] - slice[0]) : diff_up_0; double slope_0 = minmod(diff_down_0, diff_up_0); - face_flux[0] = slice[0] - 0.5 * slope_0; + face_flux[0] = slice[0] - slope_weight * slope_0; - double cfl_frac = abs_vel * rem_dt / ds; for (int m = 0; m < N; ++m) { slice[m] -= cfl_frac * (face_flux[m] - face_flux[m + 1]); } From ed9fc4e6fc47287351812ce574178281a21260a3 Mon Sep 17 00:00:00 2001 From: "Sebastian D. Eastham" Date: Wed, 19 Aug 2026 19:10:05 -0400 Subject: [PATCH 03/14] Fixed OpenMP handling in advection --- Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp | 2 +- Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp | 6 +++--- Code.v05-00/src/FVM_ANDS/FVM_Solver.cpp | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp b/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp index 772dbcc7b..fd67fe20f 100644 --- a/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp +++ b/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp @@ -106,7 +106,7 @@ namespace FVM_ANDS{ const Eigen::VectorXd& calcRHS(); void applyBoundaryCondition(); void updateBoundaryCondition(const BoundaryConditions& bc); - void semiLagrangianAdvection(double dt); + void semiLagrangianAdvection(double dt, bool parallelAdvection = false); Eigen::VectorXd forwardEulerAdvection(bool operatorSplit = false, bool parallelAdvection = false) const noexcept; // Breakup the implementation of sor_solve to allow for easy testing by inputing an arbitrary linear system to solve: // Implementation is moved outside of the class, and make class method to be used in code diff --git a/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp b/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp index abce979c5..678a5c8d3 100644 --- a/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp +++ b/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp @@ -662,9 +662,9 @@ namespace FVM_ANDS{ } } - void AdvDiffSystem::semiLagrangianAdvection(double dt) { + void AdvDiffSystem::semiLagrangianAdvection(double dt, bool parallelAdvection) { // 1. Horizontal Advection along X (row by row) - #pragma omp parallel for default(shared) schedule(static) + #pragma omp parallel for if (parallelAdvection) default(shared) schedule(static) for (int j = 0; j < ny_; ++j) { double u_j = u_double_ - yCoord_[j] * shear_; if (std::abs(u_j) > 1.0e-14) { @@ -685,7 +685,7 @@ namespace FVM_ANDS{ // 2. Vertical Advection along Y (column by column) if (std::abs(v_double_) > 1.0e-14) { - #pragma omp parallel for default(shared) schedule(static) + #pragma omp parallel for if (parallelAdvection) default(shared) schedule(static) for (int i = 0; i < nx_; ++i) { std::vector col(ny_); for (int j = 0; j < ny_; ++j) { diff --git a/Code.v05-00/src/FVM_ANDS/FVM_Solver.cpp b/Code.v05-00/src/FVM_ANDS/FVM_Solver.cpp index dd906b060..66def456b 100644 --- a/Code.v05-00/src/FVM_ANDS/FVM_Solver.cpp +++ b/Code.v05-00/src/FVM_ANDS/FVM_Solver.cpp @@ -74,7 +74,7 @@ namespace FVM_ANDS{ #endif //Step 1: Solve Advection for half timestep via Semi-Lagrangian advection - advDiffSys_.semiLagrangianAdvection(0.5 * dt_max); + advDiffSys_.semiLagrangianAdvection(0.5 * dt_max, parallelAdvection); #ifdef ENABLE_TIMING auto stop = std::chrono::high_resolution_clock::now(); @@ -127,7 +127,7 @@ namespace FVM_ANDS{ #endif //Step 3: Solve advection for second half timestep via Semi-Lagrangian advection - advDiffSys_.semiLagrangianAdvection(0.5 * dt_max); + advDiffSys_.semiLagrangianAdvection(0.5 * dt_max, parallelAdvection); #ifdef ENABLE_TIMING stop = std::chrono::high_resolution_clock::now(); From 45bf28b2843dbbcf824fb85b97c98684fe60a1c3 Mon Sep 17 00:00:00 2001 From: "Sebastian D. Eastham" Date: Wed, 19 Aug 2026 19:51:06 -0400 Subject: [PATCH 04/14] Fix variable name error in AdvDiffSystem - zero diff --- Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp b/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp index 678a5c8d3..cabe76072 100644 --- a/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp +++ b/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp @@ -649,10 +649,10 @@ namespace FVM_ANDS{ face_flux[m + 1] = slice[m + 1] - slope_weight * slope; } - // Outflow face at m = 0 - double diff_up_0 = 2.0 * (slice[0] - bc_left); - double diff_down_0 = (N >= 2) ? (slice[1] - slice[0]) : diff_up_0; - double slope_0 = minmod(diff_down_0, diff_up_0); + // Outflow face at m = 0 (leftward flow towards bc_left) + double diff_down_0 = 2.0 * (slice[0] - bc_left); + double diff_up_0 = (N >= 2) ? (slice[1] - slice[0]) : diff_down_0; + double slope_0 = minmod(diff_up_0, diff_down_0); face_flux[0] = slice[0] - slope_weight * slope_0; for (int m = 0; m < N; ++m) { From 07f5765b46dafd7b020b3862968999fbe1a7306a Mon Sep 17 00:00:00 2001 From: "Sebastian D. Eastham" Date: Wed, 19 Aug 2026 20:40:56 -0400 Subject: [PATCH 05/14] Improved parallelisation in AdvDiff - zero diff --- Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp | 56 +++++++++++----------- 1 file changed, 29 insertions(+), 27 deletions(-) diff --git a/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp b/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp index cabe76072..70a3d4451 100644 --- a/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp +++ b/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp @@ -664,40 +664,42 @@ namespace FVM_ANDS{ void AdvDiffSystem::semiLagrangianAdvection(double dt, bool parallelAdvection) { // 1. Horizontal Advection along X (row by row) - #pragma omp parallel for if (parallelAdvection) default(shared) schedule(static) - for (int j = 0; j < ny_; ++j) { - double u_j = u_double_ - yCoord_[j] * shear_; - if (std::abs(u_j) > 1.0e-14) { - std::vector row(nx_); - for (int i = 0; i < nx_; ++i) { - int idx = twoDIdx_to_vecIdx(i, j, nx_, ny_, format_); - row[i] = phi_[idx]; - } - double bc_left = bcVals_left_.empty() ? 0.0 : bcVals_left_[j]; - double bc_right = bcVals_right_.empty() ? 0.0 : bcVals_right_[j]; - semiLagrangianAdvection1D(row, u_j, dt, dx_, bc_left, bc_right); - for (int i = 0; i < nx_; ++i) { - int idx = twoDIdx_to_vecIdx(i, j, nx_, ny_, format_); - phi_[idx] = row[i]; + #pragma omp parallel if (parallelAdvection) default(shared) + { + std::vector row(nx_); + #pragma omp for schedule(static) + for (int j = 0; j < ny_; ++j) { + double u_j = u_double_ - yCoord_[j] * shear_; + if (std::abs(u_j) > 1.0e-14) { + for (int i = 0; i < nx_; ++i) { + int idx = twoDIdx_to_vecIdx(i, j, nx_, ny_, format_); + row[i] = phi_[idx]; + } + semiLagrangianAdvection1D(row, u_j, dt, dx_, bcVals_left_[j], bcVals_right_[j]); + for (int i = 0; i < nx_; ++i) { + int idx = twoDIdx_to_vecIdx(i, j, nx_, ny_, format_); + phi_[idx] = row[i]; + } } } } // 2. Vertical Advection along Y (column by column) if (std::abs(v_double_) > 1.0e-14) { - #pragma omp parallel for if (parallelAdvection) default(shared) schedule(static) - for (int i = 0; i < nx_; ++i) { + #pragma omp parallel if (parallelAdvection) default(shared) + { std::vector col(ny_); - for (int j = 0; j < ny_; ++j) { - int idx = twoDIdx_to_vecIdx(i, j, nx_, ny_, format_); - col[j] = phi_[idx]; - } - double bc_bot = bcVals_bot_.empty() ? 0.0 : bcVals_bot_[i]; - double bc_top = bcVals_top_.empty() ? 0.0 : bcVals_top_[i]; - semiLagrangianAdvection1D(col, v_double_, dt, dy_, bc_bot, bc_top); - for (int j = 0; j < ny_; ++j) { - int idx = twoDIdx_to_vecIdx(i, j, nx_, ny_, format_); - phi_[idx] = col[j]; + #pragma omp for schedule(static) + for (int i = 0; i < nx_; ++i) { + for (int j = 0; j < ny_; ++j) { + int idx = twoDIdx_to_vecIdx(i, j, nx_, ny_, format_); + col[j] = phi_[idx]; + } + semiLagrangianAdvection1D(col, v_double_, dt, dy_, bcVals_bot_[i], bcVals_top_[i]); + for (int j = 0; j < ny_; ++j) { + int idx = twoDIdx_to_vecIdx(i, j, nx_, ny_, format_); + phi_[idx] = col[j]; + } } } } From 5427df37459875b4366f79a016f9fb22fc8c4945 Mon Sep 17 00:00:00 2001 From: "Sebastian D. Eastham" Date: Thu, 20 Aug 2026 08:44:07 -0400 Subject: [PATCH 06/14] Fixed bug in definition of diff_up_0 @lrobion identified an incorrect definition for the upwind calculation at the boundary. This has no effect outside of exceptional cases but was misleading and made the code hard to understand. Now corrected - zero-diff for almost all cases so should have no effect on the user. --- Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp b/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp index 70a3d4451..b6457201c 100644 --- a/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp +++ b/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp @@ -651,7 +651,7 @@ namespace FVM_ANDS{ // Outflow face at m = 0 (leftward flow towards bc_left) double diff_down_0 = 2.0 * (slice[0] - bc_left); - double diff_up_0 = (N >= 2) ? (slice[1] - slice[0]) : diff_down_0; + double diff_up_0 = (N >= 2) ? (slice[1] - slice[0]) : (2.0 * (bc_right - slice[0])); double slope_0 = minmod(diff_up_0, diff_down_0); face_flux[0] = slice[0] - slope_weight * slope_0; From af408bbeaef5ada808d575195dd8bc417cd17097 Mon Sep 17 00:00:00 2001 From: "Sebastian D. Eastham" Date: Thu, 20 Aug 2026 13:16:47 -0400 Subject: [PATCH 07/14] Monotonicity and TVD tests for S-L advection --- Code.v05-00/tests/test_adv_diff_solver.cpp | 49 ++++++++++++++++++++++ 1 file changed, 49 insertions(+) diff --git a/Code.v05-00/tests/test_adv_diff_solver.cpp b/Code.v05-00/tests/test_adv_diff_solver.cpp index 313e82ecf..df9f5e26e 100644 --- a/Code.v05-00/tests/test_adv_diff_solver.cpp +++ b/Code.v05-00/tests/test_adv_diff_solver.cpp @@ -572,4 +572,53 @@ namespace FVM_ANDS{ REQUIRE(final_mass == Catch::Approx(initial_mass).margin(1e-10)); } } + + TEST_CASE("Semi-Lagrangian 1D Advection preserves monotonicity", "[advection]"){ + + // Monotone non-decreasing profile + const std::vector initial = {0.0, 1.0, 3.0, 4.0, 5.0, 5.0}; + const double dt = 1.0; + const double ds = 1.0; + const double bc_left = 0.0; + const double bc_right = 5.0; + + auto checkMonotonicity = [&](const std::vector& slice){ + // 1. Values should stay in the range spanned by the initial data and the BCs. + for (std::size_t m = 0; m < slice.size(); m++) { + INFO("cell " << m << " = " << slice[m]); + REQUIRE(slice[m] >= 0.0 - 1e-12); + REQUIRE(slice[m] <= 5.0 + 1e-12); + } + + // 2. A monotone profile must stay monotone + for (std::size_t m = 1; m < slice.size(); m++) { + INFO("cells " << m - 1 << ", " << m << " = " << slice[m-1] << ", " << slice[m]); + REQUIRE(slice[m] >= slice[m-1] - 1e-12); + } + + // 3. Total variation must not increase + double tv_initial = 0.0; + double tv_final = 0.0; + for (std::size_t m = 1; m < slice.size(); m++) { + tv_initial += std::abs(initial[m] - initial[m-1]); + tv_final += std::abs(slice[m] - slice[m-1]); + } + INFO("TV before = " << tv_initial << ", TV after = " << tv_final); + REQUIRE(tv_final <= tv_initial + 1e-12); + }; + + SECTION("Fractional CFL 0.80, above the 2/3 TVD limit"){ + std::vector slice = initial; + double velocity = 0.8; // dt = ds = 1 -> no integer shift, fractional CFL = 0.8 + semiLagrangianAdvection1D(slice, velocity, dt, ds, bc_left, bc_right); + checkMonotonicity(slice); + } + + SECTION("Fractional CFL 0.60, below the 2/3 TVD limit"){ + std::vector slice = initial; + double velocity = 0.6; // dt = ds = 1 -> no integer shift, fractional CFL = 0.6 + semiLagrangianAdvection1D(slice, velocity, dt, ds, bc_left, bc_right); + checkMonotonicity(slice); + } + } } \ No newline at end of file From 72d425ce488b29f40bb5a6caa8977f24ee4f2cd6 Mon Sep 17 00:00:00 2001 From: "Sebastian D. Eastham" Date: Wed, 19 Aug 2026 07:07:57 -0400 Subject: [PATCH 08/14] Implemented time-averaged diffusion and ice growth substepping --- Code.v05-00/defaults/input.yaml | 1 + Code.v05-00/include/Core/Input_Mod.hpp | 1 + Code.v05-00/include/Defaults/Input.hpp | 1 + Code.v05-00/include/Util/PlumeModelUtils.hpp | 5 +- Code.v05-00/src/Core/LAGRIDPlumeModel.cpp | 17 +++- Code.v05-00/src/Util/PlumeModelUtils.cpp | 83 +++++++++++++++++-- .../src/YamlInputReader/YamlInputReader.cpp | 5 ++ Code.v05-00/tests/test_yamlreader.cpp | 1 + 8 files changed, 105 insertions(+), 9 deletions(-) diff --git a/Code.v05-00/defaults/input.yaml b/Code.v05-00/defaults/input.yaml index f6e60f499..c4864a745 100644 --- a/Code.v05-00/defaults/input.yaml +++ b/Code.v05-00/defaults/input.yaml @@ -107,6 +107,7 @@ AEROSOL MENU: # Keep on Turn on ice growth (T/F): T Ice growth timestep [min] (double): 1 + Ice growth substep [s] (double): 60.0 # At least one of "Use met. input", "Impose moist layer depth", or "Impose lapse rate" must be true # Imposing moist layer depth will automatically calculate the lapse rate and override the imposed lapse rate diff --git a/Code.v05-00/include/Core/Input_Mod.hpp b/Code.v05-00/include/Core/Input_Mod.hpp index 470ab310e..cebee1aed 100644 --- a/Code.v05-00/include/Core/Input_Mod.hpp +++ b/Code.v05-00/include/Core/Input_Mod.hpp @@ -79,6 +79,7 @@ struct OptInput bool AEROSOL_COAGULATION_LIQUID; bool AEROSOL_ICE_GROWTH; double AEROSOL_ICE_GROWTH_TIMESTEP; + double AEROSOL_ICE_GROWTH_SUBSTEP; /* ========================================== */ /* ---- METEOROLOGY MENU -------------------- */ diff --git a/Code.v05-00/include/Defaults/Input.hpp b/Code.v05-00/include/Defaults/Input.hpp index 583eea0e4..99e4f8da5 100644 --- a/Code.v05-00/include/Defaults/Input.hpp +++ b/Code.v05-00/include/Defaults/Input.hpp @@ -108,6 +108,7 @@ AEROSOL MENU: # Keep on Turn on ice growth (T/F): T Ice growth timestep [min] (double): 1 + Ice growth substep [s] (double): 60.0 # At least one of "Use met. input", "Impose moist layer depth", or "Impose lapse rate" must be true # Imposing moist layer depth will automatically calculate the lapse rate and override the imposed lapse rate diff --git a/Code.v05-00/include/Util/PlumeModelUtils.hpp b/Code.v05-00/include/Util/PlumeModelUtils.hpp index f89314b4f..3632d9024 100644 --- a/Code.v05-00/include/Util/PlumeModelUtils.hpp +++ b/Code.v05-00/include/Util/PlumeModelUtils.hpp @@ -18,8 +18,11 @@ namespace PlumeModelUtils { double &v_x, double &v_y, \ double &dTrav_x, double &dTrav_y ); + void DiffParam( const double time, const double dt, double &d_x, double &d_y, \ + const double D_X, const double D_Y ); + void DiffParam( const double time, double &d_x, double &d_y, \ - const double D_X, const double D_Y ); + const double D_X, const double D_Y ); } diff --git a/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp b/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp index 7928aaf17..5d9d667ba 100644 --- a/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp +++ b/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp @@ -143,7 +143,17 @@ SimStatus LAGRIDPlumeModel::runFullModel() { #endif timestepVars_.lastTimeIceGrowth = timestepVars_.curr_Time_s + timestepVars_.dt; - iceAerosol_.Grow( timestepVars_.ICE_GROWTH_DT, H2O_, met_.Temp(), met_.Press()); + + // Substep ice growth in increments of <= dt_micro_target seconds (default 60s) + const double dt_micro_target = (optInput_.AEROSOL_ICE_GROWTH_SUBSTEP > 0.0) + ? optInput_.AEROSOL_ICE_GROWTH_SUBSTEP : 60.0; + const double dt_growth_total = timestepVars_.ICE_GROWTH_DT; + const int n_growth_substeps = std::max(1, static_cast(std::ceil(dt_growth_total / dt_micro_target))); + const double dt_growth_sub = dt_growth_total / n_growth_substeps; + + for (int sub = 0; sub < n_growth_substeps; ++sub) { + iceAerosol_.Grow(dt_growth_sub, H2O_, met_.Temp(), met_.Press()); + } #ifdef ENABLE_TIMING auto icegrowth_end = std::chrono::high_resolution_clock::now(); @@ -422,8 +432,9 @@ void LAGRIDPlumeModel::initH2O() { void LAGRIDPlumeModel::updateDiffVecs() { double dh_enhanced, dv_enhanced; - // Update Diffusion - PlumeModelUtils::DiffParam( timestepVars_.curr_Time_s - timestepVars_.tInitial_s + timestepVars_.TRANSPORT_DT / 2.0, + // Update Diffusion with exact time-averaged diffusivity over the transport step + const double time_start = timestepVars_.curr_Time_s - timestepVars_.tInitial_s; + PlumeModelUtils::DiffParam( time_start, timestepVars_.TRANSPORT_DT, dh_enhanced, dv_enhanced, input_.horizDiff(), input_.vertiDiff() ); auto number = iceAerosol_.TotalNumber(); auto num_max = VectorUtils::max(number); diff --git a/Code.v05-00/src/Util/PlumeModelUtils.cpp b/Code.v05-00/src/Util/PlumeModelUtils.cpp index 55e009deb..2ebfa9506 100644 --- a/Code.v05-00/src/Util/PlumeModelUtils.cpp +++ b/Code.v05-00/src/Util/PlumeModelUtils.cpp @@ -118,12 +118,87 @@ namespace PlumeModelUtils { } /* End of AdvGlobal */ + void DiffParam( const double time, const double dt, double &d_x, double &d_y, \ + const double D_X, const double D_Y ) + { + + /* DiffParam: + * Computes time-averaged diffusion parameters over [time, time + dt] + * + * INPUTS: + * (double) time: current time since simulation started in [s] + * (double) dt: transport timestep duration in [s] + * + * OUTPUTS: + * (double) d_x: average horizontal diffusion coefficient over [time, time + dt] in [m^2/s] + * (double) d_y: average vertical diffusion coefficient over [time, time + dt] in [m^2/s] + */ + + if ( dt <= 0.0 ) { + DiffParam( time, d_x, d_y, D_X, D_Y ); + return; + } + + const double t1 = std::max(0.0, time); + const double t2 = t1 + dt; + + if ( DPROF == 0 ) { + // Piecewise constant enhancement + const double dt_enh_x = std::max(0.0, std::min(t2, tH0) - std::min(t1, tH0)); + const double dt_amb_x = dt - dt_enh_x; + d_x = ( 1.13 * D_X * dt_enh_x + D_X * dt_amb_x ) / dt; + + const double dt_enh_y = std::max(0.0, std::min(t2, tV0) - std::min(t1, tV0)); + const double dt_amb_y = dt - dt_enh_y; + d_y = ( 7.0 * D_Y * dt_enh_y + D_Y * dt_amb_y ) / dt; + } + else if ( DPROF == 1 ) { + // Linear decay + auto int_linear = [](double t, double D_0, double D_amb, double t_0) { + if ( t <= 0.0 ) return 0.0; + if ( t <= t_0 ) { + return D_0 * t + 0.5 * (D_amb - D_0) * (t * t) / t_0; + } else { + double int_to_t0 = 0.5 * (D_0 + D_amb) * t_0; + return int_to_t0 + D_amb * (t - t_0); + } + }; + d_x = ( int_linear(t2, 1.13 * D_X, D_X, tH0) - int_linear(t1, 1.13 * D_X, D_X, tH0) ) / dt; + d_y = ( int_linear(t2, 7.00 * D_Y, D_Y, tV0) - int_linear(t1, 7.00 * D_Y, D_Y, tV0) ) / dt; + } + else if ( DPROF == 2 ) { + // Exponential decay + auto avg_exp = [](double t_start, double t_end, double delta_t, double D_0, double D_amb, double t_0) { + return D_amb + (D_0 - D_amb) * (t_0 / delta_t) * (std::exp(-t_start / t_0) - std::exp(-t_end / t_0)); + }; + d_x = avg_exp(t1, t2, dt, 1.13 * D_X, D_X, tH0); + d_y = avg_exp(t1, t2, dt, 7.00 * D_Y, D_Y, tV0); + } + else { + std::string const currFile("DiffParam.cpp"); + std::cout << "ERROR: In " << currFile << ": DPROF set to " << DPROF << "\n"; + } + + if ( d_x < 0.0 ) { + std::cout << "d_x is negative: d_x = " << d_x << " [m^2/s]" << "\n"; + std::cout << "Setting d_x to 0.0" << "\n"; + d_x = 0.0; + } + + if ( d_y < 0.0 ) { + std::cout << "d_y is negative: d_y = " << d_y << " [m^2/s]" << "\n"; + std::cout << "Setting d_y to 0.0" << "\n"; + d_y = 0.0; + } + + } /* End of DiffParam (time, dt) */ + void DiffParam( const double time, double &d_x, double &d_y, \ const double D_X, const double D_Y ) { /* DiffParam: - * Computes diffusion parameters + * Computes instantaneous diffusion parameters at time * * INPUTS: * (double) time: current time since simulation started in [s] @@ -132,7 +207,7 @@ namespace PlumeModelUtils { * (double) d_x: current horizontal diffusion coefficient at time in [m^2/s] * (double) d_y: current vertical diffusion coefficient at time in [m^2/s] */ - + if ( DPROF == 0 ) { if ( time <= tH0 ) d_x = 1.13 * D_X; @@ -168,8 +243,6 @@ namespace PlumeModelUtils { d_y = 0.0; } - - } /* End of DiffParam */ - + } /* End of DiffParam (time) */ } \ No newline at end of file diff --git a/Code.v05-00/src/YamlInputReader/YamlInputReader.cpp b/Code.v05-00/src/YamlInputReader/YamlInputReader.cpp index e4036f3ed..f68a956ea 100644 --- a/Code.v05-00/src/YamlInputReader/YamlInputReader.cpp +++ b/Code.v05-00/src/YamlInputReader/YamlInputReader.cpp @@ -394,6 +394,11 @@ namespace YamlInputReader{ input.AEROSOL_COAGULATION_LIQUID = parseBoolString(aeroNode["Turn on liquid coagulation (T/F)"].as(), "Turn on liquid coagulation (T/F)"); input.AEROSOL_ICE_GROWTH = parseBoolString(aeroNode["Turn on ice growth (T/F)"].as(), "Turn on ice growth (T/F)"); input.AEROSOL_ICE_GROWTH_TIMESTEP = parseDoubleString(aeroNode["Ice growth timestep [min] (double)"].as(), "Ice growth timestep [min] (double)"); + if (aeroNode["Ice growth substep [s] (double)"]) { + input.AEROSOL_ICE_GROWTH_SUBSTEP = parseDoubleString(aeroNode["Ice growth substep [s] (double)"].as(), "Ice growth substep [s] (double)"); + } else { + input.AEROSOL_ICE_GROWTH_SUBSTEP = 60.0; + } } void readMetMenu(OptInput& input, const YAML::Node& metNode){ YAML::Node metInputSubmenu = metNode["METEOROLOGICAL INPUT SUBMENU"]; diff --git a/Code.v05-00/tests/test_yamlreader.cpp b/Code.v05-00/tests/test_yamlreader.cpp index b0d7a10f7..ba54c4890 100644 --- a/Code.v05-00/tests/test_yamlreader.cpp +++ b/Code.v05-00/tests/test_yamlreader.cpp @@ -269,6 +269,7 @@ TEST_CASE("Read Yaml File"){ REQUIRE(input.AEROSOL_COAGULATION_LIQUID == true); REQUIRE(input.AEROSOL_ICE_GROWTH == true); REQUIRE(input.AEROSOL_ICE_GROWTH_TIMESTEP == 10); + REQUIRE(input.AEROSOL_ICE_GROWTH_SUBSTEP == 60.0); } SECTION("Read Met Menu"){ OptInput input; From 02fea53d711fba593912e2f55549db6dc01ebaf2 Mon Sep 17 00:00:00 2001 From: "Sebastian D. Eastham" Date: Wed, 19 Aug 2026 07:21:44 -0400 Subject: [PATCH 09/14] Implemented consistent substepping for transport and ice growth --- Code.v05-00/defaults/input.yaml | 2 +- Code.v05-00/include/Core/Input_Mod.hpp | 2 +- Code.v05-00/include/Core/LAGRIDPlumeModel.hpp | 4 +- Code.v05-00/include/Defaults/Input.hpp | 2 +- Code.v05-00/src/Core/LAGRIDPlumeModel.cpp | 80 +++++++++---------- .../src/YamlInputReader/YamlInputReader.cpp | 10 +-- Code.v05-00/tests/test_yamlreader.cpp | 2 +- 7 files changed, 47 insertions(+), 55 deletions(-) diff --git a/Code.v05-00/defaults/input.yaml b/Code.v05-00/defaults/input.yaml index c4864a745..3d25587d9 100644 --- a/Code.v05-00/defaults/input.yaml +++ b/Code.v05-00/defaults/input.yaml @@ -85,6 +85,7 @@ TRANSPORT MENU: # Outdated, not used (was used by spectral solver) Fill Negative Values (T/F): T Transport Timestep [min] (double): 1 + Transport and ice growth substep [s] (double): 60.0 # Keep off: not sure of the effect yet + met updraft is included (if met file input) PLUME UPDRAFT SUBMENU: Turn on plume updraft (T/F): F @@ -107,7 +108,6 @@ AEROSOL MENU: # Keep on Turn on ice growth (T/F): T Ice growth timestep [min] (double): 1 - Ice growth substep [s] (double): 60.0 # At least one of "Use met. input", "Impose moist layer depth", or "Impose lapse rate" must be true # Imposing moist layer depth will automatically calculate the lapse rate and override the imposed lapse rate diff --git a/Code.v05-00/include/Core/Input_Mod.hpp b/Code.v05-00/include/Core/Input_Mod.hpp index cebee1aed..921507360 100644 --- a/Code.v05-00/include/Core/Input_Mod.hpp +++ b/Code.v05-00/include/Core/Input_Mod.hpp @@ -58,6 +58,7 @@ struct OptInput bool TRANSPORT_TRANSPORT; bool TRANSPORT_FILL; double TRANSPORT_TIMESTEP; + double TRANSPORT_ICE_GROWTH_SUBSTEP; bool TRANSPORT_UPDRAFT; double TRANSPORT_UPDRAFT_TIMESCALE; double TRANSPORT_UPDRAFT_VELOCITY; @@ -79,7 +80,6 @@ struct OptInput bool AEROSOL_COAGULATION_LIQUID; bool AEROSOL_ICE_GROWTH; double AEROSOL_ICE_GROWTH_TIMESTEP; - double AEROSOL_ICE_GROWTH_SUBSTEP; /* ========================================== */ /* ---- METEOROLOGY MENU -------------------- */ diff --git a/Code.v05-00/include/Core/LAGRIDPlumeModel.hpp b/Code.v05-00/include/Core/LAGRIDPlumeModel.hpp index 260350a0e..41ed589f6 100644 --- a/Code.v05-00/include/Core/LAGRIDPlumeModel.hpp +++ b/Code.v05-00/include/Core/LAGRIDPlumeModel.hpp @@ -97,8 +97,8 @@ class LAGRIDPlumeModel { void initializeGrid(const EPM::Output &epmOut); void saveTSAerosol(); void initH2O(); - void updateDiffVecs(); - void runTransport(double timestep); + void updateDiffVecs(double time_start, double dt); + void runTransport(double timestep, double time_start); void remapAllVars(double remapTimestep, const std::vector>& mask, const VectorUtils::MaskInfo& maskInfo); std::pair remapVariable(const VectorUtils::MaskInfo& maskInfo, const BufferInfo& buffers, const Vector_2D& phi, const std::vector>& mask); double totalAirMass(); diff --git a/Code.v05-00/include/Defaults/Input.hpp b/Code.v05-00/include/Defaults/Input.hpp index 99e4f8da5..ce9c151e7 100644 --- a/Code.v05-00/include/Defaults/Input.hpp +++ b/Code.v05-00/include/Defaults/Input.hpp @@ -86,6 +86,7 @@ TRANSPORT MENU: # Outdated, not used (was used by spectral solver) Fill Negative Values (T/F): T Transport Timestep [min] (double): 1 + Transport and ice growth substep [s] (double): 60.0 # Keep off: not sure of the effect yet + met updraft is included (if met file input) PLUME UPDRAFT SUBMENU: Turn on plume updraft (T/F): F @@ -108,7 +109,6 @@ AEROSOL MENU: # Keep on Turn on ice growth (T/F): T Ice growth timestep [min] (double): 1 - Ice growth substep [s] (double): 60.0 # At least one of "Use met. input", "Impose moist layer depth", or "Impose lapse rate" must be true # Imposing moist layer depth will automatically calculate the lapse rate and override the imposed lapse rate diff --git a/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp b/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp index 5d9d667ba..c4726f3d4 100644 --- a/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp +++ b/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp @@ -104,22 +104,43 @@ SimStatus LAGRIDPlumeModel::runFullModel() { std::cout << "\n - Time step: " << timestepVars_.nTime + 1 << " out of " << timestepVars_.timeArray.size(); std::cout << "\n -> Solar time: " << std::fmod( timestepVars_.curr_Time_s/3600.0, 24.0 ) << " [hr]" << std::endl; - // Run Transport - std::cout << "Running Transport..." << std::endl; + // Interleaved Transport and Ice Growth Subcycling bool timeForTransport = (simVars_.TRANSPORT && (timestepVars_.nTime == 0 || timestepVars_.checkTimeForTransport())); - if (timeForTransport) { - timestepVars_.lastTimeTransport = timestepVars_.curr_Time_s + timestepVars_.dt; + bool timeForIceGrowth = (simVars_.ICE_GROWTH && timestepVars_.checkTimeForIceGrowth()); + + if (timeForTransport || timeForIceGrowth) { + double dt_total = timeForTransport ? timestepVars_.TRANSPORT_DT : timestepVars_.ICE_GROWTH_DT; + const double dt_sub_target = (optInput_.TRANSPORT_ICE_GROWTH_SUBSTEP > 0.0) + ? optInput_.TRANSPORT_ICE_GROWTH_SUBSTEP : 60.0; + const int n_subcycles = std::max(1, static_cast(std::ceil(dt_total / dt_sub_target))); + const double dt_sub = dt_total / n_subcycles; + + if (timeForTransport) { + std::cout << "Running Transport and ice growth subcycling (" << n_subcycles << " x " << dt_sub << " s)..." << std::endl; + timestepVars_.lastTimeTransport = timestepVars_.curr_Time_s + timestepVars_.dt; + } + if (timeForIceGrowth) { + timestepVars_.lastTimeIceGrowth = timestepVars_.curr_Time_s + timestepVars_.dt; + } #ifdef ENABLE_TIMING - auto transport_start = std::chrono::high_resolution_clock::now(); + auto subcycling_start = std::chrono::high_resolution_clock::now(); #endif - runTransport(timestepVars_.TRANSPORT_DT); + for (int sub = 0; sub < n_subcycles; ++sub) { + const double t_sub_start = (timestepVars_.curr_Time_s - timestepVars_.tInitial_s) + sub * dt_sub; + if (timeForTransport) { + runTransport(dt_sub, t_sub_start); + } + if (timeForIceGrowth) { + iceAerosol_.Grow(dt_sub, H2O_, met_.Temp(), met_.Press()); + } + } #ifdef ENABLE_TIMING - auto transport_end = std::chrono::high_resolution_clock::now(); - auto transport_duration = std::chrono::duration_cast(transport_end - transport_start); - std::cout << " Ran transport in " << transport_duration.count() << " ms" << std::endl; + auto subcycling_end = std::chrono::high_resolution_clock::now(); + auto subcycling_duration = std::chrono::duration_cast(subcycling_end - subcycling_start); + std::cout << " Ran transport & ice growth subcycling in " << subcycling_duration.count() << " ms" << std::endl; #endif } @@ -134,34 +155,6 @@ SimStatus LAGRIDPlumeModel::runFullModel() { solarTime_h_ = ( timestepVars_.curr_Time_s + timestepVars_.dt / 2.0 ) / 3600.0; simTime_h_ = ( timestepVars_.curr_Time_s + timestepVars_.dt / 2.0 - timestepVars_.timeArray[0] ) / 3600.0; - // Run Ice Growth - if (simVars_.ICE_GROWTH && timestepVars_.checkTimeForIceGrowth()) { - std::cout << "Running ice growth..." << std::endl; - - #ifdef ENABLE_TIMING - auto icegrowth_start = std::chrono::high_resolution_clock::now(); - #endif - - timestepVars_.lastTimeIceGrowth = timestepVars_.curr_Time_s + timestepVars_.dt; - - // Substep ice growth in increments of <= dt_micro_target seconds (default 60s) - const double dt_micro_target = (optInput_.AEROSOL_ICE_GROWTH_SUBSTEP > 0.0) - ? optInput_.AEROSOL_ICE_GROWTH_SUBSTEP : 60.0; - const double dt_growth_total = timestepVars_.ICE_GROWTH_DT; - const int n_growth_substeps = std::max(1, static_cast(std::ceil(dt_growth_total / dt_micro_target))); - const double dt_growth_sub = dt_growth_total / n_growth_substeps; - - for (int sub = 0; sub < n_growth_substeps; ++sub) { - iceAerosol_.Grow(dt_growth_sub, H2O_, met_.Temp(), met_.Press()); - } - - #ifdef ENABLE_TIMING - auto icegrowth_end = std::chrono::high_resolution_clock::now(); - auto icegrowth_duration = std::chrono::duration_cast(icegrowth_end - icegrowth_start); - std::cout << " Ran ice growth in " << icegrowth_duration.count() << " ms" << std::endl; - #endif - } - #ifdef ENABLE_TIMING auto tracer_start = std::chrono::high_resolution_clock::now(); #endif @@ -430,11 +423,10 @@ void LAGRIDPlumeModel::initH2O() { } } -void LAGRIDPlumeModel::updateDiffVecs() { +void LAGRIDPlumeModel::updateDiffVecs(double time_start, double dt) { double dh_enhanced, dv_enhanced; - // Update Diffusion with exact time-averaged diffusivity over the transport step - const double time_start = timestepVars_.curr_Time_s - timestepVars_.tInitial_s; - PlumeModelUtils::DiffParam( time_start, timestepVars_.TRANSPORT_DT, + // Update Diffusion with exact time-averaged diffusivity over the subcycle transport step + PlumeModelUtils::DiffParam( time_start, dt, dh_enhanced, dv_enhanced, input_.horizDiff(), input_.vertiDiff() ); auto number = iceAerosol_.TotalNumber(); auto num_max = VectorUtils::max(number); @@ -450,7 +442,7 @@ void LAGRIDPlumeModel::updateDiffVecs() { } } } -void LAGRIDPlumeModel::runTransport(double timestep) { +void LAGRIDPlumeModel::runTransport(double timestep, double time_start) { #ifdef ENABLE_TIMING auto start = std::chrono::high_resolution_clock::now(); @@ -466,7 +458,7 @@ void LAGRIDPlumeModel::runTransport(double timestep) { } shear_rep_ = met_.shear(maxIdx); - const FVM_ANDS::AdvDiffParams fvmSolverInitParams(0, 0, shear_rep_, input_.horizDiff(), input_.vertiDiff(), timestepVars_.TRANSPORT_DT); + const FVM_ANDS::AdvDiffParams fvmSolverInitParams(0, 0, shear_rep_, input_.horizDiff(), input_.vertiDiff(), timestep); const FVM_ANDS::BoundaryConditions ZERO_BC_INIT = FVM_ANDS::bcFrom2DVector(iceAerosol_.getPDF()[0], true); #ifdef ENABLE_TIMING @@ -477,7 +469,7 @@ void LAGRIDPlumeModel::runTransport(double timestep) { start = std::chrono::high_resolution_clock::now(); #endif - updateDiffVecs(); + updateDiffVecs(time_start, timestep); #ifdef ENABLE_TIMING end = std::chrono::high_resolution_clock::now(); diff --git a/Code.v05-00/src/YamlInputReader/YamlInputReader.cpp b/Code.v05-00/src/YamlInputReader/YamlInputReader.cpp index f68a956ea..f82c0161d 100644 --- a/Code.v05-00/src/YamlInputReader/YamlInputReader.cpp +++ b/Code.v05-00/src/YamlInputReader/YamlInputReader.cpp @@ -377,6 +377,11 @@ namespace YamlInputReader{ input.TRANSPORT_TRANSPORT = parseBoolString(transportNode["Turn on Transport (T/F)"].as(), "Turn on Transport (T/F)"); input.TRANSPORT_FILL = parseBoolString(transportNode["Fill Negative Values (T/F)"].as(), "Fill Negative Values (T/F)"); input.TRANSPORT_TIMESTEP = parseDoubleString(transportNode["Transport Timestep [min] (double)"].as(), "Transport Timestep [min] (double)"); + if (transportNode["Transport and ice growth substep [s] (double)"]) { + input.TRANSPORT_ICE_GROWTH_SUBSTEP = parseDoubleString(transportNode["Transport and ice growth substep [s] (double)"].as(), "Transport and ice growth substep [s] (double)"); + } else { + input.TRANSPORT_ICE_GROWTH_SUBSTEP = 60.0; + } YAML::Node updraftSubmenu = transportNode["PLUME UPDRAFT SUBMENU"]; input.TRANSPORT_UPDRAFT = parseBoolString(updraftSubmenu["Turn on plume updraft (T/F)"].as(), "Turn on plume updraft (T/F)"); @@ -394,11 +399,6 @@ namespace YamlInputReader{ input.AEROSOL_COAGULATION_LIQUID = parseBoolString(aeroNode["Turn on liquid coagulation (T/F)"].as(), "Turn on liquid coagulation (T/F)"); input.AEROSOL_ICE_GROWTH = parseBoolString(aeroNode["Turn on ice growth (T/F)"].as(), "Turn on ice growth (T/F)"); input.AEROSOL_ICE_GROWTH_TIMESTEP = parseDoubleString(aeroNode["Ice growth timestep [min] (double)"].as(), "Ice growth timestep [min] (double)"); - if (aeroNode["Ice growth substep [s] (double)"]) { - input.AEROSOL_ICE_GROWTH_SUBSTEP = parseDoubleString(aeroNode["Ice growth substep [s] (double)"].as(), "Ice growth substep [s] (double)"); - } else { - input.AEROSOL_ICE_GROWTH_SUBSTEP = 60.0; - } } void readMetMenu(OptInput& input, const YAML::Node& metNode){ YAML::Node metInputSubmenu = metNode["METEOROLOGICAL INPUT SUBMENU"]; diff --git a/Code.v05-00/tests/test_yamlreader.cpp b/Code.v05-00/tests/test_yamlreader.cpp index ba54c4890..1c837973d 100644 --- a/Code.v05-00/tests/test_yamlreader.cpp +++ b/Code.v05-00/tests/test_yamlreader.cpp @@ -250,6 +250,7 @@ TEST_CASE("Read Yaml File"){ REQUIRE(input.TRANSPORT_TRANSPORT == true); REQUIRE(input.TRANSPORT_FILL == true); REQUIRE(input.TRANSPORT_TIMESTEP == 10); + REQUIRE(input.TRANSPORT_ICE_GROWTH_SUBSTEP == 60.0); REQUIRE(input.TRANSPORT_UPDRAFT == true); REQUIRE(input.TRANSPORT_UPDRAFT_TIMESCALE == 3600); REQUIRE(input.TRANSPORT_UPDRAFT_VELOCITY == 5); @@ -269,7 +270,6 @@ TEST_CASE("Read Yaml File"){ REQUIRE(input.AEROSOL_COAGULATION_LIQUID == true); REQUIRE(input.AEROSOL_ICE_GROWTH == true); REQUIRE(input.AEROSOL_ICE_GROWTH_TIMESTEP == 10); - REQUIRE(input.AEROSOL_ICE_GROWTH_SUBSTEP == 60.0); } SECTION("Read Met Menu"){ OptInput input; From 19758197f4ca79eb96020ff634aae62d24078262 Mon Sep 17 00:00:00 2001 From: "Sebastian D. Eastham" Date: Wed, 19 Aug 2026 10:17:33 -0400 Subject: [PATCH 10/14] Fixed and simplified subcycling --- Code.v05-00/defaults/input.yaml | 5 +- .../include/Core/TimestepVarsWrapper.hpp | 5 +- Code.v05-00/include/Defaults/Input.hpp | 5 +- Code.v05-00/src/Core/LAGRIDPlumeModel.cpp | 25 +++------ Code.v05-00/src/Core/TimestepVarsWrapper.cpp | 25 ++------- .../src/YamlInputReader/YamlInputReader.cpp | 52 +++++++++++++++++-- Code.v05-00/tests/test_yamlreader.cpp | 31 ++++++++++- 7 files changed, 96 insertions(+), 52 deletions(-) diff --git a/Code.v05-00/defaults/input.yaml b/Code.v05-00/defaults/input.yaml index 3d25587d9..c03496f7c 100644 --- a/Code.v05-00/defaults/input.yaml +++ b/Code.v05-00/defaults/input.yaml @@ -84,8 +84,8 @@ TRANSPORT MENU: Turn on Transport (T/F): T # Outdated, not used (was used by spectral solver) Fill Negative Values (T/F): T - Transport Timestep [min] (double): 1 - Transport and ice growth substep [s] (double): 60.0 + Outer time step [min] (double): 1 + Inner physics time step [s] (double): 60.0 # Keep off: not sure of the effect yet + met updraft is included (if met file input) PLUME UPDRAFT SUBMENU: Turn on plume updraft (T/F): F @@ -107,7 +107,6 @@ AEROSOL MENU: Turn on liquid coagulation (T/F): F # Keep on Turn on ice growth (T/F): T - Ice growth timestep [min] (double): 1 # At least one of "Use met. input", "Impose moist layer depth", or "Impose lapse rate" must be true # Imposing moist layer depth will automatically calculate the lapse rate and override the imposed lapse rate diff --git a/Code.v05-00/include/Core/TimestepVarsWrapper.hpp b/Code.v05-00/include/Core/TimestepVarsWrapper.hpp index b3ba7a3c3..34612cb35 100644 --- a/Code.v05-00/include/Core/TimestepVarsWrapper.hpp +++ b/Code.v05-00/include/Core/TimestepVarsWrapper.hpp @@ -74,13 +74,12 @@ struct TimestepVarsWrapper } inline bool checkTimeForTransport() { - ITS_TIME_FOR_TRANSPORT = (((curr_Time_s + dt - lastTimeTransport) >= TRANSPORT_DT) || LAST_STEP); + ITS_TIME_FOR_TRANSPORT = (((curr_Time_s - lastTimeTransport) >= -1.0e-6) || LAST_STEP); return ITS_TIME_FOR_TRANSPORT; } inline bool checkTimeForIceGrowth() { - /* TODO: For now perform growth at every time step */ - ITS_TIME_FOR_ICE_GROWTH = (((curr_Time_s + dt - lastTimeIceGrowth) >= ICE_GROWTH_DT) || LAST_STEP); + ITS_TIME_FOR_ICE_GROWTH = (((curr_Time_s - lastTimeIceGrowth) >= -1.0e-6) || LAST_STEP); return ITS_TIME_FOR_ICE_GROWTH; } }; diff --git a/Code.v05-00/include/Defaults/Input.hpp b/Code.v05-00/include/Defaults/Input.hpp index ce9c151e7..cf531c0ad 100644 --- a/Code.v05-00/include/Defaults/Input.hpp +++ b/Code.v05-00/include/Defaults/Input.hpp @@ -85,8 +85,8 @@ TRANSPORT MENU: Turn on Transport (T/F): T # Outdated, not used (was used by spectral solver) Fill Negative Values (T/F): T - Transport Timestep [min] (double): 1 - Transport and ice growth substep [s] (double): 60.0 + Outer time step [min] (double): 1 + Inner physics time step [s] (double): 60.0 # Keep off: not sure of the effect yet + met updraft is included (if met file input) PLUME UPDRAFT SUBMENU: Turn on plume updraft (T/F): F @@ -108,7 +108,6 @@ AEROSOL MENU: Turn on liquid coagulation (T/F): F # Keep on Turn on ice growth (T/F): T - Ice growth timestep [min] (double): 1 # At least one of "Use met. input", "Impose moist layer depth", or "Impose lapse rate" must be true # Imposing moist layer depth will automatically calculate the lapse rate and override the imposed lapse rate diff --git a/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp b/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp index c4726f3d4..a8655d7b4 100644 --- a/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp +++ b/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp @@ -104,24 +104,15 @@ SimStatus LAGRIDPlumeModel::runFullModel() { std::cout << "\n - Time step: " << timestepVars_.nTime + 1 << " out of " << timestepVars_.timeArray.size(); std::cout << "\n -> Solar time: " << std::fmod( timestepVars_.curr_Time_s/3600.0, 24.0 ) << " [hr]" << std::endl; - // Interleaved Transport and Ice Growth Subcycling - bool timeForTransport = (simVars_.TRANSPORT && (timestepVars_.nTime == 0 || timestepVars_.checkTimeForTransport())); - bool timeForIceGrowth = (simVars_.ICE_GROWTH && timestepVars_.checkTimeForIceGrowth()); - - if (timeForTransport || timeForIceGrowth) { - double dt_total = timeForTransport ? timestepVars_.TRANSPORT_DT : timestepVars_.ICE_GROWTH_DT; + // Interleaved Transport and Ice Growth Subcycling over the outer timestep dt + if (simVars_.TRANSPORT || simVars_.ICE_GROWTH) { + const double dt_step = timestepVars_.dt; const double dt_sub_target = (optInput_.TRANSPORT_ICE_GROWTH_SUBSTEP > 0.0) ? optInput_.TRANSPORT_ICE_GROWTH_SUBSTEP : 60.0; - const int n_subcycles = std::max(1, static_cast(std::ceil(dt_total / dt_sub_target))); - const double dt_sub = dt_total / n_subcycles; + const int n_subcycles = std::max(1, static_cast(std::ceil(dt_step / dt_sub_target))); + const double dt_sub = dt_step / n_subcycles; - if (timeForTransport) { - std::cout << "Running Transport and ice growth subcycling (" << n_subcycles << " x " << dt_sub << " s)..." << std::endl; - timestepVars_.lastTimeTransport = timestepVars_.curr_Time_s + timestepVars_.dt; - } - if (timeForIceGrowth) { - timestepVars_.lastTimeIceGrowth = timestepVars_.curr_Time_s + timestepVars_.dt; - } + std::cout << "Running Transport and ice growth subcycling (" << n_subcycles << " x " << dt_sub << " s)..." << std::endl; #ifdef ENABLE_TIMING auto subcycling_start = std::chrono::high_resolution_clock::now(); @@ -129,10 +120,10 @@ SimStatus LAGRIDPlumeModel::runFullModel() { for (int sub = 0; sub < n_subcycles; ++sub) { const double t_sub_start = (timestepVars_.curr_Time_s - timestepVars_.tInitial_s) + sub * dt_sub; - if (timeForTransport) { + if (simVars_.TRANSPORT) { runTransport(dt_sub, t_sub_start); } - if (timeForIceGrowth) { + if (simVars_.ICE_GROWTH) { iceAerosol_.Grow(dt_sub, H2O_, met_.Temp(), met_.Press()); } } diff --git a/Code.v05-00/src/Core/TimestepVarsWrapper.cpp b/Code.v05-00/src/Core/TimestepVarsWrapper.cpp index 0d3c31b3c..780738e6f 100644 --- a/Code.v05-00/src/Core/TimestepVarsWrapper.cpp +++ b/Code.v05-00/src/Core/TimestepVarsWrapper.cpp @@ -33,27 +33,8 @@ totPart_lost(0), totIce_lost(0), ABORT_THRESHOLD(1.0e-3) { - /* The base (heartbeat) timestep is determined from enabled transport and ice growth processes. */ - - Vector_1D timesteps; - - if(Input_Opt.TRANSPORT_TRANSPORT) - timesteps.push_back(TRANSPORT_DT); - if(Input_Opt.AEROSOL_ICE_GROWTH) - timesteps.push_back(ICE_GROWTH_DT); - - if (timesteps.empty()) { - throw std::runtime_error("No active processes to determine timestep"); - } - - dt = *(std::min_element(timesteps.begin(), timesteps.end())); - std::cout << "Calculated Timestep: " << dt/60.0 << "[min]" << std::endl; + dt = TRANSPORT_DT; if (dt <= 0) throw std::runtime_error("Invalid Timestep"); - - for (const double &step : timesteps) { - double ratio = step / dt; - if (std::abs(ratio - std::round(ratio)) > 1.0e-6) { - throw std::runtime_error("Process timesteps must be integer multiples of the heartbeat timestep"); - } - } + std::cout << "Outer Timestep: " << dt/60.0 << " [min]" << std::endl; + std::cout << "Inner Physics Substep: " << Input_Opt.TRANSPORT_ICE_GROWTH_SUBSTEP << " [s]" << std::endl; } \ No newline at end of file diff --git a/Code.v05-00/src/YamlInputReader/YamlInputReader.cpp b/Code.v05-00/src/YamlInputReader/YamlInputReader.cpp index f82c0161d..c311db539 100644 --- a/Code.v05-00/src/YamlInputReader/YamlInputReader.cpp +++ b/Code.v05-00/src/YamlInputReader/YamlInputReader.cpp @@ -87,6 +87,7 @@ namespace YamlInputReader{ "Chemistry Timestep [min] (double)", "Coag. timestep [min] (double)", "Temp. Perturb. Timescale (min)", + "Ice growth timestep [min] (double)", }; if (deprecatedKeys.contains(key)) { std::cout << "WARNING: Deprecated option found: '" << errorPath << "'. This option is no longer used and has no effect." << std::endl; @@ -95,6 +96,14 @@ namespace YamlInputReader{ return false; } + bool checkAliasedKey(const std::string& key) { + static const std::set aliasedKeys = { + "Transport Timestep [min] (double)", + "Transport and ice growth substep [s] (double)", + }; + return aliasedKeys.contains(key); + } + // Keys that previous versions accepted and that we now reject. Checked // before the generic unknown-key error so an outdated input file gets a message // naming the option that went away instead of "Unknown key found". @@ -144,6 +153,9 @@ namespace YamlInputReader{ for (const auto& key : userKeys) { if (!defaultKeys.contains(key)) { + if (checkAliasedKey(key)) { + continue; + } // The key from the user's YAML does not exist in the default YAML. std::string errorPath = currentPath.empty() ? key : currentPath + " -> " + key; if (checkDeprecatedKey(key, errorPath)) { @@ -186,6 +198,21 @@ namespace YamlInputReader{ throw std::runtime_error("Invalid field in YAML input file '" + filename + "': " + e.what()); } INPUT_FILE_PATH = std::filesystem::path(filename); + + // Map legacy aliases if present + if (userData["TRANSPORT MENU"]) { + if (userData["TRANSPORT MENU"]["Transport Timestep [min] (double)"] && + !userData["TRANSPORT MENU"]["Outer time step [min] (double)"]) { + userData["TRANSPORT MENU"]["Outer time step [min] (double)"] = + userData["TRANSPORT MENU"]["Transport Timestep [min] (double)"]; + } + if (userData["TRANSPORT MENU"]["Transport and ice growth substep [s] (double)"] && + !userData["TRANSPORT MENU"]["Inner physics time step [s] (double)"]) { + userData["TRANSPORT MENU"]["Inner physics time step [s] (double)"] = + userData["TRANSPORT MENU"]["Transport and ice growth substep [s] (double)"]; + } + } + mergedData = mergeYamlNodes(mergedData, userData); } @@ -376,13 +403,30 @@ namespace YamlInputReader{ void readTransportMenu(OptInput& input, const YAML::Node& transportNode){ input.TRANSPORT_TRANSPORT = parseBoolString(transportNode["Turn on Transport (T/F)"].as(), "Turn on Transport (T/F)"); input.TRANSPORT_FILL = parseBoolString(transportNode["Fill Negative Values (T/F)"].as(), "Fill Negative Values (T/F)"); - input.TRANSPORT_TIMESTEP = parseDoubleString(transportNode["Transport Timestep [min] (double)"].as(), "Transport Timestep [min] (double)"); - if (transportNode["Transport and ice growth substep [s] (double)"]) { + + if (transportNode["Outer time step [min] (double)"]) { + input.TRANSPORT_TIMESTEP = parseDoubleString(transportNode["Outer time step [min] (double)"].as(), "Outer time step [min] (double)"); + } else if (transportNode["Transport Timestep [min] (double)"]) { + input.TRANSPORT_TIMESTEP = parseDoubleString(transportNode["Transport Timestep [min] (double)"].as(), "Transport Timestep [min] (double)"); + } else { + throw std::invalid_argument("In YamlInputReader::readTransportMenu: Missing 'Outer time step [min] (double)'"); + } + + if (transportNode["Inner physics time step [s] (double)"]) { + input.TRANSPORT_ICE_GROWTH_SUBSTEP = parseDoubleString(transportNode["Inner physics time step [s] (double)"].as(), "Inner physics time step [s] (double)"); + } else if (transportNode["Transport and ice growth substep [s] (double)"]) { input.TRANSPORT_ICE_GROWTH_SUBSTEP = parseDoubleString(transportNode["Transport and ice growth substep [s] (double)"].as(), "Transport and ice growth substep [s] (double)"); } else { input.TRANSPORT_ICE_GROWTH_SUBSTEP = 60.0; } + if (input.TRANSPORT_ICE_GROWTH_SUBSTEP > input.TRANSPORT_TIMESTEP * 60.0) { + throw std::invalid_argument("In YamlInputReader::readTransportMenu: 'Inner physics time step [s]' (" + + std::to_string(input.TRANSPORT_ICE_GROWTH_SUBSTEP) + + " s) cannot be greater than 'Outer time step [min]' (" + + std::to_string(input.TRANSPORT_TIMESTEP * 60.0) + " s)"); + } + YAML::Node updraftSubmenu = transportNode["PLUME UPDRAFT SUBMENU"]; input.TRANSPORT_UPDRAFT = parseBoolString(updraftSubmenu["Turn on plume updraft (T/F)"].as(), "Turn on plume updraft (T/F)"); input.TRANSPORT_UPDRAFT_TIMESCALE = parseDoubleString(updraftSubmenu["Updraft timescale [s] (double)"].as(), "Updraft timescale [s] (double)"); @@ -398,7 +442,9 @@ namespace YamlInputReader{ input.AEROSOL_COAGULATION_SOLID = parseBoolString(aeroNode["Turn on solid coagulation (T/F)"].as(), "Turn on solid coagulation (T/F)"); input.AEROSOL_COAGULATION_LIQUID = parseBoolString(aeroNode["Turn on liquid coagulation (T/F)"].as(), "Turn on liquid coagulation (T/F)"); input.AEROSOL_ICE_GROWTH = parseBoolString(aeroNode["Turn on ice growth (T/F)"].as(), "Turn on ice growth (T/F)"); - input.AEROSOL_ICE_GROWTH_TIMESTEP = parseDoubleString(aeroNode["Ice growth timestep [min] (double)"].as(), "Ice growth timestep [min] (double)"); + if (aeroNode["Ice growth timestep [min] (double)"]) { + std::cout << "WARNING: 'Ice growth timestep [min] (double)' in input.yaml is deprecated and ignored. Time stepping is now controlled by 'Outer time step [min]' and 'Inner physics time step [s]' under TRANSPORT MENU." << std::endl; + } } void readMetMenu(OptInput& input, const YAML::Node& metNode){ YAML::Node metInputSubmenu = metNode["METEOROLOGICAL INPUT SUBMENU"]; diff --git a/Code.v05-00/tests/test_yamlreader.cpp b/Code.v05-00/tests/test_yamlreader.cpp index 1c837973d..a429f3a62 100644 --- a/Code.v05-00/tests/test_yamlreader.cpp +++ b/Code.v05-00/tests/test_yamlreader.cpp @@ -254,6 +254,36 @@ TEST_CASE("Read Yaml File"){ REQUIRE(input.TRANSPORT_UPDRAFT == true); REQUIRE(input.TRANSPORT_UPDRAFT_TIMESCALE == 3600); REQUIRE(input.TRANSPORT_UPDRAFT_VELOCITY == 5); + + // Test with new "Outer time step [min]" and "Inner physics time step [s]" keys + YAML::Node customTransport = YAML::Load( + "Turn on Transport (T/F): T\n" + "Fill Negative Values (T/F): T\n" + "Outer time step [min] (double): 5\n" + "Inner physics time step [s] (double): 30.0\n" + "PLUME UPDRAFT SUBMENU:\n" + " Turn on plume updraft (T/F): F\n" + " Updraft timescale [s] (double): 3600\n" + " Updraft veloc. [cm/s] (double): 5\n" + ); + OptInput customInput; + readTransportMenu(customInput, customTransport); + REQUIRE(customInput.TRANSPORT_TIMESTEP == 5.0); + REQUIRE(customInput.TRANSPORT_ICE_GROWTH_SUBSTEP == 30.0); + + // Test validation: Inner substep > Outer step throws + YAML::Node invalidTransport = YAML::Load( + "Turn on Transport (T/F): T\n" + "Fill Negative Values (T/F): T\n" + "Outer time step [min] (double): 1\n" + "Inner physics time step [s] (double): 120.0\n" + "PLUME UPDRAFT SUBMENU:\n" + " Turn on plume updraft (T/F): F\n" + " Updraft timescale [s] (double): 3600\n" + " Updraft veloc. [cm/s] (double): 5\n" + ); + OptInput invalidInput; + REQUIRE_THROWS_AS(readTransportMenu(invalidInput, invalidTransport), std::invalid_argument); } SECTION("Read Chemistry Menu"){ OptInput input; @@ -269,7 +299,6 @@ TEST_CASE("Read Yaml File"){ REQUIRE(input.AEROSOL_COAGULATION_SOLID == true); REQUIRE(input.AEROSOL_COAGULATION_LIQUID == true); REQUIRE(input.AEROSOL_ICE_GROWTH == true); - REQUIRE(input.AEROSOL_ICE_GROWTH_TIMESTEP == 10); } SECTION("Read Met Menu"){ OptInput input; From a0e3eff7ab9d2ea7f28110ffbeed0bb7b1e9104f Mon Sep 17 00:00:00 2001 From: "Sebastian D. Eastham" Date: Wed, 19 Aug 2026 17:18:31 -0400 Subject: [PATCH 11/14] First attempt at flexible grid resolution --- Code.v05-00/defaults/input.yaml | 5 +++ Code.v05-00/include/Core/Input_Mod.hpp | 5 +++ Code.v05-00/include/Defaults/Input.hpp | 5 +++ Code.v05-00/src/Core/LAGRIDPlumeModel.cpp | 24 +++++++++--- .../src/YamlInputReader/YamlInputReader.cpp | 37 +++++++++++++++++++ 5 files changed, 70 insertions(+), 6 deletions(-) diff --git a/Code.v05-00/defaults/input.yaml b/Code.v05-00/defaults/input.yaml index c03496f7c..f2e5af07b 100644 --- a/Code.v05-00/defaults/input.yaml +++ b/Code.v05-00/defaults/input.yaml @@ -182,6 +182,11 @@ ADVANCED OPTIONS MENU: XLIM_LEFT (positive double): 1.0e+3 YLIM_UP (positive double): 300 YLIM_DOWN (positive double): 1.5e+3 + Target points in plume [-] (int): 50 + Min DX [m] (double): 20.0 + Max DX [m] (double): 50.0 + Min DY [m] (double): 5.0 + Max DY [m] (double): 7.0 INITIAL CONTRAIL SIZE SUBMENU: #Depth = BaseDepth + DepthScalingFactor * Default_Depth #Same formula for width diff --git a/Code.v05-00/include/Core/Input_Mod.hpp b/Code.v05-00/include/Core/Input_Mod.hpp index 921507360..d6e9b1b51 100644 --- a/Code.v05-00/include/Core/Input_Mod.hpp +++ b/Code.v05-00/include/Core/Input_Mod.hpp @@ -139,6 +139,11 @@ struct OptInput double ADV_GRID_XLIM_LEFT; double ADV_GRID_YLIM_UP; double ADV_GRID_YLIM_DOWN; + unsigned int ADV_GRID_TARGET_PLUME_PTS; + double ADV_GRID_MIN_DX; + double ADV_GRID_MAX_DX; + double ADV_GRID_MIN_DY; + double ADV_GRID_MAX_DY; double ADV_CSIZE_DEPTH_BASE; double ADV_CSIZE_DEPTH_SCALING_FACTOR; double ADV_CSIZE_WIDTH_BASE; diff --git a/Code.v05-00/include/Defaults/Input.hpp b/Code.v05-00/include/Defaults/Input.hpp index cf531c0ad..e205fa044 100644 --- a/Code.v05-00/include/Defaults/Input.hpp +++ b/Code.v05-00/include/Defaults/Input.hpp @@ -183,6 +183,11 @@ ADVANCED OPTIONS MENU: XLIM_LEFT (positive double): 1.0e+3 YLIM_UP (positive double): 300 YLIM_DOWN (positive double): 1.5e+3 + Target points in plume [-] (int): 50 + Min DX [m] (double): 20.0 + Max DX [m] (double): 50.0 + Min DY [m] (double): 5.0 + Max DY [m] (double): 7.0 INITIAL CONTRAIL SIZE SUBMENU: #Depth = BaseDepth + DepthScalingFactor * Default_Depth #Same formula for width diff --git a/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp b/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp index a8655d7b4..9b4d9639f 100644 --- a/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp +++ b/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp @@ -651,9 +651,15 @@ std::pair LAGRIDPlumeModel::r double dx_grid_old = xCoords_[1] - xCoords_[0]; auto boxGrid = LAGRID::rectToBoxGrid(met_.dy_vec(), dx_grid_old, xEdges_[0], yEdges_[0], phi, mask); - //Enforce at least x many points in the contrail while limiting minimum/maximum dx and dy - double dx_grid_new = std::max(20.0, std::min((maskInfo.maxX - maskInfo.minX) / 50.0, 50.0)); - double dy_grid_new = std::max(5.0, std::min((maskInfo.maxY - maskInfo.minY) / 50.0, 7.0)); + //Enforce at least target points in the contrail while limiting minimum/maximum dx and dy + const double target_pts = (optInput_.ADV_GRID_TARGET_PLUME_PTS > 0) ? static_cast(optInput_.ADV_GRID_TARGET_PLUME_PTS) : 50.0; + const double min_dx = (optInput_.ADV_GRID_MIN_DX > 0) ? optInput_.ADV_GRID_MIN_DX : 20.0; + const double max_dx = (optInput_.ADV_GRID_MAX_DX >= min_dx) ? optInput_.ADV_GRID_MAX_DX : 50.0; + const double min_dy = (optInput_.ADV_GRID_MIN_DY > 0) ? optInput_.ADV_GRID_MIN_DY : 5.0; + const double max_dy = (optInput_.ADV_GRID_MAX_DY >= min_dy) ? optInput_.ADV_GRID_MAX_DY : 7.0; + + double dx_grid_new = std::max(min_dx, std::min((maskInfo.maxX - maskInfo.minX) / target_pts, max_dx)); + double dy_grid_new = std::max(min_dy, std::min((maskInfo.maxY - maskInfo.minY) / target_pts, max_dy)); //Need 2 extra points account for the buffer int nx_new = floor((maskInfo.maxX - maskInfo.minX) / dx_grid_new) + 2; int ny_new = floor((maskInfo.maxY - maskInfo.minY) / dy_grid_new) + 2; @@ -811,9 +817,15 @@ Eigen::SparseMatrix LAGRIDPlumeModel::createRegriddingWeightsSparse(cons double dx_grid_old = xCoords_[1] - xCoords_[0]; auto boxGrid = LAGRID::rectToBoxGrid(met_.dy_vec(), dx_grid_old, xEdges_[0], yEdges_[0], phi, mask); - //Enforce at least x many points in the contrail while limiting minimum/maximum dx and dy - double dx_grid_new = std::max(20.0, std::min((maskInfo.maxX - maskInfo.minX) / 50.0, 50.0)); - double dy_grid_new = std::max(5.0, std::min((maskInfo.maxY - maskInfo.minY) / 50.0, 7.0)); + //Enforce at least target points in the contrail while limiting minimum/maximum dx and dy + const double target_pts = (optInput_.ADV_GRID_TARGET_PLUME_PTS > 0) ? static_cast(optInput_.ADV_GRID_TARGET_PLUME_PTS) : 50.0; + const double min_dx = (optInput_.ADV_GRID_MIN_DX > 0) ? optInput_.ADV_GRID_MIN_DX : 20.0; + const double max_dx = (optInput_.ADV_GRID_MAX_DX >= min_dx) ? optInput_.ADV_GRID_MAX_DX : 50.0; + const double min_dy = (optInput_.ADV_GRID_MIN_DY > 0) ? optInput_.ADV_GRID_MIN_DY : 5.0; + const double max_dy = (optInput_.ADV_GRID_MAX_DY >= min_dy) ? optInput_.ADV_GRID_MAX_DY : 7.0; + + double dx_grid_new = std::max(min_dx, std::min((maskInfo.maxX - maskInfo.minX) / target_pts, max_dx)); + double dy_grid_new = std::max(min_dy, std::min((maskInfo.maxY - maskInfo.minY) / target_pts, max_dy)); //Need 2 extra points account for the buffer int nx_new = floor((maskInfo.maxX - maskInfo.minX) / dx_grid_new) + 2; int ny_new = floor((maskInfo.maxY - maskInfo.minY) / dy_grid_new) + 2; diff --git a/Code.v05-00/src/YamlInputReader/YamlInputReader.cpp b/Code.v05-00/src/YamlInputReader/YamlInputReader.cpp index c311db539..cf0226001 100644 --- a/Code.v05-00/src/YamlInputReader/YamlInputReader.cpp +++ b/Code.v05-00/src/YamlInputReader/YamlInputReader.cpp @@ -496,6 +496,43 @@ namespace YamlInputReader{ input.ADV_GRID_XLIM_LEFT = parseDoubleString(gridSubmenu["XLIM_LEFT (positive double)"].as(), "XLIM_LEFT (positive double)"); input.ADV_GRID_YLIM_UP = parseDoubleString(gridSubmenu["YLIM_UP (positive double)"].as(), "YLIM_UP (positive double)"); input.ADV_GRID_YLIM_DOWN = parseDoubleString(gridSubmenu["YLIM_DOWN (positive double)"].as(), "YLIM_DOWN (positive double)"); + + if (gridSubmenu["Target points in plume [-] (int)"]) { + input.ADV_GRID_TARGET_PLUME_PTS = parseUIntString(gridSubmenu["Target points in plume [-] (int)"].as(), "Target points in plume [-] (int)"); + } else { + input.ADV_GRID_TARGET_PLUME_PTS = 50; + } + + if (gridSubmenu["Min DX [m] (double)"]) { + input.ADV_GRID_MIN_DX = parseDoubleString(gridSubmenu["Min DX [m] (double)"].as(), "Min DX [m] (double)"); + } else { + input.ADV_GRID_MIN_DX = 20.0; + } + + if (gridSubmenu["Max DX [m] (double)"]) { + input.ADV_GRID_MAX_DX = parseDoubleString(gridSubmenu["Max DX [m] (double)"].as(), "Max DX [m] (double)"); + } else { + input.ADV_GRID_MAX_DX = 50.0; + } + + if (gridSubmenu["Min DY [m] (double)"]) { + input.ADV_GRID_MIN_DY = parseDoubleString(gridSubmenu["Min DY [m] (double)"].as(), "Min DY [m] (double)"); + } else { + input.ADV_GRID_MIN_DY = 5.0; + } + + if (gridSubmenu["Max DY [m] (double)"]) { + input.ADV_GRID_MAX_DY = parseDoubleString(gridSubmenu["Max DY [m] (double)"].as(), "Max DY [m] (double)"); + } else { + input.ADV_GRID_MAX_DY = 7.0; + } + + if (input.ADV_GRID_MIN_DX > input.ADV_GRID_MAX_DX) { + throw std::invalid_argument("Min DX [m] cannot be greater than Max DX [m] in GRID SUBMENU!"); + } + if (input.ADV_GRID_MIN_DY > input.ADV_GRID_MAX_DY) { + throw std::invalid_argument("Min DY [m] cannot be greater than Max DY [m] in GRID SUBMENU!"); + } YAML::Node csizeSubmenu = advancedNode["INITIAL CONTRAIL SIZE SUBMENU"]; input.ADV_CSIZE_DEPTH_BASE = parseDoubleString(csizeSubmenu["Base Contrail Depth [m] (double)"].as(), "Base Contrail Depth [m] (double)"); From 2e9724a418bf99fa047041c6b97883a825149b7b Mon Sep 17 00:00:00 2001 From: "Sebastian D. Eastham" Date: Sun, 23 Aug 2026 10:44:19 -0400 Subject: [PATCH 12/14] Aerosol size bin resolution now specified in input.yaml --- Code.v05-00/include/Core/Input_Mod.hpp | 3 + Code.v05-00/include/Defaults/Input.hpp | 4 + .../src/EPM/Models/Original/Integrate.cpp | 11 ++- Code.v05-00/src/EPM/Solution.cpp | 13 ++- .../src/YamlInputReader/YamlInputReader.cpp | 23 +++++ Code.v05-00/tests/test_aerosol.cpp | 96 +++++++++++++++++++ 6 files changed, 141 insertions(+), 9 deletions(-) diff --git a/Code.v05-00/include/Core/Input_Mod.hpp b/Code.v05-00/include/Core/Input_Mod.hpp index d6e9b1b51..f8f6d2f55 100644 --- a/Code.v05-00/include/Core/Input_Mod.hpp +++ b/Code.v05-00/include/Core/Input_Mod.hpp @@ -155,6 +155,9 @@ struct OptInput bool ADV_EP_N_POSTJET_OVERRIDE; double ADV_EP_N_POSTJET; bool ADV_SAVE_PSD_GRID; + double ADV_AERO_ICE_BIN_VRAT = 1.80; + double ADV_AERO_ICE_BIN_R_LOW = 5.00E-08; + double ADV_AERO_ICE_BIN_R_HIG = 8.00E-05; }; diff --git a/Code.v05-00/include/Defaults/Input.hpp b/Code.v05-00/include/Defaults/Input.hpp index e205fa044..30a2edb6a 100644 --- a/Code.v05-00/include/Defaults/Input.hpp +++ b/Code.v05-00/include/Defaults/Input.hpp @@ -197,6 +197,10 @@ ADVANCED OPTIONS MENU: Contrail Width Scaling Factor [-] (double): 1.0 Ambient Lapse Rate [K/km] (double): -3.0 Tropopause Pressure [Pa] (double): 2.0e+4 + AEROSOL GRID SUBMENU: + Ice bin volume ratio [-] (double): 1.80 + Ice bin r_min [m] (double): 5.0e-08 + Ice bin r_max [m] (double): 8.0e-05 EARLY PLUME SUBMENU: Reference ice crystal count [#/m] (double): 3.38e12 Reference wingspan [m] (double): 60.3 diff --git a/Code.v05-00/src/EPM/Models/Original/Integrate.cpp b/Code.v05-00/src/EPM/Models/Original/Integrate.cpp index d04ca4415..1268b95c3 100644 --- a/Code.v05-00/src/EPM/Models/Original/Integrate.cpp +++ b/Code.v05-00/src/EPM/Models/Original/Integrate.cpp @@ -125,16 +125,19 @@ namespace EPM::Models /* Number of ice size distribution bins based on specified min and max radii * * and volume ratio between two adjacent bins */ + const double ice_vrat = optInput_.ADV_AERO_ICE_BIN_VRAT; + const double ice_r_low = optInput_.ADV_AERO_ICE_BIN_R_LOW; + const double ice_r_hig = optInput_.ADV_AERO_ICE_BIN_R_HIG; const UInt Ice_NBIN = static_cast( - std::floor(1 + 3.0*log(PA_R_HIG / PA_R_LOW) / log(PA_VRAT))); + std::floor(1 + 3.0*log(ice_r_hig / ice_r_low) / log(ice_vrat))); /* Adjacent bin radius ratio */ - const double PA_RRAT = cbrt( PA_VRAT ); + const double PA_RRAT = cbrt( ice_vrat ); /* Ice bin center and edge radii */ Vector_1D Ice_rJ( Ice_NBIN , 0.0 ); Vector_1D Ice_rE( Ice_NBIN + 1, 0.0 ); - Ice_rE[0] = PA_R_LOW; + Ice_rE[0] = ice_r_low; for ( UInt iBin = 1; iBin < Ice_NBIN + 1; iBin++ ) Ice_rE[iBin] = Ice_rE[iBin-1] * PA_RRAT; /* [m] */ @@ -401,7 +404,7 @@ namespace EPM::Models SO4Aer = pSO4pdf_3mins; const double expsIce = 1.15; - AIM::Aerosol solidAer( Ice_rJ, Ice_rE, Ice_den, std::max(Ice_rad * exp ( -2.5 * log(expsIce) * log(expsIce) ), 1.5 * PA_R_LOW ), expsIce, "lognormal" ); + AIM::Aerosol solidAer( Ice_rJ, Ice_rE, Ice_den, std::max(Ice_rad * exp ( -2.5 * log(expsIce) * log(expsIce) ), 1.5 * optInput_.ADV_AERO_ICE_BIN_R_LOW ), expsIce, "lognormal" ); IceAer = solidAer; /* Compute plume area */ diff --git a/Code.v05-00/src/EPM/Solution.cpp b/Code.v05-00/src/EPM/Solution.cpp index b0cde6ecc..3f0644b5c 100644 --- a/Code.v05-00/src/EPM/Solution.cpp +++ b/Code.v05-00/src/EPM/Solution.cpp @@ -143,15 +143,18 @@ void Solution::Initialize(std::string fileName, nBin_LA = 2; //dumb hardcoded Grid_Aerosol default constructor } - double pa_r_hig_low = PA_R_HIG/PA_R_LOW; - nBin_PA = std::floor( 1 + log( pa_r_hig_low * pa_r_hig_low* pa_r_hig_low ) / log( PA_VRAT ) ); + const double pa_vrat = Input_Opt.ADV_AERO_ICE_BIN_VRAT; + const double pa_r_low = Input_Opt.ADV_AERO_ICE_BIN_R_LOW; + const double pa_r_hig = Input_Opt.ADV_AERO_ICE_BIN_R_HIG; + double pa_r_hig_low = pa_r_hig / pa_r_low; + nBin_PA = std::floor( 1 + log( pa_r_hig_low * pa_r_hig_low* pa_r_hig_low ) / log( pa_vrat ) ); Vector_1D PA_rE( nBin_PA + 1, 0.0 ); /* Bin edges in m */ Vector_1D PA_rJ( nBin_PA , 0.0 ); /* Bin center radius in m */ Vector_1D PA_vJ( nBin_PA , 0.0 ); /* Bin volume centers in m^3 */ - const double PA_RRAT = cbrt( PA_VRAT ); - PA_rE[0] = PA_R_LOW; + const double PA_RRAT = cbrt( pa_vrat ); + PA_rE[0] = pa_r_low; for ( UInt iBin_PA = 1; iBin_PA < nBin_PA + 1; iBin_PA++ ) PA_rE[iBin_PA] = PA_rE[iBin_PA-1] * PA_RRAT; /* [m] */ @@ -169,7 +172,7 @@ void Solution::Initialize(std::string fileName, if ( PA_nDens >= 0.0E+00 ) { const double expsPA = 1.6; - const double rPA = std::max( RAD[0] * exp( - 2.5 * log(expsPA) * log(expsPA) ), 1.5 * PA_R_LOW ); + const double rPA = std::max( RAD[0] * exp( - 2.5 * log(expsPA) * log(expsPA) ), 1.5 * pa_r_low ); AIM::Grid_Aerosol PAAerosol( size_x, size_y, PA_rJ, PA_rE, PA_nDens, rPA, expsPA, "lognormal" ); solidAerosol = PAAerosol; diff --git a/Code.v05-00/src/YamlInputReader/YamlInputReader.cpp b/Code.v05-00/src/YamlInputReader/YamlInputReader.cpp index cf0226001..27bfd42e6 100644 --- a/Code.v05-00/src/YamlInputReader/YamlInputReader.cpp +++ b/Code.v05-00/src/YamlInputReader/YamlInputReader.cpp @@ -556,6 +556,29 @@ namespace YamlInputReader{ input.ADV_EP_N_POSTJET_OVERRIDE = parseBoolString(earlyPlumeSubmenu["Override post-jet ice crystal count (T/F)"].as(), "Override post-jet ice crystal count (T/F)"); input.ADV_EP_N_POSTJET = parseDoubleString(earlyPlumeSubmenu["Post-jet ice crystal count [#/m] (double)"].as(), "Post-jet ice crystal count [#/m] (double)"); input.ADV_SAVE_PSD_GRID = parseBoolString(advancedNode["Save gridded particle size distribution (T/F)"].as(), "Save gridded particle size distribution (T/F)"); + + if (advancedNode["AEROSOL GRID SUBMENU"]) { + YAML::Node aeroGridSubmenu = advancedNode["AEROSOL GRID SUBMENU"]; + if (aeroGridSubmenu["Ice bin volume ratio [-] (double)"]) { + input.ADV_AERO_ICE_BIN_VRAT = parseDoubleString(aeroGridSubmenu["Ice bin volume ratio [-] (double)"].as(), "Ice bin volume ratio [-] (double)"); + } + if (aeroGridSubmenu["Ice bin r_min [m] (double)"]) { + input.ADV_AERO_ICE_BIN_R_LOW = parseDoubleString(aeroGridSubmenu["Ice bin r_min [m] (double)"].as(), "Ice bin r_min [m] (double)"); + } + if (aeroGridSubmenu["Ice bin r_max [m] (double)"]) { + input.ADV_AERO_ICE_BIN_R_HIG = parseDoubleString(aeroGridSubmenu["Ice bin r_max [m] (double)"].as(), "Ice bin r_max [m] (double)"); + } + } + + if (input.ADV_AERO_ICE_BIN_VRAT <= 1.0) { + throw std::invalid_argument("Ice bin volume ratio [-] in AEROSOL GRID SUBMENU must be strictly greater than 1.0!"); + } + if (input.ADV_AERO_ICE_BIN_R_LOW <= 0.0 || input.ADV_AERO_ICE_BIN_R_HIG <= 0.0) { + throw std::invalid_argument("Ice bin radius limits in AEROSOL GRID SUBMENU must be strictly positive!"); + } + if (input.ADV_AERO_ICE_BIN_R_LOW >= input.ADV_AERO_ICE_BIN_R_HIG) { + throw std::invalid_argument("Ice bin r_min cannot be greater than or equal to r_max in AEROSOL GRID SUBMENU!"); + } } // Wrapper around parseDoubleString to have a nice rejection message for diff --git a/Code.v05-00/tests/test_aerosol.cpp b/Code.v05-00/tests/test_aerosol.cpp index ad21bdc7f..63e92c501 100644 --- a/Code.v05-00/tests/test_aerosol.cpp +++ b/Code.v05-00/tests/test_aerosol.cpp @@ -229,3 +229,99 @@ TEST_CASE ("Aerosol", "[single-file]" ) { } } + +TEST_CASE("Aerosol Flexible Bin Grid Construction", "[aerosol][bins]") { + const double r_low = 5.0e-8; // 0.05 um + const double r_hig = 8.0e-5; // 80.0 um + + SECTION("Default Grid Parameters (VRAT = 1.80)") { + const double vrat = 1.80; + const UInt expected_nbins = static_cast(std::floor(1 + 3.0 * log(r_hig / r_low) / log(vrat))); + REQUIRE(expected_nbins == 38); + + const double rrat = cbrt(vrat); + Vector_1D edges(expected_nbins + 1); + Vector_1D centers(expected_nbins); + edges[0] = r_low; + for (UInt i = 1; i <= expected_nbins; i++) { + edges[i] = edges[i - 1] * rrat; + } + for (UInt i = 0; i < expected_nbins; i++) { + centers[i] = 0.5 * (edges[i] + edges[i + 1]); + } + + REQUIRE(edges.front() == Catch::Approx(5.0e-8)); + REQUIRE(edges.back() >= 8.0e-5); + } + + SECTION("VRAT Spectrum Geometric Invariance (1.2 to 3.0)") { + const std::vector vrats = {1.2, 1.5, 1.8, 2.0, 2.1, 2.4, 2.7, 3.0}; + const std::vector expected_counts = { + static_cast(std::floor(1 + 3.0 * log(r_hig / r_low) / log(1.2))), // ~122 + static_cast(std::floor(1 + 3.0 * log(r_hig / r_low) / log(1.5))), // 55 + static_cast(std::floor(1 + 3.0 * log(r_hig / r_low) / log(1.8))), // 38 + static_cast(std::floor(1 + 3.0 * log(r_hig / r_low) / log(2.0))), // 32 + static_cast(std::floor(1 + 3.0 * log(r_hig / r_low) / log(2.1))), // 30 + static_cast(std::floor(1 + 3.0 * log(r_hig / r_low) / log(2.4))), // 26 + static_cast(std::floor(1 + 3.0 * log(r_hig / r_low) / log(2.7))), // 23 + static_cast(std::floor(1 + 3.0 * log(r_hig / r_low) / log(3.0))) // 21 + }; + + for (size_t k = 0; k < vrats.size(); k++) { + const double vrat = vrats[k]; + const UInt nbins = static_cast(std::floor(1 + 3.0 * log(r_hig / r_low) / log(vrat))); + REQUIRE(nbins == expected_counts[k]); + + const double rrat = cbrt(vrat); + Vector_1D edges(nbins + 1); + Vector_1D centers(nbins); + edges[0] = r_low; + for (UInt i = 1; i <= nbins; i++) { + edges[i] = edges[i - 1] * rrat; + // Ratio test for each adjacent edge + REQUIRE(edges[i] / edges[i - 1] == Catch::Approx(rrat).epsilon(1e-12)); + } + for (UInt i = 0; i < nbins; i++) { + centers[i] = 0.5 * (edges[i] + edges[i + 1]); + } + + REQUIRE(edges.front() == Catch::Approx(r_low)); + REQUIRE(edges.back() >= r_hig); + } + } + + SECTION("Physical Moment Conservation across VRAT Spectrum") { + const double nPart = 1.0e12; // #/cm3 + const double mu = 1.0e-6; // 1.0 um + const double sigma = 1.6; + + const std::vector vrats = {1.2, 1.5, 1.8, 2.0, 2.1, 2.4, 2.7, 3.0}; + const double analytic_m0 = nPart; + const double analytic_m3 = nPart * pow(mu, 3.0) * exp((9.0 / 2.0) * pow(log(sigma), 2.0)); + + for (double vrat : vrats) { + const UInt nbins = static_cast(std::floor(1 + 3.0 * log(r_hig / r_low) / log(vrat))); + const double rrat = cbrt(vrat); + Vector_1D edges(nbins + 1); + Vector_1D centers(nbins); + edges[0] = r_low; + for (UInt i = 1; i <= nbins; i++) { + edges[i] = edges[i - 1] * rrat; + } + for (UInt i = 0; i < nbins; i++) { + centers[i] = 0.5 * (edges[i] + edges[i + 1]); + } + + Aerosol aer(centers, edges, nPart, mu, sigma, "lognormal"); + + // Zeroth moment (total particle number) should be preserved to < 1% + double m0 = aer.Moment(0); + REQUIRE(m0 == Catch::Approx(analytic_m0).epsilon(0.01)); + + // Third moment (volume / mass proportional) should be preserved to < 5% + double m3 = aer.Moment(3); + REQUIRE(m3 == Catch::Approx(analytic_m3).epsilon(0.05)); + } + } +} + From ff9ff8cbb9788bed497c6cca7a7488c4ef5ee5d6 Mon Sep 17 00:00:00 2001 From: "Sebastian D. Eastham" Date: Sun, 23 Aug 2026 16:16:46 -0400 Subject: [PATCH 13/14] Moved timing of turbulent temperature perturbation update --- Code.v05-00/src/Core/LAGRIDPlumeModel.cpp | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp b/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp index 9b4d9639f..6b20b5d2b 100644 --- a/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp +++ b/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp @@ -104,6 +104,13 @@ SimStatus LAGRIDPlumeModel::runFullModel() { std::cout << "\n - Time step: " << timestepVars_.nTime + 1 << " out of " << timestepVars_.timeArray.size(); std::cout << "\n -> Solar time: " << std::fmod( timestepVars_.curr_Time_s/3600.0, 24.0 ) << " [hr]" << std::endl; + /* Update the temperature perturbations at the start of the transport/microphysics step + so that deposition/sublimation subcycling evaluates with the active perturbed temperature field. + */ + if (simVars_.TEMP_PERTURB){ + met_.updateTempPerturb(); + } + // Interleaved Transport and Ice Growth Subcycling over the outer timestep dt if (simVars_.TRANSPORT || simVars_.ICE_GROWTH) { const double dt_step = timestepVars_.dt; @@ -135,14 +142,6 @@ SimStatus LAGRIDPlumeModel::runFullModel() { #endif } - /* With LAGRID remapping every transport timestep, it fundamentally only makes physical sense to update - the temperature perturbations at the same interval as the transport timestep. Turbulence timestep is one - tool used to tune the intensity of the simulated turbulence, but we can also just vary the amplitude. - */ - if (simVars_.TEMP_PERTURB){ - met_.updateTempPerturb(); - } - solarTime_h_ = ( timestepVars_.curr_Time_s + timestepVars_.dt / 2.0 ) / 3600.0; simTime_h_ = ( timestepVars_.curr_Time_s + timestepVars_.dt / 2.0 - timestepVars_.timeArray[0] ) / 3600.0; From ae20628e85d218ee8dc15716ac310023d5653046 Mon Sep 17 00:00:00 2001 From: "Sebastian D. Eastham" Date: Mon, 24 Aug 2026 12:01:42 -0400 Subject: [PATCH 14/14] Fix incorrect time units in output --- Code.v05-00/src/Core/Diag_Mod.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Code.v05-00/src/Core/Diag_Mod.cpp b/Code.v05-00/src/Core/Diag_Mod.cpp index 29230281e..c5de5db9a 100644 --- a/Code.v05-00/src/Core/Diag_Mod.cpp +++ b/Code.v05-00/src/Core/Diag_Mod.cpp @@ -194,7 +194,7 @@ namespace Diag { binRadVar.putAtt("units", "m"); binRadVar.putAtt("long_name", "Ice bin center radius"); binRadVar.putVar(&(iceAer.getBinCenters())[0]); - tVar.putAtt("units", "seconds since simulation start"); + tVar.putAtt("units", "hours since simulation start"); tVar.putAtt("long_name", "time"); tVar.putVar(&(cur_time));