From fcef9553c912659370bee5194936f4c59f17f97e Mon Sep 17 00:00:00 2001 From: Marco De Lucia Date: Tue, 15 Mar 2022 18:06:39 +0100 Subject: [PATCH] Added some comparisons in app: RcppFTCS.cpp & Rcpp-interface.{R, cpp} --- app/Rcpp-interface.R | 82 ++++++++++++++++++++++++++++++++++++++++++ app/Rcpp-interface.cpp | 59 ++++++++++++++++++++++++++++++ app/RcppFTCS.cpp | 40 +++++++++++++++++++++ 3 files changed, 181 insertions(+) create mode 100644 app/Rcpp-interface.R create mode 100644 app/Rcpp-interface.cpp create mode 100644 app/RcppFTCS.cpp diff --git a/app/Rcpp-interface.R b/app/Rcpp-interface.R new file mode 100644 index 0000000..ac4999b --- /dev/null +++ b/app/Rcpp-interface.R @@ -0,0 +1,82 @@ +## Time-stamp: "Last modified 2022-03-15 18:04:27 delucia" +library(Rcpp) +library(RcppEigen) +library(ReacTran) +library(deSolve) + +options(width=110) + +setwd("app") + +## This creates the "diff1D" function with our BTCSdiffusion code +sourceCpp("Rcpp-interface.cpp") + +### FTCS explicit (same name) +sourceCpp("RcppFTCS.cpp") + + + +## Grid 101 +## Set initial conditions +N <- 501 +D.coeff <- 1E-3 +C0 <- 1 ## Initial concentration (mg/L) +X0 <- 0 ## Location of initial concentration (m) +## Yini <- c(C0, rep(0,N-1)) + +## Ode1d solution +xgrid <- setup.grid.1D(x.up = 0, x.down = 1, N = N) +x <- xgrid$x.mid +Diffusion <- function (t, Y, parms){ + tran <- tran.1D(C = Y, C.up = 0, C.down = 0, D = parms$D, dx = xgrid) + return(list(tran$dC)) +} + + +## gaussian pulse as initial condition +sigma <- 0.01 +Yini <- exp(-0.5*((x-1/2.0)**2)/sigma**2) + +## plot(x, Yini, type="l") + +parms1 <- list(D=D.coeff) +# 1 timestep, 10 s +times <- seq(from = 0, to = 1, by = 1) + +system.time({ + out1 <- ode.1D(y = Yini, times = times, func = Diffusion, + parms = parms1, dimens = N)[2,-1] +}) + +## Now with BTCS +alpha <- rep(D.coeff, N) + +system.time({ + out2 <- diff1D(n=N, length=1, field=Yini, alpha=alpha, timestep = 1, 0, 0, iterations = 1) +}) + +## plot(out1, out2) +## abline(0,1) + +## matplot(cbind(out1,out2),type="l", col=c("black","red"),lty="solid", lwd=2, +## xlab="grid element", ylab="Concentration", las=1) +## legend("topright", c("ReacTran ode1D", "BTCS 1d"), text.col=c("black","red"), bty = "n") + + + + + +system.time({ + out3 <- RcppFTCS(n=N, length=1, field=Yini, alpha=1E-3, bc_left = 0, bc_right = 0, timestep = 1) +}) + +## Poor man's +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, + 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") + + diff --git a/app/Rcpp-interface.cpp b/app/Rcpp-interface.cpp new file mode 100644 index 0000000..fc54856 --- /dev/null +++ b/app/Rcpp-interface.cpp @@ -0,0 +1,59 @@ +#include "../src/BTCSDiffusion.hpp" // for BTCSDiffusion, BTCSDiffusion::BC_DIRICHLET +#include "../src/BoundaryCondition.hpp" +#include // for copy, max +#include +#include +#include // for std +#include // for vector +#include + +using namespace std; +using namespace Diffusion; +using namespace Rcpp; +//using namespace Eigen; + +// [[Rcpp::depends(RcppEigen)]] +// [[Rcpp::plugins("cpp11")]] +// [[Rcpp::export]] +std::vector & diff1D(int n, + double length, + std::vector & field, + std::vector & alpha, + double timestep, + double bc_left, + double bc_right, + int iterations) { + // dimension of grid + int dim = 1; + + // create input + diffusion coefficients for each grid cell + // std::vector alpha(n, 1 * pow(10, -1)); + // std::vector field(n, 1 * std::pow(10, -6)); + std::vector bc(n, {0,0}); + + // create instance of diffusion module + BTCSDiffusion diffu(dim); + + diffu.setXDimensions(length, n); + + // set the boundary condition for the left ghost cell to dirichlet + bc[0] = {Diffusion::BC_CONSTANT, bc_left}; + bc[n-1] = {Diffusion::BC_CONSTANT, bc_right}; + + // set timestep for simulation to 1 second + diffu.setTimestep(timestep); + + //cout << setprecision(12); + + // loop 100 times + // output is currently generated by the method itself + for (int i = 0; i < iterations; i++) { + diffu.simulate(field.data(), alpha.data(), bc.data()); + } + + // for (auto & cell : field) { + // Rcout << cell << "\n"; + // } + + return(field); +} diff --git a/app/RcppFTCS.cpp b/app/RcppFTCS.cpp new file mode 100644 index 0000000..b1dd30a --- /dev/null +++ b/app/RcppFTCS.cpp @@ -0,0 +1,40 @@ +// Time-stamp: "Last modified 2022-03-15 17:49:54 delucia" +#include +#include // for std +#include // for vector + +using namespace std; +using namespace Rcpp; + +// [[Rcpp::plugins("cpp11")]] +// [[Rcpp::export]] +NumericVector RcppFTCS(int n, + double length, + NumericVector & field, + double alpha, + double bc_left, + double bc_right, + double timestep) +{ + // dimension of grid + NumericVector ext (clone(field)); + double dx = length / ((double) n - 1.); + double dt = 0.5*dx*dx/alpha; + + + double afac = alpha*dt/dx/dx; + int iter = (int) (timestep/dt); + + Rcout << "dt: " << dt << "; inner iterations: " << iter << endl; + + + for (int it = 0; it < iter; it++){ + for (int i = 1; i < ext.size()-1; i++) { + ext[i] = (1. - 2*afac)*ext[i] + afac*(ext[i+1]+ext[i-1]); + } + ext[0] = bc_left; + ext[n-1] = bc_right; + } + return(ext); +} +