From 378b089b1393ac1ce10943874e2f095d660d8f94 Mon Sep 17 00:00:00 2001
From: Benjamin Rodenberg <benjamin.rodenberg@in.tum.de>
Date: Mon, 20 Nov 2023 14:44:49 +0100
Subject: [PATCH] Use waveform iteration with
 partitioned-heat-conduction/fenics/heat.py (#281)

---
 partitioned-heat-conduction/fenics/heat.py    | 77 ++++++++++++++++---
 .../precice-config.xml                        | 38 ++++++---
 2 files changed, 93 insertions(+), 22 deletions(-)

diff --git a/partitioned-heat-conduction/fenics/heat.py b/partitioned-heat-conduction/fenics/heat.py
index 661d078a8..02b27eee8 100644
--- a/partitioned-heat-conduction/fenics/heat.py
+++ b/partitioned-heat-conduction/fenics/heat.py
@@ -63,7 +63,7 @@ def determine_gradient(V_g, u, flux):
 
 args = parser.parse_args()
 
-fenics_dt = .1  # time step size
+fenics_dt = .01  # time step size
 # Error is bounded by coupling accuracy. In theory we would obtain the analytical solution.
 error_tol = args.error_tol
 
@@ -112,7 +112,6 @@ def determine_gradient(V_g, u, flux):
     precice.initialize(coupling_boundary, read_function_space=W, write_object=u_D_function)
     precice_dt = precice.get_max_time_step_size()
 
-
 dt = Constant(0)
 dt.assign(np.min([fenics_dt, precice_dt]))
 
@@ -164,13 +163,32 @@ def determine_gradient(V_g, u, flux):
 
 # output solution and reference solution at t=0, n=0
 n = 0
-print('output u^%d and u_ref^%d' % (n, n))
-temperature_out << (u_n, t)
-ref_out << u_ref
+print("output u^%d and u_ref^%d" % (n, n))
 ranks << mesh_rank
 
 error_total, error_pointwise = compute_errors(u_n, u_ref, V)
-error_out << error_pointwise
+
+# create buffer for output. We need this buffer, because we only want to
+# write the converged output at the end of the window, but we also want to
+# write the samples that are resulting from substeps inside the window
+u_write = []
+ref_write = []
+error_write = []
+# copy data to buffer and rename
+uu = u_n.copy()
+uu.rename("u", "")
+u_write.append((uu, t))
+uu_ref = u_ref.copy()
+uu_ref.rename("u_ref", "")
+ref_write.append(uu_ref)
+err = error_pointwise.copy()
+err.rename("err", "")
+error_write.append(err)
+
+# set t_1 = t_0 + dt, this gives u_D^1
+# call dt(0) to evaluate FEniCS Constant. Todo: is there a better way?
+u_D.t = t + dt(0)
+f.t = t + dt(0)
 
 if problem is ProblemType.DIRICHLET:
     flux = Function(V_g)
@@ -182,6 +200,17 @@ def determine_gradient(V_g, u, flux):
     if precice.requires_writing_checkpoint():
         precice.store_checkpoint(u_n, t, n)
 
+        # output solution and reference solution at t_n+1 and substeps (read from buffer)
+        print('output u^%d and u_ref^%d' % (n, n))
+        for sample in u_write:
+            temperature_out << sample
+
+        for sample in ref_write:
+            ref_out << sample
+
+        for sample in error_write:
+            error_out << error_pointwise
+
     precice_dt = precice.get_max_time_step_size()
     dt.assign(np.min([fenics_dt, precice_dt]))
 
@@ -217,21 +246,45 @@ def determine_gradient(V_g, u, flux):
         u_n.assign(u_cp)
         t = t_cp
         n = n_cp
+        # empty buffer if window has not converged
+        u_write = []
+        ref_write = []
+        error_write = []
     else:  # update solution
         u_n.assign(u_np1)
         t += float(dt)
         n += 1
+        # copy data to buffer and rename
+        uu = u_n.copy()
+        uu.rename("u", "")
+        u_write.append((uu, t))
+        uu_ref = u_ref.copy()
+        uu_ref.rename("u_ref", "")
+        ref_write.append(uu_ref)
+        err = error_pointwise.copy()
+        err.rename("err", "")
+        error_write.append(err)
 
     if precice.is_time_window_complete():
         u_ref = interpolate(u_D, V)
         u_ref.rename("reference", " ")
         error, error_pointwise = compute_errors(u_n, u_ref, V, total_error_tol=error_tol)
-        print('n = %d, t = %.2f: L2 error on domain = %.3g' % (n, t, error))
-        # output solution and reference solution at t_n+1
-        print('output u^%d and u_ref^%d' % (n, n))
-        temperature_out << (u_n, t)
-        ref_out << u_ref
-        error_out << error_pointwise
+        print("n = %d, t = %.2f: L2 error on domain = %.3g" % (n, t, error))
+
+    # Update Dirichlet BC
+    u_D.t = t + float(dt)
+    f.t = t + float(dt)
+
+# output solution and reference solution at t_n+1 and substeps (read from buffer)
+print("output u^%d and u_ref^%d" % (n, n))
+for sample in u_write:
+    temperature_out << sample
+
+for sample in ref_write:
+    ref_out << sample
+
+for sample in error_write:
+    error_out << error_pointwise
 
 # Hold plot
 precice.finalize()
diff --git a/partitioned-heat-conduction/precice-config.xml b/partitioned-heat-conduction/precice-config.xml
index 3b20fad0d..baab5d5c3 100644
--- a/partitioned-heat-conduction/precice-config.xml
+++ b/partitioned-heat-conduction/precice-config.xml
@@ -7,8 +7,8 @@
       enabled="true" />
   </log>
 
-  <data:scalar name="Temperature" />
-  <data:scalar name="Heat-Flux" />
+  <data:scalar name="Temperature" waveform-degree="1" />
+  <data:scalar name="Heat-Flux" waveform-degree="1" />
 
   <mesh name="Dirichlet-Mesh" dimensions="2">
     <use-data name="Temperature" />
@@ -25,7 +25,11 @@
     <receive-mesh name="Neumann-Mesh" from="Neumann" />
     <write-data name="Heat-Flux" mesh="Dirichlet-Mesh" />
     <read-data name="Temperature" mesh="Dirichlet-Mesh" />
-    <mapping:nearest-neighbor direction="read" from="Neumann-Mesh" to="Dirichlet-Mesh" constraint="consistent" />
+    <mapping:nearest-neighbor
+      direction="read"
+      from="Neumann-Mesh"
+      to="Dirichlet-Mesh"
+      constraint="consistent" />
   </participant>
 
   <participant name="Neumann">
@@ -33,7 +37,11 @@
     <receive-mesh name="Dirichlet-Mesh" from="Dirichlet" />
     <write-data name="Temperature" mesh="Neumann-Mesh" />
     <read-data name="Heat-Flux" mesh="Neumann-Mesh" />
-    <mapping:nearest-neighbor direction="read" from="Dirichlet-Mesh" to="Neumann-Mesh" constraint="consistent" />
+    <mapping:nearest-neighbor
+      direction="read"
+      from="Dirichlet-Mesh"
+      to="Neumann-Mesh"
+      constraint="consistent" />
   </participant>
 
   <m2n:sockets acceptor="Dirichlet" connector="Neumann" exchange-directory=".." />
@@ -43,21 +51,31 @@
     <max-time value="1.0" />
     <time-window-size value="0.1" />
     <max-iterations value="100" />
-    <exchange data="Heat-Flux" mesh="Dirichlet-Mesh" from="Dirichlet" to="Neumann" />
+    <exchange
+      data="Heat-Flux"
+      mesh="Dirichlet-Mesh"
+      from="Dirichlet"
+      to="Neumann"
+      initialize="true"
+      substeps="true" />
     <exchange
       data="Temperature"
       mesh="Neumann-Mesh"
       from="Neumann"
       to="Dirichlet"
-      initialize="true" />
-    <relative-convergence-measure data="Heat-Flux" mesh="Dirichlet-Mesh" limit="1e-5" />
-    <relative-convergence-measure data="Temperature" mesh="Neumann-Mesh" limit="1e-5" />
-    <acceleration:IQN-ILS>
+      initialize="true"
+      substeps="true" />
+    <relative-convergence-measure data="Heat-Flux" mesh="Dirichlet-Mesh" limit="1e-10" />
+    <relative-convergence-measure data="Temperature" mesh="Neumann-Mesh" limit="1e-10" />
+    <!-- <acceleration:IQN-ILS>
       <data name="Temperature" mesh="Neumann-Mesh" />
       <initial-relaxation value="0.1" />
       <max-used-iterations value="10" />
       <time-windows-reused value="5" />
       <filter type="QR2" limit="1e-3" />
-    </acceleration:IQN-ILS>
+    </acceleration:IQN-ILS> -->
+    <acceleration:constant>
+      <relaxation value="0.5" />
+    </acceleration:constant>
   </coupling-scheme:serial-implicit>
 </precice-configuration>