From abc0a5713be79d10cb8ed4199e12a974b7ed4165 Mon Sep 17 00:00:00 2001 From: Marco De Lucia Date: Wed, 16 Mar 2022 14:02:50 +0100 Subject: [PATCH] Added Rcpp-BTCS-{1d,2d}.cpp for clarity --- app/Rcpp-BTCS-1d.cpp | 59 ++++++++++++++++++++++++++++++++++++ app/Rcpp-BTCS-2d.cpp | 60 +++++++++++++++++++++++++++++++++++++ app/Rcpp-interface.R | 71 ++++++++++++++++++++++++++++++++++++++++++-- 3 files changed, 188 insertions(+), 2 deletions(-) create mode 100644 app/Rcpp-BTCS-1d.cpp create mode 100644 app/Rcpp-BTCS-2d.cpp diff --git a/app/Rcpp-BTCS-1d.cpp b/app/Rcpp-BTCS-1d.cpp new file mode 100644 index 0000000..fc54856 --- /dev/null +++ b/app/Rcpp-BTCS-1d.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/Rcpp-BTCS-2d.cpp b/app/Rcpp-BTCS-2d.cpp new file mode 100644 index 0000000..efa3119 --- /dev/null +++ b/app/Rcpp-BTCS-2d.cpp @@ -0,0 +1,60 @@ +#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 & diff2D(int nx, + int ny, + double lenx, + double leny, + std::vector & field, + std::vector & alpha, + double timestep, + int iterations) +{ + // problem dimensionality + int dim = 2; + // total number of grid cells + int n = nx*ny; + + std::vector 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); +} diff --git a/app/Rcpp-interface.R b/app/Rcpp-interface.R index 30601b9..fc1280b 100644 --- a/app/Rcpp-interface.R +++ b/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) + + + + + + + + + + + + + +