From feb899c8c3d618292e51cebd0716b9bf1eb7e28c Mon Sep 17 00:00:00 2001
From: Benjamin Rodenberg <benjamin.rodenberg@cit.tum.de>
Date: Wed, 20 Mar 2024 13:42:06 +0100
Subject: [PATCH 01/15] Add python version of resonator.

---
 resonant-circuit/README.md                   | 83 ++++++++++++++++++++
 resonant-circuit/capacitor-python/SolverI.py | 75 ++++++++++++++++++
 resonant-circuit/coil-python/SolverU.py      | 70 +++++++++++++++++
 resonant-circuit/plot-solution.sh            | 21 +++++
 resonant-circuit/precice-config.xml          | 51 ++++++++++++
 5 files changed, 300 insertions(+)
 create mode 100644 resonant-circuit/README.md
 create mode 100644 resonant-circuit/capacitor-python/SolverI.py
 create mode 100644 resonant-circuit/coil-python/SolverU.py
 create mode 100755 resonant-circuit/plot-solution.sh
 create mode 100644 resonant-circuit/precice-config.xml

diff --git a/resonant-circuit/README.md b/resonant-circuit/README.md
new file mode 100644
index 000000000..dd381a0ed
--- /dev/null
+++ b/resonant-circuit/README.md
@@ -0,0 +1,83 @@
+## Background
+
+The purpose of this tutorial is to illustrate the usage of preCICE to couple MATLAB code. Two different MATLAB solvers will be coupled to simulate a two-element LC circuit. This type of circuit consists on a very simple system with one inductor and one capacitor:
+
+![LC circuit diagram [1]](ref_images/diagram.svg)
+
+The circuit is described by the following system of ODEs:
+
+V(t) = L dI/dt
+
+I(t) = -C dV/dt
+
+where I is the current and V the voltage of the cirucit.
+
+Each of these equations is going to be solved by a different MATLAB solver. Note that as only one scalar is solved per equation, this is a 0+1 dimensional problem.
+
+## Dependencies
+
+For running this tutorial, you have to install
+
+* **preCICE**, see [preCICE wiki](https://github.com/precice/precice/wiki/Building).
+* **MATLAB**, see [mathworks.com](https://de.mathworks.com/products.get-matlab.html).
+
+After installing both preCICE and MATLAB, you need to [build the MATLAB bindings](https://github.com/gilbertolem/precice/tree/develop/src/precice/bindings/matlab#compilation).
+
+## Running
+
+There are two different versions available, one with explicit coupling and one with implicit coupling, located in the `explicit` and `implicit` folders respectively.
+
+For running the case, first get into the desired folder (e.g. `explicit/`) and open two MATLAB instances. After adding the MATLAB bindings to the MATLAB path, run the following commands:
+
+In the first MATLAB instance one can run the solver for the current:
+```
+Solver_I
+```
+
+And in the second MATLAB instance the solver for the voltage:
+```
+Solver_U
+```
+
+The preCICE configuration file is available as `precice-config.xml`, and it is called directly in the solvers.
+
+### Running from terminal
+
+If you prefer to not open the MATLAB GUIs, one can alternatively open two shells and:
+
+Run solver for current in first shell:
+```
+# cd into desired folder (explicit or implicit)
+cd explicit
+
+# Add bindings to MATLAB path
+export MATLABPATH=<PATH OF MATLAB BINDINGS>
+
+# Run matlab code without GUI
+LD_PRELOAD=/usr/lib/x86_64-linux-gnu/libstdc++.so.6 matlab -nodisplay -nosplash -nodesktop -r "Solver_I;exit;"
+
+```
+
+Run solver for voltage in second shell:
+```
+# cd into desired folder (explicit or implicit)
+cd explicit
+
+# Add bindings to MATLAB path
+export MATLABPATH=<PATH OF MATLAB BINDINGS>
+
+# Run matlab code without GUI
+LD_PRELOAD=/usr/lib/x86_64-linux-gnu/libstdc++.so.6 matlab -nodisplay -nosplash -nodesktop -r "Solver_U;exit;"
+
+```
+
+## Visualization
+
+The solver for the current also records the current and voltage through time and at the end of the simulation saves a plot with the obtained curves, as well as the analytical solution.
+
+After successfully running the coupling, one can find the curves in their respective folder as `Curves.png`.
+
+![Example of visualization of the simulation results](ref_images/Sample_Curves.png)
+
+## References
+[1] By First Harmonic - Own work, CC BY-SA 3.0, https://commons.wikimedia.org/w/index.php?curid=21991221
diff --git a/resonant-circuit/capacitor-python/SolverI.py b/resonant-circuit/capacitor-python/SolverI.py
new file mode 100644
index 000000000..596059cdc
--- /dev/null
+++ b/resonant-circuit/capacitor-python/SolverI.py
@@ -0,0 +1,75 @@
+import precice
+import numpy as np
+from scipy.integrate import solve_ivp
+
+# Initialize and configure preCICE
+participant = precice.Participant("ParticipantI", "../precice-config.xml", 0, 1)
+
+# Geometry IDs. As it is a 0-D simulation, only one vertex is necessary.
+mesh_name ="MeshI"
+
+dimensions = participant.get_mesh_dimensions(mesh_name)
+
+vertex_ids = np.array([participant.set_mesh_vertex(mesh_name, np.zeros(dimensions))])
+
+# Data IDs
+read_data_name = "U"
+write_data_name = "I"
+
+# Simulation parameters and initial condition
+C = 2                      # Capacitance of the capacitor
+L = 1                      # Inductance of the coil
+t0 = 0                     # Initial simulation time
+t_max = 10                 # End simulation time
+Io = 1                     # Initial current
+phi = 0                    # Phase of the signal
+
+w0 = 1/np.sqrt(L*C)           # Resonant frequency
+I0 = Io*np.cos(phi)           # Initial condition for I
+U0 = -w0*L*Io*np.sin(phi)     # Initial condition for U
+
+f_I = lambda t, I, U: U/L;       # Time derivative of I; ODE determining capacitor
+
+# Initialize simulation
+if participant.requires_initial_data():
+    participant.write_data(mesh_name, write_data_name, vertex_ids, np.array([I0]))
+
+participant.initialize()
+
+dt = participant.get_max_time_step_size()
+
+# Start simulation
+t = t0 + dt
+while participant.is_coupling_ongoing():
+
+    # Record checkpoint if necessary
+    if participant.requires_writing_checkpoint():
+        I0_checkpoint = I0
+        U0_checkpoint = U0
+
+    # Make Simulation Step
+    sol = solve_ivp(lambda t, y: f_I(t, y, U0), [t0, t], np.array([I0]))
+    I0 = sol.y[-1,-1]
+    print(sol.y)
+    print(I0)
+
+    # Exchange data
+    participant.write_data(mesh_name, write_data_name, vertex_ids, np.array([I0]))
+    participant.advance(dt)
+
+    # Recover checkpoint if not converged, else finish time step
+    if participant.requires_reading_checkpoint():
+        I0 = I0_checkpoint
+        U0 = U0_checkpoint
+    else:
+        dt = participant.get_max_time_step_size()
+        U0 = participant.read_data(mesh_name, read_data_name, vertex_ids, dt)[0]
+        t0 = t
+        t = t0 + dt
+
+# Stop coupling
+participant.finalize()
+
+I_analytical = lambda t: Io*np.cos(t*w0)
+error = I0 - I_ana(10)
+print(f"error: {error}")
\ No newline at end of file
diff --git a/resonant-circuit/coil-python/SolverU.py b/resonant-circuit/coil-python/SolverU.py
new file mode 100644
index 000000000..f11a0ce61
--- /dev/null
+++ b/resonant-circuit/coil-python/SolverU.py
@@ -0,0 +1,70 @@
+import precice
+import numpy as np
+from scipy.integrate import solve_ivp
+
+# Initialize and configure preCICE
+participant = precice.Participant("ParticipantU", "../precice-config.xml", 0, 1)
+
+# Geometry IDs. As it is a 0-D simulation, only one vertex is necessary.
+mesh_name ="MeshU"
+
+dimensions = participant.get_mesh_dimensions(mesh_name)
+
+vertex_ids = np.array([participant.set_mesh_vertex(mesh_name, np.zeros(dimensions))])
+
+# Data IDs
+read_data_name = "I"
+write_data_name = "U"
+
+# Simulation parameters and initial condition
+C = 2                      # Capacitance of the capacitor
+L = 1                      # Inductance of the coil
+t0 = 0                     # Initial simulation time
+t_max = 10                 # End simulation time
+Io = 1                     # Initial current
+phi = 0                    # Phase of the signal
+
+w0 = 1/np.sqrt(L*C)           # Resonant frequency
+I0 = Io*np.cos(phi)           # Initial condition for I
+U0 = -w0*L*Io*np.sin(phi)     # Initial condition for U
+
+f_U = lambda t, U, I: -I/C      # Time derivative of U
+
+# Initialize simulation
+if participant.requires_initial_data():
+    participant.write_data(mesh_name, write_data_name, vertex_ids, np.array([U0]))
+
+participant.initialize()
+
+dt = participant.get_max_time_step_size()
+
+# Start simulation
+t = t0 + dt
+while participant.is_coupling_ongoing():
+
+    # Record checkpoint if necessary
+    if participant.requires_writing_checkpoint():
+        I0_checkpoint = I0
+        U0_checkpoint = U0
+
+    # Make Simulation Step
+    sol = solve_ivp(lambda t, y: f_U(t, y, I0), [t0, t], np.array([U0]))
+    U0 = sol.y[-1,-1]
+    print(U0)
+
+    # Exchange data
+    participant.write_data(mesh_name, write_data_name, vertex_ids, np.array([U0]))
+    participant.advance(dt)
+
+    # Recover checkpoint if not converged, else finish time step
+    if participant.requires_reading_checkpoint():
+        I0 = I0_checkpoint
+        U0 = U0_checkpoint
+    else:
+        dt = participant.get_max_time_step_size()
+        I0 = participant.read_data(mesh_name, read_data_name, vertex_ids, dt)[0]
+        t0 = t
+        t = t0 + dt
+
+# Stop coupling
+participant.finalize()
diff --git a/resonant-circuit/plot-solution.sh b/resonant-circuit/plot-solution.sh
new file mode 100755
index 000000000..30cd21a21
--- /dev/null
+++ b/resonant-circuit/plot-solution.sh
@@ -0,0 +1,21 @@
+#!/bin/sh
+
+if [ "${1:-}" = "" ]; then
+    echo "No target directory specified. Please specify the directory of the participant containing the watchpoint, e.g. ./plot-displacement.sh coil-python."
+    exit 1
+fi
+
+FILE="$1/precice-ParticipantU-watchpoint-UI.log"
+
+if [ ! -f "$FILE" ]; then
+    echo "Unable to locate the watchpoint file (precice-ParticipantU-watchpoint-UI.log) in the specified participant directory '${1}'. Make sure the specified directory matches the participant you used for the calculations."
+    exit 1
+fi
+
+gnuplot -p << EOF
+	set grid
+	set title 'Voltage and current'
+	set xlabel 'time [s]'
+	set ylabel 'Voltage / Current'
+	plot "$1/precice-ParticipantU-watchpoint-UI.log" using 1:4 with linespoints title "Voltage", "$1/precice-ParticipantU-watchpoint-UI.log" using 1:5 with linespoints title "Current"
+EOF
diff --git a/resonant-circuit/precice-config.xml b/resonant-circuit/precice-config.xml
new file mode 100644
index 000000000..f09f5a42c
--- /dev/null
+++ b/resonant-circuit/precice-config.xml
@@ -0,0 +1,51 @@
+<?xml version="1.0"?>
+
+<precice-configuration>
+
+  <data:scalar name="I"  />
+  <data:scalar name="U"  />
+
+  <mesh name="MeshI" dimensions="2">
+    <use-data name="I" />
+    <use-data name="U" />
+  </mesh>
+
+  <mesh name="MeshU" dimensions="2">
+    <use-data name="U" />
+    <use-data name="I" />
+  </mesh>
+
+  <participant name="ParticipantI">
+    <provide-mesh name="MeshI" />
+    <receive-mesh name="MeshU" from="ParticipantU"/>
+    <write-data name="I"     mesh="MeshI" />
+    <read-data  name="U" mesh="MeshI" />
+    <mapping:nearest-neighbor
+      direction="write" from="MeshI" to="MeshU" constraint="consistent"/>
+    <mapping:nearest-neighbor
+      direction="read" from="MeshU" to="MeshI" constraint="consistent"/>
+  </participant>
+
+  <participant name="ParticipantU">
+    <provide-mesh name="MeshU" />
+    <receive-mesh name="MeshI" from="ParticipantI"/>
+    <write-data name="U"     mesh="MeshU" />
+    <read-data  name="I" mesh="MeshU" />
+    <mapping:nearest-neighbor
+      direction="write" from="MeshU" to="MeshI" constraint="consistent"/>
+    <mapping:nearest-neighbor
+      direction="read" from="MeshI" to="MeshU" constraint="consistent"/>
+    <watch-point mesh="MeshU" name="UI" coordinate="0.0;0.0" />
+  </participant>
+
+  <m2n:sockets acceptor="ParticipantI" connector="ParticipantU" exchange-directory=".." />
+
+  <coupling-scheme:serial-implicit>
+    <participants first="ParticipantI" second="ParticipantU" />
+    <max-time value="10" />
+    <time-window-size value=".1" />
+    <max-iterations value="3" />
+    <exchange data="I" mesh="MeshI" from="ParticipantI" to="ParticipantU"/>
+    <exchange data="U" mesh="MeshU" from="ParticipantU" to="ParticipantI"/>
+  </coupling-scheme:serial-implicit>
+</precice-configuration>

From 18737c744048e3a165cadb211ca6a28a4c3c9a84 Mon Sep 17 00:00:00 2001
From: Benjamin Rodenberg <benjamin.rodenberg@cit.tum.de>
Date: Wed, 20 Mar 2024 13:48:13 +0100
Subject: [PATCH 02/15] Add an error estimate.

---
 resonant-circuit/capacitor-python/SolverI.py | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/resonant-circuit/capacitor-python/SolverI.py b/resonant-circuit/capacitor-python/SolverI.py
index 596059cdc..3959b663e 100644
--- a/resonant-circuit/capacitor-python/SolverI.py
+++ b/resonant-circuit/capacitor-python/SolverI.py
@@ -70,6 +70,10 @@
 # Stop coupling
 participant.finalize()
 
-I_analytical = lambda t: Io*np.cos(t*w0)
-error = I0 - I_ana(10)
+import math
+
+T_max = t
+assert(math.isclose(T_max, 10))
+I_analytical = lambda t: Io*np.cos(t*w0 + phi)
+error = I0 - I_analytical(T_max)
 print(f"error: {error}")
\ No newline at end of file

From eef2da1ca91786138e2f143e71a5876aa78c3276 Mon Sep 17 00:00:00 2001
From: Benjamin Rodenberg <benjamin.rodenberg@cit.tum.de>
Date: Wed, 20 Mar 2024 15:32:30 +0100
Subject: [PATCH 03/15] Use black-box time stepper, subcycling, and time
 interpolation.

---
 oscillator/precice-config.xml                | 12 ++---
 oscillator/python/oscillator.py              | 45 ++++++++++++-----
 oscillator/python/timeSteppers.py            | 49 +++++++------------
 resonant-circuit/capacitor-python/SolverI.py | 51 +++++++++++---------
 resonant-circuit/coil-python/SolverU.py      | 45 +++++++++--------
 resonant-circuit/precice-config.xml          | 12 +++--
 6 files changed, 116 insertions(+), 98 deletions(-)

diff --git a/oscillator/precice-config.xml b/oscillator/precice-config.xml
index 384e1a86d..62e2d9ca0 100644
--- a/oscillator/precice-config.xml
+++ b/oscillator/precice-config.xml
@@ -4,8 +4,8 @@
     <sink filter="%Severity% > debug" />
   </log>
 
-  <data:scalar name="Force-Left" waveform-degree="3" />
-  <data:scalar name="Force-Right" waveform-degree="3" />
+  <data:scalar name="Force-Left" waveform-degree="5" />
+  <data:scalar name="Force-Right" waveform-degree="5" />
 
   <mesh name="Mass-Left-Mesh" dimensions="2">
     <use-data name="Force-Left" />
@@ -45,10 +45,10 @@
   <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="Force-Left" mesh="Mass-Left-Mesh" limit="1e-10" />
-    <relative-convergence-measure data="Force-Right" mesh="Mass-Left-Mesh" limit="1e-10" />
+    <time-window-size value="1" />
+    <max-iterations value="20" />
+    <relative-convergence-measure data="Force-Left" mesh="Mass-Left-Mesh" limit="1e-7" />
+    <relative-convergence-measure data="Force-Right" mesh="Mass-Left-Mesh" limit="1e-7" />
     <exchange
       data="Force-Left"
       mesh="Mass-Left-Mesh"
diff --git a/oscillator/python/oscillator.py b/oscillator/python/oscillator.py
index bf1d6be08..274262f52 100644
--- a/oscillator/python/oscillator.py
+++ b/oscillator/python/oscillator.py
@@ -80,7 +80,7 @@ class Participant(Enum):
 
 participant.initialize()
 precice_dt = participant.get_max_time_step_size()
-my_dt = precice_dt / 4  # use my_dt < precice_dt for subcycling
+my_dt = precice_dt # / 4  # use my_dt < precice_dt for subcycling
 
 # Initial Conditions
 a0 = (f0 - stiffness * u0) / mass
@@ -118,12 +118,18 @@ class Participant(Enum):
 v_write = [v]
 t_write = [t]
 
+import pandas as pd
+
+stats = {}
+
 while participant.is_coupling_ongoing():
     if participant.requires_writing_checkpoint():
         u_cp = u
         v_cp = v
         a_cp = a
         t_cp = t
+        stats[t] = {}
+        k = 0
         # store data for plotting and postprocessing
         positions += u_write
         velocities += v_write
@@ -133,22 +139,26 @@ class Participant(Enum):
     precice_dt = participant.get_max_time_step_size()
     dt = np.min([precice_dt, my_dt])
 
-    read_times = time_stepper.rhs_eval_points(dt)
-    f = len(read_times) * [None]
-
-    for i in range(len(read_times)):
-        read_data = participant.read_data(mesh_name, read_data_name, vertex_ids, read_times[i])
-        f[i] = read_data[0]
+    f = lambda t: participant.read_data(mesh_name, read_data_name, vertex_ids, t)[0]
 
     # do generalized alpha step
-    u_new, v_new, a_new = time_stepper.do_step(u, v, a, f, dt)
+    u_new, v_new, a_new, sol = time_stepper.do_step(u, v, a, f, dt)
+    stats[t][k] = sol.ts
     t_new = t + dt
+    k += 1
 
-    write_data = [connecting_spring.k * u_new]
-
-    participant.write_data(mesh_name, write_data_name, vertex_ids, write_data)
+    print(sol.ts)
 
-    participant.advance(dt)
+    n_pseudo = 10*len(sol.ts)  # 5 times no. substeps should be on the safe side.
+    t_pseudo = 0
+    print(n_pseudo)
+    for i in range(n_pseudo):
+        # perform n_pseudo pseudosteps
+        dt_pseudo = dt / n_pseudo
+        t_pseudo += dt_pseudo
+        write_data = [connecting_spring.k * sol(t_pseudo)[0]]
+        participant.write_data(mesh_name, write_data_name, vertex_ids, write_data)
+        participant.advance(dt_pseudo)
 
     if participant.requires_reading_checkpoint():
         u = u_cp
@@ -198,3 +208,14 @@ class Participant(Enum):
     csv_write.writerow(['time', 'position', 'velocity'])
     for t, u, v in zip(times, positions, velocities):
         csv_write.writerow([t, u, v])
+
+from matplotlib import pyplot as plt
+
+print(stats)
+
+for t in stats:
+    plt.title(f"{participant_name} at {t}")
+    for k in stats[t]:
+        plt.plot(stats[t][k], k*np.ones_like(stats[t][k]),'*')
+    plt.show()
+
diff --git a/oscillator/python/timeSteppers.py b/oscillator/python/timeSteppers.py
index b82f3bc78..22d865e38 100644
--- a/oscillator/python/timeSteppers.py
+++ b/oscillator/python/timeSteppers.py
@@ -30,12 +30,8 @@ def __init__(self, stiffness, mass, alpha_f=0.4, alpha_m=0.2) -> None:
         self.stiffness = stiffness
         self.mass = mass
 
-    def rhs_eval_points(self, dt) -> List[float]:
-        return [(1 - self.alpha_f) * dt]
-
-    def do_step(self, u, v, a, f, dt) -> Tuple[float, float, float]:
-        if isinstance(f, list):  # if f is list, turn it into a number
-            f = f[0]
+    def do_step(self, u, v, a, rhs, dt) -> Tuple[float, float, float]:
+        f = rhs((1 - self.alpha_f) * dt)
 
         m = 3 * [None]
         m[0] = (1 - self.alpha_m) / (self.beta * dt**2)
@@ -71,31 +67,25 @@ def __init__(self, ode_system) -> None:
         self.ode_system = ode_system
         pass
 
-    def rhs_eval_points(self, dt) -> List[float]:
-        return [self.c[0] * dt, self.c[1] * dt, self.c[2] * dt, self.c[3] * dt]
-
-    def do_step(self, u, v, a, f, dt) -> Tuple[float, float, float]:
+    def do_step(self, u, v, a, rhs, dt) -> Tuple[float, float, float]:
         assert (isinstance(u, type(v)))
 
         n_stages = 4
 
-        if isinstance(f, numbers.Number):  # if f is number, assume constant f
-            f = n_stages * [f]
-
         if isinstance(u, np.ndarray):
             x = np.concatenate([u, v])
-            rhs = [np.concatenate([np.array([0, 0]), f[i]]) for i in range(n_stages)]
+            f = lambda t: np.concatenate([np.array([0, 0]), rhs(t)])
         elif isinstance(u, numbers.Number):
             x = np.array([u, v])
-            rhs = [np.array([0, f[i]]) for i in range(n_stages)]
+            f = lambda t: np.array([0, rhs(t)])
         else:
             raise Exception(f"Cannot handle input type {type(u)} of u and v")
 
         s = n_stages * [None]
-        s[0] = self.ode_system.dot(x) + rhs[0]
-        s[1] = self.ode_system.dot(x + self.a[1, 0] * s[0] * dt) + rhs[1]
-        s[2] = self.ode_system.dot(x + self.a[2, 1] * s[1] * dt) + rhs[2]
-        s[3] = self.ode_system.dot(x + self.a[3, 2] * s[2] * dt) + rhs[3]
+        s[0] = self.ode_system.dot(x) + f(self.c[0] * dt)
+        s[1] = self.ode_system.dot(x + self.a[1, 0] * s[0] * dt) + f(self.c[1] * dt)
+        s[2] = self.ode_system.dot(x + self.a[2, 1] * s[1] * dt) + f(self.c[2] * dt)
+        s[3] = self.ode_system.dot(x + self.a[3, 2] * s[2] * dt) + f(self.c[3] * dt)
 
         x_new = x
 
@@ -119,14 +109,9 @@ def __init__(self, ode_system) -> None:
         self.ode_system = ode_system
         pass
 
-    def rhs_eval_points(self, dt) -> List[float]:
-        return np.linspace(0, dt, 5)  # will create an interpolant from this later
-
-    def do_step(self, u, v, a, f, dt) -> Tuple[float, float, float]:
+    def do_step(self, u, v, a, rhs, dt) -> Tuple[float, float, float]:
         from brot.interpolation import do_lagrange_interpolation
 
-        ts = self.rhs_eval_points(dt)
-
         t0 = 0
 
         assert (isinstance(u, type(v)))
@@ -135,20 +120,21 @@ def do_step(self, u, v, a, f, dt) -> Tuple[float, float, float]:
             x0 = np.concatenate([u, v])
             f = np.array(f)
             assert (u.shape[0] == f.shape[1])
-            def rhs_fun(t, x): return np.concatenate([np.array([np.zeros_like(t), np.zeros_like(t)]), [
-                do_lagrange_interpolation(t, ts, f[:, i]) for i in range(u.shape[0])]])
+            def rhs_fun(t): return np.concatenate([np.array([np.zeros_like(t), np.zeros_like(t)]), rhs(t)])
         elif isinstance(u, numbers.Number):
             x0 = np.array([u, v])
-            def rhs_fun(t, x): return np.array([np.zeros_like(t), do_lagrange_interpolation(t, ts, f)])
+            def rhs_fun(t): return np.array([np.zeros_like(t), rhs(t)])
         else:
             raise Exception(f"Cannot handle input type {type(u)} of u and v")
 
         def fun(t, x):
-            return self.ode_system.dot(x) + rhs_fun(t, x)
+            return self.ode_system.dot(x) + rhs_fun(t)
 
         # use large rtol and atol to circumvent error control.
+        # ret = sp.integrate.solve_ivp(fun, [t0, t0 + dt], x0, method="Radau",
+        #                              first_step=dt, max_step=dt, rtol=10e10, atol=10e10)
         ret = sp.integrate.solve_ivp(fun, [t0, t0 + dt], x0, method="Radau",
-                                     first_step=dt, max_step=dt, rtol=10e10, atol=10e10)
+                                     dense_output=True, rtol=10e-5, atol=10e-9)
 
         a_new = None
         if isinstance(u, np.ndarray):
@@ -156,4 +142,5 @@ def fun(t, x):
         elif isinstance(u, numbers.Number):
             u_new, v_new = ret.y[:, -1]
 
-        return u_new, v_new, a_new
+
+        return u_new, v_new, a_new, ret.sol
diff --git a/resonant-circuit/capacitor-python/SolverI.py b/resonant-circuit/capacitor-python/SolverI.py
index 3959b663e..1c65d5396 100644
--- a/resonant-circuit/capacitor-python/SolverI.py
+++ b/resonant-circuit/capacitor-python/SolverI.py
@@ -21,59 +21,62 @@
 L = 1                      # Inductance of the coil
 t0 = 0                     # Initial simulation time
 t_max = 10                 # End simulation time
-Io = 1                     # Initial current
+Io = np.array([1])                     # Initial current
 phi = 0                    # Phase of the signal
 
 w0 = 1/np.sqrt(L*C)           # Resonant frequency
 I0 = Io*np.cos(phi)           # Initial condition for I
-U0 = -w0*L*Io*np.sin(phi)     # Initial condition for U
 
-f_I = lambda t, I, U: U/L;       # Time derivative of I; ODE determining capacitor
+def f_I(t_span, t):
+    if t < t_span[1]:
+        dt = t - t_span[0]
+    else:
+        dt = t_span[1] - t_span[0]
+
+    U = participant.read_data(mesh_name, read_data_name, vertex_ids, dt)
+    return U/L;       # Time derivative of I; ODE determining capacitor
 
 # Initialize simulation
 if participant.requires_initial_data():
-    participant.write_data(mesh_name, write_data_name, vertex_ids, np.array([I0]))
+    participant.write_data(mesh_name, write_data_name, vertex_ids, I0)
 
 participant.initialize()
 
-dt = participant.get_max_time_step_size()
-
 # Start simulation
-t = t0 + dt
+t = t0
 while participant.is_coupling_ongoing():
 
     # Record checkpoint if necessary
     if participant.requires_writing_checkpoint():
         I0_checkpoint = I0
-        U0_checkpoint = U0
+        t_checkpoint = t
 
     # Make Simulation Step
-    sol = solve_ivp(lambda t, y: f_I(t, y, U0), [t0, t], np.array([I0]))
-    I0 = sol.y[-1,-1]
-    print(sol.y)
-    print(I0)
+    dt = participant.get_max_time_step_size()
+    t_span = [t, t+dt]
+    print(f"solving time {t_span}")
+    sol = solve_ivp(lambda t, y: f_I(t_span, t), t_span, I0, dense_output=True,r_tol=1e-6, a_tol=1e-9)
 
     # Exchange data
-    participant.write_data(mesh_name, write_data_name, vertex_ids, np.array([I0]))
-    participant.advance(dt)
+    evals = 10
+    for i in range(evals):
+        I0 = sol.sol(t_span[0] + (i+1)*dt/evals)
+        participant.write_data(mesh_name, write_data_name, vertex_ids, np.array(I0))
+        participant.advance(dt/evals)
 
-    # Recover checkpoint if not converged, else finish time step
+    t = t + dt
+
+    # Recover checkpoint if not converged
     if participant.requires_reading_checkpoint():
         I0 = I0_checkpoint
-        U0 = U0_checkpoint
-    else:
-        dt = participant.get_max_time_step_size()
-        U0 = participant.read_data(mesh_name, read_data_name, vertex_ids, dt)[0]
-        t0 = t
-        t = t0 + dt
+        t = t_checkpoint
 
 # Stop coupling
 participant.finalize()
 
 import math
 
-T_max = t
-assert(math.isclose(T_max, 10))
+assert(math.isclose(t, t_max))
 I_analytical = lambda t: Io*np.cos(t*w0 + phi)
-error = I0 - I_analytical(T_max)
+error = I0 - I_analytical(t_max)
 print(f"error: {error}")
\ No newline at end of file
diff --git a/resonant-circuit/coil-python/SolverU.py b/resonant-circuit/coil-python/SolverU.py
index f11a0ce61..848ef1cba 100644
--- a/resonant-circuit/coil-python/SolverU.py
+++ b/resonant-circuit/coil-python/SolverU.py
@@ -21,50 +21,55 @@
 L = 1                      # Inductance of the coil
 t0 = 0                     # Initial simulation time
 t_max = 10                 # End simulation time
-Io = 1                     # Initial current
+Io = np.array([1])         # Initial current
 phi = 0                    # Phase of the signal
 
 w0 = 1/np.sqrt(L*C)           # Resonant frequency
-I0 = Io*np.cos(phi)           # Initial condition for I
 U0 = -w0*L*Io*np.sin(phi)     # Initial condition for U
 
-f_U = lambda t, U, I: -I/C      # Time derivative of U
+def f_U(t_span, t):
+    if t < t_span[1]:
+        dt = t - t_span[0]
+    else:
+        dt = t_span[1] - t_span[0]
+
+    I = participant.read_data(mesh_name, read_data_name, vertex_ids, dt)
+    return -I/C      # Time derivative of U
 
 # Initialize simulation
 if participant.requires_initial_data():
-    participant.write_data(mesh_name, write_data_name, vertex_ids, np.array([U0]))
+    participant.write_data(mesh_name, write_data_name, vertex_ids, U0)
 
 participant.initialize()
 
-dt = participant.get_max_time_step_size()
-
 # Start simulation
-t = t0 + dt
+t = t0
 while participant.is_coupling_ongoing():
 
     # Record checkpoint if necessary
     if participant.requires_writing_checkpoint():
-        I0_checkpoint = I0
         U0_checkpoint = U0
+        t_checkpoint = t
 
     # Make Simulation Step
-    sol = solve_ivp(lambda t, y: f_U(t, y, I0), [t0, t], np.array([U0]))
-    U0 = sol.y[-1,-1]
-    print(U0)
+    dt = participant.get_max_time_step_size()
+    t_span = [t, t+dt]
+    print(f"solving time {t_span}")
+    sol = solve_ivp(lambda t, y: f_U(t_span, t), t_span, U0, dense_output=True,r_tol=1e-6, a_tol=1e-9)
 
     # Exchange data
-    participant.write_data(mesh_name, write_data_name, vertex_ids, np.array([U0]))
-    participant.advance(dt)
+    evals = 10
+    for i in range(evals):
+        U0 = sol.sol(t_span[0] + (i+1)*dt/evals)
+        participant.write_data(mesh_name, write_data_name, vertex_ids, np.array(U0))
+        participant.advance(dt/evals)
 
-    # Recover checkpoint if not converged, else finish time step
+    t = t + dt
+
+    # Recover checkpoint if not converged
     if participant.requires_reading_checkpoint():
-        I0 = I0_checkpoint
         U0 = U0_checkpoint
-    else:
-        dt = participant.get_max_time_step_size()
-        I0 = participant.read_data(mesh_name, read_data_name, vertex_ids, dt)[0]
-        t0 = t
-        t = t0 + dt
+        t = t_checkpoint
 
 # Stop coupling
 participant.finalize()
diff --git a/resonant-circuit/precice-config.xml b/resonant-circuit/precice-config.xml
index f09f5a42c..7143e45e8 100644
--- a/resonant-circuit/precice-config.xml
+++ b/resonant-circuit/precice-config.xml
@@ -2,8 +2,8 @@
 
 <precice-configuration>
 
-  <data:scalar name="I"  />
-  <data:scalar name="U"  />
+  <data:scalar name="I" waveform-degree="5" />
+  <data:scalar name="U" waveform-degree="5" />
 
   <mesh name="MeshI" dimensions="2">
     <use-data name="I" />
@@ -44,8 +44,10 @@
     <participants first="ParticipantI" second="ParticipantU" />
     <max-time value="10" />
     <time-window-size value=".1" />
-    <max-iterations value="3" />
-    <exchange data="I" mesh="MeshI" from="ParticipantI" to="ParticipantU"/>
-    <exchange data="U" mesh="MeshU" from="ParticipantU" to="ParticipantI"/>
+    <max-iterations value="20" />
+    <exchange data="I" mesh="MeshI" from="ParticipantI" to="ParticipantU" substeps="true" />
+    <exchange data="U" mesh="MeshU" from="ParticipantU" to="ParticipantI" substeps="true" />
+    <relative-convergence-measure limit="1e-8" data="U" mesh="MeshU" />
+    <relative-convergence-measure limit="1e-8" data="I" mesh="MeshI" />
   </coupling-scheme:serial-implicit>
 </precice-configuration>

From 44eda69828ecd107efe8d1709fe5e963a4f8ca30 Mon Sep 17 00:00:00 2001
From: Benjamin Rodenberg <benjamin.rodenberg@cit.tum.de>
Date: Fri, 22 Mar 2024 17:24:21 +0100
Subject: [PATCH 04/15] Fix black-box adaptive time stepping

---
 resonant-circuit/capacitor-python/SolverI.py | 29 ++++++++++++--------
 resonant-circuit/coil-python/SolverU.py      | 18 ++++++------
 resonant-circuit/precice-config.xml          |  6 ++--
 3 files changed, 30 insertions(+), 23 deletions(-)

diff --git a/resonant-circuit/capacitor-python/SolverI.py b/resonant-circuit/capacitor-python/SolverI.py
index 1c65d5396..92b771729 100644
--- a/resonant-circuit/capacitor-python/SolverI.py
+++ b/resonant-circuit/capacitor-python/SolverI.py
@@ -20,18 +20,22 @@
 C = 2                      # Capacitance of the capacitor
 L = 1                      # Inductance of the coil
 t0 = 0                     # Initial simulation time
-t_max = 10                 # End simulation time
-Io = np.array([1])                     # Initial current
+t_max = 1                  # End simulation time
+Io = np.array([1])         # Initial current
 phi = 0                    # Phase of the signal
 
 w0 = 1/np.sqrt(L*C)           # Resonant frequency
 I0 = Io*np.cos(phi)           # Initial condition for I
 
-def f_I(t_span, t):
-    if t < t_span[1]:
-        dt = t - t_span[0]
-    else:
-        dt = t_span[1] - t_span[0]
+# to estimate cost
+global f_evals
+f_evals = 0
+
+def f_I(dt, max_allowed_dt):
+    global f_evals
+    f_evals += 1
+    if dt > max_allowed_dt:  # read_data will fail, if dt is outside of window
+        return np.nan
 
     U = participant.read_data(mesh_name, read_data_name, vertex_ids, dt)
     return U/L;       # Time derivative of I; ODE determining capacitor
@@ -42,6 +46,8 @@ def f_I(t_span, t):
 
 participant.initialize()
 
+solver_dt = participant.get_max_time_step_size()
+
 # Start simulation
 t = t0
 while participant.is_coupling_ongoing():
@@ -52,10 +58,10 @@ def f_I(t_span, t):
         t_checkpoint = t
 
     # Make Simulation Step
-    dt = participant.get_max_time_step_size()
+    precice_dt = participant.get_max_time_step_size()
+    dt = min([precice_dt, solver_dt])
     t_span = [t, t+dt]
-    print(f"solving time {t_span}")
-    sol = solve_ivp(lambda t, y: f_I(t_span, t), t_span, I0, dense_output=True,r_tol=1e-6, a_tol=1e-9)
+    sol = solve_ivp(lambda t, y: f_I(t-t_span[0], dt), t_span, I0, dense_output=True,r_tol=1e-12, a_tol=1e-12)
 
     # Exchange data
     evals = 10
@@ -79,4 +85,5 @@ def f_I(t_span, t):
 assert(math.isclose(t, t_max))
 I_analytical = lambda t: Io*np.cos(t*w0 + phi)
 error = I0 - I_analytical(t_max)
-print(f"error: {error}")
\ No newline at end of file
+print(f"{error=}")
+print(f"{f_evals=}")
\ No newline at end of file
diff --git a/resonant-circuit/coil-python/SolverU.py b/resonant-circuit/coil-python/SolverU.py
index 848ef1cba..15b48c8e6 100644
--- a/resonant-circuit/coil-python/SolverU.py
+++ b/resonant-circuit/coil-python/SolverU.py
@@ -27,11 +27,9 @@
 w0 = 1/np.sqrt(L*C)           # Resonant frequency
 U0 = -w0*L*Io*np.sin(phi)     # Initial condition for U
 
-def f_U(t_span, t):
-    if t < t_span[1]:
-        dt = t - t_span[0]
-    else:
-        dt = t_span[1] - t_span[0]
+def f_U(dt, max_allowed_dt):
+    if dt > max_allowed_dt:  # read_data will fail, if dt is outside of window
+        return np.nan
 
     I = participant.read_data(mesh_name, read_data_name, vertex_ids, dt)
     return -I/C      # Time derivative of U
@@ -42,6 +40,8 @@ def f_U(t_span, t):
 
 participant.initialize()
 
+solver_dt = participant.get_max_time_step_size()
+
 # Start simulation
 t = t0
 while participant.is_coupling_ongoing():
@@ -52,10 +52,10 @@ def f_U(t_span, t):
         t_checkpoint = t
 
     # Make Simulation Step
-    dt = participant.get_max_time_step_size()
+    precice_dt = participant.get_max_time_step_size()
+    dt = min([precice_dt, solver_dt])
     t_span = [t, t+dt]
-    print(f"solving time {t_span}")
-    sol = solve_ivp(lambda t, y: f_U(t_span, t), t_span, U0, dense_output=True,r_tol=1e-6, a_tol=1e-9)
+    sol = solve_ivp(lambda t, y: f_U(t-t_span[0], dt), t_span, U0, dense_output=True,r_tol=1e-12, a_tol=1e-12)
 
     # Exchange data
     evals = 10
@@ -72,4 +72,4 @@ def f_U(t_span, t):
         t = t_checkpoint
 
 # Stop coupling
-participant.finalize()
+participant.finalize()
\ No newline at end of file
diff --git a/resonant-circuit/precice-config.xml b/resonant-circuit/precice-config.xml
index 7143e45e8..777cd2890 100644
--- a/resonant-circuit/precice-config.xml
+++ b/resonant-circuit/precice-config.xml
@@ -2,8 +2,8 @@
 
 <precice-configuration>
 
-  <data:scalar name="I" waveform-degree="5" />
-  <data:scalar name="U" waveform-degree="5" />
+  <data:scalar name="I" waveform-degree="3" />
+  <data:scalar name="U" waveform-degree="3" />
 
   <mesh name="MeshI" dimensions="2">
     <use-data name="I" />
@@ -42,7 +42,7 @@
 
   <coupling-scheme:serial-implicit>
     <participants first="ParticipantI" second="ParticipantU" />
-    <max-time value="10" />
+    <max-time value="1" />
     <time-window-size value=".1" />
     <max-iterations value="20" />
     <exchange data="I" mesh="MeshI" from="ParticipantI" to="ParticipantU" substeps="true" />

From 234f9111f4e92ee7fac23981d89b6310e0f2ab47 Mon Sep 17 00:00:00 2001
From: Benjamin Rodenberg <benjamin.rodenberg@cit.tum.de>
Date: Mon, 25 Mar 2024 13:41:56 +0100
Subject: [PATCH 05/15] Remove unneeded.

---
 oscillator/precice-config.xml     | 12 ++++----
 oscillator/python/oscillator.py   | 45 ++++++++--------------------
 oscillator/python/timeSteppers.py | 49 +++++++++++++++++++------------
 3 files changed, 49 insertions(+), 57 deletions(-)

diff --git a/oscillator/precice-config.xml b/oscillator/precice-config.xml
index 62e2d9ca0..384e1a86d 100644
--- a/oscillator/precice-config.xml
+++ b/oscillator/precice-config.xml
@@ -4,8 +4,8 @@
     <sink filter="%Severity% > debug" />
   </log>
 
-  <data:scalar name="Force-Left" waveform-degree="5" />
-  <data:scalar name="Force-Right" waveform-degree="5" />
+  <data:scalar name="Force-Left" waveform-degree="3" />
+  <data:scalar name="Force-Right" waveform-degree="3" />
 
   <mesh name="Mass-Left-Mesh" dimensions="2">
     <use-data name="Force-Left" />
@@ -45,10 +45,10 @@
   <coupling-scheme:serial-implicit>
     <participants first="Mass-Left" second="Mass-Right" />
     <max-time value="1" />
-    <time-window-size value="1" />
-    <max-iterations value="20" />
-    <relative-convergence-measure data="Force-Left" mesh="Mass-Left-Mesh" limit="1e-7" />
-    <relative-convergence-measure data="Force-Right" mesh="Mass-Left-Mesh" limit="1e-7" />
+    <time-window-size value="0.01" />
+    <max-iterations value="200" />
+    <relative-convergence-measure data="Force-Left" mesh="Mass-Left-Mesh" limit="1e-10" />
+    <relative-convergence-measure data="Force-Right" mesh="Mass-Left-Mesh" limit="1e-10" />
     <exchange
       data="Force-Left"
       mesh="Mass-Left-Mesh"
diff --git a/oscillator/python/oscillator.py b/oscillator/python/oscillator.py
index 274262f52..bf1d6be08 100644
--- a/oscillator/python/oscillator.py
+++ b/oscillator/python/oscillator.py
@@ -80,7 +80,7 @@ class Participant(Enum):
 
 participant.initialize()
 precice_dt = participant.get_max_time_step_size()
-my_dt = precice_dt # / 4  # use my_dt < precice_dt for subcycling
+my_dt = precice_dt / 4  # use my_dt < precice_dt for subcycling
 
 # Initial Conditions
 a0 = (f0 - stiffness * u0) / mass
@@ -118,18 +118,12 @@ class Participant(Enum):
 v_write = [v]
 t_write = [t]
 
-import pandas as pd
-
-stats = {}
-
 while participant.is_coupling_ongoing():
     if participant.requires_writing_checkpoint():
         u_cp = u
         v_cp = v
         a_cp = a
         t_cp = t
-        stats[t] = {}
-        k = 0
         # store data for plotting and postprocessing
         positions += u_write
         velocities += v_write
@@ -139,26 +133,22 @@ class Participant(Enum):
     precice_dt = participant.get_max_time_step_size()
     dt = np.min([precice_dt, my_dt])
 
-    f = lambda t: participant.read_data(mesh_name, read_data_name, vertex_ids, t)[0]
+    read_times = time_stepper.rhs_eval_points(dt)
+    f = len(read_times) * [None]
+
+    for i in range(len(read_times)):
+        read_data = participant.read_data(mesh_name, read_data_name, vertex_ids, read_times[i])
+        f[i] = read_data[0]
 
     # do generalized alpha step
-    u_new, v_new, a_new, sol = time_stepper.do_step(u, v, a, f, dt)
-    stats[t][k] = sol.ts
+    u_new, v_new, a_new = time_stepper.do_step(u, v, a, f, dt)
     t_new = t + dt
-    k += 1
 
-    print(sol.ts)
+    write_data = [connecting_spring.k * u_new]
+
+    participant.write_data(mesh_name, write_data_name, vertex_ids, write_data)
 
-    n_pseudo = 10*len(sol.ts)  # 5 times no. substeps should be on the safe side.
-    t_pseudo = 0
-    print(n_pseudo)
-    for i in range(n_pseudo):
-        # perform n_pseudo pseudosteps
-        dt_pseudo = dt / n_pseudo
-        t_pseudo += dt_pseudo
-        write_data = [connecting_spring.k * sol(t_pseudo)[0]]
-        participant.write_data(mesh_name, write_data_name, vertex_ids, write_data)
-        participant.advance(dt_pseudo)
+    participant.advance(dt)
 
     if participant.requires_reading_checkpoint():
         u = u_cp
@@ -208,14 +198,3 @@ class Participant(Enum):
     csv_write.writerow(['time', 'position', 'velocity'])
     for t, u, v in zip(times, positions, velocities):
         csv_write.writerow([t, u, v])
-
-from matplotlib import pyplot as plt
-
-print(stats)
-
-for t in stats:
-    plt.title(f"{participant_name} at {t}")
-    for k in stats[t]:
-        plt.plot(stats[t][k], k*np.ones_like(stats[t][k]),'*')
-    plt.show()
-
diff --git a/oscillator/python/timeSteppers.py b/oscillator/python/timeSteppers.py
index 22d865e38..b82f3bc78 100644
--- a/oscillator/python/timeSteppers.py
+++ b/oscillator/python/timeSteppers.py
@@ -30,8 +30,12 @@ def __init__(self, stiffness, mass, alpha_f=0.4, alpha_m=0.2) -> None:
         self.stiffness = stiffness
         self.mass = mass
 
-    def do_step(self, u, v, a, rhs, dt) -> Tuple[float, float, float]:
-        f = rhs((1 - self.alpha_f) * dt)
+    def rhs_eval_points(self, dt) -> List[float]:
+        return [(1 - self.alpha_f) * dt]
+
+    def do_step(self, u, v, a, f, dt) -> Tuple[float, float, float]:
+        if isinstance(f, list):  # if f is list, turn it into a number
+            f = f[0]
 
         m = 3 * [None]
         m[0] = (1 - self.alpha_m) / (self.beta * dt**2)
@@ -67,25 +71,31 @@ def __init__(self, ode_system) -> None:
         self.ode_system = ode_system
         pass
 
-    def do_step(self, u, v, a, rhs, dt) -> Tuple[float, float, float]:
+    def rhs_eval_points(self, dt) -> List[float]:
+        return [self.c[0] * dt, self.c[1] * dt, self.c[2] * dt, self.c[3] * dt]
+
+    def do_step(self, u, v, a, f, dt) -> Tuple[float, float, float]:
         assert (isinstance(u, type(v)))
 
         n_stages = 4
 
+        if isinstance(f, numbers.Number):  # if f is number, assume constant f
+            f = n_stages * [f]
+
         if isinstance(u, np.ndarray):
             x = np.concatenate([u, v])
-            f = lambda t: np.concatenate([np.array([0, 0]), rhs(t)])
+            rhs = [np.concatenate([np.array([0, 0]), f[i]]) for i in range(n_stages)]
         elif isinstance(u, numbers.Number):
             x = np.array([u, v])
-            f = lambda t: np.array([0, rhs(t)])
+            rhs = [np.array([0, f[i]]) for i in range(n_stages)]
         else:
             raise Exception(f"Cannot handle input type {type(u)} of u and v")
 
         s = n_stages * [None]
-        s[0] = self.ode_system.dot(x) + f(self.c[0] * dt)
-        s[1] = self.ode_system.dot(x + self.a[1, 0] * s[0] * dt) + f(self.c[1] * dt)
-        s[2] = self.ode_system.dot(x + self.a[2, 1] * s[1] * dt) + f(self.c[2] * dt)
-        s[3] = self.ode_system.dot(x + self.a[3, 2] * s[2] * dt) + f(self.c[3] * dt)
+        s[0] = self.ode_system.dot(x) + rhs[0]
+        s[1] = self.ode_system.dot(x + self.a[1, 0] * s[0] * dt) + rhs[1]
+        s[2] = self.ode_system.dot(x + self.a[2, 1] * s[1] * dt) + rhs[2]
+        s[3] = self.ode_system.dot(x + self.a[3, 2] * s[2] * dt) + rhs[3]
 
         x_new = x
 
@@ -109,9 +119,14 @@ def __init__(self, ode_system) -> None:
         self.ode_system = ode_system
         pass
 
-    def do_step(self, u, v, a, rhs, dt) -> Tuple[float, float, float]:
+    def rhs_eval_points(self, dt) -> List[float]:
+        return np.linspace(0, dt, 5)  # will create an interpolant from this later
+
+    def do_step(self, u, v, a, f, dt) -> Tuple[float, float, float]:
         from brot.interpolation import do_lagrange_interpolation
 
+        ts = self.rhs_eval_points(dt)
+
         t0 = 0
 
         assert (isinstance(u, type(v)))
@@ -120,21 +135,20 @@ def do_step(self, u, v, a, rhs, dt) -> Tuple[float, float, float]:
             x0 = np.concatenate([u, v])
             f = np.array(f)
             assert (u.shape[0] == f.shape[1])
-            def rhs_fun(t): return np.concatenate([np.array([np.zeros_like(t), np.zeros_like(t)]), rhs(t)])
+            def rhs_fun(t, x): return np.concatenate([np.array([np.zeros_like(t), np.zeros_like(t)]), [
+                do_lagrange_interpolation(t, ts, f[:, i]) for i in range(u.shape[0])]])
         elif isinstance(u, numbers.Number):
             x0 = np.array([u, v])
-            def rhs_fun(t): return np.array([np.zeros_like(t), rhs(t)])
+            def rhs_fun(t, x): return np.array([np.zeros_like(t), do_lagrange_interpolation(t, ts, f)])
         else:
             raise Exception(f"Cannot handle input type {type(u)} of u and v")
 
         def fun(t, x):
-            return self.ode_system.dot(x) + rhs_fun(t)
+            return self.ode_system.dot(x) + rhs_fun(t, x)
 
         # use large rtol and atol to circumvent error control.
-        # ret = sp.integrate.solve_ivp(fun, [t0, t0 + dt], x0, method="Radau",
-        #                              first_step=dt, max_step=dt, rtol=10e10, atol=10e10)
         ret = sp.integrate.solve_ivp(fun, [t0, t0 + dt], x0, method="Radau",
-                                     dense_output=True, rtol=10e-5, atol=10e-9)
+                                     first_step=dt, max_step=dt, rtol=10e10, atol=10e10)
 
         a_new = None
         if isinstance(u, np.ndarray):
@@ -142,5 +156,4 @@ def fun(t, x):
         elif isinstance(u, numbers.Number):
             u_new, v_new = ret.y[:, -1]
 
-
-        return u_new, v_new, a_new, ret.sol
+        return u_new, v_new, a_new

From 19dd2c3e524fd4ebeacbd6695bffe998514d536c Mon Sep 17 00:00:00 2001
From: Benjamin Rodenberg <benjamin.rodenberg@cit.tum.de>
Date: Mon, 25 Mar 2024 15:39:12 +0100
Subject: [PATCH 06/15] Fix format.

---
 resonant-circuit/capacitor-python/SolverI.py | 26 +++++++++++---------
 resonant-circuit/coil-python/SolverU.py      | 20 ++++++++-------
 2 files changed, 26 insertions(+), 20 deletions(-)

diff --git a/resonant-circuit/capacitor-python/SolverI.py b/resonant-circuit/capacitor-python/SolverI.py
index 92b771729..91a534b91 100644
--- a/resonant-circuit/capacitor-python/SolverI.py
+++ b/resonant-circuit/capacitor-python/SolverI.py
@@ -6,7 +6,7 @@
 participant = precice.Participant("ParticipantI", "../precice-config.xml", 0, 1)
 
 # Geometry IDs. As it is a 0-D simulation, only one vertex is necessary.
-mesh_name ="MeshI"
+mesh_name = "MeshI"
 
 dimensions = participant.get_mesh_dimensions(mesh_name)
 
@@ -24,13 +24,14 @@
 Io = np.array([1])         # Initial current
 phi = 0                    # Phase of the signal
 
-w0 = 1/np.sqrt(L*C)           # Resonant frequency
-I0 = Io*np.cos(phi)           # Initial condition for I
+w0 = 1 / np.sqrt(L * C)           # Resonant frequency
+I0 = Io * np.cos(phi)           # Initial condition for I
 
 # to estimate cost
 global f_evals
 f_evals = 0
 
+
 def f_I(dt, max_allowed_dt):
     global f_evals
     f_evals += 1
@@ -38,7 +39,8 @@ def f_I(dt, max_allowed_dt):
         return np.nan
 
     U = participant.read_data(mesh_name, read_data_name, vertex_ids, dt)
-    return U/L;       # Time derivative of I; ODE determining capacitor
+    return U / L       # Time derivative of I; ODE determining capacitor
+
 
 # Initialize simulation
 if participant.requires_initial_data():
@@ -60,15 +62,15 @@ def f_I(dt, max_allowed_dt):
     # Make Simulation Step
     precice_dt = participant.get_max_time_step_size()
     dt = min([precice_dt, solver_dt])
-    t_span = [t, t+dt]
-    sol = solve_ivp(lambda t, y: f_I(t-t_span[0], dt), t_span, I0, dense_output=True,r_tol=1e-12, a_tol=1e-12)
+    t_span = [t, t + dt]
+    sol = solve_ivp(lambda t, y: f_I(t - t_span[0], dt), t_span, I0, dense_output=True, r_tol=1e-12, a_tol=1e-12)
 
     # Exchange data
     evals = 10
     for i in range(evals):
-        I0 = sol.sol(t_span[0] + (i+1)*dt/evals)
+        I0 = sol.sol(t_span[0] + (i + 1) * dt / evals)
         participant.write_data(mesh_name, write_data_name, vertex_ids, np.array(I0))
-        participant.advance(dt/evals)
+        participant.advance(dt / evals)
 
     t = t + dt
 
@@ -82,8 +84,10 @@ def f_I(dt, max_allowed_dt):
 
 import math
 
-assert(math.isclose(t, t_max))
-I_analytical = lambda t: Io*np.cos(t*w0 + phi)
+assert (math.isclose(t, t_max))
+def I_analytical(t): return Io * np.cos(t * w0 + phi)
+
+
 error = I0 - I_analytical(t_max)
 print(f"{error=}")
-print(f"{f_evals=}")
\ No newline at end of file
+print(f"{f_evals=}")
diff --git a/resonant-circuit/coil-python/SolverU.py b/resonant-circuit/coil-python/SolverU.py
index 15b48c8e6..e0b1ba28b 100644
--- a/resonant-circuit/coil-python/SolverU.py
+++ b/resonant-circuit/coil-python/SolverU.py
@@ -6,7 +6,7 @@
 participant = precice.Participant("ParticipantU", "../precice-config.xml", 0, 1)
 
 # Geometry IDs. As it is a 0-D simulation, only one vertex is necessary.
-mesh_name ="MeshU"
+mesh_name = "MeshU"
 
 dimensions = participant.get_mesh_dimensions(mesh_name)
 
@@ -24,15 +24,17 @@
 Io = np.array([1])         # Initial current
 phi = 0                    # Phase of the signal
 
-w0 = 1/np.sqrt(L*C)           # Resonant frequency
-U0 = -w0*L*Io*np.sin(phi)     # Initial condition for U
+w0 = 1 / np.sqrt(L * C)           # Resonant frequency
+U0 = -w0 * L * Io * np.sin(phi)     # Initial condition for U
+
 
 def f_U(dt, max_allowed_dt):
     if dt > max_allowed_dt:  # read_data will fail, if dt is outside of window
         return np.nan
 
     I = participant.read_data(mesh_name, read_data_name, vertex_ids, dt)
-    return -I/C      # Time derivative of U
+    return -I / C      # Time derivative of U
+
 
 # Initialize simulation
 if participant.requires_initial_data():
@@ -54,15 +56,15 @@ def f_U(dt, max_allowed_dt):
     # Make Simulation Step
     precice_dt = participant.get_max_time_step_size()
     dt = min([precice_dt, solver_dt])
-    t_span = [t, t+dt]
-    sol = solve_ivp(lambda t, y: f_U(t-t_span[0], dt), t_span, U0, dense_output=True,r_tol=1e-12, a_tol=1e-12)
+    t_span = [t, t + dt]
+    sol = solve_ivp(lambda t, y: f_U(t - t_span[0], dt), t_span, U0, dense_output=True, r_tol=1e-12, a_tol=1e-12)
 
     # Exchange data
     evals = 10
     for i in range(evals):
-        U0 = sol.sol(t_span[0] + (i+1)*dt/evals)
+        U0 = sol.sol(t_span[0] + (i + 1) * dt / evals)
         participant.write_data(mesh_name, write_data_name, vertex_ids, np.array(U0))
-        participant.advance(dt/evals)
+        participant.advance(dt / evals)
 
     t = t + dt
 
@@ -72,4 +74,4 @@ def f_U(dt, max_allowed_dt):
         t = t_checkpoint
 
 # Stop coupling
-participant.finalize()
\ No newline at end of file
+participant.finalize()

From 88733c5e5cf62e2ac45a4aa55bb36fc0c62db8a7 Mon Sep 17 00:00:00 2001
From: Benjamin Rodenberg <benjamin.rodenberg@cit.tum.de>
Date: Mon, 17 Jun 2024 13:55:29 +0200
Subject: [PATCH 07/15] Fix labels, use time interpolation, measure
 convergence.

---
 .../capacitor.py}                             | 14 ++++++------
 resonant-circuit/capacitor-python/run.sh      |  9 ++++++++
 .../SolverI.py => coil-python/coil.py}        | 22 +++++++------------
 resonant-circuit/coil-python/run.sh           |  9 ++++++++
 resonant-circuit/precice-config.xml           | 14 +++++++-----
 5 files changed, 41 insertions(+), 27 deletions(-)
 rename resonant-circuit/{coil-python/SolverU.py => capacitor-python/capacitor.py} (86%)
 create mode 100755 resonant-circuit/capacitor-python/run.sh
 rename resonant-circuit/{capacitor-python/SolverI.py => coil-python/coil.py} (86%)
 create mode 100755 resonant-circuit/coil-python/run.sh

diff --git a/resonant-circuit/coil-python/SolverU.py b/resonant-circuit/capacitor-python/capacitor.py
similarity index 86%
rename from resonant-circuit/coil-python/SolverU.py
rename to resonant-circuit/capacitor-python/capacitor.py
index e0b1ba28b..bd2639bf7 100644
--- a/resonant-circuit/coil-python/SolverU.py
+++ b/resonant-circuit/capacitor-python/capacitor.py
@@ -3,18 +3,18 @@
 from scipy.integrate import solve_ivp
 
 # Initialize and configure preCICE
-participant = precice.Participant("ParticipantU", "../precice-config.xml", 0, 1)
+participant = precice.Participant("Capacitor", "../precice-config.xml", 0, 1)
 
 # Geometry IDs. As it is a 0-D simulation, only one vertex is necessary.
-mesh_name = "MeshU"
+mesh_name = "Capacitor-Mesh"
 
 dimensions = participant.get_mesh_dimensions(mesh_name)
 
 vertex_ids = np.array([participant.set_mesh_vertex(mesh_name, np.zeros(dimensions))])
 
 # Data IDs
-read_data_name = "I"
-write_data_name = "U"
+read_data_name = "Current"
+write_data_name = "Voltage"
 
 # Simulation parameters and initial condition
 C = 2                      # Capacitance of the capacitor
@@ -57,10 +57,10 @@ def f_U(dt, max_allowed_dt):
     precice_dt = participant.get_max_time_step_size()
     dt = min([precice_dt, solver_dt])
     t_span = [t, t + dt]
-    sol = solve_ivp(lambda t, y: f_U(t - t_span[0], dt), t_span, U0, dense_output=True, r_tol=1e-12, a_tol=1e-12)
+    sol = solve_ivp(lambda t, y: f_U(t - t_span[0], dt), t_span, U0, dense_output=True, rtol=1e-12, atol=1e-12)
 
     # Exchange data
-    evals = 10
+    evals = max(len(sol.t), 3)  # at least do 3 substeps to allow cubic interpolation
     for i in range(evals):
         U0 = sol.sol(t_span[0] + (i + 1) * dt / evals)
         participant.write_data(mesh_name, write_data_name, vertex_ids, np.array(U0))
@@ -74,4 +74,4 @@ def f_U(dt, max_allowed_dt):
         t = t_checkpoint
 
 # Stop coupling
-participant.finalize()
+participant.finalize()
\ No newline at end of file
diff --git a/resonant-circuit/capacitor-python/run.sh b/resonant-circuit/capacitor-python/run.sh
new file mode 100755
index 000000000..f45652a4f
--- /dev/null
+++ b/resonant-circuit/capacitor-python/run.sh
@@ -0,0 +1,9 @@
+#!/usr/bin/env bash
+set -e -u
+
+. ../../tools/log.sh
+exec > >(tee --append "$LOGFILE") 2>&1
+
+python3 capacitor.py
+
+close_log
\ No newline at end of file
diff --git a/resonant-circuit/capacitor-python/SolverI.py b/resonant-circuit/coil-python/coil.py
similarity index 86%
rename from resonant-circuit/capacitor-python/SolverI.py
rename to resonant-circuit/coil-python/coil.py
index 91a534b91..f60988015 100644
--- a/resonant-circuit/capacitor-python/SolverI.py
+++ b/resonant-circuit/coil-python/coil.py
@@ -3,24 +3,23 @@
 from scipy.integrate import solve_ivp
 
 # Initialize and configure preCICE
-participant = precice.Participant("ParticipantI", "../precice-config.xml", 0, 1)
+participant = precice.Participant("Coil", "../precice-config.xml", 0, 1)
 
 # Geometry IDs. As it is a 0-D simulation, only one vertex is necessary.
-mesh_name = "MeshI"
+mesh_name = "Coil-Mesh"
 
 dimensions = participant.get_mesh_dimensions(mesh_name)
 
 vertex_ids = np.array([participant.set_mesh_vertex(mesh_name, np.zeros(dimensions))])
 
 # Data IDs
-read_data_name = "U"
-write_data_name = "I"
+read_data_name = "Voltage"
+write_data_name = "Current"
 
 # Simulation parameters and initial condition
 C = 2                      # Capacitance of the capacitor
 L = 1                      # Inductance of the coil
 t0 = 0                     # Initial simulation time
-t_max = 1                  # End simulation time
 Io = np.array([1])         # Initial current
 phi = 0                    # Phase of the signal
 
@@ -63,10 +62,10 @@ def f_I(dt, max_allowed_dt):
     precice_dt = participant.get_max_time_step_size()
     dt = min([precice_dt, solver_dt])
     t_span = [t, t + dt]
-    sol = solve_ivp(lambda t, y: f_I(t - t_span[0], dt), t_span, I0, dense_output=True, r_tol=1e-12, a_tol=1e-12)
+    sol = solve_ivp(lambda t, y: f_I(t - t_span[0], dt), t_span, I0, dense_output=True, rtol=1e-12, atol=1e-12)
 
     # Exchange data
-    evals = 10
+    evals = max(len(sol.t), 3)  # at least do 3 substeps to allow cubic interpolation
     for i in range(evals):
         I0 = sol.sol(t_span[0] + (i + 1) * dt / evals)
         participant.write_data(mesh_name, write_data_name, vertex_ids, np.array(I0))
@@ -82,12 +81,7 @@ def f_I(dt, max_allowed_dt):
 # Stop coupling
 participant.finalize()
 
-import math
-
-assert (math.isclose(t, t_max))
 def I_analytical(t): return Io * np.cos(t * w0 + phi)
-
-
-error = I0 - I_analytical(t_max)
+error = I0 - I_analytical(t)
 print(f"{error=}")
-print(f"{f_evals=}")
+print(f"{f_evals=}")
\ No newline at end of file
diff --git a/resonant-circuit/coil-python/run.sh b/resonant-circuit/coil-python/run.sh
new file mode 100755
index 000000000..cee172279
--- /dev/null
+++ b/resonant-circuit/coil-python/run.sh
@@ -0,0 +1,9 @@
+#!/usr/bin/env bash
+set -e -u
+
+. ../../tools/log.sh
+exec > >(tee --append "$LOGFILE") 2>&1
+
+python3 coil.py
+
+close_log
\ No newline at end of file
diff --git a/resonant-circuit/precice-config.xml b/resonant-circuit/precice-config.xml
index 0ffec9276..89667784a 100644
--- a/resonant-circuit/precice-config.xml
+++ b/resonant-circuit/precice-config.xml
@@ -1,7 +1,7 @@
 <?xml version="1.0" encoding="UTF-8" ?>
 <precice-configuration>
-  <data:scalar name="Current" />
-  <data:scalar name="Voltage" />
+  <data:scalar name="Current" waveform-degree="3" />
+  <data:scalar name="Voltage" waveform-degree="3" />
 
   <mesh name="Coil-Mesh" dimensions="2">
     <use-data name="Current" />
@@ -41,9 +41,11 @@
   <coupling-scheme:serial-implicit>
     <participants first="Coil" second="Capacitor" />
     <max-time value="10" />
-    <time-window-size value="0.01" />
-    <max-iterations value="3" />
-    <exchange data="Current" mesh="Coil-Mesh" from="Coil" to="Capacitor" />
-    <exchange data="Voltage" mesh="Coil-Mesh" from="Capacitor" to="Coil" />
+    <time-window-size value="1" />
+    <max-iterations value="100" />
+    <exchange data="Current" mesh="Coil-Mesh" from="Coil" to="Capacitor" substeps="true" />
+    <exchange data="Voltage" mesh="Coil-Mesh" from="Capacitor" to="Coil" substeps="true" />
+    <relative-convergence-measure limit="1e-3" data="Current" mesh="Coil-Mesh" />
+    <relative-convergence-measure limit="1e-3" data="Voltage" mesh="Coil-Mesh" />
   </coupling-scheme:serial-implicit>
 </precice-configuration>

From dfffd5eee2c9d559dfbd60058859289a2d19d5ba Mon Sep 17 00:00:00 2001
From: Benjamin Rodenberg <benjamin.rodenberg@cit.tum.de>
Date: Mon, 17 Jun 2024 13:58:57 +0200
Subject: [PATCH 08/15] Update README.md

---
 resonant-circuit/README.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/resonant-circuit/README.md b/resonant-circuit/README.md
index d139d64c0..5a421e809 100644
--- a/resonant-circuit/README.md
+++ b/resonant-circuit/README.md
@@ -64,7 +64,7 @@ By doing that, you can now open two shells and switch into the directories `capa
 
 ## Post-processing
 
-The solver for the current also records the current and voltage through time and at the end of the simulation saves a plot with the obtained curves, as well as the analytical solution.
+The MATLAB participant for the current also records the current and voltage through time and at the end of the simulation saves a plot with the obtained curves, as well as the analytical solution.
 
 After successfully running the coupling, one can find the curves in the folder `capacitor-matlab` as `Curves.png`.
 

From f1ecc7780b9bfa7334c41cdce194e5db441823b8 Mon Sep 17 00:00:00 2001
From: Benjamin Rodenberg <benjamin.rodenberg@cit.tum.de>
Date: Thu, 29 Aug 2024 10:56:13 +0200
Subject: [PATCH 09/15] Fix autopep8.

---
 resonant-circuit/capacitor-python/capacitor.py | 2 +-
 resonant-circuit/coil-python/coil.py           | 5 ++++-
 2 files changed, 5 insertions(+), 2 deletions(-)

diff --git a/resonant-circuit/capacitor-python/capacitor.py b/resonant-circuit/capacitor-python/capacitor.py
index bd2639bf7..da594f328 100644
--- a/resonant-circuit/capacitor-python/capacitor.py
+++ b/resonant-circuit/capacitor-python/capacitor.py
@@ -74,4 +74,4 @@ def f_U(dt, max_allowed_dt):
         t = t_checkpoint
 
 # Stop coupling
-participant.finalize()
\ No newline at end of file
+participant.finalize()
diff --git a/resonant-circuit/coil-python/coil.py b/resonant-circuit/coil-python/coil.py
index f60988015..ee1cbfa8e 100644
--- a/resonant-circuit/coil-python/coil.py
+++ b/resonant-circuit/coil-python/coil.py
@@ -81,7 +81,10 @@ def f_I(dt, max_allowed_dt):
 # Stop coupling
 participant.finalize()
 
+
 def I_analytical(t): return Io * np.cos(t * w0 + phi)
+
+
 error = I0 - I_analytical(t)
 print(f"{error=}")
-print(f"{f_evals=}")
\ No newline at end of file
+print(f"{f_evals=}")

From 4f1fe641c7547b172468513606856de7c7cb95db Mon Sep 17 00:00:00 2001
From: Gerasimos Chourdakis <chourdak@in.tum.de>
Date: Thu, 29 Aug 2024 14:22:45 +0200
Subject: [PATCH 10/15] Update resonant-circuit/README.md

Co-authored-by: Benjamin Rodenberg <benjamin.rodenberg@cit.tum.de>
---
 resonant-circuit/README.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/resonant-circuit/README.md b/resonant-circuit/README.md
index 79032a1aa..2b502acad 100644
--- a/resonant-circuit/README.md
+++ b/resonant-circuit/README.md
@@ -31,7 +31,7 @@ preCICE configuration (image generated using the [precice-config-visualizer](htt
 
 * *MATLAB* A solver using the [MATLAB bindings](https://precice.org/installation-bindings-matlab.html).
  Before running this tutorial, follow the [instructions](https://precice.org/installation-bindings-matlab.html) to correctly install the MATLAB bindings.
-* *Python* A solver using the preCICE [Python bindings](https://precice.org/installation-bindings-python.html). This solver also depends on the Python libraries `numpy` and `scipy`, which you can get from your system package manager or with `pip3 install --user <package>`.
+* *Python* A solver using the preCICE [Python bindings](https://precice.org/installation-bindings-python.html). 
 
 ## Running the simulation
 

From 033a2c1ea43fd774bcd4b77bd6bf4d6edba79c86 Mon Sep 17 00:00:00 2001
From: Niklas <niklas.vinnitchenko@tum.de>
Date: Thu, 29 Aug 2024 14:51:28 +0200
Subject: [PATCH 11/15] Add requirement.txt and modify run file to use venv

---
 resonant-circuit/capacitor-python/requirements.txt | 3 +++
 resonant-circuit/capacitor-python/run.sh           | 5 ++++-
 resonant-circuit/coil-python/requirements.txt      | 3 +++
 resonant-circuit/coil-python/run.sh                | 5 ++++-
 4 files changed, 14 insertions(+), 2 deletions(-)
 create mode 100644 resonant-circuit/capacitor-python/requirements.txt
 create mode 100644 resonant-circuit/coil-python/requirements.txt

diff --git a/resonant-circuit/capacitor-python/requirements.txt b/resonant-circuit/capacitor-python/requirements.txt
new file mode 100644
index 000000000..2b8edf2ba
--- /dev/null
+++ b/resonant-circuit/capacitor-python/requirements.txt
@@ -0,0 +1,3 @@
+numpy >1, <2
+scipy
+pyprecice~=3.0
diff --git a/resonant-circuit/capacitor-python/run.sh b/resonant-circuit/capacitor-python/run.sh
index f45652a4f..0864af70a 100755
--- a/resonant-circuit/capacitor-python/run.sh
+++ b/resonant-circuit/capacitor-python/run.sh
@@ -4,6 +4,9 @@ set -e -u
 . ../../tools/log.sh
 exec > >(tee --append "$LOGFILE") 2>&1
 
+python3 -m venv .venv
+. .venv/bin/activate
+pip install -r requirements.txt
 python3 capacitor.py
 
-close_log
\ No newline at end of file
+close_log
diff --git a/resonant-circuit/coil-python/requirements.txt b/resonant-circuit/coil-python/requirements.txt
new file mode 100644
index 000000000..2b8edf2ba
--- /dev/null
+++ b/resonant-circuit/coil-python/requirements.txt
@@ -0,0 +1,3 @@
+numpy >1, <2
+scipy
+pyprecice~=3.0
diff --git a/resonant-circuit/coil-python/run.sh b/resonant-circuit/coil-python/run.sh
index cee172279..3ca830a10 100755
--- a/resonant-circuit/coil-python/run.sh
+++ b/resonant-circuit/coil-python/run.sh
@@ -4,6 +4,9 @@ set -e -u
 . ../../tools/log.sh
 exec > >(tee --append "$LOGFILE") 2>&1
 
+python3 -m venv .venv
+. .venv/bin/activate
+pip install -r requirements.txt
 python3 coil.py
 
-close_log
\ No newline at end of file
+close_log

From 60c2b9264e6d22e1ac5880b1f395aa3a5c23861e Mon Sep 17 00:00:00 2001
From: Niklas <niklas.vinnitchenko@tum.de>
Date: Thu, 29 Aug 2024 15:19:09 +0200
Subject: [PATCH 12/15] Add description how to start the python case

---
 resonant-circuit/README.md | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/resonant-circuit/README.md b/resonant-circuit/README.md
index 2b502acad..0dcd22706 100644
--- a/resonant-circuit/README.md
+++ b/resonant-circuit/README.md
@@ -31,7 +31,7 @@ preCICE configuration (image generated using the [precice-config-visualizer](htt
 
 * *MATLAB* A solver using the [MATLAB bindings](https://precice.org/installation-bindings-matlab.html).
  Before running this tutorial, follow the [instructions](https://precice.org/installation-bindings-matlab.html) to correctly install the MATLAB bindings.
-* *Python* A solver using the preCICE [Python bindings](https://precice.org/installation-bindings-python.html). 
+* *Python* A solver using the preCICE [Python bindings](https://precice.org/installation-bindings-python.html).
 
 ## Running the simulation
 
@@ -62,6 +62,10 @@ For that, modify the path in the file `matlab-bindings-path.sh` found in the bas
 
 By doing that, you can now open two shells and switch into the directories `capacitor-matlab` and `coil-matlab` and execute the `run.sh` scripts.
 
+### Python
+
+This example can be executed by starting the `run.sh` file located in each of the directories `coil-python` and `capacitor-python`.
+
 ## Post-processing
 
 The MATLAB participant for the current also records the current and voltage through time and at the end of the simulation saves a plot with the obtained curves, as well as the analytical solution.

From 1b7f3207e238e6342bde1d019440490a24bb9325 Mon Sep 17 00:00:00 2001
From: Niklas <niklas.vinnitchenko@tum.de>
Date: Fri, 30 Aug 2024 11:09:32 +0200
Subject: [PATCH 13/15] Revert changes in README

---
 resonant-circuit/README.md | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/resonant-circuit/README.md b/resonant-circuit/README.md
index 0dcd22706..444e00773 100644
--- a/resonant-circuit/README.md
+++ b/resonant-circuit/README.md
@@ -62,10 +62,6 @@ For that, modify the path in the file `matlab-bindings-path.sh` found in the bas
 
 By doing that, you can now open two shells and switch into the directories `capacitor-matlab` and `coil-matlab` and execute the `run.sh` scripts.
 
-### Python
-
-This example can be executed by starting the `run.sh` file located in each of the directories `coil-python` and `capacitor-python`.
-
 ## Post-processing
 
 The MATLAB participant for the current also records the current and voltage through time and at the end of the simulation saves a plot with the obtained curves, as well as the analytical solution.

From 28d8010e000de9443bdfaaf966ba7b6b597c9e95 Mon Sep 17 00:00:00 2001
From: Benjamin Rodenberg <benjamin.rodenberg@cit.tum.de>
Date: Fri, 30 Aug 2024 12:53:40 +0200
Subject: [PATCH 14/15] Fix postprocessing.

---
 resonant-circuit/README.md          | 4 +++-
 resonant-circuit/plot-solution.sh   | 6 +++---
 resonant-circuit/precice-config.xml | 1 +
 3 files changed, 7 insertions(+), 4 deletions(-)

diff --git a/resonant-circuit/README.md b/resonant-circuit/README.md
index 0dcd22706..c98e666c5 100644
--- a/resonant-circuit/README.md
+++ b/resonant-circuit/README.md
@@ -68,7 +68,9 @@ This example can be executed by starting the `run.sh` file located in each of th
 
 ## Post-processing
 
-The MATLAB participant for the current also records the current and voltage through time and at the end of the simulation saves a plot with the obtained curves, as well as the analytical solution.
+As we defined a watchpoint on the 'Capacitor' participant (see `precice-config.xml`), we can plot it with gnuplot using the script `plot-solution.sh.` You need to specify the directory of the selected solid participant as a command line argument, so that the script can pick-up the desired watchpoint file, e.g. `./plot-solution.sh capacitor-python`. The resulting graph shows the voltage and current exchanged between the two participants.
+
+Additionally, the MATLAB participant for the current also records the current and voltage through time and at the end of the simulation saves a plot with the obtained curves, as well as the analytical solution.
 
 After successfully running the coupling, one can find the curves in the folder `capacitor-matlab` as `Curves.png`.
 
diff --git a/resonant-circuit/plot-solution.sh b/resonant-circuit/plot-solution.sh
index 30cd21a21..85f5b901c 100755
--- a/resonant-circuit/plot-solution.sh
+++ b/resonant-circuit/plot-solution.sh
@@ -5,10 +5,10 @@ if [ "${1:-}" = "" ]; then
     exit 1
 fi
 
-FILE="$1/precice-ParticipantU-watchpoint-UI.log"
+FILE="$1/precice-Capacitor-watchpoint-VoltageCurrent.log"
 
 if [ ! -f "$FILE" ]; then
-    echo "Unable to locate the watchpoint file (precice-ParticipantU-watchpoint-UI.log) in the specified participant directory '${1}'. Make sure the specified directory matches the participant you used for the calculations."
+    echo "Unable to locate the watchpoint file (precice-Capacitor-watchpoint-VoltageCurrent.log) in the specified participant directory '${1}'. Make sure the specified directory matches the participant you used for the calculations."
     exit 1
 fi
 
@@ -17,5 +17,5 @@ gnuplot -p << EOF
 	set title 'Voltage and current'
 	set xlabel 'time [s]'
 	set ylabel 'Voltage / Current'
-	plot "$1/precice-ParticipantU-watchpoint-UI.log" using 1:4 with linespoints title "Voltage", "$1/precice-ParticipantU-watchpoint-UI.log" using 1:5 with linespoints title "Current"
+	plot "$1/precice-Capacitor-watchpoint-VoltageCurrent.log" using 1:4 with linespoints title "Voltage", "$1/precice-Capacitor-watchpoint-VoltageCurrent.log" using 1:5 with linespoints title "Current"
 EOF
diff --git a/resonant-circuit/precice-config.xml b/resonant-circuit/precice-config.xml
index 89667784a..a4bd640eb 100644
--- a/resonant-circuit/precice-config.xml
+++ b/resonant-circuit/precice-config.xml
@@ -34,6 +34,7 @@
       from="Coil-Mesh"
       to="Capacitor-Mesh"
       constraint="consistent" />
+    <watch-point mesh="Capacitor-Mesh" name="VoltageCurrent" coordinate="0.0;0.0" />
   </participant>
 
   <m2n:sockets acceptor="Coil" connector="Capacitor" exchange-directory=".." />

From 64fdfa1a3c5f55f62f7e7904c12b1c0c6f647359 Mon Sep 17 00:00:00 2001
From: Benjamin Rodenberg <benjamin.rodenberg@cit.tum.de>
Date: Fri, 30 Aug 2024 12:57:41 +0200
Subject: [PATCH 15/15] Update resonant-circuit/README.md

---
 resonant-circuit/README.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/resonant-circuit/README.md b/resonant-circuit/README.md
index 9ddaf99dc..7af7c7c78 100644
--- a/resonant-circuit/README.md
+++ b/resonant-circuit/README.md
@@ -66,7 +66,7 @@ By doing that, you can now open two shells and switch into the directories `capa
 
 As we defined a watchpoint on the 'Capacitor' participant (see `precice-config.xml`), we can plot it with gnuplot using the script `plot-solution.sh.` You need to specify the directory of the selected solid participant as a command line argument, so that the script can pick-up the desired watchpoint file, e.g. `./plot-solution.sh capacitor-python`. The resulting graph shows the voltage and current exchanged between the two participants.
 
-Additionally, the MATLAB participant for the current also records the current and voltage through time and at the end of the simulation saves a plot with the obtained curves, as well as the analytical solution.
+Additionally, the MATLAB participant `capacitor-matlab` records the current and voltage over time. At the end of the simulation it creates a plot with the computed waveforms of current and voltage, as well as the analytical solution.
 
 After successfully running the coupling, one can find the curves in the folder `capacitor-matlab` as `Curves.png`.