One-Dimensional Advection Equation

Description

Estimates the advection term in a one-dimensional model of a liquid (volume fraction constant and equal to one) or in a porous medium (volume fraction variable and lower than one).

The interfaces between grid cells can have a variable cross-sectional area, e.g. when modelling spherical or cylindrical geometries (see example).

TVD (total variation diminishing) slope limiters ensure monotonic and positive schemes in the presence of strong gradients.

advection.1-D: uses finite differences.

This implies the use of velocity (length per time) and fluxes (mass per unit of area per unit of time).

advection.volume.1D Estimates the volumetric advection term in a one-dimensional model of an aquatic system (river, estuary). This routine is particularly suited for modelling channels (like rivers, estuaries) where the cross-sectional area changes, and hence the velocity changes.

Volumetric transport implies the use of flows (mass per unit of time).

When solved dynamically, the euler method should be used, unless the first-order upstream method is used.

Usage

advection.1D(C, C.up = NULL, C.down = NULL,
  flux.up = NULL, flux.down = NULL, v, VF = 1, A = 1, dx, 
  dt.default = 1, adv.method = c("muscl", "super", "quick", "p3", "up"),
  full.check = FALSE)

advection.volume.1D(C, C.up = C[1], C.down = C[length(C)],
  F.up = NULL, F.down = NULL, flow, V, 
  dt.default = 1, adv.method = c("muscl", "super", "quick", "p3", "up"),
  full.check = FALSE)

Arguments

Details

This implementation is based on the GOTM code

The boundary conditions are either

• zero-gradient.

• fixed concentration.

• fixed flux.

The above order also shows the priority. The default condition is the zero gradient. The fixed concentration condition overrules the zero gradient. The fixed flux overrules the other specifications.

Ensure that the boundary conditions are well defined: for instance, it does not make sense to specify an influx in a boundary cell with the advection velocity pointing outward.

Transport properties:

The advective velocity (v), the volume fraction (VF), and the interface surface (A), can either be specified as one value, a vector, or a 1D property list as generated by setup.prop.1D.

When a vector, this vector must be of length N+1, defined at all grid cell interfaces, including the upper and lower boundary.

The finite difference grid (dx) is specified either as one value, a vector or a 1D grid list, as generated by setup.grid.1D.

Several slope limiters are implemented to obtain monotonic and positive schemes also in the presence of strong gradients, i.e. to reduce the effect of numerical dispersion. The methods are (Pietrzak, 1989, Hundsdorfer and Verwer, 2003):

• "quick": third-order scheme (TVD) with ULTIMATE QUICKEST limiter (quadratic upstream interpolation for convective kinematics with estimated stream terms) (Leonard, 1988)

• "muscl": third-order scheme (TVD) with MUSCL limiter (monotonic upstream centered schemes for conservation laws) (van Leer, 1979).

• "super": third-order scheme (TVD) with Superbee limiter (method=Superbee) (Roe, 1985)

• "p3": third-order upstream-biased polynomial scheme (method=P3)

• "up": first-order upstream ( method=UPSTREAM) - this is the same method as implemented in tran.1D or tran.volume.1D

where "TVD" means a total variation diminishing scheme

Some schemes may produce artificial steepening. Scheme "p3" is not necessarily monotone (may produce negative concentrations!).

If during a certain time step the maximum Courant number is larger than one, a split iteration will be carried out which guarantees that the split step Courant numbers are just smaller than 1. The maximal number of such iterations is set to 100.

These limiters are supposed to work with explicit methods (euler). However, they will also work with implicit methods, although less effectively. Integrate ode.1D only if the model is stiff (see first example).

Value

Note

The advective equation is not checked for mass conservation. Sometimes, this is not an issue, for instance when v represents a sinking velocity of particles or a swimming velocity of organisms.

In others cases however, mass conservation needs to be accounted for.

To ensure mass conservation, the advective velocity must obey certain continuity constraints: in essence the product of the volume fraction (VF), interface surface area (A) and advective velocity (v) should be constant. In sediments, one can use setup.compaction.1D to ensure that the advective velocities for the pore water and solid phase meet these constraints.

In terms of the units of concentrations and fluxes we follow the convention in the geosciences. The concentration C, C.up, C.down as well at the rate of change of the concentration dC are always expressed per unit of phase volume (i.e. per unit volume of solid or liquid).

Total concentrations (e.g. per unit volume of bulk sediment) can be obtained by multiplication with the appropriate volume fraction. In contrast, fluxes are always expressed per unit of total interface area (so here the volume fraction is accounted for).

Author(s)

Karline Soetaert <karline.soetaert@nioz.nl>

References

Pietrzak J (1998) The use of TVD limiters for forward-in-time upstream-biased advection schemes in ocean modeling. Monthly Weather Review 126: 812 .. 830

Hundsdorfer W and Verwer JG (2003) Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations. Springer Series in Computational Mathematics, Springer-Verlag, Berlin, 471 pages

Burchard H, Bolding K, Villarreal MR (1999) GOTM, a general ocean turbulence model. Theory, applications and test cases. Tech Rep EUR 18745 EN. European Commission

Leonard BP (1988) Simple high accuracy resolution program for convective modeling of discontinuities. Int. J. Numer. Meth.Fluids 8: 1291–1318.

Roe PL (1985) Some contributions to the modeling of discontinuous flows. Lect. Notes Appl. Math. 22: 163-193.

van Leer B. (1979) Towards the ultimate conservative difference scheme V. A second order sequel to Godunov's method. J. Comput. Phys. 32: 101-136

See Also

tran.1D, for a discretisation of the general transport equations

Examples

## =============================================================================
## EXAMPLE 1: Testing the various methods - moving a square pulse  
## use of advection.1D
## The tests as in Pietrzak 
## =============================================================================

#--------------------#
# Model formulation  #
#--------------------#
model <- function (t, y, parms,...) {

  adv <- advection.1D(y, v = v, dx = dx, 
     C.up = y[n], C.down = y[1], ...)  # out on one side -> in at other
  return(list(adv$dC))

}

#--------------------#
# Parameters         #
#--------------------#

n     <- 100
dx    <- 100/n
y     <- c(rep(1, 5), rep(2, 20), rep(1, n-25))
v     <- 2 
times <- 0:300   # 3 times out and in

#--------------------#
# model solution     #
#--------------------#

## a plotting function
plotfun <- function (Out, ...) {
  plot(Out[1, -1], type = "l", col = "red", ylab = "y", xlab = "x", ...)
  lines(Out[nrow(Out), 2:(1+n)])
}

# courant number = 2
pm   <- par(mfrow = c(2, 2))

## third order TVD, quickest limiter
out <- ode.1D(y = y, times = times, func = model, parms = 0, hini = 1,
              method = "euler", nspec = 1, adv.method = "quick")

plotfun(out, main = "quickest, euler") 

## third-order ustream-biased polynomial
out2 <- ode.1D(y = y, times = times, func = model, parms = 0, hini = 1,
              method = "euler", nspec = 1, adv.method = "p3")

plotfun(out2, main = "p3, euler")

## third order TVD, superbee limiter
out3 <- ode.1D(y = y, times = times, func = model, parms = 0, hini = 1,
              method = "euler", nspec = 1, adv.method = "super")

plotfun(out3, main = "superbee, euler")

## third order TVD, muscl limiter
out4 <- ode.1D(y = y, times = times, func = model, parms = 0, hini = 1,
              method = "euler", nspec = 1, adv.method = "muscl")

plotfun(out4, main = "muscl, euler")

## =============================================================================
## upstream, different time-steps , i.e. different courant number
## =============================================================================
out <- ode.1D(y = y, times = times, func = model, parms = 0, 
              method = "euler", nspec = 1, adv.method = "up")

plotfun(out, main = "upstream, courant number = 2")

out2 <- ode.1D(y = y, times = times, func = model, parms = 0, hini = 0.5,  
               method = "euler", nspec = 1, adv.method = "up")
plotfun(out2, main = "upstream, courant number = 1")


## Now muscl scheme, velocity against x-axis
y     <- rev(c(rep(0, 5), rep(1, 20), rep(0, n-25)))
v     <- -2.0
out6 <- ode.1D(y = y, times = times, func = model, parms = 0, hini = 1,
               method = "euler", nspec = 1, adv.method = "muscl")

plotfun(out6, main = "muscl, reversed velocity, , courant number = 1")

image(out6, mfrow = NULL)
par(mfrow = pm)


## =============================================================================
## EXAMPLE 2: moving a square pulse in a widening river  
## use of advection.volume.1D
## =============================================================================

#--------------------#
# Model formulation  #
#--------------------#

river.model <- function (t=0, C, pars = NULL, ...)  {

 tran <- advection.volume.1D(C = C, C.up = 0,
                     flow = flow, V = Volume,...)
 return(list(dCdt = tran$dC, F.down = tran$F.down, F.up = tran$F.up))
}

#--------------------#
# Parameters         #
#--------------------#

# Initialising morphology river:

nbox          <- 100                # number of grid cells
lengthRiver   <- 100000             # [m]
BoxLength     <- lengthRiver / nbox # [m]

Distance      <- seq(BoxLength/2, by = BoxLength, len = nbox)   # [m]

# Cross sectional area: sigmoid function of distance [m2]
CrossArea <- 4000 + 72000 * Distance^5 /(Distance^5+50000^5)

# Volume of boxes                          (m3)
Volume    <- CrossArea*BoxLength

# Transport coefficients
flow      <- 1000*24*3600   # m3/d, main river upstream inflow

#--------------------#
# Model solution     #
#--------------------#

pm   <- par(mfrow=c(2,2))

# a square pulse
yini <-  c(rep(10, 10), rep(0, nbox-10))

## third order TVD, muscl limiter
Conc <- ode.1D(y = yini, fun = river.model,  method = "euler", hini = 1,
              parms = NULL, nspec = 1, times = 0:40, adv.method = "muscl")

image(Conc, main = "muscl", mfrow = NULL)
plot(Conc[30, 2:(1+nbox)], type = "l", lwd = 2, xlab = "x", ylab = "C", 
     main = "muscl after 30 days")

## simple upstream differencing
Conc2<- ode.1D(y = yini, fun = river.model, method = "euler", hini = 1,
              parms = NULL, nspec = 1, times = 0:40, adv.method = "up")

image(Conc2, main = "upstream", mfrow = NULL)
plot(Conc2[30, 2:(1+nbox)], type = "l", lwd = 2, xlab = "x", ylab = "C", 
     main = "upstream after 30 days")

par(mfrow = pm)


# Note: the more sophisticated the scheme, the more mass lost/created
# increase tolerances to improve this.

CC <- Conc[ , 2:(1+nbox)]
MASS <- t(CC)*Volume     
colSums(MASS)

## =============================================================================
## EXAMPLE 3: A steady-state solution 
## use of advection.volume.1D
## =============================================================================

Sink <- function (t, y, parms, ...) {
  C1 <- y[1:N]
  C2 <- y[(N+1):(2*N)]
  C3 <- y[(2*N+1):(3*N)]
  # Rate of change= Flux gradient and first-order consumption                  
  
  # upstream can be implemented in two ways:
 
  dC1  <- advection.1D(C1, v = sink, dx = dx, 
           C.up = 100, adv.method = "up", ...)$dC - decay*C1

# same, using tran.1D  
#  dC1  <- tran.1D(C1, v = sink, dx = dx, 
#        C.up = 100, D = 0)$dC -
#           decay*C1

  dC2  <- advection.1D(C2, v = sink, dx = dx, 
        C.up = 100, adv.method = "p3", ...)$dC -
           decay*C2

  dC3  <- advection.1D(C3, v = sink, dx = dx, 
        C.up = 100, adv.method = "muscl", ...)$dC -
           decay*C3

  list(c(dC1, dC2, dC3))
}

N     <- 10
L     <- 1000
dx    <- L/N                          # thickness of boxes
sink  <- 10
decay <- 0.1
out <- steady.1D(runif(3*N), func = Sink, names = c("C1", "C2", "C3"),
        parms = NULL, nspec = 3, bandwidth = 2)
matplot(out$y, 1:N, type = "l", ylim = c(10, 0), lwd = 2, 
  main = "Steady-state")          
legend("bottomright", col = 1:3, lty = 1:3, 
  c("upstream", "p3", "muscl"))

Results

R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(ReacTran)
Loading required package: rootSolve
Loading required package: deSolve

Attaching package: 'deSolve'

The following object is masked from 'package:graphics':

    matplot

Loading required package: shape
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/ReacTran/advection.Rd_%03d_medium.png", width=480, height=480)
> ### Name: advection.1D
> ### Title: One-Dimensional Advection Equation
> ### Aliases: advection.1D advection.volume.1D
> ### Keywords: utilities
> 
> ### ** Examples
> 
> 
> ## =============================================================================
> ## EXAMPLE 1: Testing the various methods - moving a square pulse  
> ## use of advection.1D
> ## The tests as in Pietrzak 
> ## =============================================================================
> 
> #--------------------#
> # Model formulation  #
> #--------------------#
> model <- function (t, y, parms,...) {
+ 
+   adv <- advection.1D(y, v = v, dx = dx, 
+      C.up = y[n], C.down = y[1], ...)  # out on one side -> in at other
+   return(list(adv$dC))
+ 
+ }
> 
> #--------------------#
> # Parameters         #
> #--------------------#
> 
> n     <- 100
> dx    <- 100/n
> y     <- c(rep(1, 5), rep(2, 20), rep(1, n-25))
> v     <- 2 
> times <- 0:300   # 3 times out and in
> 
> #--------------------#
> # model solution     #
> #--------------------#
> 
> ## a plotting function
> plotfun <- function (Out, ...) {
+   plot(Out[1, -1], type = "l", col = "red", ylab = "y", xlab = "x", ...)
+   lines(Out[nrow(Out), 2:(1+n)])
+ }
> 
> # courant number = 2
> pm   <- par(mfrow = c(2, 2))
> 
> ## third order TVD, quickest limiter
> out <- ode.1D(y = y, times = times, func = model, parms = 0, hini = 1,
+               method = "euler", nspec = 1, adv.method = "quick")
> 
> plotfun(out, main = "quickest, euler") 
> 
> ## third-order ustream-biased polynomial
> out2 <- ode.1D(y = y, times = times, func = model, parms = 0, hini = 1,
+               method = "euler", nspec = 1, adv.method = "p3")
> 
> plotfun(out2, main = "p3, euler")
> 
> ## third order TVD, superbee limiter
> out3 <- ode.1D(y = y, times = times, func = model, parms = 0, hini = 1,
+               method = "euler", nspec = 1, adv.method = "super")
> 
> plotfun(out3, main = "superbee, euler")
> 
> ## third order TVD, muscl limiter
> out4 <- ode.1D(y = y, times = times, func = model, parms = 0, hini = 1,
+               method = "euler", nspec = 1, adv.method = "muscl")
> 
> plotfun(out4, main = "muscl, euler")
> 
> ## =============================================================================
> ## upstream, different time-steps , i.e. different courant number
> ## =============================================================================
> out <- ode.1D(y = y, times = times, func = model, parms = 0, 
+               method = "euler", nspec = 1, adv.method = "up")
> 
> plotfun(out, main = "upstream, courant number = 2")
> 
> out2 <- ode.1D(y = y, times = times, func = model, parms = 0, hini = 0.5,  
+                method = "euler", nspec = 1, adv.method = "up")
> plotfun(out2, main = "upstream, courant number = 1")
> 
> 
> ## Now muscl scheme, velocity against x-axis
> y     <- rev(c(rep(0, 5), rep(1, 20), rep(0, n-25)))
> v     <- -2.0
> out6 <- ode.1D(y = y, times = times, func = model, parms = 0, hini = 1,
+                method = "euler", nspec = 1, adv.method = "muscl")
> 
> plotfun(out6, main = "muscl, reversed velocity, , courant number = 1")
> 
> image(out6, mfrow = NULL)
Warning message:
In rep(dots, length.out = n) : 'x' is NULL so the result will be NULL
> par(mfrow = pm)
> 
> 
> ## =============================================================================
> ## EXAMPLE 2: moving a square pulse in a widening river  
> ## use of advection.volume.1D
> ## =============================================================================
> 
> #--------------------#
> # Model formulation  #
> #--------------------#
> 
> river.model <- function (t=0, C, pars = NULL, ...)  {
+ 
+  tran <- advection.volume.1D(C = C, C.up = 0,
+                      flow = flow, V = Volume,...)
+  return(list(dCdt = tran$dC, F.down = tran$F.down, F.up = tran$F.up))
+ }
> 
> #--------------------#
> # Parameters         #
> #--------------------#
> 
> # Initialising morphology river:
> 
> nbox          <- 100                # number of grid cells
> lengthRiver   <- 100000             # [m]
> BoxLength     <- lengthRiver / nbox # [m]
> 
> Distance      <- seq(BoxLength/2, by = BoxLength, len = nbox)   # [m]
> 
> # Cross sectional area: sigmoid function of distance [m2]
> CrossArea <- 4000 + 72000 * Distance^5 /(Distance^5+50000^5)
> 
> # Volume of boxes                          (m3)
> Volume    <- CrossArea*BoxLength
> 
> # Transport coefficients
> flow      <- 1000*24*3600   # m3/d, main river upstream inflow
> 
> #--------------------#
> # Model solution     #
> #--------------------#
> 
> pm   <- par(mfrow=c(2,2))
> 
> # a square pulse
> yini <-  c(rep(10, 10), rep(0, nbox-10))
> 
> ## third order TVD, muscl limiter
> Conc <- ode.1D(y = yini, fun = river.model,  method = "euler", hini = 1,
+               parms = NULL, nspec = 1, times = 0:40, adv.method = "muscl")
> 
> image(Conc, main = "muscl", mfrow = NULL)
Warning message:
In rep(dots, length.out = n) : 'x' is NULL so the result will be NULL
> plot(Conc[30, 2:(1+nbox)], type = "l", lwd = 2, xlab = "x", ylab = "C", 
+      main = "muscl after 30 days")
> 
> ## simple upstream differencing
> Conc2<- ode.1D(y = yini, fun = river.model, method = "euler", hini = 1,
+               parms = NULL, nspec = 1, times = 0:40, adv.method = "up")
> 
> image(Conc2, main = "upstream", mfrow = NULL)
Warning message:
In rep(dots, length.out = n) : 'x' is NULL so the result will be NULL
> plot(Conc2[30, 2:(1+nbox)], type = "l", lwd = 2, xlab = "x", ylab = "C", 
+      main = "upstream after 30 days")
> 
> par(mfrow = pm)
> 
> 
> # Note: the more sophisticated the scheme, the more mass lost/created
> # increase tolerances to improve this.
> 
> CC <- Conc[ , 2:(1+nbox)]
> MASS <- t(CC)*Volume     
> colSums(MASS)
 [1] 400379153 400379153 400379153 400379153 400379153 400379153 400379153
 [8] 400379153 400379153 400379153 400379153 400379153 400379153 400379153
[15] 400379153 400379153 400379153 400379153 400379153 400379153 400379153
[22] 400379153 400379153 400379153 400379153 400379153 400379153 400379153
[29] 400379153 400379153 400379153 400379153 400379153 400379153 400379153
[36] 400379153 400379153 400379151 400378709 400286469 394976947
> 
> ## =============================================================================
> ## EXAMPLE 3: A steady-state solution 
> ## use of advection.volume.1D
> ## =============================================================================
> 
> Sink <- function (t, y, parms, ...) {
+   C1 <- y[1:N]
+   C2 <- y[(N+1):(2*N)]
+   C3 <- y[(2*N+1):(3*N)]
+   # Rate of change= Flux gradient and first-order consumption                  
+   
+   # upstream can be implemented in two ways:
+  
+   dC1  <- advection.1D(C1, v = sink, dx = dx, 
+            C.up = 100, adv.method = "up", ...)$dC - decay*C1
+ 
+ # same, using tran.1D  
+ #  dC1  <- tran.1D(C1, v = sink, dx = dx, 
+ #        C.up = 100, D = 0)$dC -
+ #           decay*C1
+ 
+   dC2  <- advection.1D(C2, v = sink, dx = dx, 
+         C.up = 100, adv.method = "p3", ...)$dC -
+            decay*C2
+ 
+   dC3  <- advection.1D(C3, v = sink, dx = dx, 
+         C.up = 100, adv.method = "muscl", ...)$dC -
+            decay*C3
+ 
+   list(c(dC1, dC2, dC3))
+ }
> 
> N     <- 10
> L     <- 1000
> dx    <- L/N                          # thickness of boxes
> sink  <- 10
> decay <- 0.1
> out <- steady.1D(runif(3*N), func = Sink, names = c("C1", "C2", "C3"),
+         parms = NULL, nspec = 3, bandwidth = 2)
> matplot(out$y, 1:N, type = "l", ylim = c(10, 0), lwd = 2, 
+   main = "Steady-state")          
> legend("bottomright", col = 1:3, lty = 1:3, 
+   c("upstream", "p3", "muscl"))
> 
> 
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>