Fixed use of run1D for BTCS in DiffRT

DESCRIPTION

@@ -18,5 +18,5 @@ Depends: R (>= 4.0), doParallel, phreeqc, mgcv, graphics, methods, stats, utils,
 Suggests: RcppBTCS
 License: LGPL-2.1
 Encoding: UTF-8
-RoxygenNote: 7.2.0
+RoxygenNote: 7.2.1
 LazyData: true

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 15:01:04 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]))
+                tmp <- as.numeric(c(inflow[sp], conc[, sp]))
                 cnew[,sp] <- FTCS(tmp, dx, dt, alpha)
             }
         } 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,11 @@ DiffusionScheme <- function(conc, inflow, dx, dt, alpha, scheme, transported) {
         return(cnew)
         
     } else {
-        cj <- c(inflow,conc)
         if (scheme=="FTCS") {
+            cj <- c(inflow,conc)
             cnew <- FTCS(cj, dx, dt, alpha)
         } 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,9 +382,8 @@ 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,
@@ -390,20 +391,23 @@ DiffusionScheme <- function(conc, inflow, dx, dt, alpha, scheme, transported) {
 ##' @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()}
+##' @title Call \code{RcppTUG::RcppFTCS()}
 ##' @param field NumericVector of length n+1 with Dirichlet boundary
 ##'     conditions prepended
 ##' @param dx the grid spacing in m
@@ -414,10 +418,10 @@ BTCS <- function(field, dx, dt, alpha) {
 ##' @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)
+    require(RcppTUG)
+    RcppTUG::RcppFTCS(n=length(field)-1, length=1,
+                      field=field[-1], alpha=alpha,
+                      bc_left = field[1], timestep = dt)
 }
 
 ##' @title Run a PHREEQC input and retrieve the DUMP

man/BTCS.Rd

@@ -2,13 +2,12 @@
 % Please edit documentation in R/Rphree_Diffusion.R
 \name{BTCS}
 \alias{BTCS}
-\title{Call \code{RcppBTCS::BTCS1D()}}
+\title{Call \code{RcppTUG::run1D()}}
 \usage{
-BTCS(field, dx, dt, alpha)
+BTCS(field, dx, dt, alpha, bc_left)
 }
 \arguments{
-\item{field}{NumericVector of length n+1 with Dirichlet boundary
-conditions prepended}
+\item{field}{NumericVector of length n}
 
 \item{dx}{the grid spacing in m}
 
@@ -21,7 +20,7 @@ homogeneous)}
 updated concentration vector
 }
 \description{
-Call \code{RcppBTCS::BTCS1D()}
+Call \code{RcppTUG::run1D()}
 }
 \author{
 MDL

man/FTCS.Rd

@@ -2,7 +2,7 @@
 % Please edit documentation in R/Rphree_Diffusion.R
 \name{FTCS}
 \alias{FTCS}
-\title{Call \code{RcppBTCS::RcppFTCS()}}
+\title{Call \code{RcppTUG::RcppFTCS()}}
 \usage{
 FTCS(field, dx, dt, alpha)
 }
@@ -21,7 +21,7 @@ homogeneous)}
 updated concentration vector
 }
 \description{
-Call \code{RcppBTCS::RcppFTCS()}
+Call \code{RcppTUG::RcppFTCS()}
 }
 \author{
 MDL