max/tug
1
0
Fork 0
mirror of https://git.gfz-potsdam.de/naaice/tug.git synced 2025-12-13 17:38:23 +01:00
Code Issues Packages Projects Releases Wiki Activity

Tweaked some stuff in app/Rcpp-interfaceR and FTCS

Browse Source
This commit is contained in:
Marco De Lucia 2022-03-15 18:33:17 +01:00
parent fcef9553c9
commit 3ee1cd999c
2 changed files with 21 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
app/Rcpp-interface.R
View File

@ -1,4 +1,4 @@
## Time-stamp: "Last modified 2022-03-15 18:04:27 delucia" ## Time-stamp: "Last modified 2022-03-15 18:31:55 delucia"
library(Rcpp) library(Rcpp)
library(RcppEigen) library(RcppEigen)
library(ReacTran) library(ReacTran)
@ -18,7 +18,7 @@ sourceCpp("RcppFTCS.cpp")
## Grid 101 ## Grid 101
## Set initial conditions ## Set initial conditions
N <- 501 N <- 1001
D.coeff <- 1E-3 D.coeff <- 1E-3
C0 <- 1 ## Initial concentration (mg/L) C0 <- 1 ## Initial concentration (mg/L)
X0 <- 0 ## Location of initial concentration (m) X0 <- 0 ## Location of initial concentration (m)
@ -34,25 +34,25 @@ Diffusion <- function (t, Y, parms){
## gaussian pulse as initial condition ## gaussian pulse as initial condition
sigma <- 0.01 sigma <- 0.02
Yini <- exp(-0.5*((x-1/2.0)**2)/sigma**2) Yini <- 0.5*exp(-0.5*((x-1/2.0)**2)/sigma**2)
## plot(x, Yini, type="l") ## plot(x, Yini, type="l")
parms1 <- list(D=D.coeff) parms1 <- list(D=D.coeff)
# 1 timestep, 10 s # 1 timestep, 10 s
times <- seq(from = 0, to = 1, by = 1) times <- seq(from = 0, to = 1, by = 0.1)
system.time({ system.time({
out1 <- ode.1D(y = Yini, times = times, func = Diffusion, out1 <- ode.1D(y = Yini, times = times, func = Diffusion,
parms = parms1, dimens = N)[2,-1] parms = parms1, dimens = N)[11,-1]
}) })
## Now with BTCS ## Now with BTCS
alpha <- rep(D.coeff, N) alpha <- rep(D.coeff, N)
system.time({ system.time({
out2 <- diff1D(n=N, length=1, field=Yini, alpha=alpha, timestep = 1, 0, 0, iterations = 1) out2 <- diff1D(n=N, length=1, field=Yini, alpha=alpha, timestep = 0.1, 0, 0, iterations = 10)
}) })
## plot(out1, out2) ## plot(out1, out2)
@ -75,8 +75,18 @@ mm <- colMeans(rbind(out2,out3))
matplot(cbind(Yini, out1, out2, out3, mm),type="l", col=c("black","grey","red","blue","green4"), lty=c("dashed", rep("solid",4)), lwd=2, matplot(cbind(Yini,out1, out2, out3, mm),type="l", col=c("grey","black","red","blue","green4"), lty="solid", lwd=2,
xlab="grid element", ylab="Concentration", las=1) xlab="grid element", ylab="Concentration", las=1)
legend("topright", c("Init","ReacTran ode1D", "BTCS 1d", "FTCS", "poor man's CN"), text.col=c("black","grey","red","blue","green4"), bty = "n") legend("topright", c("init","ReacTran ode1D", "BTCS 1d", "FTCS", "poor man's CN"), text.col=c("grey","black","red","blue","green4"), bty = "n")
sum(Yini)
sum(out1)
sum(out2)
sum(out3)
sum(mm)
## Yini <- 0.2*sin(pi/0.1*x)+0.2
## plot(Yini)
## plot(out3)

4
app/RcppFTCS.cpp
View File

@ -1,4 +1,4 @@
// Time-stamp: "Last modified 2022-03-15 17:49:54 delucia" // Time-stamp: "Last modified 2022-03-15 18:15:39 delucia"
#include <Rcpp.h> #include <Rcpp.h>
#include <iostream> // for std #include <iostream> // for std
#include <vector> // for vector #include <vector> // for vector
@ -19,7 +19,7 @@ NumericVector RcppFTCS(int n,
// dimension of grid // dimension of grid
NumericVector ext (clone(field)); NumericVector ext (clone(field));
double dx = length / ((double) n - 1.); double dx = length / ((double) n - 1.);
double dt = 0.5*dx*dx/alpha; double dt = 0.25*dx*dx/alpha;
double afac = alpha*dt/dx/dx; double afac = alpha*dt/dx/dx;