How to repeatedly solve FEA on procedural meshes then save the data

I/O
1

I am trying to solve a reliability analysis problem where the FEA needs to be solved many times while changing the mesh and some material parameters every time the problem is setup and solved. All the code works well when I run it once, but any time I try and run it multiple times in a loop or execute the python script multiple times in a row from a seperate python file, at the end of the loops it gives me a memory error from PETSC, specifically:
“[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range”
When I check the container memory usage its well under what it can use (~4GB vs 120 available). I tried using the python garbage collector at the end of each analysis and deleting all variables using the “del” keyword, however that caused the same PETSC error to occur, but after the first loop instead of the last. I also tried using the

Any help to solve this would be greatly appreciated. I made a simplified code based on one of the tutorials for linear elasticity and included it below, I am using Fenicsx 0.8.0 and Ubuntu 18.04.6 for the docker container.

from dolfinx import log, mesh, fem, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
import ufl
import numpy as np

for i in range(0,5):
    # physical variables
    N = 20
    L = 1.0
    D = 2*L

    E = default_scalar_type(10.0)
    nu = default_scalar_type(0.3)
    DmatPlaneStress = E/(1-nu*nu)*ufl.as_matrix([[1,nu,0],[nu,1,0],[0,0,(1-nu)]])
    Dmat = DmatPlaneStress

    # meshing variables
    elementOrder = 2

    # create mesh
    domain = mesh.create_rectangle(
        MPI.COMM_WORLD,
        [[-L, -D], [L, 0]],
        [N, int(round(N*D/L))],
        cell_type=mesh.CellType.triangle,
    )

    # find the top surface
    def top_surface(x):
        return np.isclose(x[1],0.0)

    # define other boundaries
    def bottom(x):
        return np.isclose(x[1], -D)

    # locate all exterior faces (for later boolean operations)
    tdim = domain.topology.dim
    fdim = tdim - 1
    domain.topology.create_connectivity(fdim, tdim)
    boundary_facets_all = mesh.exterior_facet_indices(domain.topology)

    # locate different regions of facets where the boundary conditions will be applied
    topFacets = mesh.locate_entities_boundary(domain, fdim, top_surface)
    bottomZ_facets = mesh.locate_entities_boundary(domain, fdim, bottom)

    # define the mesh marker for each facets on the mesh
    num_facets = domain.topology.index_map(fdim).size_local
    markers = np.zeros(num_facets, dtype=np.int32)
    top_mark,bottomZ_mark = 1,2
    markers[topFacets],markers[bottomZ_facets] = top_mark,bottomZ_mark
    facet_marker = mesh.meshtags(domain, fdim, np.arange(num_facets, dtype=np.int32), markers)

    # define the function space
    V = fem.functionspace(domain, ("Lagrange", elementOrder, (domain.geometry.dim,)))

    # define variables for energy density function
    mu = fem.Constant(domain, E / (2 * (1 + nu)))
    lmbda = fem.Constant(domain, E * nu / ((1 + nu) * (1 - 2 * nu)))

    # define the displacment boundary condition numerically, then assign them to the boundary facets
    u_fixedAll = np.array([0, 0], dtype=default_scalar_type)
    bc_bottom = fem.dirichletbc(u_fixedAll, fem.locate_dofs_topological(V, fdim, facet_marker.find(bottomZ_mark)),V)
    bcs = [bc_bottom]

    # define the traction to be 0
    T = fem.Constant(domain, default_scalar_type((0, 5)))
    # define the constant body force value
    f = fem.Constant(domain, default_scalar_type((0, 0)))

    def TensorToVec(m):
        return ufl.as_vector([m[0,0],m[1,1],.5*m[0,1]+.5*m[1,0]])

    def VecToTensor(vec):
        return ufl.as_matrix([[vec[0],vec[2]],[vec[2],vec[1]]])

    # return the strain from the displacement symbolically
    def epsilon(z):
        return ufl.sym(ufl.grad(z))

    # return the stress from the displacement symbolically
    def sigma(z):
        sig = VecToTensor(ufl.dot(Dmat,TensorToVec(epsilon(z))))
        return sig

    # define the test and solution space
    v = ufl.TestFunction(V)
    u = fem.Function(V)
    du = ufl.TrialFunction(V)

    # define the integration measure with traction integral over all facets with value 2 and set quadrature degree to 4
    ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_marker)
    dx = ufl.Measure("dx", domain=domain)

    # calculate the residual of the equations
    residual = ufl.inner(sigma(u), epsilon(v)) * dx - ufl.dot(f, v) * dx - ufl.dot(T, v) * ds(top_mark)

    # define tangent tensor
    a = ufl.derivative(residual, u, du) # for taking the derivative with respect to a function

    jit_options ={"cffi_extra_compile_args":["-O3","-ffast-math"]}
    problem = NonlinearProblem(residual, u, bcs, a,jit_options=jit_options)

    solver = NewtonSolver(MPI.COMM_WORLD, problem)
    solver.convergence_criterion = "incremental"
    solver.rtol = 1e-8
    solver.atol = 1e-8
    solver.max_it = 100
    solver.report = True

    ksp = solver.krylov_solver
    opts = PETSc.Options()
    option_prefix = ksp.getOptionsPrefix()
    opts[f"{option_prefix}ksp_type"] = "preonly"
    opts[f"{option_prefix}pc_type"] = "lu"
    opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
    ksp.setFromOptions()

    totSteps = 2
    dispHist = np.zeros(shape=[totSteps-1])
    loadHist = np.zeros(shape=[totSteps-1]) 

    # solve the problem over several timesteps increasing the penalization slowly
    log.set_log_level(log.LogLevel.INFO)
    for n in range(1, totSteps):
        try:
            (iter, converged) = solver.solve(u)
        except: # Break the loop if solver fails
            print( "Ended Early")
            break

    print("max displacement: " + str(np.max(np.abs(u.x.array))))
2

There are usually very few good reasons for putting the form definition (and solver definition) within a for loop. In your initial text you say you change the mesh, does that mean:

  • reading in a new mesh from file
  • Deforming an existing mesh?
3

I read in a new mesh with slightly modified geometry and slightly change some of the elasticity values, I used a nonlinear analysis because the real simulation is a contact FEA with a penalty approach

4

I have encountered this problem too in geometry optimization where I had to create a new mesh on every iteration. At that time, the quickest stopgap solution was to run each iteration through joblib or multiprocessing module. Of course they are used for parallelization but their feature to create and destroy isolated processes prevented the accumulation of MPI communicators that lead to this error.

5

Changing elasticity values is no good reason for putting things inside the loop.
This is for instance illustrated in: Efficient usage of the Unified Form Language — FEniCS Workshop

For changing your mesh, things become a bit more complicated, due to a few things:

  • One should not use MPI.COMM_WORLD internally in a library, but duplicate the communicator (https://www.mpich.org/static/docs/v3.3/www3/MPI_Comm_dup.html)

    . In particular, no library routine should use MPI_COMM_WORLD as the communicator; instead, a duplicate of a user-specified communicator should always be used. For more information, see Using MPI, 2nd edition

Thus we need to duplicate the communicator sent into the mesh.

Garth made several changes to dolfinx.fem.Function/dolfinx.la.Vector to address this.
For instance, on v0.10.0 (using docker), the following runs with no issue:

from mpi4py import MPI
import dolfinx.fem.petsc
import ufl

for i in range(10_000):
    mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
    V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
    uh = dolfinx.fem.Function(V)
    u = ufl.TrialFunction(V)
    v = ufl.TestFunction(V)
    a = ufl.inner(u, v)*ufl.dx
    x = ufl.SpatialCoordinate(mesh)
    f = x[0]**2
    L = ufl.inner(f, v) * ufl.dx
    problem = dolfinx.fem.petsc.LinearProblem(a, L, u=uh, petsc_options={"ksp_type":"preonly", "pc_type": "lu",
                                                                         "pc_factor_mat_solver_type": "mumps",
                                                                         "ksp_error_if_not_converged": True},
                                                                         petsc_options_prefix=f"problem_{i}")
    problem.solve()
    if MPI.COMM_WORLD.rank == 0:
        print(i, problem.solver.getConvergedReason(), flush=True)
    del problem

Edit: Another example of how to clear communicators is shown in the comment: Reopening the discussion regarding `MPI_Comm_dup` · Issue #3061 · FEniCS/dolfinx · GitHub

2 Likes
6

Awesome, thank you!

I could not get it to work in fenicsx 0.8.0, but when I updated to 0.10.0 and implemented the changes you recommended it works! I added a working version of the example code posted in the question in case anyone would find it helpful

from dolfinx import log, mesh, fem, plot, io, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
from mpi4py import MPI
from petsc4py import PETSc
import ufl
import numpy as np
from dolfinx.io import XDMFFile
import gc

import os
os.environ['OMP_NUM_THREADS'] = '2'


def solve_poisson(comm, domain, iter):
    # meshing variables
    elementOrder = 1

    # define the function space
    V = fem.functionspace(domain, ("Lagrange", 1))

    # define the pointwise values in the trial/test function space
    u = fem.Function(V)
    v = ufl.TestFunction(V)
    du = ufl.TrialFunction(V)

    f = fem.Constant(domain, 1.0)

    # find the facets (fdim) in the domain which are on the clamped boundary
    fdim = domain.topology.dim - 1
    
    domain.topology.create_connectivity(fdim, domain.topology.dim)
    boundary_facets = mesh.exterior_facet_indices(domain.topology)

    x = ufl.SpatialCoordinate(domain)

    # define the displacment boundary condition numerically, then assign them to the boundary facets
    u_bc = fem.Constant(domain, 0.0)
    bc = fem.dirichletbc(u_bc, fem.locate_dofs_topological(V, fdim, boundary_facets), V)

    # dx = ufl.Measure("dx", domain=domain)
    residual = ufl.inner(ufl.grad(u),ufl.grad(v)) * ufl.dx - f*v*ufl.dx
    # define tangent tensor
    J = ufl.derivative(residual, u, du) # for taking the derivative with respect to a function

    prefix = "sim" + str(iter)
    jit_options ={"cffi_extra_compile_args":["-O3","-ffast-math"]}
    problem = NonlinearProblem(residual, u,bcs=[bc],petsc_options_prefix=prefix, J=J,jit_options=jit_options)

    # log.set_log_level(log.LogLevel.INFO)
    log.set_log_level(log.LogLevel.WARNING)
    problem.solve()

    varU = np.var(problem.x.array)

    return varU

def Mesh(comm, meshName):
    with XDMFFile(comm, meshName+".xdmf", "r") as xdmf:
        domain = xdmf.read_mesh(name="Grid")
        xdmf.close()

    return domain

# create meshes
directery = "/.CODE/"
mshName = directery+"singleBolt"

comm = MPI.COMM_WORLD
meshCount = 100

# sole the problem multiple times (ON NEW MESHES)
for i in range(0,meshCount):
    currComm = comm.Dup()
    msh = Mesh(currComm,mshName)
    solve = solve_poisson(currComm, msh,i)
    print("Value at mesh: " + str(i+1) + " " + str(solve))

    currComm.Free()
    del msh
    gc.collect()