[BUG] - Order of MeshVariable declaration becomes important when adding B.C. using add_condition()

[jcgraciosa]

Branch tested: development
PETSc version: 3.22.0
Issue description:

When adding boundary conditions with the add_condition function, the user must declare the MeshVariables used in the solver first in order to ensure that the field numbers forwarded to PETSc are correct. Otherwise, an error occurs. As a minimal example:

import underworld3 as uw

replicate_error = True
resolution = 16
qdeg = 3
Vdeg = 2
Pdeg = Vdeg - 1
Pcont = True

height  = 1
width   = 1

minX, maxX = 0, width
minY, maxY = 0, height

meshbox = uw.meshing.StructuredQuadBox( minCoords=(minX, minY), maxCoords=(maxX, maxY), elementRes = (resolution, resolution), qdegree = qdeg)

if replicate_error: # this will cause an error since v_ana is declared first (w/c is not used in solver)
    v_ana = uw.discretisation.MeshVariable("Ua", meshbox, meshbox.dim, degree=Vdeg)
  
    v_soln = uw.discretisation.MeshVariable("U", meshbox, meshbox.dim, degree=Vdeg)
    p_soln = uw.discretisation.MeshVariable("P", meshbox, 1, degree = Pdeg, continuous = Pcont)
    p_ana = uw.discretisation.MeshVariable("Pa", meshbox, 1, degree = Pdeg, continuous = Pcont)
else: # this will have no error since v_soln and p_soln are declated first 
    v_soln = uw.discretisation.MeshVariable("U", meshbox, meshbox.dim, degree=Vdeg)
    p_soln = uw.discretisation.MeshVariable("P", meshbox, 1, degree = Pdeg, continuous = Pcont)
    v_ana = uw.discretisation.MeshVariable("Ua", meshbox, meshbox.dim, degree=Vdeg)
    p_ana = uw.discretisation.MeshVariable("Pa", meshbox, 1, degree = Pdeg, continuous = Pcont)
  
stokes = uw.systems.Stokes(
                            mesh            = meshbox, 
                            velocityField   = v_soln,
                            pressureField   = p_soln,
                            solver_name     ="stokes",
                            )

stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
stokes.constitutive_model.Parameters.viscosity = 1

stokes.add_dirichlet_bc((1., 0.), "Bottom")
stokes.add_condition(p_soln.field_id, "dirichlet", [(1.)], "Left", components = (0))

stokes.solve(zero_init_guess = True)

The PETSc error message during the solve() is:

[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR:   It appears a new error in the code was triggered after a previous error, possibly because:
[0]PETSC ERROR:   -  The first error was not properly handled via (for example) the use of
[0]PETSC ERROR:      PetscCall(TheFunctionThatErrors()); or
[0]PETSC ERROR:   -  The second error was triggered while handling the first error.
[0]PETSC ERROR:   Above is the traceback for the previous unhandled error, below the traceback for the next error
[0]PETSC ERROR:   ALL ERRORS in the PETSc libraries are fatal, you should add the appropriate error checking to the code
[0]PETSC ERROR: Argument out of range
[0]PETSC ERROR: Field number 2 must be in [0, 2)
[0]PETSC ERROR: WARNING! There are unused option(s) set! Could be the program crashed before usage or a spelling mistake, etc!

Is this expected behaviour?
I guess this will most likely be encountered by users working with more complicated models.

[lmoresi]

It is expected because the field id that you pass in to this function has nothing to do with the one that is used to store the data. It is the internal notation of the solver and really should never be exposed to the user at all. We should rewrite this because it causes confusion and I don't think it is the correct approach.

…

Sent from my iPhone On 14 Oct 2024, at 04:07, jcgraciosa ***@***.***> wrote:  Branch tested: development PETSc version: 3.22.0 Issue description: When adding boundary conditions with the add_condition function, the user must declare the MeshVariables used in the solver first in order to ensure that the field numbers forwarded to PETSc are correct. Otherwise, an error occurs. As a minimal example: import underworld3 as uw replicate_error = True resolution = 16 qdeg = 3 Vdeg = 2 Pdeg = Vdeg - 1 Pcont = True height = 1 width = 1 minX, maxX = 0, width minY, maxY = 0, height meshbox = uw.meshing.StructuredQuadBox( minCoords=(minX, minY), maxCoords=(maxX, maxY), elementRes = (resolution, resolution), qdegree = qdeg) if replicate_error: # this will cause an error since v_ana is declared first (w/c is not used in solver) v_ana = uw.discretisation.MeshVariable("Ua", meshbox, meshbox.dim, degree=Vdeg) v_soln = uw.discretisation.MeshVariable("U", meshbox, meshbox.dim, degree=Vdeg) p_soln = uw.discretisation.MeshVariable("P", meshbox, 1, degree = Pdeg, continuous = Pcont) p_ana = uw.discretisation.MeshVariable("Pa", meshbox, 1, degree = Pdeg, continuous = Pcont) else: # this will have no error since v_soln and p_soln are declated first v_soln = uw.discretisation.MeshVariable("U", meshbox, meshbox.dim, degree=Vdeg) p_soln = uw.discretisation.MeshVariable("P", meshbox, 1, degree = Pdeg, continuous = Pcont) v_ana = uw.discretisation.MeshVariable("Ua", meshbox, meshbox.dim, degree=Vdeg) p_ana = uw.discretisation.MeshVariable("Pa", meshbox, 1, degree = Pdeg, continuous = Pcont) stokes = uw.systems.Stokes( mesh = meshbox, velocityField = v_soln, pressureField = p_soln, solver_name ="stokes", ) stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel stokes.constitutive_model.Parameters.viscosity = 1 stokes.add_dirichlet_bc((1., 0.), "Bottom") stokes.add_condition(p_soln.field_id, "dirichlet", [(1.)], "Left", components = (0)) stokes.solve(zero_init_guess = True) The PETSc error message during the solve() is: [0]PETSC ERROR: --------------------- Error Message -------------------------------------------------------------- [0]PETSC ERROR: It appears a new error in the code was triggered after a previous error, possibly because: [0]PETSC ERROR: - The first error was not properly handled via (for example) the use of [0]PETSC ERROR: PetscCall(TheFunctionThatErrors()); or [0]PETSC ERROR: - The second error was triggered while handling the first error. [0]PETSC ERROR: Above is the traceback for the previous unhandled error, below the traceback for the next error [0]PETSC ERROR: ALL ERRORS in the PETSc libraries are fatal, you should add the appropriate error checking to the code [0]PETSC ERROR: Argument out of range [0]PETSC ERROR: Field number 2 must be in [0, 2) [0]PETSC ERROR: WARNING! There are unused option(s) set! Could be the program crashed before usage or a spelling mistake, etc! Is this expected behaviour? I guess this will most likely be encountered by users working with more complicated models. — Reply to this email directly, view it on GitHub<#253>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ADABPI424VELXU7MWZVSEYTZ3OQXXAVCNFSM6AAAAABP4XZWXSVHI2DSMVQWIX3LMV43ASLTON2WKOZSGU4DKNRYGM4DMOI>. You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[gthyagi]

@jcgraciosa, adding onto @lmoresi's comment, the Stokes solver takes the mesh DM and copies it to create a new DM that contains only the velocity and pressure variables. The velocity field ID is always set to 0, and the pressure field ID is set to 1.

Here is the code:

underworld3/src/underworld3/cython/petsc_generic_snes_solvers.pyx

Lines 2318 to 2338 in d1c7e0e

	if self.dm.getNumFields() == 0:
	if self.verbose:
	print(f"{uw.mpi.rank}: Building FE / quadrature for {self.name}", flush=True)
	options = PETSc.Options()
	options.setValue("private_{}_u_petscspace_degree".format(self.petsc_options_prefix), u_degree) # for private variables
	self.petsc_fe_u = PETSc.FE().createDefault(mesh.dim, mesh.dim, mesh.isSimplex, mesh.qdegree, "private_{}_u_".format(self.petsc_options_prefix), PETSc.COMM_SELF)
	self.petsc_fe_u.setName("velocity")
	self.petsc_fe_u_id = self.dm.getNumFields()
	self.dm.setField( self.petsc_fe_u_id, self.petsc_fe_u )
	options.setValue("private_{}_p_petscspace_degree".format(self.petsc_options_prefix), p_degree)
	options.setValue("private_{}_p_petscdualspace_lagrange_continuity".format(self.petsc_options_prefix), p_continous)
	options.setValue("private_{}_p_petscdualspace_lagrange_node_endpoints".format(self.petsc_options_prefix), False)
	self.petsc_fe_p = PETSc.FE().createDefault(mesh.dim, 1, mesh.isSimplex, mesh.qdegree, "private_{}_p_".format(self.petsc_options_prefix), PETSc.COMM_SELF)
	self.petsc_fe_p.setName("pressure")
	self.petsc_fe_p_id = self.dm.getNumFields()
	self.dm.setField( self.petsc_fe_p_id, self.petsc_fe_p)

[jcgraciosa]

@lmoresi and @gthyagi, thanks for the response.

I think in this case, if the function stays, we can add the information regarding the field ids used in the solver's internal variables in its Docstring (e.g. something like: "For Stokes and Navier Stokes solvers: velocity - 0, pressure - 1.").