Fix checkpointing and subcycling for solid.py

perpendicular-flap/solid-fenics/solid.py

@@ -59,7 +59,6 @@ def neumann_boundary(x, on_boundary):
 v = TestFunction(V)
 
 u_np1 = Function(V)
-saved_u_old = Function(V)
 
 # function known from previous timestep
 u_n = Function(V)
@@ -80,15 +79,18 @@ def neumann_boundary(x, on_boundary):
 
 precice_dt = precice.get_max_time_step_size()
 fenics_dt = precice_dt  # if fenics_dt == precice_dt, no subcycling is applied
-# fenics_dt = 0.02  # if fenics_dt < precice_dt, subcycling is applied
+# n_substeps = 5  # number of substeps per window
+# fenics_dt = precice_dt / n_substeps  # if fenics_dt < precice_dt, subcycling is applied
 dt = Constant(np.min([precice_dt, fenics_dt]))
 
 # clamp the beam at the bottom
 bc = DirichletBC(V, Constant((0, 0)), fixed_boundary)
 
 # alpha method parameters
-alpha_m = Constant(0)
-alpha_f = Constant(0)
+alpha_m = Constant(0.2)
+alpha_f = Constant(0.4)
+# alpha_m = Constant(0)
+# alpha_f = Constant(0)
 
 """
 Check requirements for alpha_m and alpha_f from
@@ -196,13 +198,14 @@ def avg(x_old, x_new, alpha):
 while precice.is_coupling_ongoing():
 
     if precice.requires_writing_checkpoint():  # write checkpoint
-        precice.store_checkpoint(u_n, t, n)
+        precice.store_checkpoint((u_n, v_n, a_n), t, n)
 
     precice_dt = precice.get_max_time_step_size()
     dt = Constant(np.min([precice_dt, fenics_dt]))
 
     # read data from preCICE and get a new coupling expression
-    read_data = precice.read_data(dt)
+    # sample force F at $F(t_{n+1-\alpha_f})$ (see generalized alpha paper)
+    read_data = precice.read_data((1 - float(alpha_f)) * dt)
 
     # Update the point sources on the coupling boundary with the new read data
     Forces_x, Forces_y = precice.get_point_sources(read_data)
@@ -227,17 +230,20 @@ def avg(x_old, x_new, alpha):
 
     # Either revert to old step if timestep has not converged or move to next timestep
     if precice.requires_reading_checkpoint():  # roll back to checkpoint
-        u_cp, t_cp, n_cp = precice.retrieve_checkpoint()
+        uva_cp, t_cp, n_cp = precice.retrieve_checkpoint()
+        u_cp, v_cp, a_cp = uva_cp
         u_n.assign(u_cp)
+        v_n.assign(v_cp)
+        a_n.assign(a_cp)
         t = t_cp
         n = n_cp
     else:
+        update_fields(u_np1, u_n, v_n, a_n)
         u_n.assign(u_np1)
         t += float(dt)
         n += 1
 
     if precice.is_time_window_complete():
-        update_fields(u_np1, saved_u_old, v_n, a_n)
         if n % 10 == 0:
             displacement_out << (u_n, t)