# MPI versus OMP\_NUM\_THREADS & issue with mesh not deforming in some regions **URL:** \ **Category:** dolfinx\ **Created:** [October 25, 2024, 12:45pm UTC](https://fenicsproject.discourse.group/t/mpi-versus-omp-num-threads-issue-with-mesh-not-deforming-in-some-regions/16113 "2024-10-25T12:45:16Z")\ **Posts on this page:** 8\ **Page:** 1
**Author:** ![MorenoMiguelES](https://avatars.discourse-cdn.com/v4/letter/m/35a633/32.png) [@MorenoMiguelES](https://fenicsproject.discourse.group/u/MorenoMiguelES)\ **Post date:** [October 25, 2024, 12:45pm UTC](https://fenicsproject.discourse.group/t/mpi-versus-omp-num-threads-issue-with-mesh-not-deforming-in-some-regions/16113/1 "2024-10-25T12:45:16Z")
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. ![Captura de pantalla 2024-10-25 a las 14.33.35](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/5/5fb1efde1178761753cf5201ec7e3105e73595d9.jpeg) ```auto 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. ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [October 25, 2024, 1:10pm UTC](https://fenicsproject.discourse.group/t/mpi-versus-omp-num-threads-issue-with-mesh-not-deforming-in-some-regions/16113/2 "2024-10-25T13:10:50Z")
I can’t spot anything obviously wrong in your code. Could you: 1. try with a built in box mesh? This will just tell us if it is reading in from gmsh that is the problem. 2. If it doesn’t work with the built in mesh, it will make it easier to debug this for us. Secondly, can you check that boundary conditions would be applied correctly to a dolfinx.fem.Function and export it (without solving a PDE). ---
**Author:** ![MorenoMiguelES](https://avatars.discourse-cdn.com/v4/letter/m/35a633/32.png) [@MorenoMiguelES](https://fenicsproject.discourse.group/u/MorenoMiguelES)\ **Post date:** [October 25, 2024, 2:04pm UTC](https://fenicsproject.discourse.group/t/mpi-versus-omp-num-threads-issue-with-mesh-not-deforming-in-some-regions/16113/3 "2024-10-25T14:04:53Z")
I just tried using the built in mesh (`mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD,10,10,10)`), and the problem remains. I have also observed that, when doing `mpirun -n 2 python3 file.py`, it does not work. When doing `mpirun -n 1 python3 file.py`, it does work. Given this, I would say it has not to do with the BCs (encastre on the lateral and bottom sides), or may I be wrong? Anyhow, I am trying to verify them. ---
**Author:** ![MorenoMiguelES](https://avatars.discourse-cdn.com/v4/letter/m/35a633/32.png) [@MorenoMiguelES](https://fenicsproject.discourse.group/u/MorenoMiguelES)\ **Post date:** [October 25, 2024, 3:22pm UTC](https://fenicsproject.discourse.group/t/mpi-versus-omp-num-threads-issue-with-mesh-not-deforming-in-some-regions/16113/4 "2024-10-25T15:22:54Z")
Continuation: If I simply remove (all) the Dirichlet BC and keep the penalty, the problem also remains: a half of the mesh does not move at all, and the other part (blue) moves as one would expect from free “rigid body motion”. ![Captura de pantalla 2024-10-25 a las 17.21.45](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/c/c73aea819227d184e2e5f761d852cd1c67dfda2c.jpeg) ---
**Author:** ![MorenoMiguelES](https://avatars.discourse-cdn.com/v4/letter/m/35a633/32.png) [@MorenoMiguelES](https://fenicsproject.discourse.group/u/MorenoMiguelES)\ **Post date:** [February 20, 2025, 10:37am UTC](https://fenicsproject.discourse.group/t/mpi-versus-omp-num-threads-issue-with-mesh-not-deforming-in-some-regions/16113/5 "2025-02-20T10:37:18Z")
> [@MorenoMiguelES](#): > > 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. Is there any news regarding this? **I don’t want to use MPI, just to have the same performance that I used to have in 0.5.1 but now in 0.9.0 (for a same simulation).** A summary next: In dolfinx 0.5.1 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 to 0.9.0, I have noticed that this does not work in my installation anymore. When I run a simulation, only one thread is involved and a same simulation that I used to run in 0.5.1 runs much slower in 0.9.0. I use MAC. With 0.5.1 (without mpi) I used to see a CPU use above 100%. With 0.9.0 I see that it always stays \<= 100%. Thank you and apologises for insisting. ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [February 20, 2025, 11:21am UTC](https://fenicsproject.discourse.group/t/mpi-versus-omp-num-threads-issue-with-mesh-not-deforming-in-some-regions/16113/6 "2025-02-20T11:21:54Z")
In DOLFINx we have no special threading algorithms implemented for performance increase. We rely on using MPI for parallelization. ---
**Author:** ![MorenoMiguelES](https://avatars.discourse-cdn.com/v4/letter/m/35a633/32.png) [@MorenoMiguelES](https://fenicsproject.discourse.group/u/MorenoMiguelES)\ **Post date:** [February 20, 2025, 2:40pm UTC](https://fenicsproject.discourse.group/t/mpi-versus-omp-num-threads-issue-with-mesh-not-deforming-in-some-regions/16113/7 "2025-02-20T14:40:18Z")
It seems that the issue is related to the default configuration of the `nls.petsc-NewtonSolver`, which seems to be different in 0.9.0 with respect to 0.5.1. I tried some pre-conditioners when running the code in 0.9.0 and then the computational performance does better replicate my simulation in 0.5.1. ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [February 20, 2025, 3:40pm UTC](https://fenicsproject.discourse.group/t/mpi-versus-omp-num-threads-issue-with-mesh-not-deforming-in-some-regions/16113/8 "2025-02-20T15:40:57Z")
That means that how PETSc uses threading have changed.