[R-sig-dyn-mod] a question on ODE.1D, nspec

[Daniel Reed]

Hi Dabing:

If you'd like to simulate 1D advective flow through a tube, I'd recommend using the advection.1D() function from the ReacTran package. It uses slope limiter methods to reduce numerical diffusion and, therefore, maintains nice sharp features.

Cheers,
Daniel

_______________________________
Daniel C. Reed, PhD.
Postdoctoral Fellow,
Dept. of Earth & Environmental Sciences,
University of Michigan,
Ann Arbor, MI, USA.
email: reeddc at umich.edu
web: www.danielreed.org



On Jun 27, 2013, at 2:42 PM, Dabing Chen <dabing.c at gmail.com> wrote:

> Hi  All:
>     I was trying to simulate a plug flow through a tube and found
> something very confusing to me.
>    I was hoping to see something like the advection only scenario.
>    Instead, the simulation produced a wide spread of mass, and I did not
> see a plug at all.
> 
>     what went wrong?
> 
>    Thanks,
> 
> Dabing
> 
> 
> 
> rm(list=ls())
> library (deSolve)
> 
> 
> # simulate the plug flow of solution dosing with gastric emptying
> N       <- 100   # divide the GIT into 100 segments
> dx      <- 600/N        # grid size of 6 cm, d, total length of GIT 600 cm
> v       <- 600/(4*60)       # velocity of mass movement, cm/min
> x       <- seq(dx/2, 600, by = dx)  # m, distance from river
> 
> 
> 
> 
> state <- c(100,rep(0, N-1))
> 
> plug <- function(t, state, paras) {
> 
>  mass <- state[1:N]
> 
> 
>  conc <- mass / dx  # conc of drug in mg/cm
> 
> 
>  Flux <- v * c(0, mass)   # fluxes due to water transport
> 
> #  Flux[2] <- 0.1*mass[1] #+0.2*mass[1]*(tanh(t%%240-150)-tanh(t%%240-180))
> 
> 
> 
>  ## rate of change = plug flow movement  - absorption
>  dmass        <- -diff(Flux)
> 
> 
>  return(list(c(dmass)))
> }
> 
> 
> times <- seq(0, 1000, by = 1)
> 
> 
> systime <- system.time(out <- ode(y = state, times, plug, parms = NULL))
> 
> 
> color = topo.colors
> drug <- out[,2:(N+1)]
> filled.contour(x = times, y = x, drug, color = color, nlevels = 20,
>               xlab = "time, min", ylab = "GIT distance, cm",main =
> "drug",ylim=c(0,600), xlim=c(0,200))
> 
> write.csv(drug, file = "C:/BI-Chen/Research/R script/learning/plug.csv")
> systime
> 
> 
> 
> On Wed, Jun 26, 2013 at 6:02 PM, Thomas Petzoldt <
> Thomas.Petzoldt at tu-dresden.de> wrote:
> 
>> On 6/26/2013 10:47 PM, Dabing Chen wrote:
>> 
>>> Hi Daniel: I just tried ODE, it worked. Is it true that we do not
>>> want to use ODE.1D if the length is not n*N? I am also wondering what
>>> is the benefit of ODE.1D over ODE? It seems that we can use ODE to
>>> finish the job. Thanks, Dabing
>>> 
>> 
>> Some of the solvers used by ode.1D ... ode.3D make use of the regular
>> structure of such problems specified by nspec resp. dimens, so that an
>> adequate Jacobian matrix can be constructed instead of using the full
>> Jacobian. This makes simulations faster, especially for large systems.
>> 
>> See ?lsodes for more details.
>> 
>> Another advantage are the specific plotting functions for results of
>> ode.1D.
>> 
>> 
>> Thomas
>> 
>> 
> 
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models