Skip to content

Commit efae05c

Browse files
ksagiyamksagiyam
committed
benchmark
1 parent 192d54a commit efae05c

File tree

2 files changed

+309
-0
lines changed

2 files changed

+309
-0
lines changed

‎benchmark.py‎

Lines changed: 309 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,309 @@
1+
import numpy as np
2+
import matplotlib.pyplot as plt
3+
from firedrake import *
4+
from firedrake.pyplot import triplot, plot, quiver
5+
from firedrake.cython import dmcommon
6+
7+
8+
L = 2.50
9+
H = 0.41
10+
Cx = 0.20 # center of hole: x
11+
Cy = 0.20 # y
12+
r = 0.05
13+
pointA = (0.6, 0.2)
14+
pointB = (0.2, 0.2)
15+
label_fluid = 1
16+
label_struct = 2
17+
label_left = (4, )
18+
label_right = (2, )
19+
label_bottom = (1, )
20+
label_top = (3, )
21+
label_interface = (10, 11, 12)
22+
label_struct_base = (6, )
23+
label_circle = (7, 8, 9, 5)
24+
25+
26+
def make_mesh(h):
27+
# points
28+
# 3 2
29+
# +--------------------------------------+
30+
# | __ 9 |
31+
# | 11/10 \8___________15 |
32+
# | 12+ +7__________| |
33+
# | 13\__ /6 14 |
34+
# | 4 5 |
35+
# +--------------------------------------+
36+
# 0 1
37+
# labels
38+
# 3
39+
# +--------------------------------------+
40+
# | __ 7 |
41+
# | 8 / \ ____12_____ |
42+
# 4 | | +6__________|11 | 2
43+
# | 9 \__ / 10 |
44+
# | 5 |
45+
# +--------------------------------------+
46+
# 1
47+
import netgen
48+
from netgen.geom2d import SplineGeometry
49+
comm = COMM_WORLD
50+
if comm.Get_rank() == 0:
51+
geom = SplineGeometry()
52+
pnts = [(0, 0), # 0
53+
(L, 0), # 1
54+
(L, H), # 2
55+
(0, H), # 3
56+
(0.20, 0.15), # 4
57+
(0.240824829046386, 0.15), # 5
58+
(0.248989794855664, 0.19), # 6
59+
(0.25, 0.20), # 7
60+
(0.248989794855664, 0.21), # 8
61+
(0.240824829046386, 0.25), # 9
62+
(0.20, 0.25), # 10
63+
(0.15, 0.25), # 11
64+
(0.15, 0.20), # 12
65+
(0.15, 0.15), # 13
66+
(0.60, 0.19), # 14
67+
(0.60, 0.21), # 15
68+
(0.55, 0.19), # 16
69+
(0.56, 0.15), # 17
70+
(0.60, 0.15), # 18
71+
(0.65, 0.15), # 19
72+
(0.65, 0.20), # 20
73+
(0.65, 0.25), # 21
74+
(0.60, 0.25), # 22
75+
(0.56, 0.25), # 23
76+
(0.55, 0.21), # 24
77+
(0.65, 0.25), # 25
78+
(0.65, 0.15)] # 26
79+
pind = [geom.AppendPoint(*pnt) for pnt in pnts]
80+
geom.Append(['line', pind[0], pind[1]], leftdomain=1, rightdomain=0, bc="wall")
81+
geom.Append(['line', pind[1], pind[2]], leftdomain=1, rightdomain=0, bc="outlet")
82+
geom.Append(['line', pind[2], pind[3]], leftdomain=1, rightdomain=0, bc="wall")
83+
geom.Append(['line', pind[3], pind[0]], leftdomain=1, rightdomain=0, bc="inlet")
84+
geom.Append(['spline3', pind[4], pind[5], pind[6]], leftdomain=0, rightdomain=1, bc="circ")
85+
geom.Append(['spline3', pind[6], pind[7], pind[8]], leftdomain=0, rightdomain=2, bc="circ2")
86+
geom.Append(['spline3', pind[8], pind[9], pind[10]], leftdomain=0, rightdomain=1, bc="circ")
87+
geom.Append(['spline3', pind[10], pind[11], pind[12]], leftdomain=0, rightdomain=1, bc="circ")
88+
geom.Append(['spline3', pind[12], pind[13], pind[4]], leftdomain=0, rightdomain=1, bc="circ")
89+
geom.Append(['line', pind[6], pind[14]], leftdomain=2, rightdomain=1, bc="interface")
90+
geom.Append(['line', pind[14], pind[15]], leftdomain=2, rightdomain=1, bc="interface")
91+
geom.Append(['line', pind[15], pind[8]], leftdomain=2, rightdomain=1, bc="interface")
92+
geom.SetMaterial(1, "fluid")
93+
geom.SetMaterial(2, "solid")
94+
ngmesh = geom.GenerateMesh(maxh=h)
95+
else:
96+
ngmesh = netgen.libngpy._meshing.Mesh(2)
97+
return Mesh(ngmesh)
98+
99+
100+
def _elevate_degree(mesh, degree):
101+
V = VectorFunctionSpace(mesh, "CG", degree)
102+
f = Function(V).interpolate(mesh.coordinates)
103+
bc = DirichletBC(V, 0., label_circle + label_struct_base)
104+
r_ = np.sqrt((f.dat.data_with_halos[bc.nodes, 0] - Cx) ** 2 + ((f.dat.data_with_halos[bc.nodes, 1] - Cy) ** 2))
105+
f.dat.data_with_halos[bc.nodes, 0] = (f.dat.data_with_halos[bc.nodes, 0] - Cx) * (r / r_) + Cy
106+
f.dat.data_with_halos[bc.nodes, 1] = (f.dat.data_with_halos[bc.nodes, 1] - Cy) * (r / r_) + Cy
107+
mesh = Mesh(f)
108+
return mesh
109+
110+
111+
T = 4.0 # 12.0
112+
dt = .05 #0.001
113+
ntimesteps = int(T / dt)
114+
t = Constant(0.0)
115+
dim = 2
116+
h = 0.10
117+
degree = 3 # 2 - 4
118+
nref = 2 # 2 - 5 tested
119+
mesh = make_mesh(h)
120+
if nref > 0:
121+
mh = MeshHierarchy(mesh, nref)
122+
mesh = mh[-1]
123+
mesh = _elevate_degree(mesh, degree)
124+
"""
125+
#mesh.topology_dm.viewFromOptions("-dm_view")
126+
x, y = SpatialCoordinate(mesh)
127+
v = assemble(Constant(1.0, domain=mesh) * ds(label_circle + label_struct_base))
128+
print(v - 2 * pi * r)
129+
print(assemble(x * dx(label_struct)))
130+
print((0.6 - 0.248989794855664) * (0.6 + 0.248989794855664) /2. * 0.02)
131+
print(assemble(Constant(1) * dx(domain=mesh, subdomain_id=label_fluid)) + assemble(Constant(1) * dx(domain=mesh, subdomain_id=label_struct)))
132+
print(assemble(Constant(1) * dx(domain=mesh)))
133+
print(L * H - pi * r ** 2)
134+
"""
135+
136+
mesh_f = Submesh(mesh, dmcommon.CELL_SETS_LABEL, label_fluid, mesh.topological_dimension())
137+
x_f, y_f = SpatialCoordinate(mesh_f)
138+
n_f = FacetNormal(mesh_f)
139+
mesh_s = Submesh(mesh, dmcommon.CELL_SETS_LABEL, label_struct, mesh.topological_dimension())
140+
x_s, y_s = SpatialCoordinate(mesh_s)
141+
n_s = FacetNormal(mesh_s)
142+
143+
case = "CSM1"
144+
145+
if case in ["CFD1", "CFD2", "CFD3"]:
146+
if case == "CFD1":
147+
rho_f = 1.e+3
148+
nu_f = 1.e-3
149+
Ubar = 0.2
150+
# Re = 20.
151+
elif case == "CFD2":
152+
rho_f = 1.e+3
153+
nu_f = 1.e-3
154+
Ubar = 1.
155+
# Re = 100.
156+
elif case == "CFD3":
157+
rho_f = 1.e+3
158+
nu_f = 1.e-3
159+
Ubar = 2.
160+
# Re = 200.
161+
else:
162+
raise ValueError
163+
V_0 = VectorFunctionSpace(mesh_f, "P", degree)
164+
V_1 = FunctionSpace(mesh_f, "P", degree - 2)
165+
V = V_0 * V_1
166+
solution = Function(V)
167+
solution_0 = Function(V)
168+
v_f, p_f = split(solution)
169+
v_f_0, p_f_0 = split(solution_0)
170+
dv_f, dp_f = split(TestFunction(V))
171+
residual = (
172+
rho_f * inner(v_f - v_f_0, dv_f) / dt +
173+
rho_f * inner(dot(grad(v_f), v_f), dv_f) +
174+
rho_f * nu_f * inner(2 * sym(grad(v_f)), grad(dv_f)) -
175+
inner(p_f, div(dv_f)) +
176+
inner(div(v_f), dp_f)) * dx
177+
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.)
178+
bc_inflow = DirichletBC(V.sub(0), as_vector([v_f_left, 0.]), label_left)
179+
#bc_noslip = DirichletBC(V.sub(0), Constant((0, 0)), label_bottom + label_top + label_circle + label_interface)
180+
#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))
181+
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))
182+
solver_parameters = {
183+
"mat_type": "aij",
184+
"snes_monitor": None,
185+
"ksp_type": "gmres",
186+
"pc_type": "lu",
187+
"pc_factor_mat_solver_type": "mumps"
188+
}
189+
"""
190+
solver_parameters = {
191+
#'mat_type': 'matfree',
192+
'pc_type': 'fieldsplit',
193+
'ksp_type': 'preonly',
194+
'pc_fieldsplit_type': 'schur',
195+
'fieldsplit_schur_fact_type': 'full',
196+
'fieldsplit_0': {
197+
#'ksp_type': 'cg',
198+
'ksp_type': 'gmres', # equationBC is nonsym
199+
'pc_type': 'python',
200+
'pc_python_type': 'firedrake.AssembledPC',
201+
'assembled_pc_type': 'gamg',
202+
'assembled_mg_levels_pc_type': 'sor',
203+
'assembled_mg_levels_pc_sor_diagonal_shift': 1e-100, # See https://gitlab.com/petsc/petsc/-/issues/1221
204+
'ksp_rtol': 1e-7,
205+
'ksp_converged_reason': None,
206+
},
207+
'fieldsplit_1': {
208+
'ksp_type': 'fgmres',
209+
'ksp_converged_reason': None,
210+
'pc_type': 'python',
211+
'pc_python_type': 'firedrake.MassInvPC',
212+
'Mp_pc_type': 'ksp',
213+
'Mp_ksp_ksp_type': 'cg',
214+
'Mp_ksp_pc_type': 'sor',
215+
'ksp_rtol': '1e-5',
216+
'ksp_monitor': None,
217+
},
218+
"snes_monitor": None,
219+
}
220+
"""
221+
problem = NonlinearVariationalProblem(residual, solution, bcs=[bc_inflow, bc_noslip])
222+
solver = NonlinearVariationalSolver(problem, solver_parameters=solver_parameters)
223+
for itimestep in range(ntimesteps):
224+
print(f"{itimestep} / {ntimesteps}", flush=True)
225+
t.assign((itimestep + 1) * dt)
226+
solver.solve()
227+
for subfunction, subfunction_0 in zip(solution.subfunctions, solution_0.subfunctions):
228+
subfunction_0.assign(subfunction)
229+
FD = assemble((-p_f * n_f + rho_f * nu_f * dot(2 * sym(grad(v_f)), n_f))[0] * ds(label_circle + label_interface))
230+
FL = assemble((-p_f * n_f + rho_f * nu_f * dot(2 * sym(grad(v_f)), n_f))[1] * ds(label_circle + label_interface))
231+
print(f"FD = {FD}")
232+
print(f"FL = {FL}")
233+
print("num cells = ", mesh_f.comm.allreduce(mesh_f.cell_set.size))
234+
elif case in ["CSM1", "CSM2", "CSM3"]:
235+
if case == "CSM1":
236+
rho_s = 1.e+3
237+
nu_s = 0.4
238+
E_s = 1.4 * 1.e+6
239+
elif case == "CSM2":
240+
rho_s = 1.e+3
241+
nu_s = 0.4
242+
E_s = 5.6 * 1.e+6
243+
elif case == "CSM3":
244+
rho_s = 1.e+3
245+
nu_s = 0.4
246+
E_s = 1.4 * 1.e+6
247+
else:
248+
raise ValueError
249+
g_s_float = 2.0
250+
g_s = Constant(0.)
251+
lambda_s = nu_s * E_s / (1 + nu_s) / (1 - 2 * nu_s)
252+
mu_s = E_s / 2 / (1 + nu_s)
253+
if case in ["CSM1", "CSM2"]:
254+
V = VectorFunctionSpace(mesh_s, "P", degree)
255+
u_s = Function(V)
256+
du_s = TestFunction(V)
257+
F = Identity(dim) + grad(u_s)
258+
E = 1. / 2. * (dot(transpose(F), F) - Identity(dim))
259+
S = lambda_s * tr(E) * Identity(dim) + 2.0 * mu_s * E
260+
residual = inner(dot(F, S), grad(du_s)) * dx - \
261+
rho_s * inner(as_vector([0., - g_s]), du_s) * dx
262+
bc = DirichletBC(V, as_vector([0., 0.]), label_struct_base)
263+
solver_parameters = {
264+
"mat_type": "aij",
265+
"snes_monitor": None,
266+
"ksp_monitor": None,
267+
#"ksp_view": None,
268+
"ksp_type": "gmres",
269+
"pc_type": "gamg",
270+
"mg_levels_pc_type": "sor",
271+
#'mg_levels_ksp_max_it': 3,
272+
#"pc_factor_mat_solver_type": "mumps"
273+
}
274+
near_nullspace = VectorSpaceBasis(vecs=[assemble(Function(V).interpolate(rigid_body_mode)) \
275+
for rigid_body_mode in [Constant((1, 0)), Constant((0, 1)), as_vector([y_s, -x_s])]])
276+
near_nullspace.orthonormalize()
277+
problem = NonlinearVariationalProblem(residual, u_s, bcs=[bc])
278+
solver = NonlinearVariationalSolver(problem, solver_parameters=solver_parameters, near_nullspace=near_nullspace)
279+
# Use relaxation method.
280+
nsteps = 10
281+
for g_s_temp in [g_s_float * (i + 1) / nsteps for i in range(nsteps)]:
282+
g_s.assign(g_s_temp)
283+
solver.solve()
284+
# (degree, nref) = (2-4, 2-4) with mumps work. .1 * gs
285+
# (4, 5) has 280,962 DoFs.
286+
# u_s.at(pointA) = [-0.00722496 -0.06629327]
287+
print(V.dim())
288+
print(u_s.at(pointA))
289+
assert np.allclose(u_s.at(pointA), [-0.007187, -0.06610], rtol=1e-03)
290+
else: # CSM3
291+
pass
292+
elif case in ["FSI1", "FSI2", "FSI3"]:
293+
pass
294+
else:
295+
raise ValueError
296+
297+
if mesh.comm.size == 1:
298+
fig, axes = plt.subplots()
299+
axes.axis('equal')
300+
#quiver(solution.subfunctions[0])
301+
triplot(mesh, axes=axes)
302+
axes.legend()
303+
plt.savefig('mesh_orig.pdf')
304+
fig, axes = plt.subplots()
305+
axes.axis('equal')
306+
#quiver(solution.subfunctions[0])
307+
triplot(mesh_s, axes=axes)
308+
axes.legend()
309+
plt.savefig('mesh_s.pdf')

‎mesh_s.pdf‎

49.5 KB
Binary file not shown.

0 commit comments

Comments
 (0)