Added rice cleary solution

.gitignore

@@ -1,8 +1,9 @@
+Rplots.pdf
 *.Rproj
 .Rproj.user
 .Rhistory
 .RData
 src/*.o
 src/*.so
 src/*.dll
-.DS_Store
+.DS_Store

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: deform
 Type: Package
 Title: Crustal Deformation Modeling Tools in R
-Version: 0.0.989
-Date: 2020-06-19
+Version: 0.0.990
+Date: 2021-03-30
 Author: Andrew J. Barbour; some of these tools were
     adapted from Matlab(tm) codes by Brendan Crowell (UW),
     and Francois Beauducel (IPGP).
@@ -21,8 +21,10 @@ Imports:
     RcppFaddeeva,
     rgdal (>= 0.8-16),
     raster,
+    tibble,
     maptools,
     abind,
+    dplyr,
     plyr,
     RCurl,
     Bessel

NAMESPACE

@@ -8,6 +8,7 @@ export(.surface_g_pdisk)
 export(.surface_g_pt)
 export(.surface_g_xi)
 export(D_project)
+export(Mw2M0)
 export(Tilt)
 export(Uniaxial_extension)
 export(agency.root)
@@ -33,6 +34,7 @@ export(ierfc_re)
 export(lame_first)
 export(mogi.pressure)
 export(mogi.volume)
+export(poroelastic_stress_coefficient)
 export(r96_drainage_pressure)
 export(rigid_hydraulic_diffusivity)
 export(rudnicki86)
@@ -48,6 +50,7 @@ export(timevarying_fluidmass)
 export(timevarying_porepressure)
 export(timevarying_surface_displacement)
 export(vasco98)
+import(dplyr)
 importFrom(Bessel,BesselJ)
 importFrom(RCurl,getURL)
 importFrom(RCurl,url.exists)
@@ -60,5 +63,6 @@ importFrom(stats,integrate)
 importFrom(stats,na.omit)
 importFrom(stats,pnorm)
 importFrom(stats,qnorm)
+importFrom(tibble,tibble)
 importFrom(utils,"?")
 importFrom(utils,help)

R/deform-package.R

@@ -7,7 +7,9 @@
 # @exportPattern "^[[:alpha:]]+"
 #
 # @importFrom Rcpp evalCpp
+#' @import dplyr
 #' @importFrom RCurl url.exists getURL
+#' @importFrom tibble tibble
 #' @importFrom Bessel BesselJ
 #' @importFrom graphics legend lines matplot points segments
 #' @importFrom stats integrate na.omit pnorm qnorm

R/poro-params.R

@@ -1,5 +1,9 @@
 #' Poroelastic parameters
 #' 
+#' @details 
+#' \code{\link{effstress}} Wang (2000), equation 2.52 
+#' \code{\link{poroelastic_stress_coefficient}} Wang (2000), equation 3.47
+#' 
 #' @param K numeric; the bulk modulus
 #' @param Ks numeric; the constituent bulk modulus (e.g., one from an unjacketed lab test)
 #' @param nu numeric; the drained Poisson's ratio  [0,0.5]
@@ -8,29 +12,38 @@
 #' @param mu numeric; the isotropic elastic shear modulus (keep at unity for 'normalized')
 #' @param diffusiv numeric; the isotropic hydraulic diffusivity (keep at unity for 'normalized')
 #' 
-#' @aliases poroelastic-parameters
+#' @references Wang, H.F. (2000)
 #' 
+#' @aliases poroelastic-parameters
+
 #' @export
 effstress <- function(nu, nuu, B){
-  3*(nuu - nu)/(1 - 2*nu)/(1 + nuu)/B
+  # e.g., Wang (2000) equation 2.52
+  3*(nuu - nu) / (1 - 2*nu) / (1 + nuu) / B
+}
+#' @rdname effstress
+#' @export
+poroelastic_stress_coefficient <- function(nu, nuu, B){
+  alpha <- effstress(nu, nuu, B)
+  (1 - 2*nu) * alpha / 2 / (1 - nu)
 }
 #' @rdname effstress
 #' @export
 biot_willis_effstress <- function(K, Ks){
-  1 - K/Ks
+  1 - K / Ks
 }
 #' @rdname effstress
 #' @export
 biot_compressibility <- function(nu, nuu, B, mu=1){
-  9*(1 - 2*nuu)*(nuu - nu)/2/mu/(1 - 2*nu)/(1 + nuu)^2/B^2
+  9 * (1 - 2*nuu) * (nuu - nu) / 2 / mu / (1 - 2*nu) / (1 + nuu)^2 / B^2
 }
 #' @rdname effstress
 #' @export
 darcy_conductivity <- function(nu, nuu, B, mu=1, diffusiv=1){
-  9*diffusiv*(1 - nuu)*(nuu - nu)/2/mu/(1 - nu)/(1 + nuu)^2/B^2
+  9 * diffusiv * (1 - nuu) * (nuu - nu) / 2 / mu / (1 - nu) / (1 + nuu)^2 / B^2
 }
 #' @rdname effstress
 #' @export
 lame_first <- function(nu, mu=1){
-  2*mu*nu/(1 - 2*nu)
+  2 * mu * nu / (1 - 2*nu)
 }

R/ricecleary76.R

@@ -0,0 +1,65 @@
+rc76_edge_disloc_poro <- function(Times, slip, radial_dist, angFromEast=30, 
+                                       diffusiv=1, 
+                                       nuu=0.33, nu=0.25, B = 0.6,
+                                       mu=30e9,
+                                       ...){
+  # Equation X of Rice and Cleary 1976, or equations 7.149 - 7.152 of Wang (2000)
+  data.table::inrange(radial_dist, 0, Inf)
+  stopifnot(length(slip) == 1)
+  stopifnot(length(radial_dist) == 1)
+  stopifnot(length(angFromEast) == 1)
+
+  zer <- Times <= 0
+  calc_times <- Times
+  calc_times[zer] <- NA
+
+  theta <- angFromEast * pi / 180
+  r <- radial_dist
+  geom_fac <- slip / r
+  nu_fac <- (nuu - nu) / (1 - nu) / (1 - nuu)
+
+  # scaling and temporal changes for all quantities
+  common_scaling <- (mu / 2 / pi) * nu_fac * geom_fac 
+  rsqtimes <- r^2 / 4 / diffusiv / calc_times
+  response <- 1 - exp(-1 * rsqtimes)
+
+  # pore pressure
+  eta <- poroelastic_stress_coefficient(nu, nuu, B)
+  #print(eta)
+  pp <- common_scaling * (sin(theta) / eta) * response
+
+  # stresses
+  srr <- common_scaling * sin(theta) * (-(1 - nu)/(nuu - nu) + rsqtimes * response)
+  srt <- common_scaling * cos(theta) * ((1 - nu)/(nuu - nu) - rsqtimes * response)
+  stt <- common_scaling * sin(theta) * (2*exp(-1 * rsqtimes) - (1 - nu)/(nuu - nu) - rsqtimes * response)
+
+  result <- tibble::tibble(Times=Times, sigma_tt = stt, sigma_rt = srt, sigma_rr = srr, pp = pp)
+
+  return(result)
+}
+
+library(dplyr)
+ti <- c(-100, -1, 0, 10**seq(-4, log10(6598), length.out=301))
+Mw <- 4.33
+M0 <- Mw2M0(Mw) #3.935501e+15 N*m for Mw 4.33
+bx <- 10^(-3.22 + 0.69 * Mw) #0.5
+rd <- 250
+mu <- 2e8
+Diffu <- 60
+rc <- rc76_edge_disloc_poro(ti, bx, rd, diffusiv=Diffu, mu=mu)
+nrc <- rc76_edge_disloc_poro(ti, bx, -rd, diffusiv=Diffu, mu=mu)
+all.equal(rc$pp, -nrc$pp)
+
+plot(pp / 1e3 + 5 ~ I(Times / 3600), rc, type='l', lwd=2, ylab='kPa', 
+     ylim=c(-8.8,36.9), yaxs='i', xaxs='i')
+axis(4)
+abline(v=1*60/3600, lty=3)
+rect(0,0,100,100)
+#lines(sigma_tt / 1e3 ~ Times, rc, lty=2)
+#lines(sigma_rt / 1e3 ~ Times, rc, lty=3)
+#lines(sigma_rr / 1e3 ~ Times, rc)
+#plot((sigma_rt - dplyr::first(na.omit(sigma_rt)))/1e3 ~ Times, rc, type='l')
+SrcArea <- M0 / (mu * bx)
+SrcRad <- sqrt(SrcArea / pi)
+print(c(Area=SrcArea, Rad=SrcRad))
+