[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.