diff --git a/benchmark.py b/benchmark.py new file mode 100644 index 0000000000..8ceb3a3516 --- /dev/null +++ b/benchmark.py @@ -0,0 +1,309 @@ +import numpy as np +import matplotlib.pyplot as plt +from firedrake import * +from firedrake.pyplot import triplot, plot, quiver +from firedrake.cython import dmcommon + + +L = 2.50 +H = 0.41 +Cx = 0.20 # center of hole: x +Cy = 0.20 # y +r = 0.05 +pointA = (0.6, 0.2) +pointB = (0.2, 0.2) +label_fluid = 1 +label_struct = 2 +label_left = (4, ) +label_right = (2, ) +label_bottom = (1, ) +label_top = (3, ) +label_interface = (10, 11, 12) +label_struct_base = (6, ) +label_circle = (7, 8, 9, 5) + + +def make_mesh(h): + # points + # 3 2 + # +--------------------------------------+ + # | __ 9 | + # | 11/10 \8___________15 | + # | 12+ +7__________| | + # | 13\__ /6 14 | + # | 4 5 | + # +--------------------------------------+ + # 0 1 + # labels + # 3 + # +--------------------------------------+ + # | __ 7 | + # | 8 / \ ____12_____ | + # 4 | | +6__________|11 | 2 + # | 9 \__ / 10 | + # | 5 | + # +--------------------------------------+ + # 1 + import netgen + from netgen.geom2d import SplineGeometry + comm = COMM_WORLD + if comm.Get_rank() == 0: + geom = SplineGeometry() + pnts = [(0, 0), # 0 + (L, 0), # 1 + (L, H), # 2 + (0, H), # 3 + (0.20, 0.15), # 4 + (0.240824829046386, 0.15), # 5 + (0.248989794855664, 0.19), # 6 + (0.25, 0.20), # 7 + (0.248989794855664, 0.21), # 8 + (0.240824829046386, 0.25), # 9 + (0.20, 0.25), # 10 + (0.15, 0.25), # 11 + (0.15, 0.20), # 12 + (0.15, 0.15), # 13 + (0.60, 0.19), # 14 + (0.60, 0.21), # 15 + (0.55, 0.19), # 16 + (0.56, 0.15), # 17 + (0.60, 0.15), # 18 + (0.65, 0.15), # 19 + (0.65, 0.20), # 20 + (0.65, 0.25), # 21 + (0.60, 0.25), # 22 + (0.56, 0.25), # 23 + (0.55, 0.21), # 24 + (0.65, 0.25), # 25 + (0.65, 0.15)] # 26 + pind = [geom.AppendPoint(*pnt) for pnt in pnts] + geom.Append(['line', pind[0], pind[1]], leftdomain=1, rightdomain=0, bc="wall") + geom.Append(['line', pind[1], pind[2]], leftdomain=1, rightdomain=0, bc="outlet") + geom.Append(['line', pind[2], pind[3]], leftdomain=1, rightdomain=0, bc="wall") + geom.Append(['line', pind[3], pind[0]], leftdomain=1, rightdomain=0, bc="inlet") + geom.Append(['spline3', pind[4], pind[5], pind[6]], leftdomain=0, rightdomain=1, bc="circ") + geom.Append(['spline3', pind[6], pind[7], pind[8]], leftdomain=0, rightdomain=2, bc="circ2") + geom.Append(['spline3', pind[8], pind[9], pind[10]], leftdomain=0, rightdomain=1, bc="circ") + geom.Append(['spline3', pind[10], pind[11], pind[12]], leftdomain=0, rightdomain=1, bc="circ") + geom.Append(['spline3', pind[12], pind[13], pind[4]], leftdomain=0, rightdomain=1, bc="circ") + geom.Append(['line', pind[6], pind[14]], leftdomain=2, rightdomain=1, bc="interface") + geom.Append(['line', pind[14], pind[15]], leftdomain=2, rightdomain=1, bc="interface") + geom.Append(['line', pind[15], pind[8]], leftdomain=2, rightdomain=1, bc="interface") + geom.SetMaterial(1, "fluid") + geom.SetMaterial(2, "solid") + ngmesh = geom.GenerateMesh(maxh=h) + else: + ngmesh = netgen.libngpy._meshing.Mesh(2) + return Mesh(ngmesh) + + +def _elevate_degree(mesh, degree): + V = VectorFunctionSpace(mesh, "CG", degree) + f = Function(V).interpolate(mesh.coordinates) + bc = DirichletBC(V, 0., label_circle + label_struct_base) + r_ = np.sqrt((f.dat.data_with_halos[bc.nodes, 0] - Cx) ** 2 + ((f.dat.data_with_halos[bc.nodes, 1] - Cy) ** 2)) + f.dat.data_with_halos[bc.nodes, 0] = (f.dat.data_with_halos[bc.nodes, 0] - Cx) * (r / r_) + Cy + f.dat.data_with_halos[bc.nodes, 1] = (f.dat.data_with_halos[bc.nodes, 1] - Cy) * (r / r_) + Cy + mesh = Mesh(f) + return mesh + + +T = 4.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 +mesh = make_mesh(h) +if nref > 0: + mh = MeshHierarchy(mesh, nref) + mesh = mh[-1] +mesh = _elevate_degree(mesh, degree) +""" +#mesh.topology_dm.viewFromOptions("-dm_view") +x, y = SpatialCoordinate(mesh) +v = assemble(Constant(1.0, domain=mesh) * ds(label_circle + label_struct_base)) +print(v - 2 * pi * r) +print(assemble(x * dx(label_struct))) +print((0.6 - 0.248989794855664) * (0.6 + 0.248989794855664) /2. * 0.02) +print(assemble(Constant(1) * dx(domain=mesh, subdomain_id=label_fluid)) + assemble(Constant(1) * dx(domain=mesh, subdomain_id=label_struct))) +print(assemble(Constant(1) * dx(domain=mesh))) +print(L * H - pi * r ** 2) +""" + +mesh_f = Submesh(mesh, dmcommon.CELL_SETS_LABEL, label_fluid, mesh.topological_dimension()) +x_f, y_f = SpatialCoordinate(mesh_f) +n_f = FacetNormal(mesh_f) +mesh_s = Submesh(mesh, dmcommon.CELL_SETS_LABEL, label_struct, mesh.topological_dimension()) +x_s, y_s = SpatialCoordinate(mesh_s) +n_s = FacetNormal(mesh_s) + +case = "CSM1" + +if case in ["CFD1", "CFD2", "CFD3"]: + if case == "CFD1": + rho_f = 1.e+3 + nu_f = 1.e-3 + Ubar = 0.2 + # Re = 20. + elif case == "CFD2": + rho_f = 1.e+3 + nu_f = 1.e-3 + Ubar = 1. + # Re = 100. + elif case == "CFD3": + rho_f = 1.e+3 + nu_f = 1.e-3 + Ubar = 2. + # Re = 200. + else: + raise ValueError + 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) + 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)) +elif case in ["CSM1", "CSM2", "CSM3"]: + if case == "CSM1": + rho_s = 1.e+3 + nu_s = 0.4 + E_s = 1.4 * 1.e+6 + elif case == "CSM2": + rho_s = 1.e+3 + nu_s = 0.4 + E_s = 5.6 * 1.e+6 + elif case == "CSM3": + rho_s = 1.e+3 + nu_s = 0.4 + E_s = 1.4 * 1.e+6 + else: + raise ValueError + g_s_float = 2.0 + g_s = Constant(0.) + lambda_s = nu_s * E_s / (1 + nu_s) / (1 - 2 * nu_s) + mu_s = E_s / 2 / (1 + nu_s) + if case in ["CSM1", "CSM2"]: + V = VectorFunctionSpace(mesh_s, "P", degree) + u_s = Function(V) + du_s = 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(dot(F, S), grad(du_s)) * dx - \ + rho_s * inner(as_vector([0., - g_s]), du_s) * dx + bc = DirichletBC(V, 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_factor_mat_solver_type": "mumps" + } + near_nullspace = VectorSpaceBasis(vecs=[assemble(Function(V).interpolate(rigid_body_mode)) \ + for rigid_body_mode in [Constant((1, 0)), Constant((0, 1)), as_vector([y_s, -x_s])]]) + near_nullspace.orthonormalize() + problem = NonlinearVariationalProblem(residual, u_s, bcs=[bc]) + solver = NonlinearVariationalSolver(problem, solver_parameters=solver_parameters, near_nullspace=near_nullspace) + # Use relaxation method. + nsteps = 10 + for g_s_temp in [g_s_float * (i + 1) / nsteps for i in range(nsteps)]: + g_s.assign(g_s_temp) + solver.solve() + # (degree, nref) = (2-4, 2-4) with mumps work. .1 * gs + # (4, 5) has 280,962 DoFs. + # 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) + else: # CSM3 + pass +elif case in ["FSI1", "FSI2", "FSI3"]: + pass +else: + raise ValueError + +if mesh.comm.size == 1: + fig, axes = plt.subplots() + axes.axis('equal') + #quiver(solution.subfunctions[0]) + triplot(mesh, axes=axes) + axes.legend() + plt.savefig('mesh_orig.pdf') + fig, axes = plt.subplots() + axes.axis('equal') + #quiver(solution.subfunctions[0]) + triplot(mesh_s, axes=axes) + axes.legend() + plt.savefig('mesh_s.pdf') diff --git a/mesh_s.pdf b/mesh_s.pdf index ed78bf50ac..7da77a1993 100644 Binary files a/mesh_s.pdf and b/mesh_s.pdf differ