Skip to content

cherubimro/SWIRL2D

BranchesTags

Repository files navigation

2D Axisymmetric Swirling Flow Optimizer

License: GPL v3

Shape optimization of the inner (stagnation) boundary for 2D axisymmetric swirling flow using parallel FEM solver with multiple optimizers and spline interpolation methods.

This code is designed for analyzing and optimizing draft tube geometries to minimize vortex rope formation in hydraulic turbines operating at part-load conditions.

Table of Contents

Features

  • Parallel FEM Solver: MPI-parallelized finite element solver using PETSc
  • Multiple Optimizers: DONLP2 SQP (recommended), PETSc TAO, Python SciPy
  • 11 Spline Methods: GSL, PCHIP, FITPACK, and ALGLIB with various interpolation types
  • Configurable Boundary Conditions: Independent left/right BC control for ALGLIB cubic splines
  • Automated Benchmarking: Compare all optimizer/spline combinations automatically
  • Multiple Output Formats: Tecplot, VTK (ParaView), PNG plots, GAMBIT journal

Problem Description

This code solves an optimal shape design problem for swirling flow in a pipe with:

  • Inlet: z=0, r=[0.132, 1.0] (hub to shroud)
  • Outlet: z=4, r=[0.55, 1.1] (stagnant boundary to wall)
  • Inner boundary: Optimized shape connecting inlet hub (r=0.132) to outlet stagnant region (r=0.55)
  • Outer boundary: Conical expansion from r=1.0 to r=1.1

The objective is to minimize the flow force functional by finding the optimal shape of the inner (stagnation) boundary, parameterized by control points interpolated with splines.

Requirements

Build Requirements

Requirement Purpose Notes
PETSc 3.x Parallel FEM solver, TAO optimizer Must be configured with --download-triangle
GSL Spline interpolation GNU Scientific Library
MPI Parallel computing OpenMPI, MPICH, or Intel MPI
gfortran FITPACK compilation Any Fortran compiler
C++11 compiler ALGLIB compilation gcc 4.8+, clang 3.3+
CMake 3.14+ Build system (optional) Alternative to Make

Python Requirements (Optional)

For benchmarking and plotting scripts:

pip install numpy scipy matplotlib
# Or use the provided requirements file:
pip install -r requirements.txt

PETSc Configuration

PETSc must be installed with Triangle mesh library:

cd petsc
./configure --download-triangle [other options...]
make all

After building PETSc, set environment variables (add to ~/.bashrc):

export PETSC_DIR=/path/to/petsc
export PETSC_ARCH=arch-linux-c-debug  # Your configuration name

Quick Start

# Using CMake (recommended)
mkdir build && cd build
cmake ..
make -j4
cd ..
./build/generate_rpsi_data
mpirun -np 4 ./build/swirl2d_ipopt -nparam 2 -spline_type pchip

# Or using Make
make full
./generate_rpsi_data
mpirun -np 4 ./swirl2d_ipopt -nparam 2 -spline_type pchip

# Run comprehensive benchmark
python3 benchmark_all.py --quick

Building

Three build systems are available: CMake (recommended), Autoconf, and Make.

CMake Build (Recommended)

CMake automatically downloads ALGLIB and FITPACK if not present.

# Basic build
mkdir build && cd build
cmake ..
make -j4

# With custom PETSc location
cmake -DPETSC_DIR=/path/to/petsc -DPETSC_ARCH=arch-linux-c-opt ..
make -j4

CMake Options:

Option Default Description
PETSC_DIR $PETSC_DIR env Path to PETSc installation
PETSC_ARCH $PETSC_ARCH env PETSc configuration name
DOWNLOAD_ALGLIB ON Download ALGLIB if not present
DOWNLOAD_FITPACK ON Download FITPACK if not present
ALGLIB_VERSION 4.04.0 ALGLIB version to download

CMake Targets:

make                  # Build all executables
make swirl2d_ipopt    # Build optimizer only
make generate_rpsi_data  # Build data generator only
make run              # Generate data and run optimizer
make run4             # Run with 4 MPI processes
make benchmark        # Run quick benchmark (requires Python)
make benchmark_full   # Run full benchmark
make python-deps      # Install Python dependencies via pip

Autoconf Build

Standard GNU autoconf build with ./configure.

# Generate configure script (only needed once, or after changing configure.ac)
autoconf

# Configure (uses PETSC_DIR and PETSC_ARCH from environment)
./configure

# Or specify PETSc location explicitly
./configure --with-petsc-dir=/path/to/petsc --with-petsc-arch=arch-linux-c-opt

# Build
make

Configure Options:

Option Description
--with-petsc-dir=DIR PETSc installation directory
--with-petsc-arch=ARCH PETSc architecture name
--with-alglib=DIR ALGLIB source directory (default: alglib-cpp)
--with-alglib-version=VER ALGLIB version to download (default: 4.04.0)
--with-fitpack=DIR FITPACK source directory (default: fitpack)
--enable-download-alglib Download ALGLIB if not present (default: yes)
--enable-download-fitpack Download FITPACK if not present (default: yes)
--prefix=DIR Installation prefix (default: /usr/local)

Autoconf Targets:

make                  # Build all executables
make deps             # Download/build dependencies only
make run              # Generate data and run optimizer
make run4             # Run with 4 MPI processes
make benchmark        # Run quick benchmark
make install          # Install to prefix
make clean            # Remove build files
make distclean        # Remove all generated files
make help             # Show all targets

Make Build (Direct)

Traditional Makefile build that also downloads dependencies.

make all          # Build optimizer and data generator
make full         # Build all including FITPACK and ALGLIB
make fitpack      # Download and build FITPACK library
make alglib       # Build ALGLIB C++ wrapper
make donlp2       # Build DONLP2 SQP library
make clean        # Remove object files
make distclean    # Remove all generated files (including downloads)
make print        # Show build configuration
make help         # Show all targets

Build from Clean State

Both build systems will download ALGLIB and FITPACK automatically:

# CMake - clean build
rm -rf build alglib-cpp fitpack
mkdir build && cd build
cmake ..    # Downloads ALGLIB and FITPACK
make -j4

# Make - clean build
make distclean   # Removes alglib-cpp and fitpack directories
make full        # Downloads and builds everything

Switching Build Systems

After running ./configure, the original Makefile is backed up. To switch between systems:

# Use autoconf-generated Makefile (after ./configure)
make                    # Uses ./Makefile (generated by configure)

# Use original direct Makefile
make -f Makefile.direct # Uses the original Makefile without configure

# Restore original Makefile
cp Makefile.orig Makefile

Optimizers

Three optimization backends are available:

1. DONLP2 SQP (Default - Recommended)

Sequential Quadratic Programming optimizer. Fast and reliable for this problem.

# Sequential
./swirl2d_ipopt -nparam 2

# Parallel (4 MPI ranks for FEM)
mpirun -np 4 ./swirl2d_ipopt -nparam 2

# With options
mpirun -np 4 ./swirl2d_ipopt -nparam 5 -donlp2_maxiter 100

DONLP2 Options:

Option Description
-donlp2_maxiter N Maximum iterations (default: 100)

2. PETSc TAO

PETSc's optimization solvers. Good for comparison but may be slower.

# TAO BLMVM (bounded L-BFGS)
mpirun -np 4 ./swirl2d_ipopt -nparam 2 -tao -tao_type blmvm

# TAO with monitoring
mpirun -np 4 ./swirl2d_ipopt -tao -tao_type blmvm -tao_monitor -tao_view

TAO Options:

Option Description
-tao Enable TAO optimizer (disables DONLP2)
-tao_type TYPE blmvm, bqnls, bnls, nm (default: blmvm)
-tao_max_it N Maximum iterations
-tao_monitor Print iteration progress
-tao_view Print solver summary

3. Python SQP (External)

Uses scipy.optimize for optimization, calls FEM solver via subprocess.

python3 swirl2d_sqp.py --nparam 2 --mpi 4
python3 swirl2d_sqp.py --nparam 5 --mpi 4 --method L-BFGS-B

Spline Interpolation Methods

The inner boundary shape is defined by control points interpolated with splines. Eleven different spline methods are available:

GSL Splines (-spline_type gsl)

-spline_type gsl -gsl_type cubic    # Classic cubic spline (default)
-spline_type gsl -gsl_type akima    # Akima spline (local, needs 5+ points)
-spline_type gsl -gsl_type steffen  # Steffen spline (monotonic)

PCHIP (-spline_type pchip)

Piecewise Cubic Hermite Interpolating Polynomial. Shape-preserving, monotonicity-maintaining.

-spline_type pchip

FITPACK (-spline_type fitpack)

B-spline library from netlib. Supports smoothing and constrained fitting.

-spline_type fitpack -fitpack_type curfit   # Standard B-spline
-spline_type fitpack -fitpack_type concon   # Constrained spline (needs 4+ points)
-fitpack_smooth 0.001                        # Smoothing factor (0=interpolating)
-fitpack_order 3                             # Spline order (1-5, default: 3=cubic)

ALGLIB (-spline_type alglib)

C++ spline library with multiple algorithms and configurable boundary conditions.

-spline_type alglib -alglib_type cubic     # Cubic spline with configurable BC
-spline_type alglib -alglib_type akima     # Akima spline (needs 5+ points)
-spline_type alglib -alglib_type monotone  # Monotone spline (shape-preserving)
-spline_type alglib -alglib_type hermite   # Hermite spline (auto-derivatives)
-spline_type alglib -alglib_type catmullrom # Catmull-Rom (smooth through all points)

Boundary conditions for cubic spline:

# Set both endpoints to same BC type
-alglib_bc_type 0     # Parabolically-terminated (parabolic runout)
-alglib_bc_type 1     # First derivative specified
-alglib_bc_type 2     # Second derivative specified (default: natural spline)

# Set left/right BC types independently
-alglib_bc_type_left 0   # Left: parabolic
-alglib_bc_type_right 1  # Right: first derivative

# Boundary values (for bc_type 1 or 2)
-alglib_bc_left 0.0   # Left boundary value
-alglib_bc_right 0.0  # Right boundary value

Method Comparison

Method Min Points Monotonic Best For
GSL-cubic 3 No General purpose
GSL-akima 5 No Avoiding overshoot
GSL-steffen 3 Yes Few control points
PCHIP 3 Yes Best overall results
FITPACK-curfit 3 No Configurable smoothing
FITPACK-concon 4 Constrained Convexity control
ALGLIB-cubic 3 No Configurable BC
ALGLIB-akima 5 No Avoiding overshoot
ALGLIB-monotone 3 Yes Shape-preserving
ALGLIB-hermite 3 No Auto-estimated derivatives
ALGLIB-catmullrom 3 No Smooth through all points

Recommendation: Use PCHIP (-spline_type pchip) for best (sometimes) optimization results.

Benchmark Script

Automatically test all combinations of methods:

# Quick test (subset of methods, nparam 1-5)
python3 benchmark_all.py --quick

# Full benchmark (all methods, nparam 1-10)
python3 benchmark_all.py --nparam-max 10 -o results.csv

# Custom settings
python3 benchmark_all.py --nparam-max 8 --timeout 300

Benchmark output includes:

  • Best FFeval for each (spline, optimizer) combination
  • Best FFeval by number of parameters
  • Success rate for each method
  • Detailed CSV results file

Sample Results (nparam=2)

Method Optimizer FFeval
PCHIP DONLP2 -0.130732
PCHIP TAO-blmvm -0.130731
FITPACK-concon DONLP2 -0.130719
GSL-steffen DONLP2 -0.130637

Command-Line Reference

./swirl2d_ipopt [options]

General Options:
  -nparam N              Number of control points (default: 2)
  -seq_fem               Use sequential FEM (disable parallel)

Optimizer Selection:
  (default)              Use DONLP2 SQP
  -tao                   Use PETSc TAO optimizer

DONLP2 Options:
  -donlp2_maxiter N      Maximum iterations (default: 100)

TAO Options:
  -tao_type TYPE         blmvm, bqnls, bnls, nm
  -tao_max_it N          Maximum iterations
  -tao_monitor           Print progress
  -tao_view              Print summary

Spline Options:
  -spline_type TYPE      gsl, pchip, fitpack, alglib (default: gsl)
  -gsl_type TYPE         cubic, akima, steffen (default: cubic)
  -fitpack_type TYPE     curfit, concon (default: curfit)
  -fitpack_smooth S      Smoothing factor (default: 0.001)
  -fitpack_order K       Spline order 1-5 (default: 3)
  -alglib_type TYPE      cubic, akima, monotone, hermite, catmullrom
  -alglib_bc_type N      BC type for both ends: 0=parabolic, 1=1st deriv, 2=2nd deriv
  -alglib_bc_type_left N   Left BC type (overrides -alglib_bc_type)
  -alglib_bc_type_right N  Right BC type (overrides -alglib_bc_type)
  -alglib_bc_left V      Left boundary value (for bc_type 1 or 2)
  -alglib_bc_right V     Right boundary value (for bc_type 1 or 2)

Examples

# Basic optimization
mpirun -np 4 ./swirl2d_ipopt

# PCHIP spline with 5 control points
mpirun -np 4 ./swirl2d_ipopt -nparam 5 -spline_type pchip

# FITPACK with smoothing
mpirun -np 4 ./swirl2d_ipopt -spline_type fitpack -fitpack_smooth 0.01

# TAO optimizer with monitoring
mpirun -np 4 ./swirl2d_ipopt -tao -tao_type blmvm -tao_monitor

# GSL Steffen (monotonic) spline
mpirun -np 4 ./swirl2d_ipopt -spline_type gsl -gsl_type steffen

# ALGLIB monotone spline
mpirun -np 4 ./swirl2d_ipopt -spline_type alglib -alglib_type monotone

# ALGLIB Hermite spline
mpirun -np 4 ./swirl2d_ipopt -spline_type alglib -alglib_type hermite

# ALGLIB cubic with parabolic boundary conditions
mpirun -np 4 ./swirl2d_ipopt -spline_type alglib -alglib_type cubic -alglib_bc_type 0

# ALGLIB cubic with different left/right BC types
mpirun -np 4 ./swirl2d_ipopt -spline_type alglib -alglib_type cubic \
    -alglib_bc_type_left 0 -alglib_bc_type_right 1 -alglib_bc_right 0.1

# Run benchmark
python3 benchmark_all.py --quick -o results.csv

Output Files

File Description
meshoutSQP.dat Tecplot mesh with PSI solution (DONLP2, Python SQP)
meshout.dat Tecplot mesh with PSI solution (TAO)
params.dat Optimization log with all function evaluations
gambit.jou GAMBIT journal with optimized boundary vertices
solutionSQP.vtk VTK format for ParaView visualization
solutionSQP.png PSI contour plot (generated by plot_psi.py)
inlet-r-psi.dat Inlet boundary conditions
outlet-r-psi.dat Outlet boundary conditions
inner_flux.dat Boundary flux on inner (stagnation) boundary
outer_flux.dat Boundary flux on outer (wall) boundary
benchmark_results.csv Benchmark results (from benchmark_all.py)

Note: DONLP2 and Python SQP write to meshoutSQP.dat, while TAO writes to meshout.dat.

Code Structure

.
├── CMakeLists.txt         # CMake build system (recommended)
├── configure.ac           # Autoconf input file
├── Makefile.ac.in         # Autoconf Makefile template
├── configure              # Generated configure script (run autoconf)
├── Makefile.orig          # Original Make build system
├── Makefile               # Active Makefile (from configure or original)
├── README.md              # This documentation
├── requirements.txt       # Python dependencies
├── install-sh             # Autoconf helper script
│
├── swirl2d_ipopt.c        # Main: Parallel FEM solver with DONLP2/TAO optimizers
├── generate_rpsi_data.c   # Analytical boundary condition generator
├── alglib_wrapper.cpp     # C++ wrapper for ALGLIB spline interpolation
│
├── benchmark_all.py       # Automated benchmarking of all method combinations
├── swirl2d_sqp.py         # Python SQP optimizer wrapper (uses scipy.optimize)
├── plot_psi.py            # PSI contour visualization
│
├── build/                 # CMake build directory (created by user)
│   ├── swirl2d_ipopt      # Built executable
│   └── generate_rpsi_data # Built data generator
│
├── fitpack/               # FITPACK B-spline library (auto-downloaded)
├── donlp2c/               # DONLP2 SQP optimizer library
│   ├── donlp2.c           # DONLP2 core algorithm
│   └── donlp2_compat.h    # Compatibility header
├── alglib-cpp/            # ALGLIB C++ spline library (auto-downloaded)
│   └── src/               # ALGLIB source files
└── ../src/                # Shared sources
    └── dnconf_wrapper.c   # DONLP2 configuration wrapper

Key Source Files

File Description
swirl2d_ipopt.c Main program: mesh generation, FEM assembly, optimization loop
generate_rpsi_data.c Generates inlet/outlet BCs from abc-swirl analytical solution
alglib_wrapper.cpp Exposes ALGLIB splines to C code via extern "C" interface
benchmark_all.py Tests all spline/optimizer combinations, generates CSV reports
swirl2d_sqp.py Alternative optimizer using Python scipy.optimize.minimize

Troubleshooting

Triangle Mesh Error

If you see "Internal error in finddirection()" with FITPACK, this is a numerical precision issue in the Triangle library. The code includes a workaround (tiny perturbation) but some parameter combinations may still trigger it. Use PCHIP or GSL-steffen instead.

GSL Akima Requires 5 Points

GSL and ALGLIB Akima splines require at least 5 data points (nparam >= 3). The code automatically falls back to Steffen if fewer points are provided.

FITPACK concon Requires 4 Points

FITPACK concon (constrained spline) requires m > 3, meaning nparam >= 2.

Physics

The code solves the Bragg-Hawthorne equation for axisymmetric swirling flow:

-d/dz(1/r * dpsi/dz) - d/dr(1/r * dpsi/dr) = r*H'(psi) - K(psi)*K'(psi)/r

where:

  • psi: stream function (constant on streamlines)
  • H(psi): total head distribution (Bernoulli function)
  • K(psi): circulation distribution (r * v_theta)

Boundary Conditions

Boundary conditions are derived from the abc-swirl analytical solution:

  • Inlet (z=0): Prescribed psi profile from hub to shroud
  • Outlet (z=4): Prescribed psi profile from stagnation region to wall
  • Inner boundary: Stagnation streamline (psi=0)
  • Outer boundary: Wall streamline (psi=psi_wall)

Objective Function

The optimization minimizes the flow force functional (FFeval), which represents the integrated pressure and momentum flux on the inner boundary. Lower values indicate better flow conditions with reduced recirculation and vortex rope formation.

Acknowledgements

  • PETSc team for the parallel linear algebra framework
  • Jonathan Shewchuk for the Triangle mesh generator
  • P. Spellucci for the original DONLP2 optimizer (Fortran)
  • S. Schoeffert (ASB, Bourges, France) for the DONLP2 C translation
  • Diethelm Wuertz and Ryuichi Tamura for the Rdonlp2 R package
  • Paul Dierckx for the FITPACK library
  • ALGLIB project for the spline interpolation library

Third-Party Libraries

This project uses the following external libraries:

Library Purpose License
PETSc Parallel linear algebra, TAO optimizer BSD-2-Clause
Triangle Mesh generation Non-commercial use
GSL Spline interpolation GPL-3.0
FITPACK B-spline fitting Public domain (Netlib)
ALGLIB Spline interpolation GPL-2.0+
DONLP2 SQP optimizer (C code extracted from R package Rdonlp2) Non-commercial use

Note: DONLP2 and Triangle have non-commercial use restrictions. For commercial applications, contact the respective authors for licensing terms.

License

This project is licensed under the GNU General Public License v3.0 - see the LICENSE file for details.

Citing This Work

If you use SWIRL2D in your research, you MUST cite the following paper:

Anton, A., Muntean, S., & Susan-Resiga, R. (2016). SWIRL2D: An interface tracking algorithm for computing the two-dimensional swirling flows with stagnant region. Proceedings of the Romanian Academy, Series A, 17(4), 379-386. PDF

@article{Anton2016swirl2d,
  author    = {Anton, Alin and Muntean, Sebastian and Susan-Resiga, Romeo},
  title     = {{SWIRL2D}: An interface tracking algorithm for computing the
               two-dimensional swirling flows with stagnant region},
  journal   = {Proceedings of the Romanian Academy, Series A - Mathematics,
               Physics, Technical Sciences, Information Science},
  volume    = {17},
  number    = {4},
  pages     = {379--386},
  year      = {2016},
  url       = {https://acad.ro/sectii2002/proceedings/doc2016-4/12ProcRoAcad-4-2016.pdf}
}

Recommended Additional Citations

Theoretical background (IAHR Working Group presentation, included in repository):

@inproceedings{SusanResiga2013IAHRWG,
  author    = {Susan-Resiga, Romeo and Muntean, Sebastian and Anton, Alin},
  title     = {Swirling Flow Analysis in a Francis Turbine Draft Tube},
  booktitle = {6th IAHR International Meeting of the Workgroup on
               Cavitation and Dynamic Problems in Hydraulic Machinery and Systems},
  year      = {2013},
  address   = {Ljubljana, Slovenia},
  note      = {PDF included: Resiga2013IAHRWG.pdf}
}

Theoretical background (variational model):

@article{SusanResiga2016variational,
  author  = {Susan-Resiga, R.F. and Muntean, S. and Stuparu, A. and
             Bosioc, A.I. and Tanasa, C. and Ighisan, C.},
  title   = {A variational model for swirling flow states with stagnant region},
  journal = {European Journal of Mechanics - B/Fluids},
  volume  = {55}, pages = {104--115}, year = {2016},
  doi     = {10.1016/j.euromechflu.2015.09.002}
}

PETSc (parallel solver and TAO optimizer):

@techreport{PETSc2024,
  author      = {Balay, S. and Abhyankar, S. and Adams, M.F. and others},
  title       = {{PETSc/TAO} Users Manual},
  institution = {Argonne National Laboratory},
  number      = {ANL-21/39 - Revision 3.22}, year = {2024},
  url         = {https://petsc.org/}
}

Triangle (mesh generation):

@inproceedings{Shewchuk1996triangle,
  author    = {Shewchuk, Jonathan Richard},
  title     = {Triangle: Engineering a {2D} quality mesh generator and
               {Delaunay} triangulator},
  booktitle = {Applied Computational Geometry: Towards Geometric Engineering},
  series    = {LNCS}, volume = {1148}, pages = {203--222},
  year      = {1996}, publisher = {Springer-Verlag}
}

FITPACK (B-spline fitting):

@book{Dierckx1993fitpack,
  author    = {Dierckx, Paul},
  title     = {Curve and Surface Fitting with Splines},
  publisher = {Oxford University Press}, year = {1993},
  note      = {FITPACK library: https://www.netlib.org/dierckx/}
}

DONLP2 (SQP optimizer):

@techreport{Spellucci1998donlp2,
  author      = {Spellucci, Peter},
  title       = {{DONLP2} -- Dynamic constrained nonlinear optimization},
  institution = {TU Darmstadt}, year = {1998},
  note        = {C code from Rdonlp2 R package}
}

GSL (spline interpolation):

@book{Galassi2009gsl,
  author    = {Galassi, M. and Davies, J. and Theiler, J. and others},
  title     = {{GNU} Scientific Library Reference Manual},
  publisher = {Network Theory Ltd.}, edition = {3rd}, year = {2009},
  url       = {https://www.gnu.org/software/gsl/}
}

ALGLIB (spline interpolation):

@misc{ALGLIB,
  author = {Bochkanov, Sergey},
  title  = {{ALGLIB} -- Numerical analysis library},
  year   = {2024}, url = {https://www.alglib.net/}
}

A complete BibTeX file is available in references.bib.

Contributing

Contributions are welcome! Please feel free to submit issues and pull requests.

About

2D Swirling Flow Solver using SQP Optimization with IPOPT/DONLP2

Resources

Readme

License

GPL-3.0 license
Activity

Stars

0 stars

Watchers

0 watching

Forks

0 forks
Report repository

Releases 1

Initial Release Latest
Feb 7, 2026

Packages

 
 
 

Contributors

Languages