added new default solver

festim/problem.py

@@ -24,6 +24,7 @@ def __init__(
         subdomains=None,
         boundary_conditions=None,
         settings=None,
+        petcs_options=None,
     ) -> None:
         self.mesh = mesh
         # for arguments to initialise as empty list
@@ -42,6 +43,7 @@ def __init__(
         self.formulation = None
         self.bc_forms = []
         self.show_progress_bar = True
+        self.petcs_options = petcs_options
 
     @property
     def volume_subdomains(self):
@@ -102,6 +104,13 @@ def create_solver(self):
         self.solver.rtol = self.settings.rtol
         self.solver.max_it = self.settings.max_iterations
 
+        if self.petcs_options is None:
+            ksp = self.solver.krylov_solver
+            ksp.setType("preonly")
+            ksp.getPC().setType("lu")
+            ksp.getPC().setFactorSolverType("mumps")
+            ksp.setErrorIfNotConverged(True)
+
     def run(self):
         """Runs the model"""