From 8c149d95cc7eb84415d519112bc55a1ea15ecf72 Mon Sep 17 00:00:00 2001 From: ishaandesai Date: Thu, 31 Dec 2020 13:15:38 +0100 Subject: [PATCH 1/3] Making SU2- and OpenFOAM- case of perp-flap for FEniCS consistent --- FSI/flap_perp_2D/OpenFOAM-FEniCS/Allclean | 3 -- .../OpenFOAM-FEniCS/Solid/perp-flap.py | 40 +++++++++---------- .../SU2-FEniCS/Solid/perp-flap.py | 11 ++--- 3 files changed, 23 insertions(+), 31 deletions(-) diff --git a/FSI/flap_perp_2D/OpenFOAM-FEniCS/Allclean b/FSI/flap_perp_2D/OpenFOAM-FEniCS/Allclean index 335e4274b..9b055022c 100755 --- a/FSI/flap_perp_2D/OpenFOAM-FEniCS/Allclean +++ b/FSI/flap_perp_2D/OpenFOAM-FEniCS/Allclean @@ -3,9 +3,6 @@ cd ${0%/*} || exit 1 # Run from this directory echo "Cleaning..." -# Source tutorial clean functions -. $WM_PROJECT_DIR/bin/tools/CleanFunctions - # Participant 1: Fluid Participant1="Fluid" cd ${Participant1} diff --git a/FSI/flap_perp_2D/OpenFOAM-FEniCS/Solid/perp-flap.py b/FSI/flap_perp_2D/OpenFOAM-FEniCS/Solid/perp-flap.py index 4b1891d17..ee19aefca 100644 --- a/FSI/flap_perp_2D/OpenFOAM-FEniCS/Solid/perp-flap.py +++ b/FSI/flap_perp_2D/OpenFOAM-FEniCS/Solid/perp-flap.py @@ -1,9 +1,7 @@ # Import required libs from fenics import Constant, Function, AutoSubDomain, RectangleMesh, VectorFunctionSpace, interpolate, \ - TrialFunction, TestFunction, Point, Expression, DirichletBC, nabla_grad, SubDomain, \ - Identity, inner, dx, ds, sym, grad, lhs, rhs, dot, File, solve, PointSource, assemble_system, MPI, MeshFunction -import dolfin - + TrialFunction, TestFunction, Point, Expression, DirichletBC, nabla_grad, Identity, inner, dx, ds, \ + sym, grad, lhs, rhs, dot, File, solve, PointSource, assemble_system, MPI, MeshFunction from ufl import nabla_div import numpy as np import matplotlib.pyplot as plt @@ -11,22 +9,17 @@ from enum import Enum -class clampedBoundary(SubDomain): - def inside(self, x, on_boundary): - tol = 1E-14 - if on_boundary and abs(x[1]) < tol: - return True - else: - return False +# define the two kinds of boundary: clamped and coupling Neumann Boundary +def clamped_boundary(x, on_boundary): + return on_boundary and abs(x[1]) < tol + +def neumann_boundary(x, on_boundary): + """ + determines whether a node is on the coupling boundary -class neumannBoundary(SubDomain): - def inside(self, x, on_boundary): - tol = 1E-14 - if on_boundary and ((abs(x[1] - 1) < tol) or abs(abs(x[0]) - W / 2) < tol): - return True - else: - return False + """ + return on_boundary and ((abs(x[1] - 1) < tol) or abs(abs(x[0]) - W / 2) < tol) # Geometry and material properties @@ -51,6 +44,9 @@ def inside(self, x, on_boundary): # create Function Space V = VectorFunctionSpace(mesh, 'P', 2) +# BCs +tol = 1E-14 + # Trial and Test Functions du = TrialFunction(V) v = TestFunction(V) @@ -66,19 +62,21 @@ def inside(self, x, on_boundary): f_N_function = interpolate(Expression(("1", "0"), degree=1), V) u_function = interpolate(Expression(("0", "0"), degree=1), V) +coupling_boundary = AutoSubDomain(neumann_boundary) +fixed_boundary = AutoSubDomain(clamped_boundary) + precice = Adapter(adapter_config_filename="precice-adapter-config-fsi-s.json") # 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(neumannBoundary(), read_function_space=V, write_object=V, - fixed_boundary=clampedBoundary()) +precice_dt = precice.initialize(coupling_boundary, read_function_space=V, write_object=V, fixed_boundary=fixed_boundary) 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])) # clamp the beam at the bottom -bc = DirichletBC(V, Constant((0, 0)), clampedBoundary()) +bc = DirichletBC(V, Constant((0, 0)), fixed_boundary) # alpha method parameters alpha_m = Constant(0.2) diff --git a/FSI/flap_perp_2D/SU2-FEniCS/Solid/perp-flap.py b/FSI/flap_perp_2D/SU2-FEniCS/Solid/perp-flap.py index 5376588a9..a57f922a3 100644 --- a/FSI/flap_perp_2D/SU2-FEniCS/Solid/perp-flap.py +++ b/FSI/flap_perp_2D/SU2-FEniCS/Solid/perp-flap.py @@ -2,8 +2,6 @@ from fenics import Constant, Function, AutoSubDomain, RectangleMesh, VectorFunctionSpace, interpolate, \ TrialFunction, TestFunction, Point, Expression, DirichletBC, nabla_grad, project, \ 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 @@ -63,26 +61,25 @@ def neumann_boundary(x, on_boundary): # Function to calculate displacement Deltas u_delta = Function(V) +u_ref = Function(V) f_N_function = interpolate(Expression(("1", "0"), degree=1), V) u_function = interpolate(Expression(("0", "0"), degree=1), V) coupling_boundary = AutoSubDomain(neumann_boundary) +fixed_boundary = AutoSubDomain(clamped_boundary) precice = Adapter(adapter_config_filename="precice-adapter-config-fsi-s.json") -clamped_boundary_domain = AutoSubDomain(clamped_boundary) -force_boundary = AutoSubDomain(neumann_boundary) - # Initialize the coupling interface -precice_dt = precice.initialize(coupling_boundary, read_function_space=V, write_object=V, fixed_boundary=clamped_boundary_domain) +precice_dt = precice.initialize(coupling_boundary, read_function_space=V, write_object=V, fixed_boundary=fixed_boundary) 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])) # clamp the beam at the bottom -bc = DirichletBC(V, Constant((0, 0)), clamped_boundary) +bc = DirichletBC(V, Constant((0, 0)), fixed_boundary) # alpha method parameters alpha_m = Constant(0.2) From fba42f6ae4d4a20781e92b063b3acf330a78b187 Mon Sep 17 00:00:00 2001 From: ishaandesai Date: Thu, 31 Dec 2020 14:47:18 +0100 Subject: [PATCH 2/3] Reverting wrong change in Allclean script --- FSI/flap_perp_2D/OpenFOAM-FEniCS/Allclean | 3 +++ 1 file changed, 3 insertions(+) diff --git a/FSI/flap_perp_2D/OpenFOAM-FEniCS/Allclean b/FSI/flap_perp_2D/OpenFOAM-FEniCS/Allclean index 9b055022c..335e4274b 100755 --- a/FSI/flap_perp_2D/OpenFOAM-FEniCS/Allclean +++ b/FSI/flap_perp_2D/OpenFOAM-FEniCS/Allclean @@ -3,6 +3,9 @@ cd ${0%/*} || exit 1 # Run from this directory echo "Cleaning..." +# Source tutorial clean functions +. $WM_PROJECT_DIR/bin/tools/CleanFunctions + # Participant 1: Fluid Participant1="Fluid" cd ${Participant1} From 55c7f78d866b397950b4686cd52662073f3eaaf1 Mon Sep 17 00:00:00 2001 From: Ishaan Desai Date: Thu, 3 Feb 2022 13:11:57 +0100 Subject: [PATCH 3/3] Adding deleted file again --- perpendicular-flap/solid-fenics/solid.py | 229 +++++++++++++++++++++++ 1 file changed, 229 insertions(+) create mode 100644 perpendicular-flap/solid-fenics/solid.py diff --git a/perpendicular-flap/solid-fenics/solid.py b/perpendicular-flap/solid-fenics/solid.py new file mode 100644 index 000000000..66bcb3a62 --- /dev/null +++ b/perpendicular-flap/solid-fenics/solid.py @@ -0,0 +1,229 @@ +# Import required libs +from fenics import Constant, Function, AutoSubDomain, RectangleMesh, VectorFunctionSpace, interpolate, \ + TrialFunction, TestFunction, Point, Expression, DirichletBC, nabla_grad, project, \ + Identity, inner, dx, ds, sym, grad, lhs, rhs, dot, File, solve, PointSource, assemble_system +from ufl import nabla_div +import numpy as np +import matplotlib.pyplot as plt +from fenicsprecice import Adapter +from enum import Enum + + +# define the two kinds of boundary: clamped and coupling Neumann Boundary +def clamped_boundary(x, on_boundary): + return on_boundary and abs(x[1]) < tol + + +def neumann_boundary(x, on_boundary): + """ + determines whether a node is on the coupling boundary + + """ + return on_boundary and ((abs(x[1] - 1) < tol) or abs(abs(x[0]) - W / 2) < tol) + + +# Geometry and material properties +dim = 2 # number of dimensions +H = 1 +W = 0.1 +rho = 3000 +E = 4000000 +nu = 0.3 + +mu = Constant(E / (2.0 * (1.0 + nu))) + +lambda_ = Constant(E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))) + +# create Mesh +n_x_Direction = 4 +n_y_Direction = 26 +mesh = RectangleMesh(Point(-W / 2, 0), Point(W / 2, H), n_x_Direction, n_y_Direction) + +h = Constant(H / n_y_Direction) + +# create Function Space +V = VectorFunctionSpace(mesh, 'P', 2) + +# BCs +tol = 1E-14 + +# Trial and Test Functions +du = TrialFunction(V) +v = TestFunction(V) + +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) + + +f_N_function = interpolate(Expression(("1", "0"), degree=1), V) +u_function = interpolate(Expression(("0", "0"), degree=1), V) + +coupling_boundary = AutoSubDomain(neumann_boundary) +fixed_boundary = AutoSubDomain(clamped_boundary) + +precice = Adapter(adapter_config_filename="precice-adapter-config-fsi-s.json") + +# Initialize the coupling interface +precice_dt = precice.initialize(coupling_boundary, read_function_space=V, write_object=V, fixed_boundary=fixed_boundary) + +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])) + +# clamp the beam at the bottom +bc = DirichletBC(V, Constant((0, 0)), fixed_boundary) + +# alpha method parameters +alpha_m = Constant(0) +alpha_f = Constant(0) +gamma = Constant(0.5 + alpha_f - alpha_m) +beta = Constant((gamma + 0.5)**2 / 4.) + + +# 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 + + +# External Work +def Wext(u_): + return dot(u_, p) * ds + + +# Update functions + +# Update acceleration +def update_a(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 + - (1 - 2 * beta_) / 2 / beta_ * a_old) + + +# Update velocity +def update_v(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_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 + 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_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) + +a_form = lhs(res) +L_form = rhs(res) + +# parameters for Time-Stepping +t = 0.0 +n = 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 + +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 % 10 == 0: + displacement_out << (u_n, t) + +# Plot tip displacement evolution +displacement_out << u_n + +precice.finalize()