diff --git a/festim/source.py b/festim/source.py index 6005fa21b..e7c9ad91b 100644 --- a/festim/source.py +++ b/festim/source.py @@ -48,9 +48,9 @@ def volume(self): @volume.setter def volume(self, value): - # check that volume is festim.VolumeSubdomain1D - if not isinstance(value, F.VolumeSubdomain1D): - raise TypeError("volume must be of type festim.VolumeSubdomain1D") + # check that volume is festim.VolumeSubdomain + if not isinstance(value, F.VolumeSubdomain): + raise TypeError("volume must be of type festim.VolumeSubdomain") self._volume = value @property @@ -74,9 +74,11 @@ def value_fenics(self, value): if value is None: self._value_fenics = value return - if not isinstance(value, (fem.Function, fem.Constant, np.ndarray)): + if not isinstance( + value, (fem.Function, fem.Constant, np.ndarray, ufl.core.expr.Expr) + ): raise TypeError( - f"Value must be a dolfinx.fem.Function, dolfinx.fem.Constant, or a np.ndarray not {type(value)}" + f"Value must be a dolfinx.fem.Function, dolfinx.fem.Constant, np.ndarray or a ufl.core.expr.Expr, not {type(value)}" ) self._value_fenics = value @@ -125,6 +127,9 @@ def create_value_fenics( if isinstance(self.value, (int, float)): self.value_fenics = F.as_fenics_constant(mesh=mesh, value=self.value) + elif isinstance(self.value, (fem.Function, ufl.core.expr.Expr)): + self.value_fenics = self.value + elif callable(self.value): arguments = self.value.__code__.co_varnames diff --git a/test/system_tests/__init__.py b/test/system_tests/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/test/system_tests/test_1_mobile_1_trap.py b/test/system_tests/test_1_mobile_1_trap.py new file mode 100644 index 000000000..58aeb073c --- /dev/null +++ b/test/system_tests/test_1_mobile_1_trap.py @@ -0,0 +1,443 @@ +import festim as F +import numpy as np +from dolfinx import fem +import ufl +from .tools import error_L2 +from dolfinx.mesh import meshtags, create_unit_square, create_unit_cube, locate_entities +from mpi4py import MPI + + +test_mesh_1d = F.Mesh1D(np.linspace(0, 1, 10000)) +test_mesh_2d = create_unit_square(MPI.COMM_WORLD, 50, 50) +test_mesh_3d = create_unit_cube(MPI.COMM_WORLD, 20, 20, 20) +x_1d = ufl.SpatialCoordinate(test_mesh_1d.mesh) +x_2d = ufl.SpatialCoordinate(test_mesh_2d) +x_3d = ufl.SpatialCoordinate(test_mesh_3d) + + +def test_1_mobile_1_trap_MMS_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]) + + def v_exact(mod): + return lambda x: mod.sin(3 * mod.pi * x[0]) + + mobile_analytical_ufl = u_exact(ufl) + mobile_analytical_np = u_exact(np) + trapped_analytical_ufl = v_exact(ufl) + trapped_analytical_np = v_exact(np) + + elements = ufl.FiniteElement("P", test_mesh_1d.mesh.ufl_cell(), 1) + V = fem.FunctionSpace(test_mesh_1d.mesh, elements) + T = fem.Function(V) + f = fem.Function(V) + g = fem.Function(V) + + k_0 = 2 + E_k = 1.5 + p_0 = 0.2 + E_p = 0.1 + T_expr = lambda x: 500 + 100 * x[0] + T.interpolate(T_expr) + n_trap = 3 + E_D = 0.1 + D_0 = 2 + k_B = F.k_B + D = D_0 * ufl.exp(-E_D / (k_B * T)) + k = k_0 * ufl.exp(-E_k / (k_B * T)) + p = p_0 * ufl.exp(-E_p / (k_B * T)) + + f = ( + -ufl.div(D * ufl.grad(mobile_analytical_ufl(x_1d))) + + k * mobile_analytical_ufl(x_1d) * (n_trap - trapped_analytical_ufl(x_1d)) + - p * trapped_analytical_ufl(x_1d) + ) + + g = p * trapped_analytical_ufl(x_1d) - k * mobile_analytical_ufl(x_1d) * ( + n_trap - trapped_analytical_ufl(x_1d) + ) + + my_model = F.HydrogenTransportProblem() + my_model.mesh = test_mesh_1d + my_mat = F.Material(name="mat", D_0=D_0, E_D=E_D) + vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=my_mat) + left = F.SurfaceSubdomain1D(id=2, x=0) + right = F.SurfaceSubdomain1D(id=3, x=1) + my_model.subdomains = [vol, left, right] + + mobile = F.Species("mobile") + trapped = F.Species("trapped", mobile=False) + traps = F.ImplicitSpecies(n=n_trap, others=[trapped]) + my_model.species = [mobile, trapped] + + my_model.reactions = [ + F.Reaction( + reactant1=mobile, + reactant2=traps, + product=trapped, + k_0=k_0, + E_k=E_k, + p_0=p_0, + E_p=E_p, + volume=vol, + ) + ] + + my_model.temperature = T_expr + + my_model.boundary_conditions = [ + F.DirichletBC(subdomain=left, value=mobile_analytical_ufl, species=mobile), + F.DirichletBC(subdomain=right, value=mobile_analytical_ufl, species=mobile), + F.DirichletBC(subdomain=left, value=trapped_analytical_ufl, species=trapped), + F.DirichletBC(subdomain=right, value=trapped_analytical_ufl, species=trapped), + ] + + my_model.sources = [ + F.Source(value=f, volume=vol, species=mobile), + F.Source(value=g, volume=vol, species=trapped), + ] + + my_model.settings = F.Settings(atol=1e-12, rtol=1e-12, transient=False) + + my_model.initialise() + my_model.run() + + mobile_computed = mobile.post_processing_solution + trapped_computed = trapped.post_processing_solution + + L2_error_mobile = error_L2(mobile_computed, mobile_analytical_np) + L2_error_trapped = error_L2(trapped_computed, trapped_analytical_np) + + assert L2_error_mobile < 2e-07 + assert L2_error_trapped < 1e-07 + + +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 + """ + + final_time = 0.1 + + def u_exact(mod): + return lambda x, t: 1 + mod.sin(2 * mod.pi * x[0]) + 2 * t**2 + + def u_exact_alt(mod): + return lambda x: 1 + mod.sin(2 * mod.pi * x[0]) + 2 * final_time**2 + + H_analytical_ufl = u_exact(ufl) + H_analytical_np = u_exact_alt(np) + + elements = ufl.FiniteElement("P", test_mesh_1d.mesh.ufl_cell(), 1) + V = fem.FunctionSpace(test_mesh_1d.mesh, elements) + T = fem.Function(V) + + D_0 = 1 + E_D = 0.1 + T_expr = lambda x: 600 + 50 * x[0] + T.interpolate(T_expr) + D = D_0 * ufl.exp(-E_D / (F.k_B * T)) + + # FESTIM model + + my_model = F.HydrogenTransportProblem() + my_model.mesh = test_mesh_1d + my_mat = F.Material(name="mat", D_0=D_0, E_D=E_D) + 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") + my_model.species = [H] + + my_model.temperature = T_expr + + my_model.boundary_conditions = [ + F.DirichletBC(subdomain=left, value=H_analytical_ufl, species=H), + F.DirichletBC(subdomain=right, value=H_analytical_ufl, species=H), + ] + + init_value = lambda x: 1 + ufl.sin(2 * ufl.pi * x[0]) + my_model.initial_conditions = [F.InitialCondition(value=init_value, species=H)] + + f = lambda x, t: 4 * t - ufl.div(D * ufl.grad(H_analytical_ufl(x_1d, t))) + my_model.sources = [F.Source(value=f, volume=vol, species=H)] + + my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=final_time) + my_model.settings.stepsize = final_time / 50 + + my_model.initialise() + my_model.run() + + H_computed = H.post_processing_solution + + L2_error = error_L2(H_computed, H_analytical_np) + + assert L2_error < 5e-4 + + +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""" + + def u_exact(mod): + return lambda x: 1.5 + mod.sin(3 * mod.pi * x[0]) + mod.cos(3 * mod.pi * x[1]) + + def v_exact(mod): + return lambda x: mod.sin(3 * mod.pi * x[0]) + mod.cos(3 * mod.pi * x[1]) + + mobile_analytical_ufl = u_exact(ufl) + mobile_analytical_np = u_exact(np) + trapped_analytical_ufl = v_exact(ufl) + trapped_analytical_np = v_exact(np) + + elements = ufl.FiniteElement("P", test_mesh_3d.ufl_cell(), 1) + V = fem.FunctionSpace(test_mesh_3d, elements) + T = fem.Function(V) + f = fem.Function(V) + g = fem.Function(V) + + k_0 = 2 + E_k = 1.5 + p_0 = 0.2 + E_p = 0.1 + T_expr = lambda x: 500 + 100 * x[0] + T.interpolate(T_expr) + n_trap = 3 + E_D = 0.1 + D_0 = 2 + k_B = F.k_B + D = D_0 * ufl.exp(-E_D / (k_B * T)) + k = k_0 * ufl.exp(-E_k / (k_B * T)) + p = p_0 * ufl.exp(-E_p / (k_B * T)) + + f = ( + -ufl.div(D * ufl.grad(mobile_analytical_ufl(x_3d))) + + k * mobile_analytical_ufl(x_3d) * (n_trap - trapped_analytical_ufl(x_3d)) + - p * trapped_analytical_ufl(x_3d) + ) + + g = p * trapped_analytical_ufl(x_3d) - k * mobile_analytical_ufl(x_3d) * ( + n_trap - trapped_analytical_ufl(x_3d) + ) + + my_model = F.HydrogenTransportProblem() + my_model.mesh = F.Mesh(mesh=test_mesh_3d) + + # create facet meshtags + boundaries = [ + (1, lambda x: np.isclose(x[0], 0)), + (2, lambda x: np.isclose(x[0], 1)), + ] + facet_indices, facet_markers = [], [] + fdim = test_mesh_3d.topology.dim - 1 + for marker, locator in boundaries: + facets = locate_entities(test_mesh_3d, fdim, locator) + facet_indices.append(facets) + facet_markers.append(np.full_like(facets, marker)) + facet_indices = np.hstack(facet_indices).astype(np.int32) + facet_markers = np.hstack(facet_markers).astype(np.int32) + sorted_facets = np.argsort(facet_indices) + my_facet_meshtags = meshtags( + test_mesh_3d, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets] + ) + + # create volume meshtags + vdim = test_mesh_3d.topology.dim + num_cells = test_mesh_3d.topology.index_map(vdim).size_local + mesh_cell_indices = np.arange(num_cells, dtype=np.int32) + tags_volumes = np.full(num_cells, 1, dtype=np.int32) + my_volume_meshtags = meshtags(test_mesh_3d, vdim, mesh_cell_indices, tags_volumes) + + my_model.facet_meshtags = my_facet_meshtags + my_model.volume_meshtags = my_volume_meshtags + + my_mat = F.Material(name="mat", D_0=D_0, E_D=E_D) + vol = F.VolumeSubdomain(id=1, material=my_mat) + left = F.SurfaceSubdomain(id=1) + right = F.SurfaceSubdomain(id=2) + + my_model.subdomains = [vol, left, right] + + mobile = F.Species("mobile") + trapped = F.Species("trapped", mobile=False) + traps = F.ImplicitSpecies(n=n_trap, others=[trapped]) + my_model.species = [mobile, trapped] + + my_model.reactions = [ + F.Reaction( + reactant1=mobile, + reactant2=traps, + product=trapped, + k_0=k_0, + E_k=E_k, + p_0=p_0, + E_p=E_p, + volume=vol, + ) + ] + + my_model.temperature = T_expr + + my_model.boundary_conditions = [ + F.DirichletBC(subdomain=left, value=mobile_analytical_ufl, species=mobile), + F.DirichletBC(subdomain=right, value=mobile_analytical_ufl, species=mobile), + F.DirichletBC(subdomain=left, value=trapped_analytical_ufl, species=trapped), + F.DirichletBC(subdomain=right, value=trapped_analytical_ufl, species=trapped), + ] + + my_model.sources = [ + F.Source(value=f, volume=vol, species=mobile), + F.Source(value=g, volume=vol, species=trapped), + ] + + my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False) + + my_model.initialise() + my_model.run() + + mobile_computed = mobile.post_processing_solution + trapped_computed = trapped.post_processing_solution + + L2_error_mobile = error_L2(mobile_computed, mobile_analytical_np) + L2_error_trapped = error_L2(trapped_computed, trapped_analytical_np) + + assert L2_error_mobile < 3e-02 + assert L2_error_trapped < 9e-03 + + +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""" + + def u_exact(mod): + return lambda x: 1.5 + mod.sin(3 * mod.pi * x[0]) + mod.cos(3 * mod.pi * x[1]) + + def v_exact(mod): + return lambda x: mod.sin(3 * mod.pi * x[0]) + mod.cos(3 * mod.pi * x[1]) + + mobile_analytical_ufl = u_exact(ufl) + mobile_analytical_np = u_exact(np) + trapped_analytical_ufl = v_exact(ufl) + trapped_analytical_np = v_exact(np) + + elements = ufl.FiniteElement("P", test_mesh_2d.ufl_cell(), 1) + V = fem.FunctionSpace(test_mesh_2d, elements) + T = fem.Function(V) + f = fem.Function(V) + g = fem.Function(V) + + k_0 = 2 + E_k = 1.5 + p_0 = 0.2 + E_p = 0.1 + T_expr = lambda x: 500 + 100 * x[0] + T.interpolate(T_expr) + n_trap = 3 + E_D = 0.1 + D_0 = 2 + k_B = F.k_B + D = D_0 * ufl.exp(-E_D / (k_B * T)) + k = k_0 * ufl.exp(-E_k / (k_B * T)) + p = p_0 * ufl.exp(-E_p / (k_B * T)) + + f = ( + -ufl.div(D * ufl.grad(mobile_analytical_ufl(x_2d))) + + k * mobile_analytical_ufl(x_2d) * (n_trap - trapped_analytical_ufl(x_2d)) + - p * trapped_analytical_ufl(x_2d) + ) + + g = p * trapped_analytical_ufl(x_2d) - k * mobile_analytical_ufl(x_2d) * ( + n_trap - trapped_analytical_ufl(x_2d) + ) + + my_model = F.HydrogenTransportProblem() + my_model.mesh = F.Mesh(mesh=test_mesh_2d) + + # create facet meshtags + boundaries = [ + (1, lambda x: np.isclose(x[0], 0)), + (2, lambda x: np.isclose(x[0], 1)), + ] + facet_indices, facet_markers = [], [] + fdim = test_mesh_2d.topology.dim - 1 + for marker, locator in boundaries: + facets = locate_entities(test_mesh_2d, fdim, locator) + facet_indices.append(facets) + facet_markers.append(np.full_like(facets, marker)) + facet_indices = np.hstack(facet_indices).astype(np.int32) + facet_markers = np.hstack(facet_markers).astype(np.int32) + sorted_facets = np.argsort(facet_indices) + my_facet_meshtags = meshtags( + test_mesh_2d, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets] + ) + + # create volume meshtags + vdim = test_mesh_2d.topology.dim + num_cells = test_mesh_2d.topology.index_map(vdim).size_local + mesh_cell_indices = np.arange(num_cells, dtype=np.int32) + tags_volumes = np.full(num_cells, 1, dtype=np.int32) + my_volume_meshtags = meshtags(test_mesh_2d, vdim, mesh_cell_indices, tags_volumes) + + my_model.facet_meshtags = my_facet_meshtags + my_model.volume_meshtags = my_volume_meshtags + + my_mat = F.Material(name="mat", D_0=D_0, E_D=E_D) + vol = F.VolumeSubdomain(id=1, material=my_mat) + left = F.SurfaceSubdomain(id=1) + right = F.SurfaceSubdomain(id=2) + + my_model.subdomains = [vol, left, right] + + mobile = F.Species("mobile") + trapped = F.Species("trapped", mobile=False) + traps = F.ImplicitSpecies(n=n_trap, others=[trapped]) + my_model.species = [mobile, trapped] + + my_model.reactions = [ + F.Reaction( + reactant1=mobile, + reactant2=traps, + product=trapped, + k_0=k_0, + E_k=E_k, + p_0=p_0, + E_p=E_p, + volume=vol, + ) + ] + + my_model.temperature = T_expr + + my_model.boundary_conditions = [ + F.DirichletBC(subdomain=left, value=mobile_analytical_ufl, species=mobile), + F.DirichletBC(subdomain=right, value=mobile_analytical_ufl, species=mobile), + F.DirichletBC(subdomain=left, value=trapped_analytical_ufl, species=trapped), + F.DirichletBC(subdomain=right, value=trapped_analytical_ufl, species=trapped), + ] + + my_model.sources = [ + F.Source(value=f, volume=vol, species=mobile), + F.Source(value=g, volume=vol, species=trapped), + ] + + my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False) + + my_model.initialise() + my_model.run() + + mobile_computed = mobile.post_processing_solution + trapped_computed = trapped.post_processing_solution + + L2_error_mobile = error_L2(mobile_computed, mobile_analytical_np) + L2_error_trapped = error_L2(trapped_computed, trapped_analytical_np) + + assert L2_error_mobile < 4e-03 + assert L2_error_trapped < 2e-03 diff --git a/test/system_tests/test_1_mobile_species.py b/test/system_tests/test_1_mobile_species.py new file mode 100644 index 000000000..f023b5e0f --- /dev/null +++ b/test/system_tests/test_1_mobile_species.py @@ -0,0 +1,302 @@ +import festim as F +import numpy as np +from dolfinx import fem +import ufl +from .tools import error_L2 +from dolfinx.mesh import meshtags, create_unit_square, create_unit_cube, locate_entities +from mpi4py import MPI + +test_mesh_1d = F.Mesh1D(np.linspace(0, 1, 10000)) +test_mesh_2d = create_unit_square(MPI.COMM_WORLD, 50, 50) +test_mesh_3d = create_unit_cube(MPI.COMM_WORLD, 20, 20, 20) +x_1d = ufl.SpatialCoordinate(test_mesh_1d.mesh) +x_2d = ufl.SpatialCoordinate(test_mesh_2d) +x_3d = ufl.SpatialCoordinate(test_mesh_3d) + + +def test_1_mobile_MMS_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]) + + H_analytical_ufl = u_exact(ufl) + H_analytical_np = u_exact(np) + + elements = ufl.FiniteElement("CG", test_mesh_1d.mesh.ufl_cell(), 1) + V = fem.FunctionSpace(test_mesh_1d.mesh, elements) + T = fem.Function(V) + + D_0 = 1 + E_D = 0.1 + T_expr = lambda x: 500 + 100 * x[0] + T.interpolate(T_expr) + D = D_0 * ufl.exp(-E_D / (F.k_B * T)) + + my_model = F.HydrogenTransportProblem() + my_model.mesh = test_mesh_1d + my_mat = F.Material(name="mat", D_0=D_0, E_D=E_D) + 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") + my_model.species = [H] + + my_model.temperature = T_expr + + my_model.boundary_conditions = [ + F.DirichletBC(subdomain=left, value=H_analytical_ufl, species=H), + F.DirichletBC(subdomain=right, value=H_analytical_ufl, species=H), + ] + + f = -ufl.div(D * ufl.grad(H_analytical_ufl(x_1d))) + my_model.sources = [F.Source(value=f, volume=vol, species=H)] + + my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False) + + my_model.initialise() + my_model.run() + + H_computed = H.post_processing_solution + + L2_error = error_L2(H_computed, H_analytical_np) + + assert L2_error < 1e-7 + + +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 + """ + + final_time = 0.1 + + def u_exact(mod): + return lambda x, t: 1 + mod.sin(2 * mod.pi * x[0]) + 2 * t**2 + + def u_exact_alt(mod): + return lambda x: u_exact(mod)(x, final_time) + + H_analytical_ufl = u_exact(ufl) + H_analytical_np = u_exact_alt(np) + + elements = ufl.FiniteElement("P", test_mesh_1d.mesh.ufl_cell(), 1) + V = fem.FunctionSpace(test_mesh_1d.mesh, elements) + T = fem.Function(V) + + D_0 = 1 + E_D = 0.1 + T_expr = lambda x: 600 + 50 * x[0] + T.interpolate(T_expr) + D = D_0 * ufl.exp(-E_D / (F.k_B * T)) + + # FESTIM model + + my_model = F.HydrogenTransportProblem() + my_model.mesh = test_mesh_1d + my_mat = F.Material(name="mat", D_0=D_0, E_D=E_D) + 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") + my_model.species = [H] + + my_model.temperature = T_expr + + my_model.boundary_conditions = [ + F.DirichletBC(subdomain=left, value=H_analytical_ufl, species=H), + F.DirichletBC(subdomain=right, value=H_analytical_ufl, species=H), + ] + + init_value = lambda x: 1 + ufl.sin(2 * ufl.pi * x[0]) + my_model.initial_conditions = [F.InitialCondition(value=init_value, species=H)] + + f = lambda x, t: 4 * t - ufl.div(D * ufl.grad(H_analytical_ufl(x, t))) + my_model.sources = [F.Source(value=f, volume=vol, species=H)] + + my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=final_time) + my_model.settings.stepsize = final_time / 50 + + my_model.initialise() + my_model.run() + + H_computed = H.post_processing_solution + + L2_error = error_L2(H_computed, H_analytical_np) + + assert L2_error < 5e-4 + + +def test_1_mobile_MMS_2D(): + """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]) + + H_analytical_ufl = u_exact(ufl) + H_analytical_np = u_exact(np) + + elements = ufl.FiniteElement("CG", test_mesh_2d.ufl_cell(), 1) + V = fem.FunctionSpace(test_mesh_2d, elements) + T = fem.Function(V) + + D_0 = 1 + E_D = 0.1 + T_expr = lambda x: 500 + 100 * x[0] + T.interpolate(T_expr) + D = D_0 * ufl.exp(-E_D / (F.k_B * T)) + + my_model = F.HydrogenTransportProblem() + my_model.mesh = F.Mesh(mesh=test_mesh_2d) + + # create facet meshtags + boundaries = [ + (1, lambda x: np.isclose(x[0], 0)), + (2, lambda x: np.isclose(x[0], 1)), + ] + facet_indices, facet_markers = [], [] + fdim = test_mesh_2d.topology.dim - 1 + for marker, locator in boundaries: + facets = locate_entities(test_mesh_2d, fdim, locator) + facet_indices.append(facets) + facet_markers.append(np.full_like(facets, marker)) + facet_indices = np.hstack(facet_indices).astype(np.int32) + facet_markers = np.hstack(facet_markers).astype(np.int32) + sorted_facets = np.argsort(facet_indices) + my_facet_meshtags = meshtags( + test_mesh_2d, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets] + ) + + # create volume meshtags + vdim = test_mesh_2d.topology.dim + num_cells = test_mesh_2d.topology.index_map(vdim).size_local + mesh_cell_indices = np.arange(num_cells, dtype=np.int32) + tags_volumes = np.full(num_cells, 1, dtype=np.int32) + my_volume_meshtags = meshtags(test_mesh_2d, vdim, mesh_cell_indices, tags_volumes) + + my_model.facet_meshtags = my_facet_meshtags + my_model.volume_meshtags = my_volume_meshtags + + my_mat = F.Material(name="mat", D_0=D_0, E_D=E_D) + vol = F.VolumeSubdomain(id=1, material=my_mat) + left = F.SurfaceSubdomain(id=1) + right = F.SurfaceSubdomain(id=2) + + my_model.subdomains = [vol, left, right] + + H = F.Species("H") + my_model.species = [H] + + my_model.temperature = T_expr + + my_model.boundary_conditions = [ + F.DirichletBC(subdomain=left, value=H_analytical_ufl, species=H), + F.DirichletBC(subdomain=right, value=H_analytical_ufl, species=H), + ] + + f = -ufl.div(D * ufl.grad(H_analytical_ufl(x_2d))) + my_model.sources = [F.Source(value=f, volume=vol, species=H)] + + my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False) + + my_model.initialise() + my_model.run() + + H_computed = H.post_processing_solution + + L2_error = error_L2(H_computed, H_analytical_np) + + assert L2_error < 2e-3 + + +def test_1_mobile_MMS_3D(): + """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]) + + H_analytical_ufl = u_exact(ufl) + H_analytical_np = u_exact(np) + + elements = ufl.FiniteElement("CG", test_mesh_3d.ufl_cell(), 1) + V = fem.FunctionSpace(test_mesh_3d, elements) + T = fem.Function(V) + + D_0 = 1 + E_D = 0.1 + T_expr = lambda x: 500 + 100 * x[0] + T.interpolate(T_expr) + D = D_0 * ufl.exp(-E_D / (F.k_B * T)) + + my_model = F.HydrogenTransportProblem() + my_model.mesh = F.Mesh(mesh=test_mesh_3d) + + # create facet meshtags + boundaries = [ + (1, lambda x: np.isclose(x[0], 0)), + (2, lambda x: np.isclose(x[0], 1)), + ] + facet_indices, facet_markers = [], [] + fdim = test_mesh_3d.topology.dim - 1 + for marker, locator in boundaries: + facets = locate_entities(test_mesh_3d, fdim, locator) + facet_indices.append(facets) + facet_markers.append(np.full_like(facets, marker)) + facet_indices = np.hstack(facet_indices).astype(np.int32) + facet_markers = np.hstack(facet_markers).astype(np.int32) + sorted_facets = np.argsort(facet_indices) + my_facet_meshtags = meshtags( + test_mesh_3d, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets] + ) + + # create volume meshtags + vdim = test_mesh_3d.topology.dim + num_cells = test_mesh_3d.topology.index_map(vdim).size_local + mesh_cell_indices = np.arange(num_cells, dtype=np.int32) + tags_volumes = np.full(num_cells, 1, dtype=np.int32) + my_volume_meshtags = meshtags(test_mesh_3d, vdim, mesh_cell_indices, tags_volumes) + + my_model.facet_meshtags = my_facet_meshtags + my_model.volume_meshtags = my_volume_meshtags + + my_mat = F.Material(name="mat", D_0=D_0, E_D=E_D) + vol = F.VolumeSubdomain(id=1, material=my_mat) + left = F.SurfaceSubdomain(id=1) + right = F.SurfaceSubdomain(id=2) + + my_model.subdomains = [vol, left, right] + + H = F.Species("H") + my_model.species = [H] + + my_model.temperature = T_expr + + my_model.boundary_conditions = [ + F.DirichletBC(subdomain=left, value=H_analytical_ufl, species=H), + F.DirichletBC(subdomain=right, value=H_analytical_ufl, species=H), + ] + + f = -ufl.div(D * ufl.grad(H_analytical_ufl(x_3d))) + my_model.sources = [F.Source(value=f, volume=vol, species=H)] + + my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False) + + my_model.initialise() + my_model.run() + + H_computed = H.post_processing_solution + + L2_error = error_L2(H_computed, H_analytical_np) + + assert L2_error < 1e-2 diff --git a/test/test_system.py b/test/system_tests/test_reaction_not_in_every_volume.py similarity index 100% rename from test/test_system.py rename to test/system_tests/test_reaction_not_in_every_volume.py diff --git a/test/system_tests/tools.py b/test/system_tests/tools.py new file mode 100644 index 000000000..2641904a6 --- /dev/null +++ b/test/system_tests/tools.py @@ -0,0 +1,34 @@ +import ufl +from dolfinx import fem +import numpy as np +import mpi4py.MPI as MPI + + +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) diff --git a/test/test_source.py b/test/test_source.py index 508113a3b..b3e8f0209 100644 --- a/test/test_source.py +++ b/test/test_source.py @@ -163,17 +163,18 @@ def my_value(t): [1.0], [[1]], [[F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat)]], + [[F.VolumeSubdomain(1, material=dummy_mat)]], None, ], ) def test_TypeError_is_raised_when_volume_wrong_type(volume_input): """Test that a TypeError is raised when the volume is not of type - festim.VolumeSubdomain1D""" + festim.VolumeSubdomain""" my_spe = F.Species("test") with pytest.raises( TypeError, - match="volume must be of type festim.VolumeSubdomain1D", + match="volume must be of type festim.VolumeSubdomain", ): F.Source(volume=volume_input, value=1.0, species=my_spe)