Added some comparisons in app: RcppFTCS.cpp & Rcpp-interface.{R, cpp}

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")
+
+

app/Rcpp-interface.cpp

@@ -0,0 +1,59 @@
+#include "../src/BTCSDiffusion.hpp" // for BTCSDiffusion, BTCSDiffusion::BC_DIRICHLET
+#include "../src/BoundaryCondition.hpp"
+#include <algorithm> // for copy, max
+#include <cmath>
+#include <iomanip>
+#include <iostream> // for std
+#include <vector>   // for vector
+#include <Rcpp.h>
+
+using namespace std;
+using namespace Diffusion;
+using namespace Rcpp;
+//using namespace Eigen;
+
+// [[Rcpp::depends(RcppEigen)]]
+// [[Rcpp::plugins("cpp11")]]
+// [[Rcpp::export]]
+std::vector<double> & diff1D(int n,
+                             double length,
+                             std::vector<double> & field,
+                             std::vector<double> & 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<double> alpha(n, 1 * pow(10, -1));
+    // std::vector<double> field(n, 1 * std::pow(10, -6));
+    std::vector<boundary_condition> 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);
+}

app/RcppFTCS.cpp

@@ -0,0 +1,40 @@
+//  Time-stamp: "Last modified 2022-03-15 17:49:54 delucia"
+#include <Rcpp.h>
+#include <iostream> // for std
+#include <vector>   // 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);
+}
+