From 0779e3bbcfdb966f4293fd7908946e55f36dae51 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Tue, 2 Jan 2024 16:33:06 +0100
Subject: [PATCH 01/25] save checkpoint data in single variable

---
 perpendicular-flap/fluid-nutils/fluid.py | 16 ++--------------
 1 file changed, 2 insertions(+), 14 deletions(-)

diff --git a/perpendicular-flap/fluid-nutils/fluid.py b/perpendicular-flap/fluid-nutils/fluid.py
index e0e6b15fc..9cf3a6c0e 100644
--- a/perpendicular-flap/fluid-nutils/fluid.py
+++ b/perpendicular-flap/fluid-nutils/fluid.py
@@ -174,13 +174,7 @@ def main(inflow: 'inflow velocity' = 10,
 
         # save checkpoint
         if interface.is_action_required(precice.action_write_iteration_checkpoint()):
-            lhs_checkpoint = lhs0
-            lhs00_checkpoint = lhs00
-            t_checkpoint = t
-            timestep_checkpoint = timestep
-            oldmeshdofs_checkpoint = oldmeshdofs
-            oldoldmeshdofs_checkpoint = oldoldmeshdofs
-            oldoldoldmeshdofs_checkpoint = oldoldoldmeshdofs
+            checkpoint = lhs0, lhs00, t, timestep, oldmeshdofs, oldoldmeshdofs, oldoldoldmeshdofs
             interface.mark_action_fulfilled(precice.action_write_iteration_checkpoint())
 
         # solve fluid equations
@@ -215,13 +209,7 @@ def main(inflow: 'inflow velocity' = 10,
 
         # read checkpoint if required
         if interface.is_action_required(precice.action_read_iteration_checkpoint()):
-            lhs0 = lhs_checkpoint
-            lhs00 = lhs00_checkpoint
-            t = t_checkpoint
-            timestep = timestep_checkpoint
-            oldmeshdofs = oldmeshdofs_checkpoint
-            oldoldmeshdofs = oldoldmeshdofs_checkpoint
-            oldoldoldmeshdofs = oldoldoldmeshdofs_checkpoint
+            lhs0, lhs00, t, timestep, oldmeshdofs, oldoldmeshdofs, oldoldoldmeshdofs = checkpoint
             interface.mark_action_fulfilled(precice.action_read_iteration_checkpoint())
 
         if interface.is_time_window_complete():

From 86dc5f43173f356ca257a1d30263a5acca1a76ab Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Tue, 2 Jan 2024 16:36:51 +0100
Subject: [PATCH 02/25] replace old- prefix by -0 suffix

---
 perpendicular-flap/fluid-nutils/fluid.py | 42 ++++++++++++------------
 1 file changed, 21 insertions(+), 21 deletions(-)

diff --git a/perpendicular-flap/fluid-nutils/fluid.py b/perpendicular-flap/fluid-nutils/fluid.py
index 9cf3a6c0e..1a487230e 100644
--- a/perpendicular-flap/fluid-nutils/fluid.py
+++ b/perpendicular-flap/fluid-nutils/fluid.py
@@ -15,11 +15,11 @@ def subs0(f):
     if isinstance(f, function.Argument) and f._name == 'lhs':
         return function.Argument(name='lhs0', shape=f.shape, nderiv=f._nderiv)
     if isinstance(f, function.Argument) and f._name == 'meshdofs':
-        return function.Argument(name='oldmeshdofs', shape=f.shape, nderiv=f._nderiv)
-    if isinstance(f, function.Argument) and f._name == 'oldmeshdofs':
-        return function.Argument(name='oldoldmeshdofs', shape=f.shape, nderiv=f._nderiv)
-    if isinstance(f, function.Argument) and f._name == 'oldoldmeshdofs':
-        return function.Argument(name='oldoldoldmeshdofs', shape=f.shape, nderiv=f._nderiv)
+        return function.Argument(name='meshdofs0', shape=f.shape, nderiv=f._nderiv)
+    if isinstance(f, function.Argument) and f._name == 'meshdofs0':
+        return function.Argument(name='meshdofs00', shape=f.shape, nderiv=f._nderiv)
+    if isinstance(f, function.Argument) and f._name == 'meshdofs00':
+        return function.Argument(name='meshdofs000', shape=f.shape, nderiv=f._nderiv)
 
 # some helper function to shift variables by two timesteps
 
@@ -82,7 +82,7 @@ def main(inflow: 'inflow velocity' = 10,
     ns.x0 = geom  # reference geometry
     ns.dbasis = domain.basis('std', degree=1).vector(2)
     ns.d_i = 'dbasis_ni ?meshdofs_n'
-    ns.umesh_i = 'dbasis_ni (1.5 ?meshdofs_n - 2 ?oldmeshdofs_n + 0.5 ?oldoldmeshdofs_n ) / ?dt'
+    ns.umesh_i = 'dbasis_ni (1.5 ?meshdofs_n - 2 ?meshdofs0_n + 0.5 ?meshdofs00_n ) / ?dt'
     ns.x_i = 'x0_i + d_i'  # moving geometry
     ns.ubasis, ns.pbasis = function.chain([domain.basis('std', degree=2).vector(2), domain.basis('std', degree=1), ])
     ns.F_i = 'ubasis_ni ?F_n'  # stress field
@@ -92,9 +92,9 @@ def main(inflow: 'inflow velocity' = 10,
 
     # initialization of dofs
     meshdofs = numpy.zeros(len(ns.dbasis))
-    oldmeshdofs = meshdofs
-    oldoldmeshdofs = meshdofs
-    oldoldoldmeshdofs = meshdofs
+    meshdofs0 = meshdofs
+    meshdofs00 = meshdofs
+    meshdofs000 = meshdofs
     lhs0 = numpy.zeros(len(ns.ubasis))
 
     # for visualization
@@ -156,7 +156,7 @@ def main(inflow: 'inflow velocity' = 10,
 
     # better initial guess: start from Stokes solution, comment out for comparison with other solvers
     #res_stokes = domain.integral('(ubasis_ni,j ((u_i,j + u_j,i) rho nu - p δ_ij) + pbasis_n u_k,k) d:x' @ ns, degree=4)
-    #lhs0 = solver.solve_linear('lhs', res_stokes, constrain=cons, arguments=dict(meshdofs=meshdofs, oldmeshdofs=oldmeshdofs, oldoldmeshdofs=oldoldmeshdofs, oldoldoldmeshdofs=oldoldoldmeshdofs, dt=dt))
+    #lhs0 = solver.solve_linear('lhs', res_stokes, constrain=cons, arguments=dict(meshdofs=meshdofs, meshdofs0=meshdofs0, meshdofs00=meshdofs00, meshdofs000=meshdofs000, dt=dt))
     lhs00 = lhs0
 
     timestep = 0
@@ -174,21 +174,21 @@ def main(inflow: 'inflow velocity' = 10,
 
         # save checkpoint
         if interface.is_action_required(precice.action_write_iteration_checkpoint()):
-            checkpoint = lhs0, lhs00, t, timestep, oldmeshdofs, oldoldmeshdofs, oldoldoldmeshdofs
+            checkpoint = lhs0, lhs00, t, timestep, meshdofs0, meshdofs00, meshdofs000
             interface.mark_action_fulfilled(precice.action_write_iteration_checkpoint())
 
         # solve fluid equations
         lhs1 = solver.newton('lhs', res, lhs0=lhs0, constrain=cons,
-                             arguments=dict(lhs0=lhs0, dt=dt, meshdofs=meshdofs, oldmeshdofs=oldmeshdofs,
-                                            oldoldmeshdofs=oldoldmeshdofs, oldoldoldmeshdofs=oldoldoldmeshdofs)
+                             arguments=dict(lhs0=lhs0, dt=dt, meshdofs=meshdofs, meshdofs0=meshdofs0,
+                                            meshdofs00=meshdofs00, meshdofs000=meshdofs000)
                              ).solve(tol=1e-6)
 
         # write forces to interface
         if interface.is_write_data_required(dt):
             F = solver.solve_linear('F', resF, constrain=consF,
                                     arguments=dict(lhs00=lhs00, lhs0=lhs0, lhs=lhs1, dt=dt, meshdofs=meshdofs,
-                                                   oldmeshdofs=oldmeshdofs, oldoldmeshdofs=oldoldmeshdofs,
-                                                   oldoldoldmeshdofs=oldoldoldmeshdofs))
+                                                   meshdofs0=meshdofs0, meshdofs00=meshdofs00,
+                                                   meshdofs000=meshdofs000))
             # writedata = couplingsample.eval(ns.F, F=F) # for stresses
             writedata = couplingsample.eval('F_i d:x' @ ns, F=F, meshdofs=meshdofs) * \
                 numpy.concatenate([p.weights for p in couplingsample.points])[:, numpy.newaxis]
@@ -203,18 +203,18 @@ def main(inflow: 'inflow velocity' = 10,
         t += dt
         lhs00 = lhs0
         lhs0 = lhs1
-        oldoldoldmeshdofs = oldoldmeshdofs
-        oldoldmeshdofs = oldmeshdofs
-        oldmeshdofs = meshdofs
+        meshdofs000 = meshdofs00
+        meshdofs00 = meshdofs0
+        meshdofs0 = meshdofs
 
         # read checkpoint if required
         if interface.is_action_required(precice.action_read_iteration_checkpoint()):
-            lhs0, lhs00, t, timestep, oldmeshdofs, oldoldmeshdofs, oldoldoldmeshdofs = checkpoint
+            lhs0, lhs00, t, timestep, meshdofs0, meshdofs00, meshdofs000 = checkpoint
             interface.mark_action_fulfilled(precice.action_read_iteration_checkpoint())
 
         if interface.is_time_window_complete():
-            x, u, p = bezier.eval(['x_i', 'u_i', 'p'] @ ns, lhs=lhs1, meshdofs=meshdofs, oldmeshdofs=oldmeshdofs,
-                                  oldoldmeshdofs=oldoldmeshdofs, oldoldoldmeshdofs=oldoldoldmeshdofs, dt=dt)
+            x, u, p = bezier.eval(['x_i', 'u_i', 'p'] @ ns, lhs=lhs1, meshdofs=meshdofs, meshdofs0=meshdofs0,
+                                  meshdofs00=meshdofs00, meshdofs000=meshdofs000, dt=dt)
             with treelog.add(treelog.DataLog()):
                 export.vtk('Fluid_' + str(timestep), bezier.tri, x, u=u, p=p)
 

From 13c8a8bd6e7ab4ddb43c8d9e61c87781a911566f Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Mon, 8 Jan 2024 14:34:25 +0100
Subject: [PATCH 03/25] remove time approxmation function from namespace

---
 perpendicular-flap/fluid-nutils/fluid.py | 27 ++++++++++--------------
 1 file changed, 11 insertions(+), 16 deletions(-)

diff --git a/perpendicular-flap/fluid-nutils/fluid.py b/perpendicular-flap/fluid-nutils/fluid.py
index 1a487230e..3786028ec 100644
--- a/perpendicular-flap/fluid-nutils/fluid.py
+++ b/perpendicular-flap/fluid-nutils/fluid.py
@@ -59,23 +59,18 @@ def main(inflow: 'inflow velocity' = 10,
     domain = topo.withboundary(inflow='left', wall='top,bottom', outflow='right') - \
         topo[18:20, :10].withboundary(flap='left,right,top')
 
-    # Nutils namespace
-    ns = function.Namespace()
-
     # time approximations
     # TR interpolation
-    ns._functions['t'] = lambda f: theta * f + (1 - theta) * subs0(f)
-    ns._functions_nargs['t'] = 1
+    tθ = lambda f: theta * f + (1 - theta) * subs0(f)
     # 1st order FD
-    ns._functions['δt'] = lambda f: (f - subs0(f)) / dt
-    ns._functions_nargs['δt'] = 1
+    δt = lambda f: (f - subs0(f)) / dt
     # 2nd order FD
-    ns._functions['tt'] = lambda f: (1.5 * f - 2 * subs0(f) + 0.5 * subs00(f)) / dt
-    ns._functions_nargs['tt'] = 1
+    tt = lambda f: (1.5 * f - 2 * subs0(f) + 0.5 * subs00(f)) / dt
     # extrapolation for pressure
-    ns._functions['tp'] = lambda f: (1.5 * f - 0.5 * subs0(f))
-    ns._functions_nargs['tp'] = 1
+    tp = lambda f: (1.5 * f - 0.5 * subs0(f))
 
+    # Nutils namespace
+    ns = function.Namespace()
     ns.nu = viscosity
     ns.rho = density
     ns.uin = inflow
@@ -132,16 +127,16 @@ def main(inflow: 'inflow velocity' = 10,
     cons = solver.optimize('lhs', sqr, droptol=1e-15, constrain=cons)
 
     # weak form fluid equations
-    res = domain.integral('t(ubasis_ni,j (u_i,j + u_j,i) rho nu d:x)' @ ns, degree=4)
+    res = tθ(domain.integral('ubasis_ni,j (u_i,j + u_j,i) rho nu d:x' @ ns, degree=4))
     res += domain.integral('(-ubasis_ni,j p δ_ij + pbasis_n u_k,k) d:x' @ ns, degree=4)
-    res += domain.integral('rho ubasis_ni δt(u_i d:x)' @ ns, degree=4)
-    res += domain.integral('rho ubasis_ni t(u_i,j urel_j d:x)' @ ns, degree=4)
+    res += δt(domain.integral('rho ubasis_ni u_i d:x' @ ns, degree=4))
+    res += tθ(domain.integral('rho ubasis_ni u_i,j urel_j d:x' @ ns, degree=4))
 
     # weak form for force computation
     resF = domain.integral('(ubasis_ni,j (u_i,j + u_j,i) rho nu d:x)' @ ns, degree=4)
-    resF += domain.integral('tp(-ubasis_ni,j p δ_ij d:x)' @ ns, degree=4)
+    resF += tp(domain.integral('-ubasis_ni,j p δ_ij d:x' @ ns, degree=4))
     resF += domain.integral('pbasis_n u_k,k d:x' @ ns, degree=4)
-    resF += domain.integral('rho ubasis_ni tt(u_i d:x)' @ ns, degree=4)
+    resF += tt(domain.integral('rho ubasis_ni u_i d:x' @ ns, degree=4))
     resF += domain.integral('rho ubasis_ni (u_i,j urel_j d:x)' @ ns, degree=4)
     resF += couplinginterface.sample('gauss', 4).integral('ubasis_ni F_i d:x' @ ns)
     consF = numpy.isnan(solver.optimize('F', couplinginterface.sample('gauss', 4).integral('F_i F_i' @ ns),

From c57ee0c65f7f7c6dc909cd6498174f8fe657d565 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Tue, 2 Jan 2024 16:37:35 +0100
Subject: [PATCH 04/25] remove subs00 helper function

---
 perpendicular-flap/fluid-nutils/fluid.py | 10 +---------
 1 file changed, 1 insertion(+), 9 deletions(-)

diff --git a/perpendicular-flap/fluid-nutils/fluid.py b/perpendicular-flap/fluid-nutils/fluid.py
index 3786028ec..3da8002da 100644
--- a/perpendicular-flap/fluid-nutils/fluid.py
+++ b/perpendicular-flap/fluid-nutils/fluid.py
@@ -21,14 +21,6 @@ def subs0(f):
     if isinstance(f, function.Argument) and f._name == 'meshdofs00':
         return function.Argument(name='meshdofs000', shape=f.shape, nderiv=f._nderiv)
 
-# some helper function to shift variables by two timesteps
-
-
-@function.replace
-def subs00(f):
-    if isinstance(f, function.Argument) and f._name == 'lhs':
-        return function.Argument(name='lhs00', shape=f.shape, nderiv=f._nderiv)
-
 
 def main(inflow: 'inflow velocity' = 10,
          viscosity: 'kinematic viscosity' = 1.0,
@@ -65,7 +57,7 @@ def main(inflow: 'inflow velocity' = 10,
     # 1st order FD
     δt = lambda f: (f - subs0(f)) / dt
     # 2nd order FD
-    tt = lambda f: (1.5 * f - 2 * subs0(f) + 0.5 * subs00(f)) / dt
+    tt = lambda f: (1.5 * f - 2 * subs0(f) + 0.5 * subs0(subs0(f))) / dt
     # extrapolation for pressure
     tp = lambda f: (1.5 * f - 0.5 * subs0(f))
 

From 6e05ed0471ce30c7f1622950e44af0d1e6ff7345 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Tue, 2 Jan 2024 16:41:54 +0100
Subject: [PATCH 05/25] simplify subs0 function

---
 perpendicular-flap/fluid-nutils/fluid.py | 26 ++++++++----------------
 1 file changed, 9 insertions(+), 17 deletions(-)

diff --git a/perpendicular-flap/fluid-nutils/fluid.py b/perpendicular-flap/fluid-nutils/fluid.py
index 3da8002da..0906fd29e 100644
--- a/perpendicular-flap/fluid-nutils/fluid.py
+++ b/perpendicular-flap/fluid-nutils/fluid.py
@@ -7,19 +7,9 @@
 
 # for details on this solver see https://doi.org/10.1002/nme.6443
 
-# some helper function to shift variables by one timestep
-
-
-@function.replace
 def subs0(f):
-    if isinstance(f, function.Argument) and f._name == 'lhs':
-        return function.Argument(name='lhs0', shape=f.shape, nderiv=f._nderiv)
-    if isinstance(f, function.Argument) and f._name == 'meshdofs':
-        return function.Argument(name='meshdofs0', shape=f.shape, nderiv=f._nderiv)
-    if isinstance(f, function.Argument) and f._name == 'meshdofs0':
-        return function.Argument(name='meshdofs00', shape=f.shape, nderiv=f._nderiv)
-    if isinstance(f, function.Argument) and f._name == 'meshdofs00':
-        return function.Argument(name='meshdofs000', shape=f.shape, nderiv=f._nderiv)
+    'helper function to shift variables by one timestep'
+    return function.replace_arguments(f, {arg: function.Argument(arg+'0', shape=shape, dtype=dtype) for arg, (shape, dtype) in f.arguments.items() if arg != 'dt'})
 
 
 def main(inflow: 'inflow velocity' = 10,
@@ -82,6 +72,7 @@ def main(inflow: 'inflow velocity' = 10,
     meshdofs0 = meshdofs
     meshdofs00 = meshdofs
     meshdofs000 = meshdofs
+    meshdofs0000 = meshdofs
     lhs0 = numpy.zeros(len(ns.ubasis))
 
     # for visualization
@@ -161,13 +152,13 @@ def main(inflow: 'inflow velocity' = 10,
 
         # save checkpoint
         if interface.is_action_required(precice.action_write_iteration_checkpoint()):
-            checkpoint = lhs0, lhs00, t, timestep, meshdofs0, meshdofs00, meshdofs000
+            checkpoint = lhs0, lhs00, t, timestep, meshdofs0, meshdofs00, meshdofs000, meshdofs0000
             interface.mark_action_fulfilled(precice.action_write_iteration_checkpoint())
 
         # solve fluid equations
         lhs1 = solver.newton('lhs', res, lhs0=lhs0, constrain=cons,
                              arguments=dict(lhs0=lhs0, dt=dt, meshdofs=meshdofs, meshdofs0=meshdofs0,
-                                            meshdofs00=meshdofs00, meshdofs000=meshdofs000)
+                                            meshdofs00=meshdofs00, meshdofs000=meshdofs000, meshdofs0000=meshdofs0000)
                              ).solve(tol=1e-6)
 
         # write forces to interface
@@ -175,7 +166,7 @@ def main(inflow: 'inflow velocity' = 10,
             F = solver.solve_linear('F', resF, constrain=consF,
                                     arguments=dict(lhs00=lhs00, lhs0=lhs0, lhs=lhs1, dt=dt, meshdofs=meshdofs,
                                                    meshdofs0=meshdofs0, meshdofs00=meshdofs00,
-                                                   meshdofs000=meshdofs000))
+                                                   meshdofs000=meshdofs000, meshdofs0000=meshdofs0000))
             # writedata = couplingsample.eval(ns.F, F=F) # for stresses
             writedata = couplingsample.eval('F_i d:x' @ ns, F=F, meshdofs=meshdofs) * \
                 numpy.concatenate([p.weights for p in couplingsample.points])[:, numpy.newaxis]
@@ -190,18 +181,19 @@ def main(inflow: 'inflow velocity' = 10,
         t += dt
         lhs00 = lhs0
         lhs0 = lhs1
+        meshdofs0000 = meshdofs000
         meshdofs000 = meshdofs00
         meshdofs00 = meshdofs0
         meshdofs0 = meshdofs
 
         # read checkpoint if required
         if interface.is_action_required(precice.action_read_iteration_checkpoint()):
-            lhs0, lhs00, t, timestep, meshdofs0, meshdofs00, meshdofs000 = checkpoint
+            lhs0, lhs00, t, timestep, meshdofs0, meshdofs00, meshdofs000, meshdofs0000 = checkpoint
             interface.mark_action_fulfilled(precice.action_read_iteration_checkpoint())
 
         if interface.is_time_window_complete():
             x, u, p = bezier.eval(['x_i', 'u_i', 'p'] @ ns, lhs=lhs1, meshdofs=meshdofs, meshdofs0=meshdofs0,
-                                  meshdofs00=meshdofs00, meshdofs000=meshdofs000, dt=dt)
+                                  meshdofs00=meshdofs00, meshdofs000=meshdofs000, meshdofs0000=meshdofs0000, dt=dt)
             with treelog.add(treelog.DataLog()):
                 export.vtk('Fluid_' + str(timestep), bezier.tri, x, u=u, p=p)
 

From 2a0b465e70ecc4555e5cd9396acb6e7ece66660d Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Mon, 8 Jan 2024 15:16:52 +0100
Subject: [PATCH 06/25] rename subs0 to t0 and inline

---
 perpendicular-flap/fluid-nutils/fluid.py | 15 ++++++---------
 1 file changed, 6 insertions(+), 9 deletions(-)

diff --git a/perpendicular-flap/fluid-nutils/fluid.py b/perpendicular-flap/fluid-nutils/fluid.py
index 0906fd29e..7249aa831 100644
--- a/perpendicular-flap/fluid-nutils/fluid.py
+++ b/perpendicular-flap/fluid-nutils/fluid.py
@@ -7,11 +7,6 @@
 
 # for details on this solver see https://doi.org/10.1002/nme.6443
 
-def subs0(f):
-    'helper function to shift variables by one timestep'
-    return function.replace_arguments(f, {arg: function.Argument(arg+'0', shape=shape, dtype=dtype) for arg, (shape, dtype) in f.arguments.items() if arg != 'dt'})
-
-
 def main(inflow: 'inflow velocity' = 10,
          viscosity: 'kinematic viscosity' = 1.0,
          density: 'density' = 1.0,
@@ -42,14 +37,16 @@ def main(inflow: 'inflow velocity' = 10,
         topo[18:20, :10].withboundary(flap='left,right,top')
 
     # time approximations
+    t0 = lambda f: function.replace_arguments(f, {arg: function.Argument(arg+'0', shape=shape, dtype=dtype)
+                        for arg, (shape, dtype) in f.arguments.items() if arg != 'dt'})
     # TR interpolation
-    tθ = lambda f: theta * f + (1 - theta) * subs0(f)
+    tθ = lambda f: theta * f + (1 - theta) * t0(f)
     # 1st order FD
-    δt = lambda f: (f - subs0(f)) / dt
+    δt = lambda f: (f - t0(f)) / dt
     # 2nd order FD
-    tt = lambda f: (1.5 * f - 2 * subs0(f) + 0.5 * subs0(subs0(f))) / dt
+    tt = lambda f: (1.5 * f - 2 * t0(f) + 0.5 * t0(t0(f))) / dt
     # extrapolation for pressure
-    tp = lambda f: (1.5 * f - 0.5 * subs0(f))
+    tp = lambda f: (1.5 * f - 0.5 * t0(f))
 
     # Nutils namespace
     ns = function.Namespace()

From fca2acc993ba74396d292e2f16b1a731093aff83 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Mon, 8 Jan 2024 15:55:28 +0100
Subject: [PATCH 07/25] add treelog context for every time step

---
 perpendicular-flap/fluid-nutils/fluid.py | 1 +
 perpendicular-flap/fluid-nutils/run.sh   | 2 +-
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/perpendicular-flap/fluid-nutils/fluid.py b/perpendicular-flap/fluid-nutils/fluid.py
index 7249aa831..f566ed5bd 100644
--- a/perpendicular-flap/fluid-nutils/fluid.py
+++ b/perpendicular-flap/fluid-nutils/fluid.py
@@ -138,6 +138,7 @@ def main(inflow: 'inflow velocity' = 10,
     t = 0
 
     while interface.is_coupling_ongoing():
+      with treelog.context(f'timestep {timestep}'):
 
         # read displacements from interface
         if interface.is_read_data_available():
diff --git a/perpendicular-flap/fluid-nutils/run.sh b/perpendicular-flap/fluid-nutils/run.sh
index 788f40d13..2ebf6cbf7 100755
--- a/perpendicular-flap/fluid-nutils/run.sh
+++ b/perpendicular-flap/fluid-nutils/run.sh
@@ -1,4 +1,4 @@
 #!/bin/sh
 set -e -u
 
-python3 fluid.py
+python3 fluid.py richoutput=no

From 442287cafa6701678fa9c9b312cdfe1a552f3054 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Wed, 3 Jan 2024 11:28:30 +0100
Subject: [PATCH 08/25] add matplotlib figures

---
 perpendicular-flap/fluid-nutils/fluid.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/perpendicular-flap/fluid-nutils/fluid.py b/perpendicular-flap/fluid-nutils/fluid.py
index f566ed5bd..703c91e6a 100644
--- a/perpendicular-flap/fluid-nutils/fluid.py
+++ b/perpendicular-flap/fluid-nutils/fluid.py
@@ -192,6 +192,8 @@ def main(inflow: 'inflow velocity' = 10,
         if interface.is_time_window_complete():
             x, u, p = bezier.eval(['x_i', 'u_i', 'p'] @ ns, lhs=lhs1, meshdofs=meshdofs, meshdofs0=meshdofs0,
                                   meshdofs00=meshdofs00, meshdofs000=meshdofs000, meshdofs0000=meshdofs0000, dt=dt)
+            export.triplot('velocity.jpg', x, numpy.linalg.norm(u, axis=1), tri=bezier.tri, cmap='jet')
+            export.triplot('pressure.jpg', x, p, tri=bezier.tri, cmap='jet')
             with treelog.add(treelog.DataLog()):
                 export.vtk('Fluid_' + str(timestep), bezier.tri, x, u=u, p=p)
 

From 4cbf721fe5d223a8e79fd766c2a0c85c35331199 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Thu, 4 Jan 2024 14:37:47 +0100
Subject: [PATCH 09/25] refer to dt via argument

---
 perpendicular-flap/fluid-nutils/fluid.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/perpendicular-flap/fluid-nutils/fluid.py b/perpendicular-flap/fluid-nutils/fluid.py
index 703c91e6a..38b4ab991 100644
--- a/perpendicular-flap/fluid-nutils/fluid.py
+++ b/perpendicular-flap/fluid-nutils/fluid.py
@@ -42,9 +42,9 @@ def main(inflow: 'inflow velocity' = 10,
     # TR interpolation
     tθ = lambda f: theta * f + (1 - theta) * t0(f)
     # 1st order FD
-    δt = lambda f: (f - t0(f)) / dt
+    δt = lambda f: (f - t0(f)) / function.Argument('dt', ())
     # 2nd order FD
-    tt = lambda f: (1.5 * f - 2 * t0(f) + 0.5 * t0(t0(f))) / dt
+    tt = lambda f: (1.5 * f - 2 * t0(f) + 0.5 * t0(t0(f))) / function.Argument('dt', ())
     # extrapolation for pressure
     tp = lambda f: (1.5 * f - 0.5 * t0(f))
 

From 767c19af0b8ceadb772da9378c575ca32a311cb8 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Thu, 4 Jan 2024 14:40:24 +0100
Subject: [PATCH 10/25] remove unused t variable

---
 perpendicular-flap/fluid-nutils/fluid.py | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/perpendicular-flap/fluid-nutils/fluid.py b/perpendicular-flap/fluid-nutils/fluid.py
index 38b4ab991..d089633ae 100644
--- a/perpendicular-flap/fluid-nutils/fluid.py
+++ b/perpendicular-flap/fluid-nutils/fluid.py
@@ -135,7 +135,6 @@ def main(inflow: 'inflow velocity' = 10,
     lhs00 = lhs0
 
     timestep = 0
-    t = 0
 
     while interface.is_coupling_ongoing():
       with treelog.context(f'timestep {timestep}'):
@@ -150,7 +149,7 @@ def main(inflow: 'inflow velocity' = 10,
 
         # save checkpoint
         if interface.is_action_required(precice.action_write_iteration_checkpoint()):
-            checkpoint = lhs0, lhs00, t, timestep, meshdofs0, meshdofs00, meshdofs000, meshdofs0000
+            checkpoint = lhs0, lhs00, timestep, meshdofs0, meshdofs00, meshdofs000, meshdofs0000
             interface.mark_action_fulfilled(precice.action_write_iteration_checkpoint())
 
         # solve fluid equations
@@ -176,7 +175,6 @@ def main(inflow: 'inflow velocity' = 10,
 
         # advance variables
         timestep += 1
-        t += dt
         lhs00 = lhs0
         lhs0 = lhs1
         meshdofs0000 = meshdofs000
@@ -186,7 +184,7 @@ def main(inflow: 'inflow velocity' = 10,
 
         # read checkpoint if required
         if interface.is_action_required(precice.action_read_iteration_checkpoint()):
-            lhs0, lhs00, t, timestep, meshdofs0, meshdofs00, meshdofs000, meshdofs0000 = checkpoint
+            lhs0, lhs00, timestep, meshdofs0, meshdofs00, meshdofs000, meshdofs0000 = checkpoint
             interface.mark_action_fulfilled(precice.action_read_iteration_checkpoint())
 
         if interface.is_time_window_complete():

From 628a1e67dfa2442c4ae699a5fed77df7c7f38b3b Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Thu, 4 Jan 2024 14:40:43 +0100
Subject: [PATCH 11/25] rename lhs1 to lhs

---
 perpendicular-flap/fluid-nutils/fluid.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/perpendicular-flap/fluid-nutils/fluid.py b/perpendicular-flap/fluid-nutils/fluid.py
index d089633ae..64649ad28 100644
--- a/perpendicular-flap/fluid-nutils/fluid.py
+++ b/perpendicular-flap/fluid-nutils/fluid.py
@@ -153,7 +153,7 @@ def main(inflow: 'inflow velocity' = 10,
             interface.mark_action_fulfilled(precice.action_write_iteration_checkpoint())
 
         # solve fluid equations
-        lhs1 = solver.newton('lhs', res, lhs0=lhs0, constrain=cons,
+        lhs = solver.newton('lhs', res, lhs0=lhs0, constrain=cons,
                              arguments=dict(lhs0=lhs0, dt=dt, meshdofs=meshdofs, meshdofs0=meshdofs0,
                                             meshdofs00=meshdofs00, meshdofs000=meshdofs000, meshdofs0000=meshdofs0000)
                              ).solve(tol=1e-6)
@@ -161,7 +161,7 @@ def main(inflow: 'inflow velocity' = 10,
         # write forces to interface
         if interface.is_write_data_required(dt):
             F = solver.solve_linear('F', resF, constrain=consF,
-                                    arguments=dict(lhs00=lhs00, lhs0=lhs0, lhs=lhs1, dt=dt, meshdofs=meshdofs,
+                                    arguments=dict(lhs00=lhs00, lhs0=lhs0, lhs=lhs, dt=dt, meshdofs=meshdofs,
                                                    meshdofs0=meshdofs0, meshdofs00=meshdofs00,
                                                    meshdofs000=meshdofs000, meshdofs0000=meshdofs0000))
             # writedata = couplingsample.eval(ns.F, F=F) # for stresses
@@ -176,7 +176,7 @@ def main(inflow: 'inflow velocity' = 10,
         # advance variables
         timestep += 1
         lhs00 = lhs0
-        lhs0 = lhs1
+        lhs0 = lhs
         meshdofs0000 = meshdofs000
         meshdofs000 = meshdofs00
         meshdofs00 = meshdofs0
@@ -188,7 +188,7 @@ def main(inflow: 'inflow velocity' = 10,
             interface.mark_action_fulfilled(precice.action_read_iteration_checkpoint())
 
         if interface.is_time_window_complete():
-            x, u, p = bezier.eval(['x_i', 'u_i', 'p'] @ ns, lhs=lhs1, meshdofs=meshdofs, meshdofs0=meshdofs0,
+            x, u, p = bezier.eval(['x_i', 'u_i', 'p'] @ ns, lhs=lhs, meshdofs=meshdofs, meshdofs0=meshdofs0,
                                   meshdofs00=meshdofs00, meshdofs000=meshdofs000, meshdofs0000=meshdofs0000, dt=dt)
             export.triplot('velocity.jpg', x, numpy.linalg.norm(u, axis=1), tri=bezier.tri, cmap='jet')
             export.triplot('pressure.jpg', x, p, tri=bezier.tri, cmap='jet')

From 4b7b925e400d48a0c812d5bba7f2f537c60ea14c Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Thu, 4 Jan 2024 14:52:00 +0100
Subject: [PATCH 12/25] move history update to before fluid solve

Formerly, history values were updated prior to post-processing, so that all
evaluations happened with lhs==lhs0 etc. This patch fixes this, improving VTK
and matplotlib results slightly.
---
 perpendicular-flap/fluid-nutils/fluid.py | 30 +++++++++++-------------
 1 file changed, 14 insertions(+), 16 deletions(-)

diff --git a/perpendicular-flap/fluid-nutils/fluid.py b/perpendicular-flap/fluid-nutils/fluid.py
index 64649ad28..0109d98f2 100644
--- a/perpendicular-flap/fluid-nutils/fluid.py
+++ b/perpendicular-flap/fluid-nutils/fluid.py
@@ -69,8 +69,8 @@ def main(inflow: 'inflow velocity' = 10,
     meshdofs0 = meshdofs
     meshdofs00 = meshdofs
     meshdofs000 = meshdofs
-    meshdofs0000 = meshdofs
-    lhs0 = numpy.zeros(len(ns.ubasis))
+    lhs = numpy.zeros(len(ns.ubasis))
+    lhs0 = lhs
 
     # for visualization
     bezier = domain.sample('bezier', 2)
@@ -98,7 +98,6 @@ def main(inflow: 'inflow velocity' = 10,
 
     # initialize preCICE
     precice_dt = interface.initialize()
-    dt = min(precice_dt, timestepsize)
 
     # boundary conditions for fluid equations
     sqr = domain.boundary['wall,flap'].integral('urel_k urel_k d:x0' @ ns, degree=4)
@@ -132,7 +131,6 @@ def main(inflow: 'inflow velocity' = 10,
     # better initial guess: start from Stokes solution, comment out for comparison with other solvers
     #res_stokes = domain.integral('(ubasis_ni,j ((u_i,j + u_j,i) rho nu - p δ_ij) + pbasis_n u_k,k) d:x' @ ns, degree=4)
     #lhs0 = solver.solve_linear('lhs', res_stokes, constrain=cons, arguments=dict(meshdofs=meshdofs, meshdofs0=meshdofs0, meshdofs00=meshdofs00, meshdofs000=meshdofs000, dt=dt))
-    lhs00 = lhs0
 
     timestep = 0
 
@@ -149,9 +147,19 @@ def main(inflow: 'inflow velocity' = 10,
 
         # save checkpoint
         if interface.is_action_required(precice.action_write_iteration_checkpoint()):
-            checkpoint = lhs0, lhs00, timestep, meshdofs0, meshdofs00, meshdofs000, meshdofs0000
+            checkpoint = timestep, lhs, lhs0, meshdofs, meshdofs0, meshdofs00, meshdofs000
             interface.mark_action_fulfilled(precice.action_write_iteration_checkpoint())
 
+        # advance variables
+        timestep += 1
+        lhs00 = lhs0
+        lhs0 = lhs
+        meshdofs0000 = meshdofs000
+        meshdofs000 = meshdofs00
+        meshdofs00 = meshdofs0
+        meshdofs0 = meshdofs
+        dt = min(precice_dt, timestepsize)
+
         # solve fluid equations
         lhs = solver.newton('lhs', res, lhs0=lhs0, constrain=cons,
                              arguments=dict(lhs0=lhs0, dt=dt, meshdofs=meshdofs, meshdofs0=meshdofs0,
@@ -171,20 +179,10 @@ def main(inflow: 'inflow velocity' = 10,
 
         # do the coupling
         precice_dt = interface.advance(dt)
-        dt = min(precice_dt, timestepsize)
-
-        # advance variables
-        timestep += 1
-        lhs00 = lhs0
-        lhs0 = lhs
-        meshdofs0000 = meshdofs000
-        meshdofs000 = meshdofs00
-        meshdofs00 = meshdofs0
-        meshdofs0 = meshdofs
 
         # read checkpoint if required
         if interface.is_action_required(precice.action_read_iteration_checkpoint()):
-            lhs0, lhs00, timestep, meshdofs0, meshdofs00, meshdofs000, meshdofs0000 = checkpoint
+            timestep, lhs, lhs0, meshdofs, meshdofs0, meshdofs00, meshdofs000 = checkpoint
             interface.mark_action_fulfilled(precice.action_read_iteration_checkpoint())
 
         if interface.is_time_window_complete():

From 6c70f7e41532b7b659997b5a4d16b1e347986366 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Wed, 3 Jan 2024 11:54:58 +0100
Subject: [PATCH 13/25] compute meshdofs as part of physics update

This fixes the history of meshdofs, which was updated before pushing the
previous version to history.
---
 perpendicular-flap/fluid-nutils/fluid.py | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/perpendicular-flap/fluid-nutils/fluid.py b/perpendicular-flap/fluid-nutils/fluid.py
index 0109d98f2..592487842 100644
--- a/perpendicular-flap/fluid-nutils/fluid.py
+++ b/perpendicular-flap/fluid-nutils/fluid.py
@@ -124,6 +124,8 @@ def main(inflow: 'inflow velocity' = 10,
     # boundary conditions mesh displacements
     sqr = domain.boundary['inflow,outflow,wall'].integral('d_i d_i' @ ns, degree=2)
     meshcons0 = solver.optimize('meshdofs', sqr, droptol=1e-15)
+    sqr = couplingsample.integral('d_k d_k' @ ns)
+    meshcons = solver.optimize('meshdofs', sqr, droptol=1e-15, constrain=meshcons0)
 
     # weak form mesh displacements
     meshsqr = domain.integral('d_i,x0_j d_i,x0_j d:x0' @ ns, degree=2)
@@ -143,7 +145,6 @@ def main(inflow: 'inflow velocity' = 10,
             coupledata = couplingsample.asfunction(readdata)
             sqr = couplingsample.integral(((ns.d - coupledata)**2).sum(0))
             meshcons = solver.optimize('meshdofs', sqr, droptol=1e-15, constrain=meshcons0)
-            meshdofs = solver.optimize('meshdofs', meshsqr, constrain=meshcons)
 
         # save checkpoint
         if interface.is_action_required(precice.action_write_iteration_checkpoint()):
@@ -160,6 +161,9 @@ def main(inflow: 'inflow velocity' = 10,
         meshdofs0 = meshdofs
         dt = min(precice_dt, timestepsize)
 
+        # solve mesh deformation
+        meshdofs = solver.optimize('meshdofs', meshsqr, constrain=meshcons)
+
         # solve fluid equations
         lhs = solver.newton('lhs', res, lhs0=lhs0, constrain=cons,
                              arguments=dict(lhs0=lhs0, dt=dt, meshdofs=meshdofs, meshdofs0=meshdofs0,

From b55dea8df8a76bca53ee390eb1061cfce9df9c94 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Thu, 4 Jan 2024 19:37:20 +0100
Subject: [PATCH 14/25] replace umesh by tt(d)

---
 perpendicular-flap/fluid-nutils/fluid.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/perpendicular-flap/fluid-nutils/fluid.py b/perpendicular-flap/fluid-nutils/fluid.py
index 592487842..41a0d1cb3 100644
--- a/perpendicular-flap/fluid-nutils/fluid.py
+++ b/perpendicular-flap/fluid-nutils/fluid.py
@@ -56,7 +56,7 @@ def main(inflow: 'inflow velocity' = 10,
     ns.x0 = geom  # reference geometry
     ns.dbasis = domain.basis('std', degree=1).vector(2)
     ns.d_i = 'dbasis_ni ?meshdofs_n'
-    ns.umesh_i = 'dbasis_ni (1.5 ?meshdofs_n - 2 ?meshdofs0_n + 0.5 ?meshdofs00_n ) / ?dt'
+    ns.umesh = tt(ns.d) # mesh velocity
     ns.x_i = 'x0_i + d_i'  # moving geometry
     ns.ubasis, ns.pbasis = function.chain([domain.basis('std', degree=2).vector(2), domain.basis('std', degree=1), ])
     ns.F_i = 'ubasis_ni ?F_n'  # stress field

From fd09e72d0234ad16c4df2c3eb3013c438c605348 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Wed, 3 Jan 2024 20:23:37 +0100
Subject: [PATCH 15/25] =?UTF-8?q?redefine=20tt=20via=20tp=20and=20=CE=B4t?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 perpendicular-flap/fluid-nutils/fluid.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/perpendicular-flap/fluid-nutils/fluid.py b/perpendicular-flap/fluid-nutils/fluid.py
index 41a0d1cb3..1b0c5023f 100644
--- a/perpendicular-flap/fluid-nutils/fluid.py
+++ b/perpendicular-flap/fluid-nutils/fluid.py
@@ -43,10 +43,10 @@ def main(inflow: 'inflow velocity' = 10,
     tθ = lambda f: theta * f + (1 - theta) * t0(f)
     # 1st order FD
     δt = lambda f: (f - t0(f)) / function.Argument('dt', ())
-    # 2nd order FD
-    tt = lambda f: (1.5 * f - 2 * t0(f) + 0.5 * t0(t0(f))) / function.Argument('dt', ())
     # extrapolation for pressure
     tp = lambda f: (1.5 * f - 0.5 * t0(f))
+    # 2nd order FD
+    tt = lambda f: tp(δt(f))
 
     # Nutils namespace
     ns = function.Namespace()

From bbeac373d1ca29bf2d69e45c0e3551c4d43e4854 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Wed, 3 Jan 2024 21:49:38 +0100
Subject: [PATCH 16/25] introduce arguments, fix variable time steps

This improves results in case dt varies from timestep to timestep.
---
 perpendicular-flap/fluid-nutils/fluid.py | 50 +++++++-----------------
 1 file changed, 15 insertions(+), 35 deletions(-)

diff --git a/perpendicular-flap/fluid-nutils/fluid.py b/perpendicular-flap/fluid-nutils/fluid.py
index 1b0c5023f..be0c6a5c6 100644
--- a/perpendicular-flap/fluid-nutils/fluid.py
+++ b/perpendicular-flap/fluid-nutils/fluid.py
@@ -38,7 +38,7 @@ def main(inflow: 'inflow velocity' = 10,
 
     # time approximations
     t0 = lambda f: function.replace_arguments(f, {arg: function.Argument(arg+'0', shape=shape, dtype=dtype)
-                        for arg, (shape, dtype) in f.arguments.items() if arg != 'dt'})
+                        for arg, (shape, dtype) in f.arguments.items()})
     # TR interpolation
     tθ = lambda f: theta * f + (1 - theta) * t0(f)
     # 1st order FD
@@ -64,14 +64,6 @@ def main(inflow: 'inflow velocity' = 10,
     ns.u_i = 'umesh_i + urel_i'  # total velocity
     ns.p = 'pbasis_n ?lhs_n'  # pressure
 
-    # initialization of dofs
-    meshdofs = numpy.zeros(len(ns.dbasis))
-    meshdofs0 = meshdofs
-    meshdofs00 = meshdofs
-    meshdofs000 = meshdofs
-    lhs = numpy.zeros(len(ns.ubasis))
-    lhs0 = lhs
-
     # for visualization
     bezier = domain.sample('bezier', 2)
 
@@ -135,6 +127,10 @@ def main(inflow: 'inflow velocity' = 10,
     #lhs0 = solver.solve_linear('lhs', res_stokes, constrain=cons, arguments=dict(meshdofs=meshdofs, meshdofs0=meshdofs0, meshdofs00=meshdofs00, meshdofs000=meshdofs000, dt=dt))
 
     timestep = 0
+    arguments = dict(lhs=numpy.zeros(len(ns.ubasis)), meshdofs=numpy.zeros(len(ns.dbasis)), dt=timestepsize)
+
+    nhist = 4 # history length required by the formulation
+    arguments.update((k+'0'*i, v) for k, v in tuple(arguments.items()) for i in range(nhist))
 
     while interface.is_coupling_ongoing():
       with treelog.context(f'timestep {timestep}'):
@@ -148,36 +144,21 @@ def main(inflow: 'inflow velocity' = 10,
 
         # save checkpoint
         if interface.is_action_required(precice.action_write_iteration_checkpoint()):
-            checkpoint = timestep, lhs, lhs0, meshdofs, meshdofs0, meshdofs00, meshdofs000
+            checkpoint = timestep, arguments
             interface.mark_action_fulfilled(precice.action_write_iteration_checkpoint())
 
         # advance variables
         timestep += 1
-        lhs00 = lhs0
-        lhs0 = lhs
-        meshdofs0000 = meshdofs000
-        meshdofs000 = meshdofs00
-        meshdofs00 = meshdofs0
-        meshdofs0 = meshdofs
-        dt = min(precice_dt, timestepsize)
-
-        # solve mesh deformation
-        meshdofs = solver.optimize('meshdofs', meshsqr, constrain=meshcons)
-
-        # solve fluid equations
-        lhs = solver.newton('lhs', res, lhs0=lhs0, constrain=cons,
-                             arguments=dict(lhs0=lhs0, dt=dt, meshdofs=meshdofs, meshdofs0=meshdofs0,
-                                            meshdofs00=meshdofs00, meshdofs000=meshdofs000, meshdofs0000=meshdofs0000)
-                             ).solve(tol=1e-6)
+        arguments = {k+'0'*(i+1): arguments[k+'0'*i] for k in ('lhs', 'meshdofs', 'dt') for i in range(nhist)}
+        arguments['dt'] = dt = min(precice_dt, timestepsize)
+        arguments['meshdofs'] = solver.optimize('meshdofs', meshsqr, constrain=meshcons) # solve mesh deformation
+        arguments['lhs'] = arguments['lhs0'] # initial guess for newton
+        arguments['lhs'] = solver.newton('lhs', res, arguments=arguments, constrain=cons).solve(tol=1e-6) # solve fluid equations
+        arguments['F'] = solver.solve_linear('F', resF, constrain=consF, arguments=arguments)
 
         # write forces to interface
         if interface.is_write_data_required(dt):
-            F = solver.solve_linear('F', resF, constrain=consF,
-                                    arguments=dict(lhs00=lhs00, lhs0=lhs0, lhs=lhs, dt=dt, meshdofs=meshdofs,
-                                                   meshdofs0=meshdofs0, meshdofs00=meshdofs00,
-                                                   meshdofs000=meshdofs000, meshdofs0000=meshdofs0000))
-            # writedata = couplingsample.eval(ns.F, F=F) # for stresses
-            writedata = couplingsample.eval('F_i d:x' @ ns, F=F, meshdofs=meshdofs) * \
+            writedata = couplingsample.eval('F_i d:x' @ ns, **arguments) * \
                 numpy.concatenate([p.weights for p in couplingsample.points])[:, numpy.newaxis]
             interface.write_block_vector_data(writedataID, dataIndices, writedata)
 
@@ -186,12 +167,11 @@ def main(inflow: 'inflow velocity' = 10,
 
         # read checkpoint if required
         if interface.is_action_required(precice.action_read_iteration_checkpoint()):
-            timestep, lhs, lhs0, meshdofs, meshdofs0, meshdofs00, meshdofs000 = checkpoint
+            timestep, arguments = checkpoint
             interface.mark_action_fulfilled(precice.action_read_iteration_checkpoint())
 
         if interface.is_time_window_complete():
-            x, u, p = bezier.eval(['x_i', 'u_i', 'p'] @ ns, lhs=lhs, meshdofs=meshdofs, meshdofs0=meshdofs0,
-                                  meshdofs00=meshdofs00, meshdofs000=meshdofs000, meshdofs0000=meshdofs0000, dt=dt)
+            x, u, p = bezier.eval(['x_i', 'u_i', 'p'] @ ns, **arguments)
             export.triplot('velocity.jpg', x, numpy.linalg.norm(u, axis=1), tri=bezier.tri, cmap='jet')
             export.triplot('pressure.jpg', x, p, tri=bezier.tri, cmap='jet')
             with treelog.add(treelog.DataLog()):

From 60c41f1ee2897be03808cb95dc36035157fc8150 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Thu, 4 Jan 2024 12:17:38 +0100
Subject: [PATCH 17/25] introduce evaluable quadrature weights

---
 perpendicular-flap/fluid-nutils/fluid.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/perpendicular-flap/fluid-nutils/fluid.py b/perpendicular-flap/fluid-nutils/fluid.py
index be0c6a5c6..e8ec4f87e 100644
--- a/perpendicular-flap/fluid-nutils/fluid.py
+++ b/perpendicular-flap/fluid-nutils/fluid.py
@@ -36,6 +36,9 @@ def main(inflow: 'inflow velocity' = 10,
     domain = topo.withboundary(inflow='left', wall='top,bottom', outflow='right') - \
         topo[18:20, :10].withboundary(flap='left,right,top')
 
+    couplinginterface = domain.boundary['flap']
+    couplingsample = couplinginterface.sample('gauss', degree=2)  # mesh located at Gauss points
+
     # time approximations
     t0 = lambda f: function.replace_arguments(f, {arg: function.Argument(arg+'0', shape=shape, dtype=dtype)
                         for arg, (shape, dtype) in f.arguments.items()})
@@ -63,6 +66,7 @@ def main(inflow: 'inflow velocity' = 10,
     ns.urel_i = 'ubasis_ni ?lhs_n'  # relative velocity
     ns.u_i = 'umesh_i + urel_i'  # total velocity
     ns.p = 'pbasis_n ?lhs_n'  # pressure
+    ns.qw = couplingsample.asfunction(numpy.concatenate([p.weights for p in couplingsample.points]))
 
     # for visualization
     bezier = domain.sample('bezier', 2)
@@ -77,9 +81,6 @@ def main(inflow: 'inflow velocity' = 10,
     # define coupling meshes
     meshName = "Fluid-Mesh"
     meshID = interface.get_mesh_id(meshName)
-
-    couplinginterface = domain.boundary['flap']
-    couplingsample = couplinginterface.sample('gauss', degree=2)  # mesh located at Gauss points
     dataIndices = interface.set_mesh_vertices(meshID, couplingsample.eval(ns.x0))
 
     # coupling data
@@ -158,8 +159,7 @@ def main(inflow: 'inflow velocity' = 10,
 
         # write forces to interface
         if interface.is_write_data_required(dt):
-            writedata = couplingsample.eval('F_i d:x' @ ns, **arguments) * \
-                numpy.concatenate([p.weights for p in couplingsample.points])[:, numpy.newaxis]
+            writedata = couplingsample.eval('F_i qw d:x' @ ns, **arguments)
             interface.write_block_vector_data(writedataID, dataIndices, writedata)
 
         # do the coupling

From 85f68e29fb52c5d2b35f3ae1828869e4d5c9d04c Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Wed, 3 Jan 2024 22:06:05 +0100
Subject: [PATCH 18/25] simplify mesh generation code

---
 perpendicular-flap/fluid-nutils/fluid.py | 32 ++++++++++--------------
 1 file changed, 13 insertions(+), 19 deletions(-)

diff --git a/perpendicular-flap/fluid-nutils/fluid.py b/perpendicular-flap/fluid-nutils/fluid.py
index e8ec4f87e..d1edf8838 100644
--- a/perpendicular-flap/fluid-nutils/fluid.py
+++ b/perpendicular-flap/fluid-nutils/fluid.py
@@ -14,25 +14,19 @@ def main(inflow: 'inflow velocity' = 10,
          timestepsize=0.01):
 
     # mesh and geometry definition
-    grid_x_1 = numpy.linspace(-3, -1, 7)
-    grid_x_1 = grid_x_1[:-1]
-    grid_x_2 = numpy.linspace(-1, -0.3, 8)
-    grid_x_2 = grid_x_2[:-1]
-    grid_x_3 = numpy.linspace(-0.3, 0.3, 13)
-    grid_x_3 = grid_x_3[:-1]
-    grid_x_4 = numpy.linspace(0.3, 1, 8)
-    grid_x_4 = grid_x_4[:-1]
-    grid_x_5 = numpy.linspace(1, 3, 7)
-    grid_x = numpy.concatenate((grid_x_1, grid_x_2, grid_x_3, grid_x_4, grid_x_5), axis=None)
-    grid_y_1 = numpy.linspace(0, 1.5, 16)
-    grid_y_1 = grid_y_1[:-1]
-    grid_y_2 = numpy.linspace(1.5, 2, 4)
-    grid_y_2 = grid_y_2[:-1]
-    grid_y_3 = numpy.linspace(2, 4, 7)
-    grid_y = numpy.concatenate((grid_y_1, grid_y_2, grid_y_3), axis=None)
-    grid = [grid_x, grid_y]
-
-    topo, geom = mesh.rectilinear(grid)
+    topo, geom = mesh.rectilinear([
+        numpy.concatenate((
+            numpy.linspace(-3, -1, 6, endpoint=False),
+            numpy.linspace(-1, -0.3, 7, endpoint=False),
+            numpy.linspace(-0.3, 0.3, 12, endpoint=False),
+            numpy.linspace(0.3, 1, 7, endpoint=False),
+            numpy.linspace(1, 3, 7))),
+        numpy.concatenate((
+            numpy.linspace(0, 1.5, 15, endpoint=False),
+            numpy.linspace(1.5, 2, 3, endpoint=False),
+            numpy.linspace(2, 4, 7))),
+    ])
+
     domain = topo.withboundary(inflow='left', wall='top,bottom', outflow='right') - \
         topo[18:20, :10].withboundary(flap='left,right,top')
 

From 47931253ea7de8a91287e641b158483efceb9364 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Thu, 4 Jan 2024 13:59:55 +0100
Subject: [PATCH 19/25] simplify computation of consF

---
 perpendicular-flap/fluid-nutils/fluid.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/perpendicular-flap/fluid-nutils/fluid.py b/perpendicular-flap/fluid-nutils/fluid.py
index d1edf8838..7f45f0782 100644
--- a/perpendicular-flap/fluid-nutils/fluid.py
+++ b/perpendicular-flap/fluid-nutils/fluid.py
@@ -105,8 +105,7 @@ def main(inflow: 'inflow velocity' = 10,
     resF += tt(domain.integral('rho ubasis_ni u_i d:x' @ ns, degree=4))
     resF += domain.integral('rho ubasis_ni (u_i,j urel_j d:x)' @ ns, degree=4)
     resF += couplinginterface.sample('gauss', 4).integral('ubasis_ni F_i d:x' @ ns)
-    consF = numpy.isnan(solver.optimize('F', couplinginterface.sample('gauss', 4).integral('F_i F_i' @ ns),
-                                        droptol=1e-10))
+    consF = couplingsample.integrate((ns.ubasis**2).sum(1)) == 0
 
     # boundary conditions mesh displacements
     sqr = domain.boundary['inflow,outflow,wall'].integral('d_i d_i' @ ns, degree=2)

From d89f08de059de1e9b5acc6e63a89c872eb8b4bfd Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Mon, 8 Jan 2024 14:01:35 +0100
Subject: [PATCH 20/25] simplify weak traction evaluation

---
 perpendicular-flap/fluid-nutils/fluid.py | 7 +------
 1 file changed, 1 insertion(+), 6 deletions(-)

diff --git a/perpendicular-flap/fluid-nutils/fluid.py b/perpendicular-flap/fluid-nutils/fluid.py
index 7f45f0782..28fedf721 100644
--- a/perpendicular-flap/fluid-nutils/fluid.py
+++ b/perpendicular-flap/fluid-nutils/fluid.py
@@ -99,12 +99,7 @@ def main(inflow: 'inflow velocity' = 10,
     res += tθ(domain.integral('rho ubasis_ni u_i,j urel_j d:x' @ ns, degree=4))
 
     # weak form for force computation
-    resF = domain.integral('(ubasis_ni,j (u_i,j + u_j,i) rho nu d:x)' @ ns, degree=4)
-    resF += tp(domain.integral('-ubasis_ni,j p δ_ij d:x' @ ns, degree=4))
-    resF += domain.integral('pbasis_n u_k,k d:x' @ ns, degree=4)
-    resF += tt(domain.integral('rho ubasis_ni u_i d:x' @ ns, degree=4))
-    resF += domain.integral('rho ubasis_ni (u_i,j urel_j d:x)' @ ns, degree=4)
-    resF += couplinginterface.sample('gauss', 4).integral('ubasis_ni F_i d:x' @ ns)
+    resF = res + couplinginterface.sample('gauss', 4).integral('ubasis_ni F_i d:x' @ ns)
     consF = couplingsample.integrate((ns.ubasis**2).sum(1)) == 0
 
     # boundary conditions mesh displacements

From 8459059a506eb1dd9b8ca01308af288e9bdc875a Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Thu, 4 Jan 2024 19:31:50 +0100
Subject: [PATCH 21/25] change umesh to first order accuracy

---
 perpendicular-flap/fluid-nutils/fluid.py | 8 ++------
 1 file changed, 2 insertions(+), 6 deletions(-)

diff --git a/perpendicular-flap/fluid-nutils/fluid.py b/perpendicular-flap/fluid-nutils/fluid.py
index 28fedf721..0c1b40e5e 100644
--- a/perpendicular-flap/fluid-nutils/fluid.py
+++ b/perpendicular-flap/fluid-nutils/fluid.py
@@ -40,10 +40,6 @@ def main(inflow: 'inflow velocity' = 10,
     tθ = lambda f: theta * f + (1 - theta) * t0(f)
     # 1st order FD
     δt = lambda f: (f - t0(f)) / function.Argument('dt', ())
-    # extrapolation for pressure
-    tp = lambda f: (1.5 * f - 0.5 * t0(f))
-    # 2nd order FD
-    tt = lambda f: tp(δt(f))
 
     # Nutils namespace
     ns = function.Namespace()
@@ -53,7 +49,7 @@ def main(inflow: 'inflow velocity' = 10,
     ns.x0 = geom  # reference geometry
     ns.dbasis = domain.basis('std', degree=1).vector(2)
     ns.d_i = 'dbasis_ni ?meshdofs_n'
-    ns.umesh = tt(ns.d) # mesh velocity
+    ns.umesh = δt(ns.d) # mesh velocity
     ns.x_i = 'x0_i + d_i'  # moving geometry
     ns.ubasis, ns.pbasis = function.chain([domain.basis('std', degree=2).vector(2), domain.basis('std', degree=1), ])
     ns.F_i = 'ubasis_ni ?F_n'  # stress field
@@ -118,7 +114,7 @@ def main(inflow: 'inflow velocity' = 10,
     timestep = 0
     arguments = dict(lhs=numpy.zeros(len(ns.ubasis)), meshdofs=numpy.zeros(len(ns.dbasis)), dt=timestepsize)
 
-    nhist = 4 # history length required by the formulation
+    nhist = 2 # history length required by the formulation
     arguments.update((k+'0'*i, v) for k, v in tuple(arguments.items()) for i in range(nhist))
 
     while interface.is_coupling_ongoing():

From 6c2a551e68c9e7c627785dd837e442d54a41d52d Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Thu, 4 Jan 2024 21:02:47 +0100
Subject: [PATCH 22/25] change main signature for nutils 8 compatibility

---
 perpendicular-flap/fluid-nutils/fluid.py | 6 +-----
 1 file changed, 1 insertion(+), 5 deletions(-)

diff --git a/perpendicular-flap/fluid-nutils/fluid.py b/perpendicular-flap/fluid-nutils/fluid.py
index 0c1b40e5e..55ba50919 100644
--- a/perpendicular-flap/fluid-nutils/fluid.py
+++ b/perpendicular-flap/fluid-nutils/fluid.py
@@ -7,11 +7,7 @@
 
 # for details on this solver see https://doi.org/10.1002/nme.6443
 
-def main(inflow: 'inflow velocity' = 10,
-         viscosity: 'kinematic viscosity' = 1.0,
-         density: 'density' = 1.0,
-         theta=0.5,
-         timestepsize=0.01):
+def main(inflow=10., viscosity=1., density=1., theta=.5, timestepsize=.01):
 
     # mesh and geometry definition
     topo, geom = mesh.rectilinear([

From d406c27864a8caafbe578c81e0789905d6ed992b Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Fri, 5 Jan 2024 15:57:17 +0100
Subject: [PATCH 23/25] separate sample definitions from precice setup

---
 perpendicular-flap/fluid-nutils/fluid.py | 40 +++++++++---------------
 1 file changed, 15 insertions(+), 25 deletions(-)

diff --git a/perpendicular-flap/fluid-nutils/fluid.py b/perpendicular-flap/fluid-nutils/fluid.py
index 55ba50919..91dfa7330 100644
--- a/perpendicular-flap/fluid-nutils/fluid.py
+++ b/perpendicular-flap/fluid-nutils/fluid.py
@@ -27,7 +27,9 @@ def main(inflow=10., viscosity=1., density=1., theta=.5, timestepsize=.01):
         topo[18:20, :10].withboundary(flap='left,right,top')
 
     couplinginterface = domain.boundary['flap']
-    couplingsample = couplinginterface.sample('gauss', degree=2)  # mesh located at Gauss points
+    couplingsample = couplinginterface.sample('uniform', degree=npoints_per_elem)
+
+    bezier = domain.sample('bezier', 2)
 
     # time approximations
     t0 = lambda f: function.replace_arguments(f, {arg: function.Argument(arg+'0', shape=shape, dtype=dtype)
@@ -54,30 +56,6 @@ def main(inflow=10., viscosity=1., density=1., theta=.5, timestepsize=.01):
     ns.p = 'pbasis_n ?lhs_n'  # pressure
     ns.qw = couplingsample.asfunction(numpy.concatenate([p.weights for p in couplingsample.points]))
 
-    # for visualization
-    bezier = domain.sample('bezier', 2)
-
-    # preCICE setup
-    configFileName = "../precice-config.xml"
-    participantName = "Fluid"
-    solverProcessIndex = 0
-    solverProcessSize = 1
-    interface = precice.Interface(participantName, configFileName, solverProcessIndex, solverProcessSize)
-
-    # define coupling meshes
-    meshName = "Fluid-Mesh"
-    meshID = interface.get_mesh_id(meshName)
-    dataIndices = interface.set_mesh_vertices(meshID, couplingsample.eval(ns.x0))
-
-    # coupling data
-    writeData = "Force"
-    readData = "Displacement"
-    writedataID = interface.get_data_id(writeData, meshID)
-    readdataID = interface.get_data_id(readData, meshID)
-
-    # initialize preCICE
-    precice_dt = interface.initialize()
-
     # boundary conditions for fluid equations
     sqr = domain.boundary['wall,flap'].integral('urel_k urel_k d:x0' @ ns, degree=4)
     cons = solver.optimize('lhs', sqr, droptol=1e-15)
@@ -107,6 +85,18 @@ def main(inflow=10., viscosity=1., density=1., theta=.5, timestepsize=.01):
     #res_stokes = domain.integral('(ubasis_ni,j ((u_i,j + u_j,i) rho nu - p δ_ij) + pbasis_n u_k,k) d:x' @ ns, degree=4)
     #lhs0 = solver.solve_linear('lhs', res_stokes, constrain=cons, arguments=dict(meshdofs=meshdofs, meshdofs0=meshdofs0, meshdofs00=meshdofs00, meshdofs000=meshdofs000, dt=dt))
 
+    # preCICE setup
+    solverProcessIndex = 0
+    solverProcessSize = 1
+    interface = precice.Interface("Fluid", "../precice-config.xml", solverProcessIndex, solverProcessSize)
+    meshID = interface.get_mesh_id("Fluid-Mesh")
+    dataIndices = interface.set_mesh_vertices(meshID, couplingsample.eval(ns.x0))
+    writedataID = interface.get_data_id("Force", meshID)
+    readdataID = interface.get_data_id("Displacement", meshID)
+
+    # initialize preCICE
+    precice_dt = interface.initialize()
+
     timestep = 0
     arguments = dict(lhs=numpy.zeros(len(ns.ubasis)), meshdofs=numpy.zeros(len(ns.dbasis)), dt=timestepsize)
 

From 35b8b3425367df39bb94641faa597845b9bc5a8c Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Fri, 5 Jan 2024 15:44:46 +0100
Subject: [PATCH 24/25] switch to uniform sample

---
 perpendicular-flap/fluid-nutils/fluid.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/perpendicular-flap/fluid-nutils/fluid.py b/perpendicular-flap/fluid-nutils/fluid.py
index 91dfa7330..36ed35371 100644
--- a/perpendicular-flap/fluid-nutils/fluid.py
+++ b/perpendicular-flap/fluid-nutils/fluid.py
@@ -7,7 +7,7 @@
 
 # for details on this solver see https://doi.org/10.1002/nme.6443
 
-def main(inflow=10., viscosity=1., density=1., theta=.5, timestepsize=.01):
+def main(inflow=10., viscosity=1., density=1., theta=.5, timestepsize=.01, npoints_per_elem=3):
 
     # mesh and geometry definition
     topo, geom = mesh.rectilinear([
@@ -54,7 +54,7 @@ def main(inflow=10., viscosity=1., density=1., theta=.5, timestepsize=.01):
     ns.urel_i = 'ubasis_ni ?lhs_n'  # relative velocity
     ns.u_i = 'umesh_i + urel_i'  # total velocity
     ns.p = 'pbasis_n ?lhs_n'  # pressure
-    ns.qw = couplingsample.asfunction(numpy.concatenate([p.weights for p in couplingsample.points]))
+    ns.qw = 1 / npoints_per_elem
 
     # boundary conditions for fluid equations
     sqr = domain.boundary['wall,flap'].integral('urel_k urel_k d:x0' @ ns, degree=4)

From 3e17c0dafb263fec3bf610f9fe82a6b9b930e712 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Mon, 8 Jan 2024 15:11:34 +0100
Subject: [PATCH 25/25] change to first order formulation

---
 perpendicular-flap/fluid-nutils/fluid.py | 12 +++++-------
 1 file changed, 5 insertions(+), 7 deletions(-)

diff --git a/perpendicular-flap/fluid-nutils/fluid.py b/perpendicular-flap/fluid-nutils/fluid.py
index 36ed35371..5a6e83a1d 100644
--- a/perpendicular-flap/fluid-nutils/fluid.py
+++ b/perpendicular-flap/fluid-nutils/fluid.py
@@ -33,7 +33,7 @@ def main(inflow=10., viscosity=1., density=1., theta=.5, timestepsize=.01, npoin
 
     # time approximations
     t0 = lambda f: function.replace_arguments(f, {arg: function.Argument(arg+'0', shape=shape, dtype=dtype)
-                        for arg, (shape, dtype) in f.arguments.items()})
+                        for arg, (shape, dtype) in f.arguments.items() if arg in ('lhs', 'meshdofs', 'F')})
     # TR interpolation
     tθ = lambda f: theta * f + (1 - theta) * t0(f)
     # 1st order FD
@@ -47,7 +47,7 @@ def main(inflow=10., viscosity=1., density=1., theta=.5, timestepsize=.01, npoin
     ns.x0 = geom  # reference geometry
     ns.dbasis = domain.basis('std', degree=1).vector(2)
     ns.d_i = 'dbasis_ni ?meshdofs_n'
-    ns.umesh = δt(ns.d) # mesh velocity
+    ns.umesh_i = 'dbasis_ni ?umesh_n' # mesh velocity
     ns.x_i = 'x0_i + d_i'  # moving geometry
     ns.ubasis, ns.pbasis = function.chain([domain.basis('std', degree=2).vector(2), domain.basis('std', degree=1), ])
     ns.F_i = 'ubasis_ni ?F_n'  # stress field
@@ -98,10 +98,7 @@ def main(inflow=10., viscosity=1., density=1., theta=.5, timestepsize=.01, npoin
     precice_dt = interface.initialize()
 
     timestep = 0
-    arguments = dict(lhs=numpy.zeros(len(ns.ubasis)), meshdofs=numpy.zeros(len(ns.dbasis)), dt=timestepsize)
-
-    nhist = 2 # history length required by the formulation
-    arguments.update((k+'0'*i, v) for k, v in tuple(arguments.items()) for i in range(nhist))
+    arguments = dict(lhs=numpy.zeros(len(ns.ubasis)), meshdofs=numpy.zeros(len(ns.dbasis)))
 
     while interface.is_coupling_ongoing():
       with treelog.context(f'timestep {timestep}'):
@@ -120,9 +117,10 @@ def main(inflow=10., viscosity=1., density=1., theta=.5, timestepsize=.01, npoin
 
         # advance variables
         timestep += 1
-        arguments = {k+'0'*(i+1): arguments[k+'0'*i] for k in ('lhs', 'meshdofs', 'dt') for i in range(nhist)}
+        arguments = {k+'0': arguments[k] for k in ('lhs', 'meshdofs')}
         arguments['dt'] = dt = min(precice_dt, timestepsize)
         arguments['meshdofs'] = solver.optimize('meshdofs', meshsqr, constrain=meshcons) # solve mesh deformation
+        arguments['umesh'] = (arguments['meshdofs'] - arguments['meshdofs0']) / dt
         arguments['lhs'] = arguments['lhs0'] # initial guess for newton
         arguments['lhs'] = solver.newton('lhs', res, arguments=arguments, constrain=cons).solve(tol=1e-6) # solve fluid equations
         arguments['F'] = solver.solve_linear('F', resF, constrain=consF, arguments=arguments)