createSubDM() returns inconsistent IS in petsc4py

[gthyagi]

Hi @knepley,

I've noticed that createSubDM() in petsc4py seems to return an incorrect IS. Could it be that I'm misunderstanding how createSubDM() is intended to work? Would you mind taking a look and helping us out? Thank you!

Here is the script:

import underworld3 as uw
import sympy as sp

mesh = uw.meshing.UnstructuredSimplexBox(cellSize=1.)
x, y = mesh.X

u = uw.discretisation.MeshVariable(r"mathbf{u}", mesh, mesh.dim, vtype=uw.VarType.VECTOR, degree=2)
p = uw.discretisation.MeshVariable(r"mathbf{p}", mesh, 1, vtype=uw.VarType.SCALAR, degree=1, continuous=False)
t = uw.discretisation.MeshVariable(r"mathbf{t}", mesh, 1, vtype=uw.VarType.SCALAR, degree=1, continuous=False)

stokes = uw.systems.Stokes(mesh, velocityField=u, pressureField=p)
stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
stokes.constitutive_model.Parameters.shear_viscosity_0 = 1

stokes.petsc_options["snes_monitor"] = None
stokes.petsc_options["ksp_monitor"] = None

stokes.bodyforce = 1.0e6 * sp.Matrix([0, x])

stokes.add_dirichlet_bc((0.0, 0.0), "Bottom")
stokes.add_dirichlet_bc((0.0, 0.0), "Top")

stokes.add_dirichlet_bc((0.0, sp.oo), "Left")
stokes.add_dirichlet_bc((0.0, sp.oo), "Right")

stokes.solve()

Here is the output from the section view:

stokes.dm.getLocalSection().view()

PetscSection Object: 1 MPI process
  type not yet set
2 fields
  field 0 "velocity" with 2 components
Process 0:
  (   0) dof  0 offset   0
  (   1) dof  0 offset   3
  (   2) dof  0 offset   6
  (   3) dof  0 offset   9
  (   4) dof  2 offset  12 constrained 0 1
  (   5) dof  2 offset  14 constrained 0 1
  (   6) dof  2 offset  16 constrained 0 1
  (   7) dof  2 offset  18 constrained 0 1
  (   8) dof  2 offset  20
  (   9) dof  2 offset  22 constrained 0 1
  (  10) dof  2 offset  24
  (  11) dof  2 offset  26
  (  12) dof  2 offset  28 constrained 0
  (  13) dof  2 offset  30
  (  14) dof  2 offset  32 constrained 0
  (  15) dof  2 offset  34
  (  16) dof  2 offset  36 constrained 0 1
  field 1 "pressure" with 1 components
Process 0:
  (   0) dof  3 offset   0
  (   1) dof  3 offset   3
  (   2) dof  3 offset   6
  (   3) dof  3 offset   9
  (   4) dof  0 offset  14
  (   5) dof  0 offset  16
  (   6) dof  0 offset  18
  (   7) dof  0 offset  20
  (   8) dof  0 offset  22
  (   9) dof  0 offset  24
  (  10) dof  0 offset  26
  (  11) dof  0 offset  28
  (  12) dof  0 offset  30
  (  13) dof  0 offset  32
  (  14) dof  0 offset  34
  (  15) dof  0 offset  36
  (  16) dof  0 offset  38
  PetscSectionSym Object: 1 MPI process
    type: label
    Label 'depth'
    Symmetry for stratum value 0 (0 dofs per point): no symmetries
    Symmetry for stratum value 1 (0 dofs per point): no symmetries
    Symmetry for stratum value 2 (3 dofs per point):
      Orientation range: [-3, 3)
    Symmetry for stratum value -1 (0 dofs per point): no symmetries

I got following IS for the pressure field, which is correct and I can crosscheck from the petsc section info:

stokes.dm.createSubDM(1)[0].array

array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11], dtype=int32)

For the velocity field, I got this:

stokes.dm.createSubDM(0)[0].array

array([12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23], dtype=int32)

Isn't the above line supposed to return following IS?

array([12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37], dtype=int32)

[gthyagi]

@knepley petsc documentation says this
https://petsc.org/release/manualpages/DM/DMCreateFieldDecomposition/#:~:text=Each%20IS%20contains%20the%20global%20indices%20of%20the%20dofs%20of%20the%20corresponding%20field%2C%20defined%20by%20DMAddField().%20The%20optional%20list%20of%20DMs%20define%20the%20DM%20for%20each%20subproblem.

I realize now that my understanding of createSubDM() was incorrect. With that in mind, my next question is: how can I retrieve all the indices for a given field in a DM?

[bknight1]

In the earlier version of the code it looks like we originally only got the pressure indices and assume the rest are velocity

underworld3/src/underworld3/cython/petsc_generic_snes_solvers.pyx

Lines 2739 to 2762 in 92a94f5

	pressure_is = get_local_field_is(local_section, pressure_field_num)
	# Get the total number of entries in the local vector
	size = self.dm.getLocalVec().getLocalSize()
	# Create a list of all indices
	all_indices = set(range(size))
	# Get indices of the pressure field
	pressure_indices = set(pressure_is.getIndices())
	# Compute the complement for the velocity field
	velocity_indices = sorted(list(all_indices - pressure_indices))
	# Create the index set for velocity
	velocity_is = PETSc.IS().createGeneral(velocity_indices, comm=PETSc.COMM_SELF)
	# Copy solution back into pressure and velocity variables
	with self.mesh.access(self.Unknowns.p, self.Unknowns.u):
	for name, var in self.fields.items():
	if name=='velocity':
	var.vec.array[:] = clvec.getSubVector(velocity_is).array[:]
	elif name=='pressure':
	var.vec.array[:] = clvec.getSubVector(pressure_is).array[:]

I've written an update to the code using that approach and assuming the pressure field is from the first field, i.e. self.dm.createSubDM(1)[0]:

        # pressure_is = get_local_field_is(local_section, pressure_field_num)
        # pressure_is = PETSc.IS().createGeneral(pressure_is, comm=PETSc.COMM_SELF)

        ### pressure field is field 1
        pressure_is = self.dm.createSubDM(1)[0]

        # Get the total number of entries in the local vector
        size = self.dm.getLocalVec().getLocalSize()

        # Create a list of all indices
        all_indices = set( range(size) )

        # Get indices of the pressure field
        pressure_indices = set( pressure_is.getIndices() )

        # Compute the complement for the velocity field
        velocity_indices = sorted(list(all_indices - pressure_indices))

        # Create the index set for velocity
        velocity_is = PETSc.IS().createGeneral(velocity_indices, comm=PETSc.COMM_SELF)

        # Copy solution back into pressure and velocity variables
        with self.mesh.access(self.Unknowns.p, self.Unknowns.u):
             for name, var in self.fields.items():
                 if name=='velocity':
                     var.vec.array[:] = clvec.getSubVector(velocity_is).array[:]
                 elif name=='pressure':
                     var.vec.array[:] = clvec.getSubVector(pressure_is).array[:]

[gthyagi]

Feature Request: Direct Function to Obtain Both Constrained and Unconstrained Field Indices in PETSc/petsc4py

Currently, there is no single PETSc or petsc4py function that returns both the constrained and unconstrained DOF indices for a given field in a DM. We can use DMCreateSubDM() (or DMCreateFieldIS()) to retrieve unconstrained global indices easily, but obtaining constrained indices requires manually looping over the local PetscSection and checking each point’s constrained DOFs.

Use Case:
An example scenario is separating the velocity and pressure components from the local vector after boundary values have been inserted. Having a built-in routine to fetch both types of indices would simplify this process significantly.

code:

cdef DM dm = stokes.dm
cdef Vec clvec = stokes.dm.getLocalVec()
stokes.dm.globalToLocal(gvec, clvec) # gvec is the global vector
ierr = DMPlexSNESComputeBoundaryFEM(dm.dm, <void*>clvec.vec, NULL); CHKERRQ(ierr)

if verbose and uw.mpi.rank == 0:
         print(f"SNES Compute Boundary FEM Successfull", flush=True)

# now use indice set of field with all DOFs to separately each field.

[knepley]

First I need a definition of your words. I understand the words "local" and "global". Local indices are those indexing the serial local vector defined by the local section. Global indices index the parallel global vector defined by the global section. Local vectors contain all dofs specified in the problem, some redundantly. Global vectors omit the constrained dofs and are not redundant.

I do not understand what you are asking for with "constrained and unconstrained DOF indices".

[knepley]

Maybe we start with "What kind of vector do you want to index?"

[gthyagi]

@knepley We are planning to run a few additional tests to help debug the issue. If these tests also fail, I will provide detailed information about my specific request. Thanks