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
source('tools/support/findQuakes.R')
source('tools/support/writeContours.R')
finalTime <- 1000 # s
convergenceVelocity <- 5e-5 # m/s
printlevels <- c('1','2','3','5',
'10','20','30','50',
'100','200','300','500',
'1000')
criticalVelocities <- 1e-6*as.numeric(printlevels) + convergenceVelocity
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("rtol=1e-5_diam=1e-2")) {
dir <- file.path(directories[['simulation']],
'3d-lab', basedir)
h5file <- h5::h5file(file.path(dir, 'output.h5'), 'r')
relativeTime <- h5file['relativeTime'][]
realTime <- finalTime * relativeTime
velocityProxy<- h5file['/frictionalBoundary/velocity']
## We are interested in an enlarged time range around actual events here,
## (and no other quantities!) hence we pass a rather low velocity here.
quakes <- findQuakes(200e-6 + convergenceVelocity, velocityProxy,
## Note: We only look at the last 1000 timesteps here because
## we're only interested in the last few events. This
## dramatically reduces RAM usage and runtime.
indices = seq(dim(velocityProxy)[1]-1000,
dim(velocityProxy)[1]), c(1,3))
quake <- quakes[nrow(quakes)-1,] # Q: why did we not need the -1 in julia?
weakPatchGridVelocityProxy <- h5file["/weakPatchGrid/velocity"]
stepSize <- 30 ## note: should/might differ by time/spatial resolution
ts <- seq(quake$beginningIndex, quake$endingIndex, by=stepSize)
dd <- data.frame(timeSteps = ts,
times = realTime[ts],
timeOffsets = realTime[ts] - realTime[quake$beginningIndex])
fname = pasteColon('3d-velocity-contours', basedir)
write.csv(dd, row.names = FALSE, quote = FALSE,
file = file.path(directories[['output']],
paste.(pasteColon(fname, 'times'), 'csv')))
weakPatchGridXCoordinates <- h5file["/weakPatchGrid/xCoordinates"][]
weakPatchGridZCoordinates <- h5file["/weakPatchGrid/zCoordinates"][]
for (t in ts) {
m <- weakPatchGridVelocityProxy[t,,,]
s <- sqrt(m[,,,1]^2 + m[,,,3]^2)
ret <- contourLines(weakPatchGridXCoordinates, weakPatchGridZCoordinates, s,
level=criticalVelocities)
for (i in seq(printlevels))
writeContours(ret, criticalVelocities[[i]],
file.path(directories[['output']],
paste.(pasteColon(fname,
'step', t,
'level', printlevels[[i]]),
'tex')))
}
h5::h5close(h5file)
}