2016-PippingKornhuberRosenauOncken to ed8f469d3f519f3473544c091707f37c5d774025 · podlesny / dune-tectonic

Some changes are not shown.

For a faster browsing experience, only 110 of 357 files are shown. Download one of the files below to see all changes.

CMakeLists.txt

+1 −1

Original line number	Diff line number	Diff line
		@@ -17,7 +17,7 @@ include(DuneMacros)

		# start a dune project with information from dune.module
		dune_project()

		dune_enable_all_packages()
		find_package(HDF5 COMPONENTS C REQUIRED)

		add_subdirectory("src")

data/2d-dip-contours-performance.tikz

0 → 100644

+91 −0

Original line number	Diff line number	Diff line
		\pgfplotsarraynew\heights{%
		-40\\-20\\-10\\-5\\-2.5\\%
		0\\%
		2.5\\5\\10\\20\\40\\%
		}
		\def\fnamebase{2d-dip-contours}
		\def\fname{rfpitol=100e-7}

		\newread\threequakesminread
		\openin\threequakesminread=generated/timeframe:min:threequakes:\fname.tex
		\read\threequakesminread to \threequakesmin
		\closein\threequakesminread

		\newread\threequakesmaxread
		\openin\threequakesmaxread=generated/timeframe:max:threequakes:\fname.tex
		\read\threequakesmaxread to \threequakesmax
		\closein\threequakesmaxread

		\begin{tikzpicture}[trim axis group left, trim axis group right]
		\begin{groupplot}[
		xmin=\threequakesmin, xmax=\threequakesmax, max space between ticks=40pt,
		x tick label style={ /pgf/number format/1000 sep={} },
		tick label style={font=\footnotesize},
		label style={font=\small},
		group style={
		x descriptions at = edge bottom,
		group size=1 by 4,
		vertical sep=0cm
		},
		width=12cm,
		enlargelimits=false]
		\nextgroupplot[
		semithick,
		height = 6.5cm,
		/pgf/number format/1000 sep={},
		ymax=0.7,
		colormap/jet,
		y tick label style={
		/pgf/number format/.cd,
		fixed, fixed zerofill, precision=2,
		/tikz/.cd
		},
		tick label style={font=\footnotesize},
		label style={font=\small},
		legend style={font=\small,
		at={(1.05,1)},
		anchor=north west,
		fill=none},
		ylabel = distance from trench, y unit = m,
		legend cell align=right,
		contour/labels=false,
		contour prepared]
		\pgfplotsinvokeforeach{1,...,9} { % level 0 and 10 are empty
		\pgfplotscolormapaccess[0:10]{#1}{jet}
		\def\TEMP{\definecolor{mycolor#1}{rgb}}
		\expandafter\TEMP\expandafter{\pgfmathresult}
		\def\fnameX{generated/\fnamebase:\fname:level:%
		\pgfplotsarrayvalueofelem#1\of\heights.tex}
		\addplot[contour prepared={draw color=mycolor#1}, forget plot] table \fnameX;
		\addlegendimage{line legend,color=mycolor#1}
		\addlegendentry{%
		\SI{\pgfplotsarrayvalueofelem#1\of\heights}{\micro\meter}}
		};
		\nextgroupplot[
		semithick,
		height = 3.5cm,
		ylabel style={align=center}, ylabel = time step\\size, y unit = s,
		ytick={1e-3,1e-2,1e-1},
		ymax = 1, ymin = 1e-4,
		ymode = log]
		\addplot[const plot, mark=none] table[col sep=comma, y=timeIncrement]
		{generated/2d-performance:\fname.csv};
		\nextgroupplot[
		semithick,
		height = 3.5cm,
		ylabel style={align=center}, ylabel=fixed-point\\iterations,
		ytick={2,4,6,8},
		ymin=0, ymax=10]
		\addplot[const plot, mark=none] table[col sep=comma, y=fixedPointIterations]
		{generated/2d-performance:\fname.csv};
		\nextgroupplot[
		semithick,
		height = 3.5cm,
		xlabel = time, x unit = s,
		ylabel style={align=center}, ylabel=multigrid\\iterations,
		ytick={5,10,15,20},
		ymin=0, ymax=25]
		\addplot[const plot, mark=none] table[col sep=comma, y=multiGridIterations]
		{generated/2d-performance:\fname.csv};
		\end{groupplot}
		\end{tikzpicture}

data/2d-dip-single-points.tikz

0 → 100644

+51 −0

Original line number	Diff line number	Diff line
		\def\fnamebase{dip-single-points}
		\def\fname{rfpitol=100e-7}

		\newread\threequakesminread
		\openin\threequakesminread=generated/timeframe:min:threequakes:\fname.tex
		\read\threequakesminread to \threequakesmin
		\closein\threequakesminread

		\newread\threequakesmaxread
		\openin\threequakesmaxread=generated/timeframe:max:threequakes:\fname.tex
		\read\threequakesmaxread to \threequakesmax
		\closein\threequakesmaxread

		\begin{tikzpicture}[trim axis left, trim axis right]
		\begin{axis}[
		xmin=\threequakesmin, xmax=\threequakesmax, max space between ticks=40pt,
		/pgf/number format/1000 sep={},
		x tick label style={
		/pgf/number format/.cd,
		1000 sep={},
		/tikz/.cd
		},
		y tick label style={
		/pgf/number format/.cd,
		fixed, fixed zerofill, precision=2,
		/tikz/.cd
		},
		tick label style={font=\footnotesize},
		label style={font=\small},
		legend style={font=\small, at={(1.05,1)},
		anchor=north west,
		fill=none},
		legend entries={
		\SI{15}{\centi\meter},
		\SI{30}{\centi\meter},
		\SI{45}{\centi\meter}
		},
		ylabel style={align=center}, ylabel=vertical surface\\displacement, y unit = m,
		change y base, y SI prefix=micro,
		xlabel = time, x unit=s,
		width = 12cm,
		height = 4cm,
		semithick]
		\addplot[width=2pt] table[col sep=comma, x index = 0, y index=1]
		{generated/\fnamebase:\fname.csv};
		\addplot[width=2pt, dashed] table[col sep=comma, x index = 0, y index=2]
		{generated/\fnamebase:\fname.csv};
		\addplot[width=2pt, dotted] table[col sep=comma, x index = 0, y index=3]
		{generated/\fnamebase:\fname.csv};
		\end{axis}
		\end{tikzpicture}

data/2d-effort-over-tolerance.tikz

0 → 100644

+17 −0

Original line number	Diff line number	Diff line
		\begin{tikzpicture}[trim axis left, trim axis right]
		\begin{axis}[
		semithick,
		width = 12cm,
		height = 4cm,
		xmode=log,
		%ymin = 0,
		% scaled y ticks=base 10:0,
		xlabel = fixed point tolerance,
		ylabel style={align=center}, ylabel = multigrid\\iterations, % (multigrid steps)
		extra x ticks = 1e-5,
		extra x tick labels = { time-stepping\\tolerance },
		extra x tick style = { align=center, grid = major, ticklabel pos=right }
		]
		\addplot[mark=+] table[col sep=comma, y=mg] {generated/fpi-data.csv};
		\end{axis}
		\end{tikzpicture}

data/2d-velocity-contours-threequakes.tikz

0 → 100644

+61 −0

Original line number	Diff line number	Diff line
		\pgfplotsarraynew\contourlevels{%
		1\\3\\
		10\\30\\
		100\\300\\
		1000\\3000\\
		}

		\def\fname{threequakes:rfpitol=100e-7}

		\newread\threequakesminread
		\openin\threequakesminread=generated/timeframe:min:\fname.tex
		\read\threequakesminread to \threequakesmin
		\closein\threequakesminread

		\newread\threequakesmaxread
		\openin\threequakesmaxread=generated/timeframe:max:\fname.tex
		\read\threequakesmaxread to \threequakesmax
		\closein\threequakesmaxread

		\begin{tikzpicture}[trim axis left, trim axis right]
		\begin{axis}[
		xmin=\threequakesmin, xmax=\threequakesmax, max space between ticks=40pt,
		ymax=0.85,
		colormap/jet,
		x tick label style={
		/pgf/number format/.cd,
		1000 sep={},
		/tikz/.cd
		},
		y tick label style={
		/pgf/number format/.cd,
		fixed, fixed zerofill, precision=2,
		/tikz/.cd
		},
		tick label style={font=\footnotesize},
		label style={font=\small},
		legend style={font=\small,
		at={(1.05,1)},
		anchor=north west,
		fill=none},
		ylabel = distance from trench, y unit = m,
		xlabel = time, x unit=s,
		width = 12cm, height = 6.5cm,
		legend cell align=right,
		contour/labels=false,
		enlargelimits=false,
		semithick]

		\pgfplotsinvokeforeach{0,2,4,6} {
		\pgfplotscolormapaccess[0:7]{#1}{jet}
		\def\TEMP{\definecolor{mycolor#1}{rgb}}
		\expandafter\TEMP\expandafter{\pgfmathresult}
		\def\fnameX{generated/2d-velocity-contours:\fname:level:%
		\pgfplotsarrayvalueofelem#1\of\contourlevels.tex}
		\addplot[contour prepared={draw color=mycolor#1}, forget plot] table \fnameX;
		\addlegendimage{line legend,color=mycolor#1}
		\addlegendentry{%
		\SI{\pgfplotsarrayvalueofelem#1\of\contourlevels}{\micro\meter/\second}}
		};
		\end{axis}
		\end{tikzpicture}

data/2d-velocity-contours-zoom.tikz

0 → 100644

+51 −0

Original line number	Diff line number	Diff line
		\pgfplotsarraynew\contourlevels{%
		1\\3\\
		10\\30\\
		100\\300\\
		1000\\3000\\
		}

		\def\fname{zoom:rfpitol=100e-7}
		\begin{tikzpicture}[trim axis left, trim axis right]
		\begin{axis}[
		max space between ticks=40pt,
		ymax=0.7,
		colormap/jet,
		x tick label style={
		/pgf/number format/.cd,
		fixed, fixed zerofill, precision=2,
		1000 sep={},
		/tikz/.cd
		},
		y tick label style={
		/pgf/number format/.cd,
		fixed, fixed zerofill, precision=2,
		/tikz/.cd
		},
		tick label style={font=\footnotesize},
		label style={font=\small},
		legend style={font=\small,
		at={(1.05,1)},
		anchor=north west,
		fill=none},
		ylabel = distance from trench, y unit = m,
		xlabel = time, x unit=s,
		width = 12cm, height = 6.5cm,
		legend cell align=right,
		enlargelimits=false,
		contour/labels=false,
		semithick]

		\pgfplotsinvokeforeach{0,...,6} { % level 7 is empty
		\pgfplotscolormapaccess[0:7]{#1}{jet}
		\def\TEMP{\definecolor{mycolor#1}{rgb}}
		\expandafter\TEMP\expandafter{\pgfmathresult}
		\def\fnameX{generated/2d-velocity-contours:\fname:level:%
		\pgfplotsarrayvalueofelem#1\of\contourlevels.tex}
		\addplot[contour prepared={draw color=mycolor#1}, forget plot] table \fnameX;
		\addlegendimage{line legend,color=mycolor#1}
		\addlegendentry{%
		\SI{\pgfplotsarrayvalueofelem#1\of\contourlevels}{\micro\meter/\second}}
		};
		\end{axis}
		\end{tikzpicture}

data/3d-performance.tikz

0 → 100644

+57 −0

Original line number	Diff line number	Diff line
		\def\fname{rtol=1e-5_diam=1e-2}

		\newread\threequakesminread
		\openin\threequakesminread=generated/timeframe:min:threequakes:\fname.tex
		\read\threequakesminread to \threequakesmin
		\closein\threequakesminread

		\newread\threequakesmaxread
		\openin\threequakesmaxread=generated/timeframe:max:threequakes:\fname.tex
		\read\threequakesmaxread to \threequakesmax
		\closein\threequakesmaxread

		\begin{tikzpicture}[trim axis group left, trim axis group right]
		\begin{groupplot}[
		xmin=\threequakesmin, xmax=\threequakesmax,
		max space between ticks=40pt,
		x tick label style={
		/pgf/number format/.cd,
		1000 sep={},
		/tikz/.cd
		},
		tick label style={font=\footnotesize},
		label style={font=\small},
		group style={
		x descriptions at = edge bottom,
		group size=1 by 3,
		vertical sep=0cm
		},
		height=3.5cm, width=12cm,
		enlargelimits=false]

		\nextgroupplot[
		semithick,
		ylabel style={align=center}, ylabel = time step\\size, y unit = s,
		ytick={1e-3,1e-2,1e-1},
		ymax = 1, ymin = 1e-4,
		ymode = log]
		\addplot[const plot, mark=none] table[col sep=comma, y=timeIncrement]
		{generated/3d-performance:\fname.csv};
		\nextgroupplot[
		semithick,
		ylabel style={align=center}, ylabel = fixed-point\\iterations,
		ytick={2,4,6,8},
		ymin=0, ymax=10]
		\addplot[const plot, mark=none] table[col sep=comma, y=fixedPointIterations]
		{generated/3d-performance:\fname.csv};
		\nextgroupplot[
		semithick,
		xlabel = time, x unit = s,
		ylabel style={align=center}, ylabel=multigrid\\iterations,
		ytick={5,10,15,20},
		ymin=0, ymax=25
		]
		\addplot[const plot, mark=none] table[col sep=comma, y=multiGridIterations]
		{generated/3d-performance:\fname.csv};
		\end{groupplot}
		\end{tikzpicture}

data/3d-velocity-contours.tikz

0 → 100644

+156 −0

Original line number	Diff line number	Diff line
		\def\fname{generated/3d-velocity-contours:rtol=1e-5_diam=1e-2}
		\pgfplotstableread[col sep=comma]{\fname:times.csv}\myloadedtable

		\pgfplotsarraynew\contourlevels{%
		1\\2\\3\\5\\%
		10\\20\\30\\50\\%
		100\\200\\300\\500\\%
		1000\\
		}

		\begin{tikzpicture}[trim axis group left, trim axis group right]
		\begin{groupplot}[
		group style={
		y descriptions at = edge left,
		group size=6 by 1,
		horizontal sep=0cm
		},
		ymin=-0.30, ymax= 0.30,
		enlarge x limits=false,
		colormap/jet,
		y tick label style={
		/pgf/number format/.cd,
		fixed, fixed zerofill, precision=2,
		/tikz/.cd
		},
		ytick = {-0.30,-0.20,-0.10,0.00,0.10,0.20,0.30},
		tick label style={font=\footnotesize},
		label style={font=\small},
		width = 3.5cm, height = 6cm,
		enlargelimits=false,
		contour/labels=false,
		%
		groupplot xlabel={distance from trench [\si{\meter}]},
		]

		\nextgroupplot[semithick, xlabel = {
		\pgfplotstablegetelem{0}{times}\of\myloadedtable%
		$t_0 \approx
		\SI[round-mode=places,round-precision=0]{\pgfplotsretval}{\second}$
		},
		ylabel = width, y unit = m]
		\draw[color=gray!20] (0.162533,-0.30) -- (0.162533,+0.30); % X
		\draw[color=gray!20] (0.362533+0.05,-0.30) -- (0.362533-0.05,+0.30); % Y
		\pgfplotsinvokeforeach{0,...,12} {
		\pgfplotscolormapaccess[0:14]{#1}{jet}
		\def\TEMP{\definecolor{mycolor#1}{rgb}}
		\expandafter\TEMP\expandafter{\pgfmathresult}

		\pgfplotstablegetelem{0}{timeSteps}\of\myloadedtable
		\edef\fnameX{\fname:%
		step:\pgfplotsretval:%
		level:\pgfplotsarrayvalueofelem#1\of\contourlevels.tex}
		\addplot[contour prepared={draw color=mycolor#1}] table \fnameX;
		};

		\nextgroupplot[semithick, xlabel = {
		\pgfplotstablegetelem{1}{timeOffsets}\of\myloadedtable%
		$t_0 +
		\SI[round-mode=places,round-precision=2]{\pgfplotsretval}{\second}$
		}]
		\draw[color=gray!20] (0.162533,-0.30) -- (0.162533,+0.30); % X
		\draw[color=gray!20] (0.362533+0.05,-0.30) -- (0.362533-0.05,+0.30); % Y
		\pgfplotsinvokeforeach{0,...,12} { % level 13 and 14 are empty
		\pgfplotscolormapaccess[0:14]{#1}{jet}
		\def\TEMP{\definecolor{mycolor#1}{rgb}}
		\expandafter\TEMP\expandafter{\pgfmathresult}
		\pgfplotstablegetelem{1}{timeSteps}\of\myloadedtable
		\edef\fnameX{\fname:%
		step:\pgfplotsretval:%
		level:\pgfplotsarrayvalueofelem#1\of\contourlevels.tex}
		\addplot[contour prepared={draw color=mycolor#1}] table \fnameX;
		};

		\nextgroupplot[semithick, xlabel = {
		\pgfplotstablegetelem{2}{timeOffsets}\of\myloadedtable%
		$t_0 +
		\SI[round-mode=places,round-precision=2]{\pgfplotsretval}{\second}$
		}]
		\draw[color=gray!20] (0.162533,-0.30) -- (0.162533,+0.30); % X
		\draw[color=gray!20] (0.362533+0.05,-0.30) -- (0.362533-0.05,+0.30); % Y
		\pgfplotsinvokeforeach{0,...,12} { % level 13 and 14 are empty
		\pgfplotscolormapaccess[0:14]{#1}{jet}
		\def\TEMP{\definecolor{mycolor#1}{rgb}}
		\expandafter\TEMP\expandafter{\pgfmathresult}
		\pgfplotstablegetelem{2}{timeSteps}\of\myloadedtable
		\edef\fnameX{\fname:%
		step:\pgfplotsretval:%
		level:\pgfplotsarrayvalueofelem#1\of\contourlevels.tex}
		\addplot[contour prepared={draw color=mycolor#1}] table \fnameX;
		};

		\nextgroupplot[semithick, xlabel = {
		\pgfplotstablegetelem{3}{timeOffsets}\of\myloadedtable%
		$t_0 +
		\SI[round-mode=places,round-precision=2]{\pgfplotsretval}{\second}$
		},
		legend columns=4,
		legend cell align=right,
		legend style={font=\small,
		at={(0,1.05)},
		anchor=south,
		fill=none},
		]
		\draw[color=gray!20] (0.162533,-0.30) -- (0.162533,+0.30); % X
		\draw[color=gray!20] (0.362533+0.05,-0.30) -- (0.362533-0.05,+0.30); % Y
		\pgfplotsinvokeforeach{0,...,12} { % level 13 and 14 are empty
		\pgfplotscolormapaccess[0:14]{#1}{jet}
		\def\TEMP{\definecolor{mycolor#1}{rgb}}
		\expandafter\TEMP\expandafter{\pgfmathresult}
		\pgfplotstablegetelem{3}{timeSteps}\of\myloadedtable
		\edef\fnameX{\fname:%
		step:\pgfplotsretval:%
		level:\pgfplotsarrayvalueofelem#1\of\contourlevels.tex}
		\addplot[contour prepared={draw color=mycolor#1}, forget plot] table \fnameX;
		\addlegendimage{line legend,color=mycolor#1}
		\addlegendentry{%
		\SI{\pgfplotsarrayvalueofelem#1\of\contourlevels}{\micro\meter/\second}}
		};

		\nextgroupplot[semithick, xlabel = {
		\pgfplotstablegetelem{4}{timeOffsets}\of\myloadedtable%
		$t_0 +
		\SI[round-mode=places,round-precision=2]{\pgfplotsretval}{\second}$
		}]
		\draw[color=gray!20] (0.162533,-0.30) -- (0.162533,+0.30); % X
		\draw[color=gray!20] (0.362533+0.05,-0.30) -- (0.362533-0.05,+0.30); % Y
		\pgfplotsinvokeforeach{0,...,12} { % level 13 and 14 are empty
		\pgfplotscolormapaccess[0:14]{#1}{jet}
		\def\TEMP{\definecolor{mycolor#1}{rgb}}
		\expandafter\TEMP\expandafter{\pgfmathresult}
		\pgfplotstablegetelem{4}{timeSteps}\of\myloadedtable
		\edef\fnameX{\fname:%
		step:\pgfplotsretval:%
		level:\pgfplotsarrayvalueofelem#1\of\contourlevels.tex}
		\addplot[contour prepared={draw color=mycolor#1}] table \fnameX;
		};

		\nextgroupplot[semithick, xlabel = {
		\pgfplotstablegetelem{5}{timeOffsets}\of\myloadedtable%
		$t_0 +
		\SI[round-mode=places,round-precision=2]{\pgfplotsretval}{\second}$
		}]
		\draw[color=gray!20] (0.162533,-0.30) -- (0.162533,+0.30); % X
		\draw[color=gray!20] (0.362533+0.05,-0.30) -- (0.362533-0.05,+0.30); % Y
		\pgfplotsinvokeforeach{0,...,12} { % level 13 and 14 are empty
		\pgfplotscolormapaccess[0:14]{#1}{jet}
		\def\TEMP{\definecolor{mycolor#1}{rgb}}
		\expandafter\TEMP\expandafter{\pgfmathresult}
		\pgfplotstablegetelem{5}{timeSteps}\of\myloadedtable
		\edef\fnameX{\fname:%
		step:\pgfplotsretval:%
		level:\pgfplotsarrayvalueofelem#1\of\contourlevels.tex}
		\addplot[contour prepared={draw color=mycolor#1}] table \fnameX;
		};
		\end{groupplot}
		\end{tikzpicture}

data/boxplot.tikz

0 → 100644

+43 −0

Original line number	Diff line number	Diff line
		\def\simulationtag{rfpitol=100e-7}
		\begin{tikzpicture}[trim axis group left, trim axis group right]
		\begin{groupplot}[
		tick label style={font=\footnotesize},
		label style={font=\small},
		%
		group style={
		y descriptions at = edge left,
		group size=3 by 1,
		horizontal sep=0.75cm
		},
		height=4cm,
		width=4.5cm,
		%
		ytick={1,2,3},
		yticklabels={experiment, simulation}
		]
		%
		\nextgroupplot[semithick, xlabel = recurrence time, x unit = s, xmode=log,
		log ticks with fixed point, xtick={5,10,20,40}]
		\addplot[mark=+, boxplot={ box extend = 0.5 }]
		table[y index=0] {generated/boxplot-data:lab:recurrence.tex};
		\addplot[mark=+, boxplot={ box extend = 0.5 }]
		table[y index=0] {generated/boxplot-data:simulation:\simulationtag:recurrence.tex};
		%
		\nextgroupplot[semithick, xlabel = rupture width, x unit=m, xmin=0, xmax=0.4,
		extra x ticks = 0.2,
		extra x tick labels = ,
		extra x tick style = { grid = major }
		]
		\addplot[mark=+, boxplot={ box extend = 0.5 }]
		table[y index=0] {generated/boxplot-data:lab:ruptureWidth.tex};
		\addplot[mark=+, boxplot={ box extend = 0.5 }]
		table[y index=0] {generated/boxplot-data:simulation:\simulationtag:ruptureWidth.tex};
		%
		\nextgroupplot[semithick, xlabel = peak slip, x unit=mm, xmode=log, log ticks with fixed point,
		xtick={0.03, 0.06, 0.12}]
		\addplot[mark=+, boxplot={ box extend = 0.5 }]
		table[y index=0] {generated/boxplot-data:lab:peakSlip.tex};
		\addplot[mark=+, boxplot={ box extend = 0.5 }]
		table[y index=0] {generated/boxplot-data:simulation:\simulationtag:peakSlip.tex};
		\end{groupplot}
		\end{tikzpicture}

data/includes.tex

0 → 100644

+51 −0

Original line number	Diff line number	Diff line
		\usepackage{pgfplots}
		\pgfplotsset{compat=1.11} % FIXME: 1.12 would be nice; debian:8.7 only has 1.11
		\usepackage{pgfplotstable}
		\usepgfplotslibrary{groupplots}
		\usepgfplotslibrary{statistics}
		\usepgfplotslibrary{units}


		%% Typeset the mu from '[xy] SI prefix=micro' as an upright mu
		%% From https://tex.stackexchange.com/a/224574
		\pgfplotsset{
		x SI prefix/micro/.style={/pgfplots/axis base prefix={axis x base 6 prefix \micro}},
		y SI prefix/micro/.style={/pgfplots/axis base prefix={axis y base 6 prefix \micro}},
		z SI prefix/micro/.style={/pgfplots/axis base prefix={axis z base 6 prefix \micro}},
		unit code/.code 2 args={\si{#1#2}},
		}

		%% Add support for 'groupplot [xy]label'
		%% From http://tex.stackexchange.com/a/117935/16940, see also
		%% https://sourceforge.net/p/pgfplots/feature-requests/48/
		\makeatletter
		\pgfplotsset{
		groupplot xlabel/.initial={},
		every groupplot x label/.style={
		at={($({\pgfplots@group@name\space c1r\pgfplots@group@rows.west}\|-{\pgfplots@group@name\space c1r\pgfplots@group@rows.outer south})!0.5!({\pgfplots@group@name\space c\pgfplots@group@columns r\pgfplots@group@rows.east}\|-{\pgfplots@group@name\space c\pgfplots@group@columns r\pgfplots@group@rows.outer south})$)},
		anchor=north,
		},
		groupplot ylabel/.initial={},
		every groupplot y label/.style={
		rotate=90,
		at={($({\pgfplots@group@name\space c1r1.north}-\|{\pgfplots@group@name\space c1r1.outer
		west})!0.5!({\pgfplots@group@name\space c1r\pgfplots@group@rows.south}-\|{\pgfplots@group@name\space c1r\pgfplots@group@rows.outer west})$)},
		anchor=south
		},
		execute at end groupplot/.code={%
		\node [/pgfplots/every groupplot x label]
		{\pgfkeysvalueof{/pgfplots/groupplot xlabel}};
		\node [/pgfplots/every groupplot y label]
		{\pgfkeysvalueof{/pgfplots/groupplot ylabel}};
		}
		}

		\def\endpgfplots@environment@groupplot{%
		\endpgfplots@environment@opt%
		\pgfkeys{/pgfplots/execute at end groupplot}%
		\endgroup%
		}
		\makeatother

		%% Have \includegraphics{} support tikz files
		\usepackage{tikzscale}
		No newline at end of file

data/standalone/2d-dip-contours-performance.tex

0 → 100644

+7 −0

Original line number	Diff line number	Diff line
		\documentclass[border={3cm 0cm}]{standalone}
		\usepackage{siunitx}
		\input{includes}

		\begin{document}
		\includegraphics{2d-dip-contours-performance.tikz}
		\end{document}

data/standalone/2d-dip-single-points.tex

0 → 100644

+7 −0

Original line number	Diff line number	Diff line
		\documentclass[border={3cm 0cm}]{standalone}
		\usepackage{siunitx}
		\input{includes}

		\begin{document}
		\includegraphics{2d-dip-single-points.tikz}
		\end{document}

data/standalone/2d-effort-over-tolerance.tex

0 → 100644

+7 −0

Original line number	Diff line number	Diff line
		\documentclass[border={3cm 0cm}]{standalone}
		\usepackage{siunitx}
		\input{includes}

		\begin{document}
		\includegraphics{2d-effort-over-tolerance.tikz}
		\end{document}

data/standalone/2d-velocity-contours-threequakes.tex

0 → 100644

+7 −0

Original line number	Diff line number	Diff line
		\documentclass[border={3cm 0cm}]{standalone}
		\usepackage{siunitx}
		\input{includes}

		\begin{document}
		\includegraphics{2d-velocity-contours-threequakes}
		\end{document}

data/standalone/2d-velocity-contours-zoom.tex

0 → 100644

+7 −0

Original line number	Diff line number	Diff line
		\documentclass[border={3cm 0cm}]{standalone}
		\usepackage{siunitx}
		\input{includes}

		\begin{document}
		\includegraphics{2d-velocity-contours-zoom}
		\end{document}

data/standalone/3d-performance.tex

0 → 100644

+7 −0

Original line number	Diff line number	Diff line
		\documentclass[border={3cm 0cm}]{standalone}
		\usepackage{siunitx}
		\input{includes}

		\begin{document}
		\includegraphics{3d-performance.tikz}
		\end{document}

data/standalone/3d-velocity-contours.tex

0 → 100644

+7 −0

Original line number	Diff line number	Diff line
		\documentclass[border={3cm 0cm}]{standalone}
		\usepackage{siunitx}
		\input{includes}

		\begin{document}
		\includegraphics{3d-velocity-contours.tikz}
		\end{document}

data/standalone/Makefile

0 → 100644

+18 −0

Original line number	Diff line number	Diff line
		PDFs=\
		2d-dip-contours-performance.pdf \
		2d-dip-single-points.pdf \
		2d-effort-over-tolerance.pdf \
		2d-velocity-contours-threequakes.pdf \
		2d-velocity-contours-zoom.pdf \
		3d-performance.pdf \
		3d-velocity-contours.pdf \
		boxplot.pdf
		PNGs=$(PDFs:.pdf=.png)

		pdf: $(PDFs)
		%.pdf: standalone/%.tex
		latexmk -pdf $<

		png: $(PNGs)
		%.png: %.pdf
		convert -density 300 $< -quality 100 $@

data/standalone/boxplot.tex

0 → 100644

+7 −0

Original line number	Diff line number	Diff line
		\documentclass[border={3cm 0cm}]{standalone}
		\usepackage{siunitx}
		\input{includes}

		\begin{document}
		\input{boxplot.tikz}
		\end{document}

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

0 → 100644

+108 −0

Original line number	Diff line number	Diff line
		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')))
		}

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

0 → 100644

+68 −0

Original line number	Diff line number	Diff line
		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

0 → 100644

+41 −0

Original line number	Diff line number	Diff line
		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

0 → 100644

+58 −0

Original line number	Diff line number	Diff line
		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

0 → 100644

+76 −0

Original line number	Diff line number	Diff line
		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

0 → 100644

+95 −0

Original line number	Diff line number	Diff line
		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

0 → 100644

+112 −0

Original line number	Diff line number	Diff line
		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~

0 → 100644

+112 −0

Original line number	Diff line number	Diff line
		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

0 → 100644

+62 −0

Original line number	Diff line number	Diff line
		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

0 → 100644

+67 −0

Original line number	Diff line number	Diff line
		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

0 → 100644

+42 −0

Original line number	Diff line number	Diff line
		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

0 → 100644

+9 −0

Original line number	Diff line number	Diff line
		#/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/tresca-0.6/
		experiment = ~/group/publications/2016-RosenauCorbiDominguezRudolfRitterPipping
		output = generated

data/tools/debug.py

0 → 100644

+135 −0

Original line number	Diff line number	Diff line
		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/init.py

0 → 100644

+0 −0

Original line number	Diff line number	Diff line

data/tools/debug/diffplot.py

0 → 100644

+17 −0

Original line number	Diff line number	Diff line
		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

0 → 100644

+8 −0

Original line number	Diff line number	Diff line
		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

0 → 100644

+6 −0

Original line number	Diff line number	Diff line
		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

0 → 100644

+25 −0

Original line number	Diff line number	Diff line
		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

0 → 100644

+73 −0

Original line number	Diff line number	Diff line
		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

0 → 100644

+25 −0

Original line number	Diff line number	Diff line
		#!/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

0 → 100644

+16 −0

Original line number	Diff line number	Diff line
		#!/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

0 → 100644

+13 −0

Original line number	Diff line number	Diff line
		#!/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

0 → 100644

+20 −0

Original line number	Diff line number	Diff line
		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/main.py

0 → 100644

+71 −0

Original line number	Diff line number	Diff line
		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.io import read_h5file

		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]


		NBODIES = 2
		FINAL_TIME = 150 # s
		FINAL_VELOCITY = 2e-4 # 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
		}

		h5file = read_h5file()
		print(list(h5file.keys()))

		iterations(h5file, FINAL_TIME)

		for body_ID in range(NBODIES):
		body = 'body' + str(body_ID)

		if body not in h5file:
		continue

		coords = np.array(h5file[body + '/coordinates'])
		patch = build_patch(coords, 0.05)

		friction_stats(h5file, body_ID, FINAL_TIME, [], [0, FINAL_TIME])

		plt.show()

		h5file.close()
		No newline at end of file

data/tools/support/init.py

0 → 100644

+1 −0

Original line number	Diff line number	Diff line

data/tools/support/init.pyc

0 → 100644

+106 B

File added.

No diff preview for this file type.

View file

data/tools/support/find_quakes.py

0 → 100644

+16 −0

Original line number	Diff line number	Diff line
		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.pyc

0 → 100644

+741 B

File added.

No diff preview for this file type.

View file

data/tools/support/find_quakes.py~

0 → 100644

+17 −0

Original line number	Diff line number	Diff line
		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]

data/tools/support/friction_stats.py

0 → 100644

+74 −0

Original line number	Diff line number	Diff line
		import numpy as np
		import matplotlib.pyplot as plt

		def friction_stats(h5file, body_ID, FINAL_TIME, patch=[], interval=[], TANGENTIAL_COORDS=1):
		body = 'body' + str(body_ID) # 'frictionalBoundary' 'body' + str(body_ID)

		coords = np.array(h5file[body + '/coordinates'])
		if len(patch) == 0:
		patch = np.linspace(0, len(coords)-1, len(coords), endpoint=True, dtype='int8')

		# read time
		time = np.array(h5file['relativeTime']) * FINAL_TIME
		time = np.delete(time, 0)
		if len(interval) == 0:
		interval = [0, FINAL_TIME]
		t = [i for i in range(len(time)) if time[i]>=interval[0] and time[i]<=interval[1]]
		time = time[t]

		fig = plt.figure()

		# velocity data
		v = abs(np.array(h5file[body + '/velocity']))
		v_t = v[t,:,TANGENTIAL_COORDS]
		v_tx = v_t[:,patch]
		# statistics
		avg_v = np.average(v_tx, axis=1)
		min_v = np.min(v_tx, axis=1)
		max_v = np.max(v_tx, axis=1)
		# plot
		ax_slip = fig.add_subplot(3, 1, 1)
		#ax_slip.plot(time, min_v, color='gray', linestyle='--')
		ax_slip.plot(time, avg_v, color='black', linestyle='-')
		#ax_slip.plot(time, max_v, color='gray', linestyle='--')
		ax_slip.set_ylabel('slip rate')
		#ax_slip.set_yscale('log')
		#-------------------------

		print(np.min(min_v))
		print(np.average(avg_v))
		print(np.max(max_v))

		# state
		states = np.array(h5file[body + '/state'])
		states_t = states[t,:]
		states_tx = states_t[:,patch]
		# statistics
		avg_states = np.average(states_tx, axis=1)
		min_states = np.min(states_tx, axis=1)
		max_states = np.max(states_tx, axis=1)
		# plot
		ax_state = fig.add_subplot(3, 1, 2)
		ax_state.plot(time, min_states, color='gray', linestyle='--')
		ax_state.plot(time, avg_states, color='black', linestyle='-')
		ax_state.plot(time, max_states, color='gray', linestyle='--')
		ax_state.set_ylabel('state')
		#-------------------------

		# friction coefficient
		friction_coeff = np.array(h5file[body + '/coefficient'])
		friction_coeff_t = friction_coeff[t,:]
		friction_coeff_tx = friction_coeff_t[:,patch]
		# statistics
		avg_friction_coeff = np.average(friction_coeff_tx, axis=1)
		min_friction_coeff = np.min(friction_coeff_tx, axis=1)
		max_friction_coeff = np.max(friction_coeff_tx, axis=1)
		# plot
		ax_friction = fig.add_subplot(3, 1, 3)
		ax_friction.plot(time, min_friction_coeff, color='gray', linestyle='--')
		ax_friction.plot(time, avg_friction_coeff, color='black', linestyle='-')
		ax_friction.plot(time, max_friction_coeff, color='gray', linestyle='--')
		ax_friction.set_ylabel('friction coefficient')
		#-------------------------

		fig.canvas.draw()
		No newline at end of file

data/tools/support/io.py

0 → 100644

+49 −0

Original line number	Diff line number	Diff line
		import configparser as cp
		import os
		import h5py

		def remove_comment(str, char = "#"):
		split_str = str.split(char, 1)
		return split_str[0]

		# tag = 'simulation', 'experiment', 'output'
		def read_h5file(tag = 'simulation'):
		# read config ini
		config = cp.ConfigParser()
		config_path = os.path.join('tools/config.ini')
		config.read(config_path)
		path = config.get('directories', tag)

		# read hdf5 output file
		h5path = os.path.join(path)

		return h5py.File(os.path.join(h5path, 'output.h5'), 'r')

		# tag = 'simulation', 'experiment', 'output'
		def read_params(file_name, tag = 'simulation'):
		# read config ini
		config = cp.ConfigParser()
		config_path = os.path.join('tools/config.ini')
		config.read(config_path)
		path = config.get('directories', tag)

		# read cfg parameter file
		params = cp.ConfigParser()
		params_path = os.path.join(path, file_name)
		params.read(params_path)

		#with open(params_path) as stream:
		# params.read_string("[top]\n" + stream.read())

		out = {
		'L' : float(remove_comment(params.get('boundary.friction', 'L'))),
		'V0' : float(remove_comment(params.get('boundary.friction', 'V0'))),
		'mu0': float(remove_comment(params.get('boundary.friction', 'mu0'))),
		'a' : float(remove_comment(params.get('boundary.friction.weakening', 'a'))),
		'b' : float(remove_comment(params.get('boundary.friction.weakening', 'b'))),
		'bodyCount' : int(remove_comment(params.get('problem', 'bodyCount'))),
		'finalTime' : float(remove_comment(params.get('problem', 'finalTime'))),
		'finalVelocity' : float(remove_comment(params.get('boundary.dirichlet', 'finalVelocity')))
		}

		return out
		No newline at end of file

data/tools/support/iterations.py

0 → 100644

+45 −0

Original line number	Diff line number	Diff line
		import numpy as np
		import matplotlib.pyplot as plt

		def iterations(h5file, FINAL_TIME, interval = []):

		# read time
		time = np.array(h5file['relativeTime']) * FINAL_TIME
		time = np.delete(time, 0)
		if len(interval) == 0:
		interval = [0, FINAL_TIME]
		t = [i for i in range(len(time)) if time[i]>=interval[0] and time[i]<=interval[1]]
		time = time[t]

		tau = np.array(h5file['relativeTimeIncrement']) * FINAL_TIME
		tau = np.delete(tau, 0)

		# read fpi iterations
		fpi_final = np.array(h5file['iterations/fixedPoint/final'])
		fpi_final = np.delete(fpi_final, 0)
		#fpi_total = np.array(h5file['iterations/fixedPoint/total'])

		# read multigrid iterations
		multigrid_final = np.array(h5file['iterations/multiGrid/final'])
		multigrid_final = np.delete(multigrid_final, 0)
		#multigrid_total = np.array(h5file['iterations/multiGrid/total'])

		# plot
		fig = plt.figure()

		ax_fpi = fig.add_subplot(3, 1, 1)
		ax_fpi.plot(time, fpi_final[t], color='black', linestyle='-')
		ax_fpi.set_ylabel('fpi')
		#-------------------------

		ax_mg = fig.add_subplot(3, 1, 2)
		ax_mg.plot(time, multigrid_final[t], color='black', linestyle='-')
		ax_mg.set_ylabel('multigrid')
		#-------------------------

		ax_tau = fig.add_subplot(3, 1, 3)
		ax_tau.plot(time, tau[t], color='black', linestyle='-')
		ax_tau.set_ylabel('tau')
		#-------------------------

		fig.canvas.draw()
		No newline at end of file

data/tools/support/max_distance.py

0 → 100644

+17 −0

Original line number	Diff line number	Diff line
		from itertools import combinations


		def square_distance(x, y):
		return sum([(xi-yi)**2 for xi, yi in zip(x, y)])


		def max_distance(points):
		max_square_distance = 0
		max_pair = (0, 0)
		for pair in combinations(points, 2):
		sq_dist = square_distance(*pair)

		if sq_dist > max_square_distance:
		max_square_distance = sq_dist
		max_pair = pair
		return max_pair

data/tools/support/max_distance.pyc

0 → 100644

+884 B

File added.

No diff preview for this file type.

View file

data/tools/support/max_distance.py~

0 → 100644

+17 −0

Original line number	Diff line number	Diff line
		from itertools import combinations


		def square_distance(x, y):
		return sum([(xi-yi)**2 for xi, yi in zip(x, y)])


		def max_distance(points):
		max_square_distance = 0
		max_pair = (0, 0)
		for pair in combinations(points, 2):
		sq_dist = square_distance(*pair)
		print(sq_dist)
		if sq_dist > max_square_distance:
		max_square_distance = sq_dist
		max_pair = pair
		return max_pair

data/tools/support/max_velocity.py

0 → 100644

+4 −0

Original line number	Diff line number	Diff line
		import numpy as np

		def max_velocity(x) :
		return maximum(norm(x))
		No newline at end of file

data/tools/support/max_velocity.pyc

0 → 100644

+649 B

File added.

No diff preview for this file type.

View file

data/tools/support/max_velocity.py~

0 → 100644

+11 −0

Original line number	Diff line number	Diff line
		import numpy as np

		def max_velocity(x) :
		size = x.shape
		res = np.zeros(size[0])

		for time_step in range(size[0]) :
		for vertex in range(size[1]) :
		res[time_step] = max(res[time_step], np.linalg.norm(x[time_step, vertex, :]))

		return res
		No newline at end of file

data/tools/support/maximum.py

0 → 100644

+11 −0

Original line number	Diff line number	Diff line
		import numpy as np


		def maximum(x):
		size = x.shape
		res = np.zeros(size[0])

		for time_step in range(size[0]):
		res[time_step] = max(res[time_step], max(x[time_step, :]))

		return res

data/tools/support/maximum.pyc

0 → 100644

+560 B

File added.

No diff preview for this file type.

View file

data/tools/support/maximum.py~

0 → 100644

+11 −0

Original line number	Diff line number	Diff line
		import numpy as np

		def maximum(x) :
		size = x.shape
		res = np.zeros(size[0])

		for time_step in range(size[0]) :
		for vertex in range(size[1]) :
		res[time_step] = max(res[time_step], max(x[time_step, vertex]))

		return res
		No newline at end of file

data/tools/support/norm.py

0 → 100644

+12 −0

Original line number	Diff line number	Diff line
		import numpy as np


		def norm(x):
		size = x.shape
		res = np.zeros((size[0], size[1]))

		for time_step in range(size[0]):
		for vertex in range(size[1]):
		res[time_step, vertex] = np.linalg.norm(x[time_step, vertex, :])

		return res

data/tools/support/norm.pyc

0 → 100644

+616 B

File added.

No diff preview for this file type.

View file

data/tools/support/norm.py~

0 → 100644

+11 −0

Original line number	Diff line number	Diff line
		import numpy as np

		def norm(x) :
		size = x.shape
		res = np.zeros(size[0], size[1])

		for time_step in range(size[0]) :
		for vertex in range(size[1]) :
		res[time_step, vertex] = np.linalg.norm(x[time_step, vertex, :])

		return res

data/tools/support/slip_beginnings.py

0 → 100644

+15 −0

Original line number	Diff line number	Diff line
		import array as ar


		def slip_beginnings(x):
		# returns indicies i for which x[i-1]=0 and x[i]>0
		starts = ar.array('i', (0 for i in range(len(x))))
		prev = False

		length = 0
		for i in range(len(x)):
		if (not prev and x[i]):
		starts[length] = i
		length += 1
		prev = x[i]
		return starts[:length]

data/tools/support/slip_beginnings.pyc

0 → 100644

+837 B

File added.

No diff preview for this file type.

View file

data/tools/support/slip_beginnings.py~

0 → 100644

+11 −0

Original line number	Diff line number	Diff line
		import numpy as np

		def slip_beginnings(x) :
		# returns indicies i for which x[i-1]=0 and x[i]>0
		starts = []
		prev = False
		for i in range(len(x)) :
		if (not prev and x[i]) :
		starts.append(i)
		prev = x[i]
		return starts

data/tools/support/slip_endings.py

0 → 100644

+16 −0

Original line number	Diff line number	Diff line
		import array as ar


		def slip_endings(x):
		# returns indicies i for which x[i-1]>0 and x[i]=0
		ends = ar.array('i', (0 for i in range(len(x))))
		prev = False

		length = 0
		for i in range(len(x)):
		if (prev and not x[i]):
		ends[length] = i
		length += 1
		prev = x[i]

		return ends[:length]

data/tools/support/slip_endings.pyc

0 → 100644

+823 B

File added.

No diff preview for this file type.

View file

data/tools/support/slip_endings.py~

0 → 100644

+11 −0

Original line number	Diff line number	Diff line
		import numpy as np

		def slip_endings(x) :
		# returns indicies i for which x[i-1]>0 and x[i]=0
		starts = np.array([])
		prev = False
		for i in range(len(x)) :
		if (prev and not x[i]) :
		starts.append(i)
		prev = x[i]
		return starts
		No newline at end of file

data/tools/support/trapezoidal.py

0 → 100644

+8 −0

Original line number	Diff line number	Diff line
		# x: vector of domain values
		# y: vector of function values f(x)
		def trapezoidal(x, y) :
		# returns integral of (x, y = f(x)) calculated with trapezoidal rule
		res = 0.0
		for i in range(1, len(x)) :
		res += (x[i] - x[i - 1]) * (y[i] + y[i - 1]) / 2.0
		return res
		No newline at end of file

data/tools/support/trapezoidal.py~

0 → 100644

+8 −0

Original line number	Diff line number	Diff line
		# x: vector of domain values
		# y: vector of function values f(x)
		def trapezoidal(x, y) :
		# returns integral of (x, y = f(x)) calculated with trapezoidal rule
		res = 0.0
		for i in range(1, len(x)) :
		res += (x[i] - x[i - 1]) * (y[i] + y[i - 1]) / 2.0
		return res
		No newline at end of file

Source

Target

Files

Some changes are not shown.