Hello everyone,
I have updated my DOLFINx installation from 0.5.1 to 0.9.0. So far, I had never used parallelization through MPI. Instead, I used to set the environment variable OMP_NUM_THREADS (set and export) to a value >1 and I could see sufficient increase of performance.
After updating DOLFINx, I have noticed that this does not work in my installation anymore. When I run a simulation, only one thread is involved, no matter the value I give to such an environment variable.
Consequently, I have tried to use MPI parallelization (for my first time). However, when I do mpirun -n 2 python3 file.py, I oberserve in Paraview that there is a part of the mesh that does not deform. When I do python3 file.py, all the mesh deforms as expected (and the simulation takes longer).
The particular case I am studying is an implementation of spherical contact with the penalty method. I copy the code below. May the problem be related with: (i) the implementation of mpi in my code?; (ii) the visualization? Thank you in advance.
from __future__ import print_function
import numpy as np
from numpy import array
import dolfinx
import ufl
from mpi4py import MPI
import dolfinx.io
import math
from petsc4py.PETSc import ScalarType
from petsc4py import PETSc
import petsc4py
import matplotlib.pyplot as plt
from petsc4py import PETSc
default_scalar_type = PETSc.ScalarType
from dolfinx import fem
import dolfinx.fem.petsc
import dolfinx.nls.petsc
name = "Output.xdmf" # Name of the output file
name1 = "Output.h5"
T = 0.25 # Final time
steps = 0 #Inizialization of the variable.
# Time parameters
tiempo = 0.0
tot_steps = 5;
dt = T / tot_steps
#Import mesh & subdomains from .msh file:
filename="Box"
from dolfinx.io import gmshio
mesh, cell_markers, facet_markers = gmshio.read_from_msh(filename+'.msh', MPI.COMM_WORLD, gdim=3)
metadata = {"quadrature_degree": 2}
dx = ufl.Measure('dx')(domain=mesh, subdomain_data=cell_markers, metadata=metadata) #Integration measures for each subdomain
# FUNCTION SPACES
V = dolfinx.fem.functionspace(mesh, ("CG",2,(mesh.geometry.dim, )))
VV = dolfinx.fem.functionspace(mesh,('CG',1, (mesh.geometry.dim, )))
# FUNCTIONS IN THE FUNCTION SPACE V (Function - Trial - Test).
u = dolfinx.fem.Function(V) #Unkown.
v = ufl.TestFunction(V)
########## BOUNDARY CONDITIONS ############################################
def bnd_encastre(x):
tol = 1e-4
return np.logical_or.reduce((
np.isclose(x[2], 0, atol=tol), # Fix z = 0
np.isclose(x[0], 0, atol=tol), # Fix x = 0
np.isclose(x[1], 0, atol=tol), # Fix y = 0
np.isclose(x[0], 1, atol=tol), # Fix x = 1
np.isclose(x[1], 1, atol=tol) # Fix y = 1
))
bnd_encastre_dofs0 = dolfinx.fem.locate_dofs_geometrical(V,bnd_encastre)
u_encastre_vect0 = np.array((0,) * mesh.geometry.dim, dtype=default_scalar_type)
bc_boundar = dolfinx.fem.dirichletbc(u_encastre_vect0,bnd_encastre_dofs0,V)
bc = [bc_boundar]
########## CONTACT SURFACE #####################################
def bnd_side(x):
return np.isclose(x[2],1)
upper_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim - 1, bnd_side)
facet_upper_tag = dolfinx.mesh.meshtags(mesh, mesh.topology.dim - 1, upper_facets, 1)
ds_contact = ufl.Measure('ds', domain=mesh, subdomain_data=facet_upper_tag, subdomain_id=1)#, metadata=metadata)
#######################################################################
########## OBSTACLE #####################################
x = ufl.SpatialCoordinate(mesh)
d_ = dolfinx.fem.Constant(mesh, 0.0)
def obstacle(x):
h0 = 1.1
Rad = 0.3
return -d_ + h0 + 1/2/Rad * ((x[0]-0.5)**2 + (x[1]-0.5)**2)
#######################################################################
## PRIMARY FIELDS
dd = len(u)
Id = ufl.Identity(dd)
F = Id + ufl.grad(u)
## CONSTITUTIVE MODEL:
######### NEO-HOOKEAN ##################################################################
def Piso(F,i):
j,k,l,m = ufl.indices(4)
FinvT = ufl.inv(F).T
A1 = ufl.as_tensor(Id[j,l]*Id[k,m],(j,k,l,m))
A2 = ufl.as_tensor(-1/3*FinvT[j,k]*F[l,m],(j,k,l,m))
Pfourth = A1+A2
Fiso = 1/ufl.det(F)**(1/3)*F
Piso_ = G[i]*Fiso
return 1/(ufl.det(F)**(1/3))*ufl.as_tensor(Pfourth[j,k,l,m]*Piso_[l,m],(j,k))
def Pvol(F,i):
return kappa[i]*ufl.det(F)*(ufl.det(F)-1)*ufl.inv(F).T
G = [1000,0]
poisson=[0.495,0,0]
kappa = [2*G[0]*(1+poisson[0])/3/(1-2*poisson[0]),0,0]
########################################################################################
####### VARIATIONAL PROBLEM - WEAK FORM ################################################
pen = 1e4
Res_u = (ufl.inner(Piso(F, 0), ufl.grad(v))\
+ ufl.inner(Pvol(F, 0), ufl.grad(v))) * dx\
+ pen * (((u[2]+1)-obstacle(x)) + np.abs((u[2]+1)-obstacle(x)))/2 * v[2] * ds_contact(1)
# Setup Non-linear variational problem
problem = fem.petsc.NonlinearProblem(Res_u,u,bc)
from dolfinx import nls
solver_problem = nls.petsc.NewtonSolver(mesh.comm, problem)
solver_problem.atol = 1e-8
solver_problem.rtol = 1e-8
solver_problem.convergence_criterion = "incremental"
import os
if os.path.exists(name):
os.remove(name)
os.remove(name1)
# Save results into an .xdmf
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
file_results = dolfinx.io.XDMFFile(MPI.COMM_WORLD,name,'w')
# Write the mesh and subdomains info
file_results.write_mesh(mesh)
###### TIME-INCREMENTAL SOLVING LOOP ########################
while (tiempo < T):
#Display report simulation
from dolfinx import log
log.set_log_level(log.LogLevel.INFO)
print("Iteration info", flush=True)
d_.value = tiempo
print("dt: ",dt)
print("Step: ",steps)
print("Tiempo: ",tiempo,"/",T)
solver_problem.solve(u)
log.set_log_level(log.LogLevel.OFF)
# Write to .xdmf results file
uVector=dolfinx.fem.Function(VV); uVector.name = "Displacement"; u_expr = dolfinx.fem.Expression(u, VV.element.interpolation_points()); uVector.interpolate(u_expr)
file_results.write_function(uVector,tiempo)
steps += 1
tiempo += dt
file_results.close()
print("Finished")
Note: the mesh is taken as input, but it is a simple unit square.
