Minor updates in DOLFINx fracture code

examples/03_fracture/fracture_dolfinx.py

@@ -1,7 +1,8 @@
 from argparse import ArgumentDefaultsHelpFormatter, ArgumentParser
 
 from mpi4py import MPI
-
+from typing import Literal
+from packaging.version import Version
 import basix.ufl
 import dolfinx.fem.petsc
 import numpy as np
@@ -16,6 +17,7 @@
     grad,
     inner,
     split,
+    core,
 )
 
 from lvpp import SNESProblem, SNESSolver
@@ -29,6 +31,20 @@
 class NotConvergedError(Exception):
     pass
 
+def norm(expr: core.expr.Expr, norm_type: Literal["L2", "H1", "H10"]) -> dolfinx.fem.Form:
+    """
+    Compile the norm of an UFL expression into a DOLFINx form.
+
+    Args:
+        expr: UFL expression
+        norm_type: Type of norm
+    """
+    if norm_type == "L2":
+        return dolfinx.fem.form(inner(expr, expr) * dx)
+    elif norm_type == "H1":
+        return norm(expr, "L2") + norm(grad(expr), "L2")
+    elif norm_type == "H10":
+        return norm(grad(expr), "L2")
 
 parser = ArgumentParser(
     description="Solve the obstacle problem on a unit square using Galahad.",
@@ -63,7 +79,8 @@ class NotConvergedError(Exception):
     dest="num_load_steps",
     help="Number of load steps",
 )
-parser.add_argument("--verbose", "-v", action="store_true", help="Verbose output from PETSc SNES")
+parser.add_argument("--verbose", "-v", action="store_true",
+                    help="Verbose output from PETSc SNES")
 parser.add_argument("--Tmin", type=float, default=0.0, help="Minimum load")
 parser.add_argument("--Tmax", type=float, default=5.0, help="Maximum load")
 
@@ -88,13 +105,17 @@ class NotConvergedError(Exception):
 
 # Compute maximum cell size
 W = dolfinx.fem.functionspace(mesh, ("DG", 0))
-h = dolfinx.fem.Expression(4 * Circumradius(mesh), W.element.interpolation_points())
+if Version(dolfinx.__version__) > Version("0.9.0"):
+    ip = W.element.interpolation_points
+else:
+    ip = W.element.interpolation_points()
+h = dolfinx.fem.Expression(4 * Circumradius(mesh), ip)
 diam = dolfinx.fem.Function(W)
 diam.interpolate(h)
 max_diam = mesh.comm.allreduce(np.max(diam.x.array), op=MPI.MAX)
 l = dolfinx.fem.Constant(mesh, st(max_diam))
 print(f"Using l = {float(l)}")
-
+print("System size: ", Z.dofmap.index_map.size_global*Z.dofmap.index_map_bs)
 
 z = dolfinx.fem.Function(Z)
 (u, c, psi) = split(z)
@@ -107,14 +128,19 @@ class NotConvergedError(Exception):
 z_prev = dolfinx.fem.Function(Z)
 _, c_prev, _ = split(z_prev)
 z_iter = dolfinx.fem.Function(Z)
-(_, c_iter, psi_iter) = split(z_iter)
+(u_iter, c_iter, psi_iter) = split(z_iter)
 
 
 output_space = dolfinx.fem.functionspace(mesh, ("Lagrange", 3))
 c_conform_out = dolfinx.fem.Function(output_space, name="ConformingDamage")
 alpha = dolfinx.fem.Constant(mesh, st(1.0))
 c_conform = (c_prev + exp(psi)) / (exp(psi) + 1)
-c_conform_expr = dolfinx.fem.Expression(c_conform, output_space.element.interpolation_points())
+if Version(dolfinx.__version__) > Version("0.9.0"):
+    out_ip = output_space.element.interpolation_points
+else:
+    out_ip = output_space.element.interpolation_points()
+c_conform_expr = dolfinx.fem.Expression(
+    c_conform, out_ip)
 
 eps = dolfinx.fem.Constant(mesh, 1.0e-5)
 E = (
@@ -193,8 +219,12 @@ class NotConvergedError(Exception):
 xdmf_file.write_mesh(mesh)
 xdmf_file.close()
 vtx_damage = dolfinx.io.VTXWriter(mesh.comm, "damage.bp", [c_conform_out])
-L2_c = dolfinx.fem.form(inner(c - c_iter, c - c_iter) * dx)
-L2_z = dolfinx.fem.form(inner(z - z_prev, z - z_prev) * dx)
+diff_c = c - c_iter
+diff_u = u - u_iter
+diff_z = z - z_prev
+norm_c = norm(diff_c, "L2")
+norm_u = norm(diff_u, "L2")
+norm_z = norm(diff_z, "L2")
 
 U, U_to_Z = Z.sub(0).collapse()
 u_out = dolfinx.fem.Function(U, name="u")
@@ -207,6 +237,8 @@ class NotConvergedError(Exception):
 
 
 problem = SNESProblem(F_compiled, z, bcs=bcs, J=J_reg)
+solver = SNESSolver(problem, sp)
+
 for step, T in enumerate(np.linspace(Tmin, Tmax, num_load_steps)[1:]):
     if mesh.comm.rank == 0:
         print(
@@ -224,9 +256,16 @@ class NotConvergedError(Exception):
         try:
             if mesh.comm.rank == 0:
                 print(f"Attempting {k=} alpha={float(alpha)}", flush=True)
+            del solver
             solver = SNESSolver(problem, sp)
+            import time
+            start = time.perf_counter()
             converged_reason, num_iterations = solver.solve()
-
+            end = time.perf_counter()
+            print(f"Solve time: {end-start:.2e}")
+            if converged_reason == -3:
+                # Linear solver did not converge
+                raise NotConvergedError(f"Not converged, linear solver failed with {solver._snes.getKSP().getConvergedReason()}")
             if num_iterations == 0 and converged_reason > 0:
                 # solver didn't actually get to do any work,
                 # we've just reduced alpha so much that the initial guess
@@ -258,14 +297,17 @@ class NotConvergedError(Exception):
                 continue
 
         # Termination
-        nrm = np.sqrt(mesh.comm.allreduce(dolfinx.fem.assemble_scalar(L2_c), op=MPI.SUM))
+        nrm_c = np.sqrt(mesh.comm.allreduce(dolfinx.fem.assemble_scalar(
+            norm_c), op=MPI.SUM))
+        nrm_u = np.sqrt(mesh.comm.allreduce(dolfinx.fem.assemble_scalar(
+            norm_u), op=MPI.SUM))
         if mesh.comm.rank == 0:
             print(
                 f"{_GREEN}Solved {k=} {num_iterations=} alpha={float(alpha)},"
-                f"||c_{k} - c_{k - 1}|| = {nrm}{_color_reset}",
+                f"||c_{k} - c_{k - 1}|| = {nrm_c:.3e} ||u_{k} - u_{k - 1}|| = {nrm_u:.3e} {_color_reset}",
                 flush=True,
             )
-        if nrm < 1.0e-4:
+        if nrm_c < 1.0e-4 and nrm_u < 1.0e-4:
             break
 
         # Update alpha
@@ -283,15 +325,16 @@ class NotConvergedError(Exception):
     # the failure mode of the algorithm above is that it terminates in one
     # PG iteration that does no Newton iterations, so the solution doesn't
     # change
-    norm_Z = np.sqrt(mesh.comm.allreduce(dolfinx.fem.assemble_scalar(L2_z), op=MPI.SUM))
-    if k == 1 and np.isclose(norm_Z, 0.0):
+    norm_LVPP = np.sqrt(mesh.comm.allreduce(
+        dolfinx.fem.assemble_scalar(norm_z), op=MPI.SUM))
+    if k == 1 and np.isclose(norm_LVPP, 0.0):
         break
 
     if nfail == NFAIL_MAX:
         break
-    c_conform_out.interpolate(c_conform_expr)
-
+    z_prev.interpolate(z)
     if step % write_frequency == 0:
+        c_conform_out.interpolate(c_conform_expr)
         with dolfinx.io.XDMFFile(mesh.comm, "solution.xdmf", "a") as xdmf_file:
             u_out.x.array[:] = z.x.array[U_to_Z]
             c_out.x.array[:] = z.x.array[C_to_Z]
@@ -300,6 +343,5 @@ class NotConvergedError(Exception):
             xdmf_file.write_function(c_out, T)
             xdmf_file.write_function(psi_out, T)
 
-        z_prev.interpolate(z)
         vtx_damage.write(T)
 vtx_damage.close()

examples/03_fracture/generate_mesh.py

@@ -58,7 +58,6 @@ def create_crack_mesh(comm, max_res: float = 0.05):
 
     ud = ufl.Mesh(basix.ufl.element("Lagrange", "triangle", 1, shape=(2,)))
     linear_mesh = dolfinx.mesh.create_mesh(MPI.COMM_WORLD, cells, x, ud)
-
     if comm.rank == 0:
         regions: dict[str, list[int]] = {name: [] for name in ngmesh.GetRegionNames(codim=1)}
         for i, name in enumerate(ngmesh.GetRegionNames(codim=1), 1):
@@ -80,7 +79,7 @@ def create_crack_mesh(comm, max_res: float = 0.05):
         facet_values = np.array([facet.index for facet in ng_facets], dtype=np.int32)
         regions = comm.bcast(regions, root=0)
     else:
-        facets = np.zeros((0, 3), dtype=np.int64)
+        facets = np.zeros((0, 2), dtype=np.int64)
         facet_values = np.zeros((0,), dtype=np.int32)
         regions = comm.bcast(None, root=0)
     local_entities, local_values = dolfinx.io.gmshio.distribute_entity_data(