From dfe79ccf6d77eb6f0a3b1a81e7b862db4a6da231 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Fri, 25 Oct 2024 10:04:37 +0200 Subject: [PATCH 1/4] Get original files from tutorials/cylinder-flap-3d/OpenFOAM-FEniCS/Solid/ at 0d51b33aa276d5d65a67281187cdea6b25e68cb4. --- turek-hron-fsi3/solid-fenics/cyl-flap.py | 253 ++++++++++++++++++ .../precice-adapter-config-fsi-s.json | 9 + 2 files changed, 262 insertions(+) create mode 100644 turek-hron-fsi3/solid-fenics/cyl-flap.py create mode 100644 turek-hron-fsi3/solid-fenics/precice-adapter-config-fsi-s.json diff --git a/turek-hron-fsi3/solid-fenics/cyl-flap.py b/turek-hron-fsi3/solid-fenics/cyl-flap.py new file mode 100644 index 000000000..09893e4a5 --- /dev/null +++ b/turek-hron-fsi3/solid-fenics/cyl-flap.py @@ -0,0 +1,253 @@ +# Import required libs +from fenics import Constant, Function, AutoSubDomain, RectangleMesh, VectorFunctionSpace, interpolate, \ + TrialFunction, TestFunction, Point, Expression, DirichletBC, nabla_grad, \ + Identity, inner, dx, ds, sym, grad, lhs, rhs, dot, File, solve, PointSource, assemble_system +import dolfin + +from ufl import nabla_div +import numpy as np +import matplotlib.pyplot as plt +from fenicsprecice import Adapter +from enum import Enum + +# Beam geometry +dim = 2 # number of dimensions +L = 0.35 # length +H = 0.02 # height + +y_bottom = 0.2 - 0.5 * H # y coordinate of bottom surface of beam +y_top = y_bottom + H # y coordinate of top surface of beam +x_left = 0.25 # x coordinate of left surface of beam +x_right = x_left + L # x coordinate of right surface of beam + + +# define the two inside functions to determine the boundaries: clamped Dirichlet and coupling Neumann Boundary +def left_boundary(x, on_boundary): + """ + inside-function for the clamped Dirichlet Boundary. + + Apply Dirichlet boundary on left part of the beam. + """ + return on_boundary and abs(x[0] - x_left) < tol # left boundary of beam + + +def remaining_boundary(x, on_boundary): + """ + inside-function for the coupling Neumann Boundary. + + Apply Neumann boundary on top, bottom and right (=remaining) part of the beam. + """ + return on_boundary and (abs(x[1] - y_top) < tol or # top boundary of beam + abs(x[1] - y_bottom) < tol or # bottom boundary of beam + abs(x[0] - x_right) < tol) # right boundary of beam + + +# Numerical properties +tol = 1E-14 + +# Beam material properties +rho = 1000 # density +E = 5600000.0 # Young's modulus +nu = 0.4 # Poisson's ratio +lambda_ = Constant(E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))) # first Lame constant +mu = Constant(E / (2.0 * (1.0 + nu))) # second Lame constant + +# create Mesh +n_x_Direction = 20 # DoFs in x direction +n_y_Direction = 4 # DoFs in y direction +mesh = RectangleMesh(Point(x_left, y_bottom), Point(x_right, y_top), n_x_Direction, n_y_Direction) + +# create Function Space +V = VectorFunctionSpace(mesh, 'P', 2) + +# Trial and Test Functions +du = TrialFunction(V) +v = TestFunction(V) + +# displacement fields +u_np1 = Function(V) +saved_u_old = Function(V) + +# function known from previous timestep +u_n = Function(V) +v_n = Function(V) +a_n = Function(V) + +# initial value for force and displacement field +f_N_function = interpolate(Expression(("1", "0"), degree=1), V) +u_function = interpolate(Expression(("0", "0"), degree=1), V) + +# define coupling boundary +coupling_boundary = AutoSubDomain(remaining_boundary) + +# create Adapter +precice = Adapter(adapter_config_filename="precice-adapter-config-fsi-s.json") + +# create subdomains used by the adapter +clamped_boundary_domain = AutoSubDomain(left_boundary) +force_boundary = AutoSubDomain(remaining_boundary) + +# Initialize the coupling interface +# Function space V is passed twice as both read and write functions are defined using the same space +precice_dt = precice.initialize(coupling_boundary, read_function_space=V, write_object=V, + fixed_boundary=clamped_boundary_domain) + +fenics_dt = precice_dt # if fenics_dt == precice_dt, no subcycling is applied +# fenics_dt = 0.02 # if fenics_dt < precice_dt, subcycling is applied +dt = Constant(np.min([precice_dt, fenics_dt])) + +# generalized alpha method (time stepping) parameters +alpha_m = Constant(0.2) +alpha_f = Constant(0.4) +gamma = Constant(0.5 + alpha_f - alpha_m) +beta = Constant((gamma + 0.5) ** 2 * 0.25) + +# clamp (u == 0) the beam at the left +bc = DirichletBC(V, Constant((0, 0)), left_boundary) + + +# Define strain +def epsilon(u): + return 0.5 * (nabla_grad(u) + nabla_grad(u).T) + + +# Define Stress tensor +def sigma(u): + return lambda_ * nabla_div(u) * Identity(dim) + 2 * mu * epsilon(u) + + +# Define Mass form +def m(u, v): + return rho * inner(u, v) * dx + + +# Elastic stiffness form +def k(u, v): + return inner(sigma(u), sym(grad(v))) * dx + + +def update_acceleration(u, u_old, v_old, a_old, ufl=True): + if ufl: + dt_ = dt + beta_ = beta + else: + dt_ = float(dt) + beta_ = float(beta) + + return (u - u_old - dt_ * v_old) / beta / dt_ ** 2 - .5 * (1 - 2 * beta_) / beta_ * a_old + + +def update_velocity(a, u_old, v_old, a_old, ufl=True): + if ufl: + dt_ = dt + gamma_ = gamma + else: + dt_ = float(dt) + gamma_ = float(gamma) + + return v_old + dt_ * ((1 - gamma_) * a_old + gamma_ * a) + + +def update_fields(u, u_old, v_old, a_old): + """Update all fields at the end of a timestep.""" + + u_vec, u0_vec = u.vector(), u_old.vector() + v0_vec, a0_vec = v_old.vector(), a_old.vector() + + # call update functions + a_vec = update_acceleration(u_vec, u0_vec, v0_vec, a0_vec, ufl=False) + v_vec = update_velocity(a_vec, u0_vec, v0_vec, a0_vec, ufl=False) + + # assign u->u_old + v_old.vector()[:], a_old.vector()[:] = v_vec, a_vec + u_old.vector()[:] = u.vector() + + +def avg(x_old, x_new, alpha): + return alpha * x_old + (1 - alpha) * x_new + + +# residual +a_np1 = update_acceleration(du, u_n, v_n, a_n, ufl=True) +v_np1 = update_velocity(a_np1, u_n, v_n, a_n, ufl=True) + +res = m(avg(a_n, a_np1, alpha_m), v) + k(avg(u_n, du, alpha_f), v) + +a_form = lhs(res) +L_form = rhs(res) + +# Prepare for time-stepping +t = 0.0 +n = 0 +time = [] +u_tip = [] +time.append(0.0) +u_tip.append(0.0) +E_ext = 0 + +displacement_out = File("Solid/FSI-S/u_fsi.pvd") + +u_n.rename("Displacement", "") +u_np1.rename("Displacement", "") +displacement_out << u_n + +# time loop for coupling +while precice.is_coupling_ongoing(): + + if precice.is_action_required(precice.action_write_iteration_checkpoint()): # write checkpoint + precice.store_checkpoint(u_n, t, n) + + # read data from preCICE and get a new coupling expression + read_data = precice.read_data() + + # Update the point sources on the coupling boundary with the new read data + Forces_x, Forces_y = precice.get_point_sources(read_data) + + A, b = assemble_system(a_form, L_form, bc) + + b_forces = b.copy() # b is the same for every iteration, only forces change + + for ps in Forces_x: + ps.apply(b_forces) + for ps in Forces_y: + ps.apply(b_forces) + + assert (b is not b_forces) + solve(A, u_np1.vector(), b_forces) + + dt = Constant(np.min([precice_dt, fenics_dt])) + + # Write new displacements to preCICE + precice.write_data(u_np1) + + # Call to advance coupling, also returns the optimum time step value + precice_dt = precice.advance(dt(0)) + + # Either revert to old step if timestep has not converged or move to next timestep + if precice.is_action_required(precice.action_read_iteration_checkpoint()): # roll back to checkpoint + u_cp, t_cp, n_cp = precice.retrieve_checkpoint() + u_n.assign(u_cp) + t = t_cp + n = n_cp + else: + u_n.assign(u_np1) + t += float(dt) + n += 1 + + if precice.is_time_window_complete(): + update_fields(u_np1, saved_u_old, v_n, a_n) + if n % 20 == 0: + displacement_out << (u_n, t) + + u_tip.append(u_n(0.6, 0.2)[1]) + time.append(t) + +# Plot tip displacement evolution +displacement_out << u_n +plt.figure() +plt.plot(time, u_tip) +plt.xlabel("Time") +plt.ylabel("Tip displacement") +plt.show() + +precice.finalize() diff --git a/turek-hron-fsi3/solid-fenics/precice-adapter-config-fsi-s.json b/turek-hron-fsi3/solid-fenics/precice-adapter-config-fsi-s.json new file mode 100644 index 000000000..575bcf3bb --- /dev/null +++ b/turek-hron-fsi3/solid-fenics/precice-adapter-config-fsi-s.json @@ -0,0 +1,9 @@ +{ + "participant_name": "fenics", + "config_file_name": "../precice-config.xml", + "interface": { + "coupling_mesh_name": "Solid", + "write_data_name": "Displacements0", + "read_data_name": "Forces0" + } +} From 855a6922906883f5261ea9ec102d2230bfd95916 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Fri, 25 Oct 2024 14:04:35 +0200 Subject: [PATCH 2/4] Working configuration. --- .../fluid-openfoam/system/decomposeParDict | 4 +- .../fluid-openfoam/system/preciceDict | 3 +- turek-hron-fsi3/plot-displacement.sh | 2 +- turek-hron-fsi3/precice-config.xml | 25 +++--- turek-hron-fsi3/solid-fenics/clean.sh | 6 ++ .../precice-adapter-config-fsi-s.json | 10 +-- turek-hron-fsi3/solid-fenics/requirements.txt | 11 +++ turek-hron-fsi3/solid-fenics/run.sh | 13 +++ .../solid-fenics/{cyl-flap.py => solid.py} | 83 +++++++++---------- 9 files changed, 91 insertions(+), 66 deletions(-) create mode 100755 turek-hron-fsi3/solid-fenics/clean.sh create mode 100644 turek-hron-fsi3/solid-fenics/requirements.txt create mode 100755 turek-hron-fsi3/solid-fenics/run.sh rename turek-hron-fsi3/solid-fenics/{cyl-flap.py => solid.py} (76%) diff --git a/turek-hron-fsi3/fluid-openfoam/system/decomposeParDict b/turek-hron-fsi3/fluid-openfoam/system/decomposeParDict index a9e847699..85cf083c2 100644 --- a/turek-hron-fsi3/fluid-openfoam/system/decomposeParDict +++ b/turek-hron-fsi3/fluid-openfoam/system/decomposeParDict @@ -6,12 +6,12 @@ FoamFile object decomposeParDict; } -numberOfSubdomains 25; +numberOfSubdomains 3; method hierarchical; hierarchicalCoeffs { - n (5 5 1); + n (3 1 1); delta 0.001; order xyz; } diff --git a/turek-hron-fsi3/fluid-openfoam/system/preciceDict b/turek-hron-fsi3/fluid-openfoam/system/preciceDict index 872566d38..2a2b63d9a 100644 --- a/turek-hron-fsi3/fluid-openfoam/system/preciceDict +++ b/turek-hron-fsi3/fluid-openfoam/system/preciceDict @@ -26,10 +26,9 @@ interfaces writeData ( - Stress + Force ); }; - Interface2 { mesh Fluid-Mesh-Nodes; diff --git a/turek-hron-fsi3/plot-displacement.sh b/turek-hron-fsi3/plot-displacement.sh index 141f54ec5..b44e3b759 100755 --- a/turek-hron-fsi3/plot-displacement.sh +++ b/turek-hron-fsi3/plot-displacement.sh @@ -5,5 +5,5 @@ gnuplot -p << EOF set xlabel 'time [s]' set ylabel 'y-displacement [m]' set linestyle 1 lt 2 lc 1 # red-dashed - plot "solid-dealii/precice-Solid-watchpoint-Flap-Tip.log" using 1:5 with lines notitle + plot "solid-fenics/precice-Solid-watchpoint-Flap-Tip.log" using 1:5 with lines notitle EOF diff --git a/turek-hron-fsi3/precice-config.xml b/turek-hron-fsi3/precice-config.xml index ae004917b..dfbf7f1ac 100644 --- a/turek-hron-fsi3/precice-config.xml +++ b/turek-hron-fsi3/precice-config.xml @@ -7,11 +7,11 @@ enabled="true" /> - + - + @@ -20,29 +20,28 @@ - + - + - + + + + - - + - - - @@ -51,14 +50,14 @@ - + - + - + diff --git a/turek-hron-fsi3/solid-fenics/clean.sh b/turek-hron-fsi3/solid-fenics/clean.sh new file mode 100755 index 000000000..d3f3318b4 --- /dev/null +++ b/turek-hron-fsi3/solid-fenics/clean.sh @@ -0,0 +1,6 @@ +#!/usr/bin/env sh +set -e -u + +. ../../tools/cleaning-tools.sh + +clean_fenics . diff --git a/turek-hron-fsi3/solid-fenics/precice-adapter-config-fsi-s.json b/turek-hron-fsi3/solid-fenics/precice-adapter-config-fsi-s.json index 575bcf3bb..91a042b17 100644 --- a/turek-hron-fsi3/solid-fenics/precice-adapter-config-fsi-s.json +++ b/turek-hron-fsi3/solid-fenics/precice-adapter-config-fsi-s.json @@ -1,9 +1,9 @@ { - "participant_name": "fenics", - "config_file_name": "../precice-config.xml", + "participant_name": "Solid", + "config_file_name": "../precice-config.xml", "interface": { - "coupling_mesh_name": "Solid", - "write_data_name": "Displacements0", - "read_data_name": "Forces0" + "coupling_mesh_name": "Solid-Mesh", + "write_data_name": "Displacement", + "read_data_name": "Force" } } diff --git a/turek-hron-fsi3/solid-fenics/requirements.txt b/turek-hron-fsi3/solid-fenics/requirements.txt new file mode 100644 index 000000000..fe29958de --- /dev/null +++ b/turek-hron-fsi3/solid-fenics/requirements.txt @@ -0,0 +1,11 @@ +fenicsprecice~=2.0 +numpy >1, <2 +matplotlib + +# Assuming FEniCS from ppa:fenics-packages/fenics was installed https://fenicsproject.org/download/archive/ +# Use --system-site-packages in venv +fenics-dijitso==2019.2.0.dev0 +fenics-dolfin==2019.2.0.13.dev0 +fenics-ffc==2019.2.0.dev0 +fenics-fiat==2019.2.0.dev0 +fenics-ufl-legacy==2022.3.0 \ No newline at end of file diff --git a/turek-hron-fsi3/solid-fenics/run.sh b/turek-hron-fsi3/solid-fenics/run.sh new file mode 100755 index 000000000..867c1db13 --- /dev/null +++ b/turek-hron-fsi3/solid-fenics/run.sh @@ -0,0 +1,13 @@ +#!/usr/bin/env bash +set -e -u + +. ../../tools/log.sh +exec > >(tee --append "$LOGFILE") 2>&1 + +python3 -m venv --system-site-packages .venv +. .venv/bin/activate +pip install -r requirements.txt + +python3 solid.py + +close_log diff --git a/turek-hron-fsi3/solid-fenics/cyl-flap.py b/turek-hron-fsi3/solid-fenics/solid.py similarity index 76% rename from turek-hron-fsi3/solid-fenics/cyl-flap.py rename to turek-hron-fsi3/solid-fenics/solid.py index 09893e4a5..2cafb9344 100644 --- a/turek-hron-fsi3/solid-fenics/cyl-flap.py +++ b/turek-hron-fsi3/solid-fenics/solid.py @@ -1,12 +1,9 @@ # Import required libs from fenics import Constant, Function, AutoSubDomain, RectangleMesh, VectorFunctionSpace, interpolate, \ TrialFunction, TestFunction, Point, Expression, DirichletBC, nabla_grad, \ - Identity, inner, dx, ds, sym, grad, lhs, rhs, dot, File, solve, PointSource, assemble_system -import dolfin + Identity, inner, dx, sym, grad, lhs, rhs, File, solve, assemble_system, project, div -from ufl import nabla_div import numpy as np -import matplotlib.pyplot as plt from fenicsprecice import Adapter from enum import Enum @@ -89,16 +86,26 @@ def remaining_boundary(x, on_boundary): # Initialize the coupling interface # Function space V is passed twice as both read and write functions are defined using the same space -precice_dt = precice.initialize(coupling_boundary, read_function_space=V, write_object=V, - fixed_boundary=clamped_boundary_domain) +precice.initialize(coupling_boundary, read_function_space=V, write_object=V, fixed_boundary=clamped_boundary_domain) -fenics_dt = precice_dt # if fenics_dt == precice_dt, no subcycling is applied -# fenics_dt = 0.02 # if fenics_dt < precice_dt, subcycling is applied +precice_dt = precice.get_max_time_step_size() +fenics_dt = precice_dt dt = Constant(np.min([precice_dt, fenics_dt])) -# generalized alpha method (time stepping) parameters -alpha_m = Constant(0.2) -alpha_f = Constant(0.4) +# alpha method parameters +# alpha_m = Constant(0.2) +# alpha_f = Constant(0.4) +alpha_m = Constant(0) +alpha_f = Constant(0) + +""" +Check requirements for alpha_m and alpha_f from + Chung, J., and Hulbert, G. M. (June 1, 1993). "A Time Integration Algorithm for Structural Dynamics With Improved Numerical Dissipation: + The Generalized-α Method." ASME. J. Appl. Mech. June 1993; 60(2): 371–375. https://doi.org/10.1115/1.2900803 +""" +assert (float(alpha_m) <= float(alpha_f)) +assert (float(alpha_f) <= 0.5) + gamma = Constant(0.5 + alpha_f - alpha_m) beta = Constant((gamma + 0.5) ** 2 * 0.25) @@ -113,7 +120,7 @@ def epsilon(u): # Define Stress tensor def sigma(u): - return lambda_ * nabla_div(u) * Identity(dim) + 2 * mu * epsilon(u) + return lambda_ * div(u) * Identity(dim) + 2 * mu * epsilon(u) # Define Mass form @@ -151,16 +158,14 @@ def update_velocity(a, u_old, v_old, a_old, ufl=True): def update_fields(u, u_old, v_old, a_old): """Update all fields at the end of a timestep.""" - u_vec, u0_vec = u.vector(), u_old.vector() - v0_vec, a0_vec = v_old.vector(), a_old.vector() - # call update functions - a_vec = update_acceleration(u_vec, u0_vec, v0_vec, a0_vec, ufl=False) - v_vec = update_velocity(a_vec, u0_vec, v0_vec, a0_vec, ufl=False) + a_new = update_acceleration(u, u_old, v_old, a_old) + v_new = update_velocity(a_new, u_old, v_old, a_old) # assign u->u_old - v_old.vector()[:], a_old.vector()[:] = v_vec, a_vec - u_old.vector()[:] = u.vector() + a_old.assign(project(a_new, V)) + v_old.assign(project(v_new, V)) + u_old.assign(u) def avg(x_old, x_new, alpha): @@ -179,13 +184,9 @@ def avg(x_old, x_new, alpha): # Prepare for time-stepping t = 0.0 n = 0 -time = [] -u_tip = [] -time.append(0.0) -u_tip.append(0.0) E_ext = 0 -displacement_out = File("Solid/FSI-S/u_fsi.pvd") +displacement_out = File("output/u_fsi.pvd") u_n.rename("Displacement", "") u_np1.rename("Displacement", "") @@ -194,11 +195,16 @@ def avg(x_old, x_new, alpha): # time loop for coupling while precice.is_coupling_ongoing(): - if precice.is_action_required(precice.action_write_iteration_checkpoint()): # write checkpoint - precice.store_checkpoint(u_n, t, n) + if precice.requires_writing_checkpoint(): # write checkpoint + precice.store_checkpoint((u_n, v_n, a_n), t, n) + + precice_dt = precice.get_max_time_step_size() + dt = Constant(np.min([precice_dt, fenics_dt])) # read data from preCICE and get a new coupling expression - read_data = precice.read_data() + # sample force F at $F(t_{n+1-\alpha_f})$ (see generalized alpha paper) + read_data = precice.read_data((1 - float(alpha_f)) * dt) + # read_data = precice.read_data(dt) # Update the point sources on the coupling boundary with the new read data Forces_x, Forces_y = precice.get_point_sources(read_data) @@ -215,39 +221,30 @@ def avg(x_old, x_new, alpha): assert (b is not b_forces) solve(A, u_np1.vector(), b_forces) - dt = Constant(np.min([precice_dt, fenics_dt])) - # Write new displacements to preCICE precice.write_data(u_np1) # Call to advance coupling, also returns the optimum time step value - precice_dt = precice.advance(dt(0)) + precice.advance(float(dt)) # Either revert to old step if timestep has not converged or move to next timestep - if precice.is_action_required(precice.action_read_iteration_checkpoint()): # roll back to checkpoint - u_cp, t_cp, n_cp = precice.retrieve_checkpoint() + if precice.requires_reading_checkpoint(): # roll back to checkpoint + uva_cp, t_cp, n_cp = precice.retrieve_checkpoint() + u_cp, v_cp, a_cp = uva_cp u_n.assign(u_cp) + v_n.assign(v_cp) + a_n.assign(a_cp) t = t_cp n = n_cp else: - u_n.assign(u_np1) + update_fields(u_np1, u_n, v_n, a_n) t += float(dt) n += 1 if precice.is_time_window_complete(): - update_fields(u_np1, saved_u_old, v_n, a_n) if n % 20 == 0: displacement_out << (u_n, t) - u_tip.append(u_n(0.6, 0.2)[1]) - time.append(t) - -# Plot tip displacement evolution displacement_out << u_n -plt.figure() -plt.plot(time, u_tip) -plt.xlabel("Time") -plt.ylabel("Tip displacement") -plt.show() precice.finalize() From 597d7a8e8ab8aff0705a83a7e9dabb244aae025e Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Fri, 25 Oct 2024 17:06:21 +0200 Subject: [PATCH 3/4] Apply changes from https://github.com/precice/tutorials/pull/578. --- turek-hron-fsi3/precice-config.xml | 2 +- turek-hron-fsi3/solid-fenics/solid.py | 30 ++++++++++++++------------- 2 files changed, 17 insertions(+), 15 deletions(-) diff --git a/turek-hron-fsi3/precice-config.xml b/turek-hron-fsi3/precice-config.xml index dfbf7f1ac..26f3f77a5 100644 --- a/turek-hron-fsi3/precice-config.xml +++ b/turek-hron-fsi3/precice-config.xml @@ -48,7 +48,7 @@ - + diff --git a/turek-hron-fsi3/solid-fenics/solid.py b/turek-hron-fsi3/solid-fenics/solid.py index 2cafb9344..eeae08ceb 100644 --- a/turek-hron-fsi3/solid-fenics/solid.py +++ b/turek-hron-fsi3/solid-fenics/solid.py @@ -1,7 +1,7 @@ # Import required libs from fenics import Constant, Function, AutoSubDomain, RectangleMesh, VectorFunctionSpace, interpolate, \ TrialFunction, TestFunction, Point, Expression, DirichletBC, nabla_grad, \ - Identity, inner, dx, sym, grad, lhs, rhs, File, solve, assemble_system, project, div + Identity, inner, dx, sym, grad, lhs, rhs, File, solve, assemble_system, div import numpy as np from fenicsprecice import Adapter @@ -93,10 +93,10 @@ def remaining_boundary(x, on_boundary): dt = Constant(np.min([precice_dt, fenics_dt])) # alpha method parameters -# alpha_m = Constant(0.2) -# alpha_f = Constant(0.4) -alpha_m = Constant(0) -alpha_f = Constant(0) +alpha_m = Constant(0.2) +alpha_f = Constant(0.4) +# alpha_m = Constant(0) +# alpha_f = Constant(0) """ Check requirements for alpha_m and alpha_f from @@ -133,7 +133,7 @@ def k(u, v): return inner(sigma(u), sym(grad(v))) * dx -def update_acceleration(u, u_old, v_old, a_old, ufl=True): +def update_a(u, u_old, v_old, a_old, ufl=True): if ufl: dt_ = dt beta_ = beta @@ -144,7 +144,7 @@ def update_acceleration(u, u_old, v_old, a_old, ufl=True): return (u - u_old - dt_ * v_old) / beta / dt_ ** 2 - .5 * (1 - 2 * beta_) / beta_ * a_old -def update_velocity(a, u_old, v_old, a_old, ufl=True): +def update_v(a, u_old, v_old, a_old, ufl=True): if ufl: dt_ = dt gamma_ = gamma @@ -157,15 +157,16 @@ def update_velocity(a, u_old, v_old, a_old, ufl=True): def update_fields(u, u_old, v_old, a_old): """Update all fields at the end of a timestep.""" + u_vec, u0_vec = u.vector(), u_old.vector() + v0_vec, a0_vec = v_old.vector(), a_old.vector() # call update functions - a_new = update_acceleration(u, u_old, v_old, a_old) - v_new = update_velocity(a_new, u_old, v_old, a_old) + a_vec = update_a(u_vec, u0_vec, v0_vec, a0_vec, ufl=False) + v_vec = update_v(a_vec, u0_vec, v0_vec, a0_vec, ufl=False) # assign u->u_old - a_old.assign(project(a_new, V)) - v_old.assign(project(v_new, V)) - u_old.assign(u) + v_old.vector()[:], a_old.vector()[:] = v_vec, a_vec + u_old.vector()[:] = u.vector() def avg(x_old, x_new, alpha): @@ -173,8 +174,8 @@ def avg(x_old, x_new, alpha): # residual -a_np1 = update_acceleration(du, u_n, v_n, a_n, ufl=True) -v_np1 = update_velocity(a_np1, u_n, v_n, a_n, ufl=True) +a_np1 = update_a(du, u_n, v_n, a_n, ufl=True) +v_np1 = update_v(a_np1, u_n, v_n, a_n, ufl=True) res = m(avg(a_n, a_np1, alpha_m), v) + k(avg(u_n, du, alpha_f), v) @@ -238,6 +239,7 @@ def avg(x_old, x_new, alpha): n = n_cp else: update_fields(u_np1, u_n, v_n, a_n) + u_n.assign(u_np1) t += float(dt) n += 1 From 529840436934033520e73f3120bba550f718c910 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Fri, 25 Oct 2024 19:16:28 +0200 Subject: [PATCH 4/4] Require latest fenicsadapter release. --- turek-hron-fsi3/solid-fenics/requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/turek-hron-fsi3/solid-fenics/requirements.txt b/turek-hron-fsi3/solid-fenics/requirements.txt index fe29958de..e7a847a98 100644 --- a/turek-hron-fsi3/solid-fenics/requirements.txt +++ b/turek-hron-fsi3/solid-fenics/requirements.txt @@ -1,4 +1,4 @@ -fenicsprecice~=2.0 +fenicsprecice==2.2.0rc1 numpy >1, <2 matplotlib