3 port rcpp tug

R/Rphree_Diffusion.R

 ##  Functions for RT simulations using through RcppBTCS's diffusion 
 
 ### Marco De Lucia, delucia@gfz-potsdam.de, 2009-2022
-### Time-stamp: "Last modified 2022-08-09 15:15:37 delucia"
+### Time-stamp: "Last modified 2022-10-04 17:24:58 delucia"
 
 ##' This function is somehow equivalent to
 ##' \code{phreeqc::phrGetSelectedOutput} but it won't call
@@ -60,7 +60,7 @@ DiffRT <- function(setup, init, maxtime,
                    maxsim=1, ebreak=FALSE,
                    scheme="FTCS") {
 
-    require(RcppBTCS)
+    require(RcppTUG)
 
     attach(setup)
     on.exit({
@@ -347,14 +347,16 @@ DiffusionScheme <- function(conc, inflow, dx, dt, alpha, scheme, transported) {
         cnew <- data.matrix(conc)
         if (scheme=="FTCS") {
             for (sp in totrans) {
-                tmp <- as.numeric(c(inflow[sp], cnew[, sp]))
-                cnew[,sp] <- FTCS(tmp, dx, dt, alpha)
+                tmp <- as.numeric(conc[, sp])
+                cnew[,sp] <- FTCS(tmp, dx, dt, alpha, bc_left=inflow[sp])
             }
         } else {
+            ## msg("Diffusing totrans=", paste(totrans, collapse = "; "))
             for (sp in totrans) {
-                ## msg("Diffusing names(inflow[sp])", names(inflow[sp]))
-                tmp <- as.numeric(c(inflow[sp], cnew[, sp]))
-                cnew[,sp] <- BTCS(tmp, dx, dt, alpha)
+                ## msg("Diffusing names(inflow[sp])", names(inflow[sp]), "bc_left=", inflow[sp])
+                ## msg("Colnames(cnew)=", paste(colnames(cnew), collapse = "; "))
+                ## msg("sp=", sp)
+                cnew[,sp] <- as.numeric(BTCS(as.numeric(conc[, sp]), dx, dt, alpha, bc_left=inflow[sp]))
             }
         }
 
@@ -365,11 +367,10 @@ DiffusionScheme <- function(conc, inflow, dx, dt, alpha, scheme, transported) {
         return(cnew)
         
     } else {
-        cj <- c(inflow,conc)
         if (scheme=="FTCS") {
-            cnew <- FTCS(cj, dx, dt, alpha)
+            cnew <- FTCS(conc, dx, dt, alpha, bc_left=inflow)
         } else {
-            cnew <- BTCS(cj, dx, dt, alpha)
+            cnew <- BTCS(conc, dx, dt, alpha, bc_left=inflow)
         }
         if (typeof(cnew)!="double") {
             msg("cnew is not double")
@@ -380,44 +381,48 @@ DiffusionScheme <- function(conc, inflow, dx, dt, alpha, scheme, transported) {
 }
 
 
-##' @title Call \code{RcppBTCS::BTCS1D()}
-##' @param field NumericVector of length n+1 with Dirichlet boundary
-##'     conditions prepended
+##' @title Call \code{RcppTUG::run1D()}
+##' @param field NumericVector of length n
 ##' @param dx the grid spacing in m
 ##' @param dt the required time step in s
 ##' @param alpha diffusion coefficient in m^2/s (constant,
 ##'     homogeneous)
+##' @param bc_left Dirichlet boundary condition at the left inlet
 ##' @return updated concentration vector
 ##' @author MDL
 ##' @export
-BTCS <- function(field, dx, dt, alpha) {
-    require(RcppBTCS)
-    n <- length(field)-1
-    RcppBTCS::BTCS1D(count_x=n,
-                     grid_size_x=dx*n,
-                     field=field[-1],
-                     alpha=alpha,
-                     dt = dt,
-                     dt_divide=1,
-                     bc_left = field[1],
-                     bc_right= -1)
+BTCS <- function(field, dx, dt, alpha, bc_left) {
+    require(RcppTUG)
+    ## msg("length(field):", length(field))
+    ## msg("dim(field):", dim(field))
+    n <- length(field)
+    RcppTUG::run1D(field=field,
+                   alpha=alpha,
+                   dt=dt,
+                   iter=1,
+                   dt_divide = 1,
+                   size_x = dx*n,
+                   left = bc_left,
+                   right = -1,
+                   thomas = TRUE)[[1]]
 }
 
-##' @title Call \code{RcppBTCS::RcppFTCS()}
-##' @param field NumericVector of length n+1 with Dirichlet boundary
-##'     conditions prepended
+##' @title Call \code{RcppTUG::RcppFTCS()}
+##' @param field NumericVector of length n
 ##' @param dx the grid spacing in m
 ##' @param dt the required time step in s
 ##' @param alpha diffusion coefficient in m^2/s (constant,
 ##'     homogeneous)
+##' @param bc_left Dirichlet boundary conditions at the left inlet
 ##' @return updated concentration vector
 ##' @author MDL
 ##' @export
-FTCS <- function(field, dx, dt, alpha) {
-    require(RcppBTCS)
-    RcppBTCS::RcppFTCS(n=length(field)-1, length=1,
-                       field=field[-1], alpha=alpha,
-                       bc_left = field[1], timestep = dt)
+FTCS <- function(field, dx, dt, alpha, bc_left) {
+    require(RcppTUG)
+    field_size <- length(field)*dx
+    RcppTUG::RcppFTCS(n=length(field), length=field_size,
+                      field=field, alpha=alpha,
+                      bc_left = bc_left, timestep = dt)
 }
 
 ##' @title Run a PHREEQC input and retrieve the DUMP