fixed principal stresses and fault-plane values

NOBUILD/Dev/01.test_codes.R

@@ -34,7 +34,7 @@ n.t <- ceiling(1.1*n.)
 t. <- 10**seq(-1,2,length.out=n.t)
 
 td <- 1
-hw <- td #max(h.)/5
+hw <- td
 thk <- td/2
 
 ltd <- LambertTsai2020_Diffusive(x1=0, x2=1, y1=hw, resTopDepth=td, resThickness=thk, resDiffusivity=1, Time=t.)
@@ -46,37 +46,21 @@ with(ltd, matplot(pars$Time, Def, type='l', lty=c(1,2,1,3,5,1,3,1), col=c(1,1,2,
 with(ltu, matplot(pars$x1, Def, ylim=c(-3,3), type='l', lty=c(1,2,1,3,5,1,3,1), col=c(1,1,2,2,2,3,3,4)))
 abline(v=hw*c(-1,1), col='grey')
 
+#library(magrittr)
 
-#plane stress principal stresses
-PSPS(ltd)
-PSPS(ltu)
-
-rbind(
-fault_normal_stress(ltd),
-fault_shear_stress(ltd),
-mean_stress(ltd),
-max_shear_principal(ltd),
-max_shear(ltd)
-)
-
-rbind(
-fault_normal_stress(ltu),
-fault_shear_stress(ltu),
-mean_stress(ltu),
-max_shear_principal(ltu),
-max_shear(ltu)
-)
-
-cd0 <- cfs_isoporo(ltd)
-plot(t., cd0, type='l', ylim=max(abs(cd0), na.rm=TRUE)*c(-1,1))
-lines(t., cfs_isoporo(ltd, theta=45), lty=5)
-lines(t., cfs_isoporo(ltd, theta=90), lty=2)
-cu0 <- cfs_isoporo(ltu)
-plot(h., cu0, type='l', ylim=max(abs(cu0), na.rm=TRUE)*c(-1,1))
-lines(h., cfs_isoporo(ltu, theta=45), lty=5)
-lines(h., cfs_isoporo(ltu, theta=90), lty=2)
-
-#principal_stresses(ltu); stop()
+examp <- function(){
+
+	fs <- fault_stresses(ltu, Beta=30, is.deg=TRUE)
+	cu.nf <- cfs_isoporo(ltu, Beta.deg=30, fric=fric, Skemp=skemp)
+	cu.ss <- cfs_isoporo(ltu, Beta.deg=5, fric=fric, Skemp=skemp)
+
+	matplot(h., abs(fs[,7:10]), type='l', log='xy', ylim=range(abs(c(cu.ss,cu.nf)), na.rm=TRUE))
+	lines(h., abs(cu.nf), lwd=1, type='h', col=ifelse(cu<0,4,2))
+	lines(h., abs(cu.nf), lwd=1.2)
+	lines(h., abs(cu.ss), lty=2, lwd=2)
+
+	plot(h., cu.nf / cu.ss, log='')
+}
 
 # do by depth
 do_calc_by_depth <- function(d){
@@ -88,17 +72,20 @@ d. <- d.[d. < max(d.)/1.5]
 Du <- lapply(d., do_calc_by_depth)
 
 #str(Du,1)
+library(pbapply)
+#pboptions(type = "timer")
 
-#cu <- sapply(Du, cfs_isoporo, theta=60, fric=fric, Skemp=skemp) #, verbose=FALSE)
+cu <- pbapply::pbsapply(Du, cfs_isoporo, Beta.deg=30, fric=fric, Skemp=skemp, verbose=FALSE)
 #cu <- sapply(Du, mean_stress)
 #cu <- sapply(Du, max_shear)
-cu <- sapply(Du, max_shear_principal)
+#cu <- sapply(Du, max_shear_principal)
 #cu <- sapply(Du, ext1) # OK
 #cu <- sapply(Du, ext2) # OK
 #cu <- sapply(Du, shear) # out of wack for negative distances -- added kluge
 #str(cu,1)
 
-cu[!is.finite(cu)] <- NA
+#cu[!is.finite(cu)] <- NA
+
 xy <- list(x = h./td, y = d./td)
 Cu <- c(xy, list(z = cu))
 laCu <- c(xy, list(z = log10(abs(cu))))
@@ -110,5 +97,5 @@ fields::image.plot(laCu, asp=1, ylim=rev(range(Cu[['y']], na.rm=TRUE))) #, zlim=
 #image(Cupos, add=TRUE, col=adjustcolor(c('white',NA), alpha=0.5))
 #contour(lCup, add=TRUE)
 #contour(lCun, lty=3, add=TRUE)
-contour(laCu, add=TRUE)
-rect(-hw, (td + thk)/td, hw, td/td, col='white', border='grey')
+contour(Cu, add=TRUE)
+rect(-hw/td, (td + thk)/td, hw/td, td/td, col='white', border='grey')

NOBUILD/Dev/LambertCodes.R

@@ -215,106 +215,76 @@ ext2.LT.diffusive <- ext2.LT.uniform <- function(x, angle=FALSE){
 	Def[,'S22']
 }
 
-principal_stresses <- function(x, ...) UseMethod('principal_stresses')
-principal_stresses.LT.diffusive <- principal_stresses.LT.uniform <- function(x){
-	S1 <- ext1(x)
-	S2 <- ext2(x)
-	S12 <- shear(x)
-	S <- cbind(S1, S12, S12, S2)
+ang_max_shear <- function(x, ...) UseMethod('ang_max_shear')
+ang_max_shear.LT.diffusive <- ang_max_shear.LT.uniform <- function(x){
+	s11 <- ext1(x)
+	s22 <- ext2(x)
+	s12 <- shear(x)
+	atan((s22 - s11)/(2 * s12)) / 2
+}
+
+planestress_principal_stresses <- function(x, ...) UseMethod('planestress_principal_stresses')
+planestress_principal_stresses.LT.diffusive <- planestress_principal_stresses.LT.uniform <- function(x){
+	s1 <- ext1(x)
+	s2 <- ext2(x)
+	s12 <- shear(x)
+	zer <- 0*s1
+	S <- cbind(s1, s12, zer, s12, s2, zer, zer, zer, zer)
+	S[!is.finite(S)] <- NA
 	.get_eig <- function(si){
-		sim <- matrix(si, 2)
+		sim <- matrix(si, 3)
+		navals <- rep(NA, 3)
 		vals <- if (any(is.na(sim) | !is.finite(sim))){
-			c(NA, NA)
+			navals
 		} else {
-			ei <- eigen(sim)
-			ei[['values']]
+			ei <- try(eigen(sim))
+			if (inherits(ei, 'try-error')){
+				navals
+			} else {
+				zapsmall(ei[['values']])
+			}
 		}
-		vals <- matrix(vals, ncol=2)
+		vals <- matrix(vals, ncol=3)
 
 		vals
 	}
-	s1s2 <- t(apply(S, 1, .get_eig))
-	colnames(s1s2) <- c('PS1','PS2')
-	s1s2
-}
-
-max_shear_principal <- function(x, ...) UseMethod('max_shear_principal')
-max_shear_principal.LT.diffusive <- max_shear_principal.LT.uniform <- function(x){
-	s1s2 <- principal_stresses(x)
-	max_shear_principal.default(S1=s1s2[,'PS1'], S3=s1s2[,'PS2'])
-}
-max_shear_principal.default <- function(S1, S3) (S1 - S3)/2
-
-max_shear <- function(x, ...) UseMethod('max_shear')
-max_shear.LT.diffusive <- max_shear.LT.uniform <- function(x, angle=FALSE){
-	p <- PSPS(x)
-	p[,ifelse(angle,'tau','tau_theta')]
+	s1s2s3 <- t(apply(S, 1, .get_eig))
+	colnames(s1s2s3) <- c('PS1','PS2','PS3')
+	cbind(s11=s1, s12=s12, s22=s2, s1s2s3)
 }
 
-mean_stress <- function(x, ...) UseMethod('mean_stress')
-mean_stress.LT.diffusive <- mean_stress.LT.uniform <- function(x){
-	Def <- x[['Def']]
-	mean_stress.default(S1=ext1(x), S3=ext2(x))
-}
-mean_stress.default <- function(S1, S3) (S1 + S3)/2
-
-PSPS <- function(x, ...) UseMethod('PSPS')
+MSMS <- function(x, ...) UseMethod('MSMS')
 
-PSPS.LT.diffusive <- PSPS.LT.uniform <- function(x, ...){
-	PSPS.default(s11=ext1(x), s12=shear(x), s22=ext2(x))
+MSMS.LT.diffusive <- MSMS.LT.uniform <- function(x, ...){
+	PS <- planestress_principal_stresses(x)
+	S1 <- PS[,'PS1']
+	S3 <- PS[,'PS3']
+	maxshear <- (S1 - S3) / 2
+	meanstress <- (S1 + S3) / 2
+	cbind(PS, tau_max = maxshear, Smean = meanstress)
 }
 
-PSPS.default <- function(s11, s12, s22){ # plane-stress principal stresses (i.e., 2D Mohr's circle)
-	ms <- mean_stress.default(s11, s22)
-	maxs <- max_shear_principal.default(s11, s22)
-
-	shearstress <- sqrt(maxs^2 + s12^2)
-	shearangle <- atan(-maxs / s12) / 2
-
-	S1 <- ms + sqrt(maxs^2 + s12^2)
-	S2 <- ms - sqrt(maxs^2 + s12^2)
-
-	cbind(S1=S1, S2=S2, tau = shearstress, tau_theta = shearangle)
-}
-
-#--------------------------------------------------------------
-# Fault-plane projection computations
-
-fault_normal_stress <- function(x, theta=0, ...) UseMethod('fault_normal_stress')
-
-fault_normal_stress.LT.diffusive <- fault_normal_stress.LT.uniform <- function(x, theta=0, ...){
-	fault_normal_stress.default(s11=ext1(x), s12=shear(x), s22=ext2(x), theta = theta, ...)
-}
-
-fault_normal_stress.default <- function(s11, s12, s22, theta, is.deg=TRUE){
-	# normal stress acting on a plane oriented by theta radians from the 1-2 coord system
-	ms <- mean_stress.default(s11, s22)
-	maxs <- max_shear_principal.default(s11, s22)
-	if (is.deg) theta <- theta * pi / 180
-	ms + maxs * cos(2*theta) + s12 * sin(2*theta)
-}
-
-fault_shear_stress <- function(x, theta=0, ...) UseMethod('fault_shear_stress')
-
-fault_shear_stress.LT.diffusive <- fault_shear_stress.LT.uniform <- function(x, theta=0, ...){
-	fault_shear_stress.default(s11=ext1(x), s12=shear(x), s22=ext2(x), theta = theta, ...)
-}
+fault_stresses <- function(x, theta=0, ...) UseMethod('fault_stresses')
 
-fault_shear_stress.default <- function(s11, s12, s22, theta, is.deg=TRUE){
-	# shear stress acting on a plane oriented by theta radians from the 1-2 coord system
-	maxs <- max_shear_principal.default(s11, s22)
-	if (is.deg) theta <- theta * pi / 180
-	# use double-angle identity: sin(2x) = 2*sin(x)cos(x)
-	-maxs*sin(2*theta) + s12*(cos(theta)^2 - sin(theta)^2)
+fault_stresses.LT.diffusive <- fault_stresses.LT.uniform <- function(x, Beta=45, is.deg=TRUE, ...){
+	# Beta is the angle between the plane and the greatest principal stress (S1, S3 < S2 < S1)
+	if (is.deg) Beta <- Beta * pi / 180
+	msms <- MSMS(x)
+	twopsi <- pi - 2*Beta
+	tau <- msms[,'tau_max']
+	Tau <- tau * sin(twopsi)
+	Sig <- msms[,'Smean'] + tau * cos(twopsi)
+	return(cbind(msms, Tau = Tau, Snorm = Sig))
 }
 
 # (cfs_isoporo.defaultin beeler)
-cfs_isoporo.LT.diffusive <- cfs_isoporo.LT.uniform <- function(x, theta=0, fric=0.65, Skemp=0.6, verbose=TRUE, ...){
-	if (verbose) message('CFS {S1-to-fault angle = ', theta, '°, friction = ', fric, "} for isotropic undrained poroelastic response {B = ", Skemp, "}")
-	tauf <- fault_shear_stress(x, theta)
-	Snorm <- fault_normal_stress(x, theta)
-	Smean <- mean_stress(x)
-	cfsip <- cfs_isoporo.default(tau=tauf, mu=fric, Snormal = Snorm, Smean = Smean, B=Skemp)
+cfs_isoporo.LT.diffusive <- cfs_isoporo.LT.uniform <- function(x, Beta.deg, fric, Skemp, verbose=TRUE, ...){
+	if (verbose) message('CFS {S1-to-fault angle = ', Beta.deg, '°, friction = ', fric, "} for isotropic undrained poroelastic response {B = ", Skemp, "}")
+	fs <- fault_stresses(x, Beta=Beta.deg, is.deg=TRUE)
+	Tau <- fs[,'Tau']
+	Snorm <- fs[,'Snorm']
+	Smean <- fs[,'Smean']
+	cfsip <- cfs_isoporo.default(tau=Tau, mu=fric, Snormal = Snorm, Smean = Smean, B=Skemp)
 	zapsmall(cfsip)
 }
 

NOBUILD/Dev/beeler_cfs.R

@@ -3,15 +3,15 @@ pp_isotropic <- function(Smean, B){
 	B * Smean
 }
 
-effective_normal_stress <- function(Snormal, pp) Snormal - pp
+effective_stress <- function(S, pp) S - pp
 
 cfs_basic <- function(tau, mu, Snormal, pp){
-	sn_eff <- effective_normal_stress(Snormal, pp)
+	sn_eff <- effective_stress(Snormal, pp) # effective normal stress
 	tau - mu*(sn_eff)
 }
 
 cfs_isoporo <- function(x, ...) UseMethod("cfs_isoporo")
-cfs_isoporo.default <- function(tau, mu, Snormal, Smean, B){
+cfs_isoporo.default <- function(tau, mu, Snormal, Smean, B, ...){
 	pp.iso <- pp_isotropic(Smean, B)
 	cfs_basic(tau, mu, Snormal, pp.iso)
 }