Skip to content

Commit

Permalink
benchmark
Browse files Browse the repository at this point in the history
  • Loading branch information
@ksagiyam
ksagiyam committed Jun 3, 2024
1 parent efbf351 commit 561fa5d
Show file tree
Hide file tree
Showing 2 changed files with 53 additions and 23 deletions.
76 changes: 53 additions & 23 deletions benchmark.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -193,7 +193,7 @@ def _elevate_degree(mesh, degree):
mesh_s = Submesh(mesh, dmcommon.CELL_SETS_LABEL, label_struct, mesh.topological_dimension())
mesh_s = _elevate_degree(mesh_s, degree)
else:
nref = 0
nref = 1
mesh = make_mesh(quadrilateral)
if nref > 0:
mesh = MeshHierarchy(mesh, nref)[-1]
Expand Down Expand Up @@ -657,28 +657,58 @@ def compute_elast_tensors(dim, u, lambda_s, mu_s):
E = 1. / 2. * (dot(transpose(F), F) - Identity(dim))
S = lambda_s * tr(E) * Identity(dim) + 2.0 * mu_s * E
return F, J, E, S
F_f, J_f, E_f, S_f = compute_elast_tensors(dim, (u_f + u_f_0) / 2, lambda_s, mu_s)
F_s, J_s, E_s, S_s = compute_elast_tensors(dim, (u_s + u_s_0) / 2, lambda_s, mu_s)
v_f_mid = (v_f + v_f_0) / 2
v_s_mid = (v_s + v_s_0) / 2
u_f_mid = (u_f + u_f_0) / 2
p_mid = (p + p_0) / 2
residual_f = (
inner(rho_f * J_f * (v_f - v_f_0) / dt, dv_f) +
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) +
inner(rho_f * J_f * nu_f * 2 * sym(dot(grad(v_f_mid), inv(F_f))), dot(grad(dv_f), inv(F_f))) -
J_f * inner(p_mid, tr(dot(grad(dv_f), inv(F_f)))) +
J_f * inner(tr(dot(grad(v_f_mid), inv(F_f))), dp) +
J_f * inner(dot(grad(u_f_mid), inv(F_f)), dot(grad(du_f), inv(F_f)))
) * dx_f
residual_s = (
inner(rho_s * J_s * (v_s - v_s_0) / dt, dv_s) +
inner(dot(F_s, S_s), grad(dv_s)) -
inner(rho_s * J_s * as_vector([0., - g_s]), dv_s) +
inner(J_s * (u_s - u_s_0) / dt - v_s_mid, du_s)
) * dx_s + \
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)
#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)
if True: # implicit midpoint
F_f, J_f, E_f, S_f = compute_elast_tensors(dim, (u_f + u_f_0) / 2, lambda_s, mu_s)
F_s, J_s, E_s, S_s = compute_elast_tensors(dim, (u_s + u_s_0) / 2, lambda_s, mu_s)
v_f_mid = (v_f + v_f_0) / 2
v_s_mid = (v_s + v_s_0) / 2
u_f_mid = (u_f + u_f_0) / 2
p_mid = (p + p_0) / 2
residual_f = (
inner(rho_f * J_f * (v_f - v_f_0) / dt, dv_f) +
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) +
inner(rho_f * J_f * nu_f * 2 * sym(dot(grad(v_f_mid), inv(F_f))), dot(grad(dv_f), inv(F_f))) -
J_f * inner(p_mid, tr(dot(grad(dv_f), inv(F_f)))) +
J_f * inner(tr(dot(grad(v_f_mid), inv(F_f))), dp) +
J_f * inner(dot(grad(u_f_mid), inv(F_f)), dot(grad(du_f), inv(F_f)))
) * dx_f
residual_s = (
inner(rho_s * J_s * (v_s - v_s_0) / dt, dv_s) +
inner(dot(F_s, S_s), grad(dv_s)) -
inner(rho_s * J_s * as_vector([0., - g_s]), dv_s) +
inner(J_s * ((u_s - u_s_0) / dt - v_s_mid), du_s)
) * dx_s + \
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)
#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)
else: # CN
F_f, J_f, E_f, S_f = compute_elast_tensors(dim, u_f, lambda_s, mu_s)
F_s, J_s, E_s, S_s = compute_elast_tensors(dim, u_s, lambda_s, mu_s)
F_f_0, J_f_0, E_f_0, S_f_0 = compute_elast_tensors(dim, u_f_0, lambda_s, mu_s)
F_s_0, J_s_0, E_s_0, S_s_0 = compute_elast_tensors(dim, u_s_0, lambda_s, mu_s)
theta_p = Constant(1. / 2.)
theta_m = Constant(1. / 2.)
v_f_dot = (v_f - v_f_0) / dt
u_f_dot = (u_f - u_f_0) / dt
v_s_dot = (v_s - v_s_0) / dt
u_s_dot = (u_s - u_s_0) / dt
def _fluid(v_f, u_f, F_f, J_f, p):
return (inner(rho_f * J_f * vdot, dv_f) +
inner(rho_f * J_f * dot(dot(grad(v_f), inv(F_f)), v_f - u_f_dot), dv_f) +
inner(rho_f * J_f * nu_f * 2 * sym(dot(grad(v_f), inv(F_f))), dot(grad(dv_f), inv(F_f))) +
J_f * inner(p, tr(dot(grad(dv_f), inv(F_f)))) -
J_f * inner(tr(dot(grad(v_f), inv(F_f))), dp) +
J_f * inner(dot(grad(u_f), inv(F_f)), dot(grad(du_f), inv(F_f)))) * dx_f
def _struct(F_s, S_s, J_s):
return (inner(rho_s * J_s * v_s_dot, dv_s) +
inner(dot(F_s, S_s), grad(dv_s)) -
inner(rho_s * J_s * as_vector([0., - g_s]), dv_s) +
inner(J_s * (u_s - u_s_0) / dt - v_s_mid, du_s)) * dx_s
residual_f = theta_p * _fluid(v_f, u_f, F_f, J_f, p) + \
theta_m * _fluid(v_f_0, u_f_0, F_f_0, J_f_0, p_0)
residual_s = theta_p * _struct() + \
theta_m * _struct() + \
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)
#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)
residual = residual_f + residual_s
#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.)
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.)
Expand Down
Binary file modified time_series.pdf
View file
Binary file not shown.

0 comments on commit 561fa5d

Please sign in to comment.