Error in solution of XDMF along mesh partition

I/O
1

Hello,

I’ve been running larger simulations on an HPC machine and noticed errors along the mesh partition interfaces, see

Screen Shot 2023-04-07 at 7.02.24 PM
Screen Shot 2023-04-07 at 7.02.24 PM1678×1076 59.3 KB

At each time step, I append the solution to an existing xdmf file

with io.XDMFFile(domain.comm, filename,'w') as xdmf:
    xdmf.write_mesh(domain)

via

with io.XDMFFile(domain.comm, filename,'a') as xdmf:
        xdmf.write_function(U_n, t_n)

If it’s helpful, U_n is a vector function space of size 4. The mesh domain is generated in the same python script with gmsh. In trying to debug, I came across this post that has a similar feature. I followed the marked solution to ensure that my interior mesh is indeed connected.

Has anyone else seen this before? I’ve been unsuccessful in reproducing this with a manageable/MWE, but I will keep trying.

Thank you.

2

Without a minimal example, its hard to give guidance. How long is your current m mmcode?

3

It can be stripped down to about 400 lines of code. The image I posted is the result of the first time step.

4

Looks like a scatter_forward issue to me. Are you familiar with the concept of ghost cells ? Do you currently update them at every step ?

Also ALWAYS check your mesh, especially if you automated its generation a lot. Many things can be explained by a topological defect produced by a Voronoi-like algorithm

2 Likes
5

I would suggest posting the code, as @hawkspar says, it is most likely a ghost update issue.

6

Thank you @hawkspar and @dokken. I am hesitant to post the full code but have narrowed down my problem code to an auxiliary solve of the Poisson equation, the solution of which I use in my weak form solve at each time step. In implementing this, I followed https://github.com/FEniCS/dolfinx/blob/main/python/demo/demo_stokes.py as closely as possible.

The code is:

def PoissonSolve(domain, uM, dt, n):
     
    # define finite element space
    P1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1)
    W = fem.FunctionSpace(domain, P1)

    # test function, trial solution
    psi = ufl.TrialFunction(W)
    w = ufl.TestFunction(W)

    # weak form
    lhs = -ufl.inner(ufl.grad(psi), ufl.grad(w)) *ufl.dx
    
    rhs =  ufl.inner(uM, ufl.grad(w))*ufl.dx
    rhs += ufl.inner(ufl.dot(uM, n), w)*ufl.ds
    a = fem.form(lhs)
    L =fem.form( rhs)

    # Assemble LHS matrix and RHS vector
    A = fem.petsc.assemble_matrix(a)
    A.assemble()
    b = fem.petsc.assemble_vector(L)
    b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
    
    nullspace = PETSc.NullSpace().create(constant=True)
    A.setNullSpace(nullspace)
    
    # Create and configure solver
    ksp = PETSc.KSP().create(domain.comm)
    ksp.setOperators(A)
    ksp.setType("preonly")

    # Configure MUMPS to handle nullspace
    pc = ksp.getPC()
    pc.setType("lu")
    pc.setFactorSolverType("mumps")
    pc.setFactorSetUpSolverType()
    pc.getFactorMatrix().setMumpsIcntl(icntl=24, ival=1)
    pc.getFactorMatrix().setMumpsIcntl(icntl=25, ival=0)

    # Compute the solution
    U = fem.Function(W)

    ksp.solve(b, U.vector)
    psi = U
    
    return psi

And after the mention of scatter_forward, I do see that I b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE), and never put it back forward. With that, I don’t know where is appropriate to scatter_forwarad with setting the constant nullspace.

7

You need to call scatter forward after you have solved the problem, as done in:
https://jsdokken.com/dolfinx-tutorial/chapter2/heat_code.html?highlight=scatter#solving-the-time-dependent-problem

8

Wow, thank you. I can’t believe I missed that. I was looking around the forum a bit more - do I also need nullspace = PETSc.NullSpace().create(constant=True, comm=domain.comm)? Or is that done implicitly?

9

You should create the null space, and remove it from the rhs vector if you solve a singular problem with a direct solver. See:

10

I’m sorry, I was referring to comm=domain.comm in the call
nullspace = PETSc.NullSpace().create(constant=True, comm=domain.comm)
as opposed to not including comm=domain.comm, which I was doing before:
nullspace = PETSc.NullSpace().create(constant=True) ← no comm=

11

As long as the space is the constant space, i dont think it matters. It would be relevant for non-constant nullspaces, or if you use other comms than comm_world when you define the domain.

1 Like