diff --git a/oscillator-overlap/README.md b/oscillator-overlap/README.md new file mode 100644 index 000000000..1897dae5a --- /dev/null +++ b/oscillator-overlap/README.md @@ -0,0 +1,42 @@ +--- +title: Oscillator (overlapping domain decomposition) +permalink: tutorials-oscillator-overlap.html +keywords: Python, ODE +summary: We solve an oscillator with two masses in a partitioned fashion with overlapping domain decomposition. Each mass is solved by an independent ODE. +--- + +{% note %} +Get the [case files of this tutorial](https://github.com/precice/tutorials/tree/master/oscillator-overlap). Read how in the [tutorials introduction](https://www.precice.org/tutorials.html). +{% endnote %} + +## Setup + +This tutorial solves the same problem as the [oscillator tutorial](https://precice.org/tutorials-oscillator.html), but applies a different domain decomposition strategy. See the oscillator tutorial for details on the general setup. The partitioning of the mass-spring system is shown here: + +![Schematic drawing of oscillator example with overlapping domain decomposition](images/tutorials-oscillator-overlap-dd.png) + +Note that this case applies an overlapping 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]. + +## 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 `. + +## Running the Simulation + +### Python + +Open two separate terminals and start both participants by calling: + +```bash +cd python +./run.sh -l +``` + +and + +```bash +cd python +./run.sh -r +``` diff --git a/oscillator-overlap/clean-tutorial.sh b/oscillator-overlap/clean-tutorial.sh new file mode 120000 index 000000000..4713f5092 --- /dev/null +++ b/oscillator-overlap/clean-tutorial.sh @@ -0,0 +1 @@ +../tools/clean-tutorial-base.sh \ No newline at end of file diff --git a/oscillator-overlap/images/tutorials-oscillator-overlap-dd.pdf b/oscillator-overlap/images/tutorials-oscillator-overlap-dd.pdf new file mode 100644 index 000000000..6d9cad947 Binary files /dev/null and b/oscillator-overlap/images/tutorials-oscillator-overlap-dd.pdf differ diff --git a/oscillator-overlap/images/tutorials-oscillator-overlap-dd.png b/oscillator-overlap/images/tutorials-oscillator-overlap-dd.png new file mode 100644 index 000000000..d794e16d3 Binary files /dev/null and b/oscillator-overlap/images/tutorials-oscillator-overlap-dd.png differ diff --git a/oscillator-overlap/images/tutorials-oscillator-overlap-dd.tex b/oscillator-overlap/images/tutorials-oscillator-overlap-dd.tex new file mode 100644 index 000000000..c84a3d74b --- /dev/null +++ b/oscillator-overlap/images/tutorials-oscillator-overlap-dd.tex @@ -0,0 +1,37 @@ +\documentclass{standalone} +\usepackage{tikz} +\usetikzlibrary{positioning,calc,patterns,decorations.pathmorphing} + +\begin{document} + +\begin{tikzpicture} + \tikzstyle{spring}=[thick,decorate,decoration={zigzag,pre length=0.3cm,post length=0.3cm,segment length=6}] + \tikzstyle{ground}=[fill,pattern=north east lines,draw=none,minimum width=0.75cm,minimum height=0.3cm] + \tikzstyle{mass}=[draw, fill=gray!25, thick, minimum width=1cm, minimum height=1cm] + \node (W1) [ground, rotate=-90, minimum width=3cm] {}; + \node (M1) [mass, right= 2cm of W1.north] {$m_1$}; + \node (M2fake) [circle, draw=black, fill=gray!25, right= 2cm of M1] {}; + + \draw [dashed] (M2fake.north) -- ++ (0,.5); + \draw [->] ($(M2fake.north) + (0,.25)$) -- node[above](U2ToM1){$u_2$} ++ (1,0); + \draw [dashed] (M1.north) -- ++ (0,.5); + \draw [->] ($(M1.north) + (0,.25)$) -- node[above](U1FromM1){$u_1$} ++ (1,0); + \draw [spring] (W1.north) -- node[above]{$k_1$}(M1.west); + \draw [spring] (M1.east) -- node[above]{$k_{12}$}(M2fake.west); + + \node (M1fake) [circle, draw=black, fill=gray!25, right= 1cm of M2fake] {}; + \node (M2) [mass, right= 2cm of M1fake] {$m_2$}; + \draw [dashed] (M2.north) -- ++ (0,.5); + \draw [->] ($(M2.north) + (0,.25)$) -- node[above](U2FromM2){$u_2$} ++ (1,0); + \draw [dashed] (M1fake.north) -- ++ (0,.5); + \draw [->] ($(M1fake.north) + (0,.25)$) -- node[above](U1ToM2){$u_1$} ++ (1,0); + + \node (W2) [ground, rotate=90, minimum width=3cm, right= 2cm of M2, xshift=-1.5cm] {}; + \draw [spring] (M2.east) -- node[above]{$k_2$}(W2.north); + \draw [spring] (M1fake.east) -- node[below]{$k_{12}$}(M2.west); + + \draw[](U2FromM2.north) edge[->,out=135, in=45] node[above right]{ghost layer for $m_1$} (U2ToM1.north); + \draw[](U1FromM1.north) edge[->,out=45, in=135] node[above left]{ghost layer for $m_2$} (U1ToM2.north); + +\end{tikzpicture} +\end{document} diff --git a/oscillator-overlap/mass-left-python/clean.sh b/oscillator-overlap/mass-left-python/clean.sh new file mode 100755 index 000000000..c7e3552ac --- /dev/null +++ b/oscillator-overlap/mass-left-python/clean.sh @@ -0,0 +1,8 @@ +#!/bin/sh +set -e -u + +. ../../tools/cleaning-tools.sh + +rm -rfv ./output/ + +clean_precice_logs . diff --git a/oscillator-overlap/mass-left-python/run.sh b/oscillator-overlap/mass-left-python/run.sh new file mode 100755 index 000000000..2784cfeba --- /dev/null +++ b/oscillator-overlap/mass-left-python/run.sh @@ -0,0 +1,4 @@ +#!/bin/sh +set -e -u + +python3 ../solver-python/oscillator.py Mass-Left \ No newline at end of file diff --git a/oscillator-overlap/mass-right-python/clean.sh b/oscillator-overlap/mass-right-python/clean.sh new file mode 100755 index 000000000..c7e3552ac --- /dev/null +++ b/oscillator-overlap/mass-right-python/clean.sh @@ -0,0 +1,8 @@ +#!/bin/sh +set -e -u + +. ../../tools/cleaning-tools.sh + +rm -rfv ./output/ + +clean_precice_logs . diff --git a/oscillator-overlap/mass-right-python/run.sh b/oscillator-overlap/mass-right-python/run.sh new file mode 100755 index 000000000..775d6a83e --- /dev/null +++ b/oscillator-overlap/mass-right-python/run.sh @@ -0,0 +1,4 @@ +#!/bin/sh +set -e -u + +python3 ../solver-python/oscillator.py Mass-Right \ No newline at end of file diff --git a/oscillator-overlap/precice-config.xml b/oscillator-overlap/precice-config.xml new file mode 100644 index 000000000..6bf072cb5 --- /dev/null +++ b/oscillator-overlap/precice-config.xml @@ -0,0 +1,47 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/oscillator-overlap/solver-python/oscillator.py b/oscillator-overlap/solver-python/oscillator.py new file mode 100644 index 000000000..52bc3fa4c --- /dev/null +++ b/oscillator-overlap/solver-python/oscillator.py @@ -0,0 +1,190 @@ +from __future__ import division + +import argparse +import numpy as np +from numpy.linalg import eig +import precice +from enum import Enum +import csv +import os + +import problemDefinition + + +class Scheme(Enum): + NEWMARK_BETA = "Newmark_beta" + GENERALIZED_ALPHA = "generalized_alpha" + + +class Participant(Enum): + MASS_LEFT = "Mass-Left" + MASS_RIGHT = "Mass-Right" + + +parser = argparse.ArgumentParser() +parser.add_argument("participantName", help="Name of the solver.", type=str, choices=[p.value for p in Participant]) +parser.add_argument("-ts", "--time-stepping", help="Time stepping scheme being used.", type=str, + choices=[s.value for s in Scheme], default=Scheme.NEWMARK_BETA.value) +args = parser.parse_args() + +participant_name = args.participantName + +if participant_name == Participant.MASS_LEFT.value: + write_data_name = 'Displacement-Left' + read_data_name = 'Displacement-Right' + mesh_name = 'Mass-Left-Mesh' + + this_mass = problemDefinition.MassLeft + this_spring = problemDefinition.SpringLeft + connecting_spring = problemDefinition.SpringMiddle + other_mass = problemDefinition.MassRight + +elif participant_name == Participant.MASS_RIGHT.value: + read_data_name = 'Displacement-Left' + write_data_name = 'Displacement-Right' + mesh_name = 'Mass-Right-Mesh' + + this_mass = problemDefinition.MassRight + this_spring = problemDefinition.SpringRight + connecting_spring = problemDefinition.SpringMiddle + other_mass = problemDefinition.MassLeft + +else: + raise Exception(f"wrong participant name: {participant_name}") + +mass = this_mass.m +stiffness = this_spring.k + connecting_spring.k +u0, v0, f0, d_dt_f0 = this_mass.u0, this_mass.v0, connecting_spring.k * other_mass.u0, connecting_spring.k * other_mass.v0 + +num_vertices = 1 # Number of vertices + +solver_process_index = 0 +solver_process_size = 1 + +configuration_file_name = "../precice-config.xml" + +participant = precice.Participant(participant_name, configuration_file_name, solver_process_index, solver_process_size) + +dimensions = participant.get_mesh_dimensions(mesh_name) + +vertex = np.zeros(dimensions) +read_data = np.zeros(num_vertices) +write_data = connecting_spring.k * u0 * np.ones(num_vertices) + +vertex_ids = [participant.set_mesh_vertex(mesh_name, vertex)] + +if participant.requires_initial_data(): + participant.write_data(mesh_name, write_data_name, vertex_ids, write_data) + +participant.initialize() +precice_dt = participant.get_max_time_step_size() +my_dt = precice_dt # use my_dt < precice_dt for subcycling + +# Initial Conditions +a0 = (f0 - stiffness * u0) / mass +u = u0 +v = v0 +a = a0 +t = 0 + +# Generalized Alpha Parameters +if args.time_stepping == Scheme.GENERALIZED_ALPHA.value: + alpha_f = 0.4 + alpha_m = 0.2 +elif args.time_stepping == Scheme.NEWMARK_BETA.value: + alpha_f = 0.0 + alpha_m = 0.0 +m = 3 * [None] # will be computed for each timestep depending on dt +gamma = 0.5 - alpha_m + alpha_f +beta = 0.25 * (gamma + 0.5) + +positions = [] +velocities = [] +times = [] + +u_write = [u] +v_write = [v] +t_write = [t] + +while participant.is_coupling_ongoing(): + if participant.requires_writing_checkpoint(): + u_cp = u + v_cp = v + a_cp = a + t_cp = t + # store data for plotting and postprocessing + positions += u_write + velocities += v_write + times += t_write + + # compute time step size for this time step + precice_dt = participant.get_max_time_step_size() + dt = np.min([precice_dt, my_dt]) + read_time = (1 - alpha_f) * dt + read_data = participant.read_data(mesh_name, read_data_name, vertex_ids, read_time) + f = connecting_spring.k * read_data[0] + + # do generalized alpha step + m[0] = (1 - alpha_m) / (beta * dt**2) + m[1] = (1 - alpha_m) / (beta * dt) + m[2] = (1 - alpha_m - 2 * beta) / (2 * beta) + k_bar = stiffness * (1 - alpha_f) + m[0] * mass + u_new = (f - alpha_f * stiffness * u + mass * (m[0] * u + m[1] * v + m[2] * a)) / k_bar + a_new = 1.0 / (beta * dt**2) * (u_new - u - dt * v) - (1 - 2 * beta) / (2 * beta) * a + v_new = v + dt * ((1 - gamma) * a + gamma * a_new) + t_new = t + dt + + write_data = [u_new] + + participant.write_data(mesh_name, write_data_name, vertex_ids, write_data) + + participant.advance(dt) + + if participant.requires_reading_checkpoint(): + u = u_cp + v = v_cp + a = a_cp + t = t_cp + # empty buffers for next window + u_write = [] + v_write = [] + t_write = [] + + else: + u = u_new + v = v_new + a = a_new + t = t_new + + # write data to buffers + u_write.append(u) + v_write.append(v) + t_write.append(t) + +# store final result +u = u_new +v = v_new +a = a_new +u_write.append(u) +v_write.append(v) +t_write.append(t) +positions += u_write +velocities += v_write +times += t_write + +participant.finalize() + +# print errors +error = np.max(abs(this_mass.u_analytical(np.array(times)) - np.array(positions))) +print("Error w.r.t analytical solution:") +print(f"{my_dt},{error}") + +# output trajectory +if not os.path.exists("output"): + os.makedirs("output") + +with open(f'output/trajectory-{participant_name}.csv', 'w') as file: + csv_write = csv.writer(file, delimiter=';') + csv_write.writerow(['time', 'position', 'velocity']) + for t, u, v in zip(times, positions, velocities): + csv_write.writerow([t, u, v]) diff --git a/oscillator-overlap/solver-python/problemDefinition.py b/oscillator-overlap/solver-python/problemDefinition.py new file mode 100644 index 000000000..b0ac2aea3 --- /dev/null +++ b/oscillator-overlap/solver-python/problemDefinition.py @@ -0,0 +1,68 @@ +import numpy as np +from numpy.linalg import eig + + +class SpringLeft: + k = 4 * np.pi**2 + + +class SpringMiddle: + k = 16 * (np.pi**2) + + +class SpringRight: + k = 4 * np.pi**2 + + +class MassLeft: + # mass + m = 1 + + # initial conditions (the way how we currently compute the analytical + # solution allows arbitrary u0, but requires v0 = 0) + u0 = 1.0 + v0 = 0.0 + + u_analytical, v_analytical = None, None # will be defined below + + +class MassRight: + # mass + m = 1 + + # initial conditions (the way how we currently compute the analytical + # solution allows arbitrary u0, but requires v0 = 0) + u0 = 0.0 + v0 = 0.0 + + u_analytical, v_analytical = None, None # will be defined below + + +# Mass matrix +M = np.array([ + [MassLeft.m, 0], + [0, MassRight.m] +]) +# Stiffness matrix +K = np.array([ + [SpringLeft.k + SpringMiddle.k, -SpringMiddle.k], + [-SpringMiddle.k, SpringRight.k + SpringMiddle.k] +]) + +# system: +# m ddu + k u = f +# compute analytical solution from eigenvalue ansatz + +eigenvalues, eigenvectors = eig(K) +omega = np.sqrt(eigenvalues) +A, B = eigenvectors + +c = np.linalg.solve(eigenvectors, [MassLeft.u0, MassRight.u0]) + +MassLeft.u_analytical = lambda t: c[0] * A[0] * np.cos(omega[0] * t) + c[1] * A[1] * np.cos(omega[1] * t) +MassLeft.v_analytical = lambda t: -c[0] * A[0] * omega[0] * \ + np.sin(omega[0] * t) - c[1] * A[1] * omega[1] * np.sin(omega[1] * t) + +MassRight.u_analytical = lambda t: c[0] * B[0] * np.cos(omega[0] * t) + c[1] * B[1] * np.cos(omega[1] * t) +MassRight.v_analytical = lambda t: -c[0] * B[0] * omega[0] * \ + np.sin(omega[0] * t) - c[1] * B[1] * omega[1] * np.sin(omega[1] * t) diff --git a/oscillator/precice-config.xml b/oscillator/precice-config.xml index 0c54ea3d0..80cdc7995 100644 --- a/oscillator/precice-config.xml +++ b/oscillator/precice-config.xml @@ -1,8 +1,5 @@ -