tests for the newton_solver attributes

festim/concentration/traps/extrinsic_trap.py

@@ -40,16 +40,16 @@ def __init__(
                 Defaults to None.
             preconditioner (str, optional): preconditioning method for the newton solver,
                 options can be veiwed by print(list_krylov_solver_preconditioners()).
-                Defaults to None.
+                Defaults to "default".
         """
         super().__init__(k_0, E_k, p_0, E_p, materials, density=None, id=id)
         self.absolute_tolerance = absolute_tolerance
         self.relative_tolerance = relative_tolerance
         self.maximum_iterations = maximum_iterations
         self.linear_solver = linear_solver
         self.preconditioner = preconditioner
-        self.newton_solver = None
 
+        self.newton_solver = None
         for name, val in kwargs.items():
             setattr(self, name, as_constant_or_expression(val))
         self.density_previous_solution = None

festim/nonlinear_problem.py

@@ -17,7 +17,9 @@ def __init__(self, J, F, bcs):
         self.jacobian_form = J
         self.residual_form = F
         self.bcs = bcs
-        self.assembler = f.SystemAssembler(self.jacobian_form, self.residual_form, self.bcs)
+        self.assembler = f.SystemAssembler(
+            self.jacobian_form, self.residual_form, self.bcs
+        )
         f.NonlinearProblem.__init__(self)
 
     def F(self, b, x):
@@ -26,4 +28,4 @@ def F(self, b, x):
 
     def J(self, A, x):
         """Assembles the LHS in Ax=b and applies the boundary conditions"""
-        self.assembler.assemble(A)
+        self.assembler.assemble(A)

festim/settings.py

@@ -27,7 +27,7 @@ class Settings:
             Defaults to None, for the newton solver this is: "umfpack".
         preconditioner (str, optional): preconditioning method for the newton solver,
             options can be veiwed by print(list_krylov_solver_preconditioners()).
-            Defaults to None.
+            Defaults to "default".
 
     Attributes:
         transient (bool): transient or steady state sim

festim/temperature/temperature_solver.py

@@ -23,7 +23,7 @@ class HeatTransferProblem(festim.Temperature):
             Defaults to None.
         preconditioner (str, optional): preconditioning method for the newton solver,
             options can be veiwed by print(list_krylov_solver_preconditioners()).
-            Defaults to None.
+            Defaults to "default".
 
     Attributes:
         F (fenics.Form): the variational form of the heat transfer problem
@@ -45,7 +45,7 @@ def __init__(
         relative_tolerance=1e-10,
         maximum_iterations=30,
         linear_solver=None,
-        preconditioner=None,
+        preconditioner="default",
     ) -> None:
         super().__init__()
         self.transient = transient
@@ -97,7 +97,9 @@ def create_functions(self, materials, mesh, dt=None):
 
         self.define_variational_problem(materials, mesh, dt)
         self.create_dirichlet_bcs(mesh.surface_markers)
-        self.define_newton_solver_temperature()
+
+        if not self.newton_solver:
+            self.define_newton_solver()
 
         if not self.transient:
             print("Solving stationary heat equation")
@@ -183,7 +185,7 @@ def define_variational_problem(self, materials, mesh, dt=None):
                 for surf in bc.surfaces:
                     self.F += -bc.form * self.v_T * mesh.ds(surf)
 
-    def define_newton_solver_temperature(self):
+    def define_newton_solver(self):
         """Creates the Newton solver and sets its parameters"""
         self.newton_solver = f.NewtonSolver(f.MPI.comm_world)
         self.newton_solver.parameters["error_on_nonconvergence"] = False

test/h_transport_problem/test_solving.py

@@ -98,12 +98,12 @@ def test_solve_once_returns_false():
     assert not converged
 
 
-@pytest.mark.parametrize("preconditioner",["default", "icc"])
+@pytest.mark.parametrize("preconditioner", ["default", "icc"])
 def test_solve_once_linear_solver_gmres(preconditioner):
     """
-    Checks that solve_once() works when an alternative linear solver is used 
+    Checks that solve_once() works when an alternative linear solver is used
     with/without a preconditioner rather than the default
-    
+
     Args:
         preconditioner (str): the preconditioning method
     """
@@ -136,5 +136,55 @@ def test_solve_once_linear_solver_gmres(preconditioner):
     # test
     assert converged
 
-    #def test_solve_once_with_custom_solver():
-
+
+class Test_solve_once_with_custom_solver:
+    """
+    Checks that a custom newton sovler can be used
+    """
+
+    def sim(self):
+        """Defines a model"""
+        mesh = f.UnitIntervalMesh(8)
+        V = f.FunctionSpace(mesh, "CG", 1)
+
+        my_settings = festim.Settings(
+            absolute_tolerance=1e-10,
+            relative_tolerance=1e-10,
+            maximum_iterations=50,
+        )
+        my_problem = festim.HTransportProblem(
+            festim.Mobile(), festim.Traps([]), festim.Temperature(200), my_settings, []
+        )
+        my_problem.define_newton_solver()
+        my_problem.u = f.Function(V)
+        my_problem.u_n = f.Function(V)
+        my_problem.v = f.TestFunction(V)
+        my_problem.F = (
+            (my_problem.u - my_problem.u_n) * my_problem.v * f.dx
+            + 1 * my_problem.v * f.dx
+            + f.dot(f.grad(my_problem.u), f.grad(my_problem.v)) * f.dx
+        )
+        return my_problem
+
+    def test_custom_solver(self):
+        """Solves the system using the built-in solver and using the f.NewtonSolver"""
+
+        # solve with the built-in solver
+        problem_1 = self.sim()
+        problem_1.solve_once()
+
+        # solve with the custom solver
+        problem_2 = self.sim()
+        problem_2.newton_solver = f.NewtonSolver()
+        problem_2.newton_solver.parameters[
+            "absolute_tolerance"
+        ] = problem_1.settings.absolute_tolerance
+        problem_2.newton_solver.parameters[
+            "relative_tolerance"
+        ] = problem_1.settings.relative_tolerance
+        problem_2.newton_solver.parameters[
+            "maximum_iterations"
+        ] = problem_1.settings.maximum_iterations
+        problem_2.solve_once()
+
+        assert (problem_1.u.vector() == problem_2.u.vector()).all()

test/heat_transfer_problem/test_solving.py

@@ -3,12 +3,12 @@
 import fenics as f
 
 
-@pytest.mark.parametrize("preconditioner",["default", "icc"])
+@pytest.mark.parametrize("preconditioner", ["default", "icc"])
 def test_create_functions_linear_solver_gmres(preconditioner):
     """
     Checks that the function created by create_functions() has the expected value when an
     alternative linear solver is used with/without a preconditioner rather than the default
-    
+
     Args:
         preconditioner (str): the preconditioning method
     """
@@ -35,4 +35,53 @@ def test_create_functions_linear_solver_gmres(preconditioner):
     # run
     my_problem.create_functions(materials=materials, mesh=mesh)
 
-    assert my_problem.T(0.05) == pytest.approx(1)
+    assert my_problem.T(0.05) == pytest.approx(1)
+
+
+class Test_solve_once_with_custom_solver:
+    """
+    Checks that a custom newton sovler can be used
+    """
+
+    def sim(self):
+        """Defines a model"""
+        bcs = [
+            festim.DirichletBC(surfaces=[1, 2], value=1, field="T"),
+        ]
+
+        my_problem = festim.HeatTransferProblem(
+            transient=False,
+            absolute_tolerance=1e-03,
+            relative_tolerance=1e-10,
+            maximum_iterations=30,
+        )
+        my_problem.boundary_conditions = bcs
+        return my_problem
+
+    def test_custom_solver(self):
+        """Solves the system using the built-in solver and using the f.NewtonSolver"""
+        mesh = festim.MeshFromRefinements(10, size=0.1)
+        materials = festim.Materials(
+            [festim.Material(id=1, D_0=1, E_D=0, thermal_cond=1)]
+        )
+        mesh.define_measures(materials)
+
+        # solve with the built-in solver
+        problem_1 = self.sim()
+        problem_1.create_functions(materials=materials, mesh=mesh)
+
+        # solve with the custom solver
+        problem_2 = self.sim()
+        problem_2.newton_solver = f.NewtonSolver()
+        problem_2.newton_solver.parameters[
+            "absolute_tolerance"
+        ] = problem_1.absolute_tolerance
+        problem_2.newton_solver.parameters[
+            "relative_tolerance"
+        ] = problem_1.relative_tolerance
+        problem_2.newton_solver.parameters[
+            "maximum_iterations"
+        ] = problem_1.maximum_iterations
+        problem_2.create_functions(materials=materials, mesh=mesh)
+
+        assert (problem_1.T.vector() == problem_2.T.vector()).all()

test/system/test_misc.py

@@ -1,4 +1,5 @@
 import festim as F
+import fenics as f
 import numpy as np
 import pytest
 import os
@@ -330,3 +331,46 @@ def test_materials_setter():
     test_materials = F.Materials([])
     my_model.materials = test_materials
     assert my_model.materials is test_materials
+
+
+class TestFestimProblem:
+    """Tests the methods of the festim.Problem class"""
+
+    # define the FESTIM problem
+    mesh = f.UnitIntervalMesh(8)
+    V = f.FunctionSpace(mesh, "CG", 1)
+    u = f.Function(V)
+    v = f.TestFunction(V)
+    s = u * v * f.dx + v * f.dx
+    J = f.derivative(s, u)
+    problem = F.Problem(J, s, [])
+
+    # define the fenics assembler used in festim.Problem
+    assembler = f.SystemAssembler(J, s, [])
+    x = f.PETScVector()
+
+    def test_F(self):
+        """
+        Creates two epty matrices and checks
+        that the festim.Problem.J properly assembles the RHS of Ax=b
+        """
+        b1 = f.PETScVector()
+        b2 = f.PETScVector()
+
+        self.problem.F(b1, self.x)
+        self.assembler.assemble(b2, self.x)
+
+        assert (b1 == b2).all()
+
+    def test_J(self):
+        """
+        Creates two epty matrices and checks
+        that the festim.Problem.J properly assembles the LHS of Ax=b
+        """
+        A1 = f.PETScMatrix()
+        A2 = f.PETScMatrix()
+
+        self.problem.J(A1, self.x)
+        self.assembler.assemble(A2)
+
+        assert (A1.array() == A2.array()).all()

test/unit/test_traps/test_extrinsic_trap.py

@@ -87,4 +87,5 @@ def test_solver_parameters(self):
         self.my_trap.relative_tolerance = 1
         self.my_trap.maximum_iterations = 1
         self.my_trap.linear_solver = "mumps"
-        self.my_trap.preconditioner = 'icc'
+        self.my_trap.preconditioner = "icc"
+        self.my_trap.newton_solver = f.NewtonSolver()

test/unit/test_traps/test_neutron_induced_trap.py

@@ -147,7 +147,7 @@ class TestNeutronInducedTrapSolverParameters:
         relative_tolerance=1.2e-10,
         maximum_iterations=13,
         linear_solver="gmres",
-        preconditioner="icc"
+        preconditioner="icc",
     )
 
     def test_attributes_from_instanciation(self):