steady.3D: Steady-state solver for 3-Dimensional ordinary differential equations

Description

Estimates the steady-state condition for a system of ordinary differential equations that result from 3-Dimensional partial differential equation models that have been converted to ODEs by numerical differencing.

It is assumed that exchange occurs only between adjacent layers.

Usage

steady.3D(y, time = 0, func, parms = NULL, nspec = NULL, 
          dimens, names = NULL, method = "stodes", 
          jactype = NULL, cyclicBnd = NULL, times = time, ...)

Arguments

the initial guess of (state) values for the ODE system, a vector.

time, times

time for which steady-state is wanted; the default is times=0. (note- since version 1.7, 'times' has been added as an alias to 'time').

func

either an R-function that computes the values of the derivatives in the ode system (the model defininition) at time time, or a character string giving the name of a compiled function in a dynamically loaded shared library. If func is an R-function, it must be defined as: yprime = func(t, y, parms,...). t is the current time point in the integration, y is the current estimate of the variables in the ODE system. If the initial values y has a names attribute, the names will be available inside func. parms is a vector or list of parameters; ... (optional) are any other arguments passed to the function.

The return value of func should be a list, whose first element is a vector containing the derivatives of y with respect to time, and whose next elements are global values whose steady-state value is also required.

The derivatives should be specified in the same order as the state variables y.

parms

parameters passed to func.

nspec

the number of *species* (components) in the model.

dimens

a 3-valued vector with the dimensionality of the model, i.e. the number of *boxes* in x-, y- and z- direction.

method

the solution method, one of "stodes", or "runsteady". When method = 'runsteady', then solver lsodes, the sparse solver is used by default, unless argument jactype is set to "2D", in which case lsode will be used (likely less efficient). in which case lsodes is used and the structure of the jacobian is determined by the solver.

jactype

the jacobian type - default is a regular 2-D structure where layers only interact with adjacent layers in both directions. If the structure does not comply with this, the jacobian can be set equal to 'sparse'.

cyclicBnd

if not NULL then a number or a 3-valued vector with the dimensions where a cyclic boundary is used - 1: x-dimension, 2: y-dimension; 3: z-dimension;see details.

names

the names of the components; used to label the output, and for plotting.

...

additional arguments passed to function stodes. See note.

Value

A list containing

a vector with the state variable values from the last iteration during estimation of steady-state condition of the system of equations.

...

the "global" values returned.

The output will have the attribute steady, which returns TRUE, if steady-state has been reached and the attribute

precis with the precision attained during each iteration. Another attribute, called dims returns a.o. the length of the work array actually required. This can be specified with input argument lrw. See note and first example.

Details

This is the method of choice to find the steady-state for 3-dimensional models, that are only subjected to transport between adjacent layers.

Based on the dimension of the problem, the method first calculates the sparsity pattern of the Jacobian, under the assumption that transport is only occurring between adjacent layers. Then stodes is called to find the steady-state.

In some cases, a cyclic boundary condition exists. This is when the first boxes in x-, y-, or z-direction interact with the last boxes. In this case, there will be extra non-zero fringes in the Jacobian which need to be taken into account. The occurrence of cyclic boundaries can be toggled on by specifying argument cyclicBnd. For innstance, cyclicBnd = 1 indicates that a cyclic boundary is required only for the x-direction, whereas cyclicBnd = c(1,2) imposes a cyclic boundary for both x- and y-direction. The default is no cyclic boundaries.

As stodes is used, it will probably be necessary to specify the length of the real work array, lrw.

Although a reasonable guess of lrw is made, it may occur that this will be too low. In this case, steady.3D will return with an error message telling the size of the work array actually needed. In the second try then, set lrw equal to this number.

As stodes is used, it will probably be necessary to specify the length of the real work array, lrw.

Although a reasonable guess of lrw is made, it may occur that this will be too low. In this case, steady.2D will return with an error message telling that there was insufficient storage. In the second try then, increase lrw. you may need to experiment to find suitable value. The smaller the better.

The error message that says to increase lrw may look like this:

   
 In stodes(y = y, time = time, func = func, parms = parms, nnz = c(nspec,  :
  insufficient storage in nnfc

See stodes for the additional options.

Examples

Run this code

# NOT RUN {
## =======================================================================
## Diffusion in 3-D; imposed boundary conditions
## =======================================================================
diffusion3D <- function(t, Y, par)   {
   yy    <- array(dim=c(n,n,n),data=Y)  # vector to 3-D array
   dY   <- -r*yy        # consumption
   BND   <- rep(1,n)   # boundary concentration
   for (i in 1:n) {
     y <- yy[i,,]
     #diffusion in X-direction; boundaries=imposed concentration
     Flux <- -Dy * rbind(y[1,]-BND,(y[2:n,]-y[1:(n-1),]),BND-y[n,])/dy
     dY[i,,]   <- dY[i,,] - (Flux[2:(n+1),]-Flux[1:n,])/dy

     #diffusion in Y-direction
     Flux <- -Dz * cbind(y[,1]-BND,(y[,2:n]-y[,1:(n-1)]),BND-y[,n])/dz
     dY[i,,]    <- dY[i,,] - (Flux[,2:(n+1)]-Flux[,1:n])/dz
   }
   for (j in 1:n) {
     y <- yy[,j,]
     #diffusion in X-direction; boundaries=imposed concentration
     Flux <- -Dx * rbind(y[1,]-BND,(y[2:n,]-y[1:(n-1),]),BND-y[n,])/dx
     dY[,j,]   <- dY[,j,] - (Flux[2:(n+1),]-Flux[1:n,])/dx

     #diffusion in Y-direction
     Flux <- -Dz * cbind(y[,1]-BND,(y[,2:n]-y[,1:(n-1)]),BND-y[,n])/dz
     dY[,j,]    <- dY[,j,] - (Flux[,2:(n+1)]-Flux[,1:n])/dz
   }
   for (k in 1:n) {
     y <- yy[,,k]
     #diffusion in X-direction; boundaries=imposed concentration
     Flux <- -Dx * rbind(y[1,]-BND,(y[2:n,]-y[1:(n-1),]),BND-y[n,])/dx
     dY[,,k]   <- dY[,,k] - (Flux[2:(n+1),]-Flux[1:n,])/dx

     #diffusion in Y-direction
     Flux <- -Dy * cbind(y[,1]-BND,(y[,2:n]-y[,1:(n-1)]),BND-y[,n])/dy
     dY[,,k]    <- dY[,,k] - (Flux[,2:(n+1)]-Flux[,1:n])/dy
   }
   return(list(as.vector(dY)))
}

  # parameters
  dy    <- dx <- dz <-1   # grid size
  Dy    <- Dx <- Dz <-1   # diffusion coeff, X- and Y-direction
  r     <- 0.025     # consumption rate

  n  <- 10
  y  <- array(dim=c(n, n, n), data = 10.)

  # stodes is used, so we should specify the size of the work array (lrw)
  # We take a rather large value initially

  print(system.time(
  ST3 <- steady.3D(y, func =diffusion3D, parms = NULL, pos = TRUE,
                   dimens = c(n, n, n), lrw = 100000,
                   atol = 1e-10, rtol = 1e-10, ctol = 1e-10, 
                   verbose = TRUE)
  ))

  # the actual size of lrw is in the attributes()$dims vector.     
  # it is best to set lrw as small as possible 
  attributes(ST3)     

  # image plot
  y <- array(dim=c(n, n, n), data = ST3$y)
  filled.contour(y[ , ,n/2], color.palette = terrain.colors)
    
  # rootSolve's image plot, looping over 3rd dimension
  image(ST3, mfrow = c(4,3))

  # loop over 1st dimension, contours, legends added
  image(ST3, mfrow = c(2, 2), add.contour = TRUE, legend = TRUE,
        dimselect = list(x = c(1, 4, 8, 10)))

# }

Run the code above in your browser using DataCamp Workspace