oscillator-overlap: Make consistent with oscillator (#597)

---------

Co-authored-by: Benjamin Uekermann <[email protected]>

changelog-entries/597.md

@@ -0,0 +1 @@
+- Added Runge-Kutta 4 scheme and black-box time-stepping with RadauIIA scheme as options to oscillator-overlap case. This demonstrates use of time interpolation and dense output.

oscillator-overlap/precice-config.xml

@@ -54,12 +54,14 @@
       mesh="Mass-Left-Mesh"
       from="Mass-Left"
       to="Mass-Right"
-      initialize="true" />
+      initialize="true"
+      substeps="true" />
     <exchange
       data="Displacement-Right"
       mesh="Mass-Left-Mesh"
       from="Mass-Right"
       to="Mass-Left"
-      initialize="true" />
+      initialize="true"
+      substeps="true" />
   </coupling-scheme:serial-implicit>
 </precice-configuration>

oscillator-overlap/solver-python/oscillator.py

@@ -6,13 +6,10 @@
 from enum import Enum
 import csv
 import os
+from typing import Type
 
 import problemDefinition
-
-
-class Scheme(Enum):
-    NEWMARK_BETA = "Newmark_beta"
-    GENERALIZED_ALPHA = "generalized_alpha"
+from timeSteppers import TimeStepper, TimeSteppingSchemes, GeneralizedAlpha, RungeKutta4, RadauIIA
 
 
 class Participant(Enum):
@@ -22,20 +19,30 @@ class Participant(Enum):
 
 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)
+parser.add_argument(
+    "-ts",
+    "--time-stepping",
+    help="Time stepping scheme being used.",
+    type=str,
+    choices=[
+        s.value for s in TimeSteppingSchemes],
+    default=TimeSteppingSchemes.NEWMARK_BETA.value)
 args = parser.parse_args()
 
 participant_name = args.participantName
 
+this_mass: Type[problemDefinition.Mass]
+other_mass: Type[problemDefinition.Mass]
+this_spring: Type[problemDefinition.Spring]
+connecting_spring = problemDefinition.SpringMiddle
+
 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:
@@ -45,7 +52,6 @@ class Participant(Enum):
 
     this_mass = problemDefinition.MassRight
     this_spring = problemDefinition.SpringRight
-    connecting_spring = problemDefinition.SpringMiddle
     other_mass = problemDefinition.MassLeft
 
 else:
@@ -68,7 +74,7 @@ class Participant(Enum):
 
 vertex = np.zeros(dimensions)
 read_data = np.zeros(num_vertices)
-write_data = connecting_spring.k * u0 * np.ones(num_vertices)
+write_data = u0 * np.ones(num_vertices)
 
 vertex_ids = [participant.set_mesh_vertex(mesh_name, vertex)]
 
@@ -77,7 +83,7 @@ class Participant(Enum):
 
 participant.initialize()
 precice_dt = participant.get_max_time_step_size()
-my_dt = precice_dt  # 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
@@ -86,16 +92,20 @@ class Participant(Enum):
 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)
+time_stepper: TimeStepper
+
+if args.time_stepping == TimeSteppingSchemes.GENERALIZED_ALPHA.value:
+    time_stepper = GeneralizedAlpha(stiffness=stiffness, mass=mass, alpha_f=0.4, alpha_m=0.2)
+elif args.time_stepping == TimeSteppingSchemes.NEWMARK_BETA.value:
+    time_stepper = GeneralizedAlpha(stiffness=stiffness, mass=mass, alpha_f=0.0, alpha_m=0.0)
+elif args.time_stepping == TimeSteppingSchemes.RUNGE_KUTTA_4.value:
+    time_stepper = RungeKutta4(stiffness=stiffness, mass=mass)
+elif args.time_stepping == TimeSteppingSchemes.Radau_IIA.value:
+    time_stepper = RadauIIA(stiffness=stiffness, mass=mass)
+else:
+    raise Exception(
+        f"Invalid time stepping scheme {args.time_stepping}. Please use one of {[ts.value for ts in TimeSteppingSchemes]}")
+
 
 positions = []
 velocities = []
@@ -119,25 +129,36 @@ class Participant(Enum):
     # 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]
+    def f(t: float) -> float: return connecting_spring.k * \
+        participant.read_data(mesh_name, read_data_name, vertex_ids, t)[0]
 
-    participant.write_data(mesh_name, write_data_name, vertex_ids, write_data)
+    # do time step, write data, and advance
+    u_new, v_new, a_new = time_stepper.do_step(u, v, a, f, dt)
+
+    t_new = t + dt
 
-    participant.advance(dt)
+    # RadauIIA time stepper provides dense output. Do multiple write calls per time step.
+    if isinstance(time_stepper, RadauIIA):
+        # create n samples_per_step of time stepping scheme. Degree of dense
+        # interpolating function is usually larger 1 and, therefore, we need
+        # multiple samples per step.
+        samples_per_step = 5
+        n_time_steps = len(time_stepper.dense_output.ts)  # number of time steps performed by adaptive time stepper
+        n_pseudo = samples_per_step * n_time_steps  # samples_per_step times no. time steps per window.
+        t_pseudo = 0
+        for i in range(n_pseudo):
+            # perform n_pseudo pseudosteps
+            dt_pseudo = dt / n_pseudo
+            t_pseudo += dt_pseudo
+            write_data = np.array([time_stepper.dense_output(t_pseudo)[0]])
+            participant.write_data(mesh_name, write_data_name, vertex_ids, write_data)
+            participant.advance(dt_pseudo)
+
+    else:  # simple time stepping without dense output; only a single write call per time step
+        write_data = np.array([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
@@ -161,12 +182,6 @@ class Participant(Enum):
         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

oscillator-overlap/solver-python/problemDefinition.py

@@ -1,20 +1,33 @@
 import numpy as np
 from numpy.linalg import eig
+from typing import Callable
 
 
-class SpringLeft:
+class Spring:
+    k: float
+
+
+class SpringLeft(Spring):
     k = 4 * np.pi**2
 
 
-class SpringMiddle:
+class SpringMiddle(Spring):
     k = 16 * (np.pi**2)
 
 
-class SpringRight:
+class SpringRight(Spring):
     k = 4 * np.pi**2
 
 
-class MassLeft:
+class Mass:
+    m: float
+    u0: float
+    v0: float
+    u_analytical: Callable[[float | np.ndarray], float | np.ndarray]
+    v_analytical: Callable[[float | np.ndarray], float | np.ndarray]
+
+
+class MassLeft(Mass):
     # mass
     m = 1
 
@@ -23,10 +36,8 @@ class MassLeft:
     u0 = 1.0
     v0 = 0.0
 
-    u_analytical, v_analytical = None, None  # will be defined below
-
 
-class MassRight:
+class MassRight(Mass):
     # mass
     m = 1
 
@@ -35,8 +46,6 @@ class MassRight:
     u0 = 0.0
     v0 = 0.0
 
-    u_analytical, v_analytical = None, None  # will be defined below
-
 
 # Mass matrix
 M = np.array([