From dfde3fcd69dfa0af2bc4c856d474c5509f709316 Mon Sep 17 00:00:00 2001 From: KulaginVladimir Date: Wed, 3 Apr 2024 02:52:04 +0300 Subject: [PATCH] tests for the newton_solver attributes --- festim/concentration/traps/extrinsic_trap.py | 4 +- festim/nonlinear_problem.py | 6 +- festim/settings.py | 2 +- festim/temperature/temperature_solver.py | 10 ++-- test/h_transport_problem/test_solving.py | 60 +++++++++++++++++-- test/heat_transfer_problem/test_solving.py | 55 ++++++++++++++++- test/system/test_misc.py | 44 ++++++++++++++ test/unit/test_traps/test_extrinsic_trap.py | 3 +- .../test_traps/test_neutron_induced_trap.py | 2 +- 9 files changed, 167 insertions(+), 19 deletions(-) diff --git a/festim/concentration/traps/extrinsic_trap.py b/festim/concentration/traps/extrinsic_trap.py index b63d76d78..fdefcd4b6 100644 --- a/festim/concentration/traps/extrinsic_trap.py +++ b/festim/concentration/traps/extrinsic_trap.py @@ -40,7 +40,7 @@ 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 @@ -48,8 +48,8 @@ def __init__( 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 diff --git a/festim/nonlinear_problem.py b/festim/nonlinear_problem.py index cf8234bc9..06c77fadc 100644 --- a/festim/nonlinear_problem.py +++ b/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) \ No newline at end of file + self.assembler.assemble(A) diff --git a/festim/settings.py b/festim/settings.py index 820f03f43..57b02d3b0 100644 --- a/festim/settings.py +++ b/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 diff --git a/festim/temperature/temperature_solver.py b/festim/temperature/temperature_solver.py index be54cb64d..069bd21d3 100644 --- a/festim/temperature/temperature_solver.py +++ b/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 diff --git a/test/h_transport_problem/test_solving.py b/test/h_transport_problem/test_solving.py index a93e08794..27f5d14d9 100644 --- a/test/h_transport_problem/test_solving.py +++ b/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(): - \ No newline at end of file + +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() diff --git a/test/heat_transfer_problem/test_solving.py b/test/heat_transfer_problem/test_solving.py index 7ceaf396f..b4d68ee25 100644 --- a/test/heat_transfer_problem/test_solving.py +++ b/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) \ No newline at end of file + 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() diff --git a/test/system/test_misc.py b/test/system/test_misc.py index ec3e3804e..ad73ea3e9 100644 --- a/test/system/test_misc.py +++ b/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() diff --git a/test/unit/test_traps/test_extrinsic_trap.py b/test/unit/test_traps/test_extrinsic_trap.py index a314f6f84..5eb44c92f 100644 --- a/test/unit/test_traps/test_extrinsic_trap.py +++ b/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() diff --git a/test/unit/test_traps/test_neutron_induced_trap.py b/test/unit/test_traps/test_neutron_induced_trap.py index 3ad301c32..1d6a441ac 100644 --- a/test/unit/test_traps/test_neutron_induced_trap.py +++ b/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):