diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index d78a5b0f9..7a285fd16 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/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: diff --git a/src/underworld3/discretisation.py b/src/underworld3/discretisation.py index 1a3d5886b..5f784eeae 100644 --- a/src/underworld3/discretisation.py +++ b/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,10 +191,19 @@ 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") @@ -200,15 +211,15 @@ def __init__( # 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() ## --- diff --git a/src/underworld3/meshing.py b/src/underworld3/meshing.py index 7d3814fa8..1025e7a08 100644 --- a/src/underworld3/meshing.py +++ b/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 diff --git a/src/underworld3/swarm.py b/src/underworld3/swarm.py index d3f3596b0..48f24fb88 100644 --- a/src/underworld3/swarm.py +++ b/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] diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index 7d63139d6..96c6d2c11 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -58,7 +58,6 @@ def __init__( self, mesh: uw.discretisation.Mesh, u_Field: uw.discretisation.MeshVariable = None, - solver_name: str = "", verbose=False, degree=2, DuDt: Union[SemiLagrangian_DDt, Lagrangian_DDt] = None, @@ -71,15 +70,11 @@ def __init__( mesh, u_Field, degree, - solver_name, verbose, DuDt=DuDt, DFDt=DFDt, ) - if solver_name == "": - solver_name = "Poisson_{}_".format(self.instance_number) - # default values for properties self.f = sympy.Matrix.zeros(1, 1) @@ -187,7 +182,6 @@ def __init__( h_Field: Optional[uw.discretisation.MeshVariable] = None, v_Field: Optional[uw.discretisation.MeshVariable] = None, degree: int = 2, - solver_name: str = "", verbose=False, DuDt=None, DFDt=None, @@ -197,15 +191,11 @@ def __init__( mesh, h_Field, degree, - solver_name, verbose, DuDt=DuDt, DFDt=DFDt, ) - if solver_name == "": - self.solver_name = "Darcy_{}_".format(self.instance_number) - # default values for properties self._f = sympy.Matrix([0]) self._k = 1 @@ -223,8 +213,8 @@ def __init__( self._v = self._v_projector.Unknowns.u - # If we add smoothing, it should be small relative to actual diffusion (self.viscosity) - self._v_projector.smoothing = 0 + # If we add smoothing, it should be small + self._v_projector.smoothing = 1.0e-6 ## This function is the one we will typically over-ride to build specific solvers. ## This example is a poisson-like problem with isotropic coefficients @@ -357,11 +347,11 @@ class SNES_Stokes(SNES_Stokes_SaddlePt): \nabla \cdot \color{Blue}{\underbrace{\Bigl[ \boldsymbol{\tau} - p \mathbf{I} \Bigr]}_{\mathbf{F}}} = - \color{Maroon}{\underbrace{\Bigl[ \mathbf{f} \Bigl] }_{\mathbf{f}}} + \color{Maroon}{\underbrace{\Bigl[ \mathbf{f} \Bigl] }_{\mathbf{h}}} $$ $$ - \underbrace{\Bigl[ \nabla \cdot \mathbf{u} \Bigr]}_{\mathbf{f}_p} = 0 + \underbrace{\Bigl[ \nabla \cdot \mathbf{u} \Bigr]}_{\mathbf{h}_p} = 0 $$ The flux term is a deviatoric stress ( $\boldsymbol{\tau}$ ) related to velocity gradients @@ -371,7 +361,7 @@ class SNES_Stokes(SNES_Stokes_SaddlePt): \mathbf{F}: \quad \boldsymbol{\tau} = \frac{\eta}{2}\left( \nabla \mathbf{u} + \nabla \mathbf{u}^T \right) $$ - The constraint equation, $\mathbf{f}_p = 0$ is incompressible flow by default but can be set + The constraint equation, $\mathbf{h}_p = 0$ gives incompressible flow by default but can be set to any function of the unknown $\mathbf{u}$ and $\nabla\cdot\mathbf{u}$ ## Properties @@ -417,7 +407,6 @@ def __init__( pressureField: Optional[uw.discretisation.MeshVariable] = None, degree: Optional[int] = 2, p_continuous: Optional[bool] = True, - solver_name: Optional[str] = "", verbose: Optional[bool] = False, # Not used in Stokes, but may be used in NS, VE etc DuDt: Union[SemiLagrangian_DDt, Lagrangian_DDt] = None, @@ -429,16 +418,11 @@ def __init__( pressureField, degree, p_continuous, - solver_name, verbose, DuDt=DuDt, DFDt=DFDt, ) - # change this to be less generic - if solver_name == "": - self.name = "Stokes_{}_".format(self.instance_number) - self._degree = degree # User-facing operations are matrices / vectors by preference @@ -497,9 +481,9 @@ def F1(self): def PF0(self): f0 = uw.function.expression( - r"\mathbf{f}_0\left( \mathbf{p} \right)", + r"\mathbf{h}_0\left( \mathbf{p} \right)", sympy.simplify(sympy.Matrix((self.constraints))), - "Pointwise flux term: f_0(p)", + "Pointwise force term: h_0(p)", ) # backward compatibility @@ -675,7 +659,6 @@ def __init__( degree: Optional[int] = 2, order: Optional[int] = 2, p_continuous: Optional[bool] = True, - solver_name: Optional[str] = "", verbose: Optional[bool] = False, # DuDt Not used in VE, but may be in child classes DuDt: Union[SemiLagrangian_DDt, Lagrangian_DDt] = None, @@ -689,16 +672,11 @@ def __init__( pressureField, degree, p_continuous, - solver_name, verbose, DuDt=DuDt, DFDt=DFDt, ) - # change this to be less generic - if solver_name == "": - self.name = "VE_Stokes_{}_".format(self.instance_number) - self._order = order # VE time-order if self.Unknowns.DFDt is None: @@ -819,20 +797,16 @@ def __init__( mesh: uw.discretisation.Mesh, u_Field: uw.discretisation.MeshVariable = None, degree=2, - solver_name: str = "", verbose=False, ): + super().__init__( mesh, u_Field, degree, - solver_name, verbose, ) - if solver_name == "": - self.name = "SProj_{}_".format(self.instance_number) - self.is_setup = False self._smoothing = sympy.sympify(0) self._uw_weighting_function = sympy.sympify(1) @@ -870,21 +844,6 @@ def F1(self): return F1_val - # @timing.routine_timer_decorator - # def projection_problem_description(self): - # # residual terms - defines the problem: - # # solve for a best fit to the continuous mesh - # # variable given the values in self.function - # # F0 is left in place for the user to inject - # # non-linear constraints if required - - # self._f0 = self.F0.sym - - # # F1 is left in the users control ... e.g to add other gradient constraints to the stiffness matrix - - # self._f1 = self.F1.sym - - # return @property def uw_function(self): @@ -950,20 +909,15 @@ def __init__( mesh: uw.discretisation.Mesh, u_Field: uw.discretisation.MeshVariable = None, degree=2, - solver_name: str = "", verbose=False, ): super().__init__( mesh, u_Field, degree, - solver_name, verbose, ) - if solver_name == "": - solver_name = "VProj{}_".format(self.instance_number) - self.is_setup = False self._smoothing = 0.0 self._penalty = 0.0 @@ -1095,7 +1049,6 @@ def __init__( tensor_Field: uw.discretisation.MeshVariable = None, scalar_Field: uw.discretisation.MeshVariable = None, degree=2, - solver_name: str = "", verbose=False, ): self.t_field = tensor_Field @@ -1104,13 +1057,9 @@ def __init__( mesh=mesh, u_Field=scalar_Field, degree=degree, - solver_name=solver_name, verbose=verbose, ) - if solver_name == "": - solver_name = "TProj{}_".format(self.instance_number) - return ## Need to over-ride solve method to run over all components @@ -1250,7 +1199,6 @@ def __init__( uw.discretisation.MeshVariable, sympy.Basic ], # Should be a sympy function order: int = 1, - solver_name: str = "", restore_points_func: Callable = None, verbose=False, DuDt: Union[SemiLagrangian_DDt, Lagrangian_DDt] = None, @@ -1261,15 +1209,11 @@ def __init__( mesh, u_Field, u_Field.degree, - solver_name, verbose, DuDt=DuDt, DFDt=DFDt, ) - if solver_name == "": - solver_name = "AdvDiff_slcn_{}_".format(self.instance_number) - if isinstance(V_fn, uw.discretisation._MeshVariable): self._V_fn = V_fn.sym else: @@ -1595,7 +1539,6 @@ def __init__( restore_points_func: Callable = None, order: Optional[int] = 2, p_continuous: Optional[bool] = False, - solver_name: Optional[str] = "", verbose: Optional[bool] = False, DuDt: Union[SemiLagrangian_DDt, Lagrangian_DDt] = None, DFDt: Union[SemiLagrangian_DDt, Lagrangian_DDt] = None, @@ -1607,7 +1550,6 @@ def __init__( pressureField, order, p_continuous, - solver_name, verbose, DuDt=DuDt, DFDt=DFDt, diff --git a/tests/test_1100_AdvDiffCartesian.py b/tests/test_1100_AdvDiffCartesian.py index 4344b0eca..b74f623b1 100644 --- a/tests/test_1100_AdvDiffCartesian.py +++ b/tests/test_1100_AdvDiffCartesian.py @@ -70,7 +70,6 @@ def test_advDiff_boxmesh(mesh): mesh, u_Field=T, V_fn=v, - solver_name="adv_diff", ) # ### Set up properties of the adv_diff solver diff --git a/tests/test_1110_advDiffAnnulus.py b/tests/test_1110_advDiffAnnulus.py index 2eafa0514..74891d405 100644 --- a/tests/test_1110_advDiffAnnulus.py +++ b/tests/test_1110_advDiffAnnulus.py @@ -28,7 +28,6 @@ def test_adv_diff_annulus(): mesh, u_Field=t_soln, V_fn=v_soln, - solver_name="adv_diff_annulus", ) adv_diff.constitutive_model = uw.constitutive_models.DiffusionModel