Added code

Author: BenjaminRodenberg

README.md

@@ -0,0 +1,39 @@
+# Running the code
+
+Run the code by typing ```python main.py```. You need numpy and scipy for executing the code. Output folders with experimental results are created. The folders are named in the following way: ```seriesX_SHORTNAME_DATETIME```. See below for details on the different experimental series.
+
+You can plot the results by running ```python postproc.py <foldername>```. You need matplotlib for plotting.
+
+# Convergence Study
+
+This code can be used to perform convergence studies in a partitioned multi-physics multi-scale setup. As a benchmark scenario the 1D heat transport equations is solved. For details see [1]. The convergence studies performed by this code are identical to the ones described in [1]. However, for a lower runtime, the discretization slightly differs for [1].
+
+## Experimental Series 1
+
+Order degradation to first order if classical implicit coupling is used with higher order schemes a regular domain
+
+![](./series1.png)
+
+## Experimental Series 2
+
+Second order if customized coupling schemes are used with second order schemes a regular domain
+
+![](./series2.png)
+
+## Experimental Series 3
+
+Second order if Strang splitting coupling is used with higher order schemes on a non-regular domain
+
+![](./series3.png)
+
+## Experimental Series 4
+
+High order if waveform relaxation coupling is used with higher order schemes on a non-regular domain
+
+![](./series4.png)
+
+# Reference
+
+If you want to refer to this code, please use the following reference:
+
+[1] Rüth, B., Uekermann, B., Mehl, M., & Bungartz, H.-J. (2018). Time Stepping Algorithms for Partitioned Multi-Scale Multi-Physics. In ??? Glasgow: ???

coupling_schemes.py

@@ -0,0 +1,71 @@
+import numpy as np
+import abc
+from time_integration import TimeIntegrationScheme
+
+
+class Domain():
+    __metaclass__ = abc.ABCMeta
+
+    def __init__(self, grid):
+        """
+        Simulation domain with a grid, corresponding function values u, boundary conditions and a time integration scheme
+        :param grid: grid discretizing domain
+        :type grid: OneDimensionalGrid
+        """
+        self.grid = grid
+        self.u = np.zeros(self.grid.N_gridpoints)  # :type : np.array
+        self.left_BC = {}
+        self.right_BC = {}
+        self.time_integration_scheme = None  # :type : TimeIntegrationScheme
+
+    def update_u(self, new_u):
+        """
+        updates u on the domain
+        :param new_u:
+        :type new_u: np.array
+        :return:
+        """
+        assert self.u.shape == new_u.shape
+        self.u = new_u
+
+
+class OneDimensionalGrid():
+    __metaclass__ = abc.ABCMeta
+
+    def __init__(self, x):
+        """
+        :param x: array of points on 1D grid
+        """
+        self.x = np.array(x)
+        self.N_gridpoints = self.x.__len__()  # :type : int
+        self.x_left = self.x[0]  # :type : float
+        self.x_right = self.x[-1]  # :type : float
+
+
+class RegularGrid(OneDimensionalGrid):
+
+    def __init__(self, N_gridpoints, x_left, x_right):
+        """
+        Regular grid with N_gridpoints between x_left and x_right
+        :param N_gridpoints: number of gridpoints
+        :type N_gridpoints: int
+        :param x_left: leftmost coordinate
+        :type x_left: float
+        :param x_right: rightmost coordinate
+        :type x_right: float
+        """
+        # meshwidth of grid
+        self.h = (x_right - x_left) / (N_gridpoints-1)  # :type : float
+        # vector of mesh points
+        x = np.linspace(x_left, x_right, N_gridpoints)  # :type : np.array
+        super(RegularGrid, self).__init__(x)
+
+
+class IrregularGrid(OneDimensionalGrid):
+
+    def __init__(self, x):
+        """
+        Irregular grid defined by ascending vector of gridpoints
+        :param x: vector of gridpoints
+        """
+        super(IrregularGrid, self).__init__(x)

domain.py

@@ -0,0 +1,23 @@
+import datetime
+import json
+import os
+
+
+class Experiment():
+    def __init__(self, name, h, tau, errors_left, errors_right, additional_numerical_parameters = None):
+        self.h = h
+        self.temporal_resolution = tau
+        self.numerical_parameters = additional_numerical_parameters
+        self.errors_left = errors_left
+        self.errors_right = errors_right
+        self.time = datetime.datetime.now().strftime("%Y-%m-%d--%H-%M-%S_%f")
+        self.experiment_name = name
+
+    def save(self, save_path):
+        if not os.path.exists(save_path):  # create empty directory
+            os.mkdir(save_path)
+
+        filename = save_path + "/" + self.experiment_name + "_" + self.time + ".json"
+        with open(filename, 'w') as outfile:
+            s = json.dumps(self.__dict__, outfile)
+            outfile.write(s)

experiment_documentation.py

@@ -0,0 +1,58 @@
+def central_FD_at(i, h, u, derivative=1, order=2):
+    """
+    computes the n_th central finite difference at x_i on the basis of u_i.
+    :param i:
+    :param h:
+    :param u:
+    :param derivative:
+    :param order:
+    :return:
+    """
+
+    du = 0
+
+    if derivative == 1:
+        if order == 2:
+            du = 1.0/(2.0 * h) * (-1.0 * u[i-1] + 1.0 * u[i+1])
+        else:
+            print "not implemented!"
+    else:
+        print "not implemented!"
+
+    return du
+
+
+def one_sided_forward_FD_at(i, h, u, derivative=1, order=1):
+    """
+    computes the n_th one-sided forward finite difference derivative at x_i on the basis of u_i.
+
+    Coefficients from https://en.wikipedia.org/wiki/Finite_difference_coefficient
+
+    d^n u/dn (x_i)= a*u_i + b*u_i+1 + c*u_i+2 + ...
+
+    :param i:
+    :param u:
+    :param order:
+    :param derivative:
+    :return:
+    """
+
+    du = 0
+
+    if derivative == 1:
+        if order == 1:
+            du = 1.0/h * (- 1.0 * u[i] + 1.0 * u[i+1] )
+        elif order == 2:
+            du = 1.0/h * (-1.5 * u[i] + 2.0 * u[i+1] - 0.5 * u[i+2])
+        elif order == 3:
+            du = 1.0/h * (-11.0/6.0 * u[i] + 3.0 * u[i+1] - 3.0/2.0 * u[i+2] + 1.0/3.0 * u[i+3])
+        elif order == 4:
+            du = 1.0/h * (-25.0/12.0 * u[i] + 4.0 * u[i+1] - 3.0 * u[i+2] + 4.0/3.0 * u[i+3] - 1.0/4.0 * u[i+4])
+        elif order == 5:
+            du = 1.0/h * (-137.0/60.0 * u[i] + 5.0 * u[i+1] - 5.0 * u[i+2] + 10.0/3.0 * u[i+3] - 5.0/4.0 * u[i+4] + 1.0/5.0 * u[i+5])
+        else:
+            print "not implemented!"
+    else:
+        print "not implemented!"
+
+    return du