diff --git a/DESCRIPTION b/DESCRIPTION index 93e1224..aba4c29 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -32,3 +32,4 @@ Suggests: covr LinkingTo: Rcpp RoxygenNote: 6.1.1 +Encoding: UTF-8 diff --git a/NAMESPACE b/NAMESPACE index e78a9b5..5675bce 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -3,7 +3,6 @@ S3method(plot,rudnicki.pt) export(.rudnicki_pt) export(.setleft) -export(.surface_J_ringdisk) export(.surface_g) export(.surface_g_pdisk) export(.surface_g_pt) @@ -43,6 +42,7 @@ export(surface_displacement) export(surface_displacement_point) export(surface_displacement_pressuredisk) export(surface_displacement_reservoir) +export(surface_displacement_ringdisk) export(timevarying_fluidmass) export(timevarying_porepressure) export(timevarying_surface_displacement) diff --git a/R/segall92.R b/R/segall92.R index 0bcf252..3a5e495 100644 --- a/R/segall92.R +++ b/R/segall92.R @@ -1,43 +1,95 @@ -segall92 <- function(help=FALSE){ +#' @title Elastostatic subsidence of a half space +#' +segall92 <- function(){ cat("\nThis function is simply a placeholder. See the documentation ( `?segall92` ).\n") - if (help) ?segall92 } -#' @rdname segall92 -#' @export -.surface_J_ringdisk <- function(x, ring.depth, radius, thickness){ - - z_prime <- ring.depth - r <- x - r_prime <- radius - +# Integral in the Green's function that must be integrated +# @rdname segall92 +# @export +.green_segall92 <- function(k., r., a., z_prime., radial. = TRUE){ + # r radial observation point + # a radius of disk + # z_prime depth disk + # k vector along direction of integration (k >= 0) + ka <- k. * a. + kr <- k. * r. + kz <- k. * z_prime. + J0 <- Bessel::BesselJ(kr, 0) + J1 <- Bessel::BesselJ(kr, 1) + J1a <- Bessel::BesselJ(ka, 1) + E <- exp(-1 * kz) + if (radial.){ + # radial displacements at the surface [ur(r,0)] + # Wang (2000) equation 9.29 + J1 * J1a * E + } else { + # vertical displacements at the surface [uz(r,0)] + # Wang (2000) equation 9.30 + J0 * J1 * E + } } - # Integral in the Green's function that must be integrated - .green_integ <- function(k., r., a., z_prime., radial. = TRUE){ - # r radial observation point - # a radius of disk - # z_prime depth disk - # k vector along direction of integration (k >= 0) - ka <- k. * a. - kr <- k. * r. - kz <- k. * z_prime. - J0 <- Bessel::BesselJ(kr, 0) - J1 <- Bessel::BesselJ(kr, 1) - J1a <- Bessel::BesselJ(ka, 1) - E <- exp(-1 * kz) - if (radial.){ - # radial displacements at the surface [ur(r,0)] - # Wang (2000) equation 9.29 - J1 * J1a * E - } else { - # vertical displacements at the surface [uz(r,0)] - # Wang (2000) equation 9.30 - J0 * J1 * E - } - } - -surface_displacement_ringdisk <- function(x, radius, depth, thickness, p0=1, nu=1/4, ShearModulus=15e9, BiotCoef=1, x.lim, verbose=TRUE, tol_pow){ +#' @rdname segall92 +#' @param x +#' @param radius,depth,thickness numeric; the radius, midpoint depth, and thickness of the circular-disk pressure source; the +#' midpoint-depth must be at least \code{thickness/2} +#' @param p0 numeric; the (uniform) pressure change inside the circular-disk source +#' @param nuu numeric; the undrained Poisson's ratio +#' @param ShearModulus numeric; the elastic shear modulus +#' @param BiotCoef numeric; Biot's coefficient (usually written as \eqn{\alpha}) +#' @param verbose logical; should messages be given? +#' @param x.lim numeric; the limit of integration in the radial direction +#' @param tol_pow numeric; the exponent applied to \code{.Machine$double.eps} used as the integration tolerance (see \code{\link{integrate}}) +#' @param nsub integer; the number of subdivisions in the integration (see \code{\link{integrate}}) +#' +#' @references Segall, P. (1992), Induced stresses due to fluid extraction from axisymmetric reservoirs, +#' PAGEOPH, 139: 535. doi: 10.1007/BF00879950 +#' @export +#' +#' @examples +#' +#' # Calculate surface displacements from a pressure-disk source +#' +#' # Radial distances (e.g., observation points) and pressure source strength +#' # (in this case, keep distances in km to stabilize integration +#' # and then rescale pressure to maintain unit consistency) +#' r <- seq(0, 15, length.out=151) +#' Pa <- 80e3 #Pa (N / m^2) +#' pressure <- 2.5 * Pa * (1000^2) +#' +#' # dimensions of inflating (pressure > 0) or deflating (pressure < 0) reservoir +#' disk.radius <- 4 +#' disk.thickness <- 2 +#' disk.middepth <- 3 +#' +#' # calculations... +#' res <- surface_displacement_ringdisk(r, radius = disk.radius, thickness = disk.thickness, depth = disk.middepth, p0 = pressure) +#' head(res) +#' +#' with(res,{ +#' uz_sc <- -uz +#' ur_sc <- ur +#' yl <- range(c(uz_sc, ur_sc, err), na.rm=TRUE) +#' plot(x, uz_sc, col=NA, log='', ylim=yl, xaxs='i', +#' ylab="Displacement, mm", xlab="Distance from source, km", main="Surface Displacements: Pressure-disk Source") +#' abline(h=0, col='grey') +#' +#' lines(x, uz_sc, lwd=1.5) +#' text(10, 7, "Vertical") +#' +#' lines(x, ur_sc, lty=2, lwd=1.5) +#' text(5, 5, "Radial") +#' +#' lines(x, err, lty=3, lwd=1.5) +#' text(4.5, -0.75, "Radial strain\n(mm/km)", pos=2) +#' +#' mtext(sprintf("Segall (1992) model with radius: %s, Thickness: %s, Depth: %s", disk.radius, disk.thickness, disk.middepth), font=3) +#' }) +#' +surface_displacement_ringdisk <- function(x, radius, depth, thickness, + p0=1, nuu=1/4, ShearModulus=15e9, BiotCoef=1, verbose=TRUE, + x.lim, tol_pow, nsub){ if (depth < thickness/2) stop("must be from: to ") @@ -52,18 +104,18 @@ surface_displacement_ringdisk <- function(x, radius, depth, thickness, p0=1, nu= as.numeric(x.lim) } - # Geetsma's uniaxial poroelastic coefficient + # Geetsma's uniaxial poroelastic coefficient (Young's modulus?) # Wang (2000) 3.71 c_m <- BiotCoef / ShearModulus - scaling <- 2 * c_m * (1 - nu) * p0 * thickness * radius + scaling <- 2 * c_m * (1 - nuu) * p0 * thickness * radius # Integration parameters - nsub <- 101L + nsub <- ifelse(missing(nsub), 101L, as.integer(nsub)) tol <- .Machine$double.eps^ifelse(missing(tol_pow), 0.35, tol_pow) .fun <- function(r.., ...){ - int <- integrate(.green_integ, lower=0, upper=x.lim, stop.on.error=FALSE, r. = r.., ..., subdivisions = nsub, rel.tol=tol) + int <- stats::integrate(.green_segall92, lower=0, upper=x.lim, stop.on.error=FALSE, r. = r.., ..., subdivisions = nsub, rel.tol=tol) int[['value']] } diff --git a/README.md b/README.md index 79e3aad..226896d 100644 --- a/README.md +++ b/README.md @@ -4,4 +4,14 @@ Tools for studying crustal deformation from natural and anthropogenic sources [![T-CI Build Status](https://travis-ci.org/abarbour/deform.svg?branch=master)](https://travis-ci.org/abarbour/deform) [![AppVeyor Build Status](https://ci.appveyor.com/api/projects/status/github/abarbour/deform?branch=master&svg=true)](https://ci.appveyor.com/project/abarbour/deform) [![Coverage Status](https://img.shields.io/codecov/c/github/abarbour/deform/master.svg)](https://codecov.io/github/abarbour/deform?branch=master) -## Description \ No newline at end of file +## Description + +## Installation + +This is not yet on CRAN, so: + +```{r} +if ('remotes' %in% rownames(installed.packages())) install.packages('remotes') +library(remotes) +remotes::install_github('abarbour/deform') +``` diff --git a/man/segall92.Rd b/man/segall92.Rd new file mode 100644 index 0000000..e98d7d9 --- /dev/null +++ b/man/segall92.Rd @@ -0,0 +1,83 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/segall92.R +\name{segall92} +\alias{segall92} +\alias{surface_displacement_ringdisk} +\title{Elastostatic subsidence of a half space} +\usage{ +segall92() + +surface_displacement_ringdisk(x, radius, depth, thickness, p0 = 1, + nuu = 1/4, ShearModulus = 1.5e+10, BiotCoef = 1, verbose = TRUE, + x.lim, tol_pow, nsub) +} +\arguments{ +\item{x}{} + +\item{radius, depth, thickness}{numeric; the radius, midpoint depth, and thickness of the circular-disk pressure source; the +midpoint-depth must be at least \code{thickness/2}} + +\item{p0}{numeric; the (uniform) pressure change inside the circular-disk source} + +\item{nuu}{numeric; the undrained Poisson's ratio} + +\item{ShearModulus}{numeric; the elastic shear modulus} + +\item{BiotCoef}{numeric; Biot's coefficient (usually written as \eqn{\alpha})} + +\item{verbose}{logical; should messages be given?} + +\item{x.lim}{numeric; the limit of integration in the radial direction} + +\item{tol_pow}{numeric; the exponent applied to \code{.Machine$double.eps} used as the integration tolerance (see \code{\link{integrate}})} + +\item{nsub}{integer; the number of subdivisions in the integration (see \code{\link{integrate}})} +} +\description{ +Elastostatic subsidence of a half space +} +\examples{ + +# Calculate surface displacements from a pressure-disk source + +# Radial distances (e.g., observation points) and pressure source strength +# (in this case, keep distances in km to stabilize integration +# and then rescale pressure to maintain unit consistency) +r <- seq(0, 15, length.out=151) +Pa <- 80e3 #Pa (N / m^2) +pressure <- 2.5 * Pa * (1000^2) + +# dimensions of inflating (pressure > 0) or deflating (pressure < 0) reservoir +disk.radius <- 4 +disk.thickness <- 2 +disk.middepth <- 3 + +# calculations... +res <- surface_displacement_ringdisk(r, radius = disk.radius, thickness = disk.thickness, depth = disk.middepth, p0 = pressure) +head(res) + +with(res,{ + uz_sc <- -uz + ur_sc <- ur + yl <- range(c(uz_sc, ur_sc, err), na.rm=TRUE) + plot(x, uz_sc, col=NA, log='', ylim=yl, xaxs='i', + ylab="Displacement, mm", xlab="Distance from source, km", main="Surface Displacements: Pressure-disk Source") + abline(h=0, col='grey') + + lines(x, uz_sc, lwd=1.5) + text(10, 7, "Vertical") + + lines(x, ur_sc, lty=2, lwd=1.5) + text(5, 5, "Radial") + + lines(x, err, lty=3, lwd=1.5) + text(4.5, -0.75, "Radial strain\\n(mm/km)", pos=2) + + mtext(sprintf("Segall (1992) model with radius: \%s, Thickness: \%s, Depth: \%s", disk.radius, disk.thickness, disk.middepth), font=3) +}) + +} +\references{ +Segall, P. (1992), Induced stresses due to fluid extraction from axisymmetric reservoirs, +PAGEOPH, 139: 535. doi: 10.1007/BF00879950 +}