Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add solid fake for perpendicular flap tutorial #541

Merged
merged 8 commits into from
Aug 19, 2024
Merged
Show file tree
Hide file tree
Changes from 6 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions perpendicular-flap/README.md
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,8 @@ Solid participant:

* OpenFOAM (solidDisplacementFoam). For more information, have a look at the [OpenFOAM plateHole tutorial](https://www.openfoam.com/documentation/tutorial-guide/5-stress-analysis/5.1-stress-analysis-of-a-plate-with-a-hole). The solidDisplacementFoam solver only supports linear geometry and this case is only provided for quick testing purposes, leading to outlier results. For general solid mechanics procedures in OpenFOAM, see solids4foam.

* Fake. A simple Python script that acts as a fake solver and provides an arbitrary time-dependent flap displacement in the x-direction, i.e., it performs a shear mapping on the resting flap. This solver can be used for debugging of the fluid participant and its adapter. It also technically works with implicit coupling, thus no changes to the preCICE configuration are necessary. Note that [ASTE's replay mode](https://precice.org/tooling-aste.html#replay-mode) has a similar use case and could also feed artificial or previously recorded real data, replacing an actual solver.

## Running the Simulation

All listed solvers can be used in order to run the simulation. OpenFOAM can be executed in parallel using `run.sh -parallel`. The default setting uses 4 MPI ranks. Open two separate terminals and start the desired fluid and solid participant by calling the respective run script `run.sh` located in the participant directory. For example:
Expand Down
76 changes: 76 additions & 0 deletions perpendicular-flap/solid-fake/fake.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,76 @@
from __future__ import division

import numpy as np
import precice
NiklasVin marked this conversation as resolved.
Show resolved Hide resolved


def displace_flap(x, y, t, flap_tip_y):
x_displ = np.zeros_like(x)
y_displ = np.zeros_like(y)
# first, get displacement independent of x, only dependent on y and t
# maximal displacement at flap tip should be 0.5
# initially, flap's displacement is 0
x_displ = np.minimum(((np.sin(3 * np.pi * t + np.arcsin(-0.95)) + 0.95) / 8) * y / flap_tip_y, 0.5 * y / flap_tip_y)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
# first, get displacement independent of x, only dependent on y and t
# maximal displacement at flap tip should be 0.5
# initially, flap's displacement is 0
x_displ = np.minimum(((np.sin(3 * np.pi * t + np.arcsin(-0.95)) + 0.95) / 8) * y / flap_tip_y, 0.5 * y / flap_tip_y)
# first, get displacement independent of x, only dependent on y and t
# initially, flap's displacement is 0
max_x_displacement = 0.5
x_displ = np.minimum(((np.sin(3 * np.pi * t + np.arcsin(-0.95)) + 0.95) / 8) * y / flap_tip_y, max_x_displacement * y / flap_tip_y)

I would generally try to avoid raw numbers and rather use descriptive names. This often allows you to completely save some comments.

What does the 0.95, the 3 * np.pi and the /8 do? Generally, I think if the formula you are using here works we can keep it. But a brief comment about the general rationale behind the arcsin and sin would be helpful.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I added some comments that hopefully clarify the idea behind the formula


displ = np.zeros((len(x), 2))
NiklasVin marked this conversation as resolved.
Show resolved Hide resolved
displ[:, 0] = x_displ
displ[:, 1] = y_displ

return displ


configuration_file_name = "../precice-config.xml"
participant_name = "Solid"
mesh_name = "Solid-Mesh"
write_data_name = 'Displacement'

solver_process_index = 0
solver_process_size = 1

# define mesh
H = 1
W = 0.1

interface = precice.Participant(participant_name, configuration_file_name, solver_process_index, solver_process_size)
dimensions = interface.get_mesh_dimensions(mesh_name)
assert (dimensions == 2)

x_left = 0.0 - 0.5 * W # left boundary of the flap
x_right = 0.5 * W # right boundary of the flap
y_bottom = 0.0 # bottom of the flap
y_top = y_bottom + H # top of the flap

n = 24 # Number of vertices per side
t = 0

vertices_mid = np.zeros((2 * n, dimensions))
vertices_mid[:n, 1] = np.linspace(y_bottom, y_top, n)
vertices_mid[n:, 1] = np.linspace(y_bottom, y_top, n)

vertex_ids = interface.set_mesh_vertices(mesh_name, vertices_mid)
MakisH marked this conversation as resolved.
Show resolved Hide resolved

interface.initialize()
NiklasVin marked this conversation as resolved.
Show resolved Hide resolved
# change if necessary
solver_dt = np.inf
# for checkpointing
t_cp = 0

while interface.is_coupling_ongoing():
if interface.requires_writing_checkpoint():
t_cp = t

precice_dt = interface.get_max_time_step_size()
dt = min([solver_dt, precice_dt])
# wiggle the flap
write_data = displace_flap(vertices_mid[:, 0], vertices_mid[:, 1], t, H)

interface.write_data(mesh_name, write_data_name, vertex_ids, write_data)
interface.advance(dt)

if interface.requires_reading_checkpoint():
t = t_cp
else:
# update t
t += dt

interface.finalize()