From 61a28c49368dcbf33678890cbab0909b2f59c0bc Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 22 Apr 2026 14:15:35 -0400 Subject: [PATCH 01/14] added tests --- test/system_tests/test_multi_material.py | 49 ++++++++++++++++++++++++ test/test_subdomains.py | 20 ++++++++++ 2 files changed, 69 insertions(+) diff --git a/test/system_tests/test_multi_material.py b/test/system_tests/test_multi_material.py index 86526d7b4..3a7c1551a 100644 --- a/test/system_tests/test_multi_material.py +++ b/test/system_tests/test_multi_material.py @@ -3,6 +3,7 @@ import dolfinx import dolfinx.fem.petsc import numpy as np +import pytest import ufl import festim as F @@ -396,3 +397,51 @@ def test_2_mats_particle_flux_bc(tmpdir): my_model.initialise() my_model.run() + + +def test_all_cells_are_not_tagged(): + """ + Checks that an error is raised when not all cells are tagged with a non-zero value. + This can be caused by a volume subdomain not being defined correctly, or by a mesh + that is too coarse to capture the geometry of the volume subdomains. + """ + + my_model = F.HydrogenTransportProblemDiscontinuous() + + mat = F.Material(D_0=1, E_D=0, K_S_0=1, E_K_S=0) + + vol1 = F.VolumeSubdomain1D(id=1, borders=[0, 0.5], material=mat) + vol2 = F.VolumeSubdomain1D(id=2, borders=[0.5, 1], material=mat) + + surf1 = F.SurfaceSubdomain1D(id=1, x=0) + surf2 = F.SurfaceSubdomain1D(id=2, x=1) + + my_model.subdomains = [vol1, vol2, surf1, surf2] + + my_model.surface_to_volume = { + surf1: vol1, + surf2: vol2, + } + my_model.interfaces = [F.Interface(id=3, subdomains=[vol1, vol2])] + + my_model.mesh = F.Mesh1D(np.linspace(0, 1, 10)) + + my_model.species = [F.Species("H", subdomains=[vol1, vol2])] + + my_model.boundary_conditions = [ + F.FixedConcentrationBC(species=my_model.species[0], subdomain=surf1, value=1), + ] + + my_model.temperature = 300 + + my_model.settings = F.Settings(transient=False, atol=1e-9, rtol=1e-9) + + my_model.exports = [ + F.AverageVolume(field=my_model.species[0], volume=vol1), + F.AverageVolume(field=my_model.species[0], volume=vol2), + ] + + with pytest.raises( + AssertionError, match="All cells must be tagged with a non-zero value" + ): + my_model.initialise() diff --git a/test/test_subdomains.py b/test/test_subdomains.py index 9f87ab145..757a9af01 100644 --- a/test/test_subdomains.py +++ b/test/test_subdomains.py @@ -214,3 +214,23 @@ def test_name_setter(): with pytest.raises(TypeError, match="Name must be a string"): F.VolumeSubdomain(id=1, material=None, name=1) + + +def test_all_cells_are_not_tagged(): + """ + Checks that an error is raised when not all cells are tagged with a non-zero value. + This can be caused by a volume subdomain not being defined correctly, or by a mesh + that is too coarse to capture the geometry of the volume subdomains. + """ + + mat = F.Material(D_0=1, E_D=0, K_S_0=1, E_K_S=0) + + vol1 = F.VolumeSubdomain1D(id=1, borders=[0, 0.5], material=mat) + vol2 = F.VolumeSubdomain1D(id=2, borders=[0.5, 1], material=mat) + + mesh = F.Mesh1D(np.linspace(0, 1, 10)) + + with pytest.raises( + AssertionError, match="All cells must be tagged with a non-zero value" + ): + mesh.define_meshtags(surface_subdomains=[], volume_subdomains=[vol1, vol2]) From cbf5064b34b17ec0e639e5e6a269e854d1266b2a Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 22 Apr 2026 14:16:22 -0400 Subject: [PATCH 02/14] added check in define_meshtags --- src/festim/mesh/mesh.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/festim/mesh/mesh.py b/src/festim/mesh/mesh.py index 032897c7a..d08d4cab6 100644 --- a/src/festim/mesh/mesh.py +++ b/src/festim/mesh/mesh.py @@ -180,6 +180,14 @@ def define_meshtags(self, surface_subdomains, volume_subdomains, interfaces=None ) tags_volumes[:] = vol.id + # make sure there are no zeros in the tags_volumes + assert np.all(tags_volumes != 0), ( + "All cells must be tagged with a non-zero value" + "This error can be caused by a volume subdomain not being defined" + "correctly, or by a mesh that is too coarse to capture the geometry" + "of the volume subdomains" + ) + volume_meshtags = meshtags( self._mesh, self.vdim, mesh_cell_indices, tags_volumes ) From 7e64f6375e144a2dfbf674ab83bbdcd5a06838e1 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 22 Apr 2026 14:17:11 -0400 Subject: [PATCH 03/14] added an assertion to avoid id=0 --- src/festim/subdomain/volume_subdomain.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/festim/subdomain/volume_subdomain.py b/src/festim/subdomain/volume_subdomain.py index 56bdbd806..22a735b00 100644 --- a/src/festim/subdomain/volume_subdomain.py +++ b/src/festim/subdomain/volume_subdomain.py @@ -24,7 +24,7 @@ class VolumeSubdomain: Volume subdomain class Args: - id: the id of the volume subdomain + id: the id of the volume subdomain (> 0) submesh: the submesh of the volume subdomain cell_map: the cell map of the volume subdomain parent_mesh: the parent mesh of the volume subdomain @@ -58,6 +58,7 @@ class VolumeSubdomain: def __init__( self, id, material, locator: Callable | None = None, name: str | None = None ): + assert id != 0, "Volume subdomain id cannot be 0" self.id = id self.material = material self.locator = locator From 6632f995b5dc1ec4a7561fe51b70c98b5cd5934a Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 22 Apr 2026 14:21:17 -0400 Subject: [PATCH 04/14] fixed test --- test/system_tests/test_io.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/system_tests/test_io.py b/test/system_tests/test_io.py index ffae401a8..f0eceb58a 100644 --- a/test/system_tests/test_io.py +++ b/test/system_tests/test_io.py @@ -19,7 +19,7 @@ def test_writing_and_reading_of_species_function_using_checkpoints(tmpdir): my_model.mesh = F.Mesh(mesh) my_mat = F.Material(name="mat", D_0=1, E_D=0) - vol = F.VolumeSubdomain(id=0, material=my_mat) + vol = F.VolumeSubdomain(id=1, material=my_mat) surf = F.SurfaceSubdomain(id=1) my_model.subdomains = [vol, surf] From 649e10bd5311985fc7b492b7818aab6544774d60 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 22 Apr 2026 14:24:13 -0400 Subject: [PATCH 05/14] fixed more tests --- test/test_heat_transfer_problem.py | 2 +- test/test_initial_condition.py | 4 ++-- test/test_subdomains.py | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/test/test_heat_transfer_problem.py b/test/test_heat_transfer_problem.py index 39210a320..c4e3dd5f1 100644 --- a/test/test_heat_transfer_problem.py +++ b/test/test_heat_transfer_problem.py @@ -484,7 +484,7 @@ def test_meshtags_from_xdmf(tmp_path, mesh): def test_raise_error_non_unique_vol_ids(): """Test that an error is raised if the volume ids are not unique""" my_problem = F.HeatTransferProblem() - my_problem.mesh = F.Mesh1D(vertices=np.linspace(0, 1, 2000)) + my_problem.mesh = F.Mesh1D(vertices=np.linspace(0, 1, 11)) left = F.SurfaceSubdomain1D(id=1, x=0) right = F.SurfaceSubdomain1D(id=2, x=1) mat = F.Material(D_0=1, E_D=None) diff --git a/test/test_initial_condition.py b/test/test_initial_condition.py index 6ced0b703..588785445 100644 --- a/test/test_initial_condition.py +++ b/test/test_initial_condition.py @@ -169,7 +169,7 @@ def f(x): my_problem.mesh = F.Mesh(mesh) mat = F.Material(D_0=1, E_D=0.1, name="dummy_mat") - vol = F.VolumeSubdomain(id=0, material=mat) + vol = F.VolumeSubdomain(id=1, material=mat) my_problem.subdomains = [vol] function_initial_value = F.read_function_from_file( @@ -227,7 +227,7 @@ def f2(x): my_problem.mesh = F.Mesh(mesh) mat = F.Material(D_0=1, E_D=0.1, name="dummy_mat") - vol = F.VolumeSubdomain(id=0, material=mat) + vol = F.VolumeSubdomain(id=1, material=mat) my_problem.subdomains = [vol] my_problem.initial_conditions = [ diff --git a/test/test_subdomains.py b/test/test_subdomains.py index 757a9af01..351ac0e42 100644 --- a/test/test_subdomains.py +++ b/test/test_subdomains.py @@ -125,7 +125,7 @@ def test_ValueError_raised_when_id_not_found_in_volumes_subdomains(): def test_ValueError_rasied_when_volume_ids_are_not_unique(): """Checks""" my_test_model = F.HydrogenTransportProblem( - mesh=F.Mesh1D(np.linspace(0, 2, num=10)), species=[F.Species("H")] + mesh=F.Mesh1D(np.linspace(0, 2, num=11)), species=[F.Species("H")] ) vol_1 = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=None) From 630e28106de0dd59131cd6e6ed2ef4b7d88276a7 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 22 Apr 2026 14:25:07 -0400 Subject: [PATCH 06/14] ruff fixes --- src/festim/hydrogen_transport_problem.py | 1 - src/festim/problem.py | 2 +- test/test_reaction.py | 2 +- 3 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 14d9c35b4..58e196c62 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -41,7 +41,6 @@ is_it_time_to_export, nmm_interpolate, ) -from festim.material import SolubilityLaw from .mesh import CoordinateSystem, Mesh diff --git a/src/festim/problem.py b/src/festim/problem.py index 15901c797..5694d23f0 100644 --- a/src/festim/problem.py +++ b/src/festim/problem.py @@ -1,3 +1,4 @@ +import warnings from typing import Any from mpi4py import MPI @@ -10,7 +11,6 @@ from dolfinx import fem from dolfinx.nls.petsc import NewtonSolver from packaging.version import Version -import warnings import festim as F from festim.mesh.mesh import Mesh as _Mesh diff --git a/test/test_reaction.py b/test/test_reaction.py index a34861763..0dd24e331 100644 --- a/test/test_reaction.py +++ b/test/test_reaction.py @@ -1,10 +1,10 @@ from mpi4py import MPI +import numpy as np import pytest from dolfinx.fem import Function, functionspace from dolfinx.mesh import create_unit_cube from ufl import exp -import numpy as np import festim as F From 4d22a3fe5946848329d9ba136fe61ea0e6f94a9d Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 22 Apr 2026 14:26:50 -0400 Subject: [PATCH 07/14] unsafe fixes --- .../advection_diffusion_multi_material.py | 5 +-- examples/multi_material_2d.py | 7 +-- examples/multi_material_2d_bis.py | 7 +-- src/festim/helpers.py | 10 +++-- src/festim/hydrogen_transport_problem.py | 2 +- src/festim/problem.py | 2 +- src/festim/reaction.py | 12 ++--- test/system_tests/test_1_mobile_1_trap.py | 18 +++++--- test/system_tests/test_1_mobile_species.py | 18 +++++--- .../test_advection_term_integration.py | 13 +++--- .../test_coupled_heat_hydrogen.py | 45 ++++++++++++------- .../test_cylindrical_coordinates.py | 12 +++-- test/system_tests/test_fluxes.py | 6 ++- test/system_tests/test_multi_mat_penalty.py | 15 ++++--- test/system_tests/test_multi_material.py | 29 ++++++------ .../test_multi_material_change_of_var.py | 29 ++++++------ .../test_non_homogeneous_density.py | 13 ++++-- .../test_spherical_coordinates.py | 14 +++--- test/test_average_surface.py | 2 +- test/test_average_volume.py | 2 +- test/test_h_transport_problem.py | 12 ++--- test/test_helpers.py | 3 +- test/test_initial_condition.py | 21 +++++---- test/test_maximum_surface.py | 2 +- test/test_maximum_volume.py | 2 +- test/test_minimum_surface.py | 2 +- test/test_minimum_volume.py | 2 +- test/test_species.py | 3 +- test/test_total_surface.py | 2 +- test/test_total_volume.py | 2 +- 30 files changed, 186 insertions(+), 126 deletions(-) diff --git a/examples/advection_diffusion_multi_material.py b/examples/advection_diffusion_multi_material.py index a1869335a..c56ad8389 100644 --- a/examples/advection_diffusion_multi_material.py +++ b/examples/advection_diffusion_multi_material.py @@ -27,7 +27,6 @@ def half(x): ) # Split domain in half and set an interface tag of 5 - gdim = mesh.geometry.dim tdim = mesh.topology.dim fdim = tdim - 1 top_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, top_boundary) @@ -80,8 +79,8 @@ def velocity_func(x): values[1] = 0 # Second component remains zero return values - mesh2, mt2, ct2 = generate_mesh(n=10) - submesh, cell_map, v_map = dolfinx.mesh.create_submesh(mesh2, ct2.dim, ct2.find(4))[ + mesh2, _mt2, ct2 = generate_mesh(n=10) + submesh, _cell_map, _v_map = dolfinx.mesh.create_submesh(mesh2, ct2.dim, ct2.find(4))[ 0:3 ] v_cg = basix.ufl.element( diff --git a/examples/multi_material_2d.py b/examples/multi_material_2d.py index 7a64cbe8c..314d9ef32 100644 --- a/examples/multi_material_2d.py +++ b/examples/multi_material_2d.py @@ -24,7 +24,6 @@ def half(x): ) # Split domain in half and set an interface tag of 5 - gdim = mesh.geometry.dim tdim = mesh.topology.dim fdim = tdim - 1 top_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, top_boundary) @@ -125,11 +124,13 @@ def half(x): from dolfinx.mesh import EntityMap # noqa: F401 legacy_entity_map = False - entity_map_wrapper = lambda e_map: list(e_map.values()) + def entity_map_wrapper(e_map): + return list(e_map.values()) entity_maps = [sd.cell_map for sd in my_model.volume_subdomains] except ImportError: legacy_entity_map = True - entity_map_wrapper = lambda e_map: e_map + def entity_map_wrapper(e_map): + return e_map entity_maps = { sd.submesh: sd.parent_to_submesh for sd in my_model.volume_subdomains } diff --git a/examples/multi_material_2d_bis.py b/examples/multi_material_2d_bis.py index cc6ad33da..6ea69169b 100644 --- a/examples/multi_material_2d_bis.py +++ b/examples/multi_material_2d_bis.py @@ -24,7 +24,6 @@ def half(x): ) # Split domain in half and set an interface tag of 5 - gdim = mesh.geometry.dim tdim = mesh.topology.dim fdim = tdim - 1 top_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, top_boundary) @@ -168,10 +167,12 @@ def half(x): from dolfinx.mesh import EntityMap # noqa: F401 legacy_entity_map = False - entity_map_wrapper = lambda e_map: list(e_map.values()) + def entity_map_wrapper(e_map): + return list(e_map.values()) except ImportError: legacy_entity_map = True - entity_map_wrapper = lambda e_map: e_map + def entity_map_wrapper(e_map): + return e_map form = dolfinx.fem.form( diff --git a/src/festim/helpers.py b/src/festim/helpers.py index 0a1ba5290..94b8ab41b 100644 --- a/src/festim/helpers.py +++ b/src/festim/helpers.py @@ -285,9 +285,11 @@ def update(self, t: float): # Define the appropriate method based on the version if version.parse(dolfinx_version) > version.parse("0.9.0"): - get_interpolation_points = lambda element: element.interpolation_points + def get_interpolation_points(element): + return element.interpolation_points else: - get_interpolation_points = lambda element: element.interpolation_points() + def get_interpolation_points(element): + return element.interpolation_points() def nmm_interpolate( @@ -356,7 +358,7 @@ def is_it_time_to_export( def convergenceTest(snes, it, norms): global _residual0 - xnorm, gnorm, f = norms # ||x_k||, ||x_k-x_k-1||, ||F(x_k)|| + _xnorm, gnorm, f = norms # ||x_k||, ||x_k-x_k-1||, ||F(x_k)|| rtol, atol, stol, max_its = snes.getTolerances() @@ -378,7 +380,7 @@ def convergenceTest(snes, it, norms): def SnesMonitor(snes, iter, rnorm): global _prev_xnorm if MPI.COMM_WORLD.rank == 0: - rtol, atol, stol, max_its = snes.getTolerances() + rtol, atol, stol, _max_its = snes.getTolerances() x = snes.getSolution() xnorm = x.norm() diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 58e196c62..0f30ed0f1 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -1818,7 +1818,7 @@ def iterate(self): # Solve main problem if Version(dolfinx.__version__) == Version("0.9.0"): - nb_its, converged = self.solver.solve() + nb_its, _converged = self.solver.solve() elif Version(dolfinx.__version__) > Version("0.9.0"): _ = self.solver.solve() converged_reason = self.solver.solver.getConvergedReason() diff --git a/src/festim/problem.py b/src/festim/problem.py index 5694d23f0..bfd9852cd 100644 --- a/src/festim/problem.py +++ b/src/festim/problem.py @@ -271,7 +271,7 @@ def iterate(self): # solve main problem if Version(dolfinx.__version__) == Version("0.9.0"): - nb_its, converged = self.solver.solve(self.u) + nb_its, _converged = self.solver.solve(self.u) elif Version(dolfinx.__version__) > Version("0.9.0"): _ = self.solver.solve() converged_reason = self.solver.solver.getConvergedReason() diff --git a/src/festim/reaction.py b/src/festim/reaction.py index 336b18f80..aed8a70eb 100644 --- a/src/festim/reaction.py +++ b/src/festim/reaction.py @@ -63,8 +63,8 @@ def __init__( E_k: float, volume: VS1D, product: Union[_Species, list[_Species]] | None = [], - p_0: float = None, - E_p: float = None, + p_0: float | None = None, + E_p: float | None = None, ) -> None: self.k_0 = k_0 self.E_k = E_k @@ -114,8 +114,8 @@ def __str__(self) -> str: def reaction_term( self, temperature, - reactant_concentrations: list = None, - product_concentrations: list = None, + reactant_concentrations: list | None = None, + product_concentrations: list | None = None, ) -> Expr: """Compute the reaction term at a given temperature. @@ -146,11 +146,11 @@ def reaction_term( + " when no products are present." ) else: - if self.p_0 == None: + if self.p_0 is None: raise ValueError( "p_0 cannot be None when reaction products are present." ) - elif self.E_p == None: + elif self.E_p is None: raise ValueError( "E_p cannot be None when reaction products are present." ) diff --git a/test/system_tests/test_1_mobile_1_trap.py b/test/system_tests/test_1_mobile_1_trap.py index 39d0e43ec..4942bb5b3 100644 --- a/test/system_tests/test_1_mobile_1_trap.py +++ b/test/system_tests/test_1_mobile_1_trap.py @@ -42,7 +42,8 @@ def v_exact(mod): E_k = 1.5 p_0 = 0.2 E_p = 0.1 - T_expr = lambda x: 500 + 100 * x[0] + def T_expr(x): + return 500 + 100 * x[0] T.interpolate(T_expr) n_trap = 3 E_D = 0.1 @@ -138,7 +139,8 @@ def u_exact_alt(mod): D_0 = 1 E_D = 0.1 - T_expr = lambda x: 600 + 50 * x[0] + def T_expr(x): + return 600 + 50 * x[0] T.interpolate(T_expr) D = D_0 * ufl.exp(-E_D / (F.k_B * T)) @@ -163,12 +165,14 @@ def u_exact_alt(mod): F.DirichletBC(subdomain=right, value=H_analytical_ufl, species=H), ] - init_value = lambda x: 1 + ufl.sin(2 * ufl.pi * x[0]) + def init_value(x): + return 1 + ufl.sin(2 * ufl.pi * x[0]) my_model.initial_conditions = [ F.InitialConcentration(value=init_value, species=H, volume=vol) ] - f = lambda x, t: 4 * t - ufl.div(D * ufl.grad(H_analytical_ufl(x_1d, t))) + def f(x, t): + return 4 * t - ufl.div(D * ufl.grad(H_analytical_ufl(x_1d, t))) my_model.sources = [F.ParticleSource(value=f, volume=vol, species=H)] my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=final_time) @@ -208,7 +212,8 @@ def v_exact(mod): E_k = 1.5 p_0 = 0.2 E_p = 0.1 - T_expr = lambda x: 500 + 100 * x[0] + def T_expr(x): + return 500 + 100 * x[0] T.interpolate(T_expr) n_trap = 3 E_D = 0.1 @@ -321,7 +326,8 @@ def v_exact(mod): E_k = 1.5 p_0 = 0.2 E_p = 0.1 - T_expr = lambda x: 500 + 100 * x[0] + def T_expr(x): + return 500 + 100 * x[0] T.interpolate(T_expr) n_trap = 3 E_D = 0.1 diff --git a/test/system_tests/test_1_mobile_species.py b/test/system_tests/test_1_mobile_species.py index 9168a0e3b..84925a54d 100644 --- a/test/system_tests/test_1_mobile_species.py +++ b/test/system_tests/test_1_mobile_species.py @@ -33,7 +33,8 @@ def u_exact(mod): D_0 = 1 E_D = 0.1 - T_expr = lambda x: 500 + 100 * x[0] + def T_expr(x): + return 500 + 100 * x[0] T.interpolate(T_expr) D = D_0 * ufl.exp(-E_D / (F.k_B * T)) @@ -93,7 +94,8 @@ def u_exact_alt(mod): D_0 = 1 E_D = 0.1 - T_expr = lambda x: 600 + 50 * x[0] + def T_expr(x): + return 600 + 50 * x[0] T.interpolate(T_expr) D = D_0 * ufl.exp(-E_D / (F.k_B * T)) @@ -118,12 +120,14 @@ def u_exact_alt(mod): F.DirichletBC(subdomain=right, value=H_analytical_ufl, species=H), ] - init_value = lambda x: 1 + ufl.sin(2 * ufl.pi * x[0]) + def init_value(x): + return 1 + ufl.sin(2 * ufl.pi * x[0]) my_model.initial_conditions = [ F.InitialConcentration(value=init_value, species=H, volume=vol) ] - f = lambda x, t: 4 * t - ufl.div(D * ufl.grad(H_analytical_ufl(x, t))) + def f(x, t): + return 4 * t - ufl.div(D * ufl.grad(H_analytical_ufl(x, t))) my_model.sources = [F.ParticleSource(value=f, volume=vol, species=H)] my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=final_time) @@ -154,7 +158,8 @@ def u_exact(mod): D_0 = 1 E_D = 0.1 - T_expr = lambda x: 500 + 100 * x[0] + def T_expr(x): + return 500 + 100 * x[0] T.interpolate(T_expr) D = D_0 * ufl.exp(-E_D / (F.k_B * T)) @@ -219,7 +224,8 @@ def u_exact(mod): D_0 = 1 E_D = 0.1 - T_expr = lambda x: 500 + 100 * x[0] + def T_expr(x): + return 500 + 100 * x[0] T.interpolate(T_expr) D = D_0 * ufl.exp(-E_D / (F.k_B * T)) diff --git a/test/system_tests/test_advection_term_integration.py b/test/system_tests/test_advection_term_integration.py index 3d3e728a0..2ba14e9de 100644 --- a/test/system_tests/test_advection_term_integration.py +++ b/test/system_tests/test_advection_term_integration.py @@ -38,7 +38,8 @@ def test_MMS_1_species_1_trap_with_advection(): V = fem.functionspace(test_mesh_2d, ("Lagrange", 1)) T = fem.Function(V) - T_expr = lambda x: 100 + 200 * x[0] + 100 * x[1] + def T_expr(x): + return 100 + 200 * x[0] + 100 * x[1] T.interpolate(T_expr) # create velocity field @@ -63,8 +64,10 @@ def velocity_func(x): u.interpolate(velocity_func) # define hydrogen problem - exact_mobile_solution = lambda x: 200 * x[0] ** 2 + 300 * x[1] ** 2 - exact_trapped_solution = lambda x: 10 * x[0] ** 2 + 10 * x[1] ** 2 + def exact_mobile_solution(x): + return 200 * x[0] ** 2 + 300 * x[1] ** 2 + def exact_trapped_solution(x): + return 10 * x[0] ** 2 + 10 * x[1] ** 2 D = D_0 * ufl.exp(-E_D / (k_B * T)) k = k_0 * ufl.exp(-E_k / (k_B * T)) @@ -207,8 +210,8 @@ def velocity_func(x): values[1] = 0 # Second component remains zero return values - mesh2, mt2, ct2 = generate_mesh(n=10) - submesh, cell_map, v_map = dolfinx.mesh.create_submesh( + mesh2, _mt2, ct2 = generate_mesh(n=10) + submesh, _cell_map, _v_map = dolfinx.mesh.create_submesh( mesh2, ct2.dim, ct2.find(4) )[0:3] v_cg = basix.ufl.element( diff --git a/test/system_tests/test_coupled_heat_hydrogen.py b/test/system_tests/test_coupled_heat_hydrogen.py index df6b657df..34122882b 100644 --- a/test/system_tests/test_coupled_heat_hydrogen.py +++ b/test/system_tests/test_coupled_heat_hydrogen.py @@ -38,7 +38,8 @@ def test_MMS_coupled_problem(): test_traps = F.ImplicitSpecies(n=n_trap, others=[test_mobile, test_trapped]) # define temperature sim - exact_T_solution = lambda x, t: 3 * x[0] ** 2 + 10 * t + def exact_T_solution(x, t): + return 3 * x[0] ** 2 + 10 * t dTdt = 10 mms_T_source = ( @@ -67,18 +68,25 @@ def test_MMS_coupled_problem(): ) # define hydrogen problem - exact_mobile_solution = lambda x, t: 2 * x[0] ** 2 + 15 * t - exact_trapped_solution = lambda x, t: 4 * x[0] ** 2 + 12 * t + def exact_mobile_solution(x, t): + return 2 * x[0] ** 2 + 15 * t + def exact_trapped_solution(x, t): + return 4 * x[0] ** 2 + 12 * t - exact_mobile_intial_cond = lambda x: exact_mobile_solution(x, t=0) - exact_trapped_intial_cond = lambda x: exact_trapped_solution(x, t=0) + def exact_mobile_intial_cond(x): + return exact_mobile_solution(x, t=0) + def exact_trapped_intial_cond(x): + return exact_trapped_solution(x, t=0) dmobiledt = 15 dtrappeddt = 12 - D = lambda x, t: D_0 * ufl.exp(-E_D / (k_B * exact_T_solution(x, t))) - k = lambda x, t: k_0 * ufl.exp(-E_k / (k_B * exact_T_solution(x, t))) - p = lambda x, t: p_0 * ufl.exp(-E_p / (k_B * exact_T_solution(x, t))) + def D(x, t): + return D_0 * ufl.exp(-E_D / (k_B * exact_T_solution(x, t))) + def k(x, t): + return k_0 * ufl.exp(-E_k / (k_B * exact_T_solution(x, t))) + def p(x, t): + return p_0 * ufl.exp(-E_p / (k_B * exact_T_solution(x, t))) def mms_mobile_source(x, t): return ( @@ -176,9 +184,12 @@ def mms_trapped_source(x, t): mobile_computed = test_mobile.post_processing_solution trapped_computed = test_trapped.post_processing_solution - exact_final_T = lambda x: exact_T_solution(x, t=final_time) - exact_final_mobile = lambda x: exact_mobile_solution(x, t=final_time) - exact_final_trapped = lambda x: exact_trapped_solution(x, t=final_time) + def exact_final_T(x): + return exact_T_solution(x, t=final_time) + def exact_final_mobile(x): + return exact_mobile_solution(x, t=final_time) + def exact_final_trapped(x): + return exact_trapped_solution(x, t=final_time) L2_error_T = error_L2(T_computed, exact_final_T) L2_error_mobile = error_L2(mobile_computed, exact_final_mobile) @@ -216,7 +227,8 @@ def test_coupled_problem_non_matching_mesh(): test_mobile = F.Species("mobile", mobile=True) # define temperature sim - exact_T_solution = lambda x, t: 2 * x[0] ** 2 + 5 * t + def exact_T_solution(x, t): + return 2 * x[0] ** 2 + 5 * t dTdt = 5 @@ -249,11 +261,13 @@ def test_coupled_problem_non_matching_mesh(): ) # define hydrogen problem - exact_mobile_solution = lambda x, t: 4 * x[0] ** 2 + 10 * t + def exact_mobile_solution(x, t): + return 4 * x[0] ** 2 + 10 * t dmobiledt = 10 - D = lambda x, t: D_0 * ufl.exp(-E_D / (k_B * exact_T_solution(x, t))) + def D(x, t): + return D_0 * ufl.exp(-E_D / (k_B * exact_T_solution(x, t))) def mms_mobile_source(x, t): return dmobiledt - ufl.div(D(x, t) * ufl.grad(exact_mobile_solution(x, t))) @@ -302,7 +316,8 @@ def mms_mobile_source(x, t): test_coupled_problem.initialise() test_coupled_problem.run() - exact_solution = lambda x: exact_mobile_solution(x, t=final_time) + def exact_solution(x): + return exact_mobile_solution(x, t=final_time) computed_solution = test_mobile.post_processing_solution diff --git a/test/system_tests/test_cylindrical_coordinates.py b/test/system_tests/test_cylindrical_coordinates.py index eaa313e98..432549f46 100644 --- a/test/system_tests/test_cylindrical_coordinates.py +++ b/test/system_tests/test_cylindrical_coordinates.py @@ -13,7 +13,8 @@ def test_run_MMS_cylindrical(): my_mesh = F.Mesh1D(vertices=np.linspace(1, 2, 500), coordinate_system="cylindrical") - u_exact = lambda x: 1 + x[0] ** 2 + def u_exact(x): + return 1 + x[0] ** 2 f = -4 @@ -98,7 +99,8 @@ def c_exact_left(x): def c_exact_right(x): return K_S_right / K_S_left * c_exact_left(x) - lap_c = lambda r: 4 - 2 * r_interface / r + def lap_c(r): + return 4 - 2 * r_interface / r mat_1 = F.Material(D_0=D, E_D=0, K_S_0=K_S_left, E_K_S=0, solubility_law="sievert") mat_2 = F.Material(D_0=D, E_D=0, K_S_0=K_S_right, E_K_S=0, solubility_law="sievert") @@ -130,8 +132,10 @@ def c_exact_right(x): my_model.temperature = 500 - f_left = lambda x: -D * lap_c(x[0]) - f_right = lambda x: -D * K_S_right / K_S_left * lap_c(x[0]) + def f_left(x): + return -D * lap_c(x[0]) + def f_right(x): + return -D * K_S_right / K_S_left * lap_c(x[0]) my_model.sources = [ F.ParticleSource(value=f_left, volume=vol_1, species=H), diff --git a/test/system_tests/test_fluxes.py b/test/system_tests/test_fluxes.py index 1d3237a86..5a5db61c3 100644 --- a/test/system_tests/test_fluxes.py +++ b/test/system_tests/test_fluxes.py @@ -15,7 +15,8 @@ def test_flux_bc_1_mobile_MMS_steady_state(): MMS test with a flux BC considering one mobile species at steady state """ - u_exact = lambda x: 1 + 2 * x[0] ** 2 + def u_exact(x): + return 1 + 2 * x[0] ** 2 V = fem.functionspace(test_mesh_1d.mesh, ("Lagrange", 1)) T = fem.Function(V) @@ -64,7 +65,8 @@ def test_flux_bc_heat_transfer_steady_state(): MMS test with a flux BC in a heat transfer problem at steady state """ - u_exact = lambda x: 1 + 2 * x[0] ** 2 + def u_exact(x): + return 1 + 2 * x[0] ** 2 thermal_cond = 2.3 diff --git a/test/system_tests/test_multi_mat_penalty.py b/test/system_tests/test_multi_mat_penalty.py index 116da435e..37c5aed23 100644 --- a/test/system_tests/test_multi_mat_penalty.py +++ b/test/system_tests/test_multi_mat_penalty.py @@ -25,7 +25,6 @@ def half(x): ) # Split domain in half and set an interface tag of 5 - gdim = mesh.geometry.dim tdim = mesh.topology.dim fdim = tdim - 1 top_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, top_boundary) @@ -71,12 +70,14 @@ def test_2_materials_2d_mms(tmpdir): K_S_bot = 6.0 D_top = 2.0 D_bot = 5.0 - c_exact_top_ufl = lambda x: ( - 3 + ufl.sin(ufl.pi * (2 * x[1] + 0.5)) + 0.1 * ufl.cos(2 * ufl.pi * x[0]) - ) - c_exact_top_np = lambda x: ( - 3 + np.sin(np.pi * (2 * x[1] + 0.5)) + 0.1 * np.cos(2 * np.pi * x[0]) - ) + def c_exact_top_ufl(x): + return ( + 3 + ufl.sin(ufl.pi * (2 * x[1] + 0.5)) + 0.1 * ufl.cos(2 * ufl.pi * x[0]) + ) + def c_exact_top_np(x): + return ( + 3 + np.sin(np.pi * (2 * x[1] + 0.5)) + 0.1 * np.cos(2 * np.pi * x[0]) + ) def c_exact_bot_ufl(x): return K_S_bot / K_S_top**2 * c_exact_top_ufl(x) ** 2 diff --git a/test/system_tests/test_multi_material.py b/test/system_tests/test_multi_material.py index 3a7c1551a..d21d2688f 100644 --- a/test/system_tests/test_multi_material.py +++ b/test/system_tests/test_multi_material.py @@ -26,7 +26,6 @@ def half(x): ) # Split domain in half and set an interface tag of 5 - gdim = mesh.geometry.dim tdim = mesh.topology.dim fdim = tdim - 1 top_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, top_boundary) @@ -72,16 +71,18 @@ def test_2_materials_2d_mms(tmpdir): K_S_bot = 6.0 D_top = 2.0 D_bot = 5.0 - c_exact_top_ufl = lambda x: ( - 1 + ufl.sin(ufl.pi * (2 * x[0] + 0.5)) + ufl.cos(2 * ufl.pi * x[1]) - ) + def c_exact_top_ufl(x): + return ( + 1 + ufl.sin(ufl.pi * (2 * x[0] + 0.5)) + ufl.cos(2 * ufl.pi * x[1]) + ) def c_exact_bot_ufl(x): return K_S_bot / K_S_top * c_exact_top_ufl(x) - c_exact_top_np = lambda x: ( - 1 + np.sin(np.pi * (2 * x[0] + 0.5)) + np.cos(2 * np.pi * x[1]) - ) + def c_exact_top_np(x): + return ( + 1 + np.sin(np.pi * (2 * x[0] + 0.5)) + np.cos(2 * np.pi * x[1]) + ) def c_exact_bot_np(x): return K_S_bot / K_S_top * c_exact_top_np(x) @@ -127,12 +128,14 @@ def c_exact_bot_np(x): F.FixedConcentrationBC(bottom_surface, value=c_exact_bot_ufl, species=H), ] - source_top_val = lambda x: ( - 8 * ufl.pi**2 * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1])) - ) - source_bottom_val = lambda x: ( - 40 * ufl.pi**2 * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1])) - ) + def source_top_val(x): + return ( + 8 * ufl.pi**2 * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1])) + ) + def source_bottom_val(x): + return ( + 40 * ufl.pi**2 * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1])) + ) my_model.sources = [ F.ParticleSource(volume=top_domain, species=H, value=source_top_val), F.ParticleSource(volume=bottom_domain, species=H, value=source_bottom_val), diff --git a/test/system_tests/test_multi_material_change_of_var.py b/test/system_tests/test_multi_material_change_of_var.py index 8091a2663..17d64b6e2 100644 --- a/test/system_tests/test_multi_material_change_of_var.py +++ b/test/system_tests/test_multi_material_change_of_var.py @@ -24,7 +24,6 @@ def half(x): ) # Split domain in half and set an interface tag of 5 - gdim = mesh.geometry.dim tdim = mesh.topology.dim fdim = tdim - 1 top_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, top_boundary) @@ -156,16 +155,18 @@ def test_2_materials_2d_mms(tmpdir): K_S_bot = 6.0 D_top = 2.0 D_bot = 5.0 - c_exact_top_ufl = lambda x: ( - 1 + ufl.sin(ufl.pi * (2 * x[0] + 0.5)) + ufl.cos(2 * ufl.pi * x[1]) - ) + def c_exact_top_ufl(x): + return ( + 1 + ufl.sin(ufl.pi * (2 * x[0] + 0.5)) + ufl.cos(2 * ufl.pi * x[1]) + ) def c_exact_bot_ufl(x): return K_S_bot / K_S_top * c_exact_top_ufl(x) - c_exact_top_np = lambda x: ( - 1 + np.sin(np.pi * (2 * x[0] + 0.5)) + np.cos(2 * np.pi * x[1]) - ) + def c_exact_top_np(x): + return ( + 1 + np.sin(np.pi * (2 * x[0] + 0.5)) + np.cos(2 * np.pi * x[1]) + ) def c_exact_bot_np(x): return K_S_bot / K_S_top * c_exact_top_np(x) @@ -201,12 +202,14 @@ def c_exact_bot_np(x): F.FixedConcentrationBC(bottom_surface, value=c_exact_bot_ufl, species=H), ] - source_top_val = lambda x: ( - 8 * ufl.pi**2 * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1])) - ) - source_bottom_val = lambda x: ( - 40 * ufl.pi**2 * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1])) - ) + def source_top_val(x): + return ( + 8 * ufl.pi**2 * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1])) + ) + def source_bottom_val(x): + return ( + 40 * ufl.pi**2 * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1])) + ) my_model.sources = [ F.ParticleSource(volume=top_domain, species=H, value=source_top_val), F.ParticleSource(volume=bottom_domain, species=H, value=source_bottom_val), diff --git a/test/system_tests/test_non_homogeneous_density.py b/test/system_tests/test_non_homogeneous_density.py index c021641f8..73d55a8d1 100644 --- a/test/system_tests/test_non_homogeneous_density.py +++ b/test/system_tests/test_non_homogeneous_density.py @@ -6,12 +6,17 @@ import festim as F -step_function_space = lambda x: ufl.conditional(ufl.gt(x[0], 0.5), 10, 0.0) -step_function_time_non_homogeneous = lambda x, t: ufl.conditional( + +def step_function_space(x): + return ufl.conditional(ufl.gt(x[0], 0.5), 10, 0.0) +def step_function_time_non_homogeneous(x, t): + return ufl.conditional( ufl.gt(t, 50), 10 - 10 * x[0], 0.0 ) -step_function_time_homogeneous = lambda t: 10.0 if t > 50 else 2.0 -simple_space = lambda x: 10 - 10 * x[0] +def step_function_time_homogeneous(t): + return 10.0 if t > 50 else 2.0 +def simple_space(x): + return 10 - 10 * x[0] @pytest.mark.parametrize( diff --git a/test/system_tests/test_spherical_coordinates.py b/test/system_tests/test_spherical_coordinates.py index d477f484d..956fe930d 100644 --- a/test/system_tests/test_spherical_coordinates.py +++ b/test/system_tests/test_spherical_coordinates.py @@ -13,9 +13,10 @@ def test_run_MMS_spherical(): """ my_mesh = F.Mesh1D(vertices=np.linspace(1, 2, 1000), coordinate_system="spherical") - V = fem.functionspace(my_mesh.mesh, ("Lagrange", 1)) + fem.functionspace(my_mesh.mesh, ("Lagrange", 1)) - u_exact = lambda x: 3 + 2 * x[0] ** 2 + def u_exact(x): + return 3 + 2 * x[0] ** 2 f = -12 @@ -99,7 +100,8 @@ def c_exact_left(x): def c_exact_right(x): return K_S_right / K_S_left * c_exact_left(x) - lap_c = lambda r: -(4 * r_interface / r) + 6 + def lap_c(r): + return -(4 * r_interface / r) + 6 mat_1 = F.Material(D_0=D, E_D=0, K_S_0=K_S_left, E_K_S=0, solubility_law="sievert") mat_2 = F.Material(D_0=D, E_D=0, K_S_0=K_S_right, E_K_S=0, solubility_law="sievert") @@ -131,8 +133,10 @@ def c_exact_right(x): my_model.temperature = 500 - f_left = lambda x: -D * lap_c(x[0]) - f_right = lambda x: -D * K_S_right / K_S_left * lap_c(x[0]) + def f_left(x): + return -D * lap_c(x[0]) + def f_right(x): + return -D * K_S_right / K_S_left * lap_c(x[0]) my_model.sources = [ F.ParticleSource(value=f_left, volume=vol_1, species=H), diff --git a/test/test_average_surface.py b/test/test_average_surface.py index 5ab11170e..b8c4950ea 100644 --- a/test/test_average_surface.py +++ b/test/test_average_surface.py @@ -16,7 +16,7 @@ def test_average_surface_compute_1D(): dummy_volume = F.VolumeSubdomain1D( id=1, borders=[0, L], material=F.Material(D_0=1, E_D=1, name="dummy") ) - facet_meshtags, temp = my_mesh.define_meshtags( + facet_meshtags, _temp = my_mesh.define_meshtags( surface_subdomains=[dummy_surface], volume_subdomains=[dummy_volume] ) ds = ufl.Measure("ds", domain=my_mesh.mesh, subdomain_data=facet_meshtags) diff --git a/test/test_average_volume.py b/test/test_average_volume.py index 7bb93d42e..8f71121b9 100644 --- a/test/test_average_volume.py +++ b/test/test_average_volume.py @@ -14,7 +14,7 @@ def test_average_volume_compute_1D(): L = 6.0 my_mesh = F.Mesh1D(np.linspace(0, L, 10000)) dummy_volume = F.VolumeSubdomain1D(id=1, borders=[0, L], material=dummy_mat) - temp, cell_meshtags = my_mesh.define_meshtags( + _temp, cell_meshtags = my_mesh.define_meshtags( surface_subdomains=[], volume_subdomains=[dummy_volume] ) dx = ufl.Measure("dx", domain=my_mesh.mesh, subdomain_data=cell_meshtags) diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index d33e002a2..01ffa3251 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -250,7 +250,7 @@ def test_define_D_global_different_temperatures(): is different in the volume subdomains""" D_0, E_D = 1.5, 0.1 my_mat = F.Material(D_0=D_0, E_D=E_D, name="my_mat") - surf = F.SurfaceSubdomain1D(id=1, x=0) + F.SurfaceSubdomain1D(id=1, x=0) H = F.Species("H") my_model = F.HydrogenTransportProblem( @@ -268,7 +268,7 @@ def test_define_D_global_different_temperatures(): my_model.t = 1 my_model.define_temperature() - D_computed, D_expr = my_model.define_D_global(H) + D_computed, _D_expr = my_model.define_D_global(H) computed_values = [D_computed.x.array[0], D_computed.x.array[-1]] @@ -304,7 +304,7 @@ def test_define_D_global_different_materials(): my_model.t = 0 my_model.define_temperature() - D_computed, D_expr = my_model.define_D_global(H) + D_computed, _D_expr = my_model.define_D_global(H) computed_values = [D_computed.x.array[0], D_computed.x.array[-1]] @@ -399,7 +399,7 @@ def test_define_D_global_multispecies(): my_mat = F.Material( D_0={A: D_0_A, B: D_0_B}, E_D={A: E_D_A, B: E_D_B}, name="my_mat" ) - surf = F.SurfaceSubdomain1D(id=1, x=1) + F.SurfaceSubdomain1D(id=1, x=1) my_model = F.HydrogenTransportProblem( mesh=F.Mesh1D(np.linspace(0, 1, num=101)), @@ -415,8 +415,8 @@ def test_define_D_global_multispecies(): my_model.t = 0 my_model.define_temperature() - D_A_computed, D_A_expr = my_model.define_D_global(A) - D_B_computed, D_B_expr = my_model.define_D_global(B) + D_A_computed, _D_A_expr = my_model.define_D_global(A) + D_B_computed, _D_B_expr = my_model.define_D_global(B) computed_values = [D_A_computed.x.array[-1], D_B_computed.x.array[-1]] diff --git a/test/test_helpers.py b/test/test_helpers.py index 971069e3f..38c583fe7 100644 --- a/test/test_helpers.py +++ b/test/test_helpers.py @@ -199,7 +199,8 @@ def test_input_values_of_constants_and_functions_are_accepted(value): def test_input_values_of_expressions_are_accepted(): """Test that the input values of constants and functions are accepted""" - my_func = lambda x: 1.0 + x[0] + def my_func(x): + return 1.0 + x[0] kwargs = {} kwargs["x"] = x mapped_func = my_func(**kwargs) diff --git a/test/test_initial_condition.py b/test/test_initial_condition.py index 588785445..50a3a545f 100644 --- a/test/test_initial_condition.py +++ b/test/test_initial_condition.py @@ -75,7 +75,8 @@ def test_warning_raised_when_giving_time_as_arg(): my_species = F.Species("test") my_species.prev_solution = fem.Function(V) - my_value = lambda t: 1.0 + t + def my_value(t): + return 1.0 + t init_cond = F.InitialConcentration(value=my_value, species=my_species, volume=vol) @@ -97,7 +98,8 @@ def test_warning_raised_when_giving_time_as_arg_initial_temperature(): vol = F.VolumeSubdomain(id=1, material=dummy_mat) - my_value = lambda t: 1.0 + t + def my_value(t): + return 1.0 + t init_cond = F.InitialTemperature(value=my_value, volume=vol) @@ -304,7 +306,8 @@ def test_create_initial_temperature_from_function(): vol = F.VolumeSubdomain(id=1, material=dummy_mat) # create an initial temperature from a function - T = lambda x: 300 + 10 * x[0] + def T(x): + return 300 + 10 * x[0] V = fem.functionspace(test_mesh.mesh, ("Lagrange", 1)) T_func = fem.Function(V) T_func.interpolate(T) @@ -368,10 +371,10 @@ def test_initial_condition_discontinuous(): my_model.initialise() - spe1_left, spe1_to_vol1 = vol1.u_n.function_space.sub(0).collapse() - spe2_left, spe2_to_vol1 = vol1.u_n.function_space.sub(1).collapse() - spe1_right, spe1_to_vol2 = vol2.u_n.function_space.sub(0).collapse() - spe2_right, spe2_to_vol2 = vol2.u_n.function_space.sub(1).collapse() + _spe1_left, spe1_to_vol1 = vol1.u_n.function_space.sub(0).collapse() + _spe2_left, spe2_to_vol1 = vol1.u_n.function_space.sub(1).collapse() + _spe1_right, spe1_to_vol2 = vol2.u_n.function_space.sub(0).collapse() + _spe2_right, spe2_to_vol2 = vol2.u_n.function_space.sub(1).collapse() prev_solution_spe1_left = vol1.u_n.x.array[spe1_to_vol1] prev_solution_spe2_left = vol1.u_n.x.array[spe2_to_vol1] @@ -523,8 +526,8 @@ def test_initial_condition_mixed_domain_multispecies(): my_model.initialise() - spe1_left, spe1_to_vol1 = vol1.u_n.function_space.sub(0).collapse() - spe2_left, spe2_to_vol1 = vol1.u_n.function_space.sub(1).collapse() + _spe1_left, spe1_to_vol1 = vol1.u_n.function_space.sub(0).collapse() + _spe2_left, spe2_to_vol1 = vol1.u_n.function_space.sub(1).collapse() prev_solution_spe1_left = vol1.u_n.x.array[spe1_to_vol1] prev_solution_spe2_left = vol1.u_n.x.array[spe2_to_vol1] diff --git a/test/test_maximum_surface.py b/test/test_maximum_surface.py index 54eb62de2..1636c8f66 100644 --- a/test/test_maximum_surface.py +++ b/test/test_maximum_surface.py @@ -14,7 +14,7 @@ def test_maximum_surface_compute_1D(): dummy_surface = F.SurfaceSubdomain1D(id=1, x=4) # create mesh tags - ft, ct = my_mesh.define_meshtags( + ft, _ct = my_mesh.define_meshtags( surface_subdomains=[dummy_surface], volume_subdomains=[ F.VolumeSubdomain1D(id=1, material=F.Material(D_0=1, E_D=0), borders=[0, L]) diff --git a/test/test_maximum_volume.py b/test/test_maximum_volume.py index ea3719d82..393e5b933 100644 --- a/test/test_maximum_volume.py +++ b/test/test_maximum_volume.py @@ -14,7 +14,7 @@ def test_maximum_volume_compute_1D(): dummy_volume = F.VolumeSubdomain1D(id=1, borders=[0, L], material=dummy_material) # create mesh tags - ft, ct = my_mesh.define_meshtags( + _ft, ct = my_mesh.define_meshtags( surface_subdomains=[F.SurfaceSubdomain1D(id=1, x=0)], volume_subdomains=[dummy_volume], interfaces=None, diff --git a/test/test_minimum_surface.py b/test/test_minimum_surface.py index b6a64e9aa..e7ce02808 100644 --- a/test/test_minimum_surface.py +++ b/test/test_minimum_surface.py @@ -14,7 +14,7 @@ def test_minimum_surface_export_compute_1D(): dummy_surface = F.SurfaceSubdomain1D(id=1, x=4) # create mesh tags - ft, ct = my_mesh.define_meshtags( + ft, _ct = my_mesh.define_meshtags( surface_subdomains=[dummy_surface], volume_subdomains=[ F.VolumeSubdomain1D(id=1, material=F.Material(D_0=1, E_D=0), borders=[0, L]) diff --git a/test/test_minimum_volume.py b/test/test_minimum_volume.py index 820cab40f..f9e5a985d 100644 --- a/test/test_minimum_volume.py +++ b/test/test_minimum_volume.py @@ -14,7 +14,7 @@ def test_minimum_volume_compute_1D(): dummy_volume = F.VolumeSubdomain1D(id=1, borders=[0, L], material=dummy_material) # create mesh tags - ft, ct = my_mesh.define_meshtags( + _ft, ct = my_mesh.define_meshtags( surface_subdomains=[F.SurfaceSubdomain1D(id=1, x=0)], volume_subdomains=[dummy_volume], interfaces=None, diff --git a/test/test_species.py b/test/test_species.py index e5d2af306..44f62dcf3 100644 --- a/test/test_species.py +++ b/test/test_species.py @@ -152,7 +152,8 @@ def test_implicit_species_wrong_type(): """ mesh = create_unit_cube(MPI.COMM_WORLD, 10, 10, 10) - density = lambda t: "coucou" + def density(t): + return "coucou" species = F.ImplicitSpecies(n=density) with pytest.raises( diff --git a/test/test_total_surface.py b/test/test_total_surface.py index cce43aa7a..12b97165e 100644 --- a/test/test_total_surface.py +++ b/test/test_total_surface.py @@ -16,7 +16,7 @@ def test_total_surface_compute_1D(): dummy_volume = F.VolumeSubdomain1D( id=1, borders=[0, L], material=F.Material(D_0=1, E_D=1, name="dummy") ) - facet_meshtags, temp = my_mesh.define_meshtags( + facet_meshtags, _temp = my_mesh.define_meshtags( surface_subdomains=[dummy_surface], volume_subdomains=[dummy_volume] ) ds = ufl.Measure("ds", domain=my_mesh.mesh, subdomain_data=facet_meshtags) diff --git a/test/test_total_volume.py b/test/test_total_volume.py index 2b2b43e39..e0e3716ce 100644 --- a/test/test_total_volume.py +++ b/test/test_total_volume.py @@ -17,7 +17,7 @@ def test_total_volume_export_compute(): L = 4.0 my_mesh = F.Mesh1D(np.linspace(0, L, 10000)) dummy_volume = F.VolumeSubdomain1D(id=1, borders=[0, L], material=dummy_mat) - temp, cell_meshtags = my_mesh.define_meshtags( + _temp, cell_meshtags = my_mesh.define_meshtags( surface_subdomains=[], volume_subdomains=[dummy_volume] ) dx = ufl.Measure("dx", domain=my_mesh.mesh, subdomain_data=cell_meshtags) From cb2026bb30918833f4923bca1809559327698d03 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 22 Apr 2026 14:44:58 -0400 Subject: [PATCH 08/14] fixed buidl --- test/test_reaction.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_reaction.py b/test/test_reaction.py index 0dd24e331..397d82faa 100644 --- a/test/test_reaction.py +++ b/test/test_reaction.py @@ -401,7 +401,7 @@ def test_product_setter_raise_error_E_p_no_product(): vol1 = F.VolumeSubdomain1D(id=1, borders=[0, 0.5], material=mat) vol2 = F.VolumeSubdomain1D(id=2, borders=[0.5, 1], material=mat) my_model = F.HydrogenTransportProblemDiscontinuous() -my_model.mesh = F.Mesh1D(np.linspace(0, 1, 100)) +my_model.mesh = F.Mesh1D(np.linspace(0, 1, 101)) my_model.subdomains = [vol1, vol2] From da21a4d11a1cac774b3239515add456123a9de74 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 22 Apr 2026 14:59:19 -0400 Subject: [PATCH 09/14] bunch of ruff fixes --- docs/source/conf.py | 4 +- .../advection_diffusion_multi_material.py | 6 +- examples/convergence_rates.py | 7 +- examples/multi_material_1d.py | 3 +- examples/multi_material_2d.py | 4 + examples/multi_material_2d_bis.py | 2 + examples/multi_material_transient.py | 4 +- examples/surface_reaction_example.py | 2 +- src/festim/__init__.py | 1 + src/festim/advection.py | 13 +- .../boundary_conditions/dirichlet_bc.py | 52 +++---- src/festim/boundary_conditions/flux_bc.py | 74 +++++---- src/festim/boundary_conditions/henrys_bc.py | 28 ++-- src/festim/boundary_conditions/sieverts_bc.py | 28 ++-- .../boundary_conditions/surface_reaction.py | 8 +- src/festim/coupled_heat_hydrogen_problem.py | 5 +- src/festim/exports/average_surface.py | 7 +- src/festim/exports/average_volume.py | 11 +- src/festim/exports/maximum_surface.py | 11 +- src/festim/exports/maximum_volume.py | 12 +- src/festim/exports/minimum_surface.py | 11 +- src/festim/exports/minimum_volume.py | 11 +- src/festim/exports/surface_flux.py | 4 +- src/festim/exports/surface_quantity.py | 6 +- src/festim/exports/total_surface.py | 7 +- src/festim/exports/total_volume.py | 7 +- src/festim/exports/volume_quantity.py | 6 +- src/festim/exports/vtx.py | 15 +- src/festim/exports/xdmf.py | 10 +- src/festim/heat_transfer_problem.py | 28 ++-- src/festim/helpers.py | 29 ++-- src/festim/hydrogen_transport_problem.py | 117 +++++++-------- src/festim/initial_condition.py | 12 +- src/festim/material.py | 15 +- src/festim/mesh/mesh.py | 15 +- src/festim/mesh/mesh_1d.py | 7 +- src/festim/mesh/mesh_from_xdmf.py | 10 +- src/festim/problem.py | 26 ++-- src/festim/reaction.py | 17 ++- src/festim/settings.py | 1 - src/festim/source.py | 9 +- src/festim/species.py | 20 +-- src/festim/stepsize.py | 14 +- src/festim/subdomain/interface.py | 11 +- src/festim/subdomain/surface_subdomain.py | 20 +-- src/festim/subdomain/volume_subdomain.py | 13 +- src/festim/trap.py | 13 +- test/benchmark.py | 4 +- test/system_tests/test_1_mobile_1_trap.py | 28 ++-- test/system_tests/test_1_mobile_species.py | 28 ++-- .../test_advection_term_integration.py | 5 +- .../test_coupled_heat_hydrogen.py | 10 +- .../test_cylindrical_coordinates.py | 13 +- test/system_tests/test_fluxes.py | 8 +- test/system_tests/test_heat_transfer.py | 4 +- test/system_tests/test_io.py | 12 +- test/system_tests/test_misc.py | 23 ++- test/system_tests/test_multi_mat_penalty.py | 10 +- test/system_tests/test_multi_material.py | 24 ++- .../test_multi_material_change_of_var.py | 15 +- .../test_non_homogeneous_density.py | 10 +- .../test_reaction_not_in_every_volume.py | 4 +- .../test_spherical_coordinates.py | 13 +- test/system_tests/test_surface_reactions.py | 27 ++-- test/test_advection.py | 2 +- test/test_average_surface.py | 2 +- test/test_average_volume.py | 2 +- test/test_complex_reaction.py | 20 +-- test/test_coupled_heat_hydrogen_problem.py | 20 +-- test/test_dirichlet_bc.py | 39 +++-- test/test_flux_bc.py | 34 +++-- test/test_h_transport_problem.py | 141 +++++++++--------- test/test_heat_transfer_problem.py | 44 +++--- test/test_helpers.py | 41 ++--- test/test_henrysbc.py | 5 +- test/test_initial_condition.py | 49 +++--- test/test_material.py | 56 ++++--- test/test_maximum_surface.py | 2 +- test/test_maximum_volume.py | 2 +- test/test_mesh.py | 19 ++- test/test_minimum_surface.py | 2 +- test/test_minimum_volume.py | 2 +- test/test_multispecies_problem.py | 6 +- test/test_permeation_problem.py | 10 +- test/test_reaction.py | 45 +++--- test/test_settings.py | 6 +- test/test_sievertsbc.py | 5 +- test/test_source.py | 34 +++-- test/test_species.py | 35 ++--- test/test_stepsize.py | 38 ++--- test/test_subdomains.py | 37 ++--- test/test_surface_quantity.py | 20 +-- test/test_total_surface.py | 2 +- test/test_total_volume.py | 4 +- test/test_volume_quantity.py | 7 +- test/test_vtx.py | 22 +-- test/test_xdmf.py | 16 +- 97 files changed, 876 insertions(+), 867 deletions(-) diff --git a/docs/source/conf.py b/docs/source/conf.py index 5765a588f..3e5c5c7d6 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -12,6 +12,8 @@ import os import sys +import map + sys.path.insert(0, os.path.abspath("../../src")) @@ -20,7 +22,6 @@ # Add the directory containing your Python script to the Python path sys.path.insert(0, os.path.abspath(".")) -import map m = map.generate_map() current_dir = os.path.dirname(__file__) @@ -37,7 +38,6 @@ # Add the directory containing your Python script to the Python path sys.path.insert(0, os.path.abspath(".")) -import map m = map.generate_map() current_dir = os.path.dirname(__file__) diff --git a/examples/advection_diffusion_multi_material.py b/examples/advection_diffusion_multi_material.py index c56ad8389..96850d2b8 100644 --- a/examples/advection_diffusion_multi_material.py +++ b/examples/advection_diffusion_multi_material.py @@ -80,9 +80,9 @@ def velocity_func(x): return values mesh2, _mt2, ct2 = generate_mesh(n=10) - submesh, _cell_map, _v_map = dolfinx.mesh.create_submesh(mesh2, ct2.dim, ct2.find(4))[ - 0:3 - ] + submesh, _cell_map, _v_map = dolfinx.mesh.create_submesh( + mesh2, ct2.dim, ct2.find(4) + )[0:3] v_cg = basix.ufl.element( "Lagrange", submesh.topology.cell_name(), 2, shape=(submesh.geometry.dim,) ) diff --git a/examples/convergence_rates.py b/examples/convergence_rates.py index 6286c8e5e..a6b24d702 100644 --- a/examples/convergence_rates.py +++ b/examples/convergence_rates.py @@ -112,7 +112,9 @@ def mms_source(x, t): my_problem.settings = F.Settings( atol=1e-8, rtol=1e-10, - final_time=1, # final time shouldn't be too long so that a potential error at the initial timestep is not negligible + # final time shouldn't be too long so that a potential error at the + # initial timestep is not negligible + final_time=1, ) # Forward euler isn't great so dt should be small @@ -127,7 +129,8 @@ def mms_source(x, t): my_problem.run() computed_solution = my_problem.u - # we use the exact final time of the simulation which may differ from the one specified in the settings + # we use the exact final time of the simulation which may differ from + # the one specified in the settings final_time_sim = my_problem.t.value def exact_solution_end(x): diff --git a/examples/multi_material_1d.py b/examples/multi_material_1d.py index 2fb2a815c..4810f27bb 100644 --- a/examples/multi_material_1d.py +++ b/examples/multi_material_1d.py @@ -43,7 +43,8 @@ left_surface = F.SurfaceSubdomain1D(id=1, x=vertices[0]) right_surface = F.SurfaceSubdomain1D(id=2, x=vertices[-1]) -# the ids here are arbitrary in 1D, you can put anything as long as it's not the same as the surfaces +# the ids here are arbitrary in 1D, you can put anything as long as it's +# not the same as the surfaces my_model.interfaces = [ F.Interface(6, (left_domain, middle_domain)), F.Interface(7, (middle_domain, right_domain)), diff --git a/examples/multi_material_2d.py b/examples/multi_material_2d.py index 314d9ef32..201961511 100644 --- a/examples/multi_material_2d.py +++ b/examples/multi_material_2d.py @@ -124,13 +124,17 @@ def half(x): from dolfinx.mesh import EntityMap # noqa: F401 legacy_entity_map = False + def entity_map_wrapper(e_map): return list(e_map.values()) + entity_maps = [sd.cell_map for sd in my_model.volume_subdomains] except ImportError: legacy_entity_map = True + def entity_map_wrapper(e_map): return e_map + entity_maps = { sd.submesh: sd.parent_to_submesh for sd in my_model.volume_subdomains } diff --git a/examples/multi_material_2d_bis.py b/examples/multi_material_2d_bis.py index 6ea69169b..933747c13 100644 --- a/examples/multi_material_2d_bis.py +++ b/examples/multi_material_2d_bis.py @@ -167,10 +167,12 @@ def half(x): from dolfinx.mesh import EntityMap # noqa: F401 legacy_entity_map = False + def entity_map_wrapper(e_map): return list(e_map.values()) except ImportError: legacy_entity_map = True + def entity_map_wrapper(e_map): return e_map diff --git a/examples/multi_material_transient.py b/examples/multi_material_transient.py index 8495336fa..54b30cfc0 100644 --- a/examples/multi_material_transient.py +++ b/examples/multi_material_transient.py @@ -36,8 +36,8 @@ left_surface = F.SurfaceSubdomain1D(id=1, x=vertices[0]) right_surface = F.SurfaceSubdomain1D(id=2, x=vertices[-1]) -# the ids here are arbitrary in 1D, you can put anything as long as it's not the same as the surfaces -# TODO remove mesh and meshtags from these arguments +# the ids here are arbitrary in 1D, you can put anything as long as it's not the same +# as the surfaces my_model.interfaces = [ F.Interface(6, (left_domain, middle_domain)), F.Interface(7, (middle_domain, right_domain)), diff --git a/examples/surface_reaction_example.py b/examples/surface_reaction_example.py index 2f022e0ba..14d41ea80 100644 --- a/examples/surface_reaction_example.py +++ b/examples/surface_reaction_example.py @@ -1,4 +1,5 @@ import dolfinx.fem as fem +import matplotlib.pyplot as plt import numpy as np import ufl @@ -98,7 +99,6 @@ def compute(self, ds): my_model.initialise() my_model.run() -import matplotlib.pyplot as plt plt.stackplot( H_flux_left.t, diff --git a/src/festim/__init__.py b/src/festim/__init__.py index 823d519d9..d0a824587 100644 --- a/src/festim/__init__.py +++ b/src/festim/__init__.py @@ -1,3 +1,4 @@ +# ruff: noqa: F401, E402 from importlib import metadata try: diff --git a/src/festim/advection.py b/src/festim/advection.py index 2ecf944f3..10081f823 100644 --- a/src/festim/advection.py +++ b/src/festim/advection.py @@ -10,8 +10,7 @@ class AdvectionTerm: - """ - Advection term class + """Advection term class. args: velocity: the velocity field or function @@ -22,7 +21,6 @@ class AdvectionTerm: velocity: the velocity field or function subdomain: the volume subdomain where the velocity is to be applied species: the species to which the velocity field is acting on - """ velocity: fem.Function @@ -94,9 +92,8 @@ def species(self, value): class VelocityField(Value): - """ - A class to handle input values of velocity fields from users and convert them to a - relevent fenics object + """A class to handle input values of velocity fields from users and convert them to + a relevent fenics object. Args: input_value: The value of the user input @@ -124,7 +121,7 @@ def convert_input_value( function_space: fem.FunctionSpace, t: fem.Constant | None = None, ): - """Converts a user given value to a relevent fenics object + """Converts a user given value to a relevent fenics object. Args: function_space: the function space of the fenics object @@ -165,7 +162,7 @@ def convert_input_value( nmm_interpolate(self.fenics_object, vel) def update(self, t: fem.Constant): - """Updates the velocity field + """Updates the velocity field. Args: t: the time diff --git a/src/festim/boundary_conditions/dirichlet_bc.py b/src/festim/boundary_conditions/dirichlet_bc.py index f947ab98c..4ae5604cd 100644 --- a/src/festim/boundary_conditions/dirichlet_bc.py +++ b/src/festim/boundary_conditions/dirichlet_bc.py @@ -14,9 +14,7 @@ class DirichletBCBase: - """ - Dirichlet boundary condition class - u = value + """Dirichlet boundary condition class u = value. Args: subdomain: The surface subdomain where the boundary condition is applied @@ -28,7 +26,6 @@ class DirichletBCBase: value_fenics: The value of the boundary condition in fenics format bc_expr: The expression of the boundary condition that is used to update the `value_fenics` - """ subdomain: _subdomain.SurfaceSubdomain @@ -67,14 +64,14 @@ def value_fenics(self, value: None | fem.Function | fem.Constant | np.ndarray): if not isinstance(value, (fem.Function, fem.Constant, np.ndarray)): # FIXME: Should we allow sending in a callable here? raise TypeError( - "Value must be a dolfinx.fem.Function, dolfinx.fem.Constant, or a np.ndarray not" + "Value must be a dolfinx.fem.Function, dolfinx.fem.Constant, or a np.ndarray not" # noqa: E501 + f"{type(value)}" ) self._value_fenics = value @property def time_dependent(self) -> bool: - """Returns true if the value of the boundary condition is time dependent""" + """Returns true if the value of the boundary condition is time dependent.""" if self.value is None: return False if isinstance(self.value, fem.Constant): @@ -92,17 +89,20 @@ def define_surface_subdomain_dofs( ) -> npt.NDArray[np.int32] | tuple[npt.NDArray[np.int32], npt.NDArray[np.int32]]: """Defines the facets and the degrees of freedom of the boundary condition. - Given the input meshtags, find all facets matching the boundary condition subdomain ID, + Given the input meshtags, find all facets matching the boundary + condition subdomain ID, and locate all DOFs associated with the input function space(s). Note: - For sub-spaces, a tuple of sub-spaces are expected as input, and a tuple of arrays + For sub-spaces, a tuple of sub-spaces are expected as input, and a + tuple of arrays associated to each of the function spaces are returned. Args: facet_meshtags: MeshTags describing some facets in the domain mesh: - function_space: The function space or a tuple of function spaces: (sub, collapsed) + function_space: The function space or a tuple of function spaces: + (sub, collapsed) """ mesh = ( function_space[0].mesh @@ -111,11 +111,11 @@ def define_surface_subdomain_dofs( ) if facet_meshtags.topology != mesh.topology._cpp_object: raise ValueError( - "Mesh of function-space is not the same as the one used for the meshtags" + "Mesh of function-space is not the same as the one used for the meshtags" # noqa: E501 ) if mesh.topology.dim - 1 != facet_meshtags.dim: raise ValueError( - f"Meshtags of dimension {facet_meshtags.dim}, expected {mesh.topology.dim - 1}" + f"Meshtags of dimension {facet_meshtags.dim}, expected {mesh.topology.dim - 1}" # noqa: E501 ) bc_dofs = fem.locate_dofs_topological( function_space, facet_meshtags.dim, facet_meshtags.find(self.subdomain.id) @@ -124,7 +124,7 @@ def define_surface_subdomain_dofs( return bc_dofs def update(self, t: float): - """Updates the boundary condition value + """Updates the boundary condition value. Args: t: the time @@ -144,11 +144,13 @@ class FixedConcentrationBC(DirichletBCBase): Args: subdomain (festim.Subdomain): the surface subdomain where the boundary condition is applied - value: The value of the boundary condition. It can be a function of space and/or time + value: The value of the boundary condition. It can be a function of + space and/or time species: The name of the species Attributes: - temperature_dependent (bool): True if the value of the bc is temperature dependent + temperature_dependent (bool): True if the value of the bc is + temperature dependent Examples: @@ -202,17 +204,18 @@ def create_value( K_S: fem.Function = None, ): """Creates the value of the boundary condition as a fenics object and sets it to - self.value_fenics. - If the value is a constant, it is converted to a `dolfinx.fem.Constant`. - If the value is a function of t, it is converted to `dolfinx.fem.Constant`. - Otherwise, it is converted to a `dolfinx.fem.Function`.Function and the - expression of the function is stored in `bc_expr`. + self.value_fenics. If the value is a constant, it is converted to a + `dolfinx.fem.Constant`. If the value is a function of t, it is converted to + `dolfinx.fem.Constant`. Otherwise, it is converted to a + `dolfinx.fem.Function`.Function and the expression of the function is stored in + `bc_expr`. Args: function_space: the function space temperature: The temperature t: the time - K_S: The solubility of the species. If provided, the value of the boundary condition + K_S: The solubility of the species. If provided, the value of the + boundary condition is divided by K_S (change of variable method). """ mesh = function_space.mesh @@ -313,11 +316,10 @@ def create_value( class FixedTemperatureBC(DirichletBCBase): def create_value(self, function_space: fem.FunctionSpace, t: fem.Constant): """Creates the value of the boundary condition as a fenics object and sets it to - self.value_fenics. - If the value is a constant, it is converted to a `dolfinx.fem.Constant`. - If the value is a function of t, it is converted to a `dolfinx.fem.Constant`. - Otherwise, it is converted to a` dolfinx.fem.Function` and the - expression of the function is stored in `bc_expr`. + self.value_fenics. If the value is a constant, it is converted to a + `dolfinx.fem.Constant`. If the value is a function of t, it is converted to a + `dolfinx.fem.Constant`. Otherwise, it is converted to a` dolfinx.fem.Function` + and the expression of the function is stored in `bc_expr`. Args: function_space: the function space diff --git a/src/festim/boundary_conditions/flux_bc.py b/src/festim/boundary_conditions/flux_bc.py index 5ae48c777..09eebe6a8 100644 --- a/src/festim/boundary_conditions/flux_bc.py +++ b/src/festim/boundary_conditions/flux_bc.py @@ -7,8 +7,7 @@ class FluxBCBase: - """ - Flux boundary condition class + """Flux boundary condition class. Ensuring the gradient of the solution u at a boundary: @@ -30,11 +29,12 @@ class FluxBCBase: subdomain (festim.SurfaceSubdomain): the surface subdomain where the boundary condition is applied value (float, fem.Constant, callable): the value of the boundary condition - value_fenics (fem.Function or fem.Constant): the value of the boundary condition in + value_fenics (fem.Function or fem.Constant): the value of the boundary + condition in fenics format - bc_expr (fem.Expression): the expression of the boundary condition that is used to + bc_expr (fem.Expression): the expression of the boundary condition that + is used to update the value_fenics - """ def __init__(self, subdomain: SurfaceSubdomain, value): @@ -57,7 +57,8 @@ def value_fenics(self, value): value, fem.Function | fem.Constant | np.ndarray | ufl.core.expr.Expr ): raise TypeError( - f"Value must be a dolfinx.fem.Function, dolfinx.fem.Constant, np.ndarray or ufl.core.expr.Expr not {type(value)}" # noqa: E501 + f"Value must be a dolfinx.fem.Function, dolfinx.fem.Constant," + f" np.ndarray or ufl.core.expr.Expr not {type(value)}" ) self._value_fenics = value @@ -87,10 +88,9 @@ def temperature_dependent(self): def create_value_fenics(self, mesh, temperature, t: fem.Constant): """Creates the value of the boundary condition as a fenics object and sets it to - self.value_fenics. - If the value is a constant, it is converted to a fenics.Constant. - If the value is a function of t, it is converted to a fenics.Constant. - Otherwise, it is converted to a ufl Expression + self.value_fenics. If the value is a constant, it is converted to a + fenics.Constant. If the value is a function of t, it is converted to a + fenics.Constant. Otherwise, it is converted to a ufl Expression. Args: mesh (dolfinx.mesh.Mesh) : the mesh @@ -109,7 +109,7 @@ def create_value_fenics(self, mesh, temperature, t: fem.Constant): # only t is an argument if not isinstance(self.value(t=float(t)), (float, int)): raise ValueError( - f"self.value should return a float or an int, not {type(self.value(t=float(t)))} " + f"self.value should return a float or an int, not {type(self.value(t=float(t)))} " # noqa: E501 ) self.value_fenics = F.as_fenics_constant( mesh=mesh, value=self.value(t=float(t)) @@ -126,7 +126,7 @@ def create_value_fenics(self, mesh, temperature, t: fem.Constant): self.value_fenics = self.value(**kwargs) def update(self, t): - """Updates the flux bc value + """Updates the flux bc value. Args: t (float): the time @@ -138,22 +138,27 @@ def update(self, t): class ParticleFluxBC(FluxBCBase): - """ - Particle flux boundary condition class + """Particle flux boundary condition class. + Ensuring the gradient of the solution c at a boundary: -D * grad(c) * n = f - where D is the material diffusivity, n is the outwards normal vector of the boundary, f is a function of space and time. + where D is the material diffusivity, n is the outwards normal vector of the + boundary, f is a function of space and time. Args: - subdomain (festim.SurfaceSubdomain): the surface subdomain where the particle flux + subdomain (festim.SurfaceSubdomain): the surface subdomain where the + particle flux is applied value (float, fem.Constant, callable): the value of the particle flux species (festim.Species): the species to which the flux is applied - species_dependent_value (dict): a dictionary containing the species that the value. Example: {"name": species} - where "name" is the variable name in the callable value and species is a festim.Species object. + species_dependent_value (dict): a dictionary containing the species + that the value. Example: {"name": species} + where "name" is the variable name in the callable value and species + is a festim.Species object. Attributes: - subdomain (festim.SurfaceSubdomain): the surface subdomain where the particle flux + subdomain (festim.SurfaceSubdomain): the surface subdomain where the + particle flux is applied value (float or fem.Constant): the value of the particle flux species (festim.Species): the species to which the flux is applied @@ -161,8 +166,10 @@ class ParticleFluxBC(FluxBCBase): fenics format bc_expr (fem.Expression): the expression of the particle flux that is used to update the value_fenics - species_dependent_value (dict): a dictionary containing the species that the value. Example: {"name": species} - where "name" is the variable name in the callable value and species is a festim.Species object. + species_dependent_value (dict): a dictionary containing the species + that the value. Example: {"name": species} + where "name" is the variable name in the callable value and species + is a festim.Species object. Examples: @@ -176,11 +183,14 @@ class ParticleFluxBC(FluxBCBase): .. testcode:: ParticleFluxBC ParticleFluxBC(subdomain=my_subdomain, value=1, species="H") - ParticleFluxBC(subdomain=my_subdomain, value=lambda x: 1 + x[0], species="H") + ParticleFluxBC(subdomain=my_subdomain, value=lambda x: 1 + x[0], + species="H") ParticleFluxBC(subdomain=my_subdomain, value=lambda t: 1 + t, species="H") ParticleFluxBC(subdomain=my_subdomain, value=lambda T: 1 + T, species="H") - ParticleFluxBC(subdomain=my_subdomain, value=lambda x, t: 1 + x[0] + t, species="H") - ParticleFluxBC(subdomain=my_subdomain, value=lambda c1: 2 * c1**2, species="H", species_dependent_value={"c1": species1}) + ParticleFluxBC(subdomain=my_subdomain, value=lambda x, t: 1 + x[0] + + t, species="H") + ParticleFluxBC(subdomain=my_subdomain, value=lambda c1: 2 * c1**2, + species="H", species_dependent_value={"c1": species1}) """ def __init__(self, subdomain, value, species, species_dependent_value={}): @@ -191,10 +201,9 @@ def __init__(self, subdomain, value, species, species_dependent_value={}): def create_value_fenics(self, mesh, temperature, t: fem.Constant): """Creates the value of the boundary condition as a fenics object and sets it to - self.value_fenics. - If the value is a constant, it is converted to a fenics.Constant. - If the value is a function of t, it is converted to a fenics.Constant. - Otherwise, it is converted to a ufl Expression + self.value_fenics. If the value is a constant, it is converted to a + fenics.Constant. If the value is a function of t, it is converted to a + fenics.Constant. Otherwise, it is converted to a ufl Expression. Args: mesh (dolfinx.mesh.Mesh) : the mesh @@ -213,7 +222,7 @@ def create_value_fenics(self, mesh, temperature, t: fem.Constant): # only t is an argument if not isinstance(self.value(t=float(t)), (float, int)): raise ValueError( - f"self.value should return a float or an int, not {type(self.value(t=float(t)))} " + f"self.value should return a float or an int, not {type(self.value(t=float(t)))} " # noqa: E501 ) self.value_fenics = F.as_fenics_constant( mesh=mesh, value=self.value(t=float(t)) @@ -239,11 +248,12 @@ def create_value_fenics(self, mesh, temperature, t: fem.Constant): class HeatFluxBC(FluxBCBase): - """ - Heat flux boundary condition class + """Heat flux boundary condition class. + Ensuring the gradient of the solution T at a boundary: -lambda * grad(T) * n = f - where lambda is the thermal conductivity , n is the outwards normal vector of the boundary, f is a function of space and time. + where lambda is the thermal conductivity , n is the outwards normal vector + of the boundary, f is a function of space and time. Args: subdomain (festim.SurfaceSubdomain): the surface subdomain where the heat flux diff --git a/src/festim/boundary_conditions/henrys_bc.py b/src/festim/boundary_conditions/henrys_bc.py index c92083f0f..071ca08e8 100644 --- a/src/festim/boundary_conditions/henrys_bc.py +++ b/src/festim/boundary_conditions/henrys_bc.py @@ -5,14 +5,13 @@ def henrys_law(T, H_0, E_H, pressure): - """Applies the Henry's law to compute the concentration at the boundary""" + """Applies the Henry's law to compute the concentration at the boundary.""" H = H_0 * ufl.exp(-E_H / k_B / T) return H * pressure class HenrysBC(FixedConcentrationBC): - """ - Henrys boundary condition class + """Henrys boundary condition class. c = H * pressure H = H_0 * exp(-E_H / k_B / T) @@ -21,7 +20,8 @@ class HenrysBC(FixedConcentrationBC): subdomain (festim.Subdomain): the subdomain where the boundary condition is applied species (str): the name of the species - H_0 (float or fem.Constant): the Henrys constant pre-exponential factor (H/m3/Pa) + H_0 (float or fem.Constant): the Henrys constant pre-exponential factor + (H/m3/Pa) E_H (float or fem.Constant): the Henrys constant activation energy (eV) pressure (float or callable): the pressure at the boundary (Pa) @@ -30,7 +30,8 @@ class HenrysBC(FixedConcentrationBC): condition is applied value (float or fem.Constant): the value of the boundary condition species (festim.Species or str): the name of the species - H_0 (float or fem.Constant): the Henrys constant pre-exponential factor (H/m3/Pa) + H_0 (float or fem.Constant): the Henrys constant pre-exponential factor + (H/m3/Pa) E_H (float or fem.Constant): the Henrys constant activation energy (eV) pressure (float or callable): the pressure at the boundary (Pa) @@ -43,11 +44,16 @@ class HenrysBC(FixedConcentrationBC): .. testcode:: HenrysBC - HenrysBC(subdomain=my_subdomain, H_0=1e-6, E_H=0.2, pressure=1e5, species="H") - HenrysBC(subdomain=my_subdomain, H_0=1e-6, E_H=0.2, pressure=lambda x: 1e5 + x[0], species="H") - HenrysBC(subdomain=my_subdomain, H_0=1e-6, E_H=0.2, pressure=lambda t: 1e5 + t, species="H") - HenrysBC(subdomain=my_subdomain, H_0=1e-6, E_H=0.2, pressure=lambda T: 1e5 + T, species="H") - HenrysBC(subdomain=my_subdomain, H_0=1e-6, E_H=0.2, pressure=lambda x, t: 1e5 + x[0] + t, species="H") + HenrysBC(subdomain=my_subdomain, H_0=1e-6, E_H=0.2, pressure=1e5, + species="H") + HenrysBC(subdomain=my_subdomain, H_0=1e-6, E_H=0.2, pressure=lambda + x: 1e5 + x[0], species="H") + HenrysBC(subdomain=my_subdomain, H_0=1e-6, E_H=0.2, pressure=lambda + t: 1e5 + t, species="H") + HenrysBC(subdomain=my_subdomain, H_0=1e-6, E_H=0.2, pressure=lambda + T: 1e5 + T, species="H") + HenrysBC(subdomain=my_subdomain, H_0=1e-6, E_H=0.2, pressure=lambda + x, t: 1e5 + x[0] + t, species="H") """ def __init__(self, subdomain, H_0, E_H, pressure, species) -> None: @@ -62,7 +68,7 @@ def __init__(self, subdomain, H_0, E_H, pressure, species) -> None: super().__init__(value=value, species=species, subdomain=subdomain) def create_new_value_function(self): - """Creates a new value function based on the pressure attribute + """Creates a new value function based on the pressure attribute. Raises: ValueError: if the pressure function is not supported diff --git a/src/festim/boundary_conditions/sieverts_bc.py b/src/festim/boundary_conditions/sieverts_bc.py index 911d5dcf9..d15b317b7 100644 --- a/src/festim/boundary_conditions/sieverts_bc.py +++ b/src/festim/boundary_conditions/sieverts_bc.py @@ -5,14 +5,13 @@ def sieverts_law(T, S_0, E_S, pressure): - """Applies the Sieverts law to compute the concentration at the boundary""" + """Applies the Sieverts law to compute the concentration at the boundary.""" S = S_0 * ufl.exp(-E_S / k_B / T) return S * pressure**0.5 class SievertsBC(FixedConcentrationBC): - """ - Sieverts boundary condition class + """Sieverts boundary condition class. c = S * sqrt(pressure) S = S_0 * exp(-E_S / k_B / T) @@ -21,7 +20,8 @@ class SievertsBC(FixedConcentrationBC): subdomain (festim.Subdomain): the subdomain where the boundary condition is applied species (str): the name of the species - S_0 (float or fem.Constant): the Sieverts constant pre-exponential factor (H/m3/Pa0.5) + S_0 (float or fem.Constant): the Sieverts constant pre-exponential + factor (H/m3/Pa0.5) E_S (float or fem.Constant): the Sieverts constant activation energy (eV) pressure (float or callable): the pressure at the boundary (Pa) @@ -30,7 +30,8 @@ class SievertsBC(FixedConcentrationBC): condition is applied value (float or fem.Constant): the value of the boundary condition species (festim.Species or str): the name of the species - S_0 (float or fem.Constant): the Sieverts constant pre-exponential factor (H/m3/Pa0.5) + S_0 (float or fem.Constant): the Sieverts constant pre-exponential + factor (H/m3/Pa0.5) E_S (float or fem.Constant): the Sieverts constant activation energy (eV) pressure (float or callable): the pressure at the boundary (Pa) @@ -43,11 +44,16 @@ class SievertsBC(FixedConcentrationBC): .. testcode:: SievertsBC - SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2, pressure=1e5, species="H") - SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2, pressure=lambda x: 1e5 + x[0], species="H") - SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2, pressure=lambda t: 1e5 + t, species="H") - SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2, pressure=lambda T: 1e5 + T, species="H") - SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2, pressure=lambda x, t: 1e5 + x[0] + t, species="H") + SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2, pressure=1e5, + species="H") + SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2, + pressure=lambda x: 1e5 + x[0], species="H") + SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2, + pressure=lambda t: 1e5 + t, species="H") + SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2, + pressure=lambda T: 1e5 + T, species="H") + SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2, + pressure=lambda x, t: 1e5 + x[0] + t, species="H") """ def __init__(self, subdomain, S_0, E_S, pressure, species) -> None: @@ -62,7 +68,7 @@ def __init__(self, subdomain, S_0, E_S, pressure, species) -> None: super().__init__(value=value, species=species, subdomain=subdomain) def create_new_value_function(self): - """Creates a new value function based on the pressure attribute + """Creates a new value function based on the pressure attribute. Raises: ValueError: if the pressure function is not supported diff --git a/src/festim/boundary_conditions/surface_reaction.py b/src/festim/boundary_conditions/surface_reaction.py index ca8efe92e..f56e4c56e 100644 --- a/src/festim/boundary_conditions/surface_reaction.py +++ b/src/festim/boundary_conditions/surface_reaction.py @@ -6,8 +6,8 @@ class SurfaceReactionBCpartial(ParticleFluxBC): - """Boundary condition representing a surface reaction - A + B <-> C + """Boundary condition representing a surface reaction A + B <-> C. + where A, B are the reactants and C is the product the forward reaction rate is K_r = k_r0 * exp(-E_kr / (k_B * T)) and the backward reaction rate is K_d = k_d0 * exp(-E_kd / (k_B * T)) @@ -79,8 +79,8 @@ def create_value_fenics(self, mesh, temperature, t: fem.Constant): class SurfaceReactionBC: - """Boundary condition representing a surface reaction - A + B <-> C + """Boundary condition representing a surface reaction A + B <-> C. + where A, B are the reactants and C is the product the forward reaction rate is K_r = k_r0 * exp(-E_kr / (k_B * T)) and the backward reaction rate is K_d = k_d0 * exp(-E_kd / (k_B * T)) diff --git a/src/festim/coupled_heat_hydrogen_problem.py b/src/festim/coupled_heat_hydrogen_problem.py index 64ce4cfec..d7dab6f10 100644 --- a/src/festim/coupled_heat_hydrogen_problem.py +++ b/src/festim/coupled_heat_hydrogen_problem.py @@ -11,8 +11,7 @@ class CoupledTransientHeatTransferHydrogenTransport: - """ - Coupled heat transfer and hydrogen transport transient problem + """Coupled heat transfer and hydrogen transport transient problem. Args: heat_problem: the heat transfer problem @@ -37,8 +36,6 @@ class CoupledTransientHeatTransferHydrogenTransport: heat_problem=my_heat_transfer_model, hydrogen_problem=my_h_transport_model, ) - - """ heat_problem: HeatTransferProblem diff --git a/src/festim/exports/average_surface.py b/src/festim/exports/average_surface.py index ac3354d06..da6b46b5e 100644 --- a/src/festim/exports/average_surface.py +++ b/src/festim/exports/average_surface.py @@ -6,7 +6,7 @@ class AverageSurface(SurfaceQuantity): - """Computes the average value of a field on a given surface + """Computes the average value of a field on a given surface. Args: field (festim.Species): species for which the average surface is computed @@ -25,9 +25,8 @@ def title(self): def compute( self, u: fem.Function | ufl.indexed.Indexed, ds: ufl.Measure, entity_maps=None ): - """ - Computes the average value of the field on the defined surface - subdomain, and appends it to the data list + """Computes the average value of the field on the defined surface subdomain, and + appends it to the data list. Args: u: field for which the average value is computed diff --git a/src/festim/exports/average_volume.py b/src/festim/exports/average_volume.py index f2ff01a8d..01feeae9e 100644 --- a/src/festim/exports/average_volume.py +++ b/src/festim/exports/average_volume.py @@ -4,12 +4,13 @@ class AverageVolume(VolumeQuantity): - """Computes the average value of a field in a given volume + """Computes the average value of a field in a given volume. Args: field (festim.Species): species for which the average volume is computed volume (festim.VolumeSubdomain): volume subdomain - filename (str, optional): name of the file to which the average volume is exported + filename (str, optional): name of the file to which the average volume + is exported Attributes: see `festim.VolumeQuantity` @@ -20,10 +21,8 @@ def title(self): return f"Average {self.field.name} volume {self.volume.id}" def compute(self, u, dx, entity_maps=None): - """ - Computes the average value of solution function within the defined volume - subdomain, and appends it to the data list - """ + """Computes the average value of solution function within the defined volume + subdomain, and appends it to the data list.""" self.value = assemble_scalar( u * dx(self.volume.id), entity_maps=entity_maps ) / assemble_scalar(1 * dx(self.volume.id)) diff --git a/src/festim/exports/maximum_surface.py b/src/festim/exports/maximum_surface.py index 25f45d5f7..1a29f5d7e 100644 --- a/src/festim/exports/maximum_surface.py +++ b/src/festim/exports/maximum_surface.py @@ -7,12 +7,13 @@ class MaximumSurface(SurfaceQuantity): - """Computes the maximum value of a field on a given surface + """Computes the maximum value of a field on a given surface. Args: field (festim.Species): species for which the maximum surface is computed surface (festim.SurfaceSubdomain): surface subdomain - filename (str, optional): name of the file to which the maximum surface is exported + filename (str, optional): name of the file to which the maximum surface + is exported Attributes: see `festim.SurfaceQuantity` @@ -23,10 +24,8 @@ def title(self): return f"Maximum {self.field.name} surface {self.surface.id}" def compute(self): - """ - Computes the maximum value of the field on the defined surface - subdomain, and appends it to the data list - """ + """Computes the maximum value of the field on the defined surface subdomain, and + appends it to the data list.""" solution = self.field.post_processing_solution entities = self.facet_meshtags.find(self.surface.id) if isinstance(solution, dolfinx.fem.Function): diff --git a/src/festim/exports/maximum_volume.py b/src/festim/exports/maximum_volume.py index 043b4e284..d3bb257e3 100644 --- a/src/festim/exports/maximum_volume.py +++ b/src/festim/exports/maximum_volume.py @@ -7,16 +7,16 @@ class MaximumVolume(VolumeQuantity): - """Computes the maximum value of a field in a given volume + """Computes the maximum value of a field in a given volume. Args: field (festim.Species): species for which the maximum volume is computed volume (festim.VolumeSubdomain): volume subdomain - filename (str, optional): name of the file to which the maximum volume is exported + filename (str, optional): name of the file to which the maximum volume + is exported Attributes: see `festim.VolumeQuantity` - """ @property @@ -24,10 +24,8 @@ def title(self): return f"Maximum {self.field.name} volume {self.volume.id}" def compute(self): - """ - Computes the maximum value of solution function within the defined volume - subdomain, and appends it to the data list - """ + """Computes the maximum value of solution function within the defined volume + subdomain, and appends it to the data list.""" solution = self.field.post_processing_solution entities = self.volume_meshtags.find(self.volume.id) diff --git a/src/festim/exports/minimum_surface.py b/src/festim/exports/minimum_surface.py index 454f6b763..1d2e934d5 100644 --- a/src/festim/exports/minimum_surface.py +++ b/src/festim/exports/minimum_surface.py @@ -7,12 +7,13 @@ class MinimumSurface(SurfaceQuantity): - """Computes the minimum value of a field on a given surface + """Computes the minimum value of a field on a given surface. Args: field (festim.Species): species for which the minimum surface is computed surface (festim.SurfaceSubdomain): surface subdomain - filename (str, optional): name of the file to which the minimum surface is exported + filename (str, optional): name of the file to which the minimum surface + is exported Attributes: see `festim.SurfaceQuantity` @@ -23,10 +24,8 @@ def title(self): return f"Minimum {self.field.name} surface {self.surface.id}" def compute(self): - """ - Computes the minimum value of the field on the defined surface - subdomain, and appends it to the data list - """ + """Computes the minimum value of the field on the defined surface subdomain, and + appends it to the data list.""" solution = self.field.post_processing_solution entities = self.facet_meshtags.find(self.surface.id) if isinstance(solution, dolfinx.fem.Function): diff --git a/src/festim/exports/minimum_volume.py b/src/festim/exports/minimum_volume.py index 730b3bf29..e17ce5aff 100644 --- a/src/festim/exports/minimum_volume.py +++ b/src/festim/exports/minimum_volume.py @@ -7,12 +7,13 @@ class MinimumVolume(VolumeQuantity): - """Computes the minimum value of a field in a given volume + """Computes the minimum value of a field in a given volume. Args: field (festim.Species): species for which the minimum volume is computed volume (festim.VolumeSubdomain): volume subdomain - filename (str, optional): name of the file to which the minimum volume is exported + filename (str, optional): name of the file to which the minimum volume + is exported Attributes: see `festim.VolumeQuantity` @@ -23,10 +24,8 @@ def title(self): return f"Minimum {self.field.name} volume {self.volume.id}" def compute(self): - """ - Computes the minimum value of solution function within the defined volume - subdomain, and appends it to the data list - """ + """Computes the minimum value of solution function within the defined volume + subdomain, and appends it to the data list.""" solution = self.field.post_processing_solution entities = self.volume_meshtags.find(self.volume.id) diff --git a/src/festim/exports/surface_flux.py b/src/festim/exports/surface_flux.py index fe70e8f8c..b2ff9de21 100644 --- a/src/festim/exports/surface_flux.py +++ b/src/festim/exports/surface_flux.py @@ -8,7 +8,7 @@ class SurfaceFlux(SurfaceQuantity): - """Computes the flux of a field on a given surface + """Computes the flux of a field on a given surface. Args: field: species for which the surface flux is computed @@ -39,7 +39,7 @@ def __init__( def compute( self, u: fem.Function | ufl.indexed.Indexed, ds: ufl.Measure, entity_maps=None ): - """Computes the value of the flux at the surface + """Computes the value of the flux at the surface. Args: u: field for which the flux is computed diff --git a/src/festim/exports/surface_quantity.py b/src/festim/exports/surface_quantity.py index b6f9d2726..515ab368e 100644 --- a/src/festim/exports/surface_quantity.py +++ b/src/festim/exports/surface_quantity.py @@ -5,7 +5,7 @@ class SurfaceQuantity: - """Export SurfaceQuantity + """Export SurfaceQuantity. Args: field: species for which the surface flux is computed @@ -75,8 +75,8 @@ def field(self, value): self._field = value def write(self, t): - """If the filename doesnt exist yet, create it and write the header, - then append the time and value to the file""" + """If the filename doesnt exist yet, create it and write the header, then append + the time and value to the file.""" if self.filename is not None: if self._first_time_export: diff --git a/src/festim/exports/total_surface.py b/src/festim/exports/total_surface.py index d3f51327a..21719575c 100644 --- a/src/festim/exports/total_surface.py +++ b/src/festim/exports/total_surface.py @@ -6,7 +6,7 @@ class TotalSurface(SurfaceQuantity): - """Computes the total value of a field on a given surface + """Computes the total value of a field on a given surface. Args: field (`festim.Species`): species for which the total volume is computed @@ -24,9 +24,8 @@ def title(self): def compute( self, u: fem.Function | ufl.indexed.Indexed, ds: ufl.Measure, entity_maps=None ): - """ - Computes the total value of the field on the defined surface - subdomain, and appends it to the data list + """Computes the total value of the field on the defined surface subdomain, and + appends it to the data list. Args: u: field for which the total value is computed diff --git a/src/festim/exports/total_volume.py b/src/festim/exports/total_volume.py index 312499421..194ce4ec6 100644 --- a/src/festim/exports/total_volume.py +++ b/src/festim/exports/total_volume.py @@ -5,7 +5,7 @@ class TotalVolume(VolumeQuantity): - """Computes the total value of a field in a given volume + """Computes the total value of a field in a given volume. Args: field (festim.Species): species for which the total volume is computed @@ -21,9 +21,8 @@ def title(self): return f"Total {self.field.name} volume {self.volume.id}" def compute(self, u, dx: ufl.Measure, entity_maps=None): - """ - Computes the value of the total volume of the field in the volume subdomain - and appends it to the data list + """Computes the value of the total volume of the field in the volume subdomain + and appends it to the data list. Args: u: field for which the total volume is computed diff --git a/src/festim/exports/volume_quantity.py b/src/festim/exports/volume_quantity.py index aed892a91..1daaaac6a 100644 --- a/src/festim/exports/volume_quantity.py +++ b/src/festim/exports/volume_quantity.py @@ -5,7 +5,7 @@ class VolumeQuantity: - """Export VolumeQuantity + """Export VolumeQuantity. Args: field: species for which the volume quantity is computed @@ -63,8 +63,8 @@ def field(self, value): self._field = value def write(self, t): - """If the filename doesnt exist yet, create it and write the header, - then append the time and value to the file""" + """If the filename doesnt exist yet, create it and write the header, then append + the time and value to the file.""" if self.filename is not None: if self._first_time_export: diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index 7159360fd..87ad6a631 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -8,7 +8,7 @@ class ExportBaseClass: - """Export functions to VTX file + """Export functions to VTX file. Args: filename: The name of the output file @@ -49,7 +49,7 @@ def filename(self): class VTXTemperatureExport(ExportBaseClass): - """Export temperature field functions to VTX file + """Export temperature field functions to VTX file. Args: filename: The name of the output file @@ -74,7 +74,7 @@ def __init__( class VTXSpeciesExport(ExportBaseClass): - """Export species field functions to VTX file + """Export species field functions to VTX file. Args: filename: The name of the output file @@ -117,8 +117,7 @@ def field(self) -> list[Species]: @field.setter def field(self, value: Species | list[Species]): - """ - Update the field to export. + """Update the field to export. Args: value: The species to export @@ -149,9 +148,9 @@ def field(self, value: Species | list[Species]): self._field = val def get_functions(self) -> list[fem.Function]: - """ - Returns list of species for a given subdomain. If using legacy mode, return the - whole species. + """Returns list of species for a given subdomain. + + If using legacy mode, return the whole species. """ legacy_output: bool = False diff --git a/src/festim/exports/xdmf.py b/src/festim/exports/xdmf.py index 727bf0a10..bd86171bc 100644 --- a/src/festim/exports/xdmf.py +++ b/src/festim/exports/xdmf.py @@ -10,7 +10,7 @@ class XDMFExport(ExportBaseClass): - """Export functions to XDMFfile + """Export functions to XDMFfile. Args: filename: The name of the output file @@ -43,20 +43,20 @@ def field(self, value: Species | list[Species]): for element in value: if not isinstance(element, Species): raise TypeError( - f"Each element in the list must be a species, got {type(element)}." + f"Each element in the list must be a species, got {type(element)}." # noqa: E501 ) val = value elif isinstance(value, Species): val = [value] else: raise TypeError( - f"field must be of type festim.Species or a list of festim.Species, got " + f"field must be of type festim.Species or a list of festim.Species, got " # noqa: E501 f"{type(value)}." ) self._field = val def define_writer(self, comm: mpi4py.MPI.Intracomm) -> None: - """Define the writer + """Define the writer. Args: comm (mpi4py.MPI.Intracomm): the MPI communicator @@ -64,7 +64,7 @@ def define_writer(self, comm: mpi4py.MPI.Intracomm) -> None: self._writer = XDMFFile(comm, self.filename, "w") def write(self, t: float): - """Write functions to VTX file + """Write functions to VTX file. Args: t (float): the time of export diff --git a/src/festim/heat_transfer_problem.py b/src/festim/heat_transfer_problem.py index 600b95e6b..710bdcc8b 100644 --- a/src/festim/heat_transfer_problem.py +++ b/src/festim/heat_transfer_problem.py @@ -79,10 +79,12 @@ def initialise(self): self.initialise_exports() def define_function_space(self): - """Creates the function space of the model, creates a mixed element if - model is multispecies. Creates the main solution and previous solution - function u and u_n. Create global DG function spaces of degree 0 and 1 - for the global diffusion coefficient""" + """Creates the function space of the model, creates a mixed element if model is + multispecies. + + Creates the main solution and previous solution function u and u_n. Create + global DG function spaces of degree 0 and 1 for the global diffusion coefficient + """ degree = 1 element_CG = basix.ufl.element( @@ -99,7 +101,7 @@ def define_function_space(self): self.test_function = ufl.TestFunction(self.function_space) def create_dirichletbc_form(self, bc): - """Creates a dirichlet boundary condition form + """Creates a dirichlet boundary condition form. Args: bc (festim.FixedTemperatureBC): the boundary condition @@ -128,7 +130,7 @@ def create_dirichletbc_form(self, bc): return form def create_source_values_fenics(self): - """For each source create the value_fenics""" + """For each source create the value_fenics.""" for source in self.sources: # create value_fenics for all source objects source.value.convert_input_value( @@ -138,7 +140,7 @@ def create_source_values_fenics(self): ) def create_flux_values_fenics(self): - """For each heat flux create the value_fenics""" + """For each heat flux create the value_fenics.""" for bc in self.boundary_conditions: # create value_fenics for all F.HeatFluxBC objects if isinstance(bc, boundary_conditions.HeatFluxBC): @@ -149,8 +151,8 @@ def create_flux_values_fenics(self): ) def create_initial_conditions(self): - """For each initial condition, create the value_fenics and assign it to - the previous solution of the condition's species""" + """For each initial condition, create the value_fenics and assign it to the + previous solution of the condition's species.""" if not self.initial_condition: return @@ -171,7 +173,7 @@ def create_initial_conditions(self): self.u_n.interpolate(self.initial_condition.expr_fenics) def create_formulation(self): - """Creates the formulation of the model""" + """Creates the formulation of the model.""" self.formulation = 0 @@ -216,8 +218,8 @@ def create_formulation(self): ) def initialise_exports(self): - """Defines the export writers of the model, if field is given as - a string, find species object in self.species""" + """Defines the export writers of the model, if field is given as a string, find + species object in self.species.""" for export in self.exports: if isinstance(export, exports.XDMFExport): @@ -233,7 +235,7 @@ def initialise_exports(self): ) def post_processing(self): - """Post processes the model""" + """Post processes the model.""" for export in self.exports: # skip if it isn't time to export diff --git a/src/festim/helpers.py b/src/festim/helpers.py index 94b8ab41b..104a6888c 100644 --- a/src/festim/helpers.py +++ b/src/festim/helpers.py @@ -12,7 +12,7 @@ def as_fenics_constant( value: float | int | fem.Constant, mesh: dolfinx.mesh.Mesh ) -> fem.Constant: - """Converts a value to a dolfinx.Constant + """Converts a value to a dolfinx.Constant. Args: value: the value to convert @@ -41,7 +41,7 @@ def as_mapped_function( temperature: fem.Function | fem.Constant | ufl.core.expr.Expr | None = None, ) -> ufl.core.expr.Expr: """Maps a user given callable function to the mesh, time or temperature within - festim as needed + festim as needed. Args: value: the callable to convert @@ -76,7 +76,7 @@ def as_fenics_interp_expr_and_function( ) -> tuple[fem.Expression, fem.Function]: """Takes a user given callable function, maps the function to the mesh, time or temperature within festim as needed. Then creates the fenics interpolation - expression and function objects + expression and function objects. Args: value: the callable to convert @@ -104,9 +104,8 @@ def as_fenics_interp_expr_and_function( class Value: - """ - A class to handle input values from users and convert them to a relevent fenics - object + """A class to handle input values from users and convert them to a relevent fenics + object. Args: input_value: The value of the user input @@ -119,7 +118,6 @@ class Value: explicit_time_dependent : True if the user input value is explicitly time dependent temperature_dependent : True if the user input value is temperature dependent - """ input_value: ( @@ -177,7 +175,7 @@ def input_value(self, value): @property def explicit_time_dependent(self) -> bool: - """Returns true if the value given is time dependent""" + """Returns true if the value given is time dependent.""" if self.input_value is None: return False if isinstance(self.input_value, fem.Constant | ufl.core.expr.Expr): @@ -190,7 +188,7 @@ def explicit_time_dependent(self) -> bool: @property def temperature_dependent(self) -> bool: - """Returns true if the value given is temperature dependent""" + """Returns true if the value given is temperature dependent.""" if self.input_value is None: return False if isinstance(self.input_value, fem.Constant | ufl.core.expr.Expr): @@ -208,8 +206,8 @@ def convert_input_value( temperature: fem.Function | fem.Constant | ufl.core.expr.Expr | None = None, up_to_ufl_expr: bool | None = False, ): - """Converts a user given value to a relevent fenics object depending - on the type of the value provided + """Converts a user given value to a relevent fenics object depending on the type + of the value provided. Args: function_space: the function space of the fenics object, optional @@ -264,7 +262,7 @@ def convert_input_value( ) def update(self, t: float): - """Updates the value + """Updates the value. Args: t: the time @@ -285,9 +283,11 @@ def update(self, t: float): # Define the appropriate method based on the version if version.parse(dolfinx_version) > version.parse("0.9.0"): + def get_interpolation_points(element): return element.interpolation_points else: + def get_interpolation_points(element): return element.interpolation_points() @@ -324,8 +324,7 @@ def nmm_interpolate( def is_it_time_to_export( times: list | None, current_time: float, atol=0, rtol=1.0e-5 ) -> bool: - """ - Checks if the exported field should be written to a file or not based on the + """Checks if the exported field should be written to a file or not based on the current time and the times in `export.times` After a successful match, the corresponding time is removed from the list to @@ -392,7 +391,7 @@ def SnesMonitor(snes, iter, rnorm): dolfinx.log.log( dolfinx.log.LogLevel.INFO, - f"SNES {iter=} ; {rnorm=:.5e} ({atol=}) ; {relative_residual=:.5e} ({rtol=}) ; {stepsize_rel=:.5e} ({stol=:.5e})", + f"SNES {iter=} ; {rnorm=:.5e} ({atol=}) ; {relative_residual=:.5e} ({rtol=}) ; {stepsize_rel=:.5e} ({stol=:.5e})", # noqa: E501 ) # Update previous xnorm diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 0f30ed0f1..51ceee8c3 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -51,8 +51,7 @@ class HydrogenTransportProblem(problem.ProblemBase): - """ - Hydrogen Transport Problem. + """Hydrogen Transport Problem. Args: mesh: The mesh @@ -133,7 +132,6 @@ class HydrogenTransportProblem(problem.ProblemBase): species=[F.Species(name="H"), F.Species(name="Trap")], ) my_model.initialise() - """ _temperature_as_function: fem.Function @@ -281,7 +279,7 @@ def volume_meshtags(self, value): @property def _unpacked_bcs(self): - """Returns all boundary conditions, including fluxes from surface reactions""" + """Returns all boundary conditions, including fluxes from surface reactions.""" all_boundary_conditions = [] for bc in self.boundary_conditions: if isinstance(bc, boundary_conditions.SurfaceReactionBC): @@ -330,7 +328,7 @@ def initialise(self): self.initialise_exports() def create_implicit_species_value_fenics(self): - """For each implicit species, create the value_fenics""" + """For each implicit species, create the value_fenics.""" for reaction in self.reactions: for reactant in reaction.reactant: if isinstance(reactant, _species.ImplicitSpecies): @@ -340,7 +338,7 @@ def create_implicit_species_value_fenics(self): ) def create_species_from_traps(self): - """Generate a species and reaction per trap defined in self.traps""" + """Generate a species and reaction per trap defined in self.traps.""" for trap in self.traps: trap.create_species_and_reaction() @@ -348,12 +346,13 @@ def create_species_from_traps(self): self.reactions.append(trap.reaction) def define_temperature(self): - """Sets the value of temperature_fenics_value. The type depends on - self.temperature. If self.temperature is a function on t only, create - a fem.Constant. Else, create an dolfinx.fem.Expression (stored in - self.temperature_expr) to be updated, a dolfinx.fem.Function object - is created from the Expression (stored in self.temperature_fenics_value). - Raise a ValueError if temperature is None. + """Sets the value of temperature_fenics_value. + + The type depends on self.temperature. If self.temperature is a function on t + only, create a fem.Constant. Else, create an dolfinx.fem.Expression (stored in + self.temperature_expr) to be updated, a dolfinx.fem.Function object is created + from the Expression (stored in self.temperature_fenics_value). Raise a + ValueError if temperature is None. """ # check if temperature is None if self.temperature is None: @@ -409,8 +408,8 @@ def define_temperature(self): self.temperature_fenics.interpolate(self.temperature_expr) def initialise_exports(self): - """Defines the export writers of the model, if field is given as - a string, find species object in self.species""" + """Defines the export writers of the model, if field is given as a string, find + species object in self.species.""" for export in self.exports: if isinstance(export, exports.ExportBaseClass): @@ -512,9 +511,8 @@ def initialise_exports(self): export.data = [] def _get_temperature_field_as_function(self) -> dolfinx.fem.Function: - """ - Based on the type of the temperature_fenics attribute, converts - it as a Function to be used in VTX export + """Based on the type of the temperature_fenics attribute, converts it as a + Function to be used in VTX export. Returns: the temperature field of the simulation @@ -536,7 +534,7 @@ def _get_temperature_field_as_function(self) -> dolfinx.fem.Function: return temperature_field def define_D_global(self, species): - """Defines the global diffusion coefficient for a given species + """Defines the global diffusion coefficient for a given species. Args: species (F.Species): the species @@ -635,7 +633,7 @@ def define_function_spaces(self, element_degree: int = 1): def assign_functions_to_species(self): """Creates the solution, prev solution, test function and post-processing solution for each species, as well as a collapsed function space for each - species""" + species.""" sub_solutions = list(ufl.split(self.u)) sub_prev_solution = list(ufl.split(self.u_n)) @@ -655,7 +653,7 @@ def assign_functions_to_species(self): spe.test_function = sub_test_functions[idx] def define_boundary_conditions(self): - """Defines the boundary conditions of the model""" + """Defines the boundary conditions of the model.""" for bc in self._unpacked_bcs: if isinstance(bc.species, str): @@ -671,7 +669,7 @@ def define_boundary_conditions(self): super().define_boundary_conditions() def create_dirichletbc_form(self, bc): - """Creates a dirichlet boundary condition form + """Creates a dirichlet boundary condition form. Args: bc (festim.DirichletBC): the boundary condition @@ -711,7 +709,7 @@ def create_dirichletbc_form(self, bc): return form def convert_source_input_values_to_fenics_objects(self): - """For each source create the value_fenics""" + """For each source create the value_fenics.""" for source in self.sources: # create value_fenics for all F.ParticleSource objects if isinstance(source, _source.ParticleSource): @@ -723,7 +721,7 @@ def convert_source_input_values_to_fenics_objects(self): ) def convert_advection_term_to_fenics_objects(self): - """For each advection term convert the input value""" + """For each advection term convert the input value.""" for advec_term in self.advection_terms: advec_term.velocity.convert_input_value( @@ -731,7 +729,7 @@ def convert_advection_term_to_fenics_objects(self): ) def create_flux_values_fenics(self): - """For each particle flux create the value_fenics""" + """For each particle flux create the value_fenics.""" for bc in self.boundary_conditions: # create value_fenics for all F.ParticleFluxBC objects if isinstance(bc, boundary_conditions.ParticleFluxBC): @@ -742,8 +740,8 @@ def create_flux_values_fenics(self): ) def create_initial_conditions(self): - """For each initial condition, create the value_fenics and assign it to - the previous solution of the condition's species""" + """For each initial condition, create the value_fenics and assign it to the + previous solution of the condition's species.""" if len(self.initial_conditions) > 0 and not self.settings.transient: raise ValueError( @@ -772,7 +770,7 @@ def create_initial_conditions(self): ) def create_formulation(self): - """Creates the formulation of the model""" + """Creates the formulation of the model.""" self.formulation = 0 @@ -808,7 +806,7 @@ def create_formulation(self): ) case _: raise NotImplementedError( - f"Unknown coordinate system {self.mesh.coordinate_system!s}" + f"Unknown coordinate system {self.mesh.coordinate_system!s}" # noqa: E501 ) if self.settings.transient: @@ -928,7 +926,7 @@ def update_time_dependent_values(self): advec_term.velocity.update(t=t) def update_post_processing_solutions(self): - """Updates the post-processing solutions of each species""" + """Updates the post-processing solutions of each species.""" for spe in self.species: spe.post_processing_solution.x.array[:] = self.u.x.array[ @@ -936,7 +934,7 @@ def update_post_processing_solutions(self): ] def post_processing(self): - """Post processes the model""" + """Post processes the model.""" self.update_post_processing_solutions() @@ -1057,8 +1055,8 @@ def __init__( surface_to_volume: dict | None = None, petsc_options: dict | None = None, ): - """Class for a multi-material hydrogen transport problem - For other arguments see ``festim.HydrogenTransportProblem``. + """Class for a multi-material hydrogen transport problem For other arguments see + ``festim.HydrogenTransportProblem``. Args: interfaces (list, optional): list of interfaces (``festim.Interface`` @@ -1217,8 +1215,7 @@ def define_temperature(self): subdomain.sub_T = sub_T def create_dirichletbc_form(self, bc: boundary_conditions.FixedConcentrationBC): - """ - Creates the ``value_fenics`` attribute for a given + """Creates the ``value_fenics`` attribute for a given ``festim.FixedConcentrationBC`` and returns the appropriate ``dolfinx.fem.DirichletBC`` object. @@ -1269,8 +1266,8 @@ def create_dirichletbc_form(self, bc: boundary_conditions.FixedConcentrationBC): return form def create_initial_conditions(self): - """For each intial condition, create the value_fenics and assign it to - the previous solution of the condition's species""" + """For each intial condition, create the value_fenics and assign it to the + previous solution of the condition's species.""" for condition in self.initial_conditions: idx = self.species.index(condition.species) @@ -1298,13 +1295,13 @@ def create_initial_conditions(self): def define_function_spaces( self, subdomain: _subdomain.VolumeSubdomain, element_degree=1 ): - """ - Creates appropriate function space and functions for a given subdomain (submesh) - based on the number of species existing in this subdomain. Then stores the - functionspace, the current solution (``u``) and the previous solution (``u_n``) - functions. It also populates the correspondance dicts attributes of the species - (eg. ``species.subdomain_to_solution``, ``species.subdomain_to_test_function``, - etc) for easy access to the right subfunctions, sub-testfunctions etc. + """Creates appropriate function space and functions for a given subdomain + (submesh) based on the number of species existing in this subdomain. Then stores + the functionspace, the current solution (``u``) and the previous solution + (``u_n``) functions. It also populates the correspondance dicts attributes of + the species (eg. ``species.subdomain_to_solution``, + ``species.subdomain_to_test_function``, etc) for easy access to the right + subfunctions, sub-testfunctions etc. Args: subdomain (F.VolumeSubdomain): a subdomain of the geometry @@ -1357,7 +1354,7 @@ def define_function_spaces( species.subdomain_to_post_processing_solution[subdomain].name = name def convert_source_input_values_to_fenics_objects(self): - """For each source create the value_fenics""" + """For each source create the value_fenics.""" for source in self.sources: # create value_fenics for all F.ParticleSource objects if isinstance(source, _source.ParticleSource): @@ -1372,7 +1369,7 @@ def convert_source_input_values_to_fenics_objects(self): ) def convert_advection_term_to_fenics_objects(self): - """For each advection term convert the input value""" + """For each advection term convert the input value.""" for advec_term in self.advection_terms: if isinstance(advec_term, AdvectionTerm): @@ -1388,8 +1385,7 @@ def define_boundary_conditions(self): super().define_boundary_conditions() def create_subdomain_formulation(self, subdomain: _subdomain.VolumeSubdomain): - """ - Creates the variational formulation for each subdomain and stores it in + """Creates the variational formulation for each subdomain and stores it in ``subdomain.F`` Args: @@ -1432,7 +1428,7 @@ def create_subdomain_formulation(self, subdomain: _subdomain.VolumeSubdomain): ) case _: raise ValueError( - f"Unsupported coordinate system {self.mesh.coordinate_system}" + f"Unsupported coordinate system {self.mesh.coordinate_system}" # noqa: E501 ) # add reaction terms @@ -1498,8 +1494,8 @@ def create_subdomain_formulation(self, subdomain: _subdomain.VolumeSubdomain): subdomain.F = form def create_formulation(self): - """ - Takes all the formulations for each subdomain and adds the interface conditions. + """Takes all the formulations for each subdomain and adds the interface + conditions. Finally compute the jacobian matrix and store it in the ``J`` attribute, adds the ``entity_maps`` to the forms and store them in the ``forms`` attribute @@ -1566,11 +1562,12 @@ def create_formulation(self): ) def override_solution_attributes(self, reaction: _reaction.Reaction): - """ - Reaction.reaction_term() relies on the .solution attribute of the species + """Reaction.reaction_term() relies on the .solution attribute of the species however, in the discontinuous class, this attribute doesn't really make sense since there is one solution per subdomain. - Therefore we temporarily override the .solution attribute based on the reactants, + + Therefore we temporarily override the .solution attribute based on the + reactants, products, and `others` if there are implicit species """ list_of_species_to_override = reaction.reactant + reaction.product @@ -1635,7 +1632,7 @@ def create_solver(self): del opts[f"{prefix}{k}"] def create_flux_values_fenics(self): - """For each particle flux create the ``value_fenics`` attribute""" + """For each particle flux create the ``value_fenics`` attribute.""" for bc in self._unpacked_bcs: if isinstance(bc, boundary_conditions.ParticleFluxBC): volume_subdomain = self.surface_to_volume[bc.subdomain] @@ -1807,7 +1804,7 @@ def post_processing(self): export.t.append(float(self.t)) def iterate(self): - """Iterates the model for a given time step""" + """Iterates the model for a given time step.""" if self.show_progress_bar: self.progress_bar.update( min(self.dt.value, abs(self.settings.final_time - self.t.value)) @@ -1823,7 +1820,7 @@ def iterate(self): _ = self.solver.solve() converged_reason = self.solver.solver.getConvergedReason() assert converged_reason > 0, ( - f"Non-linear solver did not converge. Reason code: {converged_reason}. \n See https://petsc.org/release/manualpages/SNES/SNESConvergedReason/ for more information." + f"Non-linear solver did not converge. Reason code: {converged_reason}. \n See https://petsc.org/release/manualpages/SNES/SNESConvergedReason/ for more information." # noqa: E501 ) nb_its = self.solver.solver.getIterationNumber() @@ -1887,7 +1884,7 @@ def initialise(self): super().initialise() def create_formulation(self): - """Creates the formulation of the model""" + """Creates the formulation of the model.""" self.formulation = 0 @@ -1962,7 +1959,7 @@ def create_formulation(self): ) def add_reaction_term(self, reaction: _reaction.Reaction): - """Adds the reaction term to the formulation""" + """Adds the reaction term to the formulation.""" products = ( reaction.product @@ -2043,7 +2040,7 @@ def override_post_processing_solution(self): ) # NOTE: do we need this line since it's in initialise? def update_post_processing_solutions(self): - """Updates the post-processing solutions after each time step""" + """Updates the post-processing solutions after each time step.""" # need to compute c = theta * K_S # this expression is stored in species.dg_expr @@ -2053,7 +2050,7 @@ def update_post_processing_solutions(self): spe.post_processing_solution.interpolate(spe.dg_expr) def create_dirichletbc_form(self, bc: boundary_conditions.FixedConcentrationBC): - """Creates a dirichlet boundary condition form + """Creates a dirichlet boundary condition form. Args: bc (festim.DirichletBC): the boundary condition diff --git a/src/festim/initial_condition.py b/src/festim/initial_condition.py index 7153dffc4..e10768a80 100644 --- a/src/festim/initial_condition.py +++ b/src/festim/initial_condition.py @@ -15,8 +15,7 @@ class InitialConditionBase: - """ - Base initial condition class + """Base initial condition class. Args: value: the value of the initial condition. @@ -67,8 +66,7 @@ def volume(self, value): class InitialConcentration(InitialConditionBase): - """ - Initial concentration class + """Initial concentration class. Args: value: the value of the initial concentration of a given species. @@ -172,8 +170,7 @@ def create_expr_fenics( class InitialTemperature(InitialConditionBase): - """ - Initial temperature class + """Initial temperature class. Args: value: the value of the initial temperature @@ -247,8 +244,7 @@ def read_function_from_file( order: int = 1, mesh: dolfinx.mesh.Mesh | None = None, ) -> fem.Function: - """ - Read a function from a file + """Read a function from a file. note:: The function is read from a file using adios4dolfinx. For more information diff --git a/src/festim/material.py b/src/festim/material.py index c66a9cb97..1ffff3dc0 100644 --- a/src/festim/material.py +++ b/src/festim/material.py @@ -29,8 +29,7 @@ def from_string(cls, s: str): class Material: - """ - Material class + """Material class. Args: D_0: the pre-exponential factor of the @@ -147,7 +146,7 @@ def D(self, value): raise TypeError("D must be of type fem.Function") def get_D_0(self, species=None): - """Returns the pre-exponential factor of the diffusion coefficient + """Returns the pre-exponential factor of the diffusion coefficient. Args: species (festim.Species or str, optional): the species we want the @@ -176,7 +175,7 @@ def get_D_0(self, species=None): raise TypeError("D_0 must be either a float, int or a dict") def get_E_D(self, species=None): - """Returns the activation energy of the diffusion coefficient + """Returns the activation energy of the diffusion coefficient. Args: species (festim.Species or str, optional): the species we want the @@ -205,7 +204,7 @@ def get_E_D(self, species=None): raise TypeError("E_D must be either a float, int or a dict") def get_K_S_0(self, species=None) -> float: - """Returns the pre-exponential factor of the solubility coefficient + """Returns the pre-exponential factor of the solubility coefficient. Args: species: the species we want the pre-exponential @@ -233,7 +232,7 @@ def get_K_S_0(self, species=None) -> float: raise TypeError("K_S_0 must be either a float, int or a dict") def get_E_K_S(self, species=None) -> float: - """Returns the activation energy of the solubility coefficient + """Returns the activation energy of the solubility coefficient. Args: species: the species we want the activation @@ -261,7 +260,7 @@ def get_E_K_S(self, species=None) -> float: raise TypeError("E_K_S must be either a float, int or a dict") def get_diffusion_coefficient(self, mesh=None, temperature=None, species=None): - """Defines the diffusion coefficient + """Defines the diffusion coefficient. Args: @@ -319,7 +318,7 @@ def get_diffusion_coefficient(self, mesh=None, temperature=None, species=None): raise ValueError("D_0 and E_D must be either floats or dicts") def get_solubility_coefficient(self, mesh, temperature, species=None): - """Defines the solubility coefficient + """Defines the solubility coefficient. Args: diff --git a/src/festim/mesh/mesh.py b/src/festim/mesh/mesh.py index d08d4cab6..76fa9003f 100644 --- a/src/festim/mesh/mesh.py +++ b/src/festim/mesh/mesh.py @@ -26,7 +26,7 @@ def from_string(cls, s: str): return cls.SPHERICAL else: raise ValueError( - "coordinate_system must be one of 'cartesian', 'cylindrical', or 'spherical'" + "coordinate_system must be one of 'cartesian', 'cylindrical', or 'spherical'" # noqa: E501 ) def __str__(self): @@ -41,8 +41,7 @@ def __str__(self): class Mesh: - """ - Mesh class + """Mesh class. Args: mesh: The mesh. Defaults to None. @@ -129,11 +128,13 @@ def n(self): return ufl.FacetNormal(self._mesh) def define_meshtags(self, surface_subdomains, volume_subdomains, interfaces=None): - """Defines the facet and volume meshtags of the mesh + """Defines the facet and volume meshtags of the mesh. Args: - surface_subdomains (list of festim.SufaceSubdomains): the surface subdomains of the model - volume_subdomains (list of festim.VolumeSubdomains): the volume subdomains of the model + surface_subdomains (list of festim.SufaceSubdomains): the surface + subdomains of the model + volume_subdomains (list of festim.VolumeSubdomains): the volume + subdomains of the model interfaces (dict, optional): the interfaces between volume subdomains {int: [VolumeSubdomain, VolumeSubdomain]}. Defaults to None. @@ -219,7 +220,7 @@ def define_meshtags(self, surface_subdomains, volume_subdomains, interfaces=None def check_mesh_dim_coords(self): """Checks if the used coordinates can be applied for geometry with the specified - dimensions""" + dimensions.""" if self.coordinate_system == CoordinateSystem.SPHERICAL and self.vdim != 1: raise AttributeError( diff --git a/src/festim/mesh/mesh_1d.py b/src/festim/mesh/mesh_1d.py index d07d2c58e..45b4b66a5 100644 --- a/src/festim/mesh/mesh_1d.py +++ b/src/festim/mesh/mesh_1d.py @@ -9,8 +9,7 @@ class Mesh1D(Mesh): - """ - 1D Mesh + """1D Mesh. Args: vertices (list or np.ndarray): the mesh x-coordinates (m) @@ -34,7 +33,7 @@ def vertices(self, value): self._vertices = np.sort(np.unique(value)).astype(np.float64) def generate_mesh(self): - """Generates a 1D mesh""" + """Generates a 1D mesh.""" if MPI.COMM_WORLD.rank == 0: mesh_points = np.reshape(self.vertices, (len(self.vertices), 1)) @@ -54,7 +53,7 @@ def generate_mesh(self): ) def check_borders(self, volume_subdomains): - """Checks that the borders of the subdomain are within the domain + """Checks that the borders of the subdomain are within the domain. Args: volume_subdomains (list of festim.VolumeSubdomain1D): the volume subdomains diff --git a/src/festim/mesh/mesh_from_xdmf.py b/src/festim/mesh/mesh_from_xdmf.py index 9d0d20a1d..a318f06bc 100644 --- a/src/festim/mesh/mesh_from_xdmf.py +++ b/src/festim/mesh/mesh_from_xdmf.py @@ -6,8 +6,7 @@ class MeshFromXDMF(Mesh): - """ - Mesh read from the XDMF files + """Mesh read from the XDMF files. Args: volume_file (str): path to the volume file @@ -23,7 +22,8 @@ class MeshFromXDMF(Mesh): volume_file (str): path to the volume file facet_file (str): path to the facet file mesh_name (str): name of the mesh in the volume XDMF file. - surface_meshtags_name (str): name of the surface meshtags in the facet XDMF file. + surface_meshtags_name (str): name of the surface meshtags in the facet + XDMF file. volume_meshtags_name (str): name of the volume meshtags in the volume XDMF file mesh (fenics.mesh.Mesh): the fenics mesh """ @@ -48,7 +48,7 @@ def __init__( super().__init__(mesh=mesh) def define_surface_meshtags(self): - """Creates the facet meshtags + """Creates the facet meshtags. Returns: dolfinx.MeshTags: the facet meshtags @@ -61,7 +61,7 @@ def define_surface_meshtags(self): return facet_meshtags def define_volume_meshtags(self): - """Creates the volume meshtags + """Creates the volume meshtags. Returns: dolfinx.MeshTags: the volume meshtags diff --git a/src/festim/problem.py b/src/festim/problem.py index bfd9852cd..675ef7474 100644 --- a/src/festim/problem.py +++ b/src/festim/problem.py @@ -21,9 +21,10 @@ class ProblemBase: - """ - Base class for :py:class:`HeatTransferProblem ` and - :py:class:`HydrogenTransportProblem `. + """Base class for :py:class:`HeatTransferProblem + ` and + :py:class:`HydrogenTransportProblem + `. Attributes: show_progress_bar: If `True` a progress bar is displayed during the simulation @@ -85,8 +86,8 @@ def timesteps(self): return self._timesteps def define_meshtags_and_measures(self): - """Defines the facet and volume meshtags of the model which are used - to define the measures fo the model, dx and ds""" + """Defines the facet and volume meshtags of the model which are used to define + the measures fo the model, dx and ds.""" if isinstance(self.mesh, F.MeshFromXDMF): # TODO: fix naming inconsistency between facet and surface meshtags @@ -119,16 +120,15 @@ def define_meshtags_and_measures(self): ) def define_boundary_conditions(self): - """Defines the dirichlet boundary conditions of the model""" + """Defines the dirichlet boundary conditions of the model.""" for bc in self.boundary_conditions: if isinstance(bc, F.DirichletBCBase): form = self.create_dirichletbc_form(bc) self.bc_forms.append(form) def get_petsc_options(self) -> dict[str, Any]: - """ - Gets the PETSc options to pass to the NewtonProblem solver. Default - options are updated with user-provided options, if any. + """Gets the PETSc options to pass to the NewtonProblem solver. Default options + are updated with user-provided options, if any. Returns: the petsc options to pass to the NewtonProblem solver. @@ -165,7 +165,7 @@ def get_petsc_options(self) -> dict[str, Any]: return petsc_options def create_solver(self): - """Creates the solver of the model""" + """Creates the solver of the model.""" if Version(dolfinx.__version__) == Version("0.9.0"): problem = fem.petsc.NonlinearProblem( @@ -226,7 +226,7 @@ def create_solver(self): del opts[f"{prefix}{k}"] def run(self): - """Runs the model""" + """Runs the model.""" if self.settings.transient: # Solve transient @@ -250,7 +250,7 @@ def run(self): self.post_processing() def iterate(self): - """Iterates the model for a given time step""" + """Iterates the model for a given time step.""" self._timesteps.append(float(self.t)) if self.show_progress_bar: @@ -276,7 +276,7 @@ def iterate(self): _ = self.solver.solve() converged_reason = self.solver.solver.getConvergedReason() assert converged_reason > 0, ( - f"Non-linear solver did not converge. Reason code: {converged_reason}. \n See https://petsc.org/release/manualpages/SNES/SNESConvergedReason/ for more information." + f"Non-linear solver did not converge. Reason code: {converged_reason}. \n See https://petsc.org/release/manualpages/SNES/SNESConvergedReason/ for more information." # noqa: E501 ) nb_its = self.solver.solver.getIterationNumber() diff --git a/src/festim/reaction.py b/src/festim/reaction.py index aed8a70eb..946bfc7e2 100644 --- a/src/festim/reaction.py +++ b/src/festim/reaction.py @@ -13,22 +13,26 @@ class Reaction: """A reaction between two species, with a forward and backward rate. Arguments: - reactant (Union[F.Species, F.ImplicitSpecies], List[Union[F.Species, F.ImplicitSpecies]]): The reactant. + reactant (Union[F.Species, F.ImplicitSpecies], List[Union[F.Species, + F.ImplicitSpecies]]): The reactant. product (Optional[Union[F.Species, List[F.Species]]]): The product. k_0 (float): The forward rate constant pre-exponential factor. E_k (float): The forward rate constant activation energy. p_0 (float): The backward rate constant pre-exponential factor. E_p (float): The backward rate constant activation energy. - volume (F.VolumeSubdomain1D): The volume subdomain where the reaction takes place. + volume (F.VolumeSubdomain1D): The volume subdomain where the reaction + takes place. Attributes: - reactant (Union[F.Species, F.ImplicitSpecies], List[Union[F.Species, F.ImplicitSpecies]]): The reactant. + reactant (Union[F.Species, F.ImplicitSpecies], List[Union[F.Species, + F.ImplicitSpecies]]): The reactant. product (Optional[Union[F.Species, List[F.Species]]]): The product. k_0 (float): The forward rate constant pre-exponential factor. E_k (float): The forward rate constant activation energy. p_0 (float): The backward rate constant pre-exponential factor. E_p (float): The backward rate constant activation energy. - volume (F.VolumeSubdomain1D): The volume subdomain where the reaction takes place. + volume (F.VolumeSubdomain1D): The volume subdomain where the reaction + takes place. Examples: @@ -53,7 +57,6 @@ class Reaction: # compute the reaction term at a given temperature temperature = 300.0 reaction_term = reaction.reaction_term(temperature) - """ def __init__( @@ -84,7 +87,7 @@ def reactant(self, value): value = [value] if len(value) == 0: raise ValueError( - "reactant must be an entry of one or more species objects, not an empty list." + "reactant must be an entry of one or more species objects, not an empty list." # noqa: E501 ) for i in value: if not isinstance(i, (_Species, _ImplicitSpecies)): @@ -101,7 +104,7 @@ def __repr__(self) -> str: products = " + ".join([str(product) for product in self.product]) else: products = self.product - return f"Reaction({reactants} <--> {products}, {self.k_0}, {self.E_k}, {self.p_0}, {self.E_p})" + return f"Reaction({reactants} <--> {products}, {self.k_0}, {self.E_k}, {self.p_0}, {self.E_p})" # noqa: E501 def __str__(self) -> str: reactants = " + ".join([str(reactant) for reactant in self.reactant]) diff --git a/src/festim/settings.py b/src/festim/settings.py index 16fdbdad6..c3652211b 100644 --- a/src/festim/settings.py +++ b/src/festim/settings.py @@ -30,7 +30,6 @@ class Settings: stepsize (festim.Stepsize): stepsize for a transient simulation. convergence_criterion: resiudal or incremental (for Newton solver) - """ def __init__( diff --git a/src/festim/source.py b/src/festim/source.py index 054bbde1e..9d8dc23f6 100644 --- a/src/festim/source.py +++ b/src/festim/source.py @@ -8,8 +8,7 @@ class SourceBase: - """ - Source base class + """Source base class. Args: value: the value of the source @@ -72,8 +71,7 @@ def volume(self, value): class ParticleSource(SourceBase): - """ - Particle source class + """Particle source class. Args: value: the value of the source @@ -119,8 +117,7 @@ def species(self, value): class HeatSource(SourceBase): - """ - Heat source class + """Heat source class. Args: value: the value of the source diff --git a/src/festim/species.py b/src/festim/species.py index 9b9aa412b..8f8d8babb 100644 --- a/src/festim/species.py +++ b/src/festim/species.py @@ -10,8 +10,7 @@ class Species: - """ - Hydrogen species class for H transport simulation. + """Hydrogen species class for H transport simulation. Args: name: a name given to the species. Defaults to None. @@ -53,8 +52,6 @@ class Species: Species(name="H") Species(name="Trap", mobile=False) - - """ name: str | None @@ -113,9 +110,7 @@ def concentration(self): @property def legacy(self) -> bool: - """ - Check if we are using FESTIM 1.0 implementation or FESTIM 2.0 - """ + """Check if we are using FESTIM 1.0 implementation or FESTIM 2.0.""" if not self.subdomain_to_solution: return True else: @@ -176,10 +171,9 @@ def concentration(self): def create_value_fenics(self, mesh, t: fem.Constant): """Creates the value of the density as a fenics object and sets it to - self.value_fenics. - If the value is a constant, it is converted to a fenics.Constant. - If the value is a function of t, it is converted to a fenics.Constant. - Otherwise, it is converted to a ufl Expression + self.value_fenics. If the value is a constant, it is converted to a + fenics.Constant. If the value is a function of t, it is converted to a + fenics.Constant. Otherwise, it is converted to a ufl Expression. Args: mesh (dolfinx.mesh.Mesh) : the mesh @@ -231,8 +225,7 @@ def update_density(self, t): def find_species_from_name(name: str, species: list): - """Returns the correct species object from a list of species - based on a string + """Returns the correct species object from a list of species based on a string. Args: name (str): the name of the species @@ -243,7 +236,6 @@ def find_species_from_name(name: str, species: list): Raises: ValueError: if the species name is not found in the list of species - """ for spe in species: if spe.name == name: diff --git a/src/festim/stepsize.py b/src/festim/stepsize.py index e555e6594..e8f92512c 100644 --- a/src/festim/stepsize.py +++ b/src/festim/stepsize.py @@ -2,8 +2,7 @@ class Stepsize: - """ - A class for evaluating the stepsize of transient simulations. + """A class for evaluating the stepsize of transient simulations. Args: initial_value (float, int): initial stepsize. @@ -113,8 +112,7 @@ def max_stepsize(self, value): self._max_stepsize = value def get_max_stepsize(self, t): - """ - Returns the maximum stepsize at time t. + """Returns the maximum stepsize at time t. Args: t (float): the current time @@ -152,9 +150,7 @@ def modify_value(self, value, nb_iterations, t=None): return updated_value def is_adapt(self, t): - """ - Methods that defines if the stepsize should be - adapted or not + """Methods that defines if the stepsize should be adapted or not. Args: t (float): the current time @@ -165,8 +161,8 @@ def is_adapt(self, t): return True def next_milestone(self, current_time: float): - """Returns the next milestone that the simulation must pass. - Returns None if there are no more milestones. + """Returns the next milestone that the simulation must pass. Returns None if + there are no more milestones. Args: current_time (float): current time. diff --git a/src/festim/subdomain/interface.py b/src/festim/subdomain/interface.py index 02060d31c..8b568d40e 100644 --- a/src/festim/subdomain/interface.py +++ b/src/festim/subdomain/interface.py @@ -54,9 +54,8 @@ def from_string(cls, s: str) -> "InterfaceMethod": class InterfaceBase(ABC): """Abstract base class for interfaces between subdomains. - Provides common functionality for handling interfaces in discontinuous - finite element problems, including integration data computation and - restriction handling. + Provides common functionality for handling interfaces in discontinuous finite + element problems, including integration data computation and restriction handling. """ def __init__( @@ -76,9 +75,9 @@ def __init__( def pad_parent_maps(self): """Pad parent-to-submesh maps for correct sparsity pattern. - This is a workaround to ensure the sparsity pattern is correct when - assembling forms with interface integrals. It pads the mapping between - parent mesh cells and submesh cells for DOLFINx versions that require it. + This is a workaround to ensure the sparsity pattern is correct when assembling + forms with interface integrals. It pads the mapping between parent mesh cells + and submesh cells for DOLFINx versions that require it. """ try: # No padding needed for latest version of DOLFINx diff --git a/src/festim/subdomain/surface_subdomain.py b/src/festim/subdomain/surface_subdomain.py index bdf673374..5a7e2803d 100644 --- a/src/festim/subdomain/surface_subdomain.py +++ b/src/festim/subdomain/surface_subdomain.py @@ -5,8 +5,7 @@ class SurfaceSubdomain: - """ - Surface subdomain class + """Surface subdomain class. Args: id: the id of the surface subdomain @@ -21,9 +20,12 @@ class SurfaceSubdomain: .. testcode:: SurfaceSubdomain SurfaceSubdomain(id=1, locator=lambda x: np.isclose(x[0], 1.0)) - SurfaceSubdomain(id=1, locator=lambda x: np.logical_or(np.isclose(x[1], 0.0), np.isclose(x[1], 1.0))) - SurfaceSubdomain(id=1, locator=lambda x: np.logical_and(np.isclose(x[0], 0.0), np.isclose(x[1], 1.0))) - SurfaceSubdomain(id=1, locator=lambda x: np.logical_and(np.isclose(x[0], 0.0), x[1] <= 0.5)) + SurfaceSubdomain(id=1, locator=lambda x: + np.logical_or(np.isclose(x[1], 0.0), np.isclose(x[1], 1.0))) + SurfaceSubdomain(id=1, locator=lambda x: + np.logical_and(np.isclose(x[0], 0.0), np.isclose(x[1], 1.0))) + SurfaceSubdomain(id=1, locator=lambda x: + np.logical_and(np.isclose(x[0], 0.0), x[1] <= 0.5)) """ id: int @@ -56,8 +58,7 @@ def locate_boundary_facet_indices(self, mesh: dolfinx.mesh.Mesh) -> np.ndarray: class SurfaceSubdomain1D(SurfaceSubdomain): - """ - Surface subdomain class for 1D cases + """Surface subdomain class for 1D cases. Args: id: the id of the surface subdomain @@ -88,8 +89,8 @@ def __init__(self, id: int, x: float) -> None: def find_surface_from_id(id: int, surfaces: list): - """Returns the correct surface subdomain object from a list of surface ids - based on an int + """Returns the correct surface subdomain object from a list of surface ids based on + an int. Args: id (int): the id of the surface subdomain @@ -100,7 +101,6 @@ def find_surface_from_id(id: int, surfaces: list): Raises: ValueError: if the surface name is not found in the list of surfaces - """ for surf in surfaces: if surf.id == id: diff --git a/src/festim/subdomain/volume_subdomain.py b/src/festim/subdomain/volume_subdomain.py index 22a735b00..3815b5ebc 100644 --- a/src/festim/subdomain/volume_subdomain.py +++ b/src/festim/subdomain/volume_subdomain.py @@ -20,8 +20,7 @@ class VolumeSubdomain: - """ - Volume subdomain class + """Volume subdomain class. Args: id: the id of the volume subdomain (> 0) @@ -116,7 +115,7 @@ def transfer_meshtag(self, mesh: dolfinx.mesh.Mesh, tag: dolfinx.mesh.MeshTags): ) def locate_subdomain_entities(self, mesh: Mesh) -> npt.NDArray[np.int32]: - """Locates all cells in subdomain borders within domain + """Locates all cells in subdomain borders within domain. Args: mesh: the mesh of the model @@ -132,8 +131,7 @@ def locate_subdomain_entities(self, mesh: Mesh) -> npt.NDArray[np.int32]: class VolumeSubdomain1D(VolumeSubdomain): - """ - Volume subdomain class for 1D cases + """Volume subdomain class for 1D cases. Args: id (int): the id of the volume subdomain @@ -167,8 +165,8 @@ def __init__(self, id, borders, material) -> None: def find_volume_from_id(id: int, volumes: list): - """Returns the correct volume subdomain object from a list of volume ids - based on an int + """Returns the correct volume subdomain object from a list of volume ids based on an + int. Args: id (int): the id of the volume subdomain @@ -179,7 +177,6 @@ def find_volume_from_id(id: int, volumes: list): Raises: ValueError: if the volume name is not found in the list of volumes - """ for vol in volumes: if vol.id == id: diff --git a/src/festim/trap.py b/src/festim/trap.py index 0e5553a32..809baf333 100644 --- a/src/festim/trap.py +++ b/src/festim/trap.py @@ -28,13 +28,15 @@ class Trap(_Species): volume (F.VolumeSubdomain1D): The volume subdomain where the trap is. trapped_concentration (_Species): The immobile trapped concentration trap_reaction (_Reaction): The reaction for trapping the mobile conc. - empty_trap_sites (F.ImplicitSpecies): The implicit species for the empty trap sites + empty_trap_sites (F.ImplicitSpecies): The implicit species for the + empty trap sites Examples: .. testsetup:: Trap - from festim import Trap, Species, VolumeSubdomain, Material, HydrogenTransportProblem + from festim import Trap, Species, VolumeSubdomain, Material, + HydrogenTransportProblem my_mat = Material(D_0=1, E_D=1, name="test_mat") my_vol = VolumeSubdomain(id=1, material=my_mat) @@ -42,7 +44,8 @@ class Trap(_Species): .. testcode:: Trap - trap = Trap(name="Trap", mobile_species=H, k_0=1.0, E_k=0.2, p_0=0.1, E_p=0.3, n=100, volume=my_vol) + trap = Trap(name="Trap", mobile_species=H, k_0=1.0, E_k=0.2, + p_0=0.1, E_p=0.3, n=100, volume=my_vol) my_model = HydrogenTransportProblem() my_model.traps = [trap] @@ -85,8 +88,6 @@ class Trap(_Species): ) my_model.species = [cm, ct] my_model.reactions = [trap_reaction] - - """ def __init__( @@ -105,7 +106,7 @@ def __init__( self.reaction = None def create_species_and_reaction(self): - """create the immobile trapped species object and the reaction for trapping""" + """Create the immobile trapped species object and the reaction for trapping.""" self.trapped_concentration = _Species(name=self.name, mobile=False) self.empty_trap_sites = _ImplicitSpecies( n=self.n, others=[self.trapped_concentration] diff --git a/test/benchmark.py b/test/benchmark.py index 64dfd99f2..387907484 100644 --- a/test/benchmark.py +++ b/test/benchmark.py @@ -187,8 +187,8 @@ def siverts_law(T, S_0, E_S, pressure): def test_festim_vs_fenics_permeation_benchmark(): - """Runs a problem with pure fenicsx and the same problem with FESTIM and - raise ValueError if difference is too high""" + """Runs a problem with pure fenicsx and the same problem with FESTIM and raise + ValueError if difference is too high.""" start = time.time() fenics_test_permeation_problem(mesh_size=20001) diff --git a/test/system_tests/test_1_mobile_1_trap.py b/test/system_tests/test_1_mobile_1_trap.py index 4942bb5b3..f08c0152d 100644 --- a/test/system_tests/test_1_mobile_1_trap.py +++ b/test/system_tests/test_1_mobile_1_trap.py @@ -18,9 +18,7 @@ def test_1_mobile_1_trap_MMS_steady_state(): - """ - MMS test with 1 mobile species and 1 trap at steady state - """ + """MMS test with 1 mobile species and 1 trap at steady state.""" def u_exact(mod): return lambda x: 1.5 + mod.sin(3 * mod.pi * x[0]) @@ -42,8 +40,10 @@ def v_exact(mod): E_k = 1.5 p_0 = 0.2 E_p = 0.1 + def T_expr(x): return 500 + 100 * x[0] + T.interpolate(T_expr) n_trap = 3 E_D = 0.1 @@ -118,10 +118,8 @@ def T_expr(x): def test_1_mobile_1_trap_MMS_transient(): - """ - MMS test with 1 mobile species and 1 trap in 0.1s transient, the value at the last time step is - compared to an analytical solution - """ + """MMS test with 1 mobile species and 1 trap in 0.1s transient, the value at the + last time step is compared to an analytical solution.""" final_time = 0.1 @@ -139,8 +137,10 @@ def u_exact_alt(mod): D_0 = 1 E_D = 0.1 + def T_expr(x): return 600 + 50 * x[0] + T.interpolate(T_expr) D = D_0 * ufl.exp(-E_D / (F.k_B * T)) @@ -167,12 +167,14 @@ def T_expr(x): def init_value(x): return 1 + ufl.sin(2 * ufl.pi * x[0]) + my_model.initial_conditions = [ F.InitialConcentration(value=init_value, species=H, volume=vol) ] def f(x, t): return 4 * t - ufl.div(D * ufl.grad(H_analytical_ufl(x_1d, t))) + my_model.sources = [F.ParticleSource(value=f, volume=vol, species=H)] my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=final_time) @@ -189,8 +191,8 @@ def f(x, t): def test_1_mobile_1_trap_MMS_3D(): - """Tests that a steady simulation can be run in a 3D domain with - 1 mobile and 1 trapped species""" + """Tests that a steady simulation can be run in a 3D domain with 1 mobile and 1 + trapped species.""" def u_exact(mod): return lambda x: 1.5 + mod.sin(3 * mod.pi * x[0]) + mod.cos(3 * mod.pi * x[1]) @@ -212,8 +214,10 @@ def v_exact(mod): E_k = 1.5 p_0 = 0.2 E_p = 0.1 + def T_expr(x): return 500 + 100 * x[0] + T.interpolate(T_expr) n_trap = 3 E_D = 0.1 @@ -303,8 +307,8 @@ def locate_boundary_facet_indices(self, mesh): def test_1_mobile_1_trap_MMS_2D(): - """Tests that a steady simulation can be run in a 2D domain with - 1 mobile and 1 trapped species""" + """Tests that a steady simulation can be run in a 2D domain with 1 mobile and 1 + trapped species.""" def u_exact(mod): return lambda x: 1.5 + mod.sin(3 * mod.pi * x[0]) + mod.cos(3 * mod.pi * x[1]) @@ -326,8 +330,10 @@ def v_exact(mod): E_k = 1.5 p_0 = 0.2 E_p = 0.1 + def T_expr(x): return 500 + 100 * x[0] + T.interpolate(T_expr) n_trap = 3 E_D = 0.1 diff --git a/test/system_tests/test_1_mobile_species.py b/test/system_tests/test_1_mobile_species.py index 84925a54d..b75f3334e 100644 --- a/test/system_tests/test_1_mobile_species.py +++ b/test/system_tests/test_1_mobile_species.py @@ -18,9 +18,7 @@ def test_1_mobile_MMS_steady_state(): - """ - MMS test with one mobile species at steady state - """ + """MMS test with one mobile species at steady state.""" def u_exact(mod): return lambda x: 1 + mod.sin(2 * mod.pi * x[0]) @@ -33,8 +31,10 @@ def u_exact(mod): D_0 = 1 E_D = 0.1 + def T_expr(x): return 500 + 100 * x[0] + T.interpolate(T_expr) D = D_0 * ufl.exp(-E_D / (F.k_B * T)) @@ -73,10 +73,8 @@ def T_expr(x): def test_1_mobile_MMS_transient(): - """ - MMS test with 1 mobile species in 0.1s transient, the value at the last time step is - compared to an analytical solution - """ + """MMS test with 1 mobile species in 0.1s transient, the value at the last time step + is compared to an analytical solution.""" final_time = 0.1 @@ -94,8 +92,10 @@ def u_exact_alt(mod): D_0 = 1 E_D = 0.1 + def T_expr(x): return 600 + 50 * x[0] + T.interpolate(T_expr) D = D_0 * ufl.exp(-E_D / (F.k_B * T)) @@ -122,12 +122,14 @@ def T_expr(x): def init_value(x): return 1 + ufl.sin(2 * ufl.pi * x[0]) + my_model.initial_conditions = [ F.InitialConcentration(value=init_value, species=H, volume=vol) ] def f(x, t): return 4 * t - ufl.div(D * ufl.grad(H_analytical_ufl(x, t))) + my_model.sources = [F.ParticleSource(value=f, volume=vol, species=H)] my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=final_time) @@ -144,8 +146,8 @@ def f(x, t): def test_1_mobile_MMS_2D(): - """Tests that a steady simulation can be run in a 2D domain with - 1 mobile species""" + """Tests that a steady simulation can be run in a 2D domain with 1 mobile + species.""" def u_exact(mod): return lambda x: 1 + mod.sin(2 * mod.pi * x[0]) + mod.cos(2 * mod.pi * x[1]) @@ -158,8 +160,10 @@ def u_exact(mod): D_0 = 1 E_D = 0.1 + def T_expr(x): return 500 + 100 * x[0] + T.interpolate(T_expr) D = D_0 * ufl.exp(-E_D / (F.k_B * T)) @@ -210,8 +214,8 @@ def locate_boundary_facet_indices(self, mesh): def test_1_mobile_MMS_3D(): - """Tests that a steady simulation can be run in a 3D domain with - 1 mobile species""" + """Tests that a steady simulation can be run in a 3D domain with 1 mobile + species.""" def u_exact(mod): return lambda x: 1 + mod.sin(2 * mod.pi * x[0]) + mod.cos(2 * mod.pi * x[1]) @@ -224,8 +228,10 @@ def u_exact(mod): D_0 = 1 E_D = 0.1 + def T_expr(x): return 500 + 100 * x[0] + T.interpolate(T_expr) D = D_0 * ufl.exp(-E_D / (F.k_B * T)) diff --git a/test/system_tests/test_advection_term_integration.py b/test/system_tests/test_advection_term_integration.py index 2ba14e9de..0eec01693 100644 --- a/test/system_tests/test_advection_term_integration.py +++ b/test/system_tests/test_advection_term_integration.py @@ -15,7 +15,7 @@ def test_MMS_1_species_1_trap_with_advection(): """MMS coupled heat and hydrogen test with 1 mobile species and 1 trap in a 1s transient, the values of the temperature, mobile and trapped solutions at the last - time step is compared to an analytical solution""" + time step is compared to an analytical solution.""" test_mesh_2d = create_unit_square(MPI.COMM_WORLD, 200, 200) x_2d = ufl.SpatialCoordinate(test_mesh_2d) @@ -38,8 +38,10 @@ def test_MMS_1_species_1_trap_with_advection(): V = fem.functionspace(test_mesh_2d, ("Lagrange", 1)) T = fem.Function(V) + def T_expr(x): return 100 + 200 * x[0] + 100 * x[1] + T.interpolate(T_expr) # create velocity field @@ -66,6 +68,7 @@ def velocity_func(x): # define hydrogen problem def exact_mobile_solution(x): return 200 * x[0] ** 2 + 300 * x[1] ** 2 + def exact_trapped_solution(x): return 10 * x[0] ** 2 + 10 * x[1] ** 2 diff --git a/test/system_tests/test_coupled_heat_hydrogen.py b/test/system_tests/test_coupled_heat_hydrogen.py index 34122882b..f9865e692 100644 --- a/test/system_tests/test_coupled_heat_hydrogen.py +++ b/test/system_tests/test_coupled_heat_hydrogen.py @@ -9,7 +9,7 @@ def test_MMS_coupled_problem(): """MMS coupled heat and hydrogen test with 1 mobile species and 1 trap in a 1s transient, the values of the temperature, mobile and trapped solutions at the last - time step is compared to an analytical solution""" + time step is compared to an analytical solution.""" # coupled simulation properties density, heat_capacity = 1.2, 2.6 @@ -70,11 +70,13 @@ def exact_T_solution(x, t): # define hydrogen problem def exact_mobile_solution(x, t): return 2 * x[0] ** 2 + 15 * t + def exact_trapped_solution(x, t): return 4 * x[0] ** 2 + 12 * t def exact_mobile_intial_cond(x): return exact_mobile_solution(x, t=0) + def exact_trapped_intial_cond(x): return exact_trapped_solution(x, t=0) @@ -83,8 +85,10 @@ def exact_trapped_intial_cond(x): def D(x, t): return D_0 * ufl.exp(-E_D / (k_B * exact_T_solution(x, t))) + def k(x, t): return k_0 * ufl.exp(-E_k / (k_B * exact_T_solution(x, t))) + def p(x, t): return p_0 * ufl.exp(-E_p / (k_B * exact_T_solution(x, t))) @@ -186,8 +190,10 @@ def mms_trapped_source(x, t): def exact_final_T(x): return exact_T_solution(x, t=final_time) + def exact_final_mobile(x): return exact_mobile_solution(x, t=final_time) + def exact_final_trapped(x): return exact_trapped_solution(x, t=final_time) @@ -204,7 +210,7 @@ def exact_final_trapped(x): def test_coupled_problem_non_matching_mesh(): """MMS coupled heat and hydrogen test with 1 mobile species in a 1s transient with mismatched meshes, the value of the mobile solution at the last time step is - compared to an analytical solution""" + compared to an analytical solution.""" # coupled simulation properties density, heat_capacity = 1.3, 2.3 diff --git a/test/system_tests/test_cylindrical_coordinates.py b/test/system_tests/test_cylindrical_coordinates.py index 432549f46..0335c8715 100644 --- a/test/system_tests/test_cylindrical_coordinates.py +++ b/test/system_tests/test_cylindrical_coordinates.py @@ -6,10 +6,8 @@ def test_run_MMS_cylindrical(): - """ - Tests that festim produces the correct concentration field in cylindrical - coordinates - """ + """Tests that festim produces the correct concentration field in cylindrical + coordinates.""" my_mesh = F.Mesh1D(vertices=np.linspace(1, 2, 500), coordinate_system="cylindrical") @@ -68,10 +66,8 @@ def u_exact(x): def test_run_MMS_cylindrical_mixed_domain(): - """ - Tests that festim produces the correct concentration field in cylindrical - coordinates in a discontinuous domain with two materials - """ + """Tests that festim produces the correct concentration field in cylindrical + coordinates in a discontinuous domain with two materials.""" my_model = F.HydrogenTransportProblemDiscontinuous() @@ -134,6 +130,7 @@ def lap_c(r): def f_left(x): return -D * lap_c(x[0]) + def f_right(x): return -D * K_S_right / K_S_left * lap_c(x[0]) diff --git a/test/system_tests/test_fluxes.py b/test/system_tests/test_fluxes.py index 5a5db61c3..70304488b 100644 --- a/test/system_tests/test_fluxes.py +++ b/test/system_tests/test_fluxes.py @@ -11,9 +11,7 @@ def test_flux_bc_1_mobile_MMS_steady_state(): - """ - MMS test with a flux BC considering one mobile species at steady state - """ + """MMS test with a flux BC considering one mobile species at steady state.""" def u_exact(x): return 1 + 2 * x[0] ** 2 @@ -61,9 +59,7 @@ def u_exact(x): def test_flux_bc_heat_transfer_steady_state(): - """ - MMS test with a flux BC in a heat transfer problem at steady state - """ + """MMS test with a flux BC in a heat transfer problem at steady state.""" def u_exact(x): return 1 + 2 * x[0] ** 2 diff --git a/test/system_tests/test_heat_transfer.py b/test/system_tests/test_heat_transfer.py index ded5284f1..96e41fcd6 100644 --- a/test/system_tests/test_heat_transfer.py +++ b/test/system_tests/test_heat_transfer.py @@ -11,9 +11,7 @@ def test_heat_transfer_steady_state(): - """ - MMS test with heat transfer simulation - """ + """MMS test with heat transfer simulation.""" def u_exact(mod): return lambda x: 1 + mod.sin(2 * mod.pi * x[0]) diff --git a/test/system_tests/test_io.py b/test/system_tests/test_io.py index f0eceb58a..4039431cc 100644 --- a/test/system_tests/test_io.py +++ b/test/system_tests/test_io.py @@ -9,9 +9,7 @@ def test_writing_and_reading_of_species_function_using_checkpoints(tmpdir): - """ - Tests that a model can write a checkpoint file and another model can read it. - """ + """Tests that a model can write a checkpoint file and another model can read it.""" mesh = dolfinx.mesh.create_unit_square( MPI.COMM_WORLD, nx=10, ny=10, cell_type=dolfinx.cpp.mesh.CellType.quadrilateral ) @@ -117,8 +115,8 @@ def test_writing_and_reading_of_species_function_using_checkpoints(tmpdir): def test_VTXExport_times_added_to_milestones(tmpdir): """Creates a HydrogenTransportProblem object and checks that, if no - stepsize.milestones are given and VTXExport.times are given, VTXExport.times are - are added to stepsize.milestones by .initialise() + stepsize.milestones are given and VTXExport.times are given, VTXExport.times are are + added to stepsize.milestones by .initialise() Args: tmpdir (os.PathLike): path to the pytest temporary folder @@ -156,8 +154,8 @@ def test_VTXExport_times_added_to_milestones(tmpdir): def test_vtx_writer_called_only_at_specified_times(tmpdir): - """test that the VTXWriter.write function is called the number of times specified in - the export.times""" + """Test that the VTXWriter.write function is called the number of times specified in + the export.times.""" filename = str(tmpdir.join("my_export.bp")) diff --git a/test/system_tests/test_misc.py b/test/system_tests/test_misc.py index 47d855942..f565b2d9a 100644 --- a/test/system_tests/test_misc.py +++ b/test/system_tests/test_misc.py @@ -44,10 +44,11 @@ def test_petsc_options(): def test_D_global_on_2d_mesh(): - """ - Test that the D_global is defined correctly on a 2D mesh with two different - materials. The D_global should be defined as a piecewise constant function - with two different values, one for each material. + """Test that the D_global is defined correctly on a 2D mesh with two different + materials. + + The D_global should be defined as a piecewise constant function with two different + values, one for each material. """ mesh, mt, ct = generate_mesh(20) @@ -98,7 +99,7 @@ def test_D_global_on_2d_mesh(): ], ) def test_min_max_vol_on_2d_mesh(species): - """Added test that catches bug #908""" + """Added test that catches bug #908.""" mesh, mt, ct = generate_mesh(10) my_model = F.HydrogenTransportProblem() @@ -183,7 +184,7 @@ def test_min_max_vol_on_2d_mesh(species): def test_temp_dependent_bc_mixed_domain_temperature_as_function(): - """Test to catch bug 986""" + """Test to catch bug 986.""" mesh, mt, ct = generate_mesh(8) V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1)) @@ -248,8 +249,8 @@ def test_del(): def test_multispecies_with_immobile(): - """test to catch bug 1036, that the diffusion coefficent for a non-mobile species - is not requested in the dict given to material""" + """Test to catch bug 1036, that the diffusion coefficent for a non-mobile species is + not requested in the dict given to material.""" my_model = F.HydrogenTransportProblem() my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, num=50)) @@ -287,9 +288,7 @@ def test_multispecies_with_immobile(): def test_implicit_species_bug_reaction(): - """ - This test catches the bug described in issue #1084 - """ + """This test catches the bug described in issue #1084.""" tungsten = F.Material( D_0=1, E_D=0, @@ -336,7 +335,7 @@ def test_implicit_species_bug_reaction(): def test_timesteps(): - """Test that the timesteps are correctly saved in the model""" + """Test that the timesteps are correctly saved in the model.""" my_model = F.HydrogenTransportProblem() my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, num=50)) diff --git a/test/system_tests/test_multi_mat_penalty.py b/test/system_tests/test_multi_mat_penalty.py index 37c5aed23..db693e58c 100644 --- a/test/system_tests/test_multi_mat_penalty.py +++ b/test/system_tests/test_multi_mat_penalty.py @@ -70,14 +70,12 @@ def test_2_materials_2d_mms(tmpdir): K_S_bot = 6.0 D_top = 2.0 D_bot = 5.0 + def c_exact_top_ufl(x): - return ( - 3 + ufl.sin(ufl.pi * (2 * x[1] + 0.5)) + 0.1 * ufl.cos(2 * ufl.pi * x[0]) - ) + return 3 + ufl.sin(ufl.pi * (2 * x[1] + 0.5)) + 0.1 * ufl.cos(2 * ufl.pi * x[0]) + def c_exact_top_np(x): - return ( - 3 + np.sin(np.pi * (2 * x[1] + 0.5)) + 0.1 * np.cos(2 * np.pi * x[0]) - ) + return 3 + np.sin(np.pi * (2 * x[1] + 0.5)) + 0.1 * np.cos(2 * np.pi * x[0]) def c_exact_bot_ufl(x): return K_S_bot / K_S_top**2 * c_exact_top_ufl(x) ** 2 diff --git a/test/system_tests/test_multi_material.py b/test/system_tests/test_multi_material.py index d21d2688f..700c8cebf 100644 --- a/test/system_tests/test_multi_material.py +++ b/test/system_tests/test_multi_material.py @@ -71,18 +71,15 @@ def test_2_materials_2d_mms(tmpdir): K_S_bot = 6.0 D_top = 2.0 D_bot = 5.0 + def c_exact_top_ufl(x): - return ( - 1 + ufl.sin(ufl.pi * (2 * x[0] + 0.5)) + ufl.cos(2 * ufl.pi * x[1]) - ) + return 1 + ufl.sin(ufl.pi * (2 * x[0] + 0.5)) + ufl.cos(2 * ufl.pi * x[1]) def c_exact_bot_ufl(x): return K_S_bot / K_S_top * c_exact_top_ufl(x) def c_exact_top_np(x): - return ( - 1 + np.sin(np.pi * (2 * x[0] + 0.5)) + np.cos(2 * np.pi * x[1]) - ) + return 1 + np.sin(np.pi * (2 * x[0] + 0.5)) + np.cos(2 * np.pi * x[1]) def c_exact_bot_np(x): return K_S_bot / K_S_top * c_exact_top_np(x) @@ -129,13 +126,13 @@ def c_exact_bot_np(x): ] def source_top_val(x): - return ( - 8 * ufl.pi**2 * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1])) - ) + return 8 * ufl.pi**2 * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1])) + def source_bottom_val(x): return ( 40 * ufl.pi**2 * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1])) ) + my_model.sources = [ F.ParticleSource(volume=top_domain, species=H, value=source_top_val), F.ParticleSource(volume=bottom_domain, species=H, value=source_bottom_val), @@ -267,8 +264,8 @@ def test_3_materials_transient(tmpdir): left_surface = F.SurfaceSubdomain1D(id=1, x=vertices[0]) right_surface = F.SurfaceSubdomain1D(id=2, x=vertices[-1]) - # the ids here are arbitrary in 1D, you can put anything as long as it's not the same as the surfaces - # TODO remove mesh and meshtags from these arguments + # the ids here are arbitrary in 1D, you can put anything as long as it's not the + # same as the surfaces my_model.interfaces = [ F.Interface(6, (left_domain, middle_domain)), F.Interface(7, (middle_domain, right_domain)), @@ -403,8 +400,9 @@ def test_2_mats_particle_flux_bc(tmpdir): def test_all_cells_are_not_tagged(): - """ - Checks that an error is raised when not all cells are tagged with a non-zero value. + """Checks that an error is raised when not all cells are tagged with a non-zero + value. + This can be caused by a volume subdomain not being defined correctly, or by a mesh that is too coarse to capture the geometry of the volume subdomains. """ diff --git a/test/system_tests/test_multi_material_change_of_var.py b/test/system_tests/test_multi_material_change_of_var.py index 17d64b6e2..14e68af7b 100644 --- a/test/system_tests/test_multi_material_change_of_var.py +++ b/test/system_tests/test_multi_material_change_of_var.py @@ -155,18 +155,15 @@ def test_2_materials_2d_mms(tmpdir): K_S_bot = 6.0 D_top = 2.0 D_bot = 5.0 + def c_exact_top_ufl(x): - return ( - 1 + ufl.sin(ufl.pi * (2 * x[0] + 0.5)) + ufl.cos(2 * ufl.pi * x[1]) - ) + return 1 + ufl.sin(ufl.pi * (2 * x[0] + 0.5)) + ufl.cos(2 * ufl.pi * x[1]) def c_exact_bot_ufl(x): return K_S_bot / K_S_top * c_exact_top_ufl(x) def c_exact_top_np(x): - return ( - 1 + np.sin(np.pi * (2 * x[0] + 0.5)) + np.cos(2 * np.pi * x[1]) - ) + return 1 + np.sin(np.pi * (2 * x[0] + 0.5)) + np.cos(2 * np.pi * x[1]) def c_exact_bot_np(x): return K_S_bot / K_S_top * c_exact_top_np(x) @@ -203,13 +200,13 @@ def c_exact_bot_np(x): ] def source_top_val(x): - return ( - 8 * ufl.pi**2 * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1])) - ) + return 8 * ufl.pi**2 * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1])) + def source_bottom_val(x): return ( 40 * ufl.pi**2 * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1])) ) + my_model.sources = [ F.ParticleSource(volume=top_domain, species=H, value=source_top_val), F.ParticleSource(volume=bottom_domain, species=H, value=source_bottom_val), diff --git a/test/system_tests/test_non_homogeneous_density.py b/test/system_tests/test_non_homogeneous_density.py index 73d55a8d1..5b2fa4176 100644 --- a/test/system_tests/test_non_homogeneous_density.py +++ b/test/system_tests/test_non_homogeneous_density.py @@ -9,12 +9,16 @@ def step_function_space(x): return ufl.conditional(ufl.gt(x[0], 0.5), 10, 0.0) + + def step_function_time_non_homogeneous(x, t): - return ufl.conditional( - ufl.gt(t, 50), 10 - 10 * x[0], 0.0 -) + return ufl.conditional(ufl.gt(t, 50), 10 - 10 * x[0], 0.0) + + def step_function_time_homogeneous(t): return 10.0 if t > 50 else 2.0 + + def simple_space(x): return 10 - 10 * x[0] diff --git a/test/system_tests/test_reaction_not_in_every_volume.py b/test/system_tests/test_reaction_not_in_every_volume.py index b1cdea67b..5b9fd8e22 100644 --- a/test/system_tests/test_reaction_not_in_every_volume.py +++ b/test/system_tests/test_reaction_not_in_every_volume.py @@ -4,8 +4,8 @@ def test_sim_reaction_not_in_every_volume(): - """Tests that a steady simulation can be run if a reaction is not defined - in every volume""" + """Tests that a steady simulation can be run if a reaction is not defined in every + volume.""" my_model = F.HydrogenTransportProblem() my_model.mesh = F.Mesh1D(np.linspace(0, 2e-04, num=1001)) diff --git a/test/system_tests/test_spherical_coordinates.py b/test/system_tests/test_spherical_coordinates.py index 956fe930d..5a59b493b 100644 --- a/test/system_tests/test_spherical_coordinates.py +++ b/test/system_tests/test_spherical_coordinates.py @@ -7,10 +7,8 @@ def test_run_MMS_spherical(): - """ - Tests that festim produces the correct concentration field in spherical - coordinates - """ + """Tests that festim produces the correct concentration field in spherical + coordinates.""" my_mesh = F.Mesh1D(vertices=np.linspace(1, 2, 1000), coordinate_system="spherical") fem.functionspace(my_mesh.mesh, ("Lagrange", 1)) @@ -69,10 +67,8 @@ def u_exact(x): def test_run_MMS_spherical_mixed_domain(): - """ - Tests that festim produces the correct concentration field in spherical - coordinates in a discontinuous domain with two materials - """ + """Tests that festim produces the correct concentration field in spherical + coordinates in a discontinuous domain with two materials.""" my_model = F.HydrogenTransportProblemDiscontinuous() @@ -135,6 +131,7 @@ def lap_c(r): def f_left(x): return -D * lap_c(x[0]) + def f_right(x): return -D * K_S_right / K_S_left * lap_c(x[0]) diff --git a/test/system_tests/test_surface_reactions.py b/test/system_tests/test_surface_reactions.py index 6754c10f1..7efaa2bad 100644 --- a/test/system_tests/test_surface_reactions.py +++ b/test/system_tests/test_surface_reactions.py @@ -22,13 +22,12 @@ def compute(self, u, ds, entity_maps=None): def test_2_isotopes_no_pressure(): - """ - Runs a simple 1D hydrogen transport problem with 2 isotopes - and 3 surface reactions on the right boundary. + """Runs a simple 1D hydrogen transport problem with 2 isotopes and 3 surface + reactions on the right boundary. - Then checks that the fluxes of the isotopes are consistent with the surface reactions - by computing the flux of H and D (from the gradient) and comparing it to the fluxes - computed from the surface reactions. + Then checks that the fluxes of the isotopes are consistent with the + surface reactions by computing the flux of H and D (from the gradient) + and comparing it to the fluxes computed from the surface reactions. Example: H + D <-> HD, H + H <-> HH, D + D <-> DD -D grad(c_H) n = 2*kr*c_H*c_H + kr*c_H*c_D @@ -149,13 +148,12 @@ def test_2_isotopes_no_pressure(): def test_2_isotopes_with_pressure(): - """ - Runs a simple 1D hydrogen transport problem with 2 isotopes - and 3 surface reactions on the right boundary. + """Runs a simple 1D hydrogen transport problem with 2 isotopes and 3 surface + reactions on the right boundary. - Then checks that the fluxes of the isotopes are consistent with the surface reactions - by computing the flux of H and D (from the gradient) and comparing it to the fluxes - computed from the surface reactions. + Then checks that the fluxes of the isotopes are consistent with the + surface reactions by computing the flux of H and D (from the gradient) and + comparing it to the fluxes computed from the surface reactions. Example: H + D <-> HD, H + H <-> HH, D + D <-> DD -D grad(c_H) n = 2*kr*c_H*c_H + kr*c_H*c_D - kd*P_H2 @@ -276,9 +274,8 @@ def test_2_isotopes_with_pressure(): def test_pressure_varies_in_time(): - """ - Runs a problem with a surface reaction and a time-dependent pressure - on the right boundary. + """Runs a problem with a surface reaction and a time-dependent pressure on the right + boundary. Then checks that the flux is consistent with the surface reaction """ diff --git a/test/test_advection.py b/test/test_advection.py index 88ed2f857..17c00e2fc 100644 --- a/test/test_advection.py +++ b/test/test_advection.py @@ -91,7 +91,7 @@ def test_subdomain_accepts_None_value(): def test_velocity_field_update(): - """Test a explicit time dependent input value is updated correctly""" + """Test a explicit time dependent input value is updated correctly.""" test_mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10) V = fem.functionspace(test_mesh, ("Lagrange", 1)) diff --git a/test/test_average_surface.py b/test/test_average_surface.py index b8c4950ea..74f875736 100644 --- a/test/test_average_surface.py +++ b/test/test_average_surface.py @@ -6,7 +6,7 @@ def test_average_surface_compute_1D(): - """Test that the average surface export computes the correct value""" + """Test that the average surface export computes the correct value.""" # BUILD L = 6.0 diff --git a/test/test_average_volume.py b/test/test_average_volume.py index 8f71121b9..de1fddd8b 100644 --- a/test/test_average_volume.py +++ b/test/test_average_volume.py @@ -8,7 +8,7 @@ def test_average_volume_compute_1D(): - """Test that the average volume export computes the correct value""" + """Test that the average volume export computes the correct value.""" # BUILD L = 6.0 diff --git a/test/test_complex_reaction.py b/test/test_complex_reaction.py index bcbbb7923..170e68c98 100644 --- a/test/test_complex_reaction.py +++ b/test/test_complex_reaction.py @@ -5,10 +5,10 @@ def concentration_A_exact(t, c_A_0, k, p): - """Analytical solution for the concentration of species A in a reaction A + B <-> C + D - assuming [A]_0 = [B]_0 and [C]_0 = [D]_0 = 0 - can be obtained by solving - d[A]/dt = -k[A][B] + p[C][D] + """Analytical solution for the concentration of species A in a reaction A + B <-> C + + D assuming [A]_0 = [B]_0 and [C]_0 = [D]_0 = 0 can be obtained by solving d[A]/dt + = -k[A][B] + p[C][D] + where [A] = [B] and [C] = [D] = 0 at t=0 Moreover, [C] = [D] = [A]_0 - [A] The equation then becomes @@ -39,8 +39,8 @@ def concentration_A_exact(t, c_A_0, k, p): def model_test_reaction(stepsize=1, k=350e-4, p=120e-4, c_A_0=1): - """Creates a festim model with a single reaction and runs it. - The reaction is A + B <-> C + D + """Creates a festim model with a single reaction and runs it. The reaction is A + B + <-> C + D. Args: stepsize (float): the stepsize @@ -142,7 +142,8 @@ def model_test_reaction(stepsize=1, k=350e-4, p=120e-4, c_A_0=1): def compute_error(model): - """Computes the relative L2 error norm between the analytical solution and the numerical solution + """Computes the relative L2 error norm between the analytical solution and the + numerical solution. Args: model (F.HydrogenTransportModel): the numerical festim model @@ -170,14 +171,15 @@ def compute_error(model): @pytest.mark.parametrize("k, p, c_A_0", [(350e-4, 120e-4, 3), (200e-4, 100e-4, 2)]) def test_reaction(k, p, c_A_0): - """Test the reaction A + B <-> C + D with a festim model and compare the results with the analytical solution""" + """Test the reaction A + B <-> C + D with a festim model and compare the results + with the analytical solution.""" model = model_test_reaction(stepsize=0.5, k=k, p=p, c_A_0=c_A_0) error = compute_error(model) assert error < 1e-2 def test_surface_reaction_in_discontinuous_case(): - """Test the SurfaceReactionBC in a discontinuous problem setup""" + """Test the SurfaceReactionBC in a discontinuous problem setup.""" my_model = F.HydrogenTransportProblemDiscontinuous() my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, num=1000)) diff --git a/test/test_coupled_heat_hydrogen_problem.py b/test/test_coupled_heat_hydrogen_problem.py index c905d74ac..b7e51f0fb 100644 --- a/test/test_coupled_heat_hydrogen_problem.py +++ b/test/test_coupled_heat_hydrogen_problem.py @@ -15,13 +15,13 @@ ) def test_error_raised_when_wrong_hydrogen_problem_given(object): """Test TypeError is raised when an object that isnt a - festim.HydrogenTransportProblem object is given to hydrogen_problem""" + festim.HydrogenTransportProblem object is given to hydrogen_problem.""" test_heat_problem = F.HeatTransferProblem() with pytest.raises( TypeError, - match="hydrogen_problem must be a festim.HydrogenTransportProblem object", + match=r"hydrogen_problem must be a festim.HydrogenTransportProblem object", ): F.CoupledTransientHeatTransferHydrogenTransport( heat_problem=test_heat_problem, hydrogen_problem=object @@ -37,7 +37,7 @@ def test_error_raised_when_wrong_hydrogen_problem_given(object): ) def test_error_raised_when_wrong_type_hydrogen_problem_given(object): """Test TypeError is raised when an object that isnt a - festim.HydrogenTransportProblem object is given to hydrogen_problem""" + festim.HydrogenTransportProblem object is given to hydrogen_problem.""" test_heat_problem = F.HeatTransferProblem() @@ -58,14 +58,14 @@ def test_error_raised_when_wrong_type_hydrogen_problem_given(object): ["coucou", 1, 1.0], ) def test_error_raised_when_wrong_heat_problem_given(object): - """Test TypeError is raised when an object that isnt a - festim.HeatTransferProblem object is given to heat_problem""" + """Test TypeError is raised when an object that isnt a festim.HeatTransferProblem + object is given to heat_problem.""" test_hydrogen_problem = F.HydrogenTransportProblem() with pytest.raises( TypeError, - match="heat_problem must be a festim.HeatTransferProblem object", + match=r"heat_problem must be a festim.HeatTransferProblem object", ): F.CoupledTransientHeatTransferHydrogenTransport( heat_problem=object, hydrogen_problem=test_hydrogen_problem @@ -74,7 +74,7 @@ def test_error_raised_when_wrong_heat_problem_given(object): def test_initial_dt_values_are_the_same(): """Test that the smallest of the stepsize intial_value values given is used in both - the heat_problem and the hydrogen_problem""" + the heat_problem and the hydrogen_problem.""" test_heat_problem = F.HeatTransferProblem( mesh=test_mesh, @@ -107,7 +107,7 @@ def test_initial_dt_values_are_the_same(): def test_dts_always_the_same(): - """Test that the correct dt value is modified""" + """Test that the correct dt value is modified.""" dt = F.Stepsize( 1000, growth_factor=1000, cutback_factor=0.9, target_nb_iterations=4 @@ -150,7 +150,7 @@ def test_dts_always_the_same(): def test_error_raised_when_final_times_not_the_same(): """Test that an error is raised when the final time values given in the heat_problem - and the hydrogen_problem are not the same""" + and the hydrogen_problem are not the same.""" test_heat_problem = F.HeatTransferProblem( mesh=test_mesh, @@ -184,7 +184,7 @@ def test_error_raised_when_final_times_not_the_same(): def test_error_raised_when_both_problems_not_transient(): """Test that an error is raised when the heat_problem and hydorgen_problem are not - set to transient""" + set to transient.""" test_heat_problem = F.HeatTransferProblem( settings=F.Settings(atol=1, rtol=0.9999, transient=False), diff --git a/test/test_dirichlet_bc.py b/test/test_dirichlet_bc.py index 0dd366644..bb632552c 100644 --- a/test/test_dirichlet_bc.py +++ b/test/test_dirichlet_bc.py @@ -15,7 +15,7 @@ def test_init(): - """Test that the attributes are set correctly""" + """Test that the attributes are set correctly.""" # create a DirichletBC object subdomain = F.SurfaceSubdomain1D(1, x=0) value = 1.0 @@ -31,9 +31,8 @@ def test_init(): def test_value_fenics(): - """Test that the value_fenics attribute can be set to a valid value - and that an invalid type throws an error - """ + """Test that the value_fenics attribute can be set to a valid value and that an + invalid type throws an error.""" # create a DirichletBC object subdomain = F.SurfaceSubdomain1D(1, x=0) value = 1.0 @@ -53,7 +52,7 @@ def test_value_fenics(): def test_callable_for_value(): - """Test that the value attribute can be a callable function of x and t""" + """Test that the value attribute can be a callable function of x and t.""" subdomain = F.SurfaceSubdomain1D(1, x=1) vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat) @@ -97,7 +96,7 @@ def value(x, t): def test_value_callable_x_t_T(): - """Test that the value attribute can be a callable function of x, t, and T""" + """Test that the value attribute can be a callable function of x, t, and T.""" subdomain = F.SurfaceSubdomain1D(1, x=1) vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat) @@ -145,7 +144,7 @@ def value(x, t, T): @pytest.mark.parametrize("value", [lambda t: t, lambda t: 1.0 + t]) def test_callable_t_only(value): - """Test that the value attribute can be a callable function of t only""" + """Test that the value attribute can be a callable function of t only.""" subdomain = F.SurfaceSubdomain1D(1, x=1) vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat) @@ -186,7 +185,7 @@ def test_callable_t_only(value): def test_callable_x_only(): - """Test that the value attribute can be a callable function of x only""" + """Test that the value attribute can be a callable function of x only.""" # BUILD subdomain = F.SurfaceSubdomain1D(1, x=1) @@ -247,8 +246,8 @@ def value(x): ], ) def test_integration_with_HTransportProblem(value): - """test that different callable functions can be applied to a dirichlet - boundary condition, asserting in each case they match an expected value""" + """Test that different callable functions can be applied to a dirichlet boundary + condition, asserting in each case they match an expected value.""" subdomain = F.SurfaceSubdomain1D(1, x=1) vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat) @@ -311,8 +310,8 @@ def test_integration_with_HTransportProblem(value): ], ) def test_define_value_error_if_ufl_conditional_t_only(value): - """Test that a ValueError is raised when the value attribute is a callable - of t only and contains a ufl conditional""" + """Test that a ValueError is raised when the value attribute is a callable of t only + and contains a ufl conditional.""" subdomain = F.SurfaceSubdomain1D(1, x=1) species = F.Species("test") @@ -322,14 +321,14 @@ def test_define_value_error_if_ufl_conditional_t_only(value): t = fem.Constant(mesh, 0.0) V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1)) with pytest.raises( - ValueError, match="self.value should return a float or an int, not " + ValueError, match=r"self.value should return a float or an int, not " ): bc.create_value(V, temperature=None, t=t) def test_species_predefined(): - """Test a ValueError is raised when the species defined in the boundary - condition is not predefined in the model""" + """Test a ValueError is raised when the species defined in the boundary condition is + not predefined in the model.""" subdomain = F.SurfaceSubdomain1D(1, x=1) vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat) @@ -362,9 +361,9 @@ def test_species_predefined(): ], ) def test_integration_with_a_multispecies_HTransportProblem(value_A, value_B): - """test that a mixture of callable functions can be applied to dirichlet - boundary conditions in a multispecies case, asserting in each case they - match an expected value""" + """Test that a mixture of callable functions can be applied to dirichlet boundary + conditions in a multispecies case, asserting in each case they match an expected + value.""" subdomain_A = F.SurfaceSubdomain1D(1, x=0) subdomain_B = F.SurfaceSubdomain1D(2, x=1) vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat) @@ -438,7 +437,7 @@ def test_integration_with_a_multispecies_HTransportProblem(value_A, value_B): ], ) def test_bc_time_dependent_attribute(input, expected_value): - """Test that the time_dependent attribute is correctly set""" + """Test that the time_dependent attribute is correctly set.""" subdomain = F.SurfaceSubdomain1D(1, x=0) species = F.Species("test") my_bc = F.DirichletBC(subdomain, input, species) @@ -463,7 +462,7 @@ def test_bc_time_dependent_attribute(input, expected_value): ], ) def test_bc_temperature_dependent_attribute(input, expected_value): - """Test that the temperature_dependent attribute is correctly set""" + """Test that the temperature_dependent attribute is correctly set.""" subdomain = F.SurfaceSubdomain1D(1, x=0) species = F.Species("test") my_bc = F.DirichletBC(subdomain, input, species) diff --git a/test/test_flux_bc.py b/test/test_flux_bc.py index 04ff92183..d048d3277 100644 --- a/test/test_flux_bc.py +++ b/test/test_flux_bc.py @@ -14,7 +14,7 @@ def test_init(): - """Test that the attributes are set correctly""" + """Test that the attributes are set correctly.""" # create a DirichletBC object subdomain = F.SurfaceSubdomain1D(1, x=0) value = 1.0 @@ -46,7 +46,7 @@ def test_init(): ], ) def test_create_value_fenics_type(value, expected_type): - """Test that""" + """Test that.""" # BUILD left = F.SurfaceSubdomain1D(1, x=0) my_species = F.Species("test") @@ -72,7 +72,7 @@ def test_create_value_fenics_type(value, expected_type): ], ) def test_create_value_fenics_value(value, expected_value): - """Test that""" + """Test that.""" # BUILD left = F.SurfaceSubdomain1D(1, x=0) my_species = F.Species("test") @@ -90,7 +90,8 @@ def test_create_value_fenics_value(value, expected_value): def test_create_value_fenics_dependent_conc(): - """Test that the value_fenics of ParticleFluxBC is set correctly when the value is dependent on the concentration""" + """Test that the value_fenics of ParticleFluxBC is set correctly when the value is + dependent on the concentration.""" # BUILD left = F.SurfaceSubdomain1D(1, x=0) my_species = F.Species("test") @@ -119,14 +120,17 @@ def test_value_fenics_setter_error(): with pytest.raises( TypeError, - match="Value must be a dolfinx.fem.Function, dolfinx.fem.Constant, np.ndarray or ufl.core.expr.Expr not ", + match=( + r"Value must be a dolfinx.fem.Function, dolfinx.fem.Constant, np.ndarray or " # noqa: E501 + r"ufl.core.expr.Expr not " + ), ): bc.value_fenics = "coucou" def test_ValueError_raised_when_callable_returns_wrong_type(): """The create_value_fenics method should raise a ValueError when the callable - returns an object which is not a float or int""" + returns an object which is not a float or int.""" surface = F.SurfaceSubdomain(id=1) species = F.Species("test") @@ -141,7 +145,10 @@ def my_value(t): with pytest.raises( ValueError, - match="self.value should return a float or an int, not ", ): my_model.species = [1, 2, 3] def test_create_initial_conditions_ValueError_raised_when_not_transient(): - """Test that ValueError is raised if initial conditions are defined in - a steady state simulation""" + """Test that ValueError is raised if initial conditions are defined in a steady + state simulation.""" my_vol = F.VolumeSubdomain1D(id=1, borders=[0, 4], material=dummy_mat) H = F.Species("H") @@ -852,8 +853,8 @@ def test_create_initial_conditions_ValueError_raised_when_not_transient(): ], ) def test_create_initial_conditions_expr_fenics(input_value, expected_value): - """Test that after calling create_initial_conditions, the prev_solution - attribute of the species has the correct value at x=4.0.""" + """Test that after calling create_initial_conditions, the prev_solution attribute of + the species has the correct value at x=4.0.""" # BUILD vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 4], material=dummy_mat) @@ -931,8 +932,8 @@ def test_create_species_from_trap(): def test_create_initial_conditions_value_fenics_multispecies( input_value_1, input_value_2, expected_value_1, expected_value_2 ): - """Test that after calling create_initial_conditions, the prev_solution - attribute of each species has the correct value at x=4.0 in a multispecies case""" + """Test that after calling create_initial_conditions, the prev_solution attribute of + each species has the correct value at x=4.0 in a multispecies case.""" # BUILD test_mesh = F.Mesh1D(vertices=np.linspace(0, 4, num=101)) @@ -970,7 +971,7 @@ def test_create_initial_conditions_value_fenics_multispecies( def test_adaptive_timestepping_grows(): - """Tests that the stepsize grows""" + """Tests that the stepsize grows.""" # BUILD my_vol = F.VolumeSubdomain1D(id=1, borders=[0, 4], material=dummy_mat) my_model = F.HydrogenTransportProblem( @@ -1006,7 +1007,7 @@ def test_adaptive_timestepping_grows(): def test_adaptive_timestepping_shrinks(): - """Tests that the stepsize shrinks""" + """Tests that the stepsize shrinks.""" # BUILD my_vol = F.VolumeSubdomain1D(id=1, borders=[0, 4], material=dummy_mat) my_model = F.HydrogenTransportProblem( @@ -1053,8 +1054,8 @@ def test_adaptive_timestepping_shrinks(): ], ) def test_reinstantiation_of_class(attribute, value): - """Test that when an attribute defaults to empty list, when the class - is reinstantiated the list is not passed to the new class object""" + """Test that when an attribute defaults to empty list, when the class is + reinstantiated the list is not passed to the new class object.""" model_1 = F.HydrogenTransportProblem() getattr(model_1, attribute).append(value) @@ -1064,8 +1065,8 @@ def test_reinstantiation_of_class(attribute, value): def test_define_meshtags_and_measures_with_custom_fenics_mesh(): - """Test that the define_meshtags_and_measures method works when the mesh is - a custom fenics mesh""" + """Test that the define_meshtags_and_measures method works when the mesh is a custom + fenics mesh.""" # BUILD mesh_1D = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, 10) @@ -1096,7 +1097,7 @@ def test_define_meshtags_and_measures_with_custom_fenics_mesh(): def test_error_raised_when_custom_fenics_mesh_wrong_facet_meshtags_type(): - """Test the facet_meshtags type hinting raises error when given as wrong type""" + """Test the facet_meshtags type hinting raises error when given as wrong type.""" # BUILD mesh_1D = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, 10) @@ -1104,12 +1105,12 @@ def test_error_raised_when_custom_fenics_mesh_wrong_facet_meshtags_type(): my_model = F.HydrogenTransportProblem(mesh=my_mesh) # TEST - with pytest.raises(TypeError, match="value must be of type dolfinx.mesh.MeshTags"): + with pytest.raises(TypeError, match=r"value must be of type dolfinx.mesh.MeshTags"): my_model.facet_meshtags = [0, 1] def test_error_raised_when_custom_fenics_mesh_wrong_volume_meshtags_type(): - """Test the volume_meshtags type hinting raises error when given as wrong type""" + """Test the volume_meshtags type hinting raises error when given as wrong type.""" # BUILD mesh_1D = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, 10) @@ -1117,7 +1118,7 @@ def test_error_raised_when_custom_fenics_mesh_wrong_volume_meshtags_type(): my_model = F.HydrogenTransportProblem(mesh=my_mesh) # TEST - with pytest.raises(TypeError, match="value must be of type dolfinx.mesh.MeshTags"): + with pytest.raises(TypeError, match=r"value must be of type dolfinx.mesh.MeshTags"): my_model.volume_meshtags = [0, 1] @@ -1135,8 +1136,8 @@ def test_error_raised_when_custom_fenics_mesh_wrong_volume_meshtags_type(): ], ) def test_update_time_dependent_values_flux(bc_value, expected_values): - """Test that time dependent fluxes are updated at each time step, - and match an expected value""" + """Test that time dependent fluxes are updated at each time step, and match an + expected value.""" # BUILD my_vol = F.VolumeSubdomain1D(id=1, borders=[0, 4], material=dummy_mat) surface = F.SurfaceSubdomain1D(id=2, x=0) @@ -1187,8 +1188,8 @@ def test_update_time_dependent_values_flux(bc_value, expected_values): def test_update_fluxes_with_time_dependent_temperature( temperature_value, bc_value, expected_values ): - """Test that temperature dependent flux terms are updated at each time step - when the temperature is time dependent, and match an expected value""" + """Test that temperature dependent flux terms are updated at each time step when the + temperature is time dependent, and match an expected value.""" # BUILD H = F.Species("H") @@ -1226,7 +1227,7 @@ def test_update_fluxes_with_time_dependent_temperature( def test_create_flux_values_fenics_multispecies(): """Test that the create_flux_values_fenics method correctly sets the value_fenics - attribute in a multispecies case""" + attribute in a multispecies case.""" # BUILD my_vol = F.VolumeSubdomain1D(id=1, borders=[0, 4], material=dummy_mat) surface = F.SurfaceSubdomain1D(id=2, x=0) @@ -1258,7 +1259,7 @@ def test_create_flux_values_fenics_multispecies(): def test_not_implemented_error_raised_with_D_as_function(): """Test that NotImplementedError is raised when D is a function and more than one - volume subdomain has been defined""" + volume subdomain has been defined.""" # BUILD test_mesh = F.Mesh1D(vertices=np.linspace(0, 1, num=101)) @@ -1293,7 +1294,7 @@ def test_not_implemented_error_raised_with_D_as_function(): def test_define_D_global_with_D_as_function(): """Test that if a function is given as diffusion coeff of a material, that it is - passed to D_global""" + passed to D_global.""" # BUILD test_mesh = F.Mesh1D(vertices=np.linspace(0, 1, num=101)) @@ -1323,7 +1324,7 @@ def test_define_D_global_with_D_as_function(): @pytest.mark.parametrize("coord_sys", ["cylindrical", "spherical"]) def test_exports_cyl_sph_coord_system_raise_not_implemented(coord_sys): """Test that NotImplementedError is raised when trying to use exports in a - cylindrical or spherical coordinate system""" + cylindrical or spherical coordinate system.""" # BUILD my_mesh = F.Mesh1D(np.linspace(0, 1, num=11), coordinate_system=coord_sys) @@ -1352,7 +1353,7 @@ def test_exports_cyl_sph_coord_system_raise_not_implemented(coord_sys): def test_traps_with_CG_elements(): """Test that when creating species from traps, if the element for traps is set to - CG, the created species has a CG element as well""" + CG, the created species has a CG element as well.""" # BUILD my_mesh = F.Mesh1D(np.linspace(0, 1, num=11)) @@ -1376,7 +1377,7 @@ def test_traps_with_CG_elements(): def test_surface_reaction_BC_discontinuous(): - """Test that surface reaction BC conserves flux""" + """Test that surface reaction BC conserves flux.""" # BUILD my_model = F.HydrogenTransportProblemDiscontinuous() @@ -1461,7 +1462,7 @@ def test_surface_reaction_BC_discontinuous(): def test_temperature_as_function_in_discontinuous(): """Test that a discontinuous problem can be initialised with a temperature field - given as a function""" + given as a function.""" # BUILD my_model = F.HydrogenTransportProblemDiscontinuous() @@ -1503,7 +1504,7 @@ def test_temperature_as_function_in_discontinuous(): def test_wrong_element_for_immobile_species(): """Test that an error is raised when an immobile species is defined on a subdomain - with an element that is incorrect""" + with an element that is incorrect.""" my_model = F.HydrogenTransportProblem() with pytest.raises(ValueError, match="element_immobile should be in"): @@ -1514,8 +1515,8 @@ def test_wrong_element_for_immobile_species(): "element_type, expected_element", [("CG", "CG"), ("DG", "DG"), ("P", "CG")] ) def test_element_for_immobile_species(element_type, expected_element): - """Test that the element for immobile species is correctly set when the user sets - it to CG or DG""" + """Test that the element for immobile species is correctly set when the user sets it + to CG or DG.""" my_model = F.HydrogenTransportProblem() my_model.element_immobile = element_type @@ -1524,7 +1525,7 @@ def test_element_for_immobile_species(element_type, expected_element): def test_nb_dofs_dg_or_cg(): """Test that the number of dofs is different when using CG or DG elements for - immobile species""" + immobile species.""" # BUILD my_model = F.HydrogenTransportProblem() diff --git a/test/test_heat_transfer_problem.py b/test/test_heat_transfer_problem.py index c4e3dd5f1..9cfcb7e73 100644 --- a/test/test_heat_transfer_problem.py +++ b/test/test_heat_transfer_problem.py @@ -110,7 +110,7 @@ def exact_solution(x): def test_MMS_T_dependent_thermal_cond(): - """MMS test with space T dependent thermal cond""" + """MMS test with space T dependent thermal cond.""" def thermal_conductivity(T): return 3 * T + 2 @@ -159,10 +159,8 @@ def mms_source(x): def test_heat_transfer_transient(tmpdir): - """ - MMS test for transient heat transfer - constant thermal conductivity density and heat capacity - """ + """MMS test for transient heat transfer constant thermal conductivity density and + heat capacity.""" density = 2 heat_capacity = 3 thermal_conductivity = 4 @@ -207,7 +205,8 @@ def exact_solution(x, t): my_problem.settings = F.Settings( atol=1e-8, rtol=1e-10, - final_time=1, # final time shouldn't be too long so that a potential error at the initial timestep is not negligible + # final time shouldn't be too long so that a potential error at the initial timestep is not negligible # noqa: E501 + final_time=1, ) # Forward euler isn't great so dt should be small @@ -222,7 +221,7 @@ def exact_solution(x, t): my_problem.run() computed_solution = my_problem.u - # we use the exact final time of the simulation which may differ from the one specified in the settings + # we use the exact final time of the simulation which may differ from the one specified in the settings # noqa: E501 final_time_sim = my_problem.t.value def exact_solution_end(x): @@ -233,8 +232,8 @@ def exact_solution_end(x): def test_MES(): - """Method of Exact Solution test for transient heat transfer - with thermal cond. k = 2, T = 0 on surfaces, and source term q = 8 k + """Method of Exact Solution test for transient heat transfer with thermal cond. k = + 2, T = 0 on surfaces, and source term q = 8 k. Analytical solution: T(x) = 4 x (1 - x) """ @@ -331,7 +330,8 @@ def mms_source(x, t): my_problem.settings = F.Settings( atol=1e-8, rtol=1e-10, - final_time=1, # final time shouldn't be too long so that a potential error at the initial timestep is not negligible + # final time shouldn't be too long so that a potential error at the initial timestep is not negligible # noqa: E501 + final_time=1, ) # Forward euler isn't great so dt should be small @@ -346,7 +346,7 @@ def mms_source(x, t): my_problem.run() computed_solution = my_problem.u - # we use the exact final time of the simulation which may differ from the one specified in the settings + # we use the exact final time of the simulation which may differ from the one specified in the settings # noqa: E501 final_time_sim = my_problem.t.value def exact_solution_end(x): @@ -357,7 +357,7 @@ def exact_solution_end(x): def test_sources(): - """Tests the sources setter of the HeatTransferProblem class""" + """Tests the sources setter of the HeatTransferProblem class.""" htp = F.HeatTransferProblem() vol = F.VolumeSubdomain1D(1, borders=[0, 1], material=None) # Test that setting valid sources works @@ -369,7 +369,7 @@ def test_sources(): assert htp.sources == valid_sources # Test that setting invalid sources raises a TypeError - with pytest.raises(TypeError, match="festim.HeatSource objects"): + with pytest.raises(TypeError, match=r"festim.HeatSource objects"): spe = F.Species("H") htp.sources = [ F.ParticleSource(1, vol, spe), @@ -378,7 +378,7 @@ def test_sources(): def test_boundary_conditions(): - """Tests the boundary_conditions setter of the HeatTransferProblem class""" + """Tests the boundary_conditions setter of the HeatTransferProblem class.""" htp = F.HeatTransferProblem() left = F.SurfaceSubdomain(1) @@ -399,7 +399,8 @@ def test_boundary_conditions(): @pytest.mark.parametrize("mesh", [mesh_1D, mesh_2D, mesh_3D]) def test_meshtags_from_xdmf(tmp_path, mesh): - """Test that the facet and volume meshtags are read correctly from the mesh XDMF files""" + """Test that the facet and volume meshtags are read correctly from the mesh XDMF + files.""" # create mesh functions fdim = mesh.topology.dim - 1 vdim = mesh.topology.dim @@ -482,7 +483,7 @@ def test_meshtags_from_xdmf(tmp_path, mesh): def test_raise_error_non_unique_vol_ids(): - """Test that an error is raised if the volume ids are not unique""" + """Test that an error is raised if the volume ids are not unique.""" my_problem = F.HeatTransferProblem() my_problem.mesh = F.Mesh1D(vertices=np.linspace(0, 1, 11)) left = F.SurfaceSubdomain1D(id=1, x=0) @@ -503,7 +504,8 @@ def test_raise_error_non_unique_vol_ids(): def test_initial_condition(): - """Test that an error is raised when initial conditions are defined for a transient simulation""" + """Test that an error is raised when initial conditions are defined for a transient + simulation.""" my_problem = F.HeatTransferProblem() my_problem.mesh = F.Mesh1D(vertices=np.linspace(0, 1, 3000)) @@ -536,7 +538,7 @@ def test_initial_condition(): def test_adaptive_timestepping_grows(): - """Tests that the stepsize grows""" + """Tests that the stepsize grows.""" # BUILD my_vol = F.VolumeSubdomain1D(id=1, borders=[0, 4], material=dummy_mat) my_model = F.HeatTransferProblem( @@ -570,7 +572,7 @@ def test_adaptive_timestepping_grows(): def test_adaptive_timestepping_shrinks(): - """Tests that the stepsize shrinks""" + """Tests that the stepsize shrinks.""" # BUILD my_vol = F.VolumeSubdomain1D(id=1, borders=[0, 4], material=dummy_mat) my_model = F.HeatTransferProblem( @@ -616,8 +618,8 @@ def test_adaptive_timestepping_shrinks(): ], ) def test_update_time_dependent_values_HeatFluxBC(bc_value, expected_values): - """Test that time dependent fluxes are updated at each time step, - and match an expected value""" + """Test that time dependent fluxes are updated at each time step, and match an + expected value.""" # BUILD my_vol = F.VolumeSubdomain1D(id=1, borders=[0, 4], material=dummy_mat) surface = F.SurfaceSubdomain1D(id=2, x=0) diff --git a/test/test_helpers.py b/test/test_helpers.py index 38c583fe7..a4fff0129 100644 --- a/test/test_helpers.py +++ b/test/test_helpers.py @@ -26,7 +26,7 @@ ], ) def test_temperature_type_and_processing(value): - """Test that the temperature type is correctly set""" + """Test that the temperature type is correctly set.""" if not isinstance(value, fem.Constant | int | float): with pytest.raises(TypeError): @@ -39,7 +39,7 @@ def test_temperature_type_and_processing(value): "input_value, expected_output_type", [(1.0, fem.Constant), (3, fem.Constant)] ) def test_value_convert_float_int_inputs(input_value, expected_output_type): - """Test that float and value is correctly converted""" + """Test that float and value is correctly converted.""" test_value = F.Value(input_value) @@ -63,7 +63,7 @@ def test_value_convert_float_int_inputs(input_value, expected_output_type): ], ) def test_value_convert_up_to_ufl_inputs(input_value, expected_output_type): - """Test that float and value is correctly converted""" + """Test that float and value is correctly converted.""" my_mesh = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, 10) V = fem.functionspace(my_mesh, ("Lagrange", 1)) @@ -97,7 +97,7 @@ def test_value_convert_up_to_ufl_inputs(input_value, expected_output_type): ], ) def test_value_convert_callable_inputs(input_value, expected_output_type): - """Test that float and value is correctly converted""" + """Test that float and value is correctly converted.""" my_mesh = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, 12) my_t = fem.Constant(my_mesh, default_scalar_type(8)) @@ -117,13 +117,13 @@ def test_value_convert_callable_inputs(input_value, expected_output_type): def test_error_raised_wehn_input_value_is_not_accepted(): - """Test that an error is raised when the input value is not accepted""" + """Test that an error is raised when the input value is not accepted.""" with pytest.raises( TypeError, match=( - "Value must be a float, int, fem.Constant, np.ndarray, fem.Expression, " - "ufl.core.expr.Expr, fem.Function, or callable not coucou" + r"Value must be a float, int, fem.Constant, np.ndarray, fem.Expression, " + r"ufl.core.expr.Expr, fem.Function, or callable not coucou" ), ): F.Value("coucou") @@ -147,7 +147,7 @@ def test_error_raised_wehn_input_value_is_not_accepted(): ], ) def test_time_dependent_values(input_value, expected_output): - """Test that the time_dependent attribute is correctly set""" + """Test that the time_dependent attribute is correctly set.""" test_value = F.Value(input_value) @@ -172,7 +172,7 @@ def test_time_dependent_values(input_value, expected_output): ], ) def test_temperature_dependent_values(input_value, expected_output): - """Test that the time_dependent attribute is correctly set""" + """Test that the time_dependent attribute is correctly set.""" test_value = F.Value(input_value) @@ -187,7 +187,7 @@ def test_temperature_dependent_values(input_value, expected_output): ], ) def test_input_values_of_constants_and_functions_are_accepted(value): - """Test that the input values of constants and functions are accepted""" + """Test that the input values of constants and functions are accepted.""" test_value = F.Value(value) @@ -197,10 +197,11 @@ def test_input_values_of_constants_and_functions_are_accepted(value): def test_input_values_of_expressions_are_accepted(): - """Test that the input values of constants and functions are accepted""" + """Test that the input values of constants and functions are accepted.""" def my_func(x): return 1.0 + x[0] + kwargs = {} kwargs["x"] = x mapped_func = my_func(**kwargs) @@ -218,7 +219,7 @@ def my_func(x): def test_ValueError_raised_when_callable_returns_wrong_type(): """The create_value_fenics method should raise a ValueError when the callable - returns an object which is not a float or int""" + returns an object which is not a float or int.""" def my_value(t): return ufl.conditional(ufl.lt(t, 0.5), 100, 0) @@ -230,7 +231,10 @@ def my_value(t): with pytest.raises( ValueError, - match="self.value should return a float or an int, not " + ), ): test_value.convert_input_value(function_space=test_function_space, t=t) diff --git a/test/test_henrysbc.py b/test/test_henrysbc.py index 998ea78e5..fae1a6e9f 100644 --- a/test/test_henrysbc.py +++ b/test/test_henrysbc.py @@ -10,13 +10,14 @@ def henrys_law(T, H_0, E_H, pressure): - """Applies the Henrys law to compute the concentration at the boundary""" + """Applies the Henrys law to compute the concentration at the boundary.""" H = H_0 * ufl.exp(-E_H / F.k_B / T) return H * pressure def test_raise_error(): - """Test that a value error is raised if the pressure function is not supported in HenrysBC""" + """Test that a value error is raised if the pressure function is not supported in + HenrysBC.""" with pytest.raises(ValueError, match="pressure function not supported"): F.HenrysBC(subdomain=None, H_0=1.0, E_H=1.0, pressure=lambda c: c, species="H") diff --git a/test/test_initial_condition.py b/test/test_initial_condition.py index 50a3a545f..d99defe38 100644 --- a/test/test_initial_condition.py +++ b/test/test_initial_condition.py @@ -15,7 +15,7 @@ def test_init(): - """Test that the attributes are set correctly""" + """Test that the attributes are set correctly.""" # create an InitialConcentration object value = 1.0 species = F.Species("test") @@ -38,8 +38,8 @@ def test_init(): ], ) def test_create_value_fenics(input_value, expected_type): - """Test that after calling .create_expr_fenics, the prev_solution - attribute of the species has the correct value at x=1.0.""" + """Test that after calling .create_expr_fenics, the prev_solution attribute of the + species has the correct value at x=1.0.""" # BUILD @@ -66,7 +66,7 @@ def test_create_value_fenics(input_value, expected_type): def test_warning_raised_when_giving_time_as_arg(): - """Test that a warning is raised if the value is given with t in its arguments""" + """Test that a warning is raised if the value is given with t in its arguments.""" vol = F.VolumeSubdomain(id=1, material=dummy_mat) @@ -83,13 +83,13 @@ def my_value(t): T = fem.Constant(test_mesh.mesh, 10.0) with pytest.raises( - ValueError, match="Initial condition cannot be a function of time." + ValueError, match=r"Initial condition cannot be a function of time." ): init_cond.create_expr_fenics(test_mesh.mesh, T, V) def test_warning_raised_when_giving_time_as_arg_initial_temperature(): - """Test that a warning is raised if the value is given with t in its arguments""" + """Test that a warning is raised if the value is given with t in its arguments.""" # give function to species V = fem.functionspace(test_mesh.mesh, ("Lagrange", 1)) @@ -104,7 +104,7 @@ def my_value(t): init_cond = F.InitialTemperature(value=my_value, volume=vol) with pytest.raises( - ValueError, match="Initial condition cannot be a function of time." + ValueError, match=r"Initial condition cannot be a function of time." ): init_cond.create_expr_fenics(test_mesh.mesh, V) @@ -118,8 +118,8 @@ def my_value(t): ], ) def test_create_value_fenics_initial_temperature(input_value, expected_type): - """Test that after calling .create_expr_fenics, the prev_solution - attribute of the species has the correct value at x=1.0.""" + """Test that after calling .create_expr_fenics, the prev_solution attribute of the + species has the correct value at x=1.0.""" # BUILD @@ -142,10 +142,8 @@ def test_create_value_fenics_initial_temperature(input_value, expected_type): def test_checkpointing_single_species(tmpdir): - """ - Writes a P1 function to a file and reads it back in as initial condition - for one species. - """ + """Writes a P1 function to a file and reads it back in as initial condition for one + species.""" # build initial condition mesh = dolfinx.mesh.create_unit_square( MPI.COMM_WORLD, nx=6, ny=6, cell_type=dolfinx.cpp.mesh.CellType.quadrilateral @@ -193,10 +191,8 @@ def f(x): def test_checkpointing_multiple_species(tmpdir): - """ - Writes two P1 functions to a file and reads them back in as initial conditions - for two species. - """ + """Writes two P1 functions to a file and reads them back in as initial conditions + for two species.""" # build initial condition mesh = dolfinx.mesh.create_unit_square( MPI.COMM_WORLD, nx=6, ny=6, cell_type=dolfinx.cpp.mesh.CellType.quadrilateral @@ -274,11 +270,11 @@ def f2(x): ], ) def test_error_raised_with_volume_setter(vol): - """Test that the volume setter works correctly""" + """Test that the volume setter works correctly.""" spe = F.Species("test") with pytest.raises( - TypeError, match="volume must be of type festim.VolumeSubdomain" + TypeError, match=r"volume must be of type festim.VolumeSubdomain" ): F.InitialConcentration(value=1.0, species=spe, volume=vol) @@ -292,15 +288,15 @@ def test_error_raised_with_volume_setter(vol): ], ) def test_error_raised_with_species_setter(species): - """Test that the volume setter works correctly""" + """Test that the volume setter works correctly.""" vol = F.VolumeSubdomain(id=1, material=dummy_mat) - with pytest.raises(TypeError, match="species must be of type festim.Species"): + with pytest.raises(TypeError, match=r"species must be of type festim.Species"): F.InitialConcentration(value=1.0, species=species, volume=vol) def test_create_initial_temperature_from_function(): - """Test that the initial temperature can be created from a function""" + """Test that the initial temperature can be created from a function.""" # create a volume subdomain vol = F.VolumeSubdomain(id=1, material=dummy_mat) @@ -308,6 +304,7 @@ def test_create_initial_temperature_from_function(): # create an initial temperature from a function def T(x): return 300 + 10 * x[0] + V = fem.functionspace(test_mesh.mesh, ("Lagrange", 1)) T_func = fem.Function(V) T_func.interpolate(T) @@ -323,7 +320,7 @@ def T(x): def test_initial_condition_discontinuous(): """Test the initial condition in a multispecies case with a discontinuous volume - subdomain""" + subdomain.""" my_model = F.HydrogenTransportProblemDiscontinuous() @@ -392,7 +389,7 @@ def test_initial_condition_discontinuous(): def test_initial_condition_continuous_multimaterial(): """Test the initial condition in multi-material continous case that the condition is - only appilied in the correct volume subdomain""" + only appilied in the correct volume subdomain.""" my_model = F.HydrogenTransportProblem() @@ -436,7 +433,7 @@ def test_initial_condition_continuous_multimaterial(): def test_initial_condition_mixed_domain(): - """Test the initial condition in a multi-material discontinous case""" + """Test the initial condition in a multi-material discontinous case.""" my_model = F.HydrogenTransportProblemDiscontinuous() @@ -481,7 +478,7 @@ def test_initial_condition_mixed_domain(): def test_initial_condition_mixed_domain_multispecies(): - """Test the initial condition in a multispecies multi-material discontinous case""" + """Test the initial condition in a multispecies multi-material discontinous case.""" my_model = F.HydrogenTransportProblemDiscontinuous() diff --git a/test/test_material.py b/test/test_material.py index 7c425b558..b6d166344 100644 --- a/test/test_material.py +++ b/test/test_material.py @@ -8,7 +8,7 @@ def test_define_diffusion_coefficient(): - """Test that the diffusion coefficient is correctly defined""" + """Test that the diffusion coefficient is correctly defined.""" T, D_0, E_D = 10, 1.2, 0.5 dum_spe = F.Species("dummy") my_mat = F.Material(D_0=D_0, E_D=E_D) @@ -21,7 +21,7 @@ def test_define_diffusion_coefficient(): def test_multispecies_dict_strings(): """Test that the diffusion coefficient is correctly defined when keys are - strings""" + strings.""" T = 500 D_0_A, D_0_B = 1, 2 E_D_A, E_D_B = 0.1, 0.2 @@ -43,7 +43,7 @@ def test_multispecies_dict_strings(): def test_multispecies_dict_objects(): """Test that the diffusion coefficient is correctly defined when keys are - festim.Species objects""" + festim.Species objects.""" T = 500 D_0_A, D_0_B = 1, 2 E_D_A, E_D_B = 0.1, 0.2 @@ -65,8 +65,8 @@ def test_multispecies_dict_objects(): def test_multispecies_dict_objects_and_strings(): - """Test that the diffusion coefficient is correctly defined when keys - are a mix of festim.Species objects and strings""" + """Test that the diffusion coefficient is correctly defined when keys are a mix of + festim.Species objects and strings.""" T = 500 D_0_A, D_0_B = 1, 2 E_D_A, E_D_B = 0.1, 0.2 @@ -88,8 +88,8 @@ def test_multispecies_dict_objects_and_strings(): def test_multispecies_dict_different_keys(): - """Test that a value error is raised when the keys of the D_0 and E_D - are not the same""" + """Test that a value error is raised when the keys of the D_0 and E_D are not the + same.""" A = F.Species("A") my_mat = F.Material(D_0={"A": 1, "B": 2}, E_D={"A": 0.1, "B": 0.2, "C": 0.3}) @@ -98,8 +98,7 @@ def test_multispecies_dict_different_keys(): def test_D_0_type_raises_error(): - """Test that a value error is raised in the get_diffusion_coefficient - function""" + """Test that a value error is raised in the get_diffusion_coefficient function.""" # TODO remove this when material class is updated A = F.Species("A") my_mat = F.Material(D_0=[1, 1], E_D=0.1) @@ -109,8 +108,8 @@ def test_D_0_type_raises_error(): def test_error_raised_when_species_not_given_with_dict(): - """Test that a value error is raised when a species has not been given in - the get_diffusion_coefficient function when using a dict for properties""" + """Test that a value error is raised when a species has not been given in the + get_diffusion_coefficient function when using a dict for properties.""" A = F.Species("A") B = F.Species("B") my_mat = F.Material(D_0={A: 1, B: 1}, E_D={A: 0.1, B: 0.1}) @@ -122,8 +121,8 @@ def test_error_raised_when_species_not_given_with_dict(): def test_error_raised_when_species_not_not_in_D_0_dict(): - """Test that a value error is raised when a species has not been given but - has no value in the dict""" + """Test that a value error is raised when a species has not been given but has no + value in the dict.""" A = F.Species("A") B = F.Species("B") J = F.Species("J") @@ -134,8 +133,7 @@ def test_error_raised_when_species_not_not_in_D_0_dict(): def test_D_0_raises_ValueError_if_species_not_provided_in_dict(): - """Test that a value error is raised in the get_diffusion_coefficient - function""" + """Test that a value error is raised in the get_diffusion_coefficient function.""" # TODO remove this when material class is updated A = F.Species("A") B = F.Species("B") @@ -146,8 +144,7 @@ def test_D_0_raises_ValueError_if_species_not_provided_in_dict(): def test_D_0_raises_ValueError_if_species_given_not_in_dict_keys(): - """Test that a value error is raised in the get_diffusion_coefficient - function""" + """Test that a value error is raised in the get_diffusion_coefficient function.""" # TODO remove this when material class is updated A = F.Species("A") B = F.Species("B") @@ -159,7 +156,7 @@ def test_D_0_raises_ValueError_if_species_given_not_in_dict_keys(): def test_raises_TypeError_when_D_0_is_not_correct_type(): - """Test that a TypeError is raised when D_0 is not a float or a dict""" + """Test that a TypeError is raised when D_0 is not a float or a dict.""" my_mat = F.Material(D_0=[1, 2], E_D=1) @@ -168,8 +165,7 @@ def test_raises_TypeError_when_D_0_is_not_correct_type(): def test_E_D_raises_ValueError_if_species_not_provided_in_dict(): - """Test that a value error is raised in the get_diffusion_coefficient - function""" + """Test that a value error is raised in the get_diffusion_coefficient function.""" # TODO remove this when material class is updated A = F.Species("A") B = F.Species("B") @@ -180,8 +176,7 @@ def test_E_D_raises_ValueError_if_species_not_provided_in_dict(): def test_E_D_raises_ValueError_if_species_given_not_in_dict_keys(): - """Test that a value error is raised in the get_diffusion_coefficient - function""" + """Test that a value error is raised in the get_diffusion_coefficient function.""" # TODO remove this when material class is updated A = F.Species("A") B = F.Species("B") @@ -193,7 +188,7 @@ def test_E_D_raises_ValueError_if_species_given_not_in_dict_keys(): def test_raises_TypeError_when_E_D_is_not_correct_type(): - """Test that a TypeError is raised when E_D is not a float or a dict""" + """Test that a TypeError is raised when E_D is not a float or a dict.""" my_mat = F.Material(D_0=1, E_D=[1, 2]) @@ -211,15 +206,14 @@ def test_raises_TypeError_when_E_D_is_not_correct_type(): ], ) def test_raises_TypeError_when_D_is_not_correct_type(input_value): - """Test that a TypeError is raised when D is not an fem.Function""" + """Test that a TypeError is raised when D is not an fem.Function.""" - with pytest.raises(TypeError, match="D must be of type fem.Function"): + with pytest.raises(TypeError, match=r"D must be of type fem.Function"): F.Material(D=input_value) def test_get_diffusion_coefficient_returns_function_when_given_to_D(): - """Test that the diffusion coefficient is correctly defined when D is a - function""" + """Test that the diffusion coefficient is correctly defined when D is a function.""" V = fem.functionspace(test_mesh.mesh, ("Lagrange", 1)) D = fem.Function(V) @@ -232,23 +226,23 @@ def test_get_diffusion_coefficient_returns_function_when_given_to_D(): def test_error_raised_when_D_and_D_0_given(): - """Test that a value error is raised when both D and D_0 are given""" + """Test that a value error is raised when both D and D_0 are given.""" V = fem.functionspace(test_mesh.mesh, ("Lagrange", 1)) D_func = fem.Function(V) with pytest.raises( ValueError, - match="D_0 and D cannot be set at the same time. Please set only one of them.", + match=r"D_0 and D cannot be set at the same time. Please set only one of them.", ): F.Material(D=D_func, D_0=1) def test_error_raised_when_D_and_D_0_both_None(): - """Test that a value error is raised when both D and D_0 are None""" + """Test that a value error is raised when both D and D_0 are None.""" with pytest.raises( ValueError, - match="D_0 and D cannot both be None. Please set one of them.", + match=r"D_0 and D cannot both be None. Please set one of them.", ): F.Material() diff --git a/test/test_maximum_surface.py b/test/test_maximum_surface.py index 1636c8f66..096c06b70 100644 --- a/test/test_maximum_surface.py +++ b/test/test_maximum_surface.py @@ -5,7 +5,7 @@ def test_maximum_surface_compute_1D(): - """Test that the maximum surface export computes the correct value""" + """Test that the maximum surface export computes the correct value.""" # BUILD L = 4.0 diff --git a/test/test_maximum_volume.py b/test/test_maximum_volume.py index 393e5b933..ef47c1d68 100644 --- a/test/test_maximum_volume.py +++ b/test/test_maximum_volume.py @@ -5,7 +5,7 @@ def test_maximum_volume_compute_1D(): - """Test that the maximum volume export computes the right value""" + """Test that the maximum volume export computes the right value.""" # BUILD L = 6 diff --git a/test/test_mesh.py b/test/test_mesh.py index f00872ef6..a087149f2 100644 --- a/test/test_mesh.py +++ b/test/test_mesh.py @@ -1,3 +1,4 @@ +# ruff: noqa: E402 import logging import os @@ -76,7 +77,8 @@ def test_vdim_changes_when_mesh_changes(): @pytest.mark.parametrize("mesh", [mesh_1D, mesh_2D, mesh_3D]) def test_meshtags_from_xdmf(tmp_path, mesh): - """Test that the facet and volume meshtags are read correctly from the mesh XDMF files""" + """Test that the facet and volume meshtags are read correctly from the mesh XDMF + files.""" # create mesh functions fdim = mesh.topology.dim - 1 vdim = mesh.topology.dim @@ -160,23 +162,24 @@ def test_meshtags_from_xdmf(tmp_path, mesh): @pytest.mark.parametrize("vertices", [[1, 2, 3, 4], [0, 0.1, 0.2, 0.3, 0.4, 0.5]]) def test_mesh_vertices_from_list(vertices): - """Check that giving vertices as a list is correctly processed and ends up as a np.ndarray for the mesh""" + """Check that giving vertices as a list is correctly processed and ends up as a + np.ndarray for the mesh.""" my_mesh = F.Mesh1D(vertices=vertices) assert isinstance(my_mesh.vertices, np.ndarray) def test_error_raised_when_mesh_is_wrong_type(): - """Test that an TypeError is raised when the mesh is not a dolfinx mesh""" + """Test that an TypeError is raised when the mesh is not a dolfinx mesh.""" - with pytest.raises(TypeError, match="Mesh must be of type dolfinx.mesh.Mesh"): + with pytest.raises(TypeError, match=r"Mesh must be of type dolfinx.mesh.Mesh"): F.Mesh( mesh="mesh", ) def test_create_1D_mesh_parallel(cluster): - """Test creating a 1D mesh in parallel using ipyparallel""" + """Test creating a 1D mesh in parallel using ipyparallel.""" def create_mesh(): import numpy as np @@ -193,7 +196,7 @@ def create_mesh(): @pytest.mark.parametrize("mesh", [mesh_2D, mesh_3D]) def test_attr_error_with_incompitable_mesh_in_spherical(mesh): """Test that an AttributeError is raised when trying to use an incompatible mesh - with a spherical coordinate system""" + with a spherical coordinate system.""" with pytest.raises( AttributeError, @@ -204,7 +207,7 @@ def test_attr_error_with_incompitable_mesh_in_spherical(mesh): def test_attr_error_with_3D_mesh_in_cylindrical(): """Test that an AttributeError is raised when trying to use an incompatible mesh - with a cylindrical coordinate system""" + with a cylindrical coordinate system.""" with pytest.raises( AttributeError, @@ -216,7 +219,7 @@ def test_attr_error_with_3D_mesh_in_cylindrical(): @pytest.mark.parametrize("system", ["cyl", "cart", 1.0, "coucou", mesh_1D]) def test_coordinate_system_setter(system): if isinstance(system, str): - err_msg = "coordinate_system must be one of 'cartesian', 'cylindrical', or 'spherical'" + err_msg = "coordinate_system must be one of 'cartesian', 'cylindrical', or 'spherical'" # noqa: E501 else: err_msg = "coordinate_system must be of type str or CoordinateSystem" with pytest.raises( diff --git a/test/test_minimum_surface.py b/test/test_minimum_surface.py index e7ce02808..416481819 100644 --- a/test/test_minimum_surface.py +++ b/test/test_minimum_surface.py @@ -5,7 +5,7 @@ def test_minimum_surface_export_compute_1D(): - """Test that the minimum surface export computes the correct value""" + """Test that the minimum surface export computes the correct value.""" # BUILD L = 4.0 diff --git a/test/test_minimum_volume.py b/test/test_minimum_volume.py index f9e5a985d..a6e42ccbd 100644 --- a/test/test_minimum_volume.py +++ b/test/test_minimum_volume.py @@ -5,7 +5,7 @@ def test_minimum_volume_compute_1D(): - """Test that the minimum volume export computes the right value""" + """Test that the minimum volume export computes the right value.""" # BUILD L = 6 diff --git a/test/test_multispecies_problem.py b/test/test_multispecies_problem.py index 4defae394..6529e81f8 100644 --- a/test/test_multispecies_problem.py +++ b/test/test_multispecies_problem.py @@ -33,9 +33,9 @@ def relative_error_computed_to_analytical( def test_multispecies_permeation_problem(): - """Test running a problem with 2 mobile species permeating through a 1D - domain, with different diffusion coefficients, checks that the computed - permeation flux matches the analytical solution""" + """Test running a problem with 2 mobile species permeating through a 1D domain, with + different diffusion coefficients, checks that the computed permeation flux matches + the analytical solution.""" # festim model L = 3e-04 diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index d5da1e391..011821539 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -36,9 +36,8 @@ def relative_error_computed_to_analytical( def test_permeation_problem(mesh_size=1001): - """Test running a problem with a mobile species permeating through a 1D - domain, checks that the computed permeation flux matches the analytical - solution""" + """Test running a problem with a mobile species permeating through a 1D domain, + checks that the computed permeation flux matches the analytical solution.""" # festim model L = 3e-04 @@ -107,7 +106,7 @@ def test_permeation_problem(mesh_size=1001): times = outgassing_flux.t flux_values = outgassing_flux.data - # -------------------------- analytical solution ------------------------------------- + # -------------------------- analytical solution ------------------------------------- # noqa: E501 D = my_mat.get_diffusion_coefficient(my_mesh.mesh, my_model.temperature) @@ -127,8 +126,7 @@ def test_permeation_problem(mesh_size=1001): def test_permeation_problem_multi_volume(tmp_path): - """Same permeation problem as above but with 4 volume subdomains instead - of 1""" + """Same permeation problem as above but with 4 volume subdomains instead of 1.""" L = 3e-04 vertices = np.linspace(0, L, num=1001) diff --git a/test/test_reaction.py b/test/test_reaction.py index 397d82faa..850214a58 100644 --- a/test/test_reaction.py +++ b/test/test_reaction.py @@ -12,7 +12,7 @@ def test_reaction_init(): - """Test that the Reaction class initialises correctly""" + """Test that the Reaction class initialises correctly.""" # create two species species1 = F.Species("A") species2 = F.Species("B") @@ -41,7 +41,7 @@ def test_reaction_init(): def test_reaction_repr(): - """Test that the Reaction __repr__ method returns the expected string""" + """Test that the Reaction __repr__ method returns the expected string.""" # create two species species1 = F.Species("A") @@ -67,7 +67,7 @@ def test_reaction_repr(): def test_reaction_repr_2_products(): - """Test that the Reaction __repr__ method returns the expected string""" + """Test that the Reaction __repr__ method returns the expected string.""" # create two species species1 = F.Species("A") @@ -94,7 +94,7 @@ def test_reaction_repr_2_products(): def test_reaction_repr_0_products(): - """Test that the Reaction __repr__ method returns the expected string""" + """Test that the Reaction __repr__ method returns the expected string.""" # create two species species1 = F.Species("A") @@ -113,7 +113,7 @@ def test_reaction_repr_0_products(): def test_reaction_str(): - """Test that the Reaction __str__ method returns the expected string""" + """Test that the Reaction __str__ method returns the expected string.""" # create two species species1 = F.Species("A") @@ -139,7 +139,8 @@ def test_reaction_str(): def test_reaction_str_2_products(): - """Test that the Reaction __str__ method returns the expected string when there are 2 products""" + """Test that the Reaction __str__ method returns the expected string when there are + 2 products.""" # create two species species1 = F.Species("A") @@ -166,7 +167,8 @@ def test_reaction_str_2_products(): def test_reaction_str_no_products(): - """Test that the Reaction __str__ method returns the expected string when there are 2 products""" + """Test that the Reaction __str__ method returns the expected string when there are + 2 products.""" # create two species species1 = F.Species("A") @@ -186,7 +188,8 @@ def test_reaction_str_no_products(): @pytest.mark.parametrize("temperature", [300.0, 350, 370, 500.0]) def test_reaction_reaction_term(temperature): - """Test that the Reaction.reaction_term method returns the expected reaction term""" + """Test that the Reaction.reaction_term method returns the expected reaction + term.""" mesh = create_unit_cube(MPI.COMM_WORLD, 10, 10, 10) V = functionspace(mesh, ("Lagrange", 1)) @@ -258,7 +261,8 @@ def arrhenius(pre, act, T): @pytest.mark.parametrize("temperature", [300.0, 350, 370, 500.0]) def test_reaction_reaction_term_2_products(temperature): - """Test that the Reaction.reaction_term method returns the expected reaction term with two products""" + """Test that the Reaction.reaction_term method returns the expected reaction term + with two products.""" mesh = create_unit_cube(MPI.COMM_WORLD, 10, 10, 10) V = functionspace(mesh, ("Lagrange", 1)) @@ -304,7 +308,10 @@ def test_reactant_setter_raises_error_with_zero_length_list(): """Test a value error is raised when the first reactant is given a wrong type.""" with pytest.raises( ValueError, - match="reactant must be an entry of one or more species objects, not an empty list.", + match=( + r"reactant must be an entry of one or more species objects, not an empty " + r"list." + ), ): F.Reaction( reactant=[], @@ -320,7 +327,7 @@ def test_reactant_setter_raises_error_with_wrong_type(): """Test a type error is raised when the first reactant is given a wrong type.""" with pytest.raises( TypeError, - match="reactant must be an F.Species or F.ImplicitSpecies, not ", + match=r"reactant must be an F.Species or F.ImplicitSpecies, not ", ): F.Reaction( reactant=["A", F.Species("B")], @@ -336,7 +343,7 @@ def test_reactant_setter_raises_error_with_wrong_type(): def test_product_setter_raise_error_p_0_no_product(): with pytest.raises( ValueError, - match="p_0 must be None, not 2 when no products are present.", + match=r"p_0 must be None, not 2 when no products are present.", ): reaction = F.Reaction( reactant=[F.Species("A")], @@ -351,7 +358,7 @@ def test_product_setter_raise_error_p_0_no_product(): def test_no_E_p_with_product(): with pytest.raises( ValueError, - match="E_p cannot be None when reaction products are present.", + match=r"E_p cannot be None when reaction products are present.", ): reaction = F.Reaction( reactant=[F.Species("A")], @@ -367,7 +374,7 @@ def test_no_E_p_with_product(): def test_no_p_0_with_product(): with pytest.raises( ValueError, - match="p_0 cannot be None when reaction products are present.", + match=r"p_0 cannot be None when reaction products are present.", ): reaction = F.Reaction( reactant=[F.Species("A")], @@ -383,7 +390,7 @@ def test_no_p_0_with_product(): def test_product_setter_raise_error_E_p_no_product(): with pytest.raises( ValueError, - match="E_p must be None, not 2 when no products are present.", + match=r"E_p must be None, not 2 when no products are present.", ): reaction = F.Reaction( reactant=[F.Species("A")], @@ -430,11 +437,9 @@ def test_product_setter_raise_error_E_p_no_product(): @pytest.mark.parametrize("reaction", [reac1, reac2]) def test_override_solution_attributes(reaction): - """ - Tests the HydrogenTransportProblemDiscontinuous.override_solution_attributes method - Checks that the .solution attribute is the expected one based on the volume of the - reaction - """ + """Tests the HydrogenTransportProblemDiscontinuous.override_solution_attributes + method Checks that the .solution attribute is the expected one based on the volume + of the reaction.""" # RUN my_model.override_solution_attributes(reaction) diff --git a/test/test_settings.py b/test/test_settings.py index 1a2275420..4d6a0bbeb 100644 --- a/test/test_settings.py +++ b/test/test_settings.py @@ -6,7 +6,7 @@ @pytest.mark.parametrize("test_type", [int, F.Stepsize, float]) def test_stepsize_value(test_type): - """Test that the stepsize is correctly set""" + """Test that the stepsize is correctly set.""" test_value = 23.0 my_settings = F.Settings(atol=1, rtol=0.1) my_settings.stepsize = test_type(test_value) @@ -16,7 +16,7 @@ def test_stepsize_value(test_type): def test_stepsize_value_wrong_type(): - """Checks that an error is raised when the wrong type is given""" + """Checks that an error is raised when the wrong type is given.""" my_settings = F.Settings(atol=1, rtol=0.1) with pytest.raises(TypeError): @@ -49,7 +49,7 @@ def test_stepsize_value_wrong_type(): # ], # ) # def test_tolerances_value(rtol, atol): -# """Tests that callable tolerances are called & return correct float before passed to fenics""" +# """Tests that callable tolerances are called & return correct float before passed to fenics""" # noqa: E501 # # BUILD # test_mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0])) diff --git a/test/test_sievertsbc.py b/test/test_sievertsbc.py index a70d697a5..992104887 100644 --- a/test/test_sievertsbc.py +++ b/test/test_sievertsbc.py @@ -10,13 +10,14 @@ def sieverts_law(T, S_0, E_S, pressure): - """Applies the Sieverts law to compute the concentration at the boundary""" + """Applies the Sieverts law to compute the concentration at the boundary.""" S = S_0 * ufl.exp(-E_S / F.k_B / T) return S * pressure**0.5 def test_raise_error(): - """Test that a value error is raised if the pressure function is not supported in SievertsBC""" + """Test that a value error is raised if the pressure function is not supported in + SievertsBC.""" with pytest.raises(ValueError, match="pressure function not supported"): F.SievertsBC( subdomain=None, S_0=1.0, E_S=1.0, pressure=lambda c: c, species="H" diff --git a/test/test_source.py b/test/test_source.py index 676d346a9..7290e9998 100644 --- a/test/test_source.py +++ b/test/test_source.py @@ -13,7 +13,7 @@ def test_init(): - """Test that the attributes are set correctly""" + """Test that the attributes are set correctly.""" # create a Source object volume = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=dummy_mat) value = 1.0 @@ -46,7 +46,7 @@ def test_init(): ], ) def test_create_fenics_object(value, expected_type): - """Test that the correct fenics object is created depending on the value input""" + """Test that the correct fenics object is created depending on the value input.""" # BUILD vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat) @@ -83,7 +83,7 @@ def test_create_fenics_object(value, expected_type): ], ) def test_source_explicit_time_dependent_attribute(input, expected_value): - """Test that the time_dependent attribute is correctly set""" + """Test that the time_dependent attribute is correctly set.""" volume = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat) species = F.Species("test") my_source = F.ParticleSource(input, volume, species) @@ -108,7 +108,7 @@ def test_source_explicit_time_dependent_attribute(input, expected_value): ], ) def test_source_temperature_dependent_attribute(input, expected_value): - """Test that the temperature_dependent attribute is correctly set""" + """Test that the temperature_dependent attribute is correctly set.""" volume = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat) species = F.Species("test") my_source = F.ParticleSource(input, volume, species) @@ -117,8 +117,8 @@ def test_source_temperature_dependent_attribute(input, expected_value): def test_ValueError_raised_when_callable_returns_wrong_type(): - """The create_value method should raise a ValueError when the callable - returns an object which is not a float or int""" + """The create_value method should raise a ValueError when the callable returns an + object which is not a float or int.""" vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat) species = F.Species("test") @@ -135,7 +135,10 @@ def my_value(t): with pytest.raises( ValueError, - match="self.value should return a float or an int, not Date: Wed, 22 Apr 2026 15:22:34 -0400 Subject: [PATCH 10/14] fixed conf.py --- docs/source/conf.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/docs/source/conf.py b/docs/source/conf.py index 3e5c5c7d6..2016a27ea 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -12,8 +12,6 @@ import os import sys -import map - sys.path.insert(0, os.path.abspath("../../src")) @@ -22,6 +20,7 @@ # Add the directory containing your Python script to the Python path sys.path.insert(0, os.path.abspath(".")) +import map m = map.generate_map() current_dir = os.path.dirname(__file__) From a2f81e229195090908110ca5c89991ead4a3d66d Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 23 Apr 2026 10:58:00 -0400 Subject: [PATCH 11/14] fixed line length --- src/festim/trap.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/festim/trap.py b/src/festim/trap.py index 809baf333..653c06fcc 100644 --- a/src/festim/trap.py +++ b/src/festim/trap.py @@ -35,8 +35,8 @@ class Trap(_Species): .. testsetup:: Trap - from festim import Trap, Species, VolumeSubdomain, Material, - HydrogenTransportProblem + from festim import Trap, Species, VolumeSubdomain, Material + from fsetim import HydrogenTransportProblem my_mat = Material(D_0=1, E_D=1, name="test_mat") my_vol = VolumeSubdomain(id=1, material=my_mat) From 1c827620785aa4549b99124b8559c2e788b146ce Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 Apr 2026 14:44:45 -0400 Subject: [PATCH 12/14] fixed typo --- src/festim/trap.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/festim/trap.py b/src/festim/trap.py index 653c06fcc..74bd514a7 100644 --- a/src/festim/trap.py +++ b/src/festim/trap.py @@ -36,7 +36,7 @@ class Trap(_Species): .. testsetup:: Trap from festim import Trap, Species, VolumeSubdomain, Material - from fsetim import HydrogenTransportProblem + from festim import HydrogenTransportProblem my_mat = Material(D_0=1, E_D=1, name="test_mat") my_vol = VolumeSubdomain(id=1, material=my_mat) From b21780f441cc172da8bff1ff6a1c59486d1b1042 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 Apr 2026 14:46:40 -0400 Subject: [PATCH 13/14] final fix --- src/festim/subdomain/volume_subdomain.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/festim/subdomain/volume_subdomain.py b/src/festim/subdomain/volume_subdomain.py index 1b33682cf..d1135b5da 100644 --- a/src/festim/subdomain/volume_subdomain.py +++ b/src/festim/subdomain/volume_subdomain.py @@ -198,7 +198,8 @@ def map_surface_to_volume_subdomains( Raises: - AssertionError: if a surface subdomain is connected to multiple volume subdomains + AssertionError: if a surface subdomain is connected to multiple volume + subdomains Args: ft: the facet meshtags of the parent mesh @@ -263,7 +264,8 @@ def map_surface_to_volume_subdomains( if s_subdomain and v_subdomain: if s_subdomain in surface_to_subdomain: assert surface_to_subdomain[s_subdomain] == v_subdomain, ( - f"Surface subdomain {s_subdomain.id} is connected to multiple volume subdomains: " + f"Surface subdomain {s_subdomain.id} is connected " + f"to multiple volume subdomains: " f"{surface_to_subdomain[s_subdomain].id} and {v_subdomain.id}" ) else: From bd9d0361915d93e239fcf0f87edb098f033b0900 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 Apr 2026 14:50:17 -0400 Subject: [PATCH 14/14] removed all exampels --- .../advection_diffusion_multi_material.py | 159 ------------- examples/convergence_rates.py | 157 ------------ examples/multi_isotope_trapping_example.py | 155 ------------ examples/multi_material_1d.py | 99 -------- examples/multi_material_2d.py | 202 ---------------- examples/multi_material_2d_bis.py | 224 ------------------ examples/multi_material_transient.py | 99 -------- examples/multi_material_with_one_material.py | 65 ----- examples/plot_tds.py | 19 -- examples/surface_reaction_example.py | 169 ------------- examples/tds_example.py | 136 ----------- examples/trapping_example.py | 80 ------- 12 files changed, 1564 deletions(-) delete mode 100644 examples/advection_diffusion_multi_material.py delete mode 100644 examples/convergence_rates.py delete mode 100644 examples/multi_isotope_trapping_example.py delete mode 100644 examples/multi_material_1d.py delete mode 100644 examples/multi_material_2d.py delete mode 100644 examples/multi_material_2d_bis.py delete mode 100644 examples/multi_material_transient.py delete mode 100644 examples/multi_material_with_one_material.py delete mode 100644 examples/plot_tds.py delete mode 100644 examples/surface_reaction_example.py delete mode 100644 examples/tds_example.py delete mode 100644 examples/trapping_example.py diff --git a/examples/advection_diffusion_multi_material.py b/examples/advection_diffusion_multi_material.py deleted file mode 100644 index dccfb0385..000000000 --- a/examples/advection_diffusion_multi_material.py +++ /dev/null @@ -1,159 +0,0 @@ -from mpi4py import MPI - -import basix -import dolfinx -import dolfinx.fem.petsc -import numpy as np - -import festim as F - - -# ---------------- Generate a mesh ---------------- -def generate_mesh(n=20): - def bottom_boundary(x): - return np.isclose(x[1], 0.0) - - def top_boundary(x): - return np.isclose(x[1], 1.0) - - def bottom_left_boundary(x): - return np.logical_and(np.isclose(x[0], 0.0), x[1] <= 0.5 + 1e-14) - - def half(x): - return x[1] <= 0.5 + 1e-14 - - mesh = dolfinx.mesh.create_unit_square( - MPI.COMM_WORLD, n, n, dolfinx.mesh.CellType.triangle - ) - - # Split domain in half and set an interface tag of 5 - tdim = mesh.topology.dim - fdim = tdim - 1 - top_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, top_boundary) - bottom_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, bottom_boundary) - bot_left_facets = dolfinx.mesh.locate_entities_boundary( - mesh, fdim, bottom_left_boundary - ) - num_facets_local = ( - mesh.topology.index_map(fdim).size_local - + mesh.topology.index_map(fdim).num_ghosts - ) - facets = np.arange(num_facets_local, dtype=np.int32) - values = np.full_like(facets, 0, dtype=np.int32) - values[top_facets] = 1 - values[bottom_facets] = 2 - values[bot_left_facets] = 3 - - bottom_cells = dolfinx.mesh.locate_entities(mesh, tdim, half) - num_cells_local = ( - mesh.topology.index_map(tdim).size_local - + mesh.topology.index_map(tdim).num_ghosts - ) - cells = np.full(num_cells_local, 5, dtype=np.int32) - cells[bottom_cells] = 4 - ct = dolfinx.mesh.meshtags( - mesh, tdim, np.arange(num_cells_local, dtype=np.int32), cells - ) - all_b_facets = dolfinx.mesh.compute_incident_entities( - mesh.topology, ct.find(4), tdim, fdim - ) - all_t_facets = dolfinx.mesh.compute_incident_entities( - mesh.topology, ct.find(5), tdim, fdim - ) - interface = np.intersect1d(all_b_facets, all_t_facets) - values[interface] = 6 - - mt = dolfinx.mesh.meshtags(mesh, mesh.topology.dim - 1, facets, values) - return mesh, mt, ct - - -mesh, mt, ct = generate_mesh() - - -# create velocity field -def create_velocity_field(): - def velocity_func(x): - values = np.zeros((2, x.shape[1])) # Initialize with zeros - scalar_value = -500 * x[1] * (x[1] - 0.5) # Compute the scalar function - values[0] = scalar_value # Assign to first component - values[1] = 0 # Second component remains zero - return values - - mesh2, _mt2, ct2 = generate_mesh(n=10) - submesh, _cell_map, _v_map = dolfinx.mesh.create_submesh( - mesh2, ct2.dim, ct2.find(4) - )[0:3] - v_cg = basix.ufl.element( - "Lagrange", submesh.topology.cell_name(), 2, shape=(submesh.geometry.dim,) - ) - V_velocity = dolfinx.fem.functionspace(submesh, v_cg) - u = dolfinx.fem.Function(V_velocity) - u.interpolate(velocity_func) - return u - - -u = create_velocity_field() - -writer = dolfinx.io.VTXWriter(MPI.COMM_WORLD, "velocity.bp", u, engine="BP5") -writer.write(t=0) - -my_model = F.HydrogenTransportProblemDiscontinuous() -my_model.mesh = F.Mesh(mesh) -my_model.volume_meshtags = ct -my_model.facet_meshtags = mt - -material_bottom = F.Material(D_0=1, E_D=0 * 0.1) -material_top = F.Material(D_0=1, E_D=0 * 0.1) - -material_bottom.K_S_0 = 2.0 -material_bottom.E_K_S = 0 -material_top.K_S_0 = 600 -material_top.E_K_S = 0 - -top_domain = F.VolumeSubdomain(5, material=material_top) -bottom_domain = F.VolumeSubdomain(4, material=material_bottom) - -top_surface = F.SurfaceSubdomain(id=1) -bottom_surface = F.SurfaceSubdomain(id=2) -inlet_surf = F.SurfaceSubdomain(id=3) -my_model.subdomains = [ - bottom_domain, - top_domain, - top_surface, - bottom_surface, - inlet_surf, -] - -# we should be able to automate this -my_model.interfaces = [F.Interface(6, (bottom_domain, top_domain))] - -H = F.Species("H", mobile=True) - -my_model.species = [H] - -for species in my_model.species: - species.subdomains = [bottom_domain, top_domain] - - -my_model.boundary_conditions = [ - F.FixedConcentrationBC(top_surface, value=0.0, species=H), - F.FixedConcentrationBC(inlet_surf, value=1.0, species=H), - F.FixedConcentrationBC(bottom_surface, value=0.0, species=H), -] - -my_model.advection_terms = [ - F.AdvectionTerm(velocity=u, species=[H], subdomain=bottom_domain) -] - - -my_model.temperature = 500.0 - -my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False) - -my_model.exports = [ - F.VTXSpeciesExport(f"u_{subdomain.id}.bp", field=H, subdomain=subdomain) - for subdomain in my_model.volume_subdomains -] - -my_model.initialise() -my_model.run() diff --git a/examples/convergence_rates.py b/examples/convergence_rates.py deleted file mode 100644 index a6b24d702..000000000 --- a/examples/convergence_rates.py +++ /dev/null @@ -1,157 +0,0 @@ -import mpi4py.MPI as MPI - -import matplotlib.pyplot as plt -import numpy as np -import ufl -from dolfinx import fem - -import festim as F - - -def source_from_exact_solution( - exact_solution, thermal_conductivity, density, heat_capacity -): - import sympy as sp - from sympy.vector import CoordSys3D, divergence, gradient - - R = CoordSys3D("R") - x = [R.x, R.y, R.z] - t = sp.symbols("t") - - u = exact_solution(x, t) - density = density(x, t) - heat_capacity = heat_capacity(x, t) - thermal_cond = thermal_conductivity(x, t) - source = density * heat_capacity * sp.diff(u, t) - divergence( - thermal_cond * gradient(u, x), x - ) - - return sp.lambdify([x, t], source) - - -def error_L2(u_computed, u_exact, degree_raise=3): - # Create higher order function space - degree = u_computed.function_space.ufl_element().degree() - family = u_computed.function_space.ufl_element().family() - mesh = u_computed.function_space.mesh - W = fem.FunctionSpace(mesh, (family, degree + degree_raise)) - # Interpolate approximate solution - u_W = fem.Function(W) - u_W.interpolate(u_computed) - - # Interpolate exact solution, special handling if exact solution - # is a ufl expression or a python lambda function - u_ex_W = fem.Function(W) - if isinstance(u_exact, ufl.core.expr.Expr): - u_expr = fem.Expression(u_exact, W.element.interpolation_points) - u_ex_W.interpolate(u_expr) - else: - u_ex_W.interpolate(u_exact) - - # Compute the error in the higher order function space - e_W = fem.Function(W) - e_W.x.array[:] = u_W.x.array - u_ex_W.x.array - - # Integrate the error - error = fem.form(ufl.inner(e_W, e_W) * ufl.dx) - error_local = fem.assemble_scalar(error) - error_global = mesh.comm.allreduce(error_local, op=MPI.SUM) - return np.sqrt(error_global) - - -def run(N): - def exact_solution(x, t): - return 2 * x[0] ** 2 + 20 * t - - def density(T): - return 0.2 * T + 2 - - def heat_capacity(T): - return 0.2 * T + 3 - - def thermal_conductivity(T): - return 0.1 * T + 4 - - mms_source_from_sp = source_from_exact_solution( - exact_solution, - density=lambda x, t: density(exact_solution(x, t)), - heat_capacity=lambda x, t: heat_capacity(exact_solution(x, t)), - thermal_conductivity=lambda x, t: thermal_conductivity(exact_solution(x, t)), - ) - - def mms_source(x, t): - return mms_source_from_sp((x[0], None, None), t) - - my_problem = F.HeatTransferProblem() - - my_problem.mesh = F.Mesh1D(vertices=np.linspace(2, 3, N)) - left = F.SurfaceSubdomain1D(id=1, x=2) - right = F.SurfaceSubdomain1D(id=2, x=3) - my_problem.surface_subdomains = [left, right] - mat = F.Material(D_0=None, E_D=None) - mat.thermal_conductivity = thermal_conductivity - mat.density = density - mat.heat_capacity = heat_capacity - - my_problem.volume_subdomains = [ - F.VolumeSubdomain1D(id=1, borders=[2, 3], material=mat) - ] - - # NOTE: it's good to check that without the IC the solution is not the exact one - my_problem.initial_condition = F.InitialTemperature(lambda x: exact_solution(x, 0)) - - my_problem.boundary_conditions = [ - F.FixedTemperatureBC(subdomain=left, value=exact_solution), - F.FixedTemperatureBC(subdomain=right, value=exact_solution), - ] - - my_problem.sources = [ - F.HeatSource(value=mms_source, volume=my_problem.volume_subdomains[0]) - ] - - my_problem.settings = F.Settings( - atol=1e-8, - rtol=1e-10, - # final time shouldn't be too long so that a potential error at the - # initial timestep is not negligible - final_time=1, - ) - - # Forward euler isn't great so dt should be small - # although it's ok here since the time derivative is constant - my_problem.settings.stepsize = F.Stepsize(0.05) - - my_problem.exports = [ - F.VTXSpeciesExport(filename="test_transient_heat_transfer.bp") - ] - - my_problem.initialise() - my_problem.run() - - computed_solution = my_problem.u - # we use the exact final time of the simulation which may differ from - # the one specified in the settings - final_time_sim = my_problem.t.value - - def exact_solution_end(x): - return exact_solution(x, final_time_sim) - - L2_error = error_L2(computed_solution, exact_solution_end) - return L2_error - - -if __name__ == "__main__": - Ns = np.geomspace(10, 3300, 10) - errors = [] - for N in Ns: - L2_error = run(int(N)) - print(f"With {int(N)} cells, L2_error = {L2_error}") - errors.append(L2_error) - - plt.loglog(Ns, errors, "-o", label="L2") - # plot a slope of 2 - plt.loglog(Ns, 1 / Ns**2, "--", color="black", label="order 2") - plt.xlabel("Number of cells") - plt.ylabel("Error") - plt.legend() - plt.show() diff --git a/examples/multi_isotope_trapping_example.py b/examples/multi_isotope_trapping_example.py deleted file mode 100644 index 209929f71..000000000 --- a/examples/multi_isotope_trapping_example.py +++ /dev/null @@ -1,155 +0,0 @@ -import numpy as np - -import festim as F - -my_model = F.HydrogenTransportProblem() - -# -------- Mesh --------- # - -L = 5e-6 -vertices = np.linspace(0, L, num=2000) -my_model.mesh = F.Mesh1D(vertices) - - -# -------- Materials and subdomains --------- # - -w_atom_density = 6.3e28 # atom/m3 - -tungsten = F.Material(D_0=4.1e-7, E_D=0.39, name="tungsten") - -my_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, L], material=tungsten) -left_surface = F.SurfaceSubdomain1D(id=1, x=0) -right_surface = F.SurfaceSubdomain1D(id=2, x=L) - -my_model.subdomains = [ - my_subdomain, - left_surface, - right_surface, -] - -# -------- Hydrogen species and reactions --------- # - -mobile_H = F.Species("H") -mobile_D = F.Species("D") - - -trapped_H1 = F.Species("trapped_H1", mobile=False) -trapped_D1 = F.Species("trapped_D1", mobile=False) -trapped_H2 = F.Species("trapped_H2", mobile=False) -trapped_D2 = F.Species("trapped_D2", mobile=False) -trapped_HD = F.Species("trapped_HD", mobile=False) - -empty_trap = F.ImplicitSpecies( - n=1e21, - others=[trapped_H1, trapped_D1, trapped_H2, trapped_HD, trapped_D2], - name="empty_trap", -) - - -my_model.species = [ - mobile_H, - mobile_D, - trapped_H1, - trapped_D1, - trapped_H2, - trapped_HD, - trapped_D2, -] - -my_model.reactions = [ - F.Reaction( - k_0=4.1e-7 / (1.1e-10**2 * 6 * w_atom_density), - E_k=0.39, - p_0=1e13, - E_p=1.2, - reactant=[mobile_H, empty_trap], - product=trapped_H1, - ), - F.Reaction( - k_0=4.1e-7 / (1.1e-10**2 * 6 * w_atom_density), - E_k=0.39, - p_0=1e13, - E_p=1.0, - reactant=[mobile_H, trapped_H1], - product=trapped_H2, - ), - F.Reaction( - k_0=4.1e-7 / (1.1e-10**2 * 6 * w_atom_density), - E_k=0.39, - p_0=1e13, - E_p=1.2, - reactant=[mobile_D, empty_trap], - product=trapped_D1, - ), - F.Reaction( - k_0=4.1e-7 / (1.1e-10**2 * 6 * w_atom_density), - E_k=0.39, - p_0=1e13, - E_p=1.2, - reactant=[mobile_D, trapped_D1], - product=trapped_D2, - ), - F.Reaction( - k_0=4.1e-7 / (1.1e-10**2 * 6 * w_atom_density), - E_k=0.39, - p_0=1e13, - E_p=1.0, - reactant=[mobile_H, trapped_D1], - product=trapped_HD, - ), - F.Reaction( - k_0=4.1e-7 / (1.1e-10**2 * 6 * w_atom_density), - E_k=0.39, - p_0=1e13, - E_p=1.0, - reactant=[mobile_D, trapped_H1], - product=trapped_HD, - ), -] - -# -------- Temperature --------- # - -my_model.temperature = 300 - -# -------- Boundary conditions --------- # - - -my_model.boundary_conditions = [ - F.DirichletBC(subdomain=left_surface, value=1e20, species=mobile_H), - F.DirichletBC(subdomain=right_surface, value=1e19, species=mobile_D), - F.DirichletBC(subdomain=right_surface, value=0, species=mobile_H), - F.DirichletBC(subdomain=left_surface, value=0, species=mobile_D), -] - -# -------- Exports --------- # - -left_flux = F.SurfaceFlux(field=mobile_H, surface=left_surface) -right_flux = F.SurfaceFlux(field=mobile_H, surface=right_surface) - -folder = "multi_isotope_trapping_example" - -my_model.exports = [ - F.XDMFExport(f"{folder}/mobile_concentration_h.xdmf", field=mobile_H), - F.XDMFExport(f"{folder}/mobile_concentration_d.xdmf", field=mobile_D), - F.XDMFExport(f"{folder}/trapped_concentration_h1.xdmf", field=trapped_H1), - F.XDMFExport(f"{folder}/trapped_concentration_h2.xdmf", field=trapped_H2), - F.XDMFExport(f"{folder}/trapped_concentration_d1.xdmf", field=trapped_D1), - F.XDMFExport(f"{folder}/trapped_concentration_d2.xdmf", field=trapped_D2), - F.XDMFExport(f"{folder}/trapped_concentration_hd.xdmf", field=trapped_HD), -] - -# -------- Settings --------- # - -my_model.settings = F.Settings( - atol=1e-10, rtol=1e-10, max_iterations=30, final_time=3000 -) - -my_model.settings.stepsize = F.Stepsize(initial_value=20) - -# -------- Run --------- # - -my_model.initialise() - -print(my_model.formulation) -# exit() -my_model.run() diff --git a/examples/multi_material_1d.py b/examples/multi_material_1d.py deleted file mode 100644 index 5c515d0b3..000000000 --- a/examples/multi_material_1d.py +++ /dev/null @@ -1,99 +0,0 @@ -import numpy as np - -import festim as F - -my_model = F.HydrogenTransportProblemDiscontinuous() - -interface_1 = 0.5 -interface_2 = 0.7 - -# for some reason if the mesh isn't fine enough then I have a random SEGV error -N = 2000 -vertices = np.concatenate( - [ - np.linspace(0, interface_1, num=N), - np.linspace(interface_1, interface_2, num=N), - np.linspace(interface_2, 1, num=N), - ] -) - -my_model.mesh = F.Mesh1D(vertices) - -material_left = F.Material(D_0=2.0, E_D=0) -material_mid = F.Material(D_0=2.0, E_D=0) -material_right = F.Material(D_0=2.0, E_D=0) - -material_left.K_S_0 = 2.0 -material_left.E_K_S = 0 -material_mid.K_S_0 = 4.0 -material_mid.E_K_S = 0 -material_right.K_S_0 = 6.0 -material_right.E_K_S = 0 - -left_domain = F.VolumeSubdomain1D( - 3, borders=[vertices[0], interface_1], material=material_left -) -middle_domain = F.VolumeSubdomain1D( - 4, borders=[interface_1, interface_2], material=material_mid -) -right_domain = F.VolumeSubdomain1D( - 5, borders=[interface_2, vertices[-1]], material=material_right -) - -left_surface = F.SurfaceSubdomain1D(id=1, x=vertices[0]) -right_surface = F.SurfaceSubdomain1D(id=2, x=vertices[-1]) - -# the ids here are arbitrary in 1D, you can put anything as long as it's -# not the same as the surfaces -my_model.interfaces = [ - F.Interface(6, (left_domain, middle_domain)), - F.Interface(7, (middle_domain, right_domain)), -] -my_model.subdomains = [ - left_domain, - middle_domain, - right_domain, - left_surface, - right_surface, -] - -H = F.Species("H", mobile=True) -trapped_H = F.Species("H_trapped", mobile=False) -empty_trap = F.ImplicitSpecies(n=0.5, others=[trapped_H]) - -my_model.species = [H, trapped_H] - -for species in my_model.species: - species.subdomains = my_model.volume_subdomains - - -my_model.reactions = [ - F.Reaction( - reactant=[H, empty_trap], - product=[trapped_H], - k_0=2, - E_k=0, - p_0=0.1, - E_p=0, - volume=domain, - ) - for domain in [left_domain, middle_domain, right_domain] -] - -my_model.boundary_conditions = [ - F.DirichletBC(left_surface, value=0.05, species=H), - F.DirichletBC(right_surface, value=0.2, species=H), -] - - -my_model.temperature = lambda x: 300 + 100 * x[0] - -my_model.settings = F.Settings(atol=None, rtol=1e-5, transient=False) - -my_model.exports = [ - F.VTXSpeciesExport(filename=f"u_{subdomain.id}.bp", field=H, subdomain=subdomain) - for subdomain in my_model.volume_subdomains -] - -my_model.initialise() -my_model.run() diff --git a/examples/multi_material_2d.py b/examples/multi_material_2d.py deleted file mode 100644 index 2fed1f6a9..000000000 --- a/examples/multi_material_2d.py +++ /dev/null @@ -1,202 +0,0 @@ -from mpi4py import MPI - -import dolfinx -import dolfinx.fem.petsc -import numpy as np -import ufl - -import festim as F - - -# ---------------- Generate a mesh ---------------- -def generate_mesh(): - def bottom_boundary(x): - return np.isclose(x[1], 0.0) - - def top_boundary(x): - return np.isclose(x[1], 1.0) - - def half(x): - return x[1] <= 0.5 + 1e-14 - - mesh = dolfinx.mesh.create_unit_square( - MPI.COMM_WORLD, 20, 20, dolfinx.mesh.CellType.triangle - ) - - # Split domain in half and set an interface tag of 5 - tdim = mesh.topology.dim - fdim = tdim - 1 - top_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, top_boundary) - bottom_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, bottom_boundary) - num_facets_local = ( - mesh.topology.index_map(fdim).size_local - + mesh.topology.index_map(fdim).num_ghosts - ) - facets = np.arange(num_facets_local, dtype=np.int32) - values = np.full_like(facets, 0, dtype=np.int32) - values[top_facets] = 1 - values[bottom_facets] = 2 - - bottom_cells = dolfinx.mesh.locate_entities(mesh, tdim, half) - num_cells_local = ( - mesh.topology.index_map(tdim).size_local - + mesh.topology.index_map(tdim).num_ghosts - ) - cells = np.full(num_cells_local, 4, dtype=np.int32) - cells[bottom_cells] = 3 - ct = dolfinx.mesh.meshtags( - mesh, tdim, np.arange(num_cells_local, dtype=np.int32), cells - ) - all_b_facets = dolfinx.mesh.compute_incident_entities( - mesh.topology, ct.find(3), tdim, fdim - ) - all_t_facets = dolfinx.mesh.compute_incident_entities( - mesh.topology, ct.find(4), tdim, fdim - ) - interface = np.intersect1d(all_b_facets, all_t_facets) - values[interface] = 5 - - mt = dolfinx.mesh.meshtags(mesh, mesh.topology.dim - 1, facets, values) - return mesh, mt, ct - - -mesh, mt, ct = generate_mesh() - -my_model = F.HydrogenTransportProblemDiscontinuous() -my_model.mesh = F.Mesh(mesh) -my_model.volume_meshtags = ct -my_model.facet_meshtags = mt - -material_bottom = F.Material(D_0=1, E_D=0 * 0.1) -material_top = F.Material(D_0=1, E_D=0 * 0.1) - -material_bottom.K_S_0 = 2.0 -material_bottom.E_K_S = 0 * 0.1 -material_top.K_S_0 = 4.0 -material_top.E_K_S = 0 * 0.12 - -top_domain = F.VolumeSubdomain(4, material=material_top) -bottom_domain = F.VolumeSubdomain(3, material=material_bottom) - -top_surface = F.SurfaceSubdomain(id=1) -bottom_surface = F.SurfaceSubdomain(id=2) -my_model.subdomains = [bottom_domain, top_domain, top_surface, bottom_surface] - -# we should be able to automate this -my_model.interfaces = [F.Interface(5, (bottom_domain, top_domain))] - -H = F.Species("H", mobile=True) - -my_model.species = [H] - -for species in my_model.species: - species.subdomains = [bottom_domain, top_domain] - - -my_model.boundary_conditions = [ - F.DirichletBC(top_surface, value=0.05, species=H), - # F.DirichletBC(bottom_surface, value=1, species=H), - F.ParticleFluxBC(bottom_surface, value=lambda x: 1.0 + x[0], species=H), -] - - -my_model.temperature = 500.0 # lambda x: 300 + 10 * x[1] + 100 * x[0] - -my_model.settings = F.Settings(atol=1e-5, rtol=1e-5, transient=False) - -my_model.exports = [ - F.VTXSpeciesExport(f"u_{subdomain.id}.bp", field=H, subdomain=subdomain) - for subdomain in my_model.volume_subdomains -] - -my_model.initialise() -my_model.run() - -# -------------------- post processing -------------------- - -# derived quantities -# Check for DOLFINx version compatibility -try: - from dolfinx.mesh import EntityMap # noqa: F401 - - legacy_entity_map = False - - def entity_map_wrapper(e_map): - return list(e_map.values()) - - entity_maps = [sd.cell_map for sd in my_model.volume_subdomains] -except ImportError: - legacy_entity_map = True - - def entity_map_wrapper(e_map): - return e_map - - entity_maps = { - sd.submesh: sd.parent_to_submesh for sd in my_model.volume_subdomains - } - - -ds_b = ufl.Measure("ds", domain=bottom_domain.submesh, subdomain_data=bottom_domain.ft) -ds_t = ufl.Measure("ds", domain=top_domain.submesh, subdomain_data=top_domain.ft) -dx_b = ufl.Measure("dx", domain=bottom_domain.submesh) -dx = ufl.Measure("dx", domain=mesh) -ds = ufl.Measure("ds", domain=mesh, subdomain_data=my_model.facet_meshtags) - -n = ufl.FacetNormal(mesh) -n_b = ufl.FacetNormal(bottom_domain.submesh) -n_t = ufl.FacetNormal(top_domain.submesh) - -D = bottom_domain.material.get_diffusion_coefficient( - my_model.mesh.mesh, my_model.temperature_fenics, H -) -form = dolfinx.fem.form( - ufl.dot( - D * ufl.grad(bottom_domain.u.sub(0)), - n, - ) - * ds(1), - entity_maps=entity_maps, -) - - -def assemble_and_print(form: dolfinx.fem.Form): - loc_val = dolfinx.fem.assemble_scalar(form) - glob_val = mesh.comm.allreduce(loc_val, op=MPI.SUM) - if MPI.COMM_WORLD.rank == 0: - print(glob_val, flush=True) - - -assemble_and_print(form) - - -form = dolfinx.fem.form(bottom_domain.u.sub(0) * dx_b) -assemble_and_print(form) - -form = dolfinx.fem.form( - my_model.temperature_fenics * dx_b, - entity_maps=entity_map_wrapper({mesh: bottom_domain.cell_map}), -) -assemble_and_print(form) - -D = bottom_domain.material.get_diffusion_coefficient( - my_model.mesh.mesh, my_model.temperature_fenics, H -) -id_interface = 5 -form = dolfinx.fem.form( - ufl.dot( - D * ufl.grad(bottom_domain.u.sub(0)), - n_b, - ) - * ds_b(id_interface), - entity_maps=entity_map_wrapper({mesh: bottom_domain.cell_map}), -) -assemble_and_print(form) -form = dolfinx.fem.form( - ufl.dot( - D * ufl.grad(top_domain.u.sub(0)), - n_t, - ) - * ds_t(id_interface), - entity_maps=entity_map_wrapper({mesh: top_domain.cell_map}), -) -assemble_and_print(form) diff --git a/examples/multi_material_2d_bis.py b/examples/multi_material_2d_bis.py deleted file mode 100644 index 0b7841ebd..000000000 --- a/examples/multi_material_2d_bis.py +++ /dev/null @@ -1,224 +0,0 @@ -from mpi4py import MPI - -import dolfinx -import dolfinx.fem.petsc -import numpy as np -import ufl - -import festim as F - - -# ---------------- Generate a mesh ---------------- -def generate_mesh(): - def bottom_boundary(x): - return np.isclose(x[1], 0.0) - - def top_boundary(x): - return np.isclose(x[1], 1.0) - - def half(x): - return x[1] <= 0.5 + 1e-14 - - mesh = dolfinx.mesh.create_unit_square( - MPI.COMM_WORLD, 20, 20, dolfinx.mesh.CellType.triangle - ) - - # Split domain in half and set an interface tag of 5 - tdim = mesh.topology.dim - fdim = tdim - 1 - top_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, top_boundary) - bottom_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, bottom_boundary) - num_facets_local = ( - mesh.topology.index_map(fdim).size_local - + mesh.topology.index_map(fdim).num_ghosts - ) - facets = np.arange(num_facets_local, dtype=np.int32) - values = np.full_like(facets, 0, dtype=np.int32) - values[top_facets] = 1 - values[bottom_facets] = 2 - - bottom_cells = dolfinx.mesh.locate_entities(mesh, tdim, half) - num_cells_local = ( - mesh.topology.index_map(tdim).size_local - + mesh.topology.index_map(tdim).num_ghosts - ) - cells = np.full(num_cells_local, 4, dtype=np.int32) - cells[bottom_cells] = 3 - ct = dolfinx.mesh.meshtags( - mesh, tdim, np.arange(num_cells_local, dtype=np.int32), cells - ) - all_b_facets = dolfinx.mesh.compute_incident_entities( - mesh.topology, ct.find(3), tdim, fdim - ) - all_t_facets = dolfinx.mesh.compute_incident_entities( - mesh.topology, ct.find(4), tdim, fdim - ) - interface = np.intersect1d(all_b_facets, all_t_facets) - values[interface] = 5 - - mt = dolfinx.mesh.meshtags(mesh, mesh.topology.dim - 1, facets, values) - return mesh, mt, ct - - -mesh, mt, ct = generate_mesh() - -my_model = F.HydrogenTransportProblemDiscontinuous() -my_model.mesh = F.Mesh(mesh) -my_model.volume_meshtags = ct -my_model.facet_meshtags = mt - -material_bottom = F.Material(D_0=2.0, E_D=0.1) -material_top = F.Material(D_0=2.0, E_D=0.1) - -material_bottom.K_S_0 = 2.0 -material_bottom.E_K_S = 0 -material_top.K_S_0 = 4.0 -material_top.E_K_S = 0 - -top_domain = F.VolumeSubdomain(4, material=material_top) -bottom_domain = F.VolumeSubdomain(3, material=material_bottom) - -top_surface = F.SurfaceSubdomain(id=1) -bottom_surface = F.SurfaceSubdomain(id=2) -my_model.subdomains = [bottom_domain, top_domain, top_surface, bottom_surface] - -# we should be able to automate this -my_model.interfaces = [F.Interface(5, (bottom_domain, top_domain))] - -H = F.Species("H", mobile=True) -trapped_H = F.Species("H_trapped", mobile=False) -empty_trap = F.ImplicitSpecies(n=0.5, others=[trapped_H]) - -my_model.species = [H, trapped_H] - -for species in my_model.species: - species.subdomains = [bottom_domain, top_domain] - - -my_model.reactions = [ - F.Reaction( - reactant=[H, empty_trap], - product=[trapped_H], - k_0=2, - E_k=0, - p_0=0.1, - E_p=0, - volume=top_domain, - ), - F.Reaction( - reactant=[H, empty_trap], - product=[trapped_H], - k_0=2, - E_k=0, - p_0=0.1, - E_p=0, - volume=bottom_domain, - ), -] - -my_model.boundary_conditions = [ - F.FixedConcentrationBC(top_surface, value=0.05, species=H), - F.FixedConcentrationBC(bottom_surface, value=0.2, species=H), -] - - -my_model.temperature = lambda x: 400 + 10 * x[1] + 100 * x[0] - -my_model.settings = F.Settings(atol=1e-10, rtol=1e-5, transient=False) - -my_model.exports = [ - F.VTXSpeciesExport(f"c_{subdomain.id}.bp", field=H, subdomain=subdomain) - for subdomain in my_model.volume_subdomains -] + [ - F.VTXSpeciesExport(f"c_t_{subdomain.id}.bp", field=trapped_H, subdomain=subdomain) - for subdomain in my_model.volume_subdomains -] - -my_model.initialise() -my_model.run() - -# -------------------- post processing -------------------- - - -# derived quantities - - -ds = ufl.Measure("ds", domain=mesh, subdomain_data=my_model.facet_meshtags) -ds_b = ufl.Measure("ds", domain=bottom_domain.submesh, subdomain_data=bottom_domain.ft) -ds_t = ufl.Measure("ds", domain=top_domain.submesh, subdomain_data=top_domain.ft) -dx_b = ufl.Measure("dx", domain=bottom_domain.submesh) -dx = ufl.Measure("dx", domain=mesh) - -n_b = ufl.FacetNormal(bottom_domain.submesh) -n_t = ufl.FacetNormal(top_domain.submesh) - -form = dolfinx.fem.form(bottom_domain.u.sub(0) * dx_b) -print(dolfinx.fem.assemble_scalar(form)) - -form = dolfinx.fem.form(bottom_domain.u.sub(1) * dx_b) -print(dolfinx.fem.assemble_scalar(form)) - -# Check for DOLFINx version compatibility -try: - from dolfinx.mesh import EntityMap # noqa: F401 - - legacy_entity_map = False - - def entity_map_wrapper(e_map): - return list(e_map.values()) -except ImportError: - legacy_entity_map = True - - def entity_map_wrapper(e_map): - return e_map - - -form = dolfinx.fem.form( - my_model.temperature_fenics * dx_b, - entity_maps=entity_map_wrapper({mesh: bottom_domain.cell_map}), -) -print(dolfinx.fem.assemble_scalar(form)) - -D = bottom_domain.material.get_diffusion_coefficient( - my_model.mesh.mesh, my_model.temperature_fenics, H -) -id_interface = 5 - -form = dolfinx.fem.form( - ufl.dot( - D * ufl.grad(bottom_domain.u.sub(0)), - n_b, - ) - * ds_b(id_interface), - entity_maps=entity_map_wrapper({mesh: bottom_domain.cell_map}), -) -print(dolfinx.fem.assemble_scalar(form)) -form = dolfinx.fem.form( - ufl.dot( - D * ufl.grad(top_domain.u.sub(0)), - n_t, - ) - * ds_t(id_interface), - entity_maps=entity_map_wrapper({mesh: top_domain.cell_map}), -) -print(dolfinx.fem.assemble_scalar(form)) - - -# using submesh function -u = H.subdomain_to_post_processing_solution[bottom_domain] -if legacy_entity_map: - entity_maps = { - sd.submesh: sd.parent_to_submesh for sd in my_model.volume_subdomains - } -else: - entity_maps = [sd.cell_map for sd in my_model.volume_subdomains] - -form = dolfinx.fem.form( - ufl.dot( - D * ufl.grad(u), - n_b, - ) - * ds(bottom_surface.id), - entity_maps=entity_maps, -) -print(dolfinx.fem.assemble_scalar(form)) diff --git a/examples/multi_material_transient.py b/examples/multi_material_transient.py deleted file mode 100644 index 00ec90271..000000000 --- a/examples/multi_material_transient.py +++ /dev/null @@ -1,99 +0,0 @@ -import numpy as np - -import festim as F - -my_model = F.HydrogenTransportProblemDiscontinuous() - -interface_1 = 0.5 -interface_2 = 0.7 - -# for some reason if the mesh isn't fine enough then I have a random SEGV error -N = 1500 -vertices = np.concatenate( - [ - np.linspace(0, interface_1, num=N), - np.linspace(interface_1, interface_2, num=N), - np.linspace(interface_2, 1, num=N), - ] -) - -my_model.mesh = F.Mesh1D(vertices) - -material_left = F.Material(D_0=2.0, E_D=0, K_S_0=2.0, E_K_S=0) -material_mid = F.Material(D_0=2.0, E_D=0, K_S_0=4.0, E_K_S=0) -material_right = F.Material(D_0=2.0, E_D=0, K_S_0=6.0, E_K_S=0) - -left_domain = F.VolumeSubdomain1D( - 3, borders=[vertices[0], interface_1], material=material_left -) -middle_domain = F.VolumeSubdomain1D( - 4, borders=[interface_1, interface_2], material=material_mid -) -right_domain = F.VolumeSubdomain1D( - 5, borders=[interface_2, vertices[-1]], material=material_right -) - -left_surface = F.SurfaceSubdomain1D(id=1, x=vertices[0]) -right_surface = F.SurfaceSubdomain1D(id=2, x=vertices[-1]) - -# the ids here are arbitrary in 1D, you can put anything as long as it's not the same -# as the surfaces -my_model.interfaces = [ - F.Interface(6, (left_domain, middle_domain)), - F.Interface(7, (middle_domain, right_domain)), -] - -my_model.subdomains = [ - left_domain, - middle_domain, - right_domain, - left_surface, - right_surface, -] - -H = F.Species("H", mobile=True) -trapped_H = F.Species("H_trapped", mobile=False) -empty_trap = F.ImplicitSpecies(n=0.5, others=[trapped_H]) - -my_model.species = [H, trapped_H] - -for species in my_model.species: - species.subdomains = my_model.volume_subdomains - - -my_model.reactions = [ - F.Reaction( - reactant=[H, empty_trap], - product=[trapped_H], - k_0=2, - E_k=0, - p_0=0.1, - E_p=0, - volume=domain, - ) - for domain in [left_domain, middle_domain, right_domain] -] - -my_model.boundary_conditions = [ - F.DirichletBC(left_surface, value=0.05, species=H), - F.DirichletBC(right_surface, value=0.2, species=H), -] - - -my_model.temperature = lambda x: 300 + 100 * x[0] - -my_model.settings = F.Settings(atol=None, rtol=1e-5, transient=True, final_time=100) -my_model.settings.stepsize = 1 - -my_model.exports = [ - F.VTXSpeciesExport(filename=f"u_{subdomain.id}.bp", field=H, subdomain=subdomain) - for subdomain in my_model.volume_subdomains -] + [ - F.VTXSpeciesExport( - filename=f"u_t_{subdomain.id}.bp", field=trapped_H, subdomain=subdomain - ) - for subdomain in my_model.volume_subdomains -] - -my_model.initialise() -my_model.run() diff --git a/examples/multi_material_with_one_material.py b/examples/multi_material_with_one_material.py deleted file mode 100644 index d58f3ffc2..000000000 --- a/examples/multi_material_with_one_material.py +++ /dev/null @@ -1,65 +0,0 @@ -import numpy as np - -import festim as F - -my_model = F.HydrogenTransportProblemDiscontinuous() - - -N = 1500 -vertices = np.linspace(0, 1, num=N) -my_model.mesh = F.Mesh1D(vertices) - -material = F.Material(D_0=2.0, E_D=0) - -material.K_S_0 = 2.0 -material.E_K_S = 0 - -subdomain = F.VolumeSubdomain1D( - 1, borders=[vertices[0], vertices[-1]], material=material -) - -left_surface = F.SurfaceSubdomain1D(id=1, x=vertices[0]) -right_surface = F.SurfaceSubdomain1D(id=2, x=vertices[-1]) - -my_model.subdomains = [subdomain, left_surface, right_surface] - -H = F.Species("H", mobile=True) -trapped_H = F.Species("H_trapped", mobile=False) -empty_trap = F.ImplicitSpecies(n=0.5, others=[trapped_H]) - -my_model.species = [H, trapped_H] - -for species in my_model.species: - species.subdomains = my_model.volume_subdomains - - -my_model.reactions = [ - F.Reaction( - reactant=[H, empty_trap], - product=[trapped_H], - k_0=2, - E_k=0, - p_0=0.1, - E_p=0, - volume=domain, - ) - for domain in my_model.volume_subdomains -] - -my_model.boundary_conditions = [ - F.DirichletBC(left_surface, value=0.05, species=H), - F.DirichletBC(right_surface, value=0.2, species=H), -] - - -my_model.temperature = lambda x: 300 + 100 * x[0] - -my_model.settings = F.Settings(atol=None, rtol=1e-5, transient=False) - -my_model.exports = [ - F.VTXSpeciesExport(filename=f"u_{subdomain.id}.bp", field=H, subdomain=subdomain) - for subdomain in my_model.volume_subdomains -] - -my_model.initialise() -my_model.run() diff --git a/examples/plot_tds.py b/examples/plot_tds.py deleted file mode 100644 index ead6811fe..000000000 --- a/examples/plot_tds.py +++ /dev/null @@ -1,19 +0,0 @@ -import matplotlib.pyplot as plt -import numpy as np - -times = np.genfromtxt("times_tds.txt") -flux = np.genfromtxt("outgassing_flux_tds.txt") -flux = np.abs(flux) - -plt.figure() -plt.plot(times, flux) -plt.xlim(450, 500) -# plt.ylim(top=1e18) -plt.ylim(bottom=0, top=1e19) -# plt.ylim(bottom=1.25e18, top=0.6e19) -plt.ylabel(r"Desorption flux (m$^{-2}$ s$^{-1}$)") -plt.xlabel(r"Time (s)") - -plt.ylabel(r"Desorption flux (m$^{-2}$ s$^{-1}$)") -plt.xlabel(r"Time (s)") -plt.show() diff --git a/examples/surface_reaction_example.py b/examples/surface_reaction_example.py deleted file mode 100644 index 14d41ea80..000000000 --- a/examples/surface_reaction_example.py +++ /dev/null @@ -1,169 +0,0 @@ -import dolfinx.fem as fem -import matplotlib.pyplot as plt -import numpy as np -import ufl - -import festim as F - - -class FluxFromSurfaceReaction(F.SurfaceFlux): - def __init__(self, reaction: F.SurfaceReactionBC): - super().__init__( - F.Species(), # just a dummy species here - reaction.subdomain, - ) - self.reaction = reaction.flux_bcs[0] - - def compute(self, ds): - self.value = fem.assemble_scalar( - fem.form(self.reaction.value_fenics * ds(self.surface.id)) - ) - self.data.append(self.value) - - -my_model = F.HydrogenTransportProblem() -my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, 1000)) -my_mat = F.Material(name="mat", D_0=1, E_D=0) -vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=my_mat) -left = F.SurfaceSubdomain1D(id=1, x=0) -right = F.SurfaceSubdomain1D(id=2, x=1) - -my_model.subdomains = [vol, left, right] - -H = F.Species("H") -D = F.Species("D") -my_model.species = [H, D] - -my_model.temperature = 500 - -surface_reaction_hd = F.SurfaceReactionBC( - reactant=[H, D], - gas_pressure=0, - k_r0=0.01, - E_kr=0, - k_d0=0, - E_kd=0, - subdomain=right, -) - -surface_reaction_hh = F.SurfaceReactionBC( - reactant=[H, H], - gas_pressure=lambda t: ufl.conditional(ufl.gt(t, 1), 2, 0), - k_r0=0.02, - E_kr=0, - k_d0=0.03, - E_kd=0, - subdomain=right, -) - -surface_reaction_dd = F.SurfaceReactionBC( - reactant=[D, D], - gas_pressure=0, - k_r0=0.01, - E_kr=0, - k_d0=0, - E_kd=0, - subdomain=right, -) - -my_model.boundary_conditions = [ - F.DirichletBC(subdomain=left, value=2, species=H), - F.DirichletBC(subdomain=left, value=2, species=D), - surface_reaction_hd, - surface_reaction_hh, - surface_reaction_dd, -] - -H_flux_right = F.SurfaceFlux(H, right) -H_flux_left = F.SurfaceFlux(H, left) -D_flux_right = F.SurfaceFlux(D, right) -D_flux_left = F.SurfaceFlux(D, left) -HD_flux = FluxFromSurfaceReaction(surface_reaction_hd) -HH_flux = FluxFromSurfaceReaction(surface_reaction_hh) -DD_flux = FluxFromSurfaceReaction(surface_reaction_dd) -my_model.exports = [ - F.XDMFExport("test.xdmf", H), - H_flux_left, - H_flux_right, - D_flux_left, - D_flux_right, - HD_flux, - HH_flux, - DD_flux, -] - -my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=5, transient=True) - -my_model.settings.stepsize = 0.1 - -my_model.initialise() -my_model.run() - - -plt.stackplot( - H_flux_left.t, - np.abs(H_flux_left.data), - np.abs(D_flux_left.data), - labels=["H_in", "D_in"], -) -plt.stackplot( - H_flux_right.t, - -np.abs(H_flux_right.data), - -np.abs(D_flux_right.data), - labels=["H_out", "D_out"], -) -plt.legend() -plt.xlabel("Time (s)") -plt.ylabel("Flux (atom/m^2/s)") -plt.figure() -plt.stackplot( - HD_flux.t, - np.abs(HH_flux.data), - np.abs(HD_flux.data), - np.abs(DD_flux.data), - labels=["HH", "HD", "DD"], -) -plt.legend(reverse=True) -plt.xlabel("Time (s)") -plt.ylabel("Flux (molecule/m^2/s)") - - -plt.figure() -plt.plot(H_flux_right.t, -np.array(H_flux_right.data), label="from gradient (H)") -plt.plot( - H_flux_right.t, - 2 * np.array(HH_flux.data) + np.array(HD_flux.data), - linestyle="--", - label="from reaction rates (H)", -) - -plt.plot(D_flux_right.t, -np.array(D_flux_right.data), label="from gradient (D)") -plt.plot( - D_flux_right.t, - 2 * np.array(DD_flux.data) + np.array(HD_flux.data), - linestyle="--", - label="from reaction rates (D)", -) -plt.xlabel("Time (s)") -plt.ylabel("Flux (atom/m^2/s)") -plt.legend() -plt.show() - -# check that H_flux_right == 2*HH_flux + HD_flux -H_flux_from_gradient = -np.array(H_flux_right.data) -H_flux_from_reac = 2 * np.array(HH_flux.data) + np.array(HD_flux.data) -assert np.allclose( - H_flux_from_gradient, - H_flux_from_reac, - rtol=0.5e-2, - atol=0.005, -) -# check that D_flux_right == 2*DD_flux + HD_flux -D_flux_from_gradient = -np.array(D_flux_right.data) -D_flux_from_reac = 2 * np.array(DD_flux.data) + np.array(HD_flux.data) -assert np.allclose( - D_flux_from_gradient, - D_flux_from_reac, - rtol=0.5e-2, - atol=0.005, -) diff --git a/examples/tds_example.py b/examples/tds_example.py deleted file mode 100644 index 7f7d7bcbb..000000000 --- a/examples/tds_example.py +++ /dev/null @@ -1,136 +0,0 @@ -import numpy as np - -import festim as F - -my_model = F.HydrogenTransportProblem() - -# -------- Mesh --------- # - -L = 20e-6 -vertices = np.concatenate( - [ - np.linspace(0, 30e-9, num=200, endpoint=False), - np.linspace(30e-9, 3e-6, num=300, endpoint=False), - np.linspace(3e-6, L, num=200), - ] -) -my_model.mesh = F.Mesh1D(vertices) - - -# -------- Materials and subdomains --------- # - -w_atom_density = 6.3e28 # atom/m3 - -tungsten = F.Material(D_0=4.1e-7, E_D=0.39, name="tungsten") - -my_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, L], material=tungsten) -left_surface = F.SurfaceSubdomain1D(id=1, x=0) -right_surface = F.SurfaceSubdomain1D(id=2, x=L) - -my_model.subdomains = [ - my_subdomain, - left_surface, - right_surface, -] - -# -------- Hydrogen species and reactions --------- # - -mobile_H = F.Species("H") -trapped_H1 = F.Species("trapped_H1", mobile=False) -empty_trap1 = F.ImplicitSpecies( - n=1.3e-3 * w_atom_density, others=[trapped_H1], name="empty_trap1" -) -trapped_H2 = F.Species("trapped_H2", mobile=False) -empty_trap2 = F.ImplicitSpecies( - n=4e-3 * w_atom_density, others=[trapped_H2], name="empty_trap2" -) -my_model.species = [mobile_H, trapped_H1, trapped_H2] - -my_model.reactions = [ - F.Reaction( - k_0=4.1e-7 / (1.1e-10**2 * 6 * w_atom_density), - E_k=0.39, - p_0=1e13, - E_p=0.87, - reactant=[mobile_H, empty_trap1], - product=trapped_H1, - volume=my_subdomain, - ), - F.Reaction( - k_0=4.1e-7 / (1.1e-10**2 * 6 * w_atom_density), - E_k=0.39, - p_0=1e13, - E_p=1.0, - reactant=[mobile_H, empty_trap2], - product=trapped_H2, - volume=my_subdomain, - ), -] - -# -------- Temperature --------- # - -implantation_time = 400 -start_tds = implantation_time + 50 -implantation_temp = 300 -temperature_ramp = 8 # K/s - - -def temp_function(t): - if t < start_tds: - return implantation_temp - else: - return implantation_temp + temperature_ramp * (t - start_tds) - - -my_model.temperature = temp_function - -# -------- Boundary conditions --------- # - - -def left_conc_value(t): - if t < implantation_time: - _T = temp_function(t) - D = tungsten.D_0 * np.exp(-tungsten.E_D / F.k_B / _T) - return 2.5e19 * 4.5e-9 / D - else: - return 0.0 - - -my_model.boundary_conditions = [ - F.DirichletBC(subdomain=left_surface, value=left_conc_value, species=mobile_H), - F.DirichletBC(subdomain=right_surface, value=0, species=mobile_H), -] - -# -------- Exports --------- # - -left_flux = F.SurfaceFlux(field=mobile_H, surface=left_surface) -right_flux = F.SurfaceFlux(field=mobile_H, surface=right_surface) - -my_model.exports = [ - # F.VTXSpeciesExport("mobile_concentration_h.bp", field=mobile_H), - # F.VTXSpeciesExport("trapped_concentration_h.bp", field=trapped_H1), - F.XDMFExport("mobile_concentration_h.xdmf", field=mobile_H), - F.XDMFExport("trapped_concentration_h1.xdmf", field=trapped_H1), - F.XDMFExport("trapped_concentration_h2.xdmf", field=trapped_H2), - left_flux, - right_flux, -] - -# -------- Settings --------- # - -my_model.settings = F.Settings(atol=1e10, rtol=1e-10, max_iterations=30, final_time=500) - -my_model.settings.stepsize = F.Stepsize(initial_value=0.5) - -# -------- Run --------- # - -my_model.initialise() -my_model.run() - -# -------- Save results --------- # - -np.savetxt( - "outgassing_flux_tds.txt", - np.array(left_flux.data) + np.array(right_flux.data), -) -np.savetxt("times_tds.txt", np.array(left_flux.t)) diff --git a/examples/trapping_example.py b/examples/trapping_example.py deleted file mode 100644 index 57d339907..000000000 --- a/examples/trapping_example.py +++ /dev/null @@ -1,80 +0,0 @@ -from petsc4py import PETSc - -import numpy as np - -import festim as F - -L = 3e-04 -vertices = np.linspace(0, L, num=1000) - -my_mesh = F.Mesh1D(vertices) - -my_model = F.HydrogenTransportProblem() -my_model.mesh = my_mesh - -my_mat = F.Material(D_0=1.9e-7, E_D=0.2, name="my_mat") -my_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, L], material=my_mat) -left_surface = F.SurfaceSubdomain1D(id=1, x=0) -right_surface = F.SurfaceSubdomain1D(id=2, x=L) -my_model.subdomains = [ - my_subdomain, - left_surface, - right_surface, -] - -mobile_H = F.Species("H") -trapped_H = F.Species("trapped_H", mobile=False) -empty_trap = F.ImplicitSpecies(n=1e19, others=[trapped_H], name="empty_trap") -my_model.species = [mobile_H, trapped_H] - -my_model.reactions = [ - F.Reaction( - p_0=1e13, - E_p=1.2, - k_0=3.8e-17, - E_k=0.2, - reactant=[mobile_H, empty_trap], - product=trapped_H, - ) -] - -temperature = 500.0 -my_model.temperature = temperature - -my_model.boundary_conditions = [ - F.DirichletBC(subdomain=right_surface, value=0, species="H"), - F.DirichletBC(subdomain=left_surface, value=1e12, species=mobile_H), -] -my_model.exports = [ - F.VTXSpeciesExport( - "mobile_concentration_h.bp", field=mobile_H - ), # produces 0 in file - F.VTXSpeciesExport( - "trapped_concentration_h.bp", field=trapped_H - ), # produces 0 in file - F.XDMFExport("mobile_concentration_h.xdmf", field=mobile_H), - F.XDMFExport("trapped_concentration_h.xdmf", field=trapped_H), -] - -my_model.settings = F.Settings( - atol=1e10, - rtol=1e-10, - max_iterations=30, - final_time=50, -) - -my_model.settings.stepsize = F.Stepsize(initial_value=1 / 20) - -my_model.initialise() - - -my_model.solver.convergence_criterion = "incremental" -ksp = my_model.solver.krylov_solver -opts = PETSc.Options() -option_prefix = ksp.getOptionsPrefix() -opts[f"{option_prefix}ksp_type"] = "cg" -opts[f"{option_prefix}pc_type"] = "gamg" -opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps" -ksp.setFromOptions() - -times, flux_values = my_model.run()