checkpoint

block_matrix.py

@@ -262,20 +262,20 @@ def inner(input_dofs_idx, take_groups, all_groups):
         I, J = np.meshgrid(row_idx, col_idx, sparse=True, indexing="ij", copy=False)
         self.mat[I, J] = value
 
-    def group_shape(self, i, j=None) -> tuple[int, int]:
-        groups_i, groups_j = self._correct_getitem_key((i, j) if j is not None else i)
+    # def group_shape(self, i, j=None) -> tuple[int, int]:
+    #     groups_i, groups_j = self._correct_getitem_key((i, j) if j is not None else i)
 
-        def inner(input_dofs_idx, take_groups, all_groups):
-            num_dofs = 0
-            for group in take_groups:
-                for dof_idx in all_groups[group]:
-                    if (dofs := input_dofs_idx[dof_idx]) is not None:
-                        num_dofs += dofs.size
-            return num_dofs
+    #     def inner(input_dofs_idx, take_groups, all_groups):
+    #         num_dofs = 0
+    #         for group in take_groups:
+    #             for dof_idx in all_groups[group]:
+    #                 if (dofs := input_dofs_idx[dof_idx]) is not None:
+    #                     num_dofs += dofs.size
+    #         return num_dofs
 
-        row_shape = inner(self.local_row_idx, groups_i, self.groups_row)
-        col_shape = inner(self.local_col_idx, groups_j, self.groups_col)
-        return row_shape, col_shape
+    #     row_shape = inner(self.local_row_idx, groups_i, self.groups_row)
+    #     col_shape = inner(self.local_col_idx, groups_j, self.groups_col)
+    #     return row_shape, col_shape
 
     def slice_domain(self, i, j) -> spmatrix:
         active_subgroups_i, active_subgroups_j = self.active_subgroups
@@ -284,23 +284,9 @@ def slice_domain(self, i, j) -> spmatrix:
         I, J = np.meshgrid(row_idx, col_idx, sparse=True, indexing="ij", copy=False)
         return self.mat[I, J]
 
-    def build_schur_complement(
-        self, keep: Sequence[int], elim: Sequence[int], only_stab: bool = False
-    ) -> "BlockMatrixStorage":
-        m_kk = self[keep].copy()
-        m_ee = self[elim].mat
-        m_ke = self[keep, elim].mat
-        m_ek = self[elim, keep].mat
-        stab = m_ke @ inv(m_ee) @ m_ek
-        if only_stab:
-            m_kk.mat = stab
-        else:
-            m_kk.mat -= stab
-        return m_kk
-
     def block_diag(self) -> "BlockMatrixStorage":
         """Returns a copy. Keeps only diagonal blocks.
-        E.g. how removes how one fracture affects the other."""
+        E.g. removes how one fracture affects the other."""
         active_idx = [x for x in self.local_row_idx if x is not None]
         bmats = []
         for active in active_idx:
@@ -340,6 +326,20 @@ def empty_container(self) -> "BlockMatrixStorage":
             group_row_names=self.groups_row_names,
         )
 
+    def make_restriction_matrix(self, to: list[int]) -> "BlockMatrixStorage":
+        self = self
+        active_cols = np.concatenate(
+            [self.local_col_idx[j] for i in to for j in self.groups_col[i]]
+        )
+        active_rows = np.arange(active_cols.size)
+        active_data = np.ones(active_rows.size)
+        result = self[to, self.active_groups[1]].empty_container()
+        proj = scipy.sparse.coo_matrix(
+            (active_data, (active_rows, active_cols)), result.shape
+        )
+        result.mat = proj.tocsr()
+        return result
+
     def local_rhs(self, global_rhs: np.ndarray) -> np.ndarray:
         row_idx = [
             self.global_row_idx[j]
@@ -583,14 +583,14 @@ def color_local_rhs(
 
 
 @dataclass
-class SolveSchema:
+class FieldSplitScheme:
     groups: list[int]
     solve: callable | Literal["direct", "use_invertor"] = "direct"
     invertor: callable | Literal["use_solve", "direct"] = "use_solve"
     invertor_type: Literal["physical", "algebraic", "operator", "test_vector"] = (
         "algebraic"
     )
-    complement: Optional["SolveSchema"] = None
+    complement: Optional["FieldSplitScheme"] = None
     factorization_type: Literal["full", "upper"] = "upper"
 
     compute_cond: bool = False
@@ -610,15 +610,20 @@ def __str__(self):
         return res
 
 
-def get_complement_groups(schema: SolveSchema):
+@dataclass
+class MultiStageScheme:
+    stages: list[FieldSplitScheme]
+
+
+def get_complement_groups(schema: FieldSplitScheme):
     res = []
     if schema.complement is not None:
         res.extend(schema.complement.groups)
         res.extend(get_complement_groups(schema=schema.complement))
     return res
 
 
-def make_solver(schema: SolveSchema, mat_orig: BlockMatrixStorage):
+def make_solver(schema: FieldSplitScheme, mat_orig: BlockMatrixStorage):
     groups_0 = schema.groups
     groups_1 = get_complement_groups(schema)
 

experiments/runscript2.py

@@ -4,7 +4,16 @@
 
 def experiment_2():
     setups = []
-    for grid_refinement in [1, 2, 3, 4, 5, 6, 10, 33]:
+    for grid_refinement in [
+        1,
+        # 2,
+        # 3,
+        # 4,
+        # 5,
+        # 6,
+        # 10,
+        # 33,
+    ]:
         setups.append(
             {
                 "physics": 1,

experiments/thermal/runscript1.py

@@ -0,0 +1,89 @@
+from thermal_models import run_model, make_model
+from plot_utils import write_dofs_info
+
+
+setup_reference = {
+    "geometry": 1,  # 1 - 2D 1 fracture, 2 - 2D 7 fractures;
+    "barton_bandis_stiffness_type": 1,  # 0 - off; 1 - small; 2 - medium, 3 - large
+    "friction_type": 1,  # 0 - small, 1 - medium, 2 - large
+    "grid_refinement": 1,  # 1 - coarsest level
+    "solver": 2,  # 0 - Direct solver, 1 - Richardson + fixed stress splitting scheme, 2 - GMRES + scalable prec.
+    "save_matrix": False,  # Save each matrix and state array. Might take space.
+}
+
+
+def experiment_1_barton_bandis_friction():
+    setups = []
+    for barton_bandis in [
+        0,
+        1,
+        2,
+        3,
+        4,
+        5,
+    ]:
+        for friction in [0, 1, 2]:
+            for solver in [1, 2]:
+                setups.append(
+                    {
+                        "physics": 0,
+                        "geometry": 1,
+                        "barton_bandis_stiffness_type": barton_bandis,
+                        "friction_type": friction,
+                        "grid_refinement": 1,
+                        "solver": solver,
+                        "save_matrix": True,
+                    }
+                )
+    for setup in setups:
+        model = make_model(setup)
+        run_model(setup)
+        write_dofs_info(model)
+        print(model.simulation_name())
+
+
+def experiment_1_grid_refinement():
+    setups = []
+    # for grid_refinement in [1, 2, 3, 4]:
+    #     for solver in [1, 11, 12]:
+    #         setups.append(
+    #             {
+    #                 "physics": 1,
+    #                 "geometry": 1,
+    #                 "barton_bandis_stiffness_type": 2,
+    #                 "friction_type": 1,
+    #                 "grid_refinement": grid_refinement,
+    #                 "solver": solver,
+    #                 "save_matrix": True,
+    #             }
+    #         )
+    for grid_refinement in [
+        1,
+        # 2,
+        # 3,
+        # 4,
+        # 5,
+        # 6,
+        # 10,
+        # 33,
+    ]:
+        setups.append(
+            {
+                "physics": 1,
+                "geometry": 1,
+                "barton_bandis_stiffness_type": 2,
+                "friction_type": 1,
+                "grid_refinement": grid_refinement,
+                "solver": 0,
+                "save_matrix": True,
+            }
+        )
+    for setup in setups:
+        model = make_model(setup)
+        run_model(setup)
+        write_dofs_info(model)
+        print(model.simulation_name())
+
+
+if __name__ == "__main__":
+    experiment_1_grid_refinement()