diff --git a/resonant-circuit/README.md b/resonant-circuit/README.md index fbbcbc4e9..ca2cb21fd 100644 --- a/resonant-circuit/README.md +++ b/resonant-circuit/README.md @@ -1,7 +1,7 @@ --- 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). --- @@ -33,6 +33,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). +* *Julia* A solver using the preCICE [Julia bindings](https://precice.org/installation-bindings-julia.html). ## Running the simulation 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/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/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