Skip to content

Implemented semi-lagrangian advection solver #125

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
sdeastham wants to merge 7 commits into
base: main
Choose a base branch
from
Open

Implemented semi-lagrangian advection solver #125

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
c116a74
Implemented semi-lagrangian advection solver in place of diffusive mu…
sdeastham Aug 18, 2026
68a143c
Correction of error in S-L scheme and addition of docstring
sdeastham Aug 19, 2026
ed9fc4e
Fixed OpenMP handling in advection
sdeastham Aug 19, 2026
45bf28b
Fix variable name error in AdvDiffSystem - zero diff
sdeastham Aug 19, 2026
07f5765
Improved parallelisation in AdvDiff - zero diff
sdeastham Aug 20, 2026
5427df3
Fixed bug in definition of diff_up_0
sdeastham Aug 20, 2026
af408bb
Monotonicity and TVD tests for S-L advection
sdeastham Aug 20, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
72 changes: 72 additions & 0 deletions Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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
Expand Down
166 changes: 166 additions & 0 deletions Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<double>& slice,
double velocity,
double dt,
double ds,
double bc_left,
double bc_right)
{
const int N = static_cast<int>(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<int>(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<double> face_flux(N + 1, 0.0);
Comment thread
lrobion marked this conversation as resolved.
face_flux[0] = bc_left;

auto minmod = [](double a, double b) -> double {

@Calebsakhtar Calebsakhtar Aug 24, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Out of curiosity, why use a Lambda function here? This lambda function is repeated in Line 640 unnecessarily. Either extract the minmod function by turning it into a helper function or define the lambda function once.

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;

@Calebsakhtar Calebsakhtar Aug 24, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These quantities appear to be independent of flow direction and may only need to be calculated once outside the if (velocity > 0.0) / else block:

const double abs_vel  = std::abs(velocity);
const double abs_disp = abs_vel * dt;

int k = static_cast<int>(abs_disp / ds);
double rem_disp = abs_disp - k * ds;
double rem_dt = rem_disp / abs_vel;

double abs_disp = -disp;
int k = static_cast<int>(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<double> face_flux(N + 1, 0.0);
face_flux[N] = bc_right;

auto minmod = [](double a, double b) -> double {

@Calebsakhtar Calebsakhtar Aug 24, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

See Lambda function comment above.

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)
Comment thread
lrobion marked this conversation as resolved.
#pragma omp parallel if (parallelAdvection) default(shared)
{
std::vector<double> 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<double> 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;
Expand Down
48 changes: 6 additions & 42 deletions Code.v05-00/src/FVM_ANDS/FVM_Solver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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();
Expand Down Expand Up @@ -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());
}

Expand Down
Loading