Clean solvers

Author: a-latyshev

.gitignore

@@ -19,4 +19,6 @@ build/
 
 *.out
 *.pkl
-*.txt
+*.txt
+
+*.core

doc/demo/demo_plasticity_mohr_coulomb.py

@@ -9,6 +9,10 @@
 #       format_name: percent
 #       format_version: '1.3'
 #       jupytext_version: 1.11.2
+#   kernelspec:
+#     display_name: dolfinx-env
+#     language: python
+#     name: python3
 # ---
 
 # %% [markdown]
@@ -113,7 +117,7 @@
 
 # %%
 L, H = (1.2, 1.0)
-Nx, Ny = (200, 200)
+Nx, Ny = (25, 25)
 gamma = 1.0 / 6778
 domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([L, H])], [Nx, Ny])
 
@@ -669,10 +673,11 @@ def F_ext(v):
 # parameters of the manual Newton method
 max_iterations, relative_tolerance = 200, 1e-8
 
-load_steps_1 = np.linspace(1.5, 21, 40)
-load_steps_2 = np.linspace(21, 22.75, 20)[1:]
+load_steps_1 = np.linspace(2, 22.9, 50)
+load_steps_2 = np.array([22.96, 22.99])
 load_steps = np.concatenate([load_steps_1, load_steps_2])
-load_steps = np.concatenate([np.linspace(1.1, 22.3, 100)[:-1]])
+
+# load_steps = np.concatenate([np.linspace(1.1, 22.3, 100)[:-1]])
 num_increments = len(load_steps)
 results = np.zeros((num_increments + 1, 2))
 
@@ -702,8 +707,6 @@ def constitutive_update():
 external_operator_problem = SNESProblem(Du, F_replaced, J_replaced, bcs=bcs, petsc_options=petsc_options, system_update=constitutive_update)
 
 # %%
-timer_total = common.Timer("Total_timer")
-timer_total.start()
 for i, load in enumerate(load_steps):
     q.value = load * np.array([0, -gamma])
 
@@ -721,14 +724,8 @@ def constitutive_update():
 
     if len(points_on_process) > 0:
         results[i + 1, :] = (-u.eval(points_on_process, cells)[0], load)
-timer_total.stop()
-total_time = timer_total.elapsed()[0]
 
 print(f"Slope stability factor: {-q.value[-1]*H/c}")
-print(f"Total time: {total_time}")
-
-# %%
-external_operator_problem.performance_monitor
 
 # %% [markdown]
 # ## Verification
@@ -768,34 +765,6 @@ def constitutive_update():
     plt.grid()
     plt.legend()
 
-# %%
-external_operator_problem.performance_monitor
-
-# %%
-import pickle
-with open("performance_monitor_100_100", "wb") as f:
-        pickle.dump(external_operator_problem.performance_monitor, f)
-
-# %%
-# pickle_file = Path("performance_monitor_200_200") 
-with open("performance_monitor_200_200", 'rb') as f:
-    performance_monitor_200_200 = pickle.load(f)
-
-# %%
-summary_monitor = performance_monitor_200_200.copy()
-
-cols = ["matrix_assembling", "vector_assembling", "linear_solver", "constitutive_model_update"]
-summary_monitor["linear_solver"] = summary_monitor["nonlinear_solver"] - summary_monitor["matrix_assembling"] - summary_monitor["vector_assembling"] - summary_monitor["constitutive_model_update"]
-
-fig, ax = plt.subplots(figsize=(10, 5))
-summary_monitor.plot(use_index=True, y=cols, kind="bar", stacked=True, ax=ax)
-
-# %%
-fig, ax = plt.subplots(figsize=(10, 5))
-for col in cols:
-    summary_monitor[col] = summary_monitor[col] / (summary_monitor["Newton_iterations"]+1)
-summary_monitor.plot(use_index=True, y=cols, kind="bar", stacked=True, ax=ax)
-
 # %% [markdown]
 # The slope profile reaching its stability limit:
 

doc/demo/solvers.py

@@ -110,17 +110,6 @@ def __init__(
 
         self.solver = self.solver_setup()
 
-        self.local_monitor = {"matrix_assembling": 0.0, "vector_assembling": 0.0, "constitutive_model_update": 0.0}
-        self.performance_monitor = pd.DataFrame({
-            # "loading_step": np.array([], dtype=np.int64),
-            "Newton_iterations": np.array([], dtype=np.int64),
-            "matrix_assembling": np.array([], dtype=np.float64),
-            "vector_assembling": np.array([], dtype=np.float64),
-            "nonlinear_solver": np.array([], dtype=np.float64),
-            "constitutive_model_update": np.array([], dtype=np.float64),
-        })
-        self.timer = common.Timer("SNES")
-
     def set_petsc_options(self):
         # Set PETSc options
         opts = PETSc.Options()
@@ -159,22 +148,16 @@ def F(self, snes: PETSc.SNES, x: PETSc.Vec, b: PETSc.Vec) -> None:
 
         #TODO: SNES makes the iteration #0, where it calculates the b norm.
         #`system_update()` can be omitted in that case
-        self.timer.start()
         self.system_update()
-        self.timer.stop()
-        self.local_monitor["constitutive_model_update"] += self.comm.allreduce(self.timer.elapsed()[0], op=MPI.SUM)
-
-        self.timer.start()
+       
         with b.localForm() as b_local:
             b_local.set(0.0)
         fem.petsc.assemble_vector(b, self.F_form)
 
         fem.petsc.apply_lifting(b, [self.J_form], [self.bcs], [x], -1.0)
         b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
         fem.petsc.set_bc(b, self.bcs, x, -1.0)
-        self.timer.stop()
-        self.local_monitor["vector_assembling"] += self.comm.allreduce(self.timer.elapsed()[0], op=MPI.SUM)
-
+       
     def J(self, snes, x: PETSc.Vec, A: PETSc.Mat, P: PETSc.Mat) -> None:
         """Assemble the Jacobian matrix.
 
@@ -183,26 +166,13 @@ def J(self, snes, x: PETSc.Vec, A: PETSc.Mat, P: PETSc.Mat) -> None:
         x: Vector containing the latest solution.
         A: Matrix to assemble the Jacobian into.
         """
-        self.timer.start()
         A.zeroEntries()
         fem.petsc.assemble_matrix(A, self.J_form, self.bcs)
         A.assemble()
-        self.timer.stop()
-        self.local_monitor["matrix_assembling"] += self.comm.allreduce(self.timer.elapsed()[0], op=MPI.SUM)
-
+        
     def solve(self,) -> Tuple[int, int]:
-        # self.local_monitor["loading_step"] = loading_step
-        self.local_monitor["vector_assembling"] = 0.0
-        self.local_monitor["matrix_assembling"] = 0.0
-        self.local_monitor["constitutive_model_update"] = 0.0
-        timer = common.Timer("nonlinear_solver")
-        self.timer.start()
         self.solver.solve(None, self.u.x.petsc_vec)
-        timer.stop()
-        self.local_monitor["nonlinear_solver"] = self.comm.allreduce(self.timer.elapsed()[0], op=MPI.SUM)
-        self.local_monitor["Newton_iterations"] = self.solver.getIterationNumber()
         self.u.x.scatter_forward()
-        self.performance_monitor.loc[len(self.performance_monitor.index)] = self.local_monitor
         return (self.solver.getIterationNumber(), self.solver.getConvergedReason())
 
     def __del__(self):

doc/demo/utilities.py

@@ -60,7 +60,7 @@ def build_cylinder_quarter(lc=0.3, R_e=1.3, R_i=1.0):
 
     # NOTE: Do not forget to check the leaks produced by the following line of code
     # [WARNING] yaksa: 2 leaked handle pool objects
-    mesh, cell_tags, facet_tags = gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, rank=model_rank, gdim=2)
+    mesh, cell_tags, facet_tags, _, _, _ = gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, rank=model_rank, gdim=2)
 
     mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
     mesh.name = "quarter_cylinder"