Add FMU participants to oscillator tutorial (#466)

Co-authored-by: Leonard Willeke <[email protected]>
Co-authored-by: LeonardWilleke <[email protected]>
Co-authored-by: Gerasimos Chourdakis <[email protected]>

oscillator/.gitignore

@@ -0,0 +1 @@
+fmi/Oscillator.fmu

oscillator/README.md

@@ -1,7 +1,7 @@
 ---
 title: Oscillator
 permalink: tutorials-oscillator.html
-keywords: Python, ODE
+keywords: Python, ODE, FMI
 summary: We solve an oscillator with two masses in a partitioned fashion. Each mass is solved by an independent ODE.
 ---
 
@@ -17,19 +17,16 @@ This tutorial solves a simple mass-spring oscillator with two masses and three s
 
 Note that this case applies a 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].
 
-This tutorial is useful to study different time stepping and coupling schemes. It uses subcycling and time interpolation: Each solver performs 4 time steps in each time window. The data of these 4 substeps is used to create a third order B-spline interpolation (`waveform-degree="3"` in `precice-config.xml`).
-
 ## 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>`. Using the option `-ts` allows to pick the time stepping scheme being used. Available choices are Newmark beta, generalized alpha, explicit Runge Kutta 4, and implicit RadauIIA.
+There are two different implementations:
 
-## Running the Simulation
+- *Python*: A 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>`. Using the option `-ts` allows you to pick the time stepping scheme being used. Available choices are Newmark beta, generalized alpha, explicit Runge Kutta 4, and implicit RadauIIA. The solver uses subcycling: Each participant performs 4 time steps in each time window. The data of these 4 substeps is then used by preCICE to create a third order B-spline interpolation (`waveform-degree="3"` in `precice-config.xml`).
+- *FMI*: A solver using the [preCICE-FMI runner](https://github.com/precice/fmi-runner) (requires at least v0.2). The Runner executes the FMU model `Oscillator.fmu` for computation. The provided run script (see below) builds this model if not already there. If you want to change the model parameters or the initial conditions of the simulation, have a look inside the setting files for [MassLeft](https://github.com/precice/tutorials/tree/master/oscillator/fmi/MassLeft) and [MassRight](https://github.com/precice/tutorials/tree/master/oscillator/fmi/MassRight). For more information, please refer to [2].
 
-### Python
+## Running the simulation
 
-Open two separate terminals and start both participants by calling:
+Open two separate terminals and start both participants. For example, you can run a simulation where the left participant is computed in Python and the right participant is computed with FMI with these commands:
 
 ```bash
 cd python
@@ -39,19 +36,21 @@ cd python
 and
 
 ```bash
-cd python
+cd fmi
 ./run.sh -r
 ```
 
+Of course, you can also use the same solver for both sides.
+
 ## Post-processing
 
-Each simulation run creates two files containing position and velocity of the two masses over time. These files are called `trajectory-Mass-Left.csv` and `trajectory-Mass-Right.csv`. You can use the script `plot-trajectory.py` for post-processing. Type `python3 plot-trajectory --help` to see available options. You can, for example plot the trajectory by running
+Each simulation run creates two files containing position and velocity of the two masses over time. These files are called `trajectory-Mass-Left.csv` and `trajectory-Mass-Right.csv`. You can use the script `plot-trajectory.py` for post-processing. Type `python3 plot-trajectory --help` to see available options. You can, for example, plot the trajectory of the Python solver by running
 
 ```bash
 python3 plot-trajectory.py python/output/trajectory-Mass-Left.csv TRAJECTORY
 ```
 
-This allows you to study the effect of different time stepping schemes on energy conservation. Newmark beta conserves energy:
+The solvers allow you to study the effect of different time stepping schemes on energy conservation. Newmark beta conserves energy:
 
 ![Trajectory for Newmark beta scheme](images/tutorials-oscillator-trajectory-newmark-beta.png)
 
@@ -64,3 +63,5 @@ For details, refer to [1].
 ## References
 
 [1] V. Schüller, B. Rodenberg, B. Uekermann and H. Bungartz, A Simple Test Case for Convergence Order in Time and Energy Conservation of Black-Box Coupling Schemes, in: WCCM-APCOM2022. [URL](https://www.scipedia.com/public/Rodenberg_2022a)
+
+[2] L. Willeke, D. Schneider and B. Uekermann, A preCICE-FMI Runner to Couple FMUs to PDE-Based Simulations, Proceedings 15th Intern. Modelica Conference, 2023. [DOI](https://doi.org/10.3384/ecp204)

oscillator/fmi/MassLeft/fmi-settings.json

@@ -0,0 +1,20 @@
+{
+    "model_params": {
+        "mass.m": 1,
+        "spring_fixed.c": 39.4784176044,
+        "spring_middle.c": 157.913670417
+    },
+    "initial_conditions": {
+        "mass.u": 1,
+        "mass.v": 0, 
+        "mass.a": -197.392088022
+    },
+    "simulation_params": {
+        "fmu_file_name": "./Oscillator.fmu",
+        "output_file_name": "./output/trajectory-Mass-Left.csv",
+        "output": ["mass.u", "mass.v"],
+        "fmu_read_data_names": ["force_in"],
+        "fmu_write_data_names": ["force_out"],
+        "fmu_instance_name": "instance_left"
+    }
+}

oscillator/fmi/MassLeft/precice-settings.json

@@ -0,0 +1,9 @@
+{
+    "coupling_params": {
+        "participant_name": "Mass-Left",
+        "config_file_name": "../precice-config.xml",
+        "mesh_name": "Mass-Left-Mesh",
+        "write_data_name": "Force-Left",
+        "read_data_name": "Force-Right"
+    }
+}

oscillator/fmi/MassRight/fmi-settings.json

@@ -0,0 +1,20 @@
+{
+    "model_params": {
+        "mass.m": 1,
+        "spring_fixed.c": 39.4784176044,
+        "spring_middle.c": 157.913670417
+    },
+    "initial_conditions": {
+        "mass.u": 0,
+        "mass.v": 0, 
+        "mass.a": 157.913670417
+    },
+    "simulation_params": {
+        "fmu_file_name": "./Oscillator.fmu",
+        "output_file_name": "./output/trajectory-Mass-Right.csv",
+        "output": ["mass.u", "mass.v"],
+        "fmu_read_data_names": ["force_in"],
+        "fmu_write_data_names": ["force_out"],
+        "fmu_instance_name": "instance_right"
+    }
+}

oscillator/fmi/MassRight/precice-settings.json

@@ -0,0 +1,9 @@
+{
+    "coupling_params": {
+        "participant_name": "Mass-Right",
+        "config_file_name": "../precice-config.xml",
+        "mesh_name": "Mass-Right-Mesh",
+        "write_data_name": "Force-Right",
+        "read_data_name": "Force-Left"
+    }
+}

oscillator/fmi/calculate-error.py

@@ -0,0 +1,93 @@
+import numpy as np
+import pandas as pd
+import json
+import os
+import argparse
+from enum import Enum
+import sys
+
+
+class Participant(Enum):
+    MASS_LEFT = "Mass-Left"
+    MASS_RIGHT = "Mass-Right"
+
+
+parser = argparse.ArgumentParser()
+parser.add_argument("fmi_setting_file_left", help="Path to the fmi setting file for MassLeft.", type=str)
+parser.add_argument("precice_setting_file_left", help="Path to the precice setting file for MassLeft.", type=str)
+parser.add_argument("fmi_setting_file_right", help="Path to the fmi setting file for MassRight.", type=str)
+parser.add_argument("precice_setting_file_right", help="Path to the precice setting file for MassRight.", type=str)
+parser.add_argument("participant_name", help="Participant for which the error should be calculated", type=str,
+                    choices=[p.value for p in Participant])
+args = parser.parse_args()
+
+# Get input files
+fmi_file_left = args.fmi_setting_file_left
+precice_file_left = args.precice_setting_file_left
+fmi_file_right = args.fmi_setting_file_right
+precice_file_right = args.precice_setting_file_right
+participant_name = args.participant_name
+
+# Read json files
+folder = os.path.dirname(os.path.join(os.getcwd(), os.path.dirname(sys.argv[0]), fmi_file_left))
+path = os.path.join(folder, os.path.basename(fmi_file_left))
+read_file = open(path, "r")
+fmi_data_left = json.load(read_file)
+
+folder = os.path.dirname(os.path.join(os.getcwd(), os.path.dirname(sys.argv[0]), precice_file_left))
+path = os.path.join(folder, os.path.basename(precice_file_left))
+read_file = open(path, "r")
+precice_data_left = json.load(read_file)
+
+folder = os.path.dirname(os.path.join(os.getcwd(), os.path.dirname(sys.argv[0]), fmi_file_right))
+path = os.path.join(folder, os.path.basename(fmi_file_right))
+read_file = open(path, "r")
+fmi_data_right = json.load(read_file)
+
+folder = os.path.dirname(os.path.join(os.getcwd(), os.path.dirname(sys.argv[0]), precice_file_right))
+path = os.path.join(folder, os.path.basename(precice_file_right))
+read_file = open(path, "r")
+precice_data_right = json.load(read_file)
+
+
+# Define variables
+k_1 = fmi_data_left["model_params"]["spring_fixed.c"]
+k_2 = fmi_data_right["model_params"]["spring_fixed.c"]
+u0_1 = fmi_data_left["initial_conditions"]["mass.u"]
+u0_2 = fmi_data_right["initial_conditions"]["mass.u"]
+k_12_left = fmi_data_left["model_params"]["spring_middle.c"]
+k_12_right = fmi_data_right["model_params"]["spring_middle.c"]
+
+if k_12_left == k_12_right:
+    k_12 = k_12_left
+else:
+    raise Exception("k_12 has to be equal in both participants. Please adjust input values.")
+
+
+# Define analytical solution and read computed results
+K = np.array([[k_1 + k_12, -k_12], [-k_12, k_2 + k_12]])
+eigenvalues, eigenvectors = np.linalg.eig(K)
+omega = np.sqrt(eigenvalues)
+A, B = eigenvectors
+c = np.linalg.solve(eigenvectors, [u0_1, u0_2])
+
+if participant_name == Participant.MASS_LEFT.value:
+    filename = fmi_data_left["simulation_params"]["output_file_name"]
+    df = pd.read_csv(filename, delimiter=',')
+
+    def u_analytical(t): return c[0] * A[0] * np.cos(omega[0] * t) + c[1] * A[1] * np.cos(omega[1] * t)
+
+elif participant_name == Participant.MASS_RIGHT.value:
+    filename = fmi_data_right["simulation_params"]["output_file_name"]
+    df = pd.read_csv(filename, delimiter=',')
+
+    def u_analytical(t): return c[0] * B[0] * np.cos(omega[0] * t) + c[1] * B[1] * np.cos(omega[1] * t)
+
+times = df.iloc[:, 0]
+positions = df.iloc[:, 1]
+
+
+# Calculate error
+error = np.max(abs(u_analytical(np.array(times)) - np.array(positions)))
+print("Error w.r.t analytical solution:")
+print(f"{error}")

oscillator/fmi/clean.sh

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