From 7e3c2480499ee7d685b303abf96e3fed2fdb796c Mon Sep 17 00:00:00 2001 From: ksagiyam Date: Wed, 15 May 2024 22:54:53 +0100 Subject: [PATCH] benchmark --- benchmark.py | 31 ++++++++++++++++++++----------- 1 file changed, 20 insertions(+), 11 deletions(-) diff --git a/benchmark.py b/benchmark.py index 9b09932da5..8099046dc0 100644 --- a/benchmark.py +++ b/benchmark.py @@ -357,20 +357,29 @@ def _elevate_degree(mesh, degree): E_s = mu_s * 2 * (1 + nu_s) lambda_s = nu_s * E_s / (1 + nu_s) / (1 - 2 * nu_s) V_0 = VectorFunctionSpace(mesh_f, "P", degree) - V_1 = FunctionSpace(mesh_f, "P", degree - 2) - V = V_0 * V_1 + V_1 = FunctionSpace(mesh_f, "P", degree - 1) + V = V_0 * V_0 * V_1 solution = Function(V) solution_0 = Function(V) - v_f, p_f = split(solution) - u_f = Function(V_0) - v_f_0, p_f_0 = split(solution_0) - dv_f, dp_f = split(TestFunction(V)) + v_f, u_f, p_f = split(solution) + v_f_0, u_f_0, p_f_0 = split(solution_0) + dv_f, du_f, dp_f = split(TestFunction(V)) + F = Identity(dim) + grad(u_f) + J = det(F) + E = 1. / 2. * (dot(transpose(F), F) - Identity(dim)) + S = lambda_s * tr(E) * Identity(dim) + 2.0 * mu_s * E + #residual = ( + # rho_f * J * inner(v_f - v_f_0, dv_f) / dt + + # rho_f * J * inner(dot(grad(v_f), dot(inv(F), v_f - (u_f - u_f_0) / 2)), dv_f) + + # rho_f * J * nu_f * inner(dot(2 * sym(dot(grad(v_f), inv(F))), transpose(inv(F))), grad(dv_f)) - + # J * inner(p_f * transpose(inv(F)), grad(dv_f)) + + # J * inner(tr(dot(grad(v_f), inv(F))), dp_f)) * dx 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 + rho_f * J * inner(v_f - v_f_0, dv_f) / dt + + rho_f * J * inner(dot(grad(v_f), dot(inv(F), v_f - (u_f - u_f_0) / 2)), dv_f) + + rho_f * J * nu_f * inner(2 * sym(dot(grad(v_f), inv(F))), dot(grad(dv_f), inv(F))) - + J * inner(p_f, tr(dot(grad(dv_f), inv(F)))) + + J * inner(tr(dot(grad(v_f), inv(F))), dp_f) + J * inner(u_f, du_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)