diff --git a/resonant-circuit/README.md b/resonant-circuit/README.md index 13fb11bec..7af7c7c78 100644 --- 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 new file mode 100644 index 000000000..da594f328 --- /dev/null +++ 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 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 new file mode 100755 index 000000000..0864af70a --- /dev/null +++ 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 new file mode 100644 index 000000000..ee1cbfa8e --- /dev/null +++ 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 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 new file mode 100755 index 000000000..3ca830a10 --- /dev/null +++ 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 new file mode 100755 index 000000000..85f5b901c --- /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-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 index 0ffec9276..a4bd640eb 100644 --- a/resonant-circuit/precice-config.xml +++ b/resonant-circuit/precice-config.xml @@ -1,7 +1,7 @@ - - + + @@ -34,6 +34,7 @@ from="Coil-Mesh" to="Capacitor-Mesh" constraint="consistent" /> + @@ -41,9 +42,11 @@ - - - - + + + + + +