Confused on parallel running on Fipy

[bladesaber]

Thinks for the great work on Fipy. I'm using Fipy on solving stoke equation based on SIMPLE method in CFD. And after using 'mpirun -np 4 python test4.py', only part of mesh has result and the rest are all full zero which is wrong. I'm conflused and what is the problem?

Result is shown below:

and here is my code:

import os
import shutil
import numpy as np
from mpi4py import MPI

from fipy import Gmsh2D
from fipy import CellVariable, FaceVariable
from fipy.variables.faceGradVariable import _FaceGradVariable
from fipy import DiffusionTerm, ConvectionTerm
from fipy.tools import numerix
from fipy import Viewer


def inflow_velocity_exp(
    pos: np.ndarray, 
    velocity: np.ndarray, 
    tdim: int, 
    face_center: np.ndarray, 
    face_norm_vector: np.ndarray,
    face_radius: float
):
    bry_values = np.zeros_like(pos.T)
        
    face_norm_vector = face_norm_vector / np.linalg.norm(face_norm_vector)
    face_norm_velocity = np.dot(velocity, face_norm_vector) * face_norm_vector
    bry_values[:, :tdim] = face_norm_velocity
        
    pos = pos.T[:, :tdim]
    dist2center = np.linalg.norm(pos - face_center, ord=2, axis=1)
    dist2plane = np.sum((pos - face_center) * face_norm_vector, axis=1)
        
    weight = np.maximum((face_radius - dist2center) / face_radius, 0.0)
    weight[np.abs(dist2plane) > 1e-6] = 0.0
    weight = -np.power(weight - 1.0, 2.0) + 1.0

    bry_values = bry_values * weight.reshape((-1, 1))
    bry_values = bry_values.T[:tdim, :]
        
    return bry_values


comm = MPI.COMM_WORLD
rank = comm.Get_rank()
size = comm.Get_size()

base_dir = "/home/admin123456/Desktop/work/gmsh_test/t38"
mesh = Gmsh2D(os.path.join(base_dir, "t38_model_partition_4.msh"))
viscosity=1.0 / 400.0

pressure = CellVariable(mesh=mesh, name="pressure")
xVelocity = CellVariable(mesh=mesh, name="X velocity")
yVelocity = CellVariable(mesh=mesh, name="Y velocity")
face_velocity = FaceVariable(mesh=mesh, rank=1)

xVelocityEq = DiffusionTerm(coeff=viscosity) - pressure.grad.dot([1., 0.])
yVelocityEq = DiffusionTerm(coeff=viscosity) - pressure.grad.dot([0., 1.])
        
# ------ define pressure correction equation ------
ap = CellVariable(mesh=mesh, value=1.)
coeff = 1./ ap.arithmeticFaceValue * mesh._faceAreas * mesh._cellDistances
pressureCorrectionEq = DiffusionTerm(coeff=coeff) - face_velocity.divergence
        
pressureCorrection = CellVariable(mesh=mesh)
pressureCorrectionSolver = pressureCorrectionEq.getDefaultSolver(tolerance=1e-10, iterations=2000)
        
# ------ define other variables ------
volume = CellVariable(mesh=mesh, value=mesh.cellVolumes, name="Volume")
control_volume = volume.arithmeticFaceValue

# ------ define boundary ------
velocity_values = inflow_velocity_exp(
    pos=mesh.faceCenters.value,
    velocity=np.array([1.0, 0.0]),
    tdim=2,
    face_center=[-0.05, 0.2],  # type: ignore
    face_norm_vector=[1.0, 0.0],  # type: ignore
    face_radius=0.1
)
xVelocity.constrain(velocity_values[0, :], where=mesh.physicalFaceMap == 22)
yVelocity.constrain(velocity_values[1, :], where=mesh.physicalFaceMap == 22)

pressure.constrain(0.0, where=mesh.physicalFaceMap == 23)
pressureCorrection.constrain(0., where=mesh.physicalFaceMap == 23)

xVelocity.constrain(0., where=mesh.physicalFaceMap == 24)
yVelocity.constrain(0., where=mesh.physicalFaceMap == 24)

# ------ define save directory ------
save_dir = os.path.join(base_dir, "cache")
if rank == 0:
    if save_dir is not None:
        if os.path.exists(save_dir):
            shutil.rmtree(save_dir)
        os.makedirs(save_dir)
        
# ------ wait for all processes to finish initialization ------
comm.Barrier()

velocityRelaxation = 0.5
pressureRelaxation = 0.8
viewer = Viewer(
    vars=(face_velocity),
    xmin=0., xmax=1., ymin=0., ymax=1., colorbar="vertical", scale=5
)
        
for iter_idx in range(60):
    # ------ solve the momentum equation ------
    xVelocityEq.cacheMatrix()
    x_residual = xVelocityEq.sweep(var=xVelocity, underRelaxation=velocityRelaxation)
            
    yVelocityEq.cacheMatrix()
    y_residual = yVelocityEq.sweep(var=yVelocity, underRelaxation=velocityRelaxation)
    # ------ solve the momentum equation ------

    # ------ update the ap coefficient from the matrix diagonal ------
    xmat = xVelocityEq.matrix
    ap_x = -numerix.asarray(xmat.takeDiagonal())

    ymat = yVelocityEq.matrix
    ap_y = -numerix.asarray(ymat.takeDiagonal())

    ap[:] = (ap_x + ap_y) * 0.5
    # ------ update the ap coefficient from the matrix diagonal ------
            
    # ------ interpolate the pressure gradient to the face based on linear pressure gradient difference
    presgrad = pressure.grad  # cell pressure gradient
    facepresgrad = _FaceGradVariable(pressure)  # face pressure gradient

    face_velocity[0] = (
        xVelocity.arithmeticFaceValue 
        + control_volume / ap.arithmeticFaceValue * (presgrad[0].arithmeticFaceValue - facepresgrad[0])
        )
    face_velocity[1] = (
        yVelocity.arithmeticFaceValue 
        + control_volume / ap.arithmeticFaceValue * (presgrad[1].arithmeticFaceValue - facepresgrad[1])
    )

    # ------ apply_boundary_to_correction_velocity ------
    velocity_values = inflow_velocity_exp(
        pos=mesh.faceCenters.value,
        velocity=np.array([1.0, 0.0]),
        tdim=2,
        face_center=[-0.05, 0.2],  # type: ignore
        face_norm_vector=[1.0, 0.0],  # type: ignore
        face_radius=0.1
    )
    face_mask = mesh.physicalFaceMap == 22
    face_velocity[0, face_mask] = velocity_values[0, face_mask]
    face_velocity[1, face_mask] = velocity_values[1, face_mask]
    
    face_mask = mesh.physicalFaceMap == 24
    face_velocity[0, face_mask] = 0.
    face_velocity[1, face_mask] = 0.
    # ------ apply_boundary_to_correction_velocity ------

    # ------ solve the pressure correction equation ------
    pressureCorrectionEq.cacheRHSvector()
    ## left bottom point must remain at pressure 0, so no correction
    p_residual = pressureCorrectionEq.sweep(var=pressureCorrection, solver=pressureCorrectionSolver)
    continuity_residual = pressureCorrectionEq.RHSvector

    # ------ update the pressure and velocity ------
    pressure.setValue(pressure + pressureRelaxation * pressureCorrection)
    xVelocity.setValue(xVelocity - pressureCorrection.grad[0] / ap * mesh.cellVolumes)
    yVelocity.setValue(yVelocity - pressureCorrection.grad[1] / ap * mesh.cellVolumes)
            
    # ------ wait for all processes to finish calculating ------
    comm.Barrier()
            
    # ------ reduce the residual to the master process ------
    global_x_residual: float = comm.reduce(x_residual, op=MPI.MAX, root=0)  # type: ignore
    global_y_residual: float = comm.reduce(y_residual, op=MPI.MAX, root=0)  # type: ignore
    global_p_residual: float = comm.reduce(p_residual, op=MPI.MAX, root=0)  # type: ignore
    global_continuity_residual: float = comm.reduce(np.max(np.abs(continuity_residual)), op=MPI.MAX, root=0)  # type: ignore
            
    if rank == 0:
        log = f"[Log: {iter_idx}]: x residual:{global_x_residual:.8f}, y residual:{global_y_residual:.8f}, "
        log += f" p residual:{global_p_residual:.8f}, continuity_residual:{global_continuity_residual:.8f}"
        print(log, flush=True)
    
    # record(save_file=os.path.join(save_dir, f"iter_{iter_idx}.vtu"))
    viewer.plot()

Looking forward to response.
Best regards.

[bladesaber]

And mesh file are uploaded:

t38_model_partition_4.zip

[bladesaber]

I find the point, when using gmsh tool, command line should be gmsh *.msh -format msh2 -part 4 -part_ghosts -o *.msh -save , option part_ghosts must be enable.

[guyer]

Oh, I see. FiPy accounts for -part_ghosts with .geo files, but not .msh. I'm in the middle of overhauling our Gmsh interaction, so I'll have to be sure to account for this.