Process subcycling and time-averaged diffusion

Code.v05-00/defaults/input.yaml

@@ -84,7 +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
+  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
@@ -106,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
@@ -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

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;
@@ -138,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;
@@ -149,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;
 
 };
 

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<std::vector<int>>& mask, const VectorUtils::MaskInfo& maskInfo);
         std::pair<LAGRID::twoDGridVariable,LAGRID::twoDGridVariable> remapVariable(const VectorUtils::MaskInfo& maskInfo, const BufferInfo& buffers, const Vector_2D& phi, const std::vector<std::vector<int>>& mask);
         double totalAirMass();

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;
     }
 };

Code.v05-00/include/Defaults/Input.hpp

@@ -85,7 +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
+  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 +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
@@ -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
@@ -192,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

Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp

@@ -11,6 +11,77 @@ namespace FVM_ANDS{
     // Separate the SOR solver for testing without having to build an AdvDiffSystem object
     void sor_solve(const Eigen::SparseMatrix<double, Eigen::RowMajor> &A, const Eigen::VectorXd &rhs, Eigen::VectorXd &phi, double omega = 1.0, double threshold = 1e-3, int n_iters = 3);
 
+    /**
+     * @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<double>& 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 +106,7 @@ namespace FVM_ANDS{
             const Eigen::VectorXd& calcRHS();
             void applyBoundaryCondition();
             void updateBoundaryCondition(const BoundaryConditions& bc);
+            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

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 );
 }
 
 

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));