data/tools/2d-event-writer.R

+source('tools/support/findQuakes.R')
+
+finalTime              <- 1000 # s
+convergenceVelocity    <- 5e-5 # m/s
+criticalVelocity       <- 1000e-6 + 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("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
+
+  velocityProxy<- h5file['/frictionalBoundary/velocity']
+
+  quakes <- findQuakes(criticalVelocity, velocityProxy,
+                       indices = 1:dim(velocityProxy)[1], 1)
+
+  basalCoordinates <- h5file['/frictionalBoundary/coordinates'][]
+
+  maximumVelocities<- maxVelocity(velocityProxy[,,1])
+  displacement     <- h5file['/frictionalBoundary/displacement']
+  numVertices      <- dim(displacement)[[2]]
+  for (quakeID in seq(nrow(quakes))) {
+    qb <- quakes[[quakeID,'beginningIndex']]
+    qe <- quakes[[quakeID,'endingIndex']]
+    localSlippingTimes <- abs(velocityProxy[qb:qe,,1][,,1]) > criticalVelocity
+
+    qd   <- displacement[qb:qe,,1]
+    slip <- vector(mode='numeric', length=numVertices)
+    for (i in seq(numVertices)) {
+      if (any(localSlippingTimes[,i])) {
+        begs <- positiveStarts(localSlippingTimes[,i])
+        ends <- negativeStarts(localSlippingTimes[,i])
+        slip[i]<- sum(qd[ends,i,] - qd[begs,i,])
+      }
+    }
+    quakes[quakeID,'peakSlip']     <- max(abs(slip))
+    quakes[quakeID,'peakSlipRate'] <- max(maximumVelocities[qb:qe])
+
+    maxRuptureWidth <- 0
+    for (ts in seq(dim(localSlippingTimes)[[1]])) {
+      st <- localSlippingTimes[ts,]
+      if (!any(st))
+        next;
+
+      slippingCoordinates <- basalCoordinates[st,1] # x-coordinate only
+      maxRuptureWidth <- max(maxRuptureWidth, diff(range(slippingCoordinates)))
+    }
+    quakes[quakeID,'ruptureWidth'] <- maxRuptureWidth
+  }
+  h5::h5close(h5file)
+
+  quakes$beginning <- realTime[quakes$beginningIndex]
+  quakes$ending    <- realTime[quakes$endingIndex]
+
+  write.csv(quakes[c('beginning','ending',
+                     'peakSlip','peakSlipRate',
+                     'ruptureWidth')],
+            row.names = FALSE, quote = FALSE,
+            file = file.path(directories[['output']],
+                             paste.(pasteColon('events', basedir), 'csv')))
+}

data/tools/2d-fpi-tolerance.R

+rfpitols <- c('1e-7', '2e-7', '3e-7', '5e-7',
+              '10e-7', '20e-7', '30e-7', '50e-7',
+              '100e-7', '200e-7', '300e-7', '500e-7',
+              '1000e-7', '2000e-7', '3000e-7', '5000e-7',
+              '10000e-7', '20000e-7', '30000e-7', '50000e-7',
+              '100000e-7')
+numEntries <- length(rfpitols)
+
+data <- data.frame(row.names = rfpitols,
+                   tol = rep(NA, numEntries),
+                   time = rep(NA, numEntries),
+                   fpi = rep(NA, numEntries),
+                   mg = rep(NA, numEntries))
+
+directories <- ini::read.ini('config.ini')$directories
+dir.create(directories[['output']], recursive=TRUE, showWarnings=FALSE)
+for (rfpitol in rfpitols) {
+  basedir <- paste('rfpitol', rfpitol, sep='=')
+  dir     <- file.path(directories[['simulation']],
+                       '2d-lab-fpi-tolerance', basedir)
+  h5file  <- h5::h5file(file.path(dir, 'output.h5'), 'r')
+
+  data[rfpitol,'tol'] <- as.numeric(rfpitol)
+
+  relativeTimeProxy <- h5file['/relativeTime']
+  relativeTimeLen <- relativeTimeProxy@dim
+  data[rfpitol,'time'] <- relativeTimeProxy[relativeTimeProxy@dim]
+
+  ## FIXME: why do we drop the first entry?
+  fixedPointIterationsProxy <- h5file["/iterations/fixedPoint/total"]
+  fixedPointIterationsLen <- fixedPointIterationsProxy@dim
+  data[rfpitol,'fpi'] <- sum(fixedPointIterationsProxy[2:fixedPointIterationsLen])
+
+  multiGridIterationsProxy <- h5file["/iterations/multiGrid/total"]
+  multiGridIterationsLen <- multiGridIterationsProxy@dim
+  data[rfpitol,'mg'] <- sum(multiGridIterationsProxy[2:multiGridIterationsLen])
+  h5::h5close(h5file)
+}
+
+write.csv(data, file.path(directories[['output']], 'fpi-data.csv'),
+          row.names = FALSE, quote = FALSE)

data/tools/2d-performance.R

+source('tools/support/findQuakes.R')
+
+finalTime           <- 1000 # s
+convergenceVelocity <- 5e-5 # m/s
+
+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
+
+  velocityProxy<- h5file['/frictionalBoundary/velocity']
+
+  ## 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, velocityProxy,
+                       indices = 1:dim(velocityProxy)[1], 1)
+  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
+  plotIndices<- which(plotMask)
+  plotIndices<- c(plotIndices[1]-1,plotIndices,plotIndices[length(plotIndices)]+1)
+
+  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')))
+
+  timeWindow = realTime[plotIndices]
+
+  relativeTau          <- h5file['relativeTimeIncrement'][]
+  fixedPointIterations <- h5file['/iterations/fixedPoint/final'][]
+  multiGridIterations  <- h5file['/iterations/multiGrid/final'][]
+  write.csv(data.frame(time = timeWindow,
+                       timeIncrement = finalTime * relativeTau[plotIndices],
+                       fixedPointIterations = fixedPointIterations[plotIndices],
+                       multiGridIterations = multiGridIterations[plotIndices]),
+            file.path(directories[['output']],
+                      paste.(pasteColon('2d-performance', basedir), 'csv')),
+            row.names = FALSE, quote = FALSE)
+  h5::h5close(h5file)
+}

data/tools/2d-performance.py

+import ConfigParser as cp
+import os
+import numpy as np
+import csv
+import h5py
+
+from support.find_quakes import find_quakes
+
+NBODIES = 2
+FINAL_TIME = 1000  # s
+FINAL_VELOCITY = 1e-5  # m/s
+THRESHOLD_VELOCITY = 0.5*FINAL_VELOCITY  # 1000e-6 + FINAL_VELOCITY
+
+TANGENTIAL_COORDS = 1
+
+# read config ini
+config = cp.ConfigParser()
+config_path = os.path.join('config.ini')
+config.read(config_path)
+sim_path = config.get('directories', 'simulation')
+exp_path = config.get('directories', 'experiment')
+out_path = config.get('directories', 'output')
+
+# create output directory
+out_dir = os.path.join(out_path)
+if not os.path.exists(out_dir):
+    os.mkdir(out_dir)
+
+h5path = os.path.join(sim_path)
+h5file = h5py.File(os.path.join(h5path, 'output.h5'), 'r')
+
+# read time
+relative_time = np.array(h5file['relativeTime'])
+real_time = relative_time * FINAL_TIME
+
+
+  velocityProxy<- h5file['/frictionalBoundary/velocity']
+
+  ## 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, velocityProxy,
+                       indices = 1:dim(velocityProxy)[1], 1)
+  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
+  plotIndices<- which(plotMask)
+  plotIndices<- c(plotIndices[1]-1,plotIndices,plotIndices[length(plotIndices)]+1)
+
+  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')))
+
+  timeWindow = realTime[plotIndices]
+
+  relativeTau          <- h5file['relativeTimeIncrement'][]
+  fixedPointIterations <- h5file['/iterations/fixedPoint/final'][]
+  multiGridIterations  <- h5file['/iterations/multiGrid/final'][]
+  write.csv(data.frame(time = timeWindow,
+                       timeIncrement = finalTime * relativeTau[plotIndices],
+                       fixedPointIterations = fixedPointIterations[plotIndices],
+                       multiGridIterations = multiGridIterations[plotIndices]),
+            file.path(directories[['output']],
+                      paste.(pasteColon('2d-performance', basedir), 'csv')),
+            row.names = FALSE, quote = FALSE)
+  h5::h5close(h5file)

data/tools/2d-velocity-contours.R

+source('tools/support/findQuakes.R')
+source('tools/support/writeContours.R')
+
+finalTime           <- 1000 # s
+convergenceVelocity <- 5e-5 # m/s
+
+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
+
+  velocityProxy<- h5file['/frictionalBoundary/velocity']
+
+  basalCoordinates <- h5file['/frictionalBoundary/coordinates'][]
+  basalTrenchDistance <- basalCoordinates[,1]
+  perm                <- order(basalTrenchDistance)
+  sortedBasalTrenchDistance <- basalTrenchDistance[perm]
+
+  {
+    ## 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, velocityProxy,
+                         indices = 1:dim(velocityProxy)[1], 1)
+    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)
+    plotMask   <- (realTime > relaxedTime[[1]]) & (realTime < relaxedTime[[2]])
+    plotIndices<- which(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')))
+
+    printlevels <- c('1000','100','10','1')
+    levels      <- 1e-6 * as.numeric(printlevels) + convergenceVelocity
+    ret <- contourLines(realTime[plotIndices],
+                        sortedBasalTrenchDistance,
+                        abs(velocityProxy[plotIndices,perm,1][,,1]),
+                        levels = levels)
+
+    for (i in seq(printlevels))
+      writeContours(ret, levels[[i]],
+                    file.path(directories[['output']],
+                              paste.(pasteColon('2d-velocity-contours',
+                                                'threequakes', basedir, 'level',
+                                                printlevels[[i]]), 'tex')))
+  }
+  {
+    ## 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(300e-6 + convergenceVelocity, velocityProxy,
+                         indices = 1:dim(velocityProxy)[1], 1)
+    quakes$beginning <- realTime[quakes$beginningIndex]
+    quakes$ending    <- realTime[quakes$endingIndex]
+    quakes$duration  <- quakes$ending - quakes$beginning
+    numQuakes        <- nrow(quakes)
+    quake            <- quakes[numQuakes,]
+    relaxedTime      <-
+      c(quake[['beginning']] - 0.9*(quake[['ending']] - quake[['beginning']]),
+        quake[['ending']]    + 0.1*(quake[['ending']] - quake[['beginning']]))
+    plotMask <- (realTime > relaxedTime[[1]]) & (realTime < relaxedTime[[2]])
+    plotIndices<- which(plotMask)
+
+    printlevels <- c('3000','1000','300','100','30','10','3','1')
+    levels      <- 1e-6 * as.numeric(printlevels) + convergenceVelocity
+    ret <- contourLines(realTime[plotIndices],
+                        sortedBasalTrenchDistance,
+                        abs(velocityProxy[plotIndices,perm,1][,,1]),
+                        levels = levels)
+
+    for (i in seq(printlevels))
+      writeContours(ret, levels[[i]],
+                    file.path(directories[['output']],
+                              paste.(pasteColon('2d-velocity-contours',
+                                                'zoom', basedir, 'level',
+                                                printlevels[[i]]), 'tex')))
+  }
+  h5::h5close(h5file)
+}

data/tools/2d_event_writer.py

+import ConfigParser as cp
+import os
+import numpy as np
+import csv
+import h5py
+
+from support.maximum import maximum
+from support.norm import norm
+from support.find_quakes import find_quakes
+from support.slip_beginnings import slip_beginnings
+from support.slip_endings import slip_endings
+from support.max_distance import max_distance
+
+
+NBODIES = 2
+FINAL_TIME = 1000  # s
+FINAL_VELOCITY = 1e-5  # m/s
+THRESHOLD_VELOCITY = 0.5*FINAL_VELOCITY  # 1000e-6 + FINAL_VELOCITY
+
+TANGENTIAL_COORDS = 1
+
+# read config ini
+config = cp.ConfigParser()
+config_path = os.path.join('config.ini')
+config.read(config_path)
+sim_path = config.get('directories', 'simulation')
+exp_path = config.get('directories', 'experiment')
+out_path = config.get('directories', 'output')
+
+# create output directory
+out_dir = os.path.join(out_path)
+if not os.path.exists(out_dir):
+    os.mkdir(out_dir)
+
+h5path = os.path.join(sim_path)
+h5file = h5py.File(os.path.join(h5path, 'output.h5'), 'r')
+
+# read time
+relative_time = np.array(h5file['relativeTime'])
+real_time = relative_time * FINAL_TIME
+
+for body_ID in range(NBODIES):
+    body = 'body' + str(body_ID)
+
+    if body not in h5file:
+        continue
+
+    # read data
+    coordinates = np.array(h5file[body + '/coordinates'])
+    displacement = np.array(h5file[body + '/displacement'])
+    velocity = np.array(h5file[body + '/velocity'])
+
+    num_vertices = displacement.shape[1]
+    velocity_norm = norm(velocity[:, :, TANGENTIAL_COORDS:])
+    maximum_velocity = maximum(velocity_norm)
+
+    [quake_starts, quake_ends] = find_quakes(THRESHOLD_VELOCITY,
+                                             maximum_velocity)
+
+    quakes = {}
+    for quake_ID in range(len(quake_starts)):
+        quake_start = int(quake_starts[quake_ID])
+        quake_end = int(quake_ends[quake_ID])
+
+        quake = {}
+        quake['start'] = real_time[quake_start]
+        quake['end'] = real_time[quake_end]
+
+        local_slipping_times = velocity_norm[quake_start:quake_end, :] \
+            > THRESHOLD_VELOCITY
+        quake_displacement = displacement[quake_start:quake_end, :, :]
+        slip = np.zeros(num_vertices)
+
+        for i in range(num_vertices):
+            if any(local_slipping_times[:, i]):
+                starts = slip_beginnings(local_slipping_times[:, i])
+                ends = slip_endings(local_slipping_times[:, i])
+                slip[i] = np.linalg.norm(quake_displacement[ends, i, :]
+                                         - quake_displacement[starts, i, :])
+
+        quake['peakSlip'] = max(slip)
+        quake['peakSlipRate'] = max(maximum_velocity[quake_start:quake_end])
+
+        max_rupture_width = 0
+        for time_step in range(local_slipping_times.shape[0]):
+            slipping_time = local_slipping_times[time_step, :]
+            if not any(slipping_time):
+                continue
+
+            slipping_coords = coordinates[slipping_time] \
+                + displacement[quake_start + time_step, slipping_time]
+            if len(slipping_coords) > 1:
+                pair = np.array(max_distance(slipping_coords))
+                max_rupture_width = max(max_rupture_width,
+                                        np.linalg.norm(pair[0] - pair[1]))
+
+        quake['ruptureWidth'] = max_rupture_width
+        quakes[quake_ID] = quake
+
+    # output quake data to csv file
+    csv_columns = quakes[0].keys()
+    csv_file_path = os.path.join(out_path, 'events' + str(body_ID) + '.csv')
+    try:
+        with open(csv_file_path, 'w') as csv_file:
+            writer = csv.DictWriter(csv_file, fieldnames=csv_columns)
+            writer.writeheader()
+            for quake_ID in quakes:
+                writer.writerow(quakes[quake_ID])
+    except IOError:
+        print('I/O error')
+
+h5file.close()

data/tools/2d_event_writer.py~

+import ConfigParser as cp
+import os
+import numpy as np
+import csv
+import h5py
+
+from support.maximum import maximum
+from support.norm import norm
+from support.find_quakes import find_quakes
+from support.slip_beginnings import slip_beginnings
+from support.slip_endings import slip_endings
+from support.max_distance import max_distance
+
+
+NBODIES = 2
+FINAL_TIME = 1000  # s
+FINAL_VELOCITY = 1e-5  # m/s
+THRESHOLD_VELOCITY = 0.5*FINAL_VELOCITY  # 1000e-6 + FINAL_VELOCITY
+
+TANGENTIAL_COORDS = 1
+
+# read config ini
+config = cp.ConfigParser()
+config_path = os.path.join('config.ini')
+config.read(config_path)
+sim_path = config.get('directories', 'simulation')
+exp_path = config.get('directories', 'experiment')
+out_path = config.get('directories', 'output')
+
+# create output directory
+out_dir = os.path.join(out_path)
+if not os.path.exists(out_dir):
+    os.mkdir(out_dir)
+
+h5path = os.path.join(sim_path)
+h5file = h5py.File(os.path.join(h5path, 'output.h5'), 'r')
+
+# read time
+relative_time = np.array(h5file['relativeTime'])
+real_time = relative_time * FINAL_TIME
+
+for body_ID in range(NBODIES):
+    body = 'body' + str(body_ID)
+
+    if body not in h5file:
+        continue
+
+    # read data
+    coordinates = np.array(h5file[body + '/coordinates'])
+    displacement = np.array(h5file[body + '/displacement'])
+    velocity = np.array(h5file[body + '/velocity'])
+
+    num_vertices = displacement.shape[1]
+    velocity_norm = norm(velocity[:, :, TANGENTIAL_COORDS:])
+    maximum_velocity = maximum(velocity_norm)
+
+    [quake_starts, quake_ends] = find_quakes(THRESHOLD_VELOCITY,
+                                             maximum_velocity)
+
+    quakes = {}
+    for quake_ID in range(len(quake_starts)):
+        quake_start = int(quake_starts[quake_ID])
+        quake_end = int(quake_ends[quake_ID])
+
+        quake = {}
+        quake['start'] = real_time[quake_start]
+        quake['end'] = real_time[quake_end]
+
+        local_slipping_times = velocity_norm[quake_start:quake_end, :] \
+            > THRESHOLD_VELOCITY
+        quake_displacement = displacement[quake_start:quake_end, :, :]
+        slip = np.zeros(num_vertices)
+
+        for i in range(num_vertices):
+            if any(local_slipping_times[:, i]):
+                starts = slip_beginnings(local_slipping_times[:, i])
+                ends = slip_endings(local_slipping_times[:, i])
+                slip[i] = np.linalg.norm(quake_displacement[ends, i, :]
+                                         - quake_displacement[starts, i, :])
+
+        quake['peakSlip'] = max(slip)
+        quake['peakSlipRate'] = max(maximum_velocity[quake_start:quake_end])
+
+        max_rupture_width = 0
+        for time_step in range(local_slipping_times.shape[0]):
+            slipping_time = local_slipping_times[time_step, :]
+            if not any(slipping_time):
+                continue
+
+            slipping_coords = coordinates[slipping_time] \
+                + displacement[quake_start + time_step, slipping_time]
+            if len(slipping_coords) > 1:
+                pair = np.array(max_distance(slipping_coords))
+                max_rupture_width = max(max_rupture_width,
+                                        np.linalg.norm(pair[0] - pair[1]))
+
+        quake['ruptureWidth'] = max_rupture_width
+        quakes[quake_ID] = quake
+
+    # output quake data to csv file
+    csv_columns = quakes[0].keys()
+    csv_file_path = os.path.join(out_path, 'events' + str(body_ID) + '.csv')
+    try:
+        with open(csv_file_path, 'w') as csv_file:
+            writer = csv.DictWriter(csv_file, fieldnames=csv_columns)
+            writer.writeheader()
+            for quake_ID in quakes:
+                writer.writerow(quakes[quake_ID])
+    except IOError:
+        print('I/O error')
+
+h5file.close()

data/tools/3d-performance.R

+source('tools/support/findQuakes.R')
+
+finalTime           <- 1000 # s
+convergenceVelocity <- 5e-5 # m/s
+
+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 very low velocity here.
+  quakes <- findQuakes(1e-6 + convergenceVelocity, velocityProxy,
+                       ## Note: We only look at the last 2000 timesteps here because
+                       ## we're only interested in the last few events. This
+                       ## dramatically reduces RAM usage and runtime.
+                       indices = seq(dim(velocityProxy)[1]-2000,
+                                     dim(velocityProxy)[1]), c(1,3))
+  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
+  plotIndices<- which(plotMask)
+  plotIndices<- c(plotIndices[1]-1,plotIndices,plotIndices[length(plotIndices)]+1)
+
+  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')))
+
+  timeWindow = realTime[plotIndices]
+
+  relativeTau          <- h5file['relativeTimeIncrement'][]
+  fixedPointIterations <- h5file['/iterations/fixedPoint/final'][]
+  multiGridIterations  <- h5file['/iterations/multiGrid/final'][]
+  write.csv(data.frame(time = timeWindow,
+                       timeIncrement = finalTime * relativeTau[plotIndices],
+                       fixedPointIterations = fixedPointIterations[plotIndices],
+                       multiGridIterations = multiGridIterations[plotIndices]),
+            file.path(directories[['output']],
+                      paste.(pasteColon('3d-performance', basedir), 'csv')),
+            row.names = FALSE, quote = FALSE)
+  h5::h5close(h5file)
+}

data/tools/3d-velocity-contours.R

+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)
+}

data/tools/comparison:lab-sim.R

+paste.    <- function(...) paste(..., sep='.')
+pasteColon<- function(...) paste(..., sep=':')
+
+basedir <- 'rfpitol=100e-7'
+
+directories <- ini::read.ini('config.ini')$directories
+labdata <- within(
+  read.table(unz(file.path(directories[['experiment']],
+                           'B_Scale-model-earthquake-data.zip'),
+                 'scaleEQs_12-31_catalogue.txt'),
+             sep='\t', header=FALSE, comment.char='%',
+             col.names=c('frame',
+                         'meanSlipNature', 'meanSlipLab',
+                         'peakSlipNature', 'peakSlipLab',
+                         'ruptureWidthNature', 'ruptureWidthLab',
+                         'ignored', 'ignored', 'ignored', 'ignored', 'ignored')),
+  {
+    recurrenceLab <- c(NA, diff(frame))/10 # (10Hz camera: 10 frames ~ 1s)
+    meanSlipLab   <- meanSlipLab / 1e6     # micro m -> m
+    peakSlipLab   <- peakSlipLab / 1e6})   # micro m -> m
+## NB: We skip the first two quakes here, for compatibility with my old code
+##     Maybe skipping the first was intentional so that each quake would have
+##     a recurrence time. But skipping the second was certainly an accident
+labdata <- labdata[3:nrow(labdata),]
+
+simdata <- read.csv(file.path(directories[['output']],
+                              paste.(pasteColon('events', basedir), 'csv')))
+
+report <- function(lab, sim, quantity) {
+  write.table(lab, file.path(directories[['output']],
+                             paste.(pasteColon('boxplot-data', 'lab',
+                                               quantity), 'tex')),
+              row.names=FALSE, col.names=FALSE)
+  write.table(sim, file.path(directories[['output']],
+                             paste.(pasteColon('boxplot-data', 'simulation',
+                                               basedir, quantity), 'tex')),
+              row.names=FALSE, col.names=FALSE)
+}
+report(labdata$recurrenceLab, diff(simdata$beginning), 'recurrence')
+report(labdata$ruptureWidthLab, simdata$ruptureWidth, 'ruptureWidth')
+# meters -> millimeters
+report(labdata$peakSlipLab * 1000, simdata$peakSlip*1000, 'peakSlip')

data/tools/config.ini

+[directories]
+simulation = /storage/mi/podlesny/software/dune/build-debug/dune-tectonic/src/strikeslip
+experiment = ~/group/publications/2016-RosenauCorbiDominguezRudolfRitterPipping
+output     = generated

data/tools/generate-2d.bash

+#!/usr/bin/env bash
+
+set -e
+
+rr() {
+    echo "$(date)" Running: Rscript $@
+    Rscript --vanilla --default-packages=grDevices,methods,stats,utils tools/$@
+}
+
+# contours
+rr 2d-velocity-contours.R
+
+# dip
+rr 2d-dip-contours.R
+
+# iterations / adaptivity
+rr 2d-performance.R
+
+# box plot (one half)
+rr 2d-event-writer.R
+
+# fpi data
+rr 2d-fpi-tolerance.R
+
+date

data/tools/generate-3d.bash

+#!/usr/bin/env bash
+
+set -e
+
+rr() {
+    echo "$(date)" Running: Rscript $@
+    Rscript --vanilla --default-packages=grDevices,methods,stats,utils tools/$@
+}
+
+# performance
+rr 3d-performance.R
+
+# velocity contours
+rr 3d-velocity-contours.R
+
+date

data/tools/generate-others.bash

+#!/usr/bin/env bash
+
+set -e
+
+rr() {
+    echo "$(date)" Running: Rscript $@
+    Rscript --vanilla --default-packages=grDevices,methods,stats,utils tools/$@
+}
+
+# box plot (one half)
+rr comparison:lab-sim.R
+
+date

data/tools/generated/events1.csv

+start,peakSlipRate,end,ruptureWidth,peakSlip
+20.086059570312422,5.7045853762396602e-06,20.086669921874922,0,0.0
+20.090332031249918,5.6219949973608878e-06,20.092163085937418,0,0.0
+20.094604492187415,5.8412336903395239e-06,20.095825195312415,0,0.0
+20.099487304687411,5.3352579051826403e-06,20.100097656249908,0,0.0
+20.143737792968619,5.4285740142335498e-06,20.144958496093619,0,0.0
+20.148620605468615,5.6370561043864526e-06,20.149536132812365,0,0.0
+20.353851318359059,6.0525078256123743e-06,20.354461669921559,0,0.0
+20.358734130859055,5.1806538408711459e-06,20.359344482421555,0,0.0
+20.371856689452795,5.6547173619088325e-06,20.373077392577795,0,0.0
+20.425262451171495,6.4545930360194781e-06,20.428314208983991,0.088385894124722503,1.4298690605481336e-08
+20.479583740233945,5.8310594282664788e-06,20.481414794921445,0,0.0
+20.483245849608942,5.2537138171449348e-06,20.485076904296442,0,6.2000101383992609e-09
+20.531158447265149,5.2316818926703734e-06,20.531768798827649,0,0.0
+20.544586181640138,5.4092184353134343e-06,20.545196533202635,0,0.0
+20.593109130858842,6.190190600109647e-06,20.595550537108842,0,0.0
+20.608673095702581,5.5618674305029011e-06,20.609283447265081,0,0.0
+20.625152587890067,5.4920875934131527e-06,20.625762939452567,0,0.0
+20.630035400390064,5.1740512343599518e-06,20.630645751952564,0,0.0
+20.63369750976506,5.7392806837644885e-06,20.63491821289006,0,0.0

data/tools/support/__init__.py

+ 

data/tools/support/find_quakes.py

+from support.slip_beginnings import slip_beginnings
+from support.slip_endings import slip_endings
+
+
+def find_quakes(threshold_velocity, maximum_velocities):
+    slipping_times = maximum_velocities > threshold_velocity
+
+    quake_starts = slip_beginnings(slipping_times)
+    quake_ends = slip_endings(slipping_times)
+
+    # remove incomplete quakes
+    min_len = min(len(quake_starts), len(quake_ends))
+    quake_ends = quake_ends[0:min_len]
+    quake_starts = quake_starts[0:min_len]
+
+    return [quake_starts, quake_ends]

data/tools/support/find_quakes.py~

+from support.slip_beginnings import slip_beginnings
+from support.slip_endings import slip_endings
+
+
+def find_quakes(threshold_velocity, maximum_velocities):
+    slipping_times = maximum_velocities > threshold_velocity
+    print(slipping_times)
+
+    quake_starts = slip_beginnings(slipping_times)
+    quake_ends = slip_endings(slipping_times)
+
+    # remove incomplete quakes
+    min_len = min(len(quake_starts), len(quake_ends))
+    quake_ends = quake_ends[0:min_len]
+    quake_starts = quake_starts[0:min_len]
+
+    return [quake_starts, quake_ends]