Add python variant of resonant circuit (#498)

* Add python version of resonator. * Add of python case in README.md * Add requirement.txt and modify run file to use venv --------- Co-authored-by: NiklasV <[email protected]> Co-authored-by: Gerasimos Chourdakis <[email protected]> Co-authored-by: Niklas <[email protected]> Co-authored-by: Benjamin Rodenberg <[email protected]>
precice · Aug 30, 2024 · b80509a · b80509a
1 parent c4eb0ec
commit b80509a
Show file tree

Hide file tree

Showing 9 changed files with 232 additions and 8 deletions.
diff --git a/resonant-circuit/README.md b/resonant-circuit/README.md
@@ -1,6 +1,6 @@
 ---
 title: Resonant Circuit
-keywords: MATLAB
+keywords: MATLAB, Python
 summary: We simulate a two-element LC circuit (one inductor and one capacitor).
 ---
 
@@ -31,6 +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).
 
 ## Running the simulation
 
@@ -63,7 +64,9 @@ 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.
+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 `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`.
 

diff --git a/resonant-circuit/capacitor-python/capacitor.py b/resonant-circuit/capacitor-python/capacitor.py
@@ -0,0 +1,77 @@
+import precice
+import numpy as np
+from scipy.integrate import solve_ivp
+
+# Initialize and configure preCICE
+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 = "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 = "Current"
+write_data_name = "Voltage"
+
+# 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 = 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
+
+
+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
+
+
+# Initialize simulation
+if participant.requires_initial_data():
+    participant.write_data(mesh_name, write_data_name, vertex_ids, U0)
+
+participant.initialize()
+
+solver_dt = participant.get_max_time_step_size()
+
+# Start simulation
+t = t0
+while participant.is_coupling_ongoing():
+
+    # Record checkpoint if necessary
+    if participant.requires_writing_checkpoint():
+        U0_checkpoint = U0
+        t_checkpoint = t
+
+    # 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, rtol=1e-12, atol=1e-12)
+
+    # Exchange data
+    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))
+        participant.advance(dt / evals)
+
+    t = t + dt
+
+    # Recover checkpoint if not converged
+    if participant.requires_reading_checkpoint():
+        U0 = U0_checkpoint
+        t = t_checkpoint
+
+# Stop coupling
+participant.finalize()
diff --git a/resonant-circuit/capacitor-python/requirements.txt 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
@@ -0,0 +1,12 @@
+#!/usr/bin/env bash
+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
diff --git a/resonant-circuit/coil-python/coil.py b/resonant-circuit/coil-python/coil.py
@@ -0,0 +1,90 @@
+import precice
+import numpy as np
+from scipy.integrate import solve_ivp
+
+# Initialize and configure preCICE
+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 = "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 = "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
+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
+
+# 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
+
+
+# Initialize simulation
+if participant.requires_initial_data():
+    participant.write_data(mesh_name, write_data_name, vertex_ids, I0)
+
+participant.initialize()
+
+solver_dt = participant.get_max_time_step_size()
+
+# Start simulation
+t = t0
+while participant.is_coupling_ongoing():
+
+    # Record checkpoint if necessary
+    if participant.requires_writing_checkpoint():
+        I0_checkpoint = I0
+        t_checkpoint = t
+
+    # 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, rtol=1e-12, atol=1e-12)
+
+    # Exchange data
+    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))
+        participant.advance(dt / evals)
+
+    t = t + dt
+
+    # Recover checkpoint if not converged
+    if participant.requires_reading_checkpoint():
+        I0 = I0_checkpoint
+        t = t_checkpoint
+
+# 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=}")
diff --git a/resonant-circuit/coil-python/requirements.txt 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
@@ -0,0 +1,12 @@
+#!/usr/bin/env bash
+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
diff --git a/resonant-circuit/plot-solution.sh 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-Capacitor-watchpoint-VoltageCurrent.log"
+
+if [ ! -f "$FILE" ]; then
+    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
+
+gnuplot -p << EOF
+	set grid
+	set title 'Voltage and current'
+	set xlabel 'time [s]'
+	set ylabel 'Voltage / 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
@@ -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" />
@@ -34,16 +34,19 @@
       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=".." />
 
   <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>