Issues regarding topology optimisation xdmf file results errors in parallel computation using MPI

dolfinx
1

Hi~ :slight_smile:
The reason for posting a question is when topology optimisation was executed with mpirun -n 2 ~~ in the previous code posted(link1), the shape of the entire domain did not appear in ParaView when the mesh result was extracted in half as an xdmf file.

Topology optimisation was performed using Euclidean distance filter & Optimality creteria.
In general, if the file was executed in a single process, the following shape was obtained.

capture_single

Even when the same code is executed with mpirun -n 2 ~~, the same shape (I think there will be an afterimage of the ghost node depending on the number of cores) should appear, but only half of it was displayed in the entire domain.

capture_2core

I left a question because I wondered if the problem was caused by a code error (an error that occurs the code is not parallelized) or an error that occurs while creating an xdmf file through write_mesh & write_function.

Is there any way to fix that error?

1 Like
2

Without a minimal code example illustrating what you have implemented, it is really hard to pinpoint where things have gone wrong. In the post you are linking to, the code is fragmented over many posts and many blocks, making it guesswork to find the error.

1 Like
3

Sorry. @dokken. Here is my code that I ran:

# Set python packages
import ufl
import numpy as np, sklearn.metrics.pairwise as sp
from dolfinx import fem, mesh, io
from mpi4py import MPI
from petsc4py import PETSc

# Dolfinx projection function
def project_func(dfx_func, func_space):
    trial_func = ufl.TrialFunction(func_space)
    test_func = ufl.TestFunction(func_space)
    a = trial_func * test_func * ufl.dx
    l = dfx_func * test_func * ufl.dx
    project_prob = fem.petsc.LinearProblem(a, l, [])
    result_sol = project_prob.solve()
    return result_sol

# Declare the input variables
nelx = 60
nely = 20
volfrac = 0.5
penal = 3.0
rmin = 3.0

# Put MPI commands
comm = MPI.COMM_WORLD

# Setting preliminaries
sigma = lambda _u: 2.0 * mu * ufl.sym(ufl.grad(_u)) + lmd * ufl.tr(ufl.sym(ufl.grad(_u))) * ufl.Identity(len(_u))
psi = lambda _u: lmd / 2 * (ufl.tr(ufl.sym(ufl.grad(_u))) ** 2) + mu * ufl.tr(ufl.sym(ufl.grad(_u)) * ufl.sym(ufl.grad(_u)))
mu, lmd = PETSc.ScalarType(0.4), PETSc.ScalarType(0.6)

# Prepare finite element analysis
msh = mesh.create_rectangle(comm, ((0.0, 0.0), (nelx, nely)), (nelx, nely), cell_type=mesh.CellType.triangle, diagonal=mesh.DiagonalType.right_left)
with io.XDMFFile(comm, "output/density.xdmf", "w") as file:
    file.write_mesh(msh)
U = fem.VectorFunctionSpace(msh, ("CG", 1))
D = fem.FunctionSpace(msh, ("DG", 0))
u, v = ufl.TrialFunction(U), ufl.TestFunction(U)
u_sol, density_old, density, sensitivity = fem.Function(U), fem.Function(D), fem.Function(D), fem.Function(D)
density.vector.array = volfrac

# Put the dimension of domain
t_dim = msh.topology.dim
f_dim = t_dim - 1
num_cells = msh.topology.index_map(t_dim).size_local

# Define support
def left_clamp(x):
    return np.isclose(x[0], 0.0)
bc_facets = mesh.locate_entities_boundary(msh, f_dim, left_clamp)
u_zero = np.array([0.0, 0.0], dtype=PETSc.ScalarType)
bc_l = fem.dirichletbc(u_zero, fem.locate_dofs_topological(U, f_dim, bc_facets), U)
bcs = [bc_l]

# Define external load
load_points = [(1, lambda x: np.logical_and(x[0]==nelx, x[1]<=2))]
facet_indices, facet_markers = [], []
for (marker, locator) in load_points:
    facets = mesh.locate_entities(msh, f_dim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full(len(facets), marker))
facet_indices = np.array(np.hstack(facet_indices), dtype=np.int32)
facet_markers = np.array(np.hstack(facet_markers), dtype=np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = mesh.meshtags(msh, f_dim, facet_indices[sorted_facets], facet_markers[sorted_facets])

ds = ufl.Measure("ds", domain=msh, subdomain_data=facet_tag)
f = ufl.dot(v, fem.Constant(msh, PETSc.ScalarType((0.0, -1.0)))) * ds(1)

# Set up the variational problem and solver (using SIMP)
k = ufl.inner(density**penal * sigma(u), ufl.grad(v)) * ufl.dx
problem = fem.petsc.LinearProblem(k, f, bcs)

# Prepare euclidean distance filter
midpoint = mesh.compute_midpoints(msh, t_dim, range(num_cells))
distance_mat = rmin - sp.euclidean_distances(midpoint, midpoint)
distance_mat[distance_mat < 0] = 0
distance_sum = distance_mat.sum(1)

# Topology optimisation using for-loop
loop, change = 0, 1
while change > 0.01 and loop < 2000:
    loop += 1
    density_old.vector.array = density.vector.array

    # FEA analysis
    u_sol = problem.solve()

    # Objective function
    objective = project_func(density**penal * psi(u_sol), D)

    # Sensitivity analysis using euclidean distance filter
    sensitivity.vector.array = -penal * (density.vector.array)**(penal-1) * project_func(psi(u_sol), D).vector.array
    sensitivity.vector.array = np.divide(distance_mat @ np.multiply(density.vector.array, sensitivity.vector.array), np.multiply(density.vector.array, distance_sum))

    # Domain update using optimality criteria method
    l1, l2, move = 0, 1e5, 0.2
    while l2 - l1 > 1e-4:
        l_mid = 0.5 * (l2 + l1)
        density_new = np.maximum(0.001, np.maximum(density.vector.array - move, np.minimum(1.0, np.minimum(density.vector.array + move, density.vector.array * np.sqrt(-sensitivity.vector.array/l_mid)))))
        l1, l2 = (l_mid, l2) if sum(density_new) - volfrac*num_cells > 0 else (l1, l_mid)
    density.vector.array = density_new

    # Print results
    change = np.linalg.norm(density.vector.array - density_old.vector.array, np.inf)
    change = comm.bcast(change, root=0)
    print("it.: {0} , obj.: {1:.3f} Vol.: {2:.3f}, ch.: {3:.3f}".format(loop, sum(objective.vector.array), sum(density.vector.array)/num_cells, change))
    file.write_function(density, loop)

I got the above result with that code.

1 Like
4

Your code is rather long and hard to parse. However, i note that you are using dolfinx.fem.Function.vector.array for your updates in the code above. You have to be very careful when doing that in parallel, at it only works on local data (not including ghosts).

Preferrably you should use dolfinx.fem.Function.x.array
For such computations.

As long as your code is this lengthy, i dont have time to go into further details.

1 Like
5

To @dokken. Sorry for checking your reply late.
As suggested, when changing dolfinx.fem.Function.vector.array to dolfinx.fem.Function.x.array, a problem occurred between ghost nodes.

When the total domain size is 2 by 2, density.x.array[:] is

[0.5 0.5 0.5 0.5 0.5 0.5]

(Unlike dolfin, dolfinx thought that the value for the ghost node was also output.)
And the error occured when I run the code (mpirun -n 2 ~~) was:

 sensitivity = np.divide(distance_mat @ np.multiply(density.x.array[:], sensitivity), np.multiply(density.x.array[:], distance_sum))
ValueError: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 6 is different from 4)
    sensitivity = np.divide(distance_mat @ np.multiply(density.x.array[:], sensitivity), np.multiply(density.x.array[:], distance_sum))
ValueError: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 6 is different from 4)
--------------------------------------------------------------------------
Primary job  terminated normally, but 1 process returned
a non-zero exit code. Per user-direction, the job has been aborted.
--------------------------------------------------------------------------
--------------------------------------------------------------------------
mpirun detected that one or more processes exited with non-zero status, thus causing
the job to be terminated. The first process to do so was:

  Process name: [[56127,1],0]
  Exit code:    1
--------------------------------------------------------------------------

As the cause of the error, the code I used to use the Euclidean distance filter

In the code, it was judged that the midpoint array was derived without considering the ghost node.
Like:

midpoint on process 1: 
[[0.33333333 1.66666667 0.        ]
 [0.66666667 1.33333333 0.        ]
 [1.33333333 1.33333333 0.        ]
 [1.66666667 1.66666667 0.        ]]

midpoint on process 0: 
[[0.33333333 0.33333333 0.        ]
 [0.66666667 0.66666667 0.        ]
 [1.33333333 0.66666667 0.        ]
 [1.66666667 0.33333333 0.        ]]

Then, how should the midpoint array derivation considering the ghost node proceed in a parallel environment?

6

This is because you are using:

instead of
num_cells = mesh.topology.index_map(t_dim).size_local + mesh.topology.index_map(t_dim).num_ghosts.

7

Thank you for reply @dokken! :slight_smile:
The problem I was worried about has been solved!

However, additional concerns arose.
After modifying the code, when the file was executed, the problem occurred that the entire result value obtained from each local process was derived the same on first iteration.
Like:

Initial density on process 0: 
[0.5 0.5 0.5 0.5 0.5 0.5]

Initial density on process 1: 
[0.5 0.5 0.5 0.5 0.5 0.5]

u_sol on process 0: 
[ 0.00000000e+00  0.00000000e+00 -1.51458886e+01 -2.69893899e+01
  0.00000000e+00  0.00000000e+00 -2.19098143e+01 -6.77188329e+01
 -4.63518113e-15 -2.81564987e+01 -1.06581410e-14 -6.44827586e+01
  1.51458886e+01 -2.69893899e+01]

u_sol on process 1: 
[ 0.00000000e+00  0.00000000e+00  1.51458886e+01 -2.69893899e+01
 -4.63518113e-15 -2.81564987e+01 -1.06581410e-14 -6.44827586e+01
  2.19098143e+01 -6.77188329e+01  0.00000000e+00  0.00000000e+00
 -1.51458886e+01 -2.69893899e+01]

Objective on process 0: 
[38.28299907  4.35108686 11.33439076 12.13231906  4.35108686 11.33439076]

Objective on process 1: 
[38.28299907  4.35108686 11.33439076 12.13231906  4.35108686 11.33439076]

Sensitivity on process 0: 
[-229.69799443  -26.10652119  -68.00634459  -72.79391433  -26.10652119
  -68.00634459]

Sensitivity on process 1: 
[-229.69799443  -26.10652119  -68.00634459  -72.79391433  -26.10652119
  -68.00634459]

The reason why Objective and Sensitivity are the same in the above result is thought to be an interpolation problem. Is there a method for U1 = fem.VectorFunctionSpace(msh, ("CG", 1)) to D0 = fem.FunctionSpace(msh, ("DG", 0)) projection (interpolation) for each parallel subdomain?

8

Please make a minimal reproducible example (as little code as possible) that gives you a similar result.
It is very hard to parse what you are doing with the amount of code you have supplied above, which does not have similar print statements

9

Thank you for reply @dokken!
The code for the result obtained above is as follows(Reduced code as much as possible)

# Set python packages
import ufl
import numpy as np
from dolfinx import fem, mesh
from mpi4py import MPI
from petsc4py import PETSc

# Dolfinx projection function
def project_func(dfx_func, func_space):
    trial_func = ufl.TrialFunction(func_space)
    test_func = ufl.TestFunction(func_space)
    a = trial_func * test_func * ufl.dx
    l = dfx_func * test_func * ufl.dx
    project_prob = fem.petsc.LinearProblem(a, l, [])
    result_sol = project_prob.solve()
    return result_sol

# Declare the input variables
nelx = 2
nely = 2
volfrac = 0.5
penal = 3.0

# Put MPI commands
comm = MPI.COMM_WORLD

# Setting preliminaries
sigma = lambda _u: 2.0 * mu * ufl.sym(ufl.grad(_u)) + lmd * ufl.tr(ufl.sym(ufl.grad(_u))) * ufl.Identity(len(_u))
psi = lambda _u: lmd / 2 * (ufl.tr(ufl.sym(ufl.grad(_u))) ** 2) + mu * ufl.tr(ufl.sym(ufl.grad(_u)) * ufl.sym(ufl.grad(_u)))
mu, lmd = PETSc.ScalarType(0.4), PETSc.ScalarType(0.6)

# Set mesh
msh = mesh.create_rectangle(comm, ((0.0, 0.0), (nelx, nely)), (nelx, nely), cell_type=mesh.CellType.triangle, diagonal=mesh.DiagonalType.right_left)
U = fem.VectorFunctionSpace(msh, ("CG", 1))
D = fem.FunctionSpace(msh, ("DG", 0))
u, v = ufl.TrialFunction(U), ufl.TestFunction(U)
u_sol, density = fem.Function(U), fem.Function(D)
density.x.array[:] = volfrac
print(f"Initial density on process {MPI.COMM_WORLD.rank}: \n{density.x.array[:]}\n")

# Put the dimension of domain
t_dim = msh.topology.dim
f_dim = t_dim - 1
num_cells = msh.topology.index_map(t_dim).size_local

# Define support
def left_clamp(x):
    return np.isclose(x[0], 0.0)
bc_facets = mesh.locate_entities_boundary(msh, f_dim, left_clamp)
u_zero = np.array([0.0, 0.0], dtype=PETSc.ScalarType)
bc_l = fem.dirichletbc(u_zero, fem.locate_dofs_topological(U, f_dim, bc_facets), U)
bcs = [bc_l]

# Define external load
load_points = [(1, lambda x: np.logical_and(x[0]==nelx, x[1]<=2))]
facet_indices, facet_markers = [], []
for (marker, locator) in load_points:
    facets = mesh.locate_entities(msh, f_dim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full(len(facets), marker))
facet_indices = np.array(np.hstack(facet_indices), dtype=np.int32)
facet_markers = np.array(np.hstack(facet_markers), dtype=np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = mesh.meshtags(msh, f_dim, facet_indices[sorted_facets], facet_markers[sorted_facets])

ds = ufl.Measure("ds", domain=msh, subdomain_data=facet_tag)
f = ufl.dot(v, fem.Constant(msh, PETSc.ScalarType((0.0, -1.0)))) * ds(1)

# Set up the variational problem and solver (using SIMP)
k = ufl.inner(density**penal * sigma(u), ufl.grad(v)) * ufl.dx
problem = fem.petsc.LinearProblem(k, f, bcs)

# solve PDE
u_sol = problem.solve()
print(f"u_sol on process {MPI.COMM_WORLD.rank}: \n{u_sol.x.array[:]}\n")

# Objective & sensitivity
objective = project_func(density**penal * psi(u_sol), D).x.array[:]
print(f"Objective on process {MPI.COMM_WORLD.rank}: \n{objective}\n")
sensitivity = -penal * (density.x.array[:])**(penal-1) * project_func(psi(u_sol), D).x.array[:]
print(f"Sensitivity on process {MPI.COMM_WORLD.rank}: \n{sensitivity}\n")

To comparing the result with the code built through dolfin in the past.
MWE is

# Set python packages
from dolfin import *
from mpi4py import MPI

# Declare the input variables
nelx = 2
nely = 2
volfrac = 0.5
penal = 3.0

# Setting preliminaries
sigma = lambda _u: 2.0 * mu * sym(grad(_u)) + lmbda * tr(sym(grad(_u))) * Identity(len(_u))
psi = lambda _u: lmbda / 2 * (tr(sym(grad(_u))) ** 2) + mu * tr(sym(grad(_u)) * sym(grad(_u)))
mu, lmbda = Constant(0.4), Constant(0.6)

# Set mesh
msh = RectangleMesh(Point(0, 0), Point(nelx, nely), nelx, nely, "right/left")
U = VectorFunctionSpace(msh, "P", 1)
D = FunctionSpace(msh, "DG", 0)
u, v = TrialFunction(U), TestFunction(U)
u_sol, density = Function(U), Function(D)
density.vector()[:] = volfrac
print(f"Initial density on process {MPI.COMM_WORLD.rank}: \n{density.vector()[:]}\n")

# Define support
support = CompiledSubDomain("near(x[0], 0.0, tol) && on_boundary", tol=1e-14)
bcs = [DirichletBC(U,Constant((0.0, 0.0)), support)]

# Define external load
load_marker = MeshFunction("size_t", msh, msh.topology().dim() - 1)
CompiledSubDomain("x[0]==l && x[1]<=2", l=nelx).mark(load_marker, 1)
ds = Measure("ds")(subdomain_data=load_marker)
f = dot(v, Constant((0.0, -1.0))) * ds(1)

# Set up the variational problem and solver (using SIMP)
k = inner(density**penal*sigma(u), grad(v)) * dx
problem = LinearVariationalProblem(k, f, u_sol, bcs=bcs)
solver = LinearVariationalSolver(problem)

# solve PDE
solver.solve()
print(f"u_sol on process {MPI.COMM_WORLD.rank}: \n{u_sol.vector()[:]}\n")

# Objective & sensitivity
objective = project(density ** penal * psi(u_sol), D).vector()[:]
print(f"Objective on process {MPI.COMM_WORLD.rank}: \n{objective}\n")
sensitivity = -penal * (density.vector()[:]) ** (penal-1) * project(psi(u_sol), D).vector()[:]
print(f"Sensitivity on process {MPI.COMM_WORLD.rank}: \n{sensitivity}\n")

When the reference code is executed through mpirun -n 2~~

Initial density on process 0: 
[0.5 0.5 0.5 0.5]

Initial density on process 1: 
[0.5 0.5 0.5 0.5]

Process 0: Solving linear variational problem.
Process 1: Solving linear variational problem.
u_sol on process 0: 
[0. 0. 0. 0. 0. 0.]

u_sol on process 1: 
[ 2.19098143e+01 -6.77188329e+01 -2.78721395e-15 -6.44827586e+01
  1.51458886e+01 -2.69893899e+01 -1.51458886e+01 -2.69893899e+01
 -2.19098143e+01 -6.77188329e+01 -9.70784783e-16 -2.81564987e+01]

Objective on process 0: 
[38.28299907  4.35108686  4.35108686 38.28299907]

Objective on process 1: 
[12.13231906 11.33439076 11.33439076 12.13231906]

Sensitivity on process 0: 
[-229.69799443  -26.10652119  -26.10652119 -229.69799443]

Sensitivity on process 1: 
[-72.79391433 -68.00634459 -68.00634459 -72.79391433]

As shown in the following results, when objective & sensitivity were compared with each local process, a difference was seen in the process of interpolation in the code built with dolfinx,
so I thought that the project_func function built in dolfinx did not perform parallel calculation.

I am wondering if there is any way to solve the problem.