Last fixes in FTCS

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-10-04 15:01:04 delucia"
+### Time-stamp: "Last modified 2022-10-04 16:00:31 delucia"
 
 ##' This function is somehow equivalent to
 ##' \code{phreeqc::phrGetSelectedOutput} but it won't call
@@ -347,8 +347,8 @@ 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], conc[, 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 = "; "))
@@ -368,8 +368,7 @@ DiffusionScheme <- function(conc, inflow, dx, dt, alpha, scheme, transported) {
         
     } else {
         if (scheme=="FTCS") {
-            cj <- c(inflow,conc)
-            cnew <- FTCS(cj, dx, dt, alpha)
+            cnew <- FTCS(conc, dx, dt, alpha, bc_left=inflow)
         } else {
             cnew <- BTCS(conc, dx, dt, alpha, bc_left=inflow)
         }
@@ -408,20 +407,21 @@ BTCS <- function(field, dx, dt, alpha, bc_left) {
 }
 
 ##' @title Call \code{RcppTUG::RcppFTCS()}
-##' @param field NumericVector of length n+1 with Dirichlet boundary
-##'     conditions prepended
+##' @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 on the left margin
 ##' @return updated concentration vector
 ##' @author MDL
 ##' @export
-FTCS <- function(field, dx, dt, alpha) {
+FTCS <- function(field, dx, dt, alpha, bc_left) {
     require(RcppTUG)
-    RcppTUG::RcppFTCS(n=length(field)-1, length=1,
-                      field=field[-1], alpha=alpha,
-                      bc_left = field[1], timestep = dt)
+    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

man/FTCS.Rd

@@ -4,11 +4,10 @@
 \alias{FTCS}
 \title{Call \code{RcppTUG::RcppFTCS()}}
 \usage{
-FTCS(field, dx, dt, alpha)
+FTCS(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}
 
@@ -16,6 +15,8 @@ conditions prepended}
 
 \item{alpha}{diffusion coefficient in m^2/s (constant,
 homogeneous)}
+
+\item{bc_left}{Dirichlet boundary conditions on the left margin}
 }
 \value{
 updated concentration vector