Add support for waveform relaxation with linear interpolation (#35)

* If subcycling is used, results that belong to the substeps may be exchanged, as well.
* Exchange of additional data that is associated with the substeps is realized by introducing additional data in the precice-config. `Temperature` becomes `Temperature1`, `Temperature2`,...
* Linear interpolation is used between substeps.

Be aware that this approach is not compatible with the original preCICE API and experimental. In `waveform_bindings.py` API functions like `readBlockScalarData` are extended with the additional parameter `time` (see [here](https://github.com/precice/fenics-adapter/blob/14cef37ca902b96bc4ab2b22e03c879b9e0d5a8c/fenicsadapter/waveform_bindings.py#L72-L84)).

For example usage of the modified adapter, refer to precice/tutorials#30.

.gitignore

@@ -3,3 +3,5 @@ build
 dist
 .idea
 *.pyc
+venv
+.pytest*

.travis.yml

@@ -25,4 +25,6 @@ matrix:
       - pip3 install scipy
 
 script: 
+  - mkdir -p precice && echo -e "from setuptools import setup\nsetup(name='precice')" > precice/setup.py
+  - $PY -m pip install ./precice/
   - $PY setup.py test 

README.md

@@ -29,7 +29,7 @@ Make sure to install
 
 ### Build and install the adapter
 
-Run ``python3 setup.py install --user`` from your shell.
+Run ``pip3 install --user .`` from your shell.
 
 ### Test the adapter
 

fenicsadapter/__init__.py

@@ -1,2 +1 @@
 from .fenicsadapter import Adapter, CustomExpression
-

fenicsadapter/checkpointing.py

@@ -0,0 +1,31 @@
+from .solverstate import SolverState
+
+class Checkpoint:
+
+    def __init__(self):
+        """
+        A checkpoint for the solver state
+        """
+        self._state = None
+
+    def get_state(self):
+        return self._state
+
+    def write(self, new_state):
+        """
+        write checkpoint from solver state.
+        :param u: function value
+        :param t: time
+        :param n: timestep
+        """
+        if self.is_empty():
+            self._state = SolverState(None, None, None)
+
+        self._state.copy(new_state)
+
+    def is_empty(self):
+        """
+        Returns whether checkpoint is empty. An empty checkpoint has no state saved.
+        :return:
+        """
+        return not self._state

fenicsadapter/config.py

@@ -13,6 +13,8 @@ class Config:
     :ivar _coupling_mesh_name: name of mesh as defined in preCICE config
     :ivar _read_data_name: name of read data as defined in preCICE config
     :ivar _write_data_name: name of write data as defined in preCICE config
+    :ivar _N_this: waveform relaxation substeps for this participant
+    :ivar _N_other: waveform relaxation substeps for other participant
     """
 
     def __init__(self, adapter_config_filename):
@@ -22,6 +24,8 @@ def __init__(self, adapter_config_filename):
         self._coupling_mesh_name = None
         self._read_data_name = None
         self._write_data_name = None
+        self._N_this = None
+        self._N_other = None
 
         self.readJSON(adapter_config_filename)
 
@@ -33,17 +37,27 @@ def readJSON(self, adapter_config_filename):
         :var data: data decoded from JSON files
         :var read_file: stores file path
         """
-
-        path = os.path.join(os.getcwd(), os.path.dirname(sys.argv[0]), adapter_config_filename)
+        folder = os.path.dirname(os.path.join(os.getcwd(), os.path.dirname(sys.argv[0]), adapter_config_filename))
+        path = os.path.join(folder, os.path.basename(adapter_config_filename))
         read_file = open(path, "r")
         data = json.load(read_file)
-
-        self._config_file_name = data["config_file_name"]
+        self._config_file_name = os.path.join(folder, data["config_file_name"])
         self._solver_name = data["solver_name"]
         self._coupling_mesh_name = data["interface"]["coupling_mesh_name"]
         self._write_data_name = data["interface"]["write_data_name"]
         self._read_data_name = data["interface"]["read_data_name"]
-
+        self._n_substeps = data["waveform"]["n_substeps"]
+        # todo check that either both keys (N_this and N_other) exist or not.
+        try:
+            self._N_this = data["interface"]["N_this"]
+            assert(self._N_this > 0)
+        except KeyError:
+            self._N_this = None
+        try:
+            self._N_other = data["interface"]["N_other"]
+            assert (self._N_other > 0)
+        except KeyError:
+            self._N_this = None
         read_file.close()
 
     def get_config_file_name(self):
@@ -60,3 +74,7 @@ def get_read_data_name(self):
 
     def get_write_data_name(self):
         return self._write_data_name
+
+    def get_n_substeps(self):
+        return self._n_substeps
+

fenicsadapter/fenicsadapter.py

@@ -4,25 +4,20 @@
 :raise ImportError: if PRECICE_ROOT is not defined
 """
 import dolfin
-from dolfin import UserExpression, SubDomain, Function
+from dolfin import UserExpression, SubDomain, Function, Constant
 from scipy.interpolate import Rbf
 from scipy.interpolate import interp1d
 import numpy as np
 from .config import Config
+from .checkpointing import Checkpoint
+from .solverstate import SolverState
 
-try:
-    import precice
-except ImportError:
-    import os
-    import sys
-    # check if PRECICE_ROOT is defined
-    if not os.getenv('PRECICE_ROOT'):
-       raise Exception("ERROR: PRECICE_ROOT not defined!")
 
-    precice_root = os.getenv('PRECICE_ROOT')
-    precice_python_adapter_root = precice_root+"/src/precice/bindings/python"
-    sys.path.insert(0, precice_python_adapter_root)
-    import precice
+import fenicsadapter.waveform_bindings
+
+import logging
+
+logging.basicConfig(level=logging.INFO)
 
 
 class CustomExpression(UserExpression):
@@ -56,11 +51,11 @@ def rbf_interpol(self, x):
             return f(x[0], x[1], x[2])
 
     def lin_interpol(self, x):
-        f = interp1d(self._coords_x, self._vals)
-        return f(x.x())
+        f = interp1d(self._coords_y, self._vals, bounds_error=False, fill_value="extrapolate")
+        return f(x[1])
 
     def eval(self, value, x):
-        value[0] = self.rbf_interpol(x)
+        value[0] = self.lin_interpol(x)
 
 
 class Adapter:
@@ -69,44 +64,41 @@ class Adapter:
 
     :ivar _config: object of class Config, which stores data from the JSON config file
     """
-    def __init__(self, adapter_config_filename='precice-adapter-config.json'):
+    def __init__(self, adapter_config_filename, other_adapter_config_filename):
 
         self._config = Config(adapter_config_filename)
 
         self._solver_name = self._config.get_solver_name()
 
-        self._interface = precice.Interface(self._solver_name, 0, 1)
+        self._interface = fenicsadapter.waveform_bindings.WaveformBindings(self._solver_name, 0, 1)
+        self._interface.configure_waveform_relaxation(adapter_config_filename, other_adapter_config_filename)
         self._interface.configure(self._config.get_config_file_name())
         self._dimensions = self._interface.get_dimensions()
 
-        self._coupling_subdomain = None # initialized later
-        self._mesh_fenics = None # initialized later
-        self._coupling_bc_expression = None # initialized later
+        self._coupling_subdomain = None  # initialized later
+        self._mesh_fenics = None  # initialized later
+        self._coupling_bc_expression = None  # initialized later
 
-        ## coupling mesh related quantities
-        self._coupling_mesh_vertices = None # initialized later
+        # coupling mesh related quantities
+        self._coupling_mesh_vertices = None  # initialized later
         self._mesh_name = self._config.get_coupling_mesh_name()
         self._mesh_id = self._interface.get_mesh_id(self._mesh_name)
         self._vertex_ids = None # initialized later
         self._n_vertices = None # initialized later
 
-        ## write data related quantities (write data is written by this solver to preCICE)
+        # write data related quantities (write data is written by this solver to preCICE)
         self._write_data_name = self._config.get_write_data_name()
-        self._write_data_id = self._interface.get_data_id(self._write_data_name, self._mesh_id)
         self._write_data = None
 
-        ## read data related quantities (read data is read by this solver from preCICE)
+        # read data related quantities (read data is read by this solver from preCICE)
         self._read_data_name = self._config.get_read_data_name()
-        self._read_data_id = self._interface.get_data_id(self._read_data_name, self._mesh_id)
         self._read_data = None
 
-        ## numerics
+        # numerics
         self._precice_tau = None
 
         ## checkpointing
-        self._u_cp = None  # checkpoint for temperature inside domain
-        self._t_cp = None  # time of the checkpoint
-        self._n_cp = None  # timestep of the checkpoint
+        self._checkpoint = Checkpoint()
 
     def convert_fenics_to_precice(self, data, mesh, subdomain):
         """Converts FEniCS data of type dolfin.Function into
@@ -186,6 +178,7 @@ def create_coupling_boundary_condition(self):
             self._coupling_bc_expression = CustomExpression()
         except (TypeError, KeyError):  # works with dolfin 2017.2.0
             self._coupling_bc_expression = CustomExpression(degree=0)
+
         self._coupling_bc_expression.set_boundary_data(self._read_data, x_vert, y_vert)
 
     def create_coupling_dirichlet_boundary_condition(self, function_space):
@@ -208,12 +201,40 @@ def create_coupling_neumann_boundary_condition(self, test_functions):
         self.create_coupling_boundary_condition()
         return self._coupling_bc_expression * test_functions * dolfin.ds  # this term has to be added to weak form to add a Neumann BC (see e.g. p. 83ff Langtangen, Hans Petter, and Anders Logg. "Solving PDEs in Python The FEniCS Tutorial Volume I." (2016).)
 
+    def _restore_solver_state_from_checkpoint(self, state):
+        """Resets the solver's state to the checkpoint's state.
+        :param state: current state of the FEniCS solver
+        """
+        logging.debug("Restore solver state")
+        state.update(self._checkpoint.get_state())
+        self._interface.fulfilled_action(fenicsadapter.waveform_bindings.action_read_iteration_checkpoint())
+
+    def _advance_solver_state(self, state, u_np1, dt):
+        """Advances the solver's state by one timestep.
+        :param state: old state
+        :param u_np1: new value
+        :param dt: timestep size
+        :return:
+        """
+        logging.debug("Advance solver state")
+        logging.debug("old state: t={time}".format(time=state.t))
+        state.update(SolverState(u_np1, state.t + dt, self._checkpoint.get_state().n + 1))
+        logging.debug("new state: t={time}".format(time=state.t))
+
+    def _save_solver_state_to_checkpoint(self, state):
+        """Writes given solver state to checkpoint.
+        :param state: state being saved as checkpoint
+        """
+        logging.debug("Save solver state")
+        self._checkpoint.write(state)
+        self._interface.fulfilled_action(fenicsadapter.waveform_bindings.action_write_iteration_checkpoint())
+
     def advance(self, write_function, u_np1, u_n, t, dt, n):
         """Calls preCICE advance function using precice and manages checkpointing.
         The solution u_n is updated by this function via call-by-reference. The corresponding values for t and n are returned.
 
         This means:
-        * either, the checkpoint self._u_cp is assigned to u_n to repeat the iteration,
+        * either, the olf value of the checkpoint is assigned to u_n to repeat the iteration,
         * or u_n+1 is assigned to u_n and the checkpoint is updated correspondingly.
 
         :param write_function: a FEniCS function being sent to the other participant as boundary condition at the coupling interface
@@ -225,46 +246,42 @@ def advance(self, write_function, u_np1, u_n, t, dt, n):
         :return: return starting time t and timestep n for next FEniCS solver iteration. u_n is updated by advance correspondingly.
         """
 
+        state = SolverState(u_n, t, n)
+
         # sample write data at interface
         x_vert, y_vert = self.extract_coupling_boundary_coordinates()
         self._write_data = self.convert_fenics_to_precice(write_function, self._mesh_fenics, self._coupling_subdomain)
-
-        # communication
-        self._interface.write_block_scalar_data(self._write_data_id, self._n_vertices, self._vertex_ids, self._write_data)
+        if True:  # todo: add self._interface.is_write_data_required(dt). We should add this check. However, it is currently not properly implemented for waveform relaxation
+            self._interface.write_block_scalar_data(self._write_data_name, self._mesh_id, self._n_vertices, self._vertex_ids, self._write_data, t + dt)
         max_dt = self._interface.advance(dt)
-        self._interface.read_block_scalar_data(self._read_data_id, self._n_vertices, self._vertex_ids, self._read_data)
-
-        # update boundary condition with read data
-        self._coupling_bc_expression.update_boundary_data(self._read_data, x_vert, y_vert)
 
         precice_step_complete = False
-
+        solver_state_has_been_restored = False
+
         # checkpointing
-        if self._interface.is_action_required(precice.action_read_iteration_checkpoint()):
-            # continue FEniCS computation from checkpoint
-            u_n.assign(self._u_cp)  # set u_n to value of checkpoint
-            t = self._t_cp
-            n = self._n_cp
-            self._interface.fulfilled_action(precice.action_read_iteration_checkpoint())
+        if self._interface.is_action_required(fenicsadapter.waveform_bindings.action_read_iteration_checkpoint()):
+            self._restore_solver_state_from_checkpoint(state)
+            solver_state_has_been_restored = True
         else:
-            u_n.assign(u_np1)
-            t = new_t = t + dt  # todo the variables new_t, new_n could be saved, by just using t and n below, however I think it improved readability.
-            n = new_n = n + 1
-
-        if self._interface.is_action_required(precice.action_write_iteration_checkpoint()):
-            # continue FEniCS computation with u_np1
-            # update checkpoint
-            self._u_cp.assign(u_np1)
-            self._t_cp = new_t
-            self._n_cp = new_n
-            self._interface.fulfilled_action(precice.action_write_iteration_checkpoint())
+            self._advance_solver_state(state, u_np1, dt)
+
+        if self._interface.is_action_required(fenicsadapter.waveform_bindings.action_write_iteration_checkpoint()):
+            assert (not solver_state_has_been_restored)  # avoids invalid control flow
+            self._save_solver_state_to_checkpoint(state)
             precice_step_complete = True
 
+        _, t, n = state.get_state()
+
+        if True:  # todo: add self._interface.is_read_data_available().  We should add this check. However, it is currently not properly implemented for waveform relaxation
+            self._interface.read_block_scalar_data(self._read_data_name, self._mesh_id, self._n_vertices, self._vertex_ids, self._read_data, t + dt)  # if precice_step_complete, we have to already use the new t for reading. Otherwise, we get a lag. Therefore, this command has to be called AFTER the state has been updated/recovered.
+        print(self._read_data)
+
+        self._coupling_bc_expression.update_boundary_data(self._read_data, x_vert, y_vert)
+
         return t, n, precice_step_complete, max_dt
 
     def initialize(self, coupling_subdomain, mesh, read_field, write_field, u_n, t=0, n=0):
         """Initializes remaining attributes. Called once, from the solver.
-
         :param read_field: function applied on the read field
         :param write_field: function applied on the write field
         """
@@ -273,22 +290,28 @@ def initialize(self, coupling_subdomain, mesh, read_field, write_field, u_n, t=0
         self.set_write_field(write_field)
         self._precice_tau = self._interface.initialize()
 
-        if self._interface.is_action_required(precice.action_write_initial_data()):
-            self._interface.write_block_scalar_data(self._write_data_id, self._n_vertices, self._vertex_ids, self._write_data)
-            self._interface.fulfilled_action(precice.action_write_initial_data())
+        dt = Constant(0)
+        self.fenics_dt = self._precice_tau / self._config.get_n_substeps()
+        dt.assign(np.min([self.fenics_dt, self._precice_tau]))
+
+        self._interface.initialize_waveforms(self._mesh_id, self._n_vertices, self._vertex_ids, self._write_data_name,
+                                             self._read_data_name)
+
+        if self._interface.is_action_required(fenicsadapter.waveform_bindings.action_write_initial_data()):
+            self._interface.write_block_scalar_data(self._write_data_name, self._mesh_id, self._n_vertices, self._vertex_ids, self._write_data, t)
+            self._interface.fulfilled_action(fenicsadapter.waveform_bindings.action_write_initial_data())
 
-        self._interface.initialize_data()
+        self._interface.initialize_data(read_zero=self._read_data, write_zero=self._write_data)
 
         if self._interface.is_read_data_available():
-            self._interface.read_block_scalar_data(self._read_data_id, self._n_vertices, self._vertex_ids, self._read_data)
+            self._interface.read_block_scalar_data(self._read_data_name, self._mesh_id, self._n_vertices, self._vertex_ids, self._read_data, t + dt(0))
+            print(self._read_data)
 
-        if self._interface.is_action_required(precice.action_write_iteration_checkpoint()):
-            self._u_cp = u_n.copy(deepcopy=True)
-            self._t_cp = t
-            self._n_cp = n
-            self._interface.fulfilled_action(precice.action_write_iteration_checkpoint())
+        if self._interface.is_action_required(fenicsadapter.waveform_bindings.action_write_iteration_checkpoint()):
+            initial_state = SolverState(u_n, t, n)
+            self._save_solver_state_to_checkpoint(initial_state)
 
-        return self._precice_tau
+        return dt
 
     def is_coupling_ongoing(self):
         """Determines whether simulation should continue. Called from the