port2python.txt

================================================================================
PORTING SWIRL2D TO 100% PYTHON - FEASIBILITY ANALYSIS
================================================================================

Author: Analysis generated from codebase review
Date: 2026-01-23

================================================================================
1. CURRENT DEPENDENCIES
================================================================================

1.1 Triangle Library
--------------------
Location: #include <triangle.h>
Purpose: 2D Delaunay triangulation and constrained mesh generation

Usage in swirl2d_ipopt.c:
  - struct triangulateio for mesh input/output
  - triangulate() function with switches "pzqYQa0.005"
    - p: PSLG (Planar Straight Line Graph) input
    - z: zero-based indexing
    - q: quality mesh generation
    - Y: no new vertices on boundary
    - Q: quiet mode
    - a0.005: maximum triangle area constraint

Data structures used:
  - numberofpoints, numberoftriangles
  - pointlist (x,y coordinates)
  - trianglelist (vertex indices per triangle)
  - pointmarkerlist (boundary markers)
  - segmentlist, segmentmarkerlist
  - pointattributelist (solution values)

1.2 PETSc Library
-----------------
Location: #include <petsctao.h>
Purpose: Parallel sparse linear algebra and optimization

Components used:

a) Sparse Matrices (Mat):
   - MatCreateAIJ() - parallel AIJ (CSR) matrix
   - MatSetValues() - element assembly
   - MatAssemblyBegin/End() - finalize assembly
   - MatZeroRows() - apply Dirichlet BC
   - MatTranspose() - for boundary flux computation
   - MatDuplicate(), MatDestroy()

b) Vectors (Vec):
   - VecCreateMPI() - parallel vector
   - VecSet(), VecSetValues()
   - VecAssemblyBegin/End()
   - VecCopy(), VecScale()
   - VecScatterCreateToZero() - gather to rank 0
   - VecGetArray/RestoreArray()

c) Linear Solver (KSP):
   - KSPCreate() - Krylov subspace solver
   - KSPSetOperators() - set matrix
   - KSPSetFromOptions() - command-line configuration
   - KSPSolve() - solve Ax=b
   - KSPGetIterationNumber()

d) TAO Optimizer:
   - TaoCreate(), TaoSetType()
   - TaoSetObjectiveAndGradientRoutine()
   - TaoSetVariableBounds()
   - TaoSolve()

e) MPI Parallelism:
   - MPI_Comm_rank/size
   - MPI_Bcast() - broadcast mesh data
   - Round-robin element distribution
   - VecScatter for solution gathering

f) Utilities:
   - PetscMalloc/Free - memory management
   - PetscFOpen/Close/Printf - file I/O
   - PetscOptionsGetInt/Real/String - command-line parsing

1.3 Other C Dependencies
------------------------
- GSL (gsl_spline.h, gsl_sf_bessel.h): Spline interpolation, Bessel functions
- FITPACK (Fortran): B-spline fitting
- ALGLIB (C++): Additional spline methods
- DONLP2 (C): SQP optimizer

================================================================================
2. PYTHON REPLACEMENT LIBRARIES
================================================================================

2.1 Mesh Generation
-------------------
Current: Triangle C library
Replace with: triangle Python package (pip install triangle)

The 'triangle' package is a direct Python wrapper around Shewchuk's Triangle.
It provides identical functionality with pythonic interface.

Example:
    import triangle

    # Define PSLG
    A = {
        'vertices': np.array([[x1,y1], [x2,y2], ...]),
        'segments': np.array([[0,1], [1,2], ...]),
        'segment_markers': np.array([marker1, marker2, ...])
    }

    # Generate mesh (equivalent to "pzqYQa0.005")
    B = triangle.triangulate(A, 'pqa0.005')

    # Access results
    vertices = B['vertices']        # (N, 2) array
    triangles = B['triangles']      # (M, 3) array
    markers = B['vertex_markers']   # (N,) array

Alternative: meshpy package (pip install meshpy)

2.2 Sparse Linear Algebra
-------------------------
Current: PETSc Mat, Vec, KSP
Replace with: scipy.sparse and scipy.sparse.linalg

Matrix creation:
    from scipy.sparse import lil_matrix, csr_matrix

    # Assembly (lil_matrix is efficient for incremental construction)
    K = lil_matrix((n, n), dtype=float)
    for elem in elements:
        ke = element_stiffness(...)
        for i, gi in enumerate(elem):
            for j, gj in enumerate(elem):
                K[gi, gj] += ke[i, j]
    K = K.tocsr()  # Convert to CSR for efficient solving

Linear solvers:
    from scipy.sparse.linalg import spsolve, gmres, cg, splu

    # Direct solver (good for N < 10000)
    x = spsolve(A, b)

    # Iterative solvers (for larger systems)
    x, info = gmres(A, b, tol=1e-10)
    x, info = cg(A, b, tol=1e-10)  # Only for SPD matrices

    # LU factorization (reusable for multiple RHS)
    lu = splu(A)
    x = lu.solve(b)

2.3 Optimization
----------------
Current: PETSc TAO, DONLP2
Replace with: scipy.optimize

    from scipy.optimize import minimize, Bounds

    # L-BFGS-B (matches TAO BLMVM)
    result = minimize(
        objective_function,
        x0=initial_guess,
        method='L-BFGS-B',
        bounds=Bounds(lb, ub),
        options={'maxiter': 100, 'disp': True}
    )

    # SLSQP (matches DONLP2 SQP with constraints)
    result = minimize(
        objective_function,
        x0=initial_guess,
        method='SLSQP',
        bounds=Bounds(lb, ub),
        constraints={'type': 'eq', 'fun': constraint_func}
    )

2.4 Spline Interpolation
------------------------
Current: GSL, PCHIP (custom), FITPACK, ALGLIB
Replace with: scipy.interpolate

    from scipy.interpolate import (
        CubicSpline,           # GSL cubic equivalent
        Akima1DInterpolator,   # GSL/ALGLIB akima
        PchipInterpolator,     # PCHIP (shape-preserving)
        UnivariateSpline,      # FITPACK equivalent
        make_interp_spline     # General B-splines
    )

    # Cubic spline with natural BC
    cs = CubicSpline(x, y, bc_type='natural')
    y_interp = cs(x_new)

    # PCHIP (monotonicity-preserving)
    pchip = PchipInterpolator(x, y)
    y_interp = pchip(x_new)

    # Akima
    akima = Akima1DInterpolator(x, y)
    y_interp = akima(x_new)

2.5 Bessel Functions
--------------------
Current: GSL gsl_sf_bessel_J1, gsl_sf_bessel_Y1
Replace with: scipy.special

    from scipy.special import jv, yv, j0, j1, y0, y1

    # First-order Bessel functions
    J1 = j1(x)   # or jv(1, x)
    Y1 = y1(x)   # or yv(1, x)
    J0 = j0(x)
    Y0 = y0(x)

2.6 Parallelism (Optional)
--------------------------
Current: PETSc + MPI
Replace with: Not needed for this problem size (~2000 DOF)

If parallelism is truly needed:
    # Option A: petsc4py (keeps PETSc dependency)
    from petsc4py import PETSc

    # Option B: mpi4py for distributed computing
    from mpi4py import MPI

    # Option C: multiprocessing for parallel function evals
    from multiprocessing import Pool

    # Option D: NumPy/SciPy vectorization (implicit parallelism via BLAS)
    # Modern NumPy uses multi-threaded BLAS automatically

================================================================================
3. CODE STRUCTURE FOR PYTHON PORT
================================================================================

Proposed file structure:

swirl2d_python/
├── __init__.py
├── mesh.py              # Mesh generation (Triangle wrapper)
├── fem.py               # FEM assembly (element matrices, BC)
├── solver.py            # Linear solver wrapper
├── optimizer.py         # Optimization interface
├── splines.py           # Spline interpolation
├── physics.py           # Boundary conditions, functional
├── io.py                # Tecplot, VTK output
├── main.py              # Main driver
└── utils.py             # Utilities

3.1 mesh.py
-----------
class Mesh:
    def __init__(self):
        self.vertices = None      # (N, 2) array
        self.triangles = None     # (M, 3) array
        self.markers = None       # (N,) boundary markers
        self.segments = None      # Boundary segments

    def generate_from_pslg(self, vertices, segments, markers, area=0.005):
        """Generate mesh using Triangle"""
        import triangle
        A = {'vertices': vertices, 'segments': segments,
             'segment_markers': markers}
        B = triangle.triangulate(A, f'pqa{area}')
        self.vertices = B['vertices']
        self.triangles = B['triangles']
        self.markers = B['vertex_markers'].flatten()

    def get_boundary_nodes(self, marker):
        """Get indices of nodes with given boundary marker"""
        return np.where(self.markers == marker)[0]

3.2 fem.py
----------
def element_matrix(coords, Acoef, Bcoef, Ccoef):
    """
    Compute element stiffness matrix and load vector
    for axisymmetric Bragg-Hawthorne equation

    coords: (3, 2) array of triangle vertex coordinates (z, r)
    Returns: ke (3,3), fe (3,)
    """
    # Quadrature points and weights (7-point rule)
    quad_w = np.array([0.225, 0.125939, 0.125939, 0.125939,
                       0.132394, 0.132394, 0.132394])
    quad_xi = np.array([1/3, 0.797427, 0.101287, 0.101287,
                        0.059716, 0.470142, 0.470142])
    quad_eta = np.array([1/3, 0.101287, 0.797427, 0.101287,
                         0.470142, 0.059716, 0.470142])

    z1, r1 = coords[0]
    z2, r2 = coords[1]
    z3, r3 = coords[2]

    area = (z2*r3 - z3*r2 + z3*r1 - z1*r3 + z1*r2 - z2*r1) / 2.0
    # ... (rest of element computation, matching C code)

    return ke, fe

def assemble_system(mesh, Acoef, Bcoef, Ccoef):
    """Assemble global stiffness matrix and load vector"""
    n = len(mesh.vertices)
    K = lil_matrix((n, n))
    F = np.zeros(n)

    for tri in mesh.triangles:
        coords = mesh.vertices[tri]
        ke, fe = element_matrix(coords, Acoef, Bcoef, Ccoef)
        for i in range(3):
            for j in range(3):
                K[tri[i], tri[j]] += ke[i, j]
            F[tri[i]] += fe[i]

    return K.tocsr(), F

def apply_dirichlet_bc(K, F, bc_nodes, bc_values):
    """Apply Dirichlet boundary conditions"""
    K = K.tolil()
    for node, val in zip(bc_nodes, bc_values):
        K[node, :] = 0
        K[node, node] = 1.0
        F[node] = val
    return K.tocsr(), F

3.3 solver.py
-------------
def solve_fem(K, F, method='direct'):
    """Solve the FEM system Kx = F"""
    if method == 'direct':
        return spsolve(K, F)
    elif method == 'gmres':
        x, info = gmres(K, F, tol=1e-10)
        if info != 0:
            raise RuntimeError(f"GMRES failed with info={info}")
        return x
    elif method == 'cg':
        x, info = cg(K, F, tol=1e-10)
        return x

3.4 physics.py
--------------
def compute_boundary_conditions(mesh, inlet_data, outlet_data, r_core):
    """
    Set up boundary conditions for all boundaries

    Returns: bc_nodes, bc_values arrays
    """
    bc_nodes = []
    bc_values = []

    # Inner boundary (stagnation streamline, psi=0)
    inner = mesh.get_boundary_nodes(MARKER_INNER)
    bc_nodes.extend(inner)
    bc_values.extend([0.0] * len(inner))

    # Outer boundary (wall streamline, psi=psi_wall)
    # ... etc

    return np.array(bc_nodes), np.array(bc_values)

def compute_functional(mesh, psi, inner_nodes, Acoef, Bcoef, Ccoef):
    """
    Compute flow force functional (objective function)
    """
    FFeval = 0.0
    for tri in mesh.triangles:
        # ... element integration (matches C functional())
        pass
    return FFeval

3.5 main.py
-----------
def objective(x, mesh_params, physics_params):
    """Objective function for optimizer"""
    # Generate inner boundary from control points
    inner_boundary = generate_inner_boundary(x, spline_type='pchip')

    # Generate mesh
    mesh = Mesh()
    mesh.generate_from_pslg(...)

    # Assemble and solve FEM
    K, F = assemble_system(mesh, ...)
    bc_nodes, bc_values = compute_boundary_conditions(...)
    K, F = apply_dirichlet_bc(K, F, bc_nodes, bc_values)
    psi = solve_fem(K, F)

    # Compute functional
    ff = compute_functional(mesh, psi, ...)

    return ff

def main():
    # Initial guess
    x0 = np.linspace(R_INLET, R_OUTLET, nparam + 2)[1:-1]

    # Bounds
    bounds = Bounds(R_INLET, R_OUTLET)

    # Optimize
    result = minimize(
        objective,
        x0,
        method='L-BFGS-B',
        bounds=bounds,
        options={'maxiter': 100, 'disp': True}
    )

    print(f"Optimal FFeval: {result.fun}")
    print(f"Optimal control points: {result.x}")

================================================================================
4. EFFORT ESTIMATE
================================================================================

Task                              Lines    Time (days)    Difficulty
----------------------------------------------------------------------
Mesh generation (mesh.py)          50         1           Easy
FEM assembly (fem.py)             200         3           Medium
Linear solver (solver.py)          30         0.5         Easy
Boundary conditions (physics.py)  150         2           Medium
Functional evaluation             100         1           Medium
Optimizer integration              50         0.5         Easy
I/O (Tecplot, VTK output)         100         1           Easy
Spline interpolation               50         0.5         Easy
Main driver                        50         0.5         Easy
Testing & validation               --         3           Medium
----------------------------------------------------------------------
TOTAL                            ~780        ~13         ~2 weeks

================================================================================
5. VALIDATION STRATEGY
================================================================================

5.1 Unit Tests
--------------
- Compare Triangle mesh output between C and Python wrappers
- Compare element stiffness matrices for identical triangles
- Compare Bessel function values (GSL vs scipy.special)
- Compare spline interpolation results

5.2 Integration Tests
---------------------
- Run C and Python versions with identical parameters
- Compare:
  - Mesh statistics (node count, element count)
  - Condition number of stiffness matrix
  - Solution PSI at specific nodes
  - Boundary flux values
  - Final FFeval (should match to ~1e-10)

5.3 Regression Tests
--------------------
- Verify optimal FFeval matches benchmark results in README
- Check that PCHIP still gives best results
- Timing comparison (Python may be 2-5x slower, acceptable for ~2000 DOF)

================================================================================
6. ADVANTAGES OF PYTHON PORT
================================================================================

1. Simplified Installation:
   - Current: PETSc (complex build), GSL, FITPACK, ALGLIB, Triangle, MPI
   - Python: pip install numpy scipy triangle matplotlib

2. Cross-Platform:
   - Works on Windows, Mac, Linux without compilation

3. Maintainability:
   - Single language (no C/C++/Fortran mix)
   - Easier debugging with Python tools

4. Extensibility:
   - Easy to add new spline types, optimizers, output formats
   - Integration with Jupyter notebooks for visualization

5. Dependencies:
   - NumPy, SciPy: stable, well-maintained, widely used
   - triangle: thin wrapper, same underlying C library

================================================================================
7. POTENTIAL ISSUES
================================================================================

1. Performance:
   - Python is slower than C for loops
   - Mitigation: Use NumPy vectorization, ~2000 DOF is small enough
   - Expected slowdown: 2-5x (still seconds, not minutes)

2. Numerical Differences:
   - Different BLAS/LAPACK implementations may give slightly different results
   - Mitigation: Use relative tolerance comparisons (1e-8)

3. Triangle Wrapper Limitations:
   - Python 'triangle' package may have fewer options than direct C API
   - Mitigation: Core functionality (quality meshing) is fully supported

4. No Native MPI:
   - If parallelism is needed, requires mpi4py or petsc4py
   - Mitigation: For 2000 DOF, sequential is fast enough (<1 second/solve)

================================================================================
8. MINIMAL VIABLE PYTHON VERSION
================================================================================

For quick prototyping, a minimal version needs only:

import numpy as np
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import spsolve
from scipy.interpolate import PchipInterpolator
from scipy.optimize import minimize, Bounds
from scipy.special import j1, y1
import triangle

This single-file implementation (~500 lines) could:
- Generate mesh with Triangle
- Assemble FEM system
- Solve with direct solver
- Optimize with L-BFGS-B
- Output Tecplot format

================================================================================
END OF DOCUMENT
================================================================================