diff --git a/partitioned-heat-conduction/nutils/clean.sh b/partitioned-heat-conduction/dirichlet-nutils/clean.sh
similarity index 100%
rename from partitioned-heat-conduction/nutils/clean.sh
rename to partitioned-heat-conduction/dirichlet-nutils/clean.sh
diff --git a/partitioned-heat-conduction/nutils/heat.py b/partitioned-heat-conduction/dirichlet-nutils/heat.py
similarity index 53%
rename from partitioned-heat-conduction/nutils/heat.py
rename to partitioned-heat-conduction/dirichlet-nutils/heat.py
index 09f631fb3..dc80100e6 100644
--- a/partitioned-heat-conduction/nutils/heat.py
+++ b/partitioned-heat-conduction/dirichlet-nutils/heat.py
@@ -7,20 +7,14 @@
 import precice
 
 
-def main(side='Dirichlet', n=10, degree=1, timestep=.1, alpha=3., beta=1.2):
-
-    if side == 'Dirichlet':
-        x_grid = np.linspace(0, 1, n)
-    elif side == 'Neumann':
-        x_grid = np.linspace(1, 2, n)
-    else:
-        raise Exception('invalid side {!r}'.format(side))
+def main(n=10, degree=1, timestep=.1, alpha=3., beta=1.2):
+
+    x_grid = np.linspace(0, 1, n)
     y_grid = np.linspace(0, 1, n)
 
     # define the Nutils mesh
     domain, geom = mesh.rectilinear([x_grid, y_grid])
-    coupling_boundary = domain.boundary['right' if side ==
-                                        'Dirichlet' else 'left']
+    coupling_boundary = domain.boundary['right']
     coupling_sample = coupling_boundary.sample('gauss', degree=degree * 2)
 
     # Nutils namespace
@@ -44,56 +38,51 @@ def main(side='Dirichlet', n=10, degree=1, timestep=.1, alpha=3., beta=1.2):
 
     # set boundary conditions at non-coupling boundaries
     # top and bottom boundary are non-coupling for both sides
-    sqr = domain.boundary['top,bottom,left' if side == 'Dirichlet' else 'top,bottom,right'].integral(
-        '(u - uexact)^2 d:x' @ ns, degree=degree * 2)
+    sqr = domain.boundary['top,bottom,left'].integral('(u - uexact)^2 d:x' @ ns, degree=degree * 2)
 
-    if side == 'Dirichlet':
-        sqr += coupling_sample.integral('(u - readfunc)^2 d:x' @ ns)
-    else:
-        res += coupling_sample.integral('basis_n readfunc d:x' @ ns)
+    sqr += coupling_sample.integral('(u - readfunc)^2 d:x' @ ns)
 
     # preCICE setup
-    participant = precice.Participant(side, "../precice-config.xml", 0, 1)
-    mesh_name = side + "-Mesh"
+    participant = precice.Participant('Dirichlet', "../precice-config.xml", 0, 1)
+    mesh_name = "Dirichlet-Mesh"
     vertex_ids = participant.set_mesh_vertices(
         mesh_name, coupling_sample.eval(ns.x))
     precice_write = functools.partial(
         participant.write_data,
         mesh_name,
-        "Temperature" if side == "Neumann" else "Heat-Flux",
+        "Heat-Flux",
         vertex_ids)
     precice_read = functools.partial(
         participant.read_data,
         mesh_name,
-        "Heat-Flux" if side == "Neumann" else "Temperature",
+        "Temperature",
         vertex_ids)
 
     # helper functions to project heat flux to coupling boundary
-    if side == 'Dirichlet':
-        # To communicate the flux to the Neumann side we should not simply
-        # evaluate u_,i n_i as this is an unbounded term leading to suboptimal
-        # convergence. Instead we project ∀ v: ∫_Γ v flux = ∫_Γ v u_,i n_i and
-        # evaluate flux. While the right-hand-side contains the same unbounded
-        # term, we can use the strong identity du/dt - u_,ii = f to rewrite it
-        # to ∫_Ω [v (du/dt - f) + v_,i u_,i] - ∫_∂Ω\Γ v u_,k n_k, in which we
-        # recognize the residual and an integral over the exterior boundary.
-        # While the latter still contains the problematic unbounded term, we
-        # can use the fact that the flux is a known value at the top and bottom
-        # via the Dirichlet boundary condition, and impose it as constraints.
-        rightsqr = domain.boundary['right'].integral(
-            'flux^2 d:x' @ ns, degree=degree * 2)
-        rightcons = solver.optimize('fluxdofs', rightsqr, droptol=1e-10)
-        # rightcons is NaN in dofs that are NOT supported on the right boundary
-        fluxsqr = domain.boundary['right'].boundary['top,bottom'].integral(
-            '(flux - uexact_,0)^2 d:x' @ ns, degree=degree * 2)
-        fluxcons = solver.optimize('fluxdofs',
-                                   fluxsqr,
-                                   droptol=1e-10,
-                                   constrain=np.choose(np.isnan(rightcons),
-                                                       [np.nan,
-                                                        0.]))
-        # fluxcons is NaN in dofs that are supported on ONLY the right boundary
-        fluxres = coupling_sample.integral('basis_n flux d:x' @ ns) - res
+    # To communicate the flux to the Neumann side we should not simply
+    # evaluate u_,i n_i as this is an unbounded term leading to suboptimal
+    # convergence. Instead we project ∀ v: ∫_Γ v flux = ∫_Γ v u_,i n_i and
+    # evaluate flux. While the right-hand-side contains the same unbounded
+    # term, we can use the strong identity du/dt - u_,ii = f to rewrite it
+    # to ∫_Ω [v (du/dt - f) + v_,i u_,i] - ∫_∂Ω\Γ v u_,k n_k, in which we
+    # recognize the residual and an integral over the exterior boundary.
+    # While the latter still contains the problematic unbounded term, we
+    # can use the fact that the flux is a known value at the top and bottom
+    # via the Dirichlet boundary condition, and impose it as constraints.
+    rightsqr = domain.boundary['right'].integral(
+        'flux^2 d:x' @ ns, degree=degree * 2)
+    rightcons = solver.optimize('fluxdofs', rightsqr, droptol=1e-10)
+    # rightcons is NaN in dofs that are NOT supported on the right boundary
+    fluxsqr = domain.boundary['right'].boundary['top,bottom'].integral(
+        '(flux - uexact_,0)^2 d:x' @ ns, degree=degree * 2)
+    fluxcons = solver.optimize('fluxdofs',
+                               fluxsqr,
+                               droptol=1e-10,
+                               constrain=np.choose(np.isnan(rightcons),
+                                                   [np.nan,
+                                                    0.]))
+    # fluxcons is NaN in dofs that are supported on ONLY the right boundary
+    fluxres = coupling_sample.integral('basis_n flux d:x' @ ns) - res
 
     # write initial data
     if participant.requires_initial_data():
@@ -118,8 +107,7 @@ def main(side='Dirichlet', n=10, degree=1, timestep=.1, alpha=3., beta=1.2):
         x, u, uexact = bezier.eval(['x_i', 'u', 'uexact'] @ ns, lhs=lhs, t=t)
         with treelog.add(treelog.DataLog()):
             export.vtk(
-                side +
-                "-" +
+                "Dirichlet-" +
                 str(istep),
                 bezier.tri,
                 x,
@@ -157,14 +145,12 @@ def main(side='Dirichlet', n=10, degree=1, timestep=.1, alpha=3., beta=1.2):
                 lhs0=lhs0, dt=dt, t=t, readdata=readdata))
 
         # write data to participant
-        if side == 'Dirichlet':
-            fluxdofs = solver.solve_linear(
-                'fluxdofs', fluxres, arguments=dict(
-                    lhs0=lhs0, lhs=lhs, dt=dt, t=t), constrain=fluxcons)
-            write_data = coupling_sample.eval(
-                'flux' @ ns, fluxdofs=fluxdofs)
-        else:
-            write_data = coupling_sample.eval('u' @ ns, lhs=lhs)
+        fluxdofs = solver.solve_linear(
+            'fluxdofs', fluxres, arguments=dict(
+                lhs0=lhs0, lhs=lhs, dt=dt, t=t), constrain=fluxcons)
+        write_data = coupling_sample.eval(
+            'flux' @ ns, fluxdofs=fluxdofs)
+
         precice_write(write_data)
 
         # do the coupling
diff --git a/partitioned-heat-conduction/nutils/requirements.txt b/partitioned-heat-conduction/dirichlet-nutils/requirements.txt
similarity index 100%
rename from partitioned-heat-conduction/nutils/requirements.txt
rename to partitioned-heat-conduction/dirichlet-nutils/requirements.txt
diff --git a/partitioned-heat-conduction/dirichlet-nutils/run.sh b/partitioned-heat-conduction/dirichlet-nutils/run.sh
new file mode 100755
index 000000000..e5743798a
--- /dev/null
+++ b/partitioned-heat-conduction/dirichlet-nutils/run.sh
@@ -0,0 +1,14 @@
+#!/bin/bash
+set -e -u
+
+. ../../tools/log.sh
+exec > >(tee --append "$LOGFILE") 2>&1
+
+python3 -m venv .venv
+. .venv/bin/activate
+pip install -r requirements.txt
+
+rm -rf Dirichlet-*.vtk
+NUTILS_RICHOUTPUT=no python3 heat.py
+
+close_log
diff --git a/partitioned-heat-conduction/neumann-nutils/clean.sh b/partitioned-heat-conduction/neumann-nutils/clean.sh
new file mode 100755
index 000000000..6893c1ea5
--- /dev/null
+++ b/partitioned-heat-conduction/neumann-nutils/clean.sh
@@ -0,0 +1,6 @@
+#!/bin/sh
+set -e -u
+
+. ../../tools/cleaning-tools.sh
+
+clean_nutils .
diff --git a/partitioned-heat-conduction/neumann-nutils/heat.py b/partitioned-heat-conduction/neumann-nutils/heat.py
new file mode 100644
index 000000000..4f4f06774
--- /dev/null
+++ b/partitioned-heat-conduction/neumann-nutils/heat.py
@@ -0,0 +1,136 @@
+#! /usr/bin/env python3
+
+from nutils import cli, mesh, function, solver, export
+import functools
+import treelog
+import numpy as np
+import precice
+
+
+def main(n=10, degree=1, timestep=.1, alpha=3., beta=1.2):
+
+    x_grid = np.linspace(1, 2, n)
+    y_grid = np.linspace(0, 1, n)
+
+    # define the Nutils mesh
+    domain, geom = mesh.rectilinear([x_grid, y_grid])
+    coupling_boundary = domain.boundary['left']
+    coupling_sample = coupling_boundary.sample('gauss', degree=degree * 2)
+
+    # Nutils namespace
+    ns = function.Namespace()
+    ns.x = geom
+    ns.basis = domain.basis('std', degree=degree)
+    ns.alpha = alpha  # parameter of problem
+    ns.beta = beta  # parameter of problem
+    ns.u = 'basis_n ?lhs_n'  # solution
+    ns.dudt = 'basis_n (?lhs_n - ?lhs0_n) / ?dt'  # time derivative
+    ns.flux = 'basis_n ?fluxdofs_n'  # heat flux
+    ns.f = 'beta - 2 - 2 alpha'  # rhs
+    ns.uexact = '1 + x_0 x_0 + alpha x_1 x_1 + beta ?t'  # analytical solution
+    ns.readbasis = coupling_sample.basis()
+    ns.readfunc = 'readbasis_n ?readdata_n'
+
+    # define the weak form
+    res = domain.integral(
+        '(basis_n dudt - basis_n f + basis_n,i u_,i) d:x' @ ns,
+        degree=degree * 2)
+
+    # set boundary conditions at non-coupling boundaries
+    # top and bottom boundary are non-coupling for both sides
+    sqr = domain.boundary['top,bottom,right'].integral('(u - uexact)^2 d:x' @ ns, degree=degree * 2)
+
+    res += coupling_sample.integral('basis_n readfunc d:x' @ ns)
+
+    # preCICE setup
+    participant = precice.Participant("Neumann", "../precice-config.xml", 0, 1)
+    mesh_name = "Neumann-Mesh"
+    vertex_ids = participant.set_mesh_vertices(
+        mesh_name, coupling_sample.eval(ns.x))
+    precice_write = functools.partial(
+        participant.write_data,
+        mesh_name,
+        "Temperature",
+        vertex_ids)
+    precice_read = functools.partial(
+        participant.read_data,
+        mesh_name,
+        "Heat-Flux",
+        vertex_ids)
+
+    # write initial data
+    if participant.requires_initial_data():
+        precice_write(coupling_sample.eval(0.))
+
+    participant.initialize()
+    precice_dt = participant.get_max_time_step_size()
+    solver_dt = timestep
+    dt = min(precice_dt, solver_dt)
+
+    t = 0.
+    istep = 0
+
+    # initial condition
+    sqr0 = domain.integral('(u - uexact)^2' @ ns, degree=degree * 2)
+    lhs = solver.optimize('lhs', sqr0, arguments=dict(t=t))
+    bezier = domain.sample('bezier', degree * 2)
+
+    while True:
+
+        # generate output
+        x, u, uexact = bezier.eval(['x_i', 'u', 'uexact'] @ ns, lhs=lhs, t=t)
+        with treelog.add(treelog.DataLog()):
+            export.vtk(
+                "Neumann-" +
+                str(istep),
+                bezier.tri,
+                x,
+                Temperature=u,
+                reference=uexact)
+
+        if not participant.is_coupling_ongoing():
+            break
+
+        # save checkpoint
+        if participant.requires_writing_checkpoint():
+            checkpoint = lhs, t, istep
+
+        # prepare next timestep
+        precice_dt = participant.get_max_time_step_size()
+        dt = min(solver_dt, precice_dt)
+        lhs0 = lhs
+        istep += 1
+        # read data from participant
+        readdata = precice_read(dt)
+        t += dt
+
+        # update (time-dependent) boundary condition
+        cons = solver.optimize(
+            'lhs',
+            sqr,
+            droptol=1e-15,
+            arguments=dict(
+                t=t,
+                readdata=readdata))
+
+        # solve nutils timestep
+        lhs = solver.solve_linear(
+            'lhs', res, constrain=cons, arguments=dict(
+                lhs0=lhs0, dt=dt, t=t, readdata=readdata))
+
+        # write data to participant
+        write_data = coupling_sample.eval('u' @ ns, lhs=lhs)
+        precice_write(write_data)
+
+        # do the coupling
+        participant.advance(dt)
+
+        # read checkpoint if required
+        if participant.requires_reading_checkpoint():
+            lhs, t, istep = checkpoint
+
+    participant.finalize()
+
+
+if __name__ == '__main__':
+    cli.run(main)
diff --git a/partitioned-heat-conduction/neumann-nutils/requirements.txt b/partitioned-heat-conduction/neumann-nutils/requirements.txt
new file mode 100644
index 000000000..cfd521c81
--- /dev/null
+++ b/partitioned-heat-conduction/neumann-nutils/requirements.txt
@@ -0,0 +1,2 @@
+nutils==7
+pyprecice==3
diff --git a/partitioned-heat-conduction/neumann-nutils/run.sh b/partitioned-heat-conduction/neumann-nutils/run.sh
new file mode 100755
index 000000000..021225a41
--- /dev/null
+++ b/partitioned-heat-conduction/neumann-nutils/run.sh
@@ -0,0 +1,14 @@
+#!/bin/bash
+set -e -u
+
+. ../../tools/log.sh
+exec > >(tee --append "$LOGFILE") 2>&1
+
+python3 -m venv .venv
+. .venv/bin/activate
+pip install -r requirements.txt
+
+rm -rf Neumann-*.vtk
+NUTILS_RICHOUTPUT=no python3 heat.py
+
+close_log
diff --git a/partitioned-heat-conduction/nutils/run.sh b/partitioned-heat-conduction/nutils/run.sh
deleted file mode 100755
index 11bcd5318..000000000
--- a/partitioned-heat-conduction/nutils/run.sh
+++ /dev/null
@@ -1,29 +0,0 @@
-#!/bin/bash
-set -e -u
-
-usage() { echo "Usage: cmd [-d] [-n]" 1>&2; exit 1; }
-
-# Check if no input argument was provided
-if [ -z "$*" ] ; then
-        usage
-fi
-
-python3 -m venv .venv
-. .venv/bin/activate
-pip install -r requirements.txt
-
-while getopts ":dn" opt; do
-  case ${opt} in
-  d)
-    rm -rf Dirichlet-*.vtk
-    NUTILS_RICHOUTPUT=no python3 heat.py side=Dirichlet
-    ;;
-  n)
-    rm -rf Neumann-*.vtk
-    NUTILS_RICHOUTPUT=no python3 heat.py side=Neumann
-    ;;
-  *)
-    usage
-    ;;
-  esac
-done