Use overlapping schwarz for oscillator tutorial

oscillator-overlap/README.md

@@ -0,0 +1,42 @@
+---
+title: Oscillator (overlapping domain decomposition)
+permalink: tutorials-oscillator-overlap.html
+keywords: Python, ODE
+summary: We solve an oscillator with two masses in a partitioned fashion with overlapping domain decomposition. Each mass is solved by an independent ODE.
+---
+
+{% note %}
+Get the [case files of this tutorial](https://github.com/precice/tutorials/tree/master/oscillator-overlap). Read how in the [tutorials introduction](https://www.precice.org/tutorials.html).
+{% endnote %}
+
+## Setup
+
+This tutorial solves the same problem as the [oscillator tutorial](https://precice.org/tutorials-oscillator.html), but applies a different domain decomposition strategy. See the oscillator tutorial for details on the general setup. The partitioning of the mass-spring system is shown here:
+
+![Schematic drawing of oscillator example with overlapping domain decomposition](images/tutorials-oscillator-overlap-dd.png)
+
+Note that this case applies an overlapping Schwarz-type coupling method and not (like most other tutorials in this repository) a Dirichlet-Neumann coupling. This results in a symmetric setup of the solvers. We will refer to the solver computing the trajectory of $m_1$ as `Mass-Left` and to the solver computing the trajectory of $m_2$ as `Mass-Right`. For more information, please refer to [1].
+
+## Available solvers
+
+This tutorial is only available in Python. You need to have preCICE and the Python bindings installed on your system.
+
+- *Python*: An example solver using the preCICE [Python bindings](https://www.precice.org/installation-bindings-python.html). This solver also depends on the Python libraries `numpy`, which you can get from your system package manager or with `pip3 install --user <package>`.
+
+## Running the Simulation
+
+### Python
+
+Open two separate terminals and start both participants by calling:
+
+```bash
+cd python
+./run.sh -l
+```
+
+and
+
+```bash
+cd python
+./run.sh -r
+```

oscillator-overlap/clean-tutorial.sh

@@ -0,0 +1 @@
+../tools/clean-tutorial-base.sh

oscillator-overlap/images/tutorials-oscillator-overlap-dd.tex

@@ -0,0 +1,37 @@
+\documentclass{standalone}
+\usepackage{tikz}
+\usetikzlibrary{positioning,calc,patterns,decorations.pathmorphing}
+
+\begin{document}
+
+\begin{tikzpicture}
+    \tikzstyle{spring}=[thick,decorate,decoration={zigzag,pre length=0.3cm,post length=0.3cm,segment length=6}]
+    \tikzstyle{ground}=[fill,pattern=north east lines,draw=none,minimum width=0.75cm,minimum height=0.3cm]
+    \tikzstyle{mass}=[draw, fill=gray!25, thick, minimum width=1cm, minimum height=1cm]
+    \node (W1) [ground, rotate=-90, minimum width=3cm] {};
+    \node (M1) [mass, right= 2cm of W1.north] {$m_1$};
+    \node (M2fake) [circle, draw=black, fill=gray!25, right= 2cm of M1] {};
+
+    \draw [dashed] (M2fake.north) -- ++ (0,.5);
+    \draw [->] ($(M2fake.north) + (0,.25)$) -- node[above](U2ToM1){$u_2$} ++ (1,0);
+    \draw [dashed] (M1.north) -- ++ (0,.5);
+    \draw [->] ($(M1.north) + (0,.25)$) -- node[above](U1FromM1){$u_1$} ++ (1,0);
+    \draw [spring] (W1.north) -- node[above]{$k_1$}(M1.west);
+    \draw [spring] (M1.east) -- node[above]{$k_{12}$}(M2fake.west);
+
+    \node (M1fake) [circle, draw=black, fill=gray!25, right= 1cm of M2fake] {};
+    \node (M2) [mass, right= 2cm of M1fake] {$m_2$};
+    \draw [dashed] (M2.north) -- ++ (0,.5);
+    \draw [->] ($(M2.north) + (0,.25)$) -- node[above](U2FromM2){$u_2$} ++ (1,0);
+    \draw [dashed] (M1fake.north) -- ++ (0,.5);
+    \draw [->] ($(M1fake.north) + (0,.25)$) -- node[above](U1ToM2){$u_1$} ++ (1,0);
+
+    \node (W2) [ground, rotate=90, minimum width=3cm, right= 2cm of M2, xshift=-1.5cm] {};
+    \draw [spring] (M2.east) -- node[above]{$k_2$}(W2.north);
+    \draw [spring] (M1fake.east) -- node[below]{$k_{12}$}(M2.west);
+
+    \draw[](U2FromM2.north) edge[->,out=135, in=45] node[above right]{ghost layer for $m_1$} (U2ToM1.north);
+    \draw[](U1FromM1.north) edge[->,out=45, in=135] node[above left]{ghost layer for $m_2$} (U1ToM2.north);
+
+\end{tikzpicture}
+\end{document}

oscillator-overlap/mass-left-python/clean.sh

@@ -0,0 +1,8 @@
+#!/bin/sh
+set -e -u
+
+. ../../tools/cleaning-tools.sh
+
+rm -rfv ./output/
+
+clean_precice_logs .

oscillator-overlap/mass-left-python/run.sh

@@ -0,0 +1,4 @@
+#!/bin/sh
+set -e -u
+
+python3 ../solver-python/oscillator.py Mass-Left

oscillator-overlap/mass-right-python/clean.sh

@@ -0,0 +1,8 @@
+#!/bin/sh
+set -e -u
+
+. ../../tools/cleaning-tools.sh
+
+rm -rfv ./output/
+
+clean_precice_logs .

oscillator-overlap/mass-right-python/run.sh

@@ -0,0 +1,4 @@
+#!/bin/sh
+set -e -u
+
+python3 ../solver-python/oscillator.py Mass-Right

oscillator-overlap/precice-config.xml

@@ -0,0 +1,47 @@
+<?xml version="1.0" encoding="UTF-8" ?>
+<precice-configuration>
+  <log enabled="1">
+    <sink filter="%Severity% > debug" />
+  </log>
+
+  <data:scalar name="Displacement-Left" />
+  <data:scalar name="Displacement-Right" />
+
+  <mesh name="Mass-Left-Mesh" dimensions="2">
+    <use-data name="Displacement-Left" />
+    <use-data name="Displacement-Right" />
+  </mesh>
+
+  <mesh name="Mass-Right-Mesh" dimensions="2">
+    <use-data name="Displacement-Left" />
+    <use-data name="Displacement-Right" />
+  </mesh>
+
+  <participant name="Mass-Left">
+    <provide-mesh name="Mass-Left-Mesh" />
+    <write-data name="Displacement-Left" mesh="Mass-Left-Mesh" />
+    <read-data  name="Displacement-Right" mesh="Mass-Left-Mesh" />
+  </participant>
+
+  <participant name="Mass-Right">
+    <receive-mesh name="Mass-Left-Mesh" from="Mass-Left" />
+    <provide-mesh name="Mass-Right-Mesh" />
+    <write-data name="Displacement-Right" mesh="Mass-Right-Mesh" />
+    <read-data  name="Displacement-Left" mesh="Mass-Right-Mesh" />
+    <mapping:nearest-neighbor direction="write" from="Mass-Right-Mesh" to="Mass-Left-Mesh" constraint="consistent" />
+    <mapping:nearest-neighbor direction="read"  from="Mass-Left-Mesh" to="Mass-Right-Mesh" constraint="consistent" />
+  </participant>
+
+  <m2n:sockets acceptor="Mass-Left" connector="Mass-Right" exchange-directory=".." />
+
+  <coupling-scheme:serial-implicit>
+    <participants first="Mass-Left" second="Mass-Right" />
+    <max-time value="1" />
+    <time-window-size value="0.01" />
+    <max-iterations value="200" />
+    <relative-convergence-measure data="Displacement-Left" mesh="Mass-Left-Mesh" limit="1e-10"/>
+    <relative-convergence-measure data="Displacement-Right" mesh="Mass-Left-Mesh" limit="1e-10"/>
+    <exchange data="Displacement-Left" mesh="Mass-Left-Mesh" from="Mass-Left" to="Mass-Right" initialize="true"/>
+    <exchange data="Displacement-Right" mesh="Mass-Left-Mesh" from="Mass-Right" to="Mass-Left" initialize="true"/>
+  </coupling-scheme:serial-implicit>
+</precice-configuration>

oscillator-overlap/solver-python/oscillator.py

@@ -0,0 +1,190 @@
+from __future__ import division
+
+import argparse
+import numpy as np
+from numpy.linalg import eig
+import precice
+from enum import Enum
+import csv
+import os
+
+import problemDefinition
+
+
+class Scheme(Enum):
+    NEWMARK_BETA = "Newmark_beta"
+    GENERALIZED_ALPHA = "generalized_alpha"
+
+
+class Participant(Enum):
+    MASS_LEFT = "Mass-Left"
+    MASS_RIGHT = "Mass-Right"
+
+
+parser = argparse.ArgumentParser()
+parser.add_argument("participantName", help="Name of the solver.", type=str, choices=[p.value for p in Participant])
+parser.add_argument("-ts", "--time-stepping", help="Time stepping scheme being used.", type=str,
+                    choices=[s.value for s in Scheme], default=Scheme.NEWMARK_BETA.value)
+args = parser.parse_args()
+
+participant_name = args.participantName
+
+if participant_name == Participant.MASS_LEFT.value:
+    write_data_name = 'Displacement-Left'
+    read_data_name = 'Displacement-Right'
+    mesh_name = 'Mass-Left-Mesh'
+
+    this_mass = problemDefinition.MassLeft
+    this_spring = problemDefinition.SpringLeft
+    connecting_spring = problemDefinition.SpringMiddle
+    other_mass = problemDefinition.MassRight
+
+elif participant_name == Participant.MASS_RIGHT.value:
+    read_data_name = 'Displacement-Left'
+    write_data_name = 'Displacement-Right'
+    mesh_name = 'Mass-Right-Mesh'
+
+    this_mass = problemDefinition.MassRight
+    this_spring = problemDefinition.SpringRight
+    connecting_spring = problemDefinition.SpringMiddle
+    other_mass = problemDefinition.MassLeft
+
+else:
+    raise Exception(f"wrong participant name: {participant_name}")
+
+mass = this_mass.m
+stiffness = this_spring.k + connecting_spring.k
+u0, v0, f0, d_dt_f0 = this_mass.u0, this_mass.v0, connecting_spring.k * other_mass.u0, connecting_spring.k * other_mass.v0
+
+num_vertices = 1  # Number of vertices
+
+solver_process_index = 0
+solver_process_size = 1
+
+configuration_file_name = "../precice-config.xml"
+
+participant = precice.Participant(participant_name, configuration_file_name, solver_process_index, solver_process_size)
+
+dimensions = participant.get_mesh_dimensions(mesh_name)
+
+vertex = np.zeros(dimensions)
+read_data = np.zeros(num_vertices)
+write_data = connecting_spring.k * u0 * np.ones(num_vertices)
+
+vertex_ids = [participant.set_mesh_vertex(mesh_name, vertex)]
+
+if participant.requires_initial_data():
+    participant.write_data(mesh_name, write_data_name, vertex_ids, write_data)
+
+participant.initialize()
+precice_dt = participant.get_max_time_step_size()
+my_dt = precice_dt  # use my_dt < precice_dt for subcycling
+
+# Initial Conditions
+a0 = (f0 - stiffness * u0) / mass
+u = u0
+v = v0
+a = a0
+t = 0
+
+# Generalized Alpha Parameters
+if args.time_stepping == Scheme.GENERALIZED_ALPHA.value:
+    alpha_f = 0.4
+    alpha_m = 0.2
+elif args.time_stepping == Scheme.NEWMARK_BETA.value:
+    alpha_f = 0.0
+    alpha_m = 0.0
+m = 3 * [None]  # will be computed for each timestep depending on dt
+gamma = 0.5 - alpha_m + alpha_f
+beta = 0.25 * (gamma + 0.5)
+
+positions = []
+velocities = []
+times = []
+
+u_write = [u]
+v_write = [v]
+t_write = [t]
+
+while participant.is_coupling_ongoing():
+    if participant.requires_writing_checkpoint():
+        u_cp = u
+        v_cp = v
+        a_cp = a
+        t_cp = t
+        # store data for plotting and postprocessing
+        positions += u_write
+        velocities += v_write
+        times += t_write
+
+    # compute time step size for this time step
+    precice_dt = participant.get_max_time_step_size()
+    dt = np.min([precice_dt, my_dt])
+    read_time = (1 - alpha_f) * dt
+    read_data = participant.read_data(mesh_name, read_data_name, vertex_ids, read_time)
+    f = connecting_spring.k * read_data[0]
+
+    # do generalized alpha step
+    m[0] = (1 - alpha_m) / (beta * dt**2)
+    m[1] = (1 - alpha_m) / (beta * dt)
+    m[2] = (1 - alpha_m - 2 * beta) / (2 * beta)
+    k_bar = stiffness * (1 - alpha_f) + m[0] * mass
+    u_new = (f - alpha_f * stiffness * u + mass * (m[0] * u + m[1] * v + m[2] * a)) / k_bar
+    a_new = 1.0 / (beta * dt**2) * (u_new - u - dt * v) - (1 - 2 * beta) / (2 * beta) * a
+    v_new = v + dt * ((1 - gamma) * a + gamma * a_new)
+    t_new = t + dt
+
+    write_data = [u_new]
+
+    participant.write_data(mesh_name, write_data_name, vertex_ids, write_data)
+
+    participant.advance(dt)
+
+    if participant.requires_reading_checkpoint():
+        u = u_cp
+        v = v_cp
+        a = a_cp
+        t = t_cp
+        # empty buffers for next window
+        u_write = []
+        v_write = []
+        t_write = []
+
+    else:
+        u = u_new
+        v = v_new
+        a = a_new
+        t = t_new
+
+        # write data to buffers
+        u_write.append(u)
+        v_write.append(v)
+        t_write.append(t)
+
+# store final result
+u = u_new
+v = v_new
+a = a_new
+u_write.append(u)
+v_write.append(v)
+t_write.append(t)
+positions += u_write
+velocities += v_write
+times += t_write
+
+participant.finalize()
+
+# print errors
+error = np.max(abs(this_mass.u_analytical(np.array(times)) - np.array(positions)))
+print("Error w.r.t analytical solution:")
+print(f"{my_dt},{error}")
+
+# output trajectory
+if not os.path.exists("output"):
+    os.makedirs("output")
+
+with open(f'output/trajectory-{participant_name}.csv', 'w') as file:
+    csv_write = csv.writer(file, delimiter=';')
+    csv_write.writerow(['time', 'position', 'velocity'])
+    for t, u, v in zip(times, positions, velocities):
+        csv_write.writerow([t, u, v])