diff --git a/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp b/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp index 136ff7335..fd67fe20f 100644 --- a/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp +++ b/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 &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& 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 diff --git a/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp b/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp index d4b6c823d..b6457201c 100644 --- a/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp +++ b/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp @@ -541,6 +541,172 @@ 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, + 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) { + 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; + + 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] + 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] + slope_weight * slope_last; + + 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) { + 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; + + 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] - slope_weight * slope; + } + + // 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]) : (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; + + for (int m = 0; m < N; ++m) { + slice[m] -= cfl_frac * (face_flux[m] - face_flux[m + 1]); + } + } + } + } + + void AdvDiffSystem::semiLagrangianAdvection(double dt, bool parallelAdvection) { + // 1. Horizontal Advection along X (row by row) + #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 if (parallelAdvection) default(shared) + { + std::vector col(ny_); + #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]; + } + } + } + } + + 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..66def456b 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, parallelAdvection); #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, parallelAdvection); #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..df9f5e26e 100644 --- a/Code.v05-00/tests/test_adv_diff_solver.cpp +++ b/Code.v05-00/tests/test_adv_diff_solver.cpp @@ -505,4 +505,120 @@ 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)); + } + } + + 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