benchmark

benchmark.py

@@ -151,7 +151,7 @@ def _elevate_degree(mesh, degree):
     mesh_s = Submesh(mesh, dmcommon.CELL_SETS_LABEL, label_struct, mesh.topological_dimension())
     mesh_s = _elevate_degree(mesh_s, degree)
 else:
-    nref = 0
+    nref = 1
     mesh = make_mesh(quadrilateral)
     if nref > 0:
         mesh = MeshHierarchy(mesh, nref)[-1]
@@ -626,14 +626,14 @@ def v_f_left(t_):
     print(f"Time: {end - start}")
 elif case in ["FSI1_2", "FSI2_2", "FSI3_2"]:
     T = 20 # 10.0 # 12.0
-    dt = Constant(0.002)  #0.001
+    dt = Constant(0.001)  #0.001
     chkstride = int(0.01 / float(dt))  # save every 0.01 sec
     t = Constant(0.0)
     CNshift = 1
     elast = True
     linear_elast = True
     serendipity = True
-    slip_bc = False
+    slip_bc = True
     if use_netgen:
         fname_checkpoint = f"dumbdata/fsi3_P4_P2_nref{nref}_0.002_shift{CNshift}_{elast}_{linear_elast}_netgen_new"
         fname_FD = f"time_series_FD_P4_P2_nref{nref}_0.002_shift{CNshift}_{elast}_{linear_elast}_netgen_new.dat"
@@ -643,11 +643,11 @@ def v_f_left(t_):
             #fname_checkpoint = f"dumbdata/fsi3_Q4_Q3_nref{nref}_0.002_shift{CNshift}_{elast}_{linear_elast}_ALEparamvar1_S"
             #fname_FD = f"time_series_FD_Q4_Q3_nref{nref}_0.002_shift{CNshift}_{elast}_{linear_elast}_ALEparamvar1_S.dat"
             #fname_FL = f"time_series_FL_Q4_Q3_nref{nref}_0.002_shift{CNshift}_{elast}_{linear_elast}_ALEparamvar1_S.dat"
-            fname_checkpoint = f"dumbdata/fsi3_Q4_Q3_nref{nref}_0.002_shift{CNshift}_{elast}_{linear_elast}_{slip_bc}_ALEparamvar1_S"
-            fname_FD = f"time_series_FD_Q4_Q3_nref{nref}_0.002_shift{CNshift}_{elast}_{linear_elast}_{slip_bc}_ALEparamvar1_S.dat"
-            fname_FL = f"time_series_FL_Q4_Q3_nref{nref}_0.002_shift{CNshift}_{elast}_{linear_elast}_{slip_bc}_ALEparamvar1_S.dat"
-            fname_ux = f"time_series_ux_Q4_Q3_nref{nref}_0.002_shift{CNshift}_{elast}_{linear_elast}_{slip_bc}_ALEparamvar1_S.dat"
-            fname_uy = f"time_series_uy_Q4_Q3_nref{nref}_0.002_shift{CNshift}_{elast}_{linear_elast}_{slip_bc}_ALEparamvar1_S.dat"
+            fname_checkpoint = f"dumbdata/fsi3_Q4_Q3_nref{nref}_0.001_shift{CNshift}_{elast}_{linear_elast}_{slip_bc}_ALEparam1_S"
+            fname_FD = f"time_series_FD_Q4_Q3_nref{nref}_0.001_shift{CNshift}_{elast}_{linear_elast}_{slip_bc}_ALEparam1_S.dat"
+            fname_FL = f"time_series_FL_Q4_Q3_nref{nref}_0.001_shift{CNshift}_{elast}_{linear_elast}_{slip_bc}_ALEparam1_S.dat"
+            fname_ux = f"time_series_ux_Q4_Q3_nref{nref}_0.001_shift{CNshift}_{elast}_{linear_elast}_{slip_bc}_ALEparam1_S.dat"
+            fname_uy = f"time_series_uy_Q4_Q3_nref{nref}_0.001_shift{CNshift}_{elast}_{linear_elast}_{slip_bc}_ALEparam1_S.dat"
         else:
             fname_checkpoint = f"dumbdata/fsi3_P4_P2_nref{nref}_0.002_shift{CNshift}_{elast}_{linear_elast}_ALEparamtest1"
             fname_FD = f"time_series_FD_P4_P2_nref{nref}_0.002_shift{CNshift}_{elast}_{linear_elast}_ALEparamtest1.dat"
@@ -682,11 +682,11 @@ def v_f_left(t_):
     E_s = mu_s * 2 * (1 + nu_s)
     lambda_s = nu_s * E_s / (1 + nu_s) / (1 - 2 * nu_s)
     # ALE constants
-    #nu_ale = Constant(0.49)
+    nu_ale = Constant(0.49)
     mu_ale = Constant(float(mu_s))
-    decay_r = Constant(0.5)  # decay ratio
-    profile = (1. - decay_r * abs(x_f - pointA[0]) / (L - pointA[0])) * (1. - decay_r * abs(y_f - pointA[1]) / (L - pointA[0]))
-    nu_ale = Constant(0.49) * profile
+    #decay_r = Constant(0.5)  # decay ratio
+    #profile = (1. - decay_r * abs(x_f - pointA[0]) / (L - pointA[0])) * (1. - decay_r * abs(y_f - pointA[1]) / (L - pointA[0]))
+    #nu_ale = Constant(0.49) * profile
     #mu_ale = Constant(2 * float(mu_s)) * profile
     E_ale = mu_ale * 2 * (1 + nu_ale)
     lambda_ale = nu_ale * E_ale / (1 + nu_ale) / (1 - 2 * nu_ale)
@@ -933,7 +933,7 @@ def v_f_left(t_):
                  outfile.write("t val" + "\n")
             with open(fname_uy, 'w') as outfile:
                  outfile.write("t val" + "\n")
-    if True:
+    if False:
         Vplot = FunctionSpace(mesh_f, "P", degree)
         vfplot = Function(Vplot, name="mesh_f").assign(Constant(1))
         #vf = solution.subfunctions[0]

firedrake/variational_solver.py

@@ -366,7 +366,7 @@ def __init__(self, a, L, u, bcs=None, aP=None,
         else:
             if not isinstance(L, (ufl.BaseForm, slate.slate.TensorBase)):
                 raise TypeError("Provided RHS is a '%s', not a Form or Slate Tensor" % type(L).__name__)
-            if len(L.arguments()) != 1:
+            if len(L.arguments()) != 1 and not L.empty():
                 raise ValueError("Provided RHS is not a linear form")
             F = ufl_expr.action(J, u) - L