Merge branch 'development' of github.com:underworldcode/underworld3 i…

…nto development

src/underworld3/cython/petsc_generic_snes_solvers.pyx

@@ -28,8 +28,8 @@ class SolverBaseClass(uw_object):
 
     def __init__(self, mesh):
 
-        super().__init__()
 
+        super().__init__()
 
         self.mesh = mesh
         self.mesh_dm_coordinate_hash = None
@@ -39,9 +39,13 @@ class SolverBaseClass(uw_object):
 
         self._order = 0
         self._constitutive_model = None
-
         self._rebuild_after_mesh_update = self._build
 
+        self.name = "Solver_{}_".format(self.instance_number)
+        self.petsc_options_prefix = self.name
+        self.petsc_options = PETSc.Options(self.petsc_options_prefix)
+
+
         return
 
     class _Unknowns:
@@ -132,7 +136,6 @@ class SolverBaseClass(uw_object):
 
         display(Markdown(fr"This solver is formulated in {self.mesh.dim} dimensions"))
 
-
         return
 
     def _reset(self):
@@ -172,6 +175,13 @@ class SolverBaseClass(uw_object):
                             # can insert new functions in template (surface integrals problematic in
                             # the current implementation )
 
+        if len(self.natural_bcs) > 0:
+            if not "Null_Boundary" in self.natural_bcs:
+                bc = (0,)*self.Unknowns.u.shape[1]
+                print(f"Adding Natural BC {bc }- {self.Unknowns.u.shape[1] }", flush=True)
+                self.add_natural_bc(bc, "Null_Boundary")
+
+
 
         self._setup_pointwise_functions(verbose, debug=debug, debug_name=debug_name)
         self._setup_discretisation(verbose)
@@ -400,7 +410,6 @@ class SNES_Scalar(SolverBaseClass):
                  mesh     : uw.discretisation.Mesh,
                  u_Field  : uw.discretisation.MeshVariable = None,
                  degree: int = 2,
-                 solver_name: str = "",
                  verbose    = False,
                  DuDt          : Union[uw.systems.ddt.SemiLagrangian, uw.systems.ddt.Lagrangian] = None,
                  DFDt          : Union[uw.systems.ddt.SemiLagrangian, uw.systems.ddt.Lagrangian] = None,
@@ -420,20 +429,9 @@ class SNES_Scalar(SolverBaseClass):
         self.Unknowns.DuDt = DuDt
         self.Unknowns.DFDt = DFDt
 
-        self.name = solver_name
         self.verbose = verbose
         self._tolerance = 1.0e-4
 
-        ## Todo: this is obviously not particularly robust
-        if solver_name != "" and not solver_name.endswith("_"):
-            self.petsc_options_prefix = solver_name+"_"
-        else:
-            self.petsc_options_prefix = solver_name
-
-        options = PETSc.Options()
-        self.petsc_options = PETSc.Options(self.petsc_options_prefix)
-
-
         # Here we can set some defaults for this set of KSP / SNES solvers
 
         ## FAST as possible for simple problems:
@@ -543,6 +541,9 @@ class SNES_Scalar(SolverBaseClass):
         options.setValue("private_{}_petscdualspace_lagrange_continuity".format(self.petsc_options_prefix), self.u.continuous)
         options.setValue("private_{}_petscdualspace_lagrange_node_endpoints".format(self.petsc_options_prefix), False)
 
+        # Should del these when finished
+
+
         self.petsc_fe_u = PETSc.FE().createDefault(mesh.dim, 1, mesh.isSimplex, mesh.qdegree, "private_{}_".format(self.petsc_options_prefix), PETSc.COMM_SELF,)
         self.petsc_fe_u_id = self.dm.getNumFields()
         self.dm.setField( self.petsc_fe_u_id, self.petsc_fe_u )
@@ -857,21 +858,6 @@ class SNES_Scalar(SolverBaseClass):
 
         self._build(verbose, debug, debug_name)
 
-        # if (not self.is_setup):
-        #     if self.dm is not None:
-        #         self.dm.destroy()
-        #         self.dm = None  # Should be able to avoid nuking this if we
-        #                     # can insert new functions in template (surface integrals problematic in
-        #                     # the current implementation )
-
-        #     self._setup_pointwise_functions(verbose, debug=debug, debug_name=debug_name)
-        #     self._setup_discretisation(verbose)
-        #     self._setup_solver(verbose)
-        # else:
-        #     # If the mesh has changed, this will rebuild (and do nothing if unchanged)
-        #     self._setup_discretisation(verbose)
-
-
         gvec = self.dm.getGlobalVec()
 
         if not zero_init_guess:
@@ -991,14 +977,12 @@ class SNES_Vector(SolverBaseClass):
                  mesh     : uw.discretisation.Mesh,
                  u_Field  : uw.discretisation.MeshVariable = None,
                  degree     = 2,
-                 solver_name: str = "",
                  verbose    = False,
                  DuDt          : Union[uw.systems.ddt.SemiLagrangian, uw.systems.ddt.Lagrangian] = None,
                  DFDt          : Union[uw.systems.ddt.SemiLagrangian, uw.systems.ddt.Lagrangian] = None,
                  ):
 
 
-
         super().__init__(mesh)
 
         self.Unknowns.u = u_Field
@@ -1011,22 +995,14 @@ class SNES_Vector(SolverBaseClass):
 
         ## Keep track
 
-        self.name = solver_name
         self.verbose = verbose
         self._tolerance = 1.0e-4
 
         ## Todo: this is obviously not particularly robust
 
-        if solver_name != "" and not solver_name.endswith("_"):
-            self.petsc_options_prefix = solver_name+"_"
-        else:
-            self.petsc_options_prefix = solver_name
-
-        options = PETSc.Options()
+        # options = PETSc.Options()
         # options["dm_adaptor"]= "pragmatic"
 
-        self.petsc_options = PETSc.Options(self.petsc_options_prefix)
-
         # Here we can set some defaults for this set of KSP / SNES solvers
         self.petsc_options["snes_type"] = "newtonls"
         self.petsc_options["ksp_rtol"] = 1.0e-3
@@ -1690,7 +1666,6 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
                  pressureField : Optional[underworld3.discretisation.MeshVariable] = None,
                  degree        : Optional[int] = 2,
                  p_continuous  : Optional[bool] = True,
-                 solver_name   : Optional[str]                           ="stokes_pt_",
                  verbose       : Optional[bool]                           =False,
                  DuDt          : Union[uw.systems.ddt.SemiLagrangian, uw.systems.ddt.Lagrangian] = None,
                  DFDt          : Union[uw.systems.ddt.SemiLagrangian, uw.systems.ddt.Lagrangian] = None,
@@ -1699,7 +1674,6 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
 
         super().__init__(mesh)
 
-        self.name = solver_name
         self.verbose = verbose
         self.dm = None
 
@@ -1723,16 +1697,8 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
 
         # I expect the following to break for anyone who wants to name their solver _stokes__ etc etc (LM)
 
-        if solver_name != "" and not solver_name.endswith("_"):
-            self.petsc_options_prefix = solver_name+"_"
-        else:
-            self.petsc_options_prefix = solver_name
-
         # options = PETSc.Options()
         # options["dm_adaptor"]= "pragmatic"
-
-        self.petsc_options = PETSc.Options(self.petsc_options_prefix)
-
         # Here we can set some defaults for this set of KSP / SNES solvers
 
         if self.verbose == True:

src/underworld3/discretisation.py

@@ -1,4 +1,6 @@
 from typing import Optional, Tuple, Union
+from enum import Enum
+
 import os
 import numpy
 import sympy
@@ -189,26 +191,35 @@ def __init__(
                     % (plex_or_meshfile, ext[1:])
                 )
 
+        if boundaries is None:
+            class replacement_boundaries(Enum):
+                Null_Boundary = 666
+                All_Boundaries = 1001
+
+            boundaries = replacement_boundaries
+
         self.filename = filename
         self.boundaries = boundaries
         self.boundary_normals = boundary_normals
 
+
+
         # options.delValue("dm_plex_gmsh_mark_vertices")
         # options.delValue("dm_plex_gmsh_multiple_tags")
         # options.delValue("dm_plex_gmsh_use_regions")
         self.dm.setFromOptions()
 
         # uw.adaptivity._dm_stack_bcs(self.dm, self.boundaries, "UW_Boundaries")
 
-        # all_edges_label_dm = self.dm.getLabel("depth")
-        # if all_edges_label_dm:
-        #     all_edges_IS_dm = all_edges_label_dm.getStratumIS(1)
-        #     # all_edges_IS_dm.view()
+        all_edges_label_dm = self.dm.getLabel("depth")
+        if all_edges_label_dm:
+            all_edges_IS_dm = all_edges_label_dm.getStratumIS(0)
+            # all_edges_IS_dm.view()
 
-        # self.dm.createLabel("All_Edges")
-        # all_edges_label = self.dm.getLabel("All_Edges")
-        # if all_edges_label and all_edges_IS_dm:
-        #     all_edges_label.setStratumIS(boundaries.All_Edges.value, all_edges_IS_dm)
+        self.dm.createLabel("Null_Boundary")
+        all_edges_label = self.dm.getLabel("Null_Boundary")
+        if all_edges_label and all_edges_IS_dm:
+           all_edges_label.setStratumIS(boundaries.Null_Boundary.value, all_edges_IS_dm)
 
         ## --- UW_Boundaries label
         if self.boundaries is not None:
@@ -230,6 +241,8 @@ def __init__(
                     if label_is:
                         stacked_bc_label.setStratumIS(b.value, label_is)
 
+
+
             uw.mpi.barrier()
 
         ## ---

src/underworld3/meshing.py

@@ -46,12 +46,14 @@ def UnstructuredSimplexBox(
     """
 
     class boundaries_2D(Enum):
+        Null_Boundary = 666
         Bottom = 11
         Top = 12
         Right = 13
         Left = 14
 
     class boundaries_3D(Enum):
+        Null_Boundary = 666
         Bottom = 11
         Top = 12
         Right = 13
@@ -281,12 +283,14 @@ def StructuredQuadBox(
     # boundaries = {"Bottom": 1, "Top": 2, "Right": 3, "Left": 4, "Front": 5, "Back": 6}
 
     class boundaries_2D(Enum):
+        Null_Boundary = 666
         Bottom = 11
         Top = 12
         Right = 13
         Left = 14
 
     class boundaries_3D(Enum):
+        Null_Boundary = 666
         Bottom = 11
         Top = 12
         Right = 13
@@ -576,6 +580,7 @@ def SphericalShell(
     verbose=False,
 ):
     class boundaries(Enum):
+        Null_Boundary = 666
         Lower = 11
         Upper = 12
         Centre = 1
@@ -793,6 +798,7 @@ def SphericalShellInternalBoundary(
     """
 
     class boundaries(Enum):
+        Null_Boundary = 666
         Centre = 1
         Lower = 11
         Internal = 12
@@ -1035,6 +1041,7 @@ def SegmentofSphere(
     """
 
     class boundaries(Enum):
+        Null_Boundary = 666
         Lower = 11
         Upper = 12
         East = 13
@@ -1256,6 +1263,7 @@ def QuarterAnnulus(
     verbose=False,
 ):
     class boundaries(Enum):
+        Null_Boundary = 666
         Lower = 1
         Upper = 2
         Left = 3
@@ -1426,6 +1434,7 @@ def Annulus(
     verbose=False,
 ):
     class boundaries(Enum):
+        Null_Boundary = 666
         Lower = 1
         Upper = 2
         Centre = 10
@@ -1661,6 +1670,7 @@ def SegmentofAnnulus(
     """
 
     class boundaries(Enum):
+        Null_Boundary = 666
         Lower = 1
         Upper = 2
         Left = 3
@@ -1860,6 +1870,7 @@ def AnnulusWithSpokes(
     verbose=False,
 ):
     class boundaries(Enum):
+        Null_Boundary = 666
         Lower = 10
         LowerPlus = 11
         Upper = 20
@@ -2124,7 +2135,7 @@ class boundaries(Enum):
         Internal = 2
         Upper = 3
         Centre = 10
-        All_Edges = 1000
+        Null_Boundary = 666
 
     if cellSize_Inner is None:
         cellSize_Inner = cellSize
@@ -2338,6 +2349,8 @@ class boundaries(Enum):
         Internal = 2
         Upper = 3
         Centre = 10
+        Null_Boundary = 666
+
 
     if cellSize_Lower is None:
         cellSize_Lower = cellSize
@@ -2547,6 +2560,8 @@ def CubedSphere(
     class boundaries(Enum):
         Lower = 1
         Upper = 2
+        Null_Boundary = 666
+
 
     r1 = radiusInner / np.sqrt(3)
     r2 = radiusOuter / np.sqrt(3)
@@ -2767,6 +2782,8 @@ class boundaries(Enum):
         South = 4
         East = 5
         West = 6
+        Null_Boundary = 666
+
 
     r1 = radiusInner / np.sqrt(3)
     r2 = radiusOuter / np.sqrt(3)
@@ -3149,6 +3166,8 @@ class boundaries(Enum):
         UpperPlus = 31
         Centre = 1
         Slices = 40
+        Null_Boundary = 666
+
 
     meshRes = cellSize
     num_segments = numSegments
@@ -3552,6 +3571,8 @@ class boundaries(Enum):
         UpperPlus = 31
         Centre = 1
         Slices = 40
+        Null_Boundary = 666
+
 
     meshRes = cellSize
     num_segments = numSegments

src/underworld3/swarm.py

@@ -631,7 +631,7 @@ def sym(self):
     def sym_1d(self):
         return self._MaskArray
 
-    # We could also add a __getitem__ call to access each mask
+    # We can  also add a __getitem__ call to access each mask
 
     def __getitem__(self, index):
         return self.sym[index]