Unable to use discontinuous pressure field #277

bknight1 · 2024-12-09T05:25:06Z

UW3 version - c0b15d0

It appears that UW creates a discontinuous variable with degree = 0 (n = 242) whilst the one supplied in the script is with degree = 1 (n = 726), leading to a shape mismatch.

I can't see where in petsc_generic_snes_solvers.pyx this occurs, so any help would be appreciated.

Error message:

ValueError                                Traceback (most recent call last)
Cell In[5], line 10
      6 stokes.add_dirichlet_bc((0.0, sp.oo), "Left")
      7 stokes.add_dirichlet_bc((0.0, sp.oo), "Right")
---> 10 stokes.solve()

File [~/Documents/Research/GitHub/underworld3/src/underworld3/cython/petsc_generic_snes_solvers.pyx:2757](http://localhost:8888/lab/tree/Scripts/Stokes_models/mesh_based/~/Documents/Research/GitHub/underworld3/src/underworld3/cython/petsc_generic_snes_solvers.pyx#line=2756), in underworld3.cython.generic_solvers.SNES_Stokes_SaddlePt.solve()

File [~/Documents/Research/GitHub/underworld3/src/underworld3/cython/petsc_generic_snes_solvers.pyx:2762](http://localhost:8888/lab/tree/Scripts/Stokes_models/mesh_based/~/Documents/Research/GitHub/underworld3/src/underworld3/cython/petsc_generic_snes_solvers.pyx#line=2761), in underworld3.cython.generic_solvers.SNES_Stokes_SaddlePt.solve()

ValueError: could not broadcast input array from shape (242,) into shape (726,)

Code to reproduce:

# %%
import underworld3 as uw
import sympy as sp


# mesh = uw.meshing.StructuredQuadBox(elementRes=(3,) * 2)
mesh = uw.meshing.UnstructuredSimplexBox()
# -

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
)

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()
# -

The text was updated successfully, but these errors were encountered:

gthyagi · 2024-12-09T14:00:22Z

I identified the root cause of the issue.

Can someone clarify this?

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  3 offset  12

@knepley Why does the pressure field have a degree of freedom (dof) of 3 when the field is discontinuous?

Following info gets printed when I do this:

stokes.dm.getLocalSection().view()

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

knepley · 2024-12-09T14:14:57Z

That is discontinuous P1. I suspect the degree is provided as 1, and not 0.

gthyagi · 2024-12-09T14:50:40Z

@bknight1 change this line https://github.com/gthyagi/underworld3/blob/53dbf6d07cc7c56b4ce42d957059e34fe6b9da2a/src/underworld3/cython/petsc_generic_snes_solvers.pyx#L2724
to

if not unconstrained and self.Unknowns.p.continuous:

This will fix the array shape broadcasting issue

bknight1 · 2024-12-09T23:15:59Z

That is discontinuous P1. I suspect the degree is provided as 1, and not 0.

Yeah this is the case. I specify Deg=1 but Deg=0 is returned.

@bknight1 change this line https://github.com/gthyagi/underworld3/blob/53dbf6d07cc7c56b4ce42d957059e34fe6b9da2a/src/underworld3/cython/petsc_generic_snes_solvers.pyx#L2724 to
if not unconstrained and self.Unknowns.p.continuous:
This will fix the array shape broadcasting issue

I will try this now, thanks for figuring it out. Not sure what changes caused this issue as it was fine a couple of months ago.

gthyagi · 2024-12-10T03:53:27Z

The issue arose from the way I gathered the index set for the pressure field to separate it from the velocity field. I had assumed that the pressure field would always have a degree of freedom (DOF) of 1. However, when the field is discontinuous, all three points on a face are treated as DOFs for the pressure field when the degree is 1.

Now PETSC's section view makes much more sense to me. I better understand what PETSc's section structure is conveying.

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

Fixes #277

julesghub · 2024-12-10T23:55:20Z

Question @bknight1 and @gthyagi.
Can the fix be made without asking for the 'Unknowns'? Can we grab that information from the petsc field/vec only?

knepley · 2024-12-11T00:20:14Z

@julesghub What info do you want about the fields or layout?

julesghub · 2024-12-11T01:29:14Z

Hello @knepley ,

This code is for extracting the velocity and pressure values from the SNES vector after solving Stokes, copying the values into separate Underworld meshvariables with boundary values included.

The relevant code starts at this line:
https://github.com/gthyagi/underworld3/blob/53dbf6d07cc7c56b4ce42d957059e34fe6b9da2a/src/underworld3/cython/petsc_generic_snes_solvers.pyx#L2696

The algorithm follows these steps:

Retrieve the local vector (lvec) from the solved SNES.
Apply boundary conditions to the lvec. DMPlexSNESComputeBoundaryFEM
Decompose the lvec into velocity and pressure subvectors. (currently generating indexSets to perform this)
Copy these velocity and pressure arrays into Underworld's representation of mesh variables uw_v and uw_p - including the boundary conditions where appropriate.

Steps 3 and 4 were implemented by @gthyagi. Decomposing the lvec requires knowledge of the layout of the fields and dof constraints, which is were the original bug was.
If you could propose a cleaner solution to the steps above that would be awesome. Thanks Matt!

knepley · 2024-12-11T01:34:34Z

If you would like a global vector with the BC in it, we also have code to do that. This is what we use for output for things like VTK.
Yes, there is an easier way to do this. You can call DMCreateSubDM which returns a DM that makes the right size vectors for the field and also an IS which you can use for a call to VecISCopy. Thus I think you can do this in 2 or 3 lines.

gthyagi · 2024-12-11T01:38:01Z

Question @bknight1 and @gthyagi. Can the fix be made without asking for the 'Unknowns'? Can we grab that information from the petsc field/vec only?

@julesghub
This is alternative approach to get p_continuous

options = PETSc.Options()
p_continuous = options.getBool("private_{}_p_petscdualspace_lagrange_continuity".format(self.petsc_options_prefix))
    if not unconstrained and p_continuous:

julesghub · 2024-12-11T01:41:33Z

If you would like a global vector with the BC in it, we also have code to do that. This is what we use for output for things like VTK.

Yes, there is an easier way to do this. You can call DMCreateSubDM which returns a DM that makes the right size vectors for the field and also an IS which you can use for a call to VecISCopy. Thus I think you can do this in 2 or 3 lines.

Sounds great. @gthyagi let's go for this.

bknight1 · 2024-12-11T09:54:00Z

Question @bknight1 and @gthyagi. Can the fix be made without asking for the 'Unknowns'? Can we grab that information from the petsc field/vec only?

@julesghub This is alternative approach to get p_continuous
options = PETSc.Options()
p_continuous = options.getBool("private_{}_p_petscdualspace_lagrange_continuity".format(self.petsc_options_prefix))
    if not unconstrained and p_continuous:

p_continuous is already a variable, but is set from the Unknowns

Noting the updates for this also fixed the visualisation issue experienced in #236

gthyagi · 2024-12-13T04:23:27Z

@knepley I am observing a discrepancy in the index sets of the pressure field depending on whether I use dm.getLocalSection() or dm.createSubDM(). Do you have any suggestions on where things might be going wrong?

Here is the simple example i tried:

import underworld3 as uw
import sympy as sp

# mesh = uw.meshing.StructuredQuadBox(elementRes=(3,) * 2)
mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.9)

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
)

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()

output

0 SNES Function norm 1.131004041684e+05
    Residual norms for Solver_3_ solve.
    0 KSP Residual norm 7.479313126782e+05
    1 KSP Residual norm 5.773599984736e+00
  1 SNES Function norm 2.387153796436e-01

This is from dm.getLocalSection(), which is correct IS to separately fields from clvec.
This is what I am doing

underworld3/src/underworld3/cython/petsc_generic_snes_solvers.pyx

Lines 2696 to 2763 in 7c08857

    
                   cdef DM dm = self.dm 
        
                   cdef Vec clvec = self.dm.getLocalVec() 
        
                   self.dm.globalToLocal(gvec, clvec) 
        
                   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) 
        
                   # get index set of pressure and velocity to separate solution from localvec 
        
                   # get local section 
        
                   local_section = self.dm.getLocalSection() 
        
                   # Get the index sets for velocity and pressure fields 
        
                   # Field numbers (adjust based on your setup) 
        
                   velocity_field_num = 0 
        
                   pressure_field_num = 1 
        
                   # Function to get index set for a field 
        
                   def get_local_field_is(section, field, unconstrained=False): 
        
                       """ 
        
                       This function returns the index set of unconstrained points if True, or all points if False. 
        
                       """ 
        
                       pStart, pEnd = section.getChart() 
        
                       indices = [] 
        
                       for p in range(pStart, pEnd): 
        
                           dof = section.getFieldDof(p, field) 
        
                           if dof > 0: 
        
                               offset = section.getFieldOffset(p, field) 
        
                               if not unconstrained: 
        
                                   indices.append(offset) 
        
                               else: 
        
                                   cind = section.getFieldConstraintIndices(p, field) 
        
                                   constrained = set(cind) if cind is not None else set() 
        
                                   for i in range(dof): 
        
                                       if i not in constrained: 
        
                                           index = offset + i 
        
                                           indices.append(index) 
        
                       is_field = PETSc.IS().createGeneral(indices, comm=PETSc.COMM_SELF) 
        
                       return is_field 
        
                   # Get index sets for pressure (both constrained and unconstrained points) 
        
                   # we need indexset of pressure field to separate the solution from localvec. 
        
                   # so we don't care whether a point is constrained by bc or not 
        
                   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[:]

IS Object: 1 MPI process
  type: general
Number of indices in set 12
0 2
1 5
2 8
3 11
4 14
5 17
6 20
7 23
8 26
9 29
10 32
11 35

This is from dm.createSubDM(), which is incorrect

IS Object: 1 MPI process
  type: general
Number of indices in set 12
0 0
1 1
2 2
3 3
4 4
5 5
6 7
7 9
8 12
9 15
10 18
11 21

julesghub · 2024-12-13T05:55:33Z

@gthyagi, it's unclear where you have called createSubDM().
Also createSubDM takes input and returns a DM and IS objects, which can then be use on a Vec object.

As discussed I think a pure petsc4py script to demonstrate this would be useful.

gthyagi · 2024-12-16T03:47:50Z

@knepley, thank you for the suggestions—they worked! The issue was caused by passing the wrong dm to obtain the IS. It's all sorted out now.

gthyagi · 2024-12-16T03:50:14Z

@bknight1, please try this: #279. Could you also verify if this PR resolves #236? Thanks!

knepley · 2024-12-16T12:09:08Z

@gthyagi Yes! I will keep reading all the issues.

bknight1 · 2024-12-17T00:55:33Z

Getting the following error with #279

[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: Argument out of range
[0]PETSC ERROR: Scatter indices in ix are out of range: found [0,39422), expected in [0,24832)
[0]PETSC ERROR: WARNING! There are unused option(s) set! Could be the program crashed before usage or a spelling mistake, etc!

gthyagi · 2024-12-17T01:10:38Z

# %%
import underworld3 as uw
import sympy as sp


# mesh = uw.meshing.StructuredQuadBox(elementRes=(3,) * 2)
mesh = uw.meshing.UnstructuredSimplexBox()
# -

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
)

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()
# -

@bknight1 I am running this script without any error. Which one are you running?

bknight1 · 2024-12-17T01:14:23Z

I was running a different script. That one works, but if you add the following mesh variable it'll reproduce the error:

mat = uw.discretisation.MeshVariable("mat", mesh, 1, degree=0, continuous=False)

gthyagi · 2024-12-17T01:18:37Z

Can you paste the code snippet where you create all mesh variables? I think the order how you create mesh variables is important. That determines the field id.

bknight1 · 2024-12-17T01:19:51Z

I've tried creating the mat variable first and last in the code you posted and get the same result either way

gthyagi · 2024-12-17T01:26:58Z

Yesterday we were discussing about this case, let me investigate more on this.

gthyagi · 2024-12-17T11:10:04Z

@bknight1 I know the source of the issue. But not quite sure how to fix it. We will discuss in Thursday dev meeting.

julesghub · 2025-01-15T02:39:09Z

If you would like a global vector with the BC in it, we also have code to do that. This is what we use for output for things like VTK.

@knepley -happy new year! Could you point out the code that performs this - much appreciated.

Update petsc_generic_snes_solvers.pyx Fixes #277 to is not yet optimised. More work to come on that.

knepley · 2025-01-15T12:04:10Z

@julesghub Sure, it is https://gitlab.com/petsc/petsc/-/blob/main/src/dm/interface/dm.c?ref_type=heads#L6411

knepley · 2025-01-15T12:04:33Z

I use it in the output routines in plex.c and plexhdf5.c

bknight1 added a commit that referenced this issue Dec 10, 2024

Update petsc_generic_snes_solvers.pyx

b94ea7f

Fixes #277

bknight1 mentioned this issue Dec 10, 2024

fix pressure variable #278

Merged

julesghub pushed a commit that referenced this issue Jan 15, 2025

fix pressure variable (#278)

1256716

Update petsc_generic_snes_solvers.pyx Fixes #277 to is not yet optimised. More work to come on that.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Unable to use discontinuous pressure field #277

Unable to use discontinuous pressure field #277

bknight1 commented Dec 9, 2024

gthyagi commented Dec 9, 2024 •

edited

Loading

knepley commented Dec 9, 2024

gthyagi commented Dec 9, 2024

bknight1 commented Dec 9, 2024

gthyagi commented Dec 10, 2024

julesghub commented Dec 10, 2024

knepley commented Dec 11, 2024

julesghub commented Dec 11, 2024

knepley commented Dec 11, 2024

gthyagi commented Dec 11, 2024 •

edited

Loading

julesghub commented Dec 11, 2024

bknight1 commented Dec 11, 2024

gthyagi commented Dec 13, 2024 •

edited

Loading

julesghub commented Dec 13, 2024

gthyagi commented Dec 16, 2024

gthyagi commented Dec 16, 2024

knepley commented Dec 16, 2024

bknight1 commented Dec 17, 2024

gthyagi commented Dec 17, 2024

bknight1 commented Dec 17, 2024

gthyagi commented Dec 17, 2024

bknight1 commented Dec 17, 2024 •

edited

Loading

gthyagi commented Dec 17, 2024

gthyagi commented Dec 17, 2024

julesghub commented Jan 15, 2025

knepley commented Jan 15, 2025

knepley commented Jan 15, 2025

Unable to use discontinuous pressure field #277

Unable to use discontinuous pressure field #277

Comments

bknight1 commented Dec 9, 2024

gthyagi commented Dec 9, 2024 • edited Loading

knepley commented Dec 9, 2024

gthyagi commented Dec 9, 2024

bknight1 commented Dec 9, 2024

gthyagi commented Dec 10, 2024

julesghub commented Dec 10, 2024

knepley commented Dec 11, 2024

julesghub commented Dec 11, 2024

knepley commented Dec 11, 2024

gthyagi commented Dec 11, 2024 • edited Loading

julesghub commented Dec 11, 2024

bknight1 commented Dec 11, 2024

gthyagi commented Dec 13, 2024 • edited Loading

julesghub commented Dec 13, 2024

gthyagi commented Dec 16, 2024

gthyagi commented Dec 16, 2024

knepley commented Dec 16, 2024

bknight1 commented Dec 17, 2024

gthyagi commented Dec 17, 2024

bknight1 commented Dec 17, 2024

gthyagi commented Dec 17, 2024

bknight1 commented Dec 17, 2024 • edited Loading

gthyagi commented Dec 17, 2024

gthyagi commented Dec 17, 2024

julesghub commented Jan 15, 2025

knepley commented Jan 15, 2025

knepley commented Jan 15, 2025

gthyagi commented Dec 9, 2024 •

edited

Loading

gthyagi commented Dec 11, 2024 •

edited

Loading

gthyagi commented Dec 13, 2024 •

edited

Loading

bknight1 commented Dec 17, 2024 •

edited

Loading