benchmark

firedrakeproject · May 15, 2024 · ec32033 · ec32033
1 parent efae05c
commit ec32033
Showing 1 changed file with 145 additions and 8 deletions.
diff --git a/benchmark.py b/benchmark.py
@@ -108,14 +108,14 @@ def _elevate_degree(mesh, degree):
     return mesh
 
 
-T = 4.0 # 12.0
+T = 10.0 # 12.0
 dt = .05  #0.001
 ntimesteps = int(T / dt)
 t = Constant(0.0)
 dim = 2
 h = 0.10
-degree = 3  # 2 - 4
-nref = 2  # 2 - 5 tested
+degree = 2 #3  # 2 - 4
+nref = 1 #  # 2 - 5 tested for CSM 1 and 2
 mesh = make_mesh(h)
 if nref > 0:
     mh = MeshHierarchy(mesh, nref)
@@ -140,7 +140,7 @@ def _elevate_degree(mesh, degree):
 x_s, y_s = SpatialCoordinate(mesh_s)
 n_s = FacetNormal(mesh_s)
 
-case = "CSM1"
+case = "FSI1"
 
 if case in ["CFD1", "CFD2", "CFD3"]:
     if case == "CFD1":
@@ -268,7 +268,8 @@ def _elevate_degree(mesh, degree):
             "ksp_type": "gmres",
             "pc_type": "gamg",
             "mg_levels_pc_type": "sor",
-            #'mg_levels_ksp_max_it': 3,
+            'mg_levels_ksp_max_it': 3,
+            #"pc_type": "lu",
             #"pc_factor_mat_solver_type": "mumps"
         }
         near_nullspace = VectorSpaceBasis(vecs=[assemble(Function(V).interpolate(rigid_body_mode)) \
@@ -286,11 +287,147 @@ def _elevate_degree(mesh, degree):
         # u_s.at(pointA) = [-0.00722496 -0.06629327]
         print(V.dim())
         print(u_s.at(pointA))
-        assert np.allclose(u_s.at(pointA), [-0.007187, -0.06610], rtol=1e-03)
+        #assert np.allclose(u_s.at(pointA), [-0.007187, -0.06610], rtol=1e-03)
     else:  # CSM3
-        pass
+        V0 = VectorFunctionSpace(mesh_s, "P", degree)
+        V = V0 * V0
+        solution = Function(V)
+        solution_0 = Function(V)
+        v_s, u_s = split(solution)
+        v_s_0, u_s_0 = split(solution_0)
+        dv_s, du_s = split(TestFunction(V))
+        F = Identity(dim) + grad(u_s)
+        E = 1. / 2. * (dot(transpose(F), F) - Identity(dim))
+        S = lambda_s * tr(E) * Identity(dim) + 2.0 * mu_s * E
+        residual = inner((u_s - u_s_0) / dt - v_s, dv_s) * dx + \
+                   rho_s * inner(v_s - v_s_0, du_s) / dt * dx + \
+                   inner(dot(F, S), grad(du_s)) * dx - \
+                   rho_s * inner(as_vector([0., - g_s_float]), du_s) * dx
+        bc_v = DirichletBC(V.sub(0), as_vector([0., 0.]), label_struct_base)
+        bc_u = DirichletBC(V.sub(1), as_vector([0., 0.]), label_struct_base)
+        solver_parameters = {
+            "mat_type": "aij",
+            "snes_monitor": None,
+            "ksp_monitor": None,
+            #"ksp_view": None,
+            "ksp_type": "gmres",
+            #"pc_type": "gamg",
+            #"mg_levels_pc_type": "sor",
+            #'mg_levels_ksp_max_it': 3,
+            "pc_type": "lu",
+            "pc_factor_mat_solver_type": "mumps"
+        }
+        problem = NonlinearVariationalProblem(residual, solution, bcs=[bc_v, bc_u])
+        solver = NonlinearVariationalSolver(problem, solver_parameters=solver_parameters)
+        for itimestep in range(ntimesteps):
+            print(f"{itimestep} / {ntimesteps}", flush=True)
+            t.assign((itimestep + 1) * dt)
+            solver.solve()
+            for subfunction, subfunction_0 in zip(solution.subfunctions, solution_0.subfunctions):
+                subfunction_0.assign(subfunction)
+        print(solution.subfunctions[1].at(pointA))
 elif case in ["FSI1", "FSI2", "FSI3"]:
-    pass
+    if case == "FSI1":
+        rho_s = 1.e+3
+        nu_s = 0.4
+        mu_s = 0.5 * 1.e+6
+        rho_f = 1.e+3
+        nu_f = 1.e-3
+        Ubar = 0.2
+        # Re = 20.
+    elif case == "FSI2":
+        rho_s = 10. * 1.e+3
+        nu_s = 0.4
+        mu_s = 0.5 * 1.e+6
+        rho_f = 1.e+3
+        nu_f = 1.e-3
+        Ubar = 1.0
+        # Re = 100.
+    elif case == "FSI3":
+        rho_s = 1. * 1.e+3
+        nu_s = 0.4
+        mu_s = 2.0 * 1.e+6
+        rho_f = 1.e+3
+        nu_f = 1.e-3
+        Ubar = 2.0
+        # Re = 200.
+    else:
+        raise ValueError
+    g_s = 2.0
+    E_s = mu_s * 2 * (1 + nu_s)
+    lambda_s = nu_s * E_s / (1 + nu_s) / (1 - 2 * nu_s)
+    V_0 = VectorFunctionSpace(mesh_f, "P", degree)
+    V_1 = FunctionSpace(mesh_f, "P", degree - 2)
+    V = V_0 * V_1
+    solution = Function(V)
+    solution_0 = Function(V)
+    v_f, p_f = split(solution)
+    u_f = Function(V_0)
+    v_f_0, p_f_0 = split(solution_0)
+    dv_f, dp_f = split(TestFunction(V))
+    residual = (
+                rho_f * inner(v_f - v_f_0, dv_f) / dt +
+                rho_f * inner(dot(grad(v_f), v_f), dv_f) +
+                rho_f * nu_f * inner(2 * sym(grad(v_f)), grad(dv_f)) -
+                inner(p_f, div(dv_f)) +
+                inner(div(v_f), dp_f)) * dx
+    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.)
+    bc_inflow = DirichletBC(V.sub(0), as_vector([v_f_left, 0.]), label_left)
+    bc_noslip = DirichletBC(V.sub(0), Constant((0, 0)), label_bottom + label_top + label_circle + label_interface)
+    #bc_inflow = EquationBC(inner(v_f - as_vector([v_f_left, 0.]), dv_f) * ds(label_left) == 0, solution, label_left, V=V.sub(0))
+    #bc_noslip = EquationBC(inner(v_f, dv_f) * ds(label_bottom + label_top + label_circle + label_interface) == 0, solution, label_bottom + label_top + label_circle + label_interface, V=V.sub(0))
+    solver_parameters = {
+        "mat_type": "aij",
+        "snes_monitor": None,
+        "ksp_type": "gmres",
+        "pc_type": "lu",
+        "pc_factor_mat_solver_type": "mumps"
+    }
+    """
+    solver_parameters = {
+        #'mat_type': 'matfree',
+        'pc_type': 'fieldsplit',
+        'ksp_type': 'preonly',
+        'pc_fieldsplit_type': 'schur',
+        'fieldsplit_schur_fact_type': 'full',
+        'fieldsplit_0': {
+            #'ksp_type': 'cg',
+            'ksp_type': 'gmres',  # equationBC is nonsym
+            'pc_type': 'python',
+            'pc_python_type': 'firedrake.AssembledPC',
+            'assembled_pc_type': 'gamg',
+            'assembled_mg_levels_pc_type': 'sor',
+            'assembled_mg_levels_pc_sor_diagonal_shift': 1e-100,  # See https://gitlab.com/petsc/petsc/-/issues/1221
+            'ksp_rtol': 1e-7,
+            'ksp_converged_reason': None,
+        },
+        'fieldsplit_1': {
+            'ksp_type': 'fgmres',
+            'ksp_converged_reason': None,
+            'pc_type': 'python',
+            'pc_python_type': 'firedrake.MassInvPC',
+            'Mp_pc_type': 'ksp',
+            'Mp_ksp_ksp_type': 'cg',
+            'Mp_ksp_pc_type': 'sor',
+            'ksp_rtol': '1e-5',
+            'ksp_monitor': None,
+        },
+        "snes_monitor": None,
+    }
+    """
+    problem = NonlinearVariationalProblem(residual, solution, bcs=[bc_inflow, bc_noslip])
+    solver = NonlinearVariationalSolver(problem, solver_parameters=solver_parameters)
+    for itimestep in range(ntimesteps):
+        print(f"{itimestep} / {ntimesteps}", flush=True)
+        t.assign((itimestep + 1) * dt)
+        solver.solve()
+        for subfunction, subfunction_0 in zip(solution.subfunctions, solution_0.subfunctions):
+            subfunction_0.assign(subfunction)
+    FD = assemble((-p_f * n_f + rho_f * nu_f * dot(2 * sym(grad(v_f)), n_f))[0] * ds(label_circle + label_interface))
+    FL = assemble((-p_f * n_f + rho_f * nu_f * dot(2 * sym(grad(v_f)), n_f))[1] * ds(label_circle + label_interface))
+    print(f"FD = {FD}")
+    print(f"FL = {FL}")
+    print("num cells = ", mesh_f.comm.allreduce(mesh_f.cell_set.size))
 else:
     raise ValueError