Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
source('tools/support/findQuakes.R')
source('tools/support/writeContours.R')
Rcpp::sourceCpp('tools/support/trapezoidal.cpp')
finalTime <- 1000 # s
specialTrenchDistances <- c(0.15,0.30,0.45) # m
convergenceVelocity <- 5e-5 # m/s
last <- function(x) x[[length(x)]]
paste. <- function(...) paste(..., sep='.')
pasteColon<- function(...) paste(..., sep=':')
directories <- ini::read.ini('config.ini')$directories
dir.create(directories[['output']], recursive=TRUE, showWarnings=FALSE)
for (basedir in c("rfpitol=100e-7")) {
dir <- file.path(directories[['simulation']],
'2d-lab-fpi-tolerance', basedir)
h5file <- h5::h5file(file.path(dir, 'output.h5'), 'r')
relativeTime <- h5file['relativeTime'][]
realTime <- finalTime * relativeTime
calcMask <- relativeTime >= 0.7
calcTime <- realTime[calcMask]
calcDuration <- last(calcTime) - calcTime[[1]]
calcRange <- range(which(calcMask))
calcLength <- sum(calcMask)
## We are interested in an enlarged time range around actual events here,
## (and no other quantities!) hence we pass a very low velocity here.
quakes <- findQuakes(1e-6 + convergenceVelocity,
h5file['/frictionalBoundary/velocity'][,,],
indices = calcRange[1]:calcRange[2])
quakes$beginning <- realTime[quakes$beginningIndex]
quakes$ending <- realTime[quakes$endingIndex]
quakes$duration <- quakes$ending - quakes$beginning
numQuakes <- nrow(quakes)
relaxedTime <- extendrange(c(quakes[[numQuakes-2,'beginning']],
quakes[[numQuakes, 'ending']]), f=0.02)
threeQuakeTimeMask <- (realTime > relaxedTime[[1]]) & (realTime < relaxedTime[[2]])
plotMask <- threeQuakeTimeMask
plotTime <- realTime[plotMask]
plotRange <- range(which(plotMask))
plotLength <- sum(plotMask)
write(relaxedTime[[1]],
file.path(directories[['output']],
paste.(pasteColon('timeframe', 'min', 'threequakes', basedir),
'tex')))
write(relaxedTime[[2]],
file.path(directories[['output']],
paste.(pasteColon('timeframe', 'max', 'threequakes', basedir),
'tex')))
surfaceCoordinates <- h5file['/surface/coordinates'][]
surfaceTrenchDistance <- surfaceCoordinates[,1] / cos(atan(.27))
perm <- order(surfaceTrenchDistance)
displacement <- h5file['/surface/displacement']
# This is the displacement imposed through the Dirichlet condition at the last node
displacementOffset <- displacement[calcRange[1]:calcRange[2],last(perm),1]
balancedSurfaceDisplacement <- matrix(nrow = plotLength, ncol = displacement@dim[2])
for (k in 1:displacement@dim[2]) {
d <- displacement[calcRange[1]:calcRange[2],k,1:2]
d[,,1] <- d[,,1] - displacementOffset
# We're in a tilted coordinate system
dv <- -sin(atan(.27))*d[,,1] + cos(atan(.27))*d[,,2]
meanVertSurfaceDisplacement <- trapezoidal(calcTime, dv) / calcDuration
balancedSurfaceDisplacement[,k] <- (dv - meanVertSurfaceDisplacement)[
(plotRange[1] - calcRange[1] + 1):(plotRange[2] - calcRange[1] + 1)
]
}
## Interpolate velocity to special points from surrounding velocities
{
data <- matrix(nrow=plotLength, ncol=1+length(specialTrenchDistances))
data[,1] <- plotTime
for (k in seq(plotLength)) {
interpolant <- approxfun(surfaceTrenchDistance[perm],
balancedSurfaceDisplacement[k,perm])
for (i in seq(specialTrenchDistances)) {
specialTrenchDistance <- specialTrenchDistances[[i]]
data[k,i+1] <- interpolant(specialTrenchDistance)
}
}
singleDataName <- file.path(directories[['output']],
paste.(pasteColon('dip-single-points', basedir),
'csv'))
write.csv(data, singleDataName, row.names = FALSE)
}
h5::h5close(h5file)
printlevels <- c('-40','-20','-10','-5','-2.5','0','2.5','5','10','20','40')
levels <- as.numeric(printlevels) * 1e-6
ret <- contourLines(plotTime,
surfaceTrenchDistance[perm],
balancedSurfaceDisplacement[,perm],
levels=levels)
for (i in seq(printlevels))
writeContours(ret, levels[[i]],
file.path(directories[['output']],
paste.(pasteColon('2d-dip-contours', basedir,
'level', printlevels[[i]]),
'tex')))
}