benchmark

Author: ksagiyam

benchmark.py

@@ -173,7 +173,7 @@ def _elevate_degree(mesh, degree):
     return Mesh(f)
 
 
-use_netgen = False
+use_netgen = True
 quadrilateral = True
 
 T = 20 # 10.0 # 12.0
@@ -183,17 +183,17 @@ def _elevate_degree(mesh, degree):
 ntimesteps = int(T / dt_float)
 t = Constant(0.0)
 dim = 2
-degree = 3  # 2 - 4
+degree = 2  # 2 - 4
 if use_netgen:
-    nref = 2 #  # 2 - 5 tested for CSM 1 and 2
+    nref = 1 #  # 2 - 5 tested for CSM 1 and 2
     mesh  = make_mesh_netgen(0.1 / 2 ** nref)
     mesh = _elevate_degree(mesh, degree)
     mesh_f = Submesh(mesh, dmcommon.CELL_SETS_LABEL, label_fluid, mesh.topological_dimension())
     mesh_f = _elevate_degree(mesh_f, degree)
     mesh_s = Submesh(mesh, dmcommon.CELL_SETS_LABEL, label_struct, mesh.topological_dimension())
     mesh_s = _elevate_degree(mesh_s, degree)
 else:
-    nref = 1
+    nref = 0
     mesh = make_mesh(quadrilateral)
     if nref > 0:
         mesh = MeshHierarchy(mesh, nref)[-1]
@@ -658,6 +658,8 @@ def compute_elast_tensors(dim, u, lambda_s, mu_s):
         S = lambda_s * tr(E) * Identity(dim) + 2.0 * mu_s * E
         return F, J, E, S
     if True:  # implicit midpoint
+        theta_p = Constant(1. / 2.)
+        theta_m = Constant(1. / 2.)
         F_f, J_f, E_f, S_f = compute_elast_tensors(dim, (u_f + u_f_0) / 2, lambda_s, mu_s)
         F_s, J_s, E_s, S_s = compute_elast_tensors(dim, (u_s + u_s_0) / 2, lambda_s, mu_s)
         v_f_mid = (v_f + v_f_0) / 2
@@ -681,38 +683,65 @@ def compute_elast_tensors(dim, u, lambda_s, mu_s):
         inner(dot(- p_mid * Identity(dim) + rho_f * nu_f * 2 * sym(dot(grad(v_f_mid), inv(F_f))), dot(J_f * transpose(inv(F_f)), n_f)), dv_s('|')) * ds_s(label_interface)
         #inner(dot(- p('|') * Identity(dim) + rho_f * nu_f * 2 * sym(dot(grad(v_f('|')), inv(F_f))), dot(J_f * transpose(inv(F_f)), n_f)), dv_s('|')) * ds_s(label_interface)
     else:  # CN
-        F_f, J_f, E_f, S_f = compute_elast_tensors(dim, u_f, lambda_s, mu_s)
-        F_s, J_s, E_s, S_s = compute_elast_tensors(dim, u_s, lambda_s, mu_s)
-        F_f_0, J_f_0, E_f_0, S_f_0 = compute_elast_tensors(dim, u_f_0, lambda_s, mu_s)
-        F_s_0, J_s_0, E_s_0, S_s_0 = compute_elast_tensors(dim, u_s_0, lambda_s, mu_s)
-        theta_p = Constant(1. / 2.)
-        theta_m = Constant(1. / 2.)
+        #theta_p = Constant(1. / 2. + 100 * float(dt))
+        #theta_m = Constant(1. / 2. - 100 * float(dt))
+        theta_p = Constant(1.)
+        theta_m = Constant(0.)
         v_f_dot = (v_f - v_f_0) / dt
         u_f_dot = (u_f - u_f_0) / dt
         v_s_dot = (v_s - v_s_0) / dt
         u_s_dot = (u_s - u_s_0) / dt
-        def _fluid(v_f, u_f, F_f, J_f, p):
-            return (inner(rho_f * J_f * vdot, dv_f) +
+        def _fluid(v_f, u_f, p):
+            F_f, J_f, E_f, S_f = compute_elast_tensors(dim, u_f, lambda_s, mu_s)
+            return (inner(rho_f * J_f * v_f_dot, dv_f) +
                     inner(rho_f * J_f * dot(dot(grad(v_f), inv(F_f)), v_f - u_f_dot), dv_f) +
                     inner(rho_f * J_f * nu_f * 2 * sym(dot(grad(v_f), inv(F_f))), dot(grad(dv_f), inv(F_f))) +
                     J_f * inner(p, tr(dot(grad(dv_f), inv(F_f)))) -
                     J_f * inner(tr(dot(grad(v_f), inv(F_f))), dp) +
                     J_f * inner(dot(grad(u_f), inv(F_f)), dot(grad(du_f), inv(F_f)))) * dx_f
-        def _struct(F_s, S_s, J_s):
+        def _struct(v_f, u_f, p, v_s, u_s):
+            F_f, J_f, E_f, S_f = compute_elast_tensors(dim, u_f, lambda_s, mu_s)
+            F_s, J_s, E_s, S_s = compute_elast_tensors(dim, u_s, lambda_s, mu_s)
             return (inner(rho_s * J_s * v_s_dot, dv_s) +
                     inner(dot(F_s, S_s), grad(dv_s)) -
                     inner(rho_s * J_s * as_vector([0., - g_s]), dv_s) +
-                    inner(J_s * (u_s - u_s_0) / dt - v_s_mid, du_s)) * dx_s
-        residual_f = theta_p * _fluid(v_f, u_f, F_f, J_f, p) + \
-                     theta_m * _fluid(v_f_0, u_f_0, F_f_0, J_f_0, p_0)
-        residual_s = theta_p * _struct() + \
-                     theta_m * _struct() + \
-                     inner(dot(- p_mid * Identity(dim) + rho_f * nu_f * 2 * sym(dot(grad(v_f_mid), inv(F_f))), dot(J_f * transpose(inv(F_f)), n_f)), dv_s('|')) * ds_s(label_interface)
-        #inner(dot(- p('|') * Identity(dim) + rho_f * nu_f * 2 * sym(dot(grad(v_f('|')), inv(F_f))), dot(J_f * transpose(inv(F_f)), n_f)), dv_s('|')) * ds_s(label_interface)
+                    inner(J_s * (u_s_dot - v_s), du_s)) * dx_s + \
+                   inner(dot(- p('|') * Identity(dim) + rho_f * nu_f * 2 * sym(dot(grad(v_f('|')), inv(F_f))), dot(J_f * transpose(inv(F_f)), n_f)), dv_s('|')) * ds_s(label_interface)
+        #residual_f = theta_p * _fluid(v_f, u_f, p) + \
+        #             theta_m * _fluid(v_f_0, u_f_0, p_0)
+        #residual_s = theta_p * _struct(v_f, u_f, p, v_s, u_s) + \
+        #             theta_m * _struct(v_f_0, u_f_0, p_0, v_s_0, u_s_0)
+        F_f, J_f, E_f, S_f = compute_elast_tensors(dim, u_f, lambda_s, mu_s)
+        F_s, J_s, E_s, S_s = compute_elast_tensors(dim, u_s, lambda_s, mu_s)
+        v_f_mid = v_f
+        v_s_mid = v_s
+        u_f_mid = u_f
+        p_mid = p
+        """
+        residual_f = (
+            inner(rho_f * J_f * (v_f - v_f_0) / dt, dv_f) +
+            inner(rho_f * J_f * dot(dot(grad(v_f_mid), inv(F_f)), v_f_mid - (u_f - u_f_0) / dt), dv_f) +
+            inner(rho_f * J_f * nu_f * 2 * sym(dot(grad(v_f_mid), inv(F_f))), dot(grad(dv_f), inv(F_f))) -
+            J_f * inner(p_mid, tr(dot(grad(dv_f), inv(F_f)))) +
+            J_f * inner(tr(dot(grad(v_f_mid), inv(F_f))), dp) +
+            J_f * inner(dot(grad(u_f_mid), inv(F_f)), dot(grad(du_f), inv(F_f)))
+        ) * dx_f
+        residual_s = (
+            inner(rho_s * J_s * (v_s - v_s_0) / dt, dv_s) +
+            inner(dot(F_s, S_s), grad(dv_s)) -
+            inner(rho_s * J_s * as_vector([0., - g_s]), dv_s) +
+            inner(J_s * ((u_s - u_s_0) / dt - v_s_mid), du_s)
+        ) * dx_s + \
+        inner(dot(- p('|') * Identity(dim) + rho_f * nu_f * 2 * sym(dot(grad(v_f('|')), inv(F_f))), dot(J_f * transpose(inv(F_f)), n_f)), dv_s('|')) * ds_s(label_interface)
+        #inner(dot(- p_mid * Identity(dim) + rho_f * nu_f * 2 * sym(dot(grad(v_f_mid), inv(F_f))), dot(J_f * transpose(inv(F_f)), n_f)), dv_s('|')) * ds_s(label_interface)
+        """
+        residual_f = _fluid(v_f, u_f, p)
+        residual_s = _struct(v_f, u_f, p, v_s, u_s)
     residual = residual_f + residual_s
-    #v_f_left = 1.5 * Ubar * y_f * (H - y_f) / ((H / 2) ** 2) * conditional(t < 2.0, (1 - cos(pi / 2 * t)) / 2., 1.)
-    v_f_left = 1.5 * Ubar * y_f * (H - y_f) / ((H / 2) ** 2) * conditional(t < 2.0 + dt / 10., (1 - cos(pi / 2 * (t - dt / 2))) / 2., 1.)
-    bc_v_f_inflow = DirichletBC(V.sub(0), as_vector([v_f_left, 0.]), (label_left, ))
+    def v_f_left(t_):
+        return 1.5 * Ubar * y_f * (H - y_f) / ((H / 2) ** 2) * conditional(t_ < 2.0 + dt / 10., (1 - cos(pi / 2 * t_)) / 2., 1.)
+    bc_v_f_inflow = DirichletBC(V.sub(0), as_vector([theta_p * v_f_left(t - dt + theta_p * dt) +
+                                                     theta_m * v_f_left(t - dt + theta_m * dt), 0.]), (label_left, ))
     bc_v_f_zero = DirichletBC(V.sub(0), Constant((0, 0)), (label_bottom, label_top, label_circle))
     bbc_v_f_noslip = DirichletBC(V.sub(0), Constant((0, 0)), ((label_circle, label_interface), ))
     bc_v_f_noslip = EquationBC(inner(v_f - v_s, dv_f) * ds_f(label_interface) == 0, solution, label_interface, bcs=[bbc_v_f_noslip], V=V.sub(0))