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

+#/home/joscha/software/dune/build-release/dune-tectonic/src/foam/output/test/
+#/home/joscha/software/dune/build-debug/dune-tectonic/src/foam/output/test/
+#/home/joscha/Desktop/strikeslip/viscosity-1e3/
+ #/home/joscha/Downloads/
+
+[directories]
+simulation = /home/joscha/software/dune/build-release/dune-tectonic/src/foam/output/pipping-2013-euler/
+experiment = ~/group/publications/2016-RosenauCorbiDominguezRudolfRitterPipping
+output     = generated

data/tools/debug.py

+import configparser as cp
+import os
+import numpy as np
+import csv
+import h5py
+import matplotlib.pyplot as plt
+
+from debug.outliers import outliers
+from debug.friction import truncated_friction
+from debug.state import aging_law
+from debug.diffplot import diffplot
+
+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
+
+from support.iterations import iterations
+
+NBODIES = 2
+FINAL_TIME = 15  # s
+FINAL_VELOCITY = 1e-5  # m/s
+THRESHOLD_VELOCITY = 0.5*FINAL_VELOCITY  # 1000e-6 + FINAL_VELOCITY
+
+TANGENTIAL_COORDS = 1
+
+# friction params
+params = {
+    'L'  : 1e-5,
+    'V0' : 1e-6,
+    'mu0': 0.6,
+    'a'  : 0.010,
+    'b'  : 0.015
+}
+
+# read config ini
+config = cp.ConfigParser()
+config_path = os.path.join('tools/config.ini')
+config.read(config_path)
+sim_path = config.get('directories', 'simulation')
+exp_path = config.get('directories', 'experiment')
+out_path = config.get('directories', 'output')
+
+# read hdf5 output file
+h5path = os.path.join(sim_path)
+h5file = h5py.File(os.path.join(h5path, 'output.h5'), 'r')
+
+print(list(h5file.keys()))
+print(list(h5file['body1'].keys()))
+
+# read time
+relative_time = np.array(h5file['relativeTime'])
+relative_tau = np.array(h5file['relativeTimeIncrement'])
+relative_tau = np.delete(relative_tau, 0)
+
+real_time = relative_time * FINAL_TIME
+real_tau = relative_tau * FINAL_TIME
+
+print(len(relative_time))
+
+for body_ID in range(NBODIES):
+    body = 'body' + str(body_ID)
+
+    if body not in h5file:
+        continue
+
+    # velocity data
+    v = abs(np.array(h5file[body + '/velocity']))
+
+    # statistics
+    avg_v = np.average(v[:,:,TANGENTIAL_COORDS], axis=1)
+    min_v = np.min(v[:,:,TANGENTIAL_COORDS], axis=1)
+    max_v = np.max(v[:,:,TANGENTIAL_COORDS], axis=1)
+
+    # plot
+    plt.figure(1)
+    plt.subplot(311)
+    plt.plot(min_v, color='gray', linestyle='--')
+    plt.plot(avg_v, color='black', linestyle='-')
+    plt.plot(max_v, color='gray', linestyle='--')
+    plt.ylabel('slip rate')
+    #-------------------------
+
+    # state
+    states = np.array(h5file[body + '/state'])
+    states_calc = aging_law(params, states[0], v[:,:,TANGENTIAL_COORDS], real_tau)
+
+    # statistics
+    avg_states = np.average(states, axis=1)
+    avg_states_calc  = np.average(states_calc, axis=1)
+    min_states = np.min(states, axis=1)
+    max_states = np.max(states, axis=1)
+
+    # plot
+    plt.subplot(312)
+    plt.plot(min_states, color='gray', linestyle='--')
+    plt.plot(avg_states, color='black', linestyle='-')
+    plt.plot(max_states, color='gray', linestyle='--')
+    plt.plot(avg_states_calc, color='red', linestyle='-')
+    plt.ylabel('state')
+    #-------------------------
+
+    # friction coefficient
+    friction_coeff = np.array(h5file[body + '/coefficient'])
+
+    #weighted_normal_stress = np.array(h5file[body + '/weightedNormalStress'])
+    friction_coeff_calc = truncated_friction(params, v[:,:,TANGENTIAL_COORDS], states_calc)
+
+    # statistics
+    avg_friction_coeff = np.average(friction_coeff, axis=1)
+    avg_friction_coeff_calc = np.average(friction_coeff_calc, axis=1)
+    min_friction_coeff = np.min(friction_coeff, axis=1)
+    max_friction_coeff = np.max(friction_coeff, axis=1)
+
+    outliers_friction_coeff = outliers(avg_friction_coeff)
+
+    # plot
+    plt.subplot(313)
+    plt.plot(min_friction_coeff, color='gray', linestyle='--')
+    plt.plot(avg_friction_coeff, color='black', linestyle='-')
+    plt.plot(max_friction_coeff, color='gray', linestyle='--')
+    plt.plot(avg_friction_coeff_calc, color='red', linestyle='-')
+    plt.plot(outliers_friction_coeff[0], outliers_friction_coeff[1], color='red', marker='+')
+    plt.ylabel('friction coefficient')
+    #-------------------------
+
+    diffplot(friction_coeff[1], friction_coeff_calc[1], 'diff friction coeff')
+
+    iterations(h5file, FINAL_TIME)
+
+    plt.show()
+
+h5file.close()
\ No newline at end of file

data/tools/debug/diffplot.py

+import numpy as np
+import matplotlib.pyplot as plt
+
+def diffplot(v1, v2, title):
+    v_diff = v1 - v2
+    t = np.linspace(0, len(v_diff)-1, len(v_diff), endpoint=True)
+
+    # plot
+    fig = plt.figure()
+    ax = fig.add_subplot(1, 1, 1)
+    ax.plot(t, v1, color='black', linestyle='--')
+    ax.plot(t, v2, color='gray', linestyle='--')
+    ax.plot(t, v_diff, color='red', linestyle='-')
+    ax.set_ylabel(title)
+    #-------------------------
+
+    fig.canvas.draw()
\ No newline at end of file

data/tools/debug/friction.py

+import numpy as np
+
+def truncated_friction(params, v, alpha):
+   vmin = params['V0'] / np.exp( ( params['mu0'] + params['b'] * np.array(alpha)) / params['a'] )
+   clipped_v = (v/vmin).clip(1)
+   return params['a'] * np.log(clipped_v)
+
+   #return (params['a'] * np.log(v/params['V0']) - params['mu0'] - (params['b']/params['a'])*np.array(alpha) ).clip(0)
\ No newline at end of file

data/tools/debug/outliers.py

+import numpy as np
+
+def outliers(data):
+    index = [i for i,x in enumerate(data) if np.abs(x)==np.inf]
+    val = np.zeros(len(index))
+    return [index, val]
\ No newline at end of file

data/tools/debug/state.py

+import numpy as np
+
+def lift_singularity(c, x):
+    def local_lift(y): 
+        if (y <= 0):
+            return -c
+        else:
+            return -np.expm1(c * y) / y
+    
+    return np.array([local_lift(y) for y in x])
+
+
+def aging_law(params, initial_alpha, v, time_steps):
+        #auto tangentVelocity = velocity_field[localToGlobal_[i]];
+        #tangentVelocity[0] = 0.0;
+        #double const V = tangentVelocity.two_norm();
+        
+    alpha = [initial_alpha]
+
+    for i,tau in enumerate(time_steps):
+        mtoL = -tau / params['L']
+        next_alpha = np.log(np.exp(alpha[-1] + v[i] * mtoL) + params['V0'] * lift_singularity(mtoL, v[i]))
+        alpha.append(next_alpha)
+
+    return alpha
\ No newline at end of file

data/tools/elias.py

+import configparser as cp
+import os
+import numpy as np
+import csv
+import h5py
+import matplotlib.pyplot as plt
+
+from debug.outliers import outliers
+from debug.friction import truncated_friction
+from debug.state import aging_law
+from debug.diffplot import diffplot
+
+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
+
+from support.iterations import iterations
+from support.friction_stats import friction_stats
+
+def build_patch(coords, percentage):
+    x_coords = coords[:, 0]
+    xmin = np.min(x_coords)
+    xmax = np.max(x_coords)
+    delta_x = (1 - percentage)*(xmax - xmin)/2
+
+    xmin = xmin + delta_x
+    xmax = xmax - delta_x
+
+    return [i for i in range(len(x_coords)) if x_coords[i]>=xmin and x_coords[i]<=xmax]
+
+FINAL_TIME = 1000  # s
+FINAL_VELOCITY = 5e-5  # m/s
+THRESHOLD_VELOCITY = 0.5*FINAL_VELOCITY  # 1000e-6 + FINAL_VELOCITY
+
+TANGENTIAL_COORDS = 0
+
+# friction params
+params = {
+    'L'  : 1e-5,
+    'V0' : 1e-6,
+    'mu0': 0.6,
+    'a'  : 0.010,
+    'b'  : 0.015
+}
+
+# read config ini
+config = cp.ConfigParser()
+config_path = os.path.join('tools/config.ini')
+config.read(config_path)
+sim_path = config.get('directories', 'simulation')
+exp_path = config.get('directories', 'experiment')
+out_path = config.get('directories', 'output')
+
+# read hdf5 output file
+h5path = os.path.join(sim_path)
+h5file = h5py.File(os.path.join(h5path, 'output.h5'), 'r')
+
+print(list(h5file.keys()))
+print(list(h5file['frictionalBoundary'].keys()))
+
+iterations(h5file, FINAL_TIME)
+
+coords = np.array(h5file['frictionalBoundary/coordinates'])
+patch = build_patch(coords, 0.05)
+
+friction_stats(h5file, 0, FINAL_TIME, [], [0, 50], TANGENTIAL_COORDS)
+    
+plt.show()
+
+h5file.close()
\ No newline at end of file

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