Skip to content

Commit fc6545c

Browse files
ksagiyamksagiyam
committed
benchmark
1 parent 3243047 commit fc6545c

File tree

5 files changed

+91
-33
lines changed

5 files changed

+91
-33
lines changed

‎benchmark.py‎

Lines changed: 91 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -173,29 +173,30 @@ def _elevate_degree(mesh, degree):
173173
return Mesh(f)
174174

175175

176-
use_netgen = True
176+
use_netgen = False
177177
quadrilateral = True
178178

179179
T = 20 # 10.0 # 12.0
180-
dt_float = 0.004 #.002
180+
dt_float = 0.001 #.002
181181
dt = Constant(dt_float) #0.001
182182
dt_plot = 0.01
183183
ntimesteps = int(T / dt_float)
184184
t = Constant(0.0)
185185
dim = 2
186186
degree = 3 # 2 - 4
187187
if use_netgen:
188-
nref = 4 # # 2 - 5 tested for CSM 1 and 2
188+
nref = 1 # # 2 - 5 tested for CSM 1 and 2
189189
mesh = make_mesh_netgen(0.1 / 2 ** nref)
190190
mesh = _elevate_degree(mesh, degree)
191191
mesh_f = Submesh(mesh, dmcommon.CELL_SETS_LABEL, label_fluid, mesh.topological_dimension())
192192
mesh_f = _elevate_degree(mesh_f, degree)
193193
mesh_s = Submesh(mesh, dmcommon.CELL_SETS_LABEL, label_struct, mesh.topological_dimension())
194194
mesh_s = _elevate_degree(mesh_s, degree)
195195
else:
196-
nref = 1
196+
nref = 0
197197
mesh = make_mesh(quadrilateral)
198-
mesh = MeshHierarchy(mesh, nref)[-1]
198+
if nref > 0:
199+
mesh = MeshHierarchy(mesh, nref)[-1]
199200
mesh_f = Submesh(mesh, dmcommon.CELL_SETS_LABEL, label_fluid, mesh.topological_dimension())
200201
mesh_s = Submesh(mesh, dmcommon.CELL_SETS_LABEL, label_struct, mesh.topological_dimension())
201202
mesh = _finalise_mesh(mesh, degree)
@@ -656,32 +657,89 @@ def compute_elast_tensors(dim, u, lambda_s, mu_s):
656657
E = 1. / 2. * (dot(transpose(F), F) - Identity(dim))
657658
S = lambda_s * tr(E) * Identity(dim) + 2.0 * mu_s * E
658659
return F, J, E, S
659-
F_f, J_f, E_f, S_f = compute_elast_tensors(dim, (u_f + u_f_0) / 2, lambda_s, mu_s)
660-
F_s, J_s, E_s, S_s = compute_elast_tensors(dim, (u_s + u_s_0) / 2, lambda_s, mu_s)
661-
v_f_mid = (v_f + v_f_0) / 2
662-
v_s_mid = (v_s + v_s_0) / 2
663-
u_f_mid = (u_f + u_f_0) / 2
664-
p_mid = (p + p_0) / 2
665-
residual_f = (
666-
inner(rho_f * J_f * (v_f - v_f_0) / dt, dv_f) +
667-
inner(rho_f * J_f * dot(dot(grad(v_f_mid), inv(F_f)), v_f_mid - (u_f - u_f_0) / dt), dv_f) +
668-
inner(rho_f * J_f * nu_f * 2 * sym(dot(grad(v_f_mid), inv(F_f))), dot(grad(dv_f), inv(F_f))) -
669-
J_f * inner(p_mid, tr(dot(grad(dv_f), inv(F_f)))) +
670-
J_f * inner(tr(dot(grad(v_f_mid), inv(F_f))), dp) +
671-
J_f * inner(dot(grad(u_f_mid), inv(F_f)), dot(grad(du_f), inv(F_f)))
672-
) * dx_f
673-
residual_s = (
674-
inner(rho_s * J_s * (v_s - v_s_0) / dt, dv_s) +
675-
inner(dot(F_s, S_s), grad(dv_s)) -
676-
inner(rho_s * J_s * as_vector([0., - g_s]), dv_s) +
677-
inner(J_s * (u_s - u_s_0) / dt - v_s_mid, du_s)
678-
) * dx_s + \
679-
inner(dot(- p_mid * Identity(dim) + rho_f * nu_f * 2 * sym(dot(grad(v_f_mid), inv(F_f))), dot(J_f * transpose(inv(F_f)), n_f)), dv_s('|')) * ds_s(label_interface)
680-
#inner(dot(- p('|') * Identity(dim) + rho_f * nu_f * 2 * sym(dot(grad(v_f('|')), inv(F_f))), dot(J_f * transpose(inv(F_f)), n_f)), dv_s('|')) * ds_s(label_interface)
660+
if True: # implicit midpoint
661+
theta_p = Constant(1. / 2.)
662+
theta_m = Constant(1. / 2.)
663+
F_f, J_f, E_f, S_f = compute_elast_tensors(dim, (u_f + u_f_0) / 2, lambda_s, mu_s)
664+
F_s, J_s, E_s, S_s = compute_elast_tensors(dim, (u_s + u_s_0) / 2, lambda_s, mu_s)
665+
v_f_mid = (v_f + v_f_0) / 2
666+
v_s_mid = (v_s + v_s_0) / 2
667+
u_f_mid = (u_f + u_f_0) / 2
668+
p_mid = (p + p_0) / 2
669+
residual_f = (
670+
inner(rho_f * J_f * (v_f - v_f_0) / dt, dv_f) +
671+
inner(rho_f * J_f * dot(dot(grad(v_f_mid), inv(F_f)), v_f_mid - (u_f - u_f_0) / dt), dv_f) +
672+
inner(rho_f * J_f * nu_f * 2 * sym(dot(grad(v_f_mid), inv(F_f))), dot(grad(dv_f), inv(F_f))) -
673+
J_f * inner(p_mid, tr(dot(grad(dv_f), inv(F_f)))) +
674+
J_f * inner(tr(dot(grad(v_f_mid), inv(F_f))), dp) +
675+
J_f * inner(dot(grad(u_f_mid), inv(F_f)), dot(grad(du_f), inv(F_f)))
676+
) * dx_f
677+
residual_s = (
678+
inner(rho_s * J_s * (v_s - v_s_0) / dt, dv_s) +
679+
inner(dot(F_s, S_s), grad(dv_s)) -
680+
inner(rho_s * J_s * as_vector([0., - g_s]), dv_s) +
681+
inner(J_s * ((u_s - u_s_0) / dt - v_s_mid), du_s)
682+
) * dx_s + \
683+
inner(dot(- p_mid * Identity(dim) + rho_f * nu_f * 2 * sym(dot(grad(v_f_mid), inv(F_f))), dot(J_f * transpose(inv(F_f)), n_f)), dv_s('|')) * ds_s(label_interface)
684+
#inner(dot(- p('|') * Identity(dim) + rho_f * nu_f * 2 * sym(dot(grad(v_f('|')), inv(F_f))), dot(J_f * transpose(inv(F_f)), n_f)), dv_s('|')) * ds_s(label_interface)
685+
else: # CN
686+
theta_p = Constant(1. / 2. + 100 * float(dt))
687+
theta_m = Constant(1. / 2. - 100 * float(dt))
688+
#theta_p = Constant(1.)
689+
#theta_m = Constant(0.)
690+
v_f_dot = (v_f - v_f_0) / dt
691+
u_f_dot = (u_f - u_f_0) / dt
692+
v_s_dot = (v_s - v_s_0) / dt
693+
u_s_dot = (u_s - u_s_0) / dt
694+
def _fluid(v_f, u_f, p):
695+
F_f, J_f, E_f, S_f = compute_elast_tensors(dim, u_f, lambda_s, mu_s)
696+
return (inner(rho_f * J_f * v_f_dot, dv_f) +
697+
inner(rho_f * J_f * dot(dot(grad(v_f), inv(F_f)), v_f - u_f_dot), dv_f) +
698+
inner(rho_f * J_f * nu_f * 2 * sym(dot(grad(v_f), inv(F_f))), dot(grad(dv_f), inv(F_f))) -
699+
J_f * inner(p, tr(dot(grad(dv_f), inv(F_f)))) +
700+
J_f * inner(tr(dot(grad(v_f), inv(F_f))), dp) +
701+
J_f * inner(dot(grad(u_f), inv(F_f)), dot(grad(du_f), inv(F_f)))) * dx_f
702+
def _struct(v_f, u_f, p, v_s, u_s):
703+
F_f, J_f, E_f, S_f = compute_elast_tensors(dim, u_f, lambda_s, mu_s)
704+
F_s, J_s, E_s, S_s = compute_elast_tensors(dim, u_s, lambda_s, mu_s)
705+
return (inner(rho_s * J_s * v_s_dot, dv_s) +
706+
inner(dot(F_s, S_s), grad(dv_s)) -
707+
inner(rho_s * J_s * as_vector([0., - g_s]), dv_s) +
708+
inner(J_s * (u_s_dot - v_s), du_s)) * dx_s + \
709+
inner(dot(- p('|') * Identity(dim) + rho_f * nu_f * 2 * sym(dot(grad(v_f('|')), inv(F_f))), dot(J_f * transpose(inv(F_f)), n_f)), dv_s('|')) * ds_s(label_interface)
710+
residual_f = theta_p * _fluid(v_f, u_f, p) + \
711+
theta_m * _fluid(v_f_0, u_f_0, p_0)
712+
residual_s = theta_p * _struct(v_f, u_f, p, v_s, u_s) + \
713+
theta_m * _struct(v_f_0, u_f_0, p_0, v_s_0, u_s_0)
714+
"""
715+
F_f, J_f, E_f, S_f = compute_elast_tensors(dim, u_f, lambda_s, mu_s)
716+
F_s, J_s, E_s, S_s = compute_elast_tensors(dim, u_s, lambda_s, mu_s)
717+
v_f_mid = v_f
718+
v_s_mid = v_s
719+
u_f_mid = u_f
720+
p_mid = p
721+
residual_f = (
722+
inner(rho_f * J_f * (v_f - v_f_0) / dt, dv_f) +
723+
inner(rho_f * J_f * dot(dot(grad(v_f_mid), inv(F_f)), v_f_mid - (u_f - u_f_0) / dt), dv_f) +
724+
inner(rho_f * J_f * nu_f * 2 * sym(dot(grad(v_f_mid), inv(F_f))), dot(grad(dv_f), inv(F_f))) -
725+
J_f * inner(p_mid, tr(dot(grad(dv_f), inv(F_f)))) +
726+
J_f * inner(tr(dot(grad(v_f_mid), inv(F_f))), dp) +
727+
J_f * inner(dot(grad(u_f_mid), inv(F_f)), dot(grad(du_f), inv(F_f)))
728+
) * dx_f
729+
residual_s = (
730+
inner(rho_s * J_s * (v_s - v_s_0) / dt, dv_s) +
731+
inner(dot(F_s, S_s), grad(dv_s)) -
732+
inner(rho_s * J_s * as_vector([0., - g_s]), dv_s) +
733+
inner(J_s * ((u_s - u_s_0) / dt - v_s_mid), du_s)
734+
) * dx_s + \
735+
inner(dot(- p('|') * Identity(dim) + rho_f * nu_f * 2 * sym(dot(grad(v_f('|')), inv(F_f))), dot(J_f * transpose(inv(F_f)), n_f)), dv_s('|')) * ds_s(label_interface)
736+
#inner(dot(- p_mid * Identity(dim) + rho_f * nu_f * 2 * sym(dot(grad(v_f_mid), inv(F_f))), dot(J_f * transpose(inv(F_f)), n_f)), dv_s('|')) * ds_s(label_interface)
737+
"""
681738
residual = residual_f + residual_s
682-
#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.)
683-
v_f_left = 1.5 * Ubar * y_f * (H - y_f) / ((H / 2) ** 2) * conditional(t < 2.0 + dt / 10., (1 - cos(pi / 2 * (t - dt / 2))) / 2., 1.)
684-
bc_v_f_inflow = DirichletBC(V.sub(0), as_vector([v_f_left, 0.]), (label_left, ))
739+
def v_f_left(t_):
740+
return 1.5 * Ubar * y_f * (H - y_f) / ((H / 2) ** 2) * conditional(t_ < 2.0 + dt / 10., (1 - cos(pi / 2 * t_)) / 2., 1.)
741+
bc_v_f_inflow = DirichletBC(V.sub(0), as_vector([theta_p * v_f_left(t - dt + theta_p * dt) +
742+
theta_m * v_f_left(t - dt + theta_m * dt), 0.]), (label_left, ))
685743
bc_v_f_zero = DirichletBC(V.sub(0), Constant((0, 0)), (label_bottom, label_top, label_circle))
686744
bbc_v_f_noslip = DirichletBC(V.sub(0), Constant((0, 0)), ((label_circle, label_interface), ))
687745
bc_v_f_noslip = EquationBC(inner(v_f - v_s, dv_f) * ds_f(label_interface) == 0, solution, label_interface, bcs=[bbc_v_f_noslip], V=V.sub(0))
@@ -746,8 +804,8 @@ def compute_elast_tensors(dim, u, lambda_s, mu_s):
746804
sample_t = np.arange(0.0, T, dt_plot) + dt_plot
747805
sample_FD = np.empty_like(sample_t)
748806
sample_FL = np.empty_like(sample_t)
749-
print("num cells = ", mesh.comm.allreduce(mesh.cell_set.size))
750-
print("num DoFs = ", V.dim())
807+
print("num cells = ", mesh.comm.allreduce(mesh.cell_set.size), flush=True)
808+
print("num DoFs = ", V.dim(), flush=True)
751809
if mesh.comm.rank == 0:
752810
with open("time_series_FD.dat", 'w') as outfile:
753811
outfile.write("t val" + "\n")
@@ -770,7 +828,7 @@ def compute_elast_tensors(dim, u, lambda_s, mu_s):
770828
print(f"FD = {FD}")
771829
print(f"FL = {FL}")
772830
print(f"uA = {u_A}")
773-
if itimestep % (ntimesteps // nsample) == 0:
831+
if (itimestep + 1) % (ntimesteps // nsample) == 0:
774832
with open("time_series_FD.dat", 'a') as outfile:
775833
outfile.write(f"{float(t)} {FD}" + "\n")
776834
with open("time_series_FL.dat", 'a') as outfile:

‎mesh_f.pdf‎

-53.5 KB
Binary file not shown.

‎mesh_orig.pdf‎

-53.5 KB
Binary file not shown.

‎mesh_s.pdf‎

-629 Bytes
Binary file not shown.

‎time_series.pdf‎

4.47 KB
Binary file not shown.

0 commit comments

Comments
 (0)