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

Improve elastic-tube-1d code and documentation #209

Merged
merged 12 commits into from
Apr 28, 2021
16 changes: 6 additions & 10 deletions elastic-tube-1d/README.md
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,8 @@ The following parameters have been chosen:
- Fluid density: $$ \rho = 1 $$
- Young modulus: E = 10000

Additionally the solvers use the parameters `N = 100`, `tau = 0.01`, `kappa = 100` by default. These values can be modified directly in each solver.

## Available solvers

Both fluid and solid participant are supported in:
Expand Down Expand Up @@ -70,8 +72,6 @@ cd solid-cpp
./run.sh
```

The solvers use the parameters `N = 100`, `tau = 0.01`, `kappa = 100` by default and can be modified in the solver.

### Python

Open two separate terminals and start each participant by calling the respective run script. Only serial run is possible:
Expand All @@ -88,25 +88,21 @@ cd solid-python
./run.sh
```

Parameters such as `N` can be modified directly at the `FluidSolver.py` and at the `SolidSolver.py`. The parameters must be consistent between the different solvers and participants.

**Optional:** Visualization and video output of the fluid participant can be triggered via the options `--enable-plot` and `--write-video` of `FluidSolver.py`. To generate .vtk files during execution, you need to add the flag `--write-vtk`.
IshaanDesai marked this conversation as resolved.
Show resolved Hide resolved

![Elastic tube animation](images/tutorials-elastic-tube-1d-animation.gif)
**Optional:** A run-time plot visualization can be trigged by passing `--enable-plot` in `run.sh` of `FluidSolver.py`. Additionally a video of the run-time plot visualization can be generated by additionally passing `--write-video`

{% include warning.html content= "The C++ and Python solvers lead to different results. Please consider the Python results as the correct ones and refer to this [open issue](https://github.com/precice/tutorials/issues/195) for more insight. Contributions are particularly welcome here." %}

## Post-processing

![Elastic tube animation](images/tutorials-elastic-tube-1d-animation.gif)

The results from each simulation are stored in each `fluid-<participant>/output/` folder. You can visualize these VTK files using the provided `plot-diameter.sh` script

```bash
./plot-diameter.sh
```

which will try to visualize the results from both fluid cases, if available.

This script calls the more flexible `plot-vtk.py` Python script, which you can use as
which will try to visualize the results from both fluid cases, if available. This script calls the more flexible `plot-vtk.py` Python script, which you can use as

```bash
python3 plot-vtk.py <quantity> <case>/output/<prefix>
Expand Down
7 changes: 2 additions & 5 deletions elastic-tube-1d/fluid-cpp/src/FluidSolver.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -72,11 +72,9 @@ int main(int argc, char **argv)

interface.setMeshVertices(meshID, chunkLength, grid.data(), vertexIDs.data());

std::cout << "Initialize preCICE..." << std::endl;
interface.initialize();

double t = 0.0;
double dt = 0.01;
MakisH marked this conversation as resolved.
Show resolved Hide resolved
std::cout << "Initialize preCICE..." << std::endl;
double dt = interface.initialize();

if (interface.isActionRequired(actionWriteInitialData())) {
interface.writeBlockScalarData(pressureID, chunkLength, vertexIDs.data(), pressure.data());
Expand All @@ -102,7 +100,6 @@ int main(int argc, char **argv)
if (interface.isActionRequired(actionWriteIterationCheckpoint())) {
interface.markActionFulfilled(actionWriteIterationCheckpoint());
}


if (interface.isReadDataAvailable()) {
interface.readBlockScalarData(crossSectionLengthID, chunkLength, vertexIDs.data(), crossSectionLength.data());
Expand Down
75 changes: 32 additions & 43 deletions elastic-tube-1d/fluid-python/FluidSolver.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,20 +6,40 @@
from thetaScheme import perform_partitioned_implicit_trapezoidal_rule_step, perform_partitioned_implicit_euler_step
import numpy as np
import tubePlotting

import matplotlib.pyplot as plt
import matplotlib.animation as manimation

from output import writeOutputToVTK

import precice
from precice import *
from precice import action_write_initial_data, action_write_iteration_checkpoint, \
action_read_iteration_checkpoint

# physical properties of the tube
r0 = 1 / np.sqrt(np.pi) # radius of the tube
a0 = r0**2 * np.pi # cross sectional area
u0 = 10 # mean velocity
ampl = 3 # amplitude of varying velocity
frequency = 10 # frequency of variation
t_shift = 0 # temporal shift of variation
p0 = 0 # pressure at outlet
kappa = 100

L = 10 # length of tube/simulation domain
N = 100
dx = L / kappa
# helper function to create constant cross section


def velocity_in(t): return u0 + ampl * np.sin(frequency *
(t + t_shift) * np.pi) # inflow velocity


def crossSection0(N):
return a0 * np.ones(N + 1)


parser = argparse.ArgumentParser()
parser.add_argument("configurationFileName", help="Name of the xml precice configuration file.",
nargs='?', type=str, default="../precice-config.xml")
parser.add_argument(
"--write-vtk", help="Save vtk files of each timestep in the 'VTK' folder.", action='store_true')
parser.add_argument(
"--enable-plot", help="Show a continuously updated plot of the tube while simulating.", action='store_true')
parser.add_argument("--write-video", help="Save a video of the simulation as 'writer_test.mp4'. \
Expand All @@ -33,7 +53,6 @@
print("Try '$ python FluidSolver.py precice-config.xml'")
quit()

output_mode = config.OutputModes.VTK if args.write_vtk else config.OutputModes.OFF
plotting_mode = config.PlottingModes.VIDEO if args.enable_plot else config.PlottingModes.OFF
if args.write_video and not args.enable_plot:
print("")
Expand All @@ -42,41 +61,14 @@
quit()
writeVideoToFile = True if args.write_video else False

print("Starting Fluid Solver...")

configFileName = args.configurationFileName

# physical properties of the tube
r0 = 1 / np.sqrt(np.pi) # radius of the tube
a0 = r0**2 * np.pi # cross sectional area
u0 = 10 # mean velocity
ampl = 3 # amplitude of varying velocity
frequency = 10 # frequency of variation
t_shift = 0 # temporal shift of variation
p0 = 0 # pressure at outlet
kappa = 100


def velocity_in(t): return u0 + ampl * np.sin(frequency *
(t + t_shift) * np.pi) # inflow velocity


L = 10 # length of tube/simulation domain
N = 100
dx = L / kappa
# helper function to create constant cross section


def crossSection0(N):
return a0 * np.ones(N + 1)
print("Plotting Mode: {}".format(plotting_mode))

print("Starting Fluid Solver...")

print("N: " + str(N))

solverName = "Fluid"

print("Configure preCICE...")
interface = precice.Interface(solverName, configFileName, 0, 1)
interface = precice.Interface("Fluid", args.configurationFileName, 0, 1)
print("preCICE configured...")

dimensions = interface.get_dimensions()
Expand All @@ -88,7 +80,6 @@ def crossSection0(N):
crossSectionLength = a0 * np.ones(N + 1)
crossSectionLength_old = a0 * np.ones(N + 1)


if plotting_mode == config.PlottingModes.VIDEO:
fig, ax = plt.subplots(1)
if writeVideoToFile:
Expand All @@ -110,11 +101,10 @@ def crossSection0(N):
vertexIDs = interface.set_mesh_vertices(meshID, grid)

t = 0
precice_dt = 0.01

print("Fluid: init precice...")
# preCICE defines timestep size of solver via precice-config.xml
interface.initialize()
precice_dt = interface.initialize()

if interface.is_action_required(action_write_initial_data()):
interface.write_block_scalar_data(pressureID, vertexIDs, pressure)
Expand Down Expand Up @@ -161,9 +151,8 @@ def crossSection0(N):
velocity_old = np.copy(velocity)
pressure_old = np.copy(pressure)
crossSectionLength_old = np.copy(crossSectionLength)
if output_mode is config.OutputModes.VTK:
writeOutputToVTK(time_it, "out_fluid_", dx, datanames=["velocity", "pressure", "diameter"], data=[
velocity_old, pressure_old, crossSectionLength_old])
writeOutputToVTK(time_it, "out_fluid_", dx, datanames=["velocity", "pressure", "diameter"], data=[
velocity_old, pressure_old, crossSectionLength_old])
time_it += 1

print("Exiting FluidSolver")
Expand Down
2 changes: 1 addition & 1 deletion elastic-tube-1d/fluid-python/run.sh
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
#!/bin/sh
set -e -u

python3 ./FluidSolver.py ../precice-config.xml --write-vtk
python3 ./FluidSolver.py ../precice-config.xml
6 changes: 4 additions & 2 deletions elastic-tube-1d/plot-diameter.sh
Original file line number Diff line number Diff line change
@@ -1,15 +1,17 @@
#!/bin/sh
set -e -u

# File check solution from: https://stackoverflow.com/questions/91368/checking-from-shell-script-if-a-directory-contains-files

# Plot diameter from fluid-cpp
if [ -d "./fluid-cpp/output/" ]; then
if [ -n "$(ls -A ./fluid-cpp/output/*.vtk 2>/dev/null)" ]; then
MakisH marked this conversation as resolved.
Show resolved Hide resolved
python3 plot-vtk.py diameter fluid-cpp/output/out_fluid_ &
else
echo "No results to plot from fluid-cpp."
fi

# Plot diameter from fluid-python
if [ -d "./fluid-python/output/" ]; then
if [ -n "$(ls -A ./fluid-python/output/*.vtk 2>/dev/null)" ]; then
python3 plot-vtk.py diameter fluid-python/output/out_fluid_ &
else
echo "No results to plot from fluid-python."
Expand Down
5 changes: 2 additions & 3 deletions elastic-tube-1d/solid-cpp/src/SolidSolver.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -45,11 +45,10 @@ int main(int argc, char **argv)

std::vector<int> vertexIDs(chunkLength);
interface.setMeshVertices(meshID, chunkLength, grid.data(), vertexIDs.data());
std::cout << "Initialize preCICE..." << std::endl;
interface.initialize();

double t = 0;
double dt = 0.01;
std::cout << "Initialize preCICE..." << std::endl;
double dt = interface.initialize();

if (interface.isActionRequired(actionWriteInitialData())) {
interface.writeBlockScalarData(crossSectionLengthID, chunkLength, vertexIDs.data(), crossSectionLength.data());
Expand Down
11 changes: 3 additions & 8 deletions elastic-tube-1d/solid-python/SolidSolver.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,8 @@
import argparse
import numpy as np
import precice
from precice import *
from precice import action_write_initial_data, action_read_iteration_checkpoint, \
action_write_iteration_checkpoint


r0 = 1 / np.sqrt(np.pi) # radius of the tube
Expand All @@ -21,8 +22,6 @@
def crossSection0(N):
return a0 * np.ones(N + 1)

###############


print("Starting Solid Solver...")

Expand All @@ -38,14 +37,10 @@ def crossSection0(N):
print("Try '$ python SolidSolver.py precice-config.xml'")
quit()

configFileName = args.configurationFileName

print("N: " + str(N))

solverName = "Solid"

print("Configure preCICE...")
interface = precice.Interface(solverName, configFileName, 0, 1)
interface = precice.Interface("Solid", args.configurationFileName, 0, 1)
print("preCICE configured...")

dimensions = interface.get_dimensions()
Expand Down