Added Rcpp-BTCS-{1d,2d}.cpp for clarity

app/Rcpp-BTCS-1d.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/Rcpp-BTCS-2d.cpp

@@ -0,0 +1,60 @@
+#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> & diff2D(int nx,
+                             int ny,
+                             double lenx,
+                             double leny,
+                             std::vector<double> & field,
+                             std::vector<double> & alpha,
+                             double timestep,
+                             int iterations)
+{
+    // problem dimensionality
+    int dim = 2;
+    // total number of grid cells
+    int n = nx*ny;
+    
+    std::vector<boundary_condition> bc(n, {0,0});
+    
+    // create instance of diffusion module
+    BTCSDiffusion diffu(dim);
+    
+    diffu.setXDimensions(lenx, nx);
+    diffu.setXDimensions(leny, ny);
+    
+    // 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/Rcpp-interface.R

@@ -1,4 +1,4 @@
-## Time-stamp: "Last modified 2022-03-15 18:31:55 delucia"
+## Time-stamp: "Last modified 2022-03-16 14:01:11 delucia"
 library(Rcpp)
 library(RcppEigen)
 library(ReacTran)
@@ -9,7 +9,7 @@ options(width=110)
 setwd("app")
 
 ## This creates the "diff1D" function with our BTCSdiffusion code
-sourceCpp("Rcpp-interface.cpp")
+sourceCpp("Rcpp-BTCS-1d.cpp")
 
 ### FTCS explicit (same name)
 sourceCpp("RcppFTCS.cpp")
@@ -90,3 +90,70 @@ sum(mm)
 ## plot(Yini)
 
 ## plot(out3)
+
+Fun <- function(dx) {
+    tmp <- diff1D(n=N, length=1, field=Yini, alpha=alpha, timestep = dx, 0, 0, iterations = floor(1/dx))
+    sqrt(sum({out1-tmp}^2))
+}
+
+reso <- optimise(f=Fun, interval=c(1E-5, 1E-1), maximum = FALSE)
+
+
+dx <- 0.0006038284
+floor(1/dx)
+
+1/dx
+
+system.time({
+    out2o <- diff1D(n=N, length=1, field=Yini, alpha=alpha, timestep = dx, 0, 0, iterations = 1656)
+})
+
+matplot(cbind(out1, out2o),type="l", col=c("black","red"), lty="solid", lwd=2,
+        xlab="grid element", ylab="Concentration", las=1)
+legend("topright", c("ReacTran ode1D", "BTCS 1d dx=0.0006"), text.col=c("black","red"), bty = "n")
+
+
+dx <- 0.05
+
+system.time({
+    out2o <- diff1D(n=N, length=1, field=Yini, alpha=alpha, timestep = dx, 0, 0, iterations = 1/dx)
+})
+
+matplot(cbind(out1, out2o),type="l", col=c("black","red"), lty="solid", lwd=2,
+        xlab="grid element", ylab="Concentration", las=1)
+legend("topright", c("ReacTran ode1D", "BTCS 1d dx=0.0006"), text.col=c("black","red"), bty = "n")
+
+Matplot
+
+
+
+
+## This creates the "diff1D" function with our BTCSdiffusion code
+sourceCpp("Rcpp-BTCS-2d.cpp")
+
+n <- 256
+a2d <- rep(1E-3, n^2)
+
+init2d <- readRDS("gs1.rds")
+
+ll <- {init2d - min(init2d)}/diff(range(init2d))
+
+system.time({
+    res1 <- diff2D(nx=N, ny=N, lenx=1, leny=1, field=ll, alpha=a2d, timestep = 0.1, iterations = 10)
+})
+
+hist(ll,32)
+
+
+
+
+
+
+
+
+
+
+
+
+
+