diff --git a/.gitignore b/.gitignore index 7f2055b46..4b7dafe69 100644 --- a/.gitignore +++ b/.gitignore @@ -40,6 +40,10 @@ __pycache__/ Cargo.lock target/ +# Julia +Manifest.toml +Project.toml + # OpenFOAM 0.*/ [1-9]*/ diff --git a/changelog-entries/658.md b/changelog-entries/658.md new file mode 100644 index 000000000..cf9e728c7 --- /dev/null +++ b/changelog-entries/658.md @@ -0,0 +1 @@ +- Added Julia-based participants to the resonant-circuit tutorial [#658](https://github.com/precice/tutorials/pull/658) \ No newline at end of file diff --git a/resonant-circuit/README.md b/resonant-circuit/README.md index fbbcbc4e9..4f81fa45b 100644 --- a/resonant-circuit/README.md +++ b/resonant-circuit/README.md @@ -1,14 +1,14 @@ --- -title: Resonant Circuit +title: Resonant circuit permalink: tutorials-resonant-circuit.html -keywords: MATLAB, Python +keywords: MATLAB, Python, Julia summary: We simulate a two-element LC circuit (one inductor and one capacitor). --- ## Setup -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 of a very simple system with one inductor and one capacitor: +Two different solvers are coupled to simulate a two-element LC circuit. This type of circuit consists of a simple system with one inductor and one capacitor: ![LC circuit diagram [1]](images/tutorials-resonant-circuit-diagram.svg) @@ -20,7 +20,7 @@ $I(t) = -C \frac{\text{d}U}{\text{d}t}$ where $I$ is the current and $U$ the voltage of the circuit. -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. +Each of these equations is solved by a different solver. Note that, as only one scalar is solved per equation, this is a 0+1 dimensional problem. ## Configuration @@ -31,43 +31,42 @@ preCICE configuration (image generated using the [precice-config-visualizer](htt ## Available solvers * *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). +* *Julia* A solver using the preCICE [Julia bindings](https://precice.org/installation-bindings-julia.html). ## Running the simulation -### MATLAB +All listed solvers can be used to run the simulation. Open two separate terminals and start the desired capacitor and coil participants by calling the respective run script. For example: -For running this example, first get into one of the solver folders and open a MATLAB instance. -Afterward, do the same for the second solver. -After adding the MATLAB bindings to the MATLAB path (in both instances), run the following commands: - -In the first MATLAB instance, one can run the solver for the current: - -```MATLAB -coil +```bash +cd capacitor-python +./run.sh ``` -And in the second MATLAB instance, the solver for the voltage: +and -```MATLAB -capacitor +```bash +cd coil-julia +./run.sh ``` -The preCICE configuration file is hard-coded as `precice-config.xml` in the solvers. +### Running the MATLAB participants + +For running this example in the MATLAB GUI, first get into each of the solver folders and open a MATLAB instance for each. +After adding the MATLAB bindings to the MATLAB path (in both instances), run the `coil` and `capacitor` commands in the two windows. -#### Running from terminal +The path to the preCICE configuration file is hard-coded as `precice-config.xml` in the solvers. -If you prefer to not open the MATLAB GUIs, you can alternatively use two shells instead. +If you prefer not to use the MATLAB GUI, you can alternatively use two shells instead. For that, modify the path in the file `matlab-bindings-path.sh` found in the base directory of this tutorial to the path to your MATLAB bindings. 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. ## Post-processing -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. +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. +Additionally, the `capacitor-matlab` case 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-julia/capacitor.jl b/resonant-circuit/capacitor-julia/capacitor.jl new file mode 100644 index 000000000..9219c5cae --- /dev/null +++ b/resonant-circuit/capacitor-julia/capacitor.jl @@ -0,0 +1,82 @@ +using PreCICE +using DifferentialEquations + +# Initialize and configure preCICE +participant = PreCICE.createParticipant("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 = PreCICE.getMeshDimensions(mesh_name) + +vertex_ids = PreCICE.setMeshVertices(mesh_name, zeros((1,dimensions))) + +let + # 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 = 1 # Initial current + phi = 0 # Phase of the signal + + w0 = 1 / sqrt(L * C) # Resonant frequency + U0 = -w0 * L * Io * sin(phi) # Initial condition for U + + function f_U(u, p, t) + (dt, C, mesh_name, read_data_name, vertex_ids) = p + I = PreCICE.readData(mesh_name, read_data_name, vertex_ids, dt) + return -I[1] / C # Time derivative of U + end + + # Initialize simulation + if PreCICE.requiresInitialData() + PreCICE.writeData(mesh_name, write_data_name, vertex_ids, I0) + end + PreCICE.initialize() + + solver_dt = PreCICE.getMaxTimeStepSize() + + # Start simulation + t = t0 + U0_checkpoint = U0 + t_checkpoint = t + while PreCICE.isCouplingOngoing() + + # Record checkpoint if necessary + if PreCICE.requiresWritingCheckpoint() + U0_checkpoint = U0 + t_checkpoint = t + end + + # Make Simulation Step + precice_dt = PreCICE.getMaxTimeStepSize() + dt = min(precice_dt, solver_dt) + t_span = (t, t + dt) + params = (dt, C, mesh_name, read_data_name, vertex_ids) + prob = ODEProblem(f_U, U0, t_span, params) + # https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts + sol = solve(prob, Tsit5(), reltol=1e-12, abstol=1e-12) + + # Exchange data + evals = max(length(sol.t), 3) # at least do 3 substeps to allow cubic interpolation + for i in range(1,evals) + U0 = sol(t_span[1] + i * dt / evals) + PreCICE.writeData(mesh_name, write_data_name, vertex_ids, fill(U0, (1,1))) + PreCICE.advance(dt / evals) + end + t = t + dt + + # Recover checkpoint if not converged + if PreCICE.requiresReadingCheckpoint() + U0 = U0_checkpoint + t = t_checkpoint + end + end + # Stop coupling + PreCICE.finalize() +end \ No newline at end of file diff --git a/resonant-circuit/capacitor-julia/clean.sh b/resonant-circuit/capacitor-julia/clean.sh new file mode 100755 index 000000000..8c3dba08a --- /dev/null +++ b/resonant-circuit/capacitor-julia/clean.sh @@ -0,0 +1,7 @@ +#!/bin/sh +set -e -u + +. ../../tools/cleaning-tools.sh + +clean_precice_logs . +clean_case_logs . \ No newline at end of file diff --git a/resonant-circuit/capacitor-julia/run.sh b/resonant-circuit/capacitor-julia/run.sh new file mode 100755 index 000000000..a61d1fa67 --- /dev/null +++ b/resonant-circuit/capacitor-julia/run.sh @@ -0,0 +1,10 @@ +#!/usr/bin/env bash +set -e -u + +. ../../tools/log.sh +exec > >(tee --append "$LOGFILE") 2>&1 + +julia --project=Project.toml -e 'using Pkg; Pkg.add("DifferentialEquations"); Pkg.add("PreCICE"); Pkg.instantiate();' +julia --project=Project.toml capacitor.jl + +close_log \ No newline at end of file diff --git a/resonant-circuit/capacitor-matlab/clean.sh b/resonant-circuit/capacitor-matlab/clean.sh index 6a77c6d6c..8b1311d03 100755 --- a/resonant-circuit/capacitor-matlab/clean.sh +++ b/resonant-circuit/capacitor-matlab/clean.sh @@ -1,8 +1,7 @@ #!/bin/sh set -e -u -rm ./*.png ./*.mat - . ../../tools/cleaning-tools.sh clean_matlab . +rm -f ./*.png ./*.mat diff --git a/resonant-circuit/capacitor-python/clean.sh b/resonant-circuit/capacitor-python/clean.sh new file mode 100755 index 000000000..8c3dba08a --- /dev/null +++ b/resonant-circuit/capacitor-python/clean.sh @@ -0,0 +1,7 @@ +#!/bin/sh +set -e -u + +. ../../tools/cleaning-tools.sh + +clean_precice_logs . +clean_case_logs . \ No newline at end of file diff --git a/resonant-circuit/coil-julia/clean.sh b/resonant-circuit/coil-julia/clean.sh new file mode 100755 index 000000000..8c3dba08a --- /dev/null +++ b/resonant-circuit/coil-julia/clean.sh @@ -0,0 +1,7 @@ +#!/bin/sh +set -e -u + +. ../../tools/cleaning-tools.sh + +clean_precice_logs . +clean_case_logs . \ No newline at end of file diff --git a/resonant-circuit/coil-julia/coil.jl b/resonant-circuit/coil-julia/coil.jl new file mode 100644 index 000000000..6cd949ebb --- /dev/null +++ b/resonant-circuit/coil-julia/coil.jl @@ -0,0 +1,85 @@ +using PreCICE +using DifferentialEquations + +# Initialize and configure preCICE +participant = PreCICE.createParticipant("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 = PreCICE.getMeshDimensions(mesh_name) + +vertex_ids = PreCICE.setMeshVertices(mesh_name, zeros((1,dimensions))) + +let + # 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 = 1 # Initial current + phi = 0 # Phase of the signal + + w0 = 1 / sqrt(L * C) # Resonant frequency + I0 = Io * cos(phi) # Initial condition for I + + # to estimate cost + f_evals = 0 + + function f_I(u, p, t) + f_evals += 1 + (dt, L, mesh_name, read_data_name, vertex_ids) = p + U = PreCICE.readData(mesh_name, read_data_name, vertex_ids, dt) + return U[1] / L # Time derivative of I; ODE determining capacitor + end + + # Initialize simulation + if PreCICE.requiresInitialData() + PreCICE.writeData(mesh_name, write_data_name, vertex_ids, I0) + end + PreCICE.initialize() + + solver_dt = PreCICE.getMaxTimeStepSize() + + # Start simulation + t = t0 + I0_checkpoint = I0 + t_checkpoint = t + while PreCICE.isCouplingOngoing() + + # Record checkpoint if necessary + if PreCICE.requiresWritingCheckpoint() + I0_checkpoint = I0 + t_checkpoint = t + end + + # Make Simulation Step + precice_dt = PreCICE.getMaxTimeStepSize() + dt = min(precice_dt, solver_dt) + t_span = (t, t + dt) + params = (dt, L, mesh_name, read_data_name, vertex_ids) + prob = ODEProblem(f_I, I0, t_span, params) + # https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts + sol = solve(prob, Tsit5(), reltol=1e-12, abstol=1e-12) + + # Exchange data + evals = max(length(sol.t), 3) # at least do 3 substeps to allow cubic interpolation + for i in range(1,evals) + I0 = sol(t_span[1] + i * dt / evals) + PreCICE.writeData(mesh_name, write_data_name, vertex_ids, fill(I0, (1,1))) + PreCICE.advance(dt / evals) + end + t = t + dt + + # Recover checkpoint if not converged + if PreCICE.requiresReadingCheckpoint() + I0 = I0_checkpoint + t = t_checkpoint + end + end + # Stop coupling + PreCICE.finalize() +end \ No newline at end of file diff --git a/resonant-circuit/coil-julia/run.sh b/resonant-circuit/coil-julia/run.sh new file mode 100755 index 000000000..6340410ae --- /dev/null +++ b/resonant-circuit/coil-julia/run.sh @@ -0,0 +1,10 @@ +#!/usr/bin/env bash +set -e -u + +. ../../tools/log.sh +exec > >(tee --append "$LOGFILE") 2>&1 + +julia --project=Project.toml -e 'using Pkg; Pkg.add("DifferentialEquations"); Pkg.add("PreCICE"); Pkg.instantiate();' +julia --project=Project.toml coil.jl + +close_log \ No newline at end of file diff --git a/resonant-circuit/coil-python/clean.sh b/resonant-circuit/coil-python/clean.sh new file mode 100755 index 000000000..8c3dba08a --- /dev/null +++ b/resonant-circuit/coil-python/clean.sh @@ -0,0 +1,7 @@ +#!/bin/sh +set -e -u + +. ../../tools/cleaning-tools.sh + +clean_precice_logs . +clean_case_logs . \ No newline at end of file diff --git a/resonant-circuit/precice-config.xml b/resonant-circuit/precice-config.xml index cda70ee19..fb72e5772 100644 --- a/resonant-circuit/precice-config.xml +++ b/resonant-circuit/precice-config.xml @@ -1,5 +1,12 @@ + + + +