support reusing the solver with different guesses

Author: pbrubeck

firedrake/adaptive_variational_solver.py

@@ -1,3 +1,4 @@
+from firedrake import dmhooks
 from firedrake.petsc import PETSc
 from firedrake.function import Function
 from firedrake.cofunction import Cofunction
@@ -91,7 +92,19 @@ def __init__(self, problem: NonlinearVariationalProblem, goal_functional, tolera
         self.degree = self.element.degree()
         self.nullspace = nullspace
         self.atm = AdaptiveTransferManager()
-        self.amh = AdaptiveMeshHierarchy(V.mesh())
+
+        mesh = V.mesh().unique()
+        amh, level = utils.get_level(mesh)
+        if amh is None:
+            amh = AdaptiveMeshHierarchy(mesh)
+        if not isinstance(amh, AdaptiveMeshHierarchy):
+            raise ValueError("Problem needs to be defined on an AdaptiveMeshHierarchy")
+
+        # FIXME we should construct a subset of the original hierarchy instead
+        amh, level = utils.get_level(mesh)
+        while len(amh) > level+1:
+            amh.pop_mesh()
+        self.amh = amh
 
         # Data storage and writing
         self.output_dir = Path(self.options.output_dir)

firedrake/mg/adaptive_hierarchy.py

@@ -21,7 +21,7 @@ class AdaptiveMeshHierarchy(HierarchyBase):
     """
     def __init__(self, base_mesh: MeshGeometry, nested: bool = True):
         self.meshes = []
-        self._meshes = []
+        self._meshes = self.meshes
         self.nested = nested
         self.add_mesh(base_mesh)
 
@@ -35,10 +35,20 @@ def add_mesh(self, mesh: MeshGeometry):
             The mesh to be added to the finest level.
         """
         level = len(self.meshes)
-        self._meshes.append(mesh)
         self.meshes.append(mesh)
         set_level(mesh, self, level)
 
+    def pop_mesh(self):
+        """
+        Adds a mesh into the hierarchy.
+
+        Parameters
+        ----------
+        mesh
+            The mesh to be added to the finest level.
+        """
+        self.meshes.pop()
+
     def adapt(self, eta: Function | Cofunction, theta: float):
         """
         Adds a new mesh to the hierarchy by locally refining the finest mesh

firedrake/mg/netgen.py

@@ -246,8 +246,6 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None):
     logger.info(f"\tOrder of the hierarchy: {order}")
     if isinstance(order, int):
         order = [order]*(levs+1)
-    permutation_tol = flags.get("permutation_tol", 1e-8)
-    location_tol = flags.get("location_tol", 1e-8)
     refType = flags.get("refinement_type", "uniform")
     logger.info(f"\tRefinement type: {refType}")
     optMoves = flags.get("optimisation_moves", False)
@@ -270,7 +268,7 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None):
         temp_flags = dict(flags)
         temp_flags['degree'] = order[0]
         temp = fd.Mesh(mesh.netgen_mesh, distribution_parameters=parameters,
-                      netgen_flags=temp_flags, comm=comm)
+                       netgen_flags=temp_flags, comm=comm)
         mesh = temp
     # Make a plex (cdm) without overlap.
     dm_cell_type, = mesh.dm_cell_types
@@ -321,6 +319,7 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None):
                     netgen_flags=temp_flags,
                     comm=comm)
             elif snap == "coarse":
+                ho_field = meshes[0].coordinates
                 mesh = snapToCoarse(ho_field, mesh, order[l+1], snap_smoothing, cg)
             toc = time.time()
             logger.info(f"\t\t\tMeshed curved. Time taken: {toc-tic}")

firedrake/mg/ufl_utils.py

@@ -145,8 +145,9 @@ def coarsen_equation_bc(ebc, self, coefficient_mapping=None):
 
 @coarsen.register(firedrake.functionspaceimpl.WithGeometryBase)
 def coarsen_function_space(V, self, coefficient_mapping=None):
-    if hasattr(V, "_coarse") and self == coarsen:
-        return V._coarse
+    # FIXME support multiple hierarchies
+    # if hasattr(V, "_coarse") and self == coarsen:
+    #     return V._coarse
     V_fine = V
     # Handle MixedFunctionSpace : V_fine.reconstruct requires MeshSequence.
     fine_mesh = V_fine.mesh() if V_fine.index is None else V_fine.parent.mesh()