diff --git a/benchmark.py b/benchmark.py index d9b93b0265..eba7b61a88 100644 --- a/benchmark.py +++ b/benchmark.py @@ -173,7 +173,7 @@ def _elevate_degree(mesh, degree): return Mesh(f) -use_netgen = False +use_netgen = True quadrilateral = True T = 20 # 10.0 # 12.0 @@ -183,9 +183,9 @@ def _elevate_degree(mesh, degree): ntimesteps = int(T / dt_float) t = Constant(0.0) dim = 2 -degree = 3 # 2 - 4 +degree = 2 # 2 - 4 if use_netgen: - nref = 2 # # 2 - 5 tested for CSM 1 and 2 + nref = 1 # # 2 - 5 tested for CSM 1 and 2 mesh = make_mesh_netgen(0.1 / 2 ** nref) mesh = _elevate_degree(mesh, degree) mesh_f = Submesh(mesh, dmcommon.CELL_SETS_LABEL, label_fluid, mesh.topological_dimension()) @@ -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 = 1 + nref = 0 mesh = make_mesh(quadrilateral) if nref > 0: mesh = MeshHierarchy(mesh, nref)[-1] @@ -658,6 +658,8 @@ def compute_elast_tensors(dim, u, lambda_s, mu_s): S = lambda_s * tr(E) * Identity(dim) + 2.0 * mu_s * E return F, J, E, S if True: # implicit midpoint + theta_p = Constant(1. / 2.) + theta_m = Constant(1. / 2.) 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 @@ -681,38 +683,65 @@ def compute_elast_tensors(dim, u, lambda_s, mu_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.) + #theta_p = Constant(1. / 2. + 100 * float(dt)) + #theta_m = Constant(1. / 2. - 100 * float(dt)) + theta_p = Constant(1.) + theta_m = Constant(0.) 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) + + def _fluid(v_f, u_f, p): + F_f, J_f, E_f, S_f = compute_elast_tensors(dim, u_f, lambda_s, mu_s) + return (inner(rho_f * J_f * v_f_dot, 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): + def _struct(v_f, u_f, p, v_s, u_s): + 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) 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) + inner(J_s * (u_s_dot - v_s), du_s)) * dx_s + \ + 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_f = theta_p * _fluid(v_f, u_f, p) + \ + # theta_m * _fluid(v_f_0, u_f_0, p_0) + #residual_s = theta_p * _struct(v_f, u_f, p, v_s, u_s) + \ + # theta_m * _struct(v_f_0, u_f_0, p_0, v_s_0, u_s_0) + 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) + v_f_mid = v_f + v_s_mid = v_s + u_f_mid = u_f + p_mid = p + """ + 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('|') * 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) + #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) + """ + residual_f = _fluid(v_f, u_f, p) + residual_s = _struct(v_f, u_f, p, v_s, u_s) 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.) - bc_v_f_inflow = DirichletBC(V.sub(0), as_vector([v_f_left, 0.]), (label_left, )) + def v_f_left(t_): + 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.) + bc_v_f_inflow = DirichletBC(V.sub(0), as_vector([theta_p * v_f_left(t - dt + theta_p * dt) + + theta_m * v_f_left(t - dt + theta_m * dt), 0.]), (label_left, )) bc_v_f_zero = DirichletBC(V.sub(0), Constant((0, 0)), (label_bottom, label_top, label_circle)) bbc_v_f_noslip = DirichletBC(V.sub(0), Constant((0, 0)), ((label_circle, label_interface), )) 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)) diff --git a/time_series.pdf b/time_series.pdf index e03f78500a..7eba2d278d 100644 Binary files a/time_series.pdf and b/time_series.pdf differ