PETSc produces inconsistent surface normals in parallel

[gthyagi]

Hi @knepley,

I've observed that when running in parallel, PETSc produces surface normals with flipped signs compared to serial execution, where all normals are consistently oriented. This inconsistency seems to be causing anomalous pressure and velocity fields. Could you please take a look?

Setup: This is a basic Stokes problem where a natural BC is applied on an internal boundary. Natural BC aligns in the surface normal orientation.

Thanks!

Serial output:
All normals (N) orient in same direction. Arrow represent velocity field and background color represent pressure field.

bdres - equation 18 X / N (2.96e-01, 6.70e-01, 0.00e+00 / 0e+00, 1e+00, 0.00e+00 ) -> 0.000000e+00, -1.000000e+00, 
bdres - equation 18 X / N (1.67e-01, 6.70e-01, 0.00e+00 / 0e+00, 1e+00, 0.00e+00 ) -> 0.000000e+00, -1.000000e+00, 
bdres - equation 18 X / N (3.76e-02, 6.70e-01, 0.00e+00 / 0e+00, 1e+00, 0.00e+00 ) -> 0.000000e+00, -1.000000e+00, 
bdres - equation 18 X / N (6.29e-01, 6.70e-01, 0.00e+00 / 0e+00, 1e+00, 0.00e+00 ) -> 0.000000e+00, -1.000000e+00, 
bdres - equation 18 X / N (5.00e-01, 6.70e-01, 0.00e+00 / 0e+00, 1e+00, 0.00e+00 ) -> 0.000000e+00, -1.000000e+00, 
bdres - equation 18 X / N (3.71e-01, 6.70e-01, 0.00e+00 / 0e+00, 1e+00, 0.00e+00 ) -> 0.000000e+00, -1.000000e+00, 
bdres - equation 18 X / N (9.62e-01, 6.70e-01, 0.00e+00 / 0e+00, 1e+00, 0.00e+00 ) -> 0.000000e+00, -1.000000e+00, 
bdres - equation 18 X / N (8.33e-01, 6.70e-01, 0.00e+00 / 0e+00, 1e+00, 0.00e+00 ) -> 0.000000e+00, -1.000000e+00, 
bdres - equation 18 X / N (7.04e-01, 6.70e-01, 0.00e+00 / 0e+00, 1e+00, 0.00e+00 ) -> 0.000000e+00, -1.000000e+00,

Parallel (2 cpus) output:

bdres - equation 18 X / N (6.29e-01, 6.70e-01, 0.00e+00 / 0e+00, 1e+00, 0.00e+00 ) -> 0.000000e+00, -1.000000e+00, 
bdres - equation 18 X / N (2.96e-01, 6.70e-01, 0.00e+00 / 0e+00, 1e+00, 0.00e+00 ) -> 0.000000e+00, -1.000000e+00, 
bdres - equation 18 X / N (5.00e-01, 6.70e-01, 0.00e+00 / 0e+00, 1e+00, 0.00e+00 ) -> 0.000000e+00, -1.000000e+00, 
bdres - equation 18 X / N (1.67e-01, 6.70e-01, 0.00e+00 / 0e+00, 1e+00, 0.00e+00 ) -> 0.000000e+00, -1.000000e+00, 
bdres - equation 18 X / N (3.71e-01, 6.70e-01, 0.00e+00 / 0e+00, 1e+00, 0.00e+00 ) -> 0.000000e+00, -1.000000e+00, 
bdres - equation 18 X / N (3.76e-02, 6.70e-01, 0.00e+00 / 0e+00, 1e+00, 0.00e+00 ) -> 0.000000e+00, -1.000000e+00, 
_bdres - equation 18 X / N (9.62e-01, 6.70e-01, 0.00e+00 / -0e+00, -1e+00, 0.00e+00 ) -> 0.000000e+00, -1.000000e+00,_ 
bdres - equation 18 X / N (9.62e-01, 6.70e-01, 0.00e+00 / 0e+00, 1e+00, 0.00e+00 ) -> 0.000000e+00, -1.000000e+00, 
_bdres - equation 18 X / N (8.33e-01, 6.70e-01, 0.00e+00 / -0e+00, -1e+00, 0.00e+00 ) -> 0.000000e+00, -1.000000e+00,_ 
bdres - equation 18 X / N (8.33e-01, 6.70e-01, 0.00e+00 / 0e+00, 1e+00, 0.00e+00 ) -> 0.000000e+00, -1.000000e+00, 
_bdres - equation 18 X / N (7.04e-01, 6.70e-01, 0.00e+00 / -0e+00, -1e+00, 0.00e+00 ) -> 0.000000e+00, -1.000000e+00,_ 
bdres - equation 18 X / N (7.04e-01, 6.70e-01, 0.00e+00 / 0e+00, 1e+00, 0.00e+00 ) -> 0.000000e+00, -1.000000e+00,

UW script:

import underworld3 as uw
from underworld3.systems import Stokes
import sympy

mesh = uw.discretisation.Mesh(mesh_fname, boundaries=boundaries, 
                              boundary_normals=boundary_normals_2D,
                              markVertices=True, useMultipleTags=True, useRegions=True,
                              coordinate_system_type=uw.coordinates.CoordinateSystemType.CARTESIAN, )

v_soln = uw.discretisation.MeshVariable(r"u", mesh, 2, degree=2)
p_soln = uw.discretisation.MeshVariable(r"p", mesh, 1, degree=1, continuous=True)

# Create Stokes object
stokes = Stokes(mesh, velocityField=v_soln, pressureField=p_soln,) 
stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
stokes.constitutive_model.Parameters.viscosity = 1.0
stokes.saddle_preconditioner = 1.0

# bc's
stokes.add_natural_bc(sympy.Matrix([0.0, -1.0]), "Internal")
stokes.add_essential_bc(sympy.Matrix([sympy.oo, 0.0]), "Top")
stokes.add_essential_bc(sympy.Matrix([sympy.oo, 0.0]), "Bottom")
stokes.add_essential_bc(sympy.Matrix([0.0,sympy.oo]), "Left")
stokes.add_essential_bc(sympy.Matrix([0.0,sympy.oo]), "Right")

stokes.bodyforce = sympy.Matrix([0, 0])

stokes.tolerance = 1e-6
stokes.solve(verbose=False, debug=True, debug_name=f'Natural_BC_Bug_{uw.mpi.size}')

[bknight1]

There seems to be additional values in list on multiple CPU's, and all the ones with the wrong normals are '_bdres', whilst the correct ones are 'bdres'.

Your code snippet also won't work as mesh_fname, boundaries, boundary_normals_2D aren't defined.

[gthyagi]

I intentionally did not include those details to explain issue to Matt with less details.

Here are the script that will reproduce the error. I am using latest dev branch.
Stokes problem script:
https://github.com/gthyagi/UW3_Annulus_Spherical_Benchmarks/blob/main/Annulus_Benchmarks/Natural_BC_Bug_with_Int_BD.py
Create mesh for the problem:
https://github.com/gthyagi/UW3_Annulus_Spherical_Benchmarks/blob/main/Annulus_Benchmarks/create_natural_bc_mesh.py

Run Natural_BC_Bug_with_Int_BD.py and it automatically calls create_natural_bc_mesh.py

[gthyagi]

@bknight1 I added _ before bdres to highlight the incorrect normals.

[gthyagi]

Here are the full output files from the jit compiler.

AMVXD_debug_1_cpus_n_petsc_quad.txt

CWFWE_debug_2_cpus_n_petsc_quad.txt

[knepley]

This had happened in the past, but August 30th, 2024 I fixed this bug for PyLith. Is it possible that you have PETSc older than this?

Also, what function are you using to extract the internal boundary? If you just label it, rather than creating a submesh, then it is never consistently oriented.

[gthyagi]

I was originally using version 3.21.3, which didn’t include your fix. Now I'm on 3.22.2 with your fix, but I'm still encountering the same issue. As you mentioned, the problem is due to merely labeling the internal boundary. I create my mesh with an internal boundary label in gmsh. How can I create a submesh for the internal boundary? Is there any example script to do that?Thanks for your suggestions.

[gthyagi]

I intentionally did not include those details to explain issue to Matt with less details.

Here are the script that will reproduce the error. I am using latest dev branch. Stokes problem script: https://github.com/gthyagi/UW3_Annulus_Spherical_Benchmarks/blob/main/Annulus_Benchmarks/Natural_BC_Bug_with_Int_BD.py Create mesh for the problem: https://github.com/gthyagi/UW3_Annulus_Spherical_Benchmarks/blob/main/Annulus_Benchmarks/create_natural_bc_mesh.py

Run Natural_BC_Bug_with_Int_BD.py and it automatically calls create_natural_bc_mesh.py

@lmoresi Here are the scripts.
Your method might work if we retrieve GammaNorm from each MPI rank. In shallow space, the same point can end up with different normals depending on which processor owns it.

[gthyagi]

This how normals are oriented in each rank. Red points are degree two mesh variable coordinates and blue points are surface integration points.

Rank 0 (normals)

Rank 1 (normals)

Generated normals are not consistently with rank decomposition.

[lmoresi]

@knepley - I understand the idea of the sub mesh construction, but I am not sure how to stuff this back into the existing mesh so that SNES can see it as a surface to build surface integrals on. This is a surface that was handed to DMPlex in the first place via gmsh (an embedded surface). Can you expand on your thoughts here ?

[knepley]

Sure. Here is what is happening. We mark a series of faces in the some volumetric mesh. The faces are oriented randomly, determined by the mesh generator. If we want that marked surface to have a consistent normal, we have to go into the mesh data structure and change the canonical orientation for those marked faces (updating the cells they are connected to and the points in their closure in the process). You should be able to do this by calling DMPlexOrientLabel() on your label first. Note that this changes your mesh. Can you try it out and see if it works for you?

[gthyagi]

@knepley I attempted to call DMPlexOrientLabel(), which requires two arguments—the DM and the DMLabel. Here’s how I obtained them:

label_int = self.dm.getLabel("Internal")
DMPlexOrientLabel(self.dm, label_int)

The label_int contains both vertices and edges of the internal boundary. However, this results in the error:

[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: Invalid argument
[0]PETSC ERROR: Wrong type of object: Parameter # 1
[0]PETSC ERROR: #1 DMLabelGetValueIS() at [/Users/tgol0006/manual_install_pkg/venv_uw3_git_gthyagi_41224/petsc/src/dm/label/dmlabel.c:1212](http://localhost:8888/lab/tree/uw_folder/UW3_Annulus_Spherical_Benchmarks/Annulus_Benchmarks/manual_install_pkg/venv_uw3_git_gthyagi_41224/petsc/src/dm/label/dmlabel.c#line=1211)
[0]PETSC ERROR: #2 CreateCellAndFaceIS_Private() at [/Users/tgol0006/manual_install_pkg/venv_uw3_git_gthyagi_41224/petsc/src/dm/impls/plex/plexorient.c:701](http://localhost:8888/lab/tree/uw_folder/UW3_Annulus_Spherical_Benchmarks/Annulus_Benchmarks/manual_install_pkg/venv_uw3_git_gthyagi_41224/petsc/src/dm/impls/plex/plexorient.c#line=700)
[0]PETSC ERROR: #3 DMPlexOrientLabel() at [/Users/tgol0006/manual_install_pkg/venv_uw3_git_gthyagi_41224/petsc/src/dm/impls/plex/plexorient.c:726](http://localhost:8888/lab/tree/uw_folder/UW3_Annulus_Spherical_Benchmarks/Annulus_Benchmarks/manual_install_pkg/venv_uw3_git_gthyagi_41224/petsc/src/dm/impls/plex/plexorient.c#line=725)

Which corresponds to the following PETSc calls:

PetscCall(CreateCellAndFaceIS_Private(dm, label, &cellIS, &faceIS)); 726
PetscCall(DMLabelGetValueIS(label, &valueIS)); 701
PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 1); 1212

I even tried removing the vertices so that only edges remain in the label, but I still encounter the same error. Any insights or suggestions to resolve this would be greatly appreciated. What information does DMLabel should contain here? Thanks!

I am using petsc 3.22.2.

[knepley]

The error is saying that label_int is not a DMLabel. Is it?

[gthyagi]

yes, that is true. label_int is not a DMLabel. However, petsc4py documentation says getlabel() returns DMLabel.
https://petsc.org/release/petsc4py/reference/petsc4py.PETSc.DM.html#petsc4py.PETSc.DM.getLabel:~:text=(str)-,Return%20type,%3A,-DMLabel

[knepley]

If a label of that name does not exist, I believe it returns None

[gthyagi]

it does exist.

cdef DMLabel label_int = <DMLabel> self.dm.getLabel("Internal")
print(label_int.view())

Label 'Internal':
[0]: 21 (5)
[0]: 22 (5)
[0]: 30 (5)
[0]: 36 (5)
[0]: 41 (5)