IpoptSolver issue with mpirun

[Nlanglet]

Hello,
In order to carry out topology optimisation simulations in fluid mechanics using the density based method, i have to use an optimizer. I've decided to use IpoptSolver as in many tutorials I've found, but when I run the calculation with mpirun with more than 1 processor, ipopt stops at the first iteration. However, other tutorials or online code use this method, without problem.

Ipopt solver output :

This is Ipopt version 3.14.16, running with linear solver MUMPS 5.7.3.
Number of nonzeros in equality constraint Jacobian...:        0
Number of nonzeros in inequality constraint Jacobian.:    67876
Number of nonzeros in Lagrangian Hessian.............:        0
Process 0: Solving nonlinear variational problem.
  0 SNES Function norm 8.532842353045e+01
  1 SNES Function norm 3.440327740165e-03
  2 SNES Function norm 1.039981982021e-08
  Process 0: PETSc SNES solver converged in 2 iterations with convergence reason CONVERGED_FNORM_RELATIVE.
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  6.1117734e+01 0.00e+00 6.56e-02   0.0 0.00e+00    -  0.00e+00 0.00e+00   0

import os
import numpy as np
import matplotlib
matplotlib.use('Agg')  # Utilisation du backend Agg pour un environnement sans GUI
import matplotlib.pyplot as plt
from fenics import *
from xml.etree.ElementTree import Element, SubElement, tostring
from xml.dom.minidom import parseString
from dolfin_adjoint import * 
import pyadjoint
from pyadjoint.tape import no_annotations
import mpi4py
from fenics import XDMFFile

# Some flags for FEniCS
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["cpp_optimize_flags"] = "-O3 -ffast-math -march=native"
parameters['allow_extrapolation'] = True # Allow small numerical differences in the boundary definition.

# Quadrature degree in FEniCS  (sometimes, the "automatic" determination of the quadrature degree becomes excessively high, meaning that it should be manually reduced)
parameters['form_compiler']['quadrature_degree'] = 5 

parameters["std_out_all_processes"] = False
import time
start_time = time.time()

# Fluid flow setup
Probtype = "DoublePipe" # "DoublePipe" ou "PipeBend" ou "Diffuser2D"

rho_ = 1.0; mu_ = 1.0 # Density and dynamic viscosity
if Probtype == "DoublePipe":
    width_inlet_outlet = 1.0/6.0 # Inlet/outlet width

x_min = 0.0; x_max = 1.5 # x dimensions
y_min = 0.0; y_max = 1.0 # y dimensions
v_max_inlet = 1      # Inlet velocity
f_V = 1/3 # Volume fraction
xinit = f_V - 0.01
flow_regime = 'laminar' # Flow regime: 'laminar' or 'turbulent (Spalart-Allmaras)'
w_ud = 2000

# Topology optimization setup
k_max = 2.5*mu_/(0.01**2); k_min = 2.5*mu_/(100**2)
#k_max = 25000

def k(alpha):
    return k_max + (k_min-k_max)*alpha*(1+q)/(alpha+q)

obj = []
ité = []

# Output folders
output_folder = 'Thèse/FEniCS/Doublepipe/Re1/output'
problem_folder = "%s/foam_problem" %(output_folder)
if (MPI.comm_world.Get_size() == 1) or (MPI.comm_world.Get_rank() == 0):
	if not os.path.exists(output_folder):
		os.makedirs(output_folder) # Create the output folder if it still does not exist


# Create the 2D mesh and plot it
N_mesh = 150
delta_x = x_max - x_min; delta_y = y_max - y_min
mesh = RectangleMesh(Point(x_min, y_min), Point(x_max, y_max), int(N_mesh*delta_x/delta_y), N_mesh, diagonal = "crossed")
File('%s/mesh.pvd' %(output_folder)) << mesh

# Function spaces -> MINI element (2D)
V1_element =  FiniteElement('Lagrange', mesh.ufl_cell(), 1)
B_element = FiniteElement('Bubble', mesh.ufl_cell(), 3)
V_element = VectorElement(NodalEnrichedElement(V1_element, B_element)) # Velocity
P_element = FiniteElement('Lagrange', mesh.ufl_cell(), 1) # Pressure
if flow_regime == 'laminar':
	U_element = MixedElement([V_element, P_element])
U = FunctionSpace(mesh, U_element) # Mixed function space
A_element = FiniteElement('Lagrange', mesh.ufl_cell(), 1)
A = FunctionSpace(mesh, A_element) # Design variable function space (nodal)


# Prepare the boundary definition

if Probtype == "DoublePipe":
    class Inlet2(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary and (x[0] == x_min and ((y_min + 3.0/4*delta_y - width_inlet_outlet/2) < x[1] < (y_min + 3.0/4*delta_y + width_inlet_outlet/2))) 

    class Inlet1(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary and (x[0] == x_min and ((y_min + 1.0/4*delta_y - width_inlet_outlet/2) < x[1] < (y_min + 1.0/4*delta_y + width_inlet_outlet/2)))
    
    class Outlet1(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary and (x[0] == x_max and ((y_min + 1.0/4*delta_y - width_inlet_outlet/2) < x[1] < (y_min + 1.0/4*delta_y + width_inlet_outlet/2)))
    
    class Outlet2(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary and (x[0] == x_max and ((y_min + 3.0/4*delta_y - width_inlet_outlet/2) < x[1] < (y_min + 3.0/4*delta_y + width_inlet_outlet/2))) 

    class Walls(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary # * It will be set before the other boundaries  

    marker_numbers = {'unset' : 0, 'wall' : 1, 'inlet1' : 2, 'inlet2' : 3, 'outlet1' : 4, 'outlet2' : 5}
    boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
    boundary_markers.set_all(marker_numbers['unset'])
    Walls().mark(boundary_markers, marker_numbers['wall'])
    Inlet1().mark(boundary_markers, marker_numbers['inlet1'])
    Inlet2().mark(boundary_markers, marker_numbers['inlet2'])
    Outlet1().mark(boundary_markers, marker_numbers['outlet1'])
    Outlet2().mark(boundary_markers, marker_numbers['outlet2'])
    File("%s/markers.pvd" %(output_folder)) << boundary_markers
    


if Probtype =="DoublePipe":
    class InletVelocity1(UserExpression):
        def eval(self, values, x):
            for i in range(len(values)):
                values[i] = 0.0 # Initialize all values with zeros
            if (x[0] == x_min or x[0] == x_max) and (1.0/4*delta_y - width_inlet_outlet/2) < x[1] < (1.0/4*delta_y + width_inlet_outlet/2):
                y_local = x[1] - 1.0/4*delta_y; values[0] = v_max_inlet*(1 - (2*y_local/width_inlet_outlet)**2)
        def value_shape(self):
            return (2,)

    class InletVelocity2(UserExpression):
        def eval(self, values, x):
            for i in range(len(values)):
                values[i] = 0.0 # Initialize all values with zeros
            if (x[0] == x_min or x[0] == x_max) and (3.0/4*delta_y - width_inlet_outlet/2) < x[1] < (3.0/4*delta_y + width_inlet_outlet/2):
                y_local = x[1] - 3.0/4*delta_y; values[0] = v_max_inlet*(1 - (2*y_local/width_inlet_outlet)**2)
        def value_shape(self):
            return (2,)

    inlet_velocity_expression1 = InletVelocity1(element = V_element)
    inlet_velocity_expression2 = InletVelocity2(element = V_element)
    wall_velocity_value = Constant((0,0))


# Function to solve the forward problem
def solve_forward_problem(alpha):
	# Set the state vector and test functions
    u = Function(U); u.rename("StateVariable", "StateVariable")
    u_split = split(u); v = u_split[0]; p = u_split[1]
    w = TestFunction(U)
    w_split = split(w); w_v = w_split[0]; w_p = w_split[1]

	# Weak form of the pressure-velocity formulation
    I_ = as_tensor(np.eye(2)); T = -p*I_ + (mu_)*(grad(v) + grad(v).T); v_mat = v
    F = div(v) * w_p*dx + inner( grad(w_v), T )*dx + rho_ * inner( dot(grad(v), v), w_v )*dx + inner(k(alpha) * v_mat, w_v)*dx

	# Boundary conditions
    
    if Probtype == "DoublePipe":
        bcs = [DirichletBC(U.sub(0), wall_velocity_value, boundary_markers, marker_numbers['wall']),
            DirichletBC(U.sub(0), inlet_velocity_expression1, boundary_markers, marker_numbers['inlet1']),
            DirichletBC(U.sub(0), inlet_velocity_expression2, boundary_markers, marker_numbers['inlet2']),
            DirichletBC(U.sub(0), inlet_velocity_expression1, boundary_markers, marker_numbers['outlet1']),
            DirichletBC(U.sub(0), inlet_velocity_expression2, boundary_markers, marker_numbers['outlet2'])]

	# Solve the simulation
    dF = derivative(F, u); problem = NonlinearVariationalProblem(F, u, bcs, dF)
    solver = NonlinearVariationalSolver(problem)
    params = {'nonlinear_solver': 'snes',
           'snes_solver':
            {'linear_solver'     : 'mumps',
                'preconditioner'  : 'hypre_euclid' }}
    solver.parameters.update(params)
    solver.solve()
    return u

global current_iteration
current_iteration = 0
# Initial setup for topology optimization
alpha = interpolate(Constant(xinit), A) # Initial guess for the design variable
for idx, q in enumerate([0.01, 0.1]):
    print("Valeur de q:",q)
    set_working_tape(Tape())           # Clear all annotations and restart the adjoint model
    file_prefix = f"q{idx:02d}_"
    # Solve the simulation
    u = solve_forward_problem(alpha)

    alpha_pvd_file = File(f"{output_folder}/{file_prefix}alpha_iterations.pvd")
    alpha_viz = Function(A, name="AlphaVisualisation")
    dj_pvd_file = File(f"{output_folder}/{file_prefix}dj_iterations.pvd")
    dj_viz = Function(A, name="dJVisualisation")

    # Callback during topology optimization

    def derivative_cb_post(j, dj, current_alpha):
        global current_iteration
        print("\n [Iteration: %d] J = %1.7e, q = %1.2f\n" %(current_iteration, j, q))
        obj.append(j)
        ité.append(current_iteration)
        current_iteration += 1
        # Save for visualization
        alpha_viz.assign(current_alpha); alpha_pvd_file << (alpha_viz,current_iteration)
        dj_viz.assign(dj); dj_pvd_file << dj_viz

    # Objective function
    u_split = split(u); v = u_split[0]
    J = assemble((1/2.*(mu_)*inner(grad(v) + grad(v).T, grad(v) + grad(v).T) + inner(k(alpha) * v, v))*dx)
    #J = assemble(w_ud*dot(v - u_target, v - u_target) * dx(6) + inner(k(alpha) * v, v)*dx)

    print(" Current objective function value: %1.7e" %(J))

    # Set the topology optimization problem and solver
    alpha_C = Control(alpha)
    Jhat = ReducedFunctional(J, alpha_C, derivative_cb_post = derivative_cb_post)

    problem_min = MinimizationProblem(Jhat, bounds = (0.0, 1.0), constraints = [UFLInequalityConstraint((f_V - alpha)*dx, alpha_C)])
    problem_min.options["print_options_documentation"] = "yes"
    max_iterations = 50 if idx == 0 else 100
    tol = 1e-3 if idx == 0 else 1e-8
    solver_opt = IPOPTSolver(problem_min, parameters = {'maximum_iterations': max_iterations, 'tol':tol, 'print_level' : 5})

    # Perform topology optimization
    alpha.assign(alpha_opt); alpha_viz.assign(alpha)
    alpha_pvd_file << alpha_viz

# Plot a simulation
u = solve_forward_problem(alpha)
u_split_deepcopy = u.split(deepcopy = True)
v_plot = u_split_deepcopy[0]
p_plot = u_split_deepcopy[1]
File("%s/simulation_final_v.pvd" %(output_folder)) << v_plot
File("%s/simulation_final_p.pvd" %(output_folder)) << p_plot

objective_image_path = os.path.join(output_folder, "objective_evolution.png")
plt.figure(figsize=(6, 5), tight_layout=True)
plt.plot(ité,obj)
plt.xlabel('Itération',fontsize=20)
plt.ylabel('Objective',fontsize=20)
plt.grid()
plt.savefig(objective_image_path, dpi=300)
plt.close()

def merge_pvd_files(output_filename, pvd_files):
    """
    Fusionne plusieurs fichiers PVD en un seul, avec des chemins corrects.
    """
    root = Element("VTKFile", type="Collection", version="0.1", byte_order="LittleEndian")
    collection = SubElement(root, "Collection")

    timestep = 0
    for pvd_file in pvd_files:
        if not os.path.exists(pvd_file):
            print(f"Fichier manquant : {pvd_file}")
            continue
        
        # Trouver le répertoire du fichier PVD pour résoudre les chemins
        pvd_dir = os.path.dirname(pvd_file)

        # Lire le contenu du fichier PVD
        with open(pvd_file, "r") as f:
            lines = f.readlines()
            for line in lines:
                if "DataSet" in line:  # Rechercher les entrées DataSet
                    # Extraire le chemin du fichier VTU
                    relative_file = line.split('file="')[1].split('"')[0].strip()
                    full_path = os.path.join(pvd_dir, relative_file)
                    if not os.path.exists(full_path):
                        print(f"Fichier .vtu manquant : {full_path}")
                        continue
                    
                    # Ajouter au fichier fusionné avec le nouveau timestep
                    dataset = SubElement(collection, "DataSet", {
                        "timestep": str(timestep),
                        "group": "",
                        "part": "0",
                        "file": relative_file,
                    })
                    timestep += 1

    # Sauvegarder le fichier fusionné
    xml_string = tostring(root)
    pretty_xml = parseString(xml_string).toprettyxml()
    with open(output_filename, "w") as out_f:
        out_f.write(pretty_xml)
# Liste des fichiers générés à fusionner
pvd_files = []
for idx, q in enumerate([0.01, 0.1]):
    file_prefix = f"q{idx:02d}_"
    pvd_files.append(f"{output_folder}/{file_prefix}alpha_iterations.pvd")
output_filename = f"{output_folder}/fused_alpha_iterations.pvd"

merge_pvd_files(output_filename, pvd_files)

end_time = time.time()
elapsed_time = end_time - start_time
print(f"Temps écoulé : {elapsed_time:.2f} secondes")

Thanks !