eig shift in pararmeters

pyproject.toml

@@ -38,4 +38,4 @@ irrevolutions = "irrevolutions:main"
 
 # Optional include files like configuration, license, etc.
 [tool.setuptools.package-data]
-irrevolutions = ["models/default_parameters.yml"]
+irrevolutions = ["models/default_parameters.yml", "test/parameters.yml"]

src/irrevolutions/algorithms/so.py

@@ -209,7 +209,7 @@ def setup_eigensolver(self, eigen):
 
         st = eigen.eps.getST()
         st.setType("sinvert")
-        st.setShift(-1.0e-3)
+        st.setShift(self.parameters["eigen"]["shift"])
 
         eigen.eps.setTolerances(
             self.parameters["eigen"]["eig_rtol"], self.parameters["eigen"]["eps_max_it"]
@@ -333,12 +333,16 @@ def normalise_eigen(self, u, mode="max-beta"):
             logging.debug(f"{rank}, |β|_infty {beta.x.petsc_vec.norm(3):.3f}")
 
         elif mode == "unit":
-            coeff_glob = np.sqrt(sum(n**2 for n in [v_i.x.petsc_vec.norm() for v_i in u]))
+            coeff_glob = np.sqrt(
+                sum(n**2 for n in [v_i.x.petsc_vec.norm() for v_i in u])
+            )
             logging.debug(f"rank {rank}, coeff_glob {coeff_glob:.3f}")
             logging.debug(f"{rank}, |(v, β)^*|_2 {coeff_glob:.3f}")
 
         if coeff_glob == 0.0:
-            logging.error(f"Damage eigenvector is null i.e. |β|={beta.x.petsc_vec.norm()}")
+            logging.error(
+                f"Damage eigenvector is null i.e. |β|={beta.x.petsc_vec.norm()}"
+            )
             return 0.0
 
         for v_i in u:
@@ -400,7 +404,7 @@ def solve(self, alpha_old: dolfinx.fem.function.Function):
             # Sort eigenmodes by eigenvalues
             spectrum.sort(key=lambda item: item.get("lambda"))
             # unstable_spectrum = list(filter(lambda item: item.get("lambda") <= 0, spectrum))
-            
+
             # Store the results
             stable = self.store_results(eigen, spectrum)
 
@@ -635,28 +639,35 @@ def __init__(
             stability_parameters=bifurcation_parameters,
         )
 
-    def log(self, logger = logger):
+    def log(self, logger=logger):
         # Check if spectrum is available
         if not self._spectrum:
             logger.info("No negative spectrum.")
 
         # Find the minimum eigenvalue
-        min_eigenvalue = min(entry["lambda"] for entry in self.spectrum if "lambda" in entry)
+        min_eigenvalue = min(
+            entry["lambda"] for entry in self.spectrum if "lambda" in entry
+        )
         # Determine if the evolution is unique (example condition)
         unique_evolution = all(entry["lambda"] > 0 for entry in self.spectrum)
 
         # Size of the computed spectrum
         spectrum_size = len(self.spectrum)
 
         # Log the information
-        logger.info(f"Processed eigenvalues: {self.get_number_of_process_eigenvalues(self.eigen)}")
+        logger.info(
+            f"Processed eigenvalues: {self.get_number_of_process_eigenvalues(self.eigen)}"
+        )
         logger.info(f"Inertia: {self.get_inertia()}")
-        logger.info(f"Eigenvalues: {' '.join([f'{value:.1e}' for value in self.data['eigs']])}")
-
+        logger.info(
+            f"Eigenvalues: {' '.join([f'{value:.1e}' for value in self.data['eigs']])}"
+        )
+
         logger.info(f"Minimum eigenvalue: {min_eigenvalue:.2f}")
         logger.info(f"Unique evolution: {unique_evolution}")
         logger.info(f"Size of computed spectrum: {spectrum_size}")
 
+
 class StabilitySolver(SecondOrderSolver):
     """Base class for a minimal implementation of the solution of eigenvalue
     problems bound to a cone. Based on numerical recipe SPA and KR existence result
@@ -1015,8 +1026,8 @@ def _convergenceTest(self, x, errors):
             )
             if self.iterations > 0:
                 _logger.critical(
-                f"     [i={self.iterations}] lambda_k = {self.data['lambda_k'].pop():.2e}, atol = {_atol}, res = {self._residual_norm}"
-            )
+                    f"     [i={self.iterations}] lambda_k = {self.data['lambda_k'].pop():.2e}, atol = {_atol}, res = {self._residual_norm}"
+                )
 
         # self.data["iterations"].append(self.iterations)
         # self.data["error_x_L2"].append(error_x_L2)
@@ -1167,13 +1178,15 @@ def log(self, logger=logger):
         # for key, value in self.data.items():
         #     logger.info(f"{key}: {value}")
         #  = {"lambda_t": lmbda_t, "xt": xt, "yt": yt}
-        if self.solution['lambda_t'] is not np.nan:
+        if self.solution["lambda_t"] is not np.nan:
             logger.info(f"Restricted Eigenvalue: {self.solution['lambda_t']}")
-            logger.info(f"Restricted Eigenfunction is in cone 🍦 ? {self._isin_cone(self.solution['xt'])}")
+            logger.info(
+                f"Restricted Eigenfunction is in cone 🍦 ? {self._isin_cone(self.solution['xt'])}"
+            )
             logger.info(f"Restricted Error {self.error:.4e}")
 
-        return 
-    
+        return
+
     def save_input_data(self, filename="data/input_data.xdmf"):
         """
         Save input data to a file.

src/irrevolutions/solvers/__init__.py

@@ -1,11 +1,18 @@
 from mpi4py import MPI
 import sys
 import petsc4py
+
 petsc4py.init(sys.argv)
 from petsc4py import PETSc
 
-from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
-                               create_matrix, create_vector, set_bc)
+from dolfinx.fem.petsc import (
+    apply_lifting,
+    assemble_matrix,
+    assemble_vector,
+    create_matrix,
+    create_vector,
+    set_bc,
+)
 from dolfinx.cpp.log import LogLevel, log
 
 import dolfinx
@@ -19,6 +26,7 @@
 
 comm = MPI.COMM_WORLD
 
+
 class SNESSolver:
     """
     Problem class for elasticity, compatible with PETSC.SNES solvers.

test/parameters.yml

@@ -7,7 +7,7 @@ loading:
 #  === Geometry === #
 geometry:
     geometric_dimension: 2
-    geom_type: "bar"
+    geom_type: 'bar'
     Lx: 1.
     Ly: .1
     lc: 0.02
@@ -20,7 +20,7 @@ model:
     nu: 0.3
     ell: 0.1
     model_dimension: 2
-    model_type: "2D"
+    model_type: '2D'
     # could be "2D"/ "3D" / "plane_stress" / "plane_strain"
 # === Solver === #
 solvers:
@@ -33,7 +33,7 @@ solvers:
             snes_rtol: 1e-8
             snes_max_it: 200
             # snes_divergence_tolerance: -1.0
-            snes_monitor: ""
+            snes_monitor: ''
             ksp_type: preonly
             pc_type: lu
             pc_factor_mat_solver_type: mumps
@@ -54,7 +54,7 @@ solvers:
             # snes_stol: 0.0
             snes_max_it: 50
             # snes_divergence_tolerance: -1.0
-            snes_monitor: ""
+            snes_monitor: ''
         tao:
             # Options in the case of TAO solver
             tao_type: tron
@@ -73,56 +73,57 @@ solvers:
             tao_ls_stepmin: 1e-8
             tao_ls_stepmax: 1e6
             pc_type: lu
-            tao_monitor: ""
+            tao_monitor: ''
 
     # Damage Elasticity Solver parameters
     damage_elasticity:
         max_it: 200
         alpha_rtol: 1.0e-4
-        criterion: "alpha_H1"
-    
+        criterion: 'alpha_H1'
+
     newton:
-        snes_type: "vinewtonrsls"
-        snes_linesearch_type: "basic"
+        snes_type: 'vinewtonrsls'
+        snes_linesearch_type: 'basic'
         snes_rtol: 1.0e-8
         snes_atol: 1.0e-8
         snes_max_it: 30
-        snes_monitor: ""
+        snes_monitor: ''
         linesearch_damping: 0.5
 
 stability:
     order: 3
     maxmodes: 10
-    checkstability: "True"
-    continuation: "False"
+    checkstability: 'True'
+    continuation: 'False'
     cont_rtol: 1.0e-10
     inactiveset_gatol: 1.e-6
     inactiveset_pwtol: 1.e-6
     is_elastic_tol: 1.e-6
 
     inertia:
         # MUMPS
-        ksp_type: "preonly"
-        pc_type: "cholesky"
-        pc_factor_mat_solver_type: "mumps"
+        ksp_type: 'preonly'
+        pc_type: 'cholesky'
+        pc_factor_mat_solver_type: 'mumps'
         mat_mumps_icntl_24: 1
         mat_mumps_icntl_13: 1
 
     eigen:
-        eps_type: "krylovschur"
+        eps_type: 'krylovschur'
         # eps_type: "lanczos"
         # eps_monitor: ""
         eps_tol: 1.e-5
         eig_rtol: 1.e-8
         eps_max_it: 100
+        shift: -1.0e-3
 
     cone:
         cone_atol: 1.0e-06
         cone_max_it: 3000
         cone_rtol: 1.0e-06
         maxmodes: 3
         scaling: 1.e-5
-    
+
     linesearch:
         order: 4
         method: 'min'