johanneshaubner/shapeopt

ShapeOpt

ShapeOpt is based on FEniCS, dolfin-adjoint, IPOPT, cyipopt, pygmsh.

Code repository for the manuscript

J. Haubner, M. Ulbrich: Numerical Methods for Shape Optimal Design of Fluid-Structure Interaction Problems.

based on the implementation for the PhD thesis

J. Haubner: Shape Optimization for Fluid-Structure Interaction, Doctoral Dissertation, Technische Universität München, 2020

Usage/Examples

Using Docker image

The Dockerfile can be used by running:

docker build -t shapeopt .
docker run -it shapeopt

or

docker pull ghcr.io/johanneshaubner/shapeopt:latest
docker run -ti -v ${PWD}:/root/shared -w /root/shared --entrypoint=/bin/bash --rm ghcr.io/johanneshaubner/shapeopt:latest

IPOPT with HSL

For runs in parallel, IPOPT needs to be installed with HSL. To do so, the Coin-HSL sources archive needs to be added such that it becomes shapeopt/hsl/coinhsl. Afterwards we build the docker image via Dockerfile_hsl.

If not ran from Docker image: Requires a recent master version of dolfin with MeshView support. Requires the changes propsed in https://bitbucket.org/fenics-project/dolfin/issues/1123/assemble-on-mixed-meshview-forms-returns.

Running Example from Paper

To create the mesh use

python3 example/FSI/create_mesh_FSI.py

To obtain the results presented in the paper we used an IPOPT installation with HSL and ran

mpiexec -n 4 python3 example/FSI/main.py >> example/FSI/mesh/Output/terminal.txt

(FFC Compiler seems to get stuck in the beginning, cancelling when this happens and restarting the simulation seems to fix this.)

To obtain the tables and figures we ran (with setting the corresponding options to True in the script)

mpiexec -n 4 python3 example/FSI/visualize.py

Running Tests

To run tests, run the following command

pytest

License

GNU GPLv3

Structure

Control_to_Trafo

We apply a method of mappings approach and solve everything on a reference domain. In order not to have too much redundancy in the representation of the transformations, we choose a scalar valued function on the design boundary as control for the optimization problem and define the deformation field using an extension operator from the design boundary to the whole domain. In the example code, we define the Control_to_Trafo-mapping as a composition of boundary operators and extension operators.

Reduced_Objective

After having defined the transformation, the reduced objective and its gradient (computed using dolfin-adjoint) can be evaluated. An example for a simple Stokes flow and Fluid-Structure Interaction is given.

Constraints

Collection of constraint for the optimization problem: volume, barycenter, determinant.

Mesh_Postprocessing

Postprocessing step for improving the mesh quality after each optimization problem solve.

Ipopt

Solving the constraint optimization problem is performed using Ipopt.

Tools

Here, different tools are collect, e.g. initialization of function spaces, save and load objects, etc.

Authors

• Johannes Haubner
• Michael Ulbrich

Acknowledgement

We would like to acknowledge Jørgen S. Dokken, Henrik N. Finsberg and Simon W. Funke for the help, support and discussions on reproducibility.

shapeopt's Issues

Speed-up linear solvers

instead of solve(a== L, u, bcs) define the solver (e.g., boundary operators)
self.lb_off should not be a float but a Constant and be updated with .values = .. (then, the form is not compiled again)

check extension operator

ShapeOpt/src/shapeopt/Control_to_Trafo/Extension_Operator/ElasticExtension.py

Line 80 in c4ea234

bc = [bc1, bc2, bc3]

seems to be wrong;
"bc" as in eval?

docker image not working

TODO: docker image uses installation with different blas and lapack for dolfin and IPOPT. write new docker image based on https://github.com/scientificcomputing/packages/blob/main/fenics/Dockerfile

Change evaluation routine for IPOPT objective and gradient

evaluation routine for objective/gradient:

• first evaluate the maximal value of the determinant of the gradient of (identity + deformation)
• if mesh degeneration; do not compute FSI solution (will fail), return large value instead (heuristic, but with additional penalization term for determinant of gradient of (identity + deformation) it should be OK; will probably mostly be needed in initial stage of algorithm)

fix path in save_output

ShapeOpt/src/shapeopt/Reduced_Objective/FluidStructure.py

Line 351 in 0e35cbb

J += assemble(0.5*Constant(param["gammaP"]) * smoothmax(Constant(param["etaP"]) - tJhat)**2*dx(mesh))

replace by 1.0/* under the integral:
assemble(0.5*Constant(param["gammaP"]) * 1.0/((tJhat- Constant(param["etaP"]) )**2)*dx(mesh))

ShapeOpt/src/shapeopt/Ipopt/ipopt_solver.py

Line 127 in 0e35cbb

if dgt.vector().min() < 0.5*self.param['det_lb']:

replace by <= 1.0*self.param['det_lb']

set gammaP to smaller number (1e-1/1e-2 e.g.)

Fix run.sh-file to run docker-image with singularity