Skip to content

Commit

Permalink
benchmark
Browse files Browse the repository at this point in the history
  • Loading branch information
@ksagiyam
ksagiyam committed May 13, 2024
1 parent eaaf3a8 commit 76e89ef
Show file tree
Hide file tree
Showing 2 changed files with 309 additions and 0 deletions.
309 changes: 309 additions & 0 deletions benchmark.py
View file
Original file line number Diff line number Diff line change
@@ -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')
Binary file modified mesh_s.pdf
View file
Binary file not shown.

0 comments on commit 76e89ef

Please sign in to comment.