ReacTran (version 1.4.3.1)

fiadeiro: Advective Finite Difference Weights

Description

Weighing coefficients used in the finite difference scheme for advection calculated according to Fiadeiro and Veronis (1977).

This particular AFDW (advective finite difference weights) scheme switches from backward differencing (in advection dominated conditions; large Peclet numbers) to central differencing (under diffusion dominated conditions; small Peclet numbers).

This way it forms a compromise between stability, accuracy and reduced numerical dispersion.

Usage

fiadeiro(v, D, dx.aux = NULL, grid = list(dx.aux = dx.aux))

Arguments

v

advective velocity; either one value or a vector of length N+1, with N the number of grid cells [L/T]

D

diffusion coefficient; either one value or a vector of length N+1 [L2/T]

dx.aux

auxiliary vector containing the distances between the locations where the concentration is defined (i.e. the grid cell centers and the two outer interfaces); either one value or a vector of length N+1 [L]

grid

discretization grid as calculated by setup.grid.1D

Value

the Advective Finite Difference Weighing (AFDW) coefficients as used in the transport routines tran.1D and tran.volume.1D; either one value or a vector of length N+1 [-]

Details

The Fiadeiro and Veronis (1977) scheme adapts the differencing method to the local situation (checks for advection or diffusion dominance).

Finite difference schemes are based on following rationale:

• When using forward differences (AFDW = 0), the scheme is first order accurate, creates a low level of (artificial) numerical dispersion, but is highly unstable (state variables may become negative).

• When using backward differences (AFDW = 1), the scheme is first order accurate, is universally stable (state variables always remain positive), but the scheme creates high levels of numerical dispersion.

• When using central differences (AFDW = 0.5), the scheme is second order accurate, is not universally stable, and has a moderate level of numerical dispersion, but state variables may become negative.

Because of the instability issue, forward schemes should be avoided. Because of the higher accuracy, the central scheme is preferred over the backward scheme.

The central scheme is stable when sufficient physical dispersion is present, it may become unstable when advection is the only transport process.

The Fiadeiro and Veronis (1977) scheme takes this into account: it uses central differencing when possible (when physical dispersion is high enough), and switches to backward differing when needed (when advection dominates). The switching is determined by the Peclet number

Pe = abs(v)*dx.aux/D

• the higher the diffusion D (Pe > 1), the closer the AFDW coefficients are to 0.5 (central differencing)

• the higher the advection v (Pe < 1), the closer the AFDW coefficients are to 1 (backward differencing)

References

• Fiadeiro ME and Veronis G (1977) Weighted-mean schemes for finite-difference approximation to advection-diffusion equation. Tellus 29, 512-522.

• Boudreau (1997) Diagnetic models and their implementation. Chapter 8: Numerical Methods. Springer.

Examples

# NOT RUN {
#============================================
# Model formulation (differential equations)
#============================================

# This is a test model to evaluate the different finite difference schemes 
# and evaluate their effect on munerical diffusion. The model describes the
# decay of organic carbon (OC) as it settles through the ocean water column.

model <- function (time, OC, pars, AFDW = 1) {
  dOC <- tran.1D(OC, flux.up = F_OC, D = D.eddy, 
                 v = v.sink, AFDW = AFDW, dx = dx)$dC - k*OC
  return(list(dOC))
}
#============================================
# Parameter set
#============================================

L <- 1000         # water depth model domain [m]
x.att <- 200      # attenuation depth of the sinking velocity [m]
v.sink.0 <- 10    # sinking velocity at the surface [m d-1]
D.eddy <- 10      # eddy diffusion coefficient [m2 d-1]
F_OC <- 10        # particle flux [mol m-2 d-1]
k <- 0.1          # decay coefficient [d-1]

## =============================================================================
## Model solution for a coarse grid (10 grid cells)
## =============================================================================

# Setting up the grid
N <- 10                               # number of grid layers 
dx <- L/N                             # thickness of boxes [m]
dx.aux <- rep(dx, N+1)                # auxilliary grid vector
x.int <- seq(from = 0, to = L, by = dx)    # water depth at box interfaces [m]
x.mid <- seq(from = dx/2, to = L, by = dx) # water depth at box centres [m]

# Exponentially declining sink velocity
v.sink <- v.sink.0 * exp(-x.int/x.att) # sink velocity [m d-1]
Pe <- v.sink * dx/D.eddy               # Peclet number

# Calculate the weighing coefficients
AFDW <- fiadeiro(v = v.sink, D = D.eddy, dx.aux = dx.aux)

par(mfrow = c(2, 1), cex.main = 1.2, cex.lab = 1.2)

# Plot the Peclet number over the grid 

plot(Pe, x.int, log = "x", pch = 19, ylim = c(L,0), xlim = c(0.1, 1000), 
     xlab = "", ylab = "depth [m]", 
     main = "Peclet number", axes = FALSE)
abline(h = 0)
axis(pos = NA, side = 2)
axis(pos = 0, side = 3)

# Plot the AFDW coefficients over the grid 

plot(AFDW, x.int, pch = 19, ylim = c(L, 0), xlim = c(0.5, 1), 
     xlab = "", ylab = "depth [m]", main = "AFDW coefficient", axes = FALSE)
abline(h = 0)
axis(pos = NA, side = 2)
axis(pos = 0, side = 3)

# Three steady-state solutions for a coarse grid based on:
# (1) backward differences (BD)
# (2) central differences (CD)
# (3) Fiadeiro & Veronis scheme (FV)

BD <- steady.1D(y = runif(N), func = model, AFDW = 1.0, nspec = 1)
CD <- steady.1D(y = runif(N), func = model, AFDW = 0.5, nspec = 1)
FV <- steady.1D(y = runif(N), func = model, AFDW = AFDW, nspec = 1)

# Plotting output - use rootSolve's plot method
plot(BD, CD, FV, grid = x.mid, xyswap = TRUE, mfrow = c(1,2), 
     xlab = "", ylab = "depth [m]", main = "conc (Low resolution grid)")

legend("bottomright", col = 1:3, lty = 1:3,
       legend = c("backward diff", "centred diff", "Fiadeiro&Veronis"))


## =============================================================================
## Model solution for a fine grid (1000 grid cells)
## =============================================================================

# Setting up the grid
N <- 1000                            # number of grid layers 
dx <- L/N                            # thickness of boxes[m]
dx.aux <- rep(dx, N+1)              # auxilliary grid vector
x.int <- seq(from = 0, to = L, by = dx)      # water depth at box interfaces [m]
x.mid <- seq(from = dx/2, to = L, by = dx)   # water depth at box centres [m]

# Exponetially declining sink velocity
v.sink <- v.sink.0 * exp(-x.int/x.att) # sink velocity [m d-1]
Pe <- v.sink * dx/D.eddy               # Peclet number

# Calculate the weighing coefficients
AFDW <- fiadeiro(v = v.sink, D = D.eddy, dx.aux = dx.aux)

# Three steady-state solutions for a coarse grid based on:
# (1) backward differences (BD)
# (2) centered differences (CD)
# (3) Fiadeiro & Veronis scheme (FV)

BD <- steady.1D(y = runif(N), func = model, AFDW = 1.0, nspec = 1)
CD <- steady.1D(y = runif(N), func = model, AFDW = 0.5, nspec = 1)
FV <- steady.1D(y = runif(N), func = model, AFDW = AFDW, nspec = 1)

# Plotting output
plot(BD, CD, FV, grid = x.mid, xyswap = TRUE, mfrow = NULL, 
     xlab = "", ylab = "depth [m]", main = "conc (High resolution grid)")

legend("bottomright", col = 1:3, lty = 1:3,
       legend = c("backward diff", "centred diff", "Fiadeiro&Veronis"))

# Results and conclusions:
# - For the fine grid, all three solutions are identical
# - For the coarse grid, the BD and FV solutions show numerical dispersion
#   while the CD provides more accurate results
# }