Dear community,
I am trying to launch my simulations in parallel but I get some issue with boundary condition.
I attached a test code of a pure traction bar in dynamic explicit, loading in displacement. I am solving the problem :
M_lumped a = - F_int with F_int = K*u
a = (M_lumped)⁻1 * (-F_int)
I followed the guidelines regarding a previous post here. It seems ok with 1 proc but when I am used more than 1 proc the problem appears.
This is the code of the test :
#!/usr/bin/env python
# coding: utf-8
from dolfinx import log, io, mesh, fem, cpp, geometry
from petsc4py import PETSc
from mpi4py import MPI
import ufl
import numpy as np
import matplotlib.pyplot as plt
savedir = "result_mpi_1proc/"
comm = MPI.COMM_WORLD
def mpi_print(s):
print(f"Rank {comm.rank}: {s}")
#loading
t_start, t_end, nsteps = 0.0, 3., 5*800
dt = t_end/nsteps
steps = np.linspace(t_start, t_end, nsteps)
# Geometry & Mesh
Lx , Ly = 1., 0.1
msh = mesh.create_rectangle(comm, [[0., 0.], [Lx,Ly]], [10,1], mesh.CellType.quadrilateral, ghost_mode=cpp.mesh.GhostMode.none) # cpp.mesh.GhostMode.shared_facet)
ndim = msh.topology.dim
fdim = ndim - 1
# Material Parameters
E, nu = fem.Constant(msh, PETSc.ScalarType(10)), fem.Constant(msh, PETSc.ScalarType(0.0)) #[MPa]
rho = fem.Constant(msh, PETSc.ScalarType(1.))
mu = E / (2.0*(1.0 + nu)) # [MPa]
lmbda = E * nu / ((1.0 + nu)* (1.0 - 2.0*nu))
K0 = lmbda+2./3.*mu
# Space function
# index_map(n) n = 0 (vertex), n=1 (edge), n=2 (cell)
V_element=ufl.VectorElement("CG",msh.ufl_cell(),1)
V=fem.FunctionSpace(msh,V_element)
nb_cell_loc = msh.topology.index_map(2).size_local
nb_vertex_loc = msh.topology.index_map(0).size_local
nb_vertex_glo = msh.topology.index_map(0).size_global
nb_node_x = len(msh.geometry.x)
# Facets / Geometry
def right(x):
return np.isclose(x[0], Lx)
def left(x):
return np.isclose(x[0], 0)
def bottom(x):
return np.isclose(x[1], 0)
facet_bottom = mesh.locate_entities(msh, fdim, bottom)
facet_right = mesh.locate_entities(msh, fdim, right)
facet_left = mesh.locate_entities(msh, fdim, left)
facet_left_dofs = fem.locate_dofs_geometrical((V.sub(0),V.sub(0).collapse()[0]), left)
facet_right_dofs = fem.locate_dofs_geometrical((V.sub(0), V.sub(0).collapse()[0]), right)
facet_bottom_dofs = fem.locate_dofs_geometrical((V.sub(1),V.sub(1).collapse()[0]), bottom)
facet_indices, facet_markers = [], []
facet_indices.append(facet_right)
facet_indices.append(facet_left)
facet_indices.append(facet_bottom)
facet_indices = np.array(np.hstack(facet_indices), dtype=np.int32)
sorted_facets = np.argsort(facet_indices)
boundary_indices = {"left": 0, "right": 1, "top": 2, "bottom": 3}
facet_markers.append(np.full(len(facet_right ), boundary_indices["right"]))
facet_markers.append(np.full(len(facet_left ), boundary_indices["left"]))
facet_markers.append(np.full(len(facet_bottom ), boundary_indices["bottom"]))
facet_markers = np.array(np.hstack(facet_markers), dtype=np.int32)
facet_tag = mesh.meshtags(msh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
metadata = {"quadrature_degree": 4}
ds = ufl.Measure("ds", domain=msh, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=msh, metadata=metadata)
# BOUNDARY CONDITIONS
u_zero = np.array((0. ,)*msh.geometry.dim, dtype=PETSc.ScalarType)
class Load():
v0 = 3.0 #m/s
t0 = 0.00001 #s
def __init__(self, t):
self.t = t
def __call__(self,x):
values = np.zeros((ndim-1, x.shape[1]),dtype=PETSc.ScalarType)
if self.t < Load.t0:
values[0] = 0.5 * (Load.v0 / Load.t0) * self.t * self.t
else:
values[0] = Load.v0 * (self.t-Load.t0) + 0.5 * Load.v0 * Load.t0
return values
t = 0.
trac_func = fem.Function(V.sub(0).collapse()[0])
s_app = Load(t)
trac_func.interpolate(s_app)
zero_u_bottom = fem.Function(V.sub(1).collapse()[0])
with zero_u_bottom.vector.localForm() as bc_local:
bc_local.set(0.0)
zero_u_left = fem.Function(V.sub(0).collapse()[0])
with zero_u_left.vector.localForm() as bc_local:
bc_local.set(0.0)
bc_right = fem.dirichletbc(trac_func , facet_right_dofs, V.sub(0))
bc_left = fem.dirichletbc(zero_u_left, facet_left_dofs ,V.sub(0))
bc_bottom = fem.dirichletbc( zero_u_bottom,facet_bottom_dofs, V.sub(1))
bcu = [bc_left, bc_bottom ,bc_right]
u, u_p = fem.Function(V, name='Displacement'),fem.Function(V, name='Velocity') # solution u, velocity
u_pp = fem.Function(V, name='Acceleration') # acceleration
u_pp_old = fem.Function(V, name='Previous Acceleration')
v, u_trial = ufl.TestFunction(V), ufl.TrialFunction(V)
# Linear elastic model
def epsilon(u):
return ufl.sym(ufl.grad(u))
def sigma(u):
return lmbda * ufl.tr(epsilon(u)) * ufl.Identity(u.geometric_dimension()) + 2*mu*epsilon(u)
def mass(u, v):
return rho*ufl.inner(u, v)*dx
def stiffness(u, v):
return ufl.inner(sigma(u),epsilon(v))*dx
def p_ext(u):
return ufl.dot(u, trac_func)*ds(1)
mass_form= fem.form(mass(u_trial, v))
M = fem.petsc.create_matrix(mass_form)
M = fem.petsc.assemble_matrix(mass_form)
M.assemble()
# Matrix Lumping
# Compute the sum of the rows of the mass matrix
ones = fem.Function(V)
with ones.vector.localForm() as loc:
loc.set(1.0)
m_action= ufl.action(mass(u_trial, v), ones)
M_lumped = fem.petsc.assemble_vector(fem.form(m_action))
M_lumped_inv = 1./M_lumped
with io.XDMFFile(msh.comm, savedir+ "displacement.xdmf", "w") as xdmf0:
xdmf0.write_mesh(msh)
xdmf0.write_function(u,t)
E_elas = 0
inc = 0
for t in steps[:50]:
s_app.t = t
trac_func.interpolate(s_app)
if comm.rank == 0 :
print("\nIncrement: %i"%(inc))
print("Time: %.8g ; Load : %.8g"%(t,trac_func.x.array[0] ))
u.x.array[:] = u.x.array[:] + dt * u_p.x.array[:] + 0.5 * dt**2 * u_pp_old.x.array[:]
fem.set_bc(u.vector.array, bcu)
# f = fext-fint
f_int_form = fem.form(stiffness(u, v))
f_int_vec = fem.petsc.create_vector(f_int_form)
with f_int_vec.localForm() as f_local:
f_local.set(0.0)
fem.petsc.assemble_vector(f_int_vec, f_int_form)
fem.petsc.apply_lifting(f_int_vec, [mass_form], [bcu])
f_int_vec.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
fem.petsc.set_bc(f_int_vec, bcu)
# u_pp.vector.setArray(-M_lumped_inv*f_int_vec)
u_pp.vector.array = -M_lumped_inv*f_int_vec
mpi_print(f" u_pp.x.array[:] : { u_pp.x.array[:]}")
u_pp.x.scatter_forward()
u_p.x.array[:] = u_p.x.array + 0.5*dt * (u_pp.x.array + u_pp_old.x.array)
u_pp_old.x.array[:] = u_pp.x.array
E_elas = fem.assemble_scalar(fem.form(0.5*stiffness(u, u)))
#mpi_print(f"Elas: {E_elas}")
xdmf0.write_function(u, t)
inc +=1
Some results :
Thank you for your help,
Best regards,
Lamia
