From ec320332164d991486899d1192efa47f0a2cef2e Mon Sep 17 00:00:00 2001 From: ksagiyam Date: Tue, 14 May 2024 16:40:52 +0100 Subject: [PATCH] benchmark --- benchmark.py | 153 ++++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 145 insertions(+), 8 deletions(-) diff --git a/benchmark.py b/benchmark.py index 8ceb3a3516..9b09932da5 100644 --- 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