From c6dbf07664ea97fbaf3475b299aa502cf60bb6c3 Mon Sep 17 00:00:00 2001 From: Andrew Barbour Date: Fri, 7 Nov 2014 13:25:26 -0800 Subject: [PATCH] Added time-varying pore pressure field --- DESCRIPTION | 3 +- NAMESPACE | 1 + NOBUILD/Sandbox/expandgrid.R | 11 +++++++ NOBUILD/Sandbox/segall.R | 45 +++++++++++++++++++++++--- R/segall85.R | 62 ++++++++++++++++++++++++++++++++++++ man/segall85.Rd | 5 +++ 6 files changed, 121 insertions(+), 6 deletions(-) create mode 100644 NOBUILD/Sandbox/expandgrid.R diff --git a/DESCRIPTION b/DESCRIPTION index b1d3c78..a86a593 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -17,7 +17,8 @@ Depends: sp, maptools, plyr, - reshape2 + reshape2, + abind Suggests: maps LinkingTo: Rcpp diff --git a/NAMESPACE b/NAMESPACE index ad34876..a3d7461 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -18,6 +18,7 @@ export(rigid_hydraulic_diffusivity) export(segall85) export(surface_displacement) export(timevarying_fluidmass) +export(timevarying_porepressure) export(timevarying_surface_displacement) exportPattern("^[[:alpha:]]+") importFrom(Rcpp,evalCpp) diff --git a/NOBUILD/Sandbox/expandgrid.R b/NOBUILD/Sandbox/expandgrid.R new file mode 100644 index 0000000..7445e5e --- /dev/null +++ b/NOBUILD/Sandbox/expandgrid.R @@ -0,0 +1,11 @@ +x <- 1:3 +z <- c(10,20,30,40) +grd <- as.matrix(expand.grid(VarX=x,VarZ=z)) + +testfun <- function(XY){ + sum(XY) +} + +Res <- apply(grd, 1, testfun) + +matrix(Res, ncol=length(z)) diff --git a/NOBUILD/Sandbox/segall.R b/NOBUILD/Sandbox/segall.R index 0c0acfb..838f122 100644 --- a/NOBUILD/Sandbox/segall.R +++ b/NOBUILD/Sandbox/segall.R @@ -1,4 +1,5 @@ library(deform) +library(kook) library(TeachingDemos) .x. <- sort(unique(c(-7:-3, seq(-3.90,3.90,by=0.15), 3:7))) @@ -33,6 +34,7 @@ F2 <- function(){ # Figs 7,8 mxx <- 50 .x.km. <- sort(unique(c((-1*mxx):-3, seq(-2.90,2.90,by=0.1), 3:mxx))) +.z.km. <- sort(unique(c(seq(0,3,by=0.25),seq(3,12,by=0.75)))) yr <- 365*86400 .time. <- seq(2,10,by=2)*10*yr .Vdot. <- 2e6/yr # volume rate m^3/yr to m^3/s @@ -41,9 +43,12 @@ yr <- 365*86400 .B. <- 0.6 # Skemptons coeff .c. <- 0.1 # hydraulic diffusivity m^2/s .Sources.x. <- 1e3*c(0) +.TwoSources.x. <- 1e3*c(0,20) # for mass computation .t. <- 100 # thickness .phi. <- 0.2 #porosity +# for pressure computation +.mu. <- 5.6 #GPa -- shear modulus zz2 <- timevarying_surface_displacement(.x.km.*1e3, .time., .Vdot., .B., .L., .D., .c., Pt.Sources.x=.Sources.x.) zz2t <- apply(zz2, 2, function(.z.) matrix(Tilt(.x.km.*1e3, z=.z.)$ztilt)) @@ -53,6 +58,15 @@ F3 <- function(){ matplot(.x.km., zz2*1e3, type="l", col="black", main="Subsidence, mm, Segall 1985, Fig 8B", sub=Sys.time()) } +F3t <- function(){ + #matplot(.time./yr, t(zz2t), type="l") + matplot(.x.km., zz2t, type="l", col="black") +} + +#try(F3()) +#try(F3t()) + + zz3 <- timevarying_fluidmass(.x.km.*1e3, .time., .Vdot., .L., .t., .c., phi.=.phi.) F4 <- function(){ @@ -62,12 +76,33 @@ F4 <- function(){ #try(F4()) +zzp <- timevarying_porepressure(.x.km.*1e3, .z.km.*1e3, .time., .Vdot., .B., .L., .D., .c., .t., .mu., Pt.Sources.x=.TwoSources.x.) -F3t <- function(){ - #matplot(.time./yr, t(zz2t), type="l") - matplot(.x.km., zz2t, type="l", col="black") +F5 <- function(do.log=FALSE){ + #matplot(.time./yr, t(zz3)*1e2, type="l") + X<- zzp[,,length(.time.)] + if (do.log) X <- log10(abs(X)) + matplot(x=.x.km., X, col=NA) + aaply(zzp, 3, .fun = function(X) { + if (do.log) X <- log10(abs(X)) + matplot(x=.x.km., X, type="l", add=TRUE) + return("x") + }) + invisible() } -try(F3()) -#try(F3t()) +try(F5()) + + +F5c <- function(){ + #matplot(.time./yr, t(zz3)*1e2, type="l") + #contour(x=.x.km., y=.z.km., zzp[,,length(.time.)], col=NA) + aaply(zzp, 3, .fun = function(X) { + image(x=.x.km., y=.z.km., X, ylim=c(12,0), col = brewerRamp()) + contour(x=.x.km., y=.z.km., X, ylim=c(12,0), add = TRUE) + return("x") + }) + invisible() +} +#try(F5c()) diff --git a/R/segall85.R b/R/segall85.R index 7d66288..c964f3d 100644 --- a/R/segall85.R +++ b/R/segall85.R @@ -88,6 +88,68 @@ timevarying_surface_displacement <- function(x, Time, Vdot., B., L., D., HD., nu apply(X=matrix(Time), MARGIN=1, FUN=.FUN, Xi.sources=Sources) } +#' @rdname segall85 +#' @export +timevarying_porepressure <- function(x, z, Time, Vdot., B., L., D., HD., t., mu.gpa., nu.=1/4, nuu.=1/3, Pt.Sources.x=0, x.lim=round(max(abs(x),na.rm=TRUE))/1){ + # Time varying p.p. associated with fluid extraction + # segall85 eq 28 + # + sc <- 1e9 + mu. <- sc * mu.gpa. + # + Sources <- matrix(Pt.Sources.x, ncol=1) + # + .FUN <- function(ti, zi=0, Xi.sources){ + message(paste("Time:", ti/365/86400, "years")) + # time varying component (scaling) + tvar1 <- -2 * mu. * (1 + nuu.)^2 * B.^2 * Vdot. * sqrt(ti / HD.) / (9 * L. * t.) + # + # Function that returns the integral value at + # a position Xi away from the source at dXi + .fun <- function(Xi., dXi.=0, X.=x, Z.=0, ti.=ti){ + # from the spatio-temporal component + ti.d <- Xi.^2 / (4 * HD. * ti.) + C1 <- ierfc2(sqrt(ti.d)) + # + xmxi <- (X. - Xi. - dXi.)^2 + zpd <- (Z. + D.) + C2A <- zpd / (zpd^2 + xmxi^2) + zpdt <- zpd + t. + C2B <- -1 * zpdt / (zpdt^2 + xmxi^2) + sc * (C2A + C2B) * C1 + } + # + # function used to integrate over Xi at different sources + .xi.integ2D <- function(dxi, ti.=ti){ + message(paste("\t2D integration of source at:", dxi/1e3, "km")) + grd <- as.matrix(expand.grid(VarX=x, VarZ=z)) + ires <- apply(grd, 1, FUN=function(xxzz){ + #message(paste(xxzz, collapse="//")) + zi <- try(integrate(.fun, -1*x.lim, x.lim, X.=xxzz[1], Z.=xxzz[2], dXi.=dxi, subdivisions = 150L, rel.tol=.Machine$double.eps^0.5)) + if (inherits(zi,"try-error")){ + #message(paste("\t==== ==== ====> shrink integration limits to below",x.lim)) + message(paste(xxzz, collapse=" x // z ")) + NA + } else { + #message(paste("integration",zi$message)) + ti.d <- xxzz[1]^2 / (4 * HD. * ti.) + C1 <- (1 - 2*nu.) * ierfc2(sqrt(ti.d)) / ((nuu. - nu.) * (1 - 2*nuu.)) + C2 <- -2/(pi * (1 - nuu.)) + C1 + C2 * zi$value + } + }) + matrix(ires, ncol=length(z), byrow = FALSE) + } + # superpostion of all sources + xvar.t2s <- apply(X = abind(lapply(X = Sources, FUN = .xi.integ2D), along=3), MARGIN = c(1,2), FUN = sum, na.rm = TRUE) + # combine and return + res <- tvar1 * xvar.t2s + return(-1*res/sc) + } + # apply FUN through time + abind(lapply(X = Time, FUN = .FUN, Xi.sources=Sources), along=3) +} + #' Simple numerical deformation estimates #' @name Simple-deformation #' @param x numeric; the spatial vector in real units (not an index!) diff --git a/man/segall85.Rd b/man/segall85.Rd index 7229318..018e64c 100644 --- a/man/segall85.Rd +++ b/man/segall85.Rd @@ -4,6 +4,7 @@ \alias{segall85} \alias{surface_displacement} \alias{timevarying_fluidmass} +\alias{timevarying_porepressure} \alias{timevarying_surface_displacement} \title{Surface deformation associated with fluid withdrawl} \usage{ @@ -17,6 +18,10 @@ timevarying_fluidmass(x, Time, Vdot., L., t., HD., phi.) timevarying_surface_displacement(x, Time, Vdot., B., L., D., HD., nuu. = 1/3, Pt.Sources.x = 0, x.lim = 2 * round(max(abs(x), na.rm = TRUE))) + +timevarying_porepressure(x, z, Time, Vdot., B., L., D., HD., t., mu.gpa., + nu. = 1/4, nuu. = 1/3, Pt.Sources.x = 0, x.lim = round(max(abs(x), + na.rm = TRUE))/1) } \description{ Surface deformation associated with fluid withdrawl