From d5c4feffef2b987dd3d78dec110d5a7783022ea0 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 25 Oct 2023 16:17:19 -0400 Subject: [PATCH] updates based on comments --- festim/hydrogen_transport_problem.py | 3 + festim/species.py | 9 +- test/test_multispecies_problem.py | 120 +++++++++++++-------------- test/test_permeation_problem.py | 68 +++++++-------- 4 files changed, 98 insertions(+), 102 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index e2aff9a16..ff171bca3 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -47,6 +47,7 @@ class HydrogenTransportProblem: model formulation (ufl.form.Form): the formulation of the model solver (dolfinx.nls.newton.NewtonSolver): the solver of the model + multispecies (bool): True if the model has more than one species Usage: >>> import festim as F @@ -380,6 +381,8 @@ def run(self): self.solver.solve(self.u) + # post processing + # TODO remove this if not self.multispecies: D_D = self.subdomains[0].material.get_diffusion_coefficient( self.mesh.mesh, self.temperature, self.species[0] diff --git a/festim/species.py b/festim/species.py index be5e528a6..663436172 100644 --- a/festim/species.py +++ b/festim/species.py @@ -15,10 +15,11 @@ class Species: timestep test_function (ufl.Argument): the testfunction associated with this species - sub_function_space (dolfinx.fem.FunctionSpace): the subspace of the - function space - collapsed_function_space (dolfinx.fem.FunctionSpace): the collapsed - function space for a species in the function space + sub_function_space (dolfinx.fem.function.FunctionSpaceBase): the + subspace of the function space + collapsed_function_space (dolfinx.fem.function.FunctionSpaceBase): the + collapsed function space for a species in the function space. In + case single species case, this is None. post_processing_solution (dolfinx.fem.Function): the solution for post processing concentration (dolfinx.fem.Function): the concentration of the species diff --git a/test/test_multispecies_problem.py b/test/test_multispecies_problem.py index f7153a881..ebae0dcf5 100644 --- a/test/test_multispecies_problem.py +++ b/test/test_multispecies_problem.py @@ -3,14 +3,35 @@ from petsc4py import PETSc from dolfinx.fem import Constant from ufl import exp -import os + + +def relative_error_computed_to_analytical( + D, permeability, computed_flux, L, times, P_up +): + n_array = np.arange(1, 10000)[:, np.newaxis] + summation = np.sum( + (-1) ** n_array * np.exp(-((np.pi * n_array) ** 2) * float(D) / L**2 * times), + axis=0, + ) + analytical_flux = P_up**0.5 * permeability / L * (2 * summation + 1) + + # post processing + analytical_flux = np.abs(analytical_flux) + indices = np.where(analytical_flux > 0.1 * np.max(analytical_flux)) + analytical_flux = analytical_flux[indices] + computed_flux = computed_flux[indices] + + # evaulate relative error compared to analytical solution + relative_error = np.abs((computed_flux - analytical_flux) / analytical_flux) + error = relative_error.mean() + + return error def test_multispecies_permeation_problem(): """Test running a problem with 2 mobile species permeating through a 1D - 0.3mm domain, with different diffusion coefficients, asserting that the - resulting concentration fields are less than 1% different from their - respecitive analytical solutions""" + domain, with different diffusion coefficients, checks that the computed + permeation flux matches the analytical solution""" # festim model L = 3e-04 @@ -38,15 +59,24 @@ def test_multispecies_permeation_problem(): temperature = Constant(my_mesh.mesh, 500.0) my_model.temperature = temperature + pressure = 100 my_model.boundary_conditions = [ F.DirichletBC(subdomain=right_surface, value=0, species="spe_1"), F.DirichletBC(subdomain=right_surface, value=0, species="spe_2"), F.SievertsBC( - subdomain=left_surface, S_0=4.02e21, E_S=1.04, pressure=100, species="spe_1" + subdomain=left_surface, + S_0=4.02e21, + E_S=1.04, + pressure=pressure, + species="spe_1", ), F.SievertsBC( - subdomain=left_surface, S_0=4.02e21, E_S=1.04, pressure=100, species="spe_2" + subdomain=left_surface, + S_0=5.0e21, + E_S=1.2, + pressure=pressure, + species="spe_2", ), ] my_model.settings = F.Settings( @@ -71,69 +101,35 @@ def test_multispecies_permeation_problem(): times, flux_values = my_model.run() # ---------------------- analytical solutions ----------------------------- - D_spe_1 = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature, spe_1) - D_spe_2 = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature, spe_2) - S_0 = float(my_model.boundary_conditions[-1].S_0) - E_S = float(my_model.boundary_conditions[-1].E_S) - P_up = float(my_model.boundary_conditions[-1].pressure) - S = S_0 * exp(-E_S / F.k_B / float(temperature)) - permeability_spe_1 = float(D_spe_1) * S - permeability_spe_2 = float(D_spe_2) * S + + # common values times = np.array(times) + P_up = float(my_model.boundary_conditions[-1].pressure) + flux_values_spe_1, flux_values_spe_2 = flux_values[0], flux_values[1] # ##### compute analyical solution for species 1 ##### # - n_array = np.arange(1, 10000)[:, np.newaxis] - - summation_spe_1 = np.sum( - (-1) ** n_array - * np.exp(-((np.pi * n_array) ** 2) * float(D_spe_1) / L**2 * times), - axis=0, - ) - analytical_flux_spe_1 = ( - P_up**0.5 * permeability_spe_1 / L * (2 * summation_spe_1 + 1) - ) - - # post processing - analytical_flux_spe_1 = np.abs(analytical_flux_spe_1) - flux_values_spe_1 = np.array(np.abs(flux_values[0])) - indices_spe_1 = np.where( - analytical_flux_spe_1 > 0.1 * np.max(analytical_flux_spe_1) - ) - - analytical_flux_spe_1 = analytical_flux_spe_1[indices_spe_1] - flux_values_spe_1 = flux_values_spe_1[indices_spe_1] - - # evaluate relative error compared to analytical solution - relative_error_spe_1 = np.abs( - (flux_values_spe_1 - analytical_flux_spe_1) / analytical_flux_spe_1 + D_spe_1 = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature, spe_1) + S_0_spe_1 = float(my_model.boundary_conditions[-2].S_0) + E_S_spe_1 = float(my_model.boundary_conditions[-2].E_S) + S_spe_1 = S_0_spe_1 * exp(-E_S_spe_1 / F.k_B / float(temperature)) + permeability_spe_1 = float(D_spe_1) * S_spe_1 + flux_values_spe_1 = np.array(np.abs(flux_values_spe_1)) + + error_spe_1 = relative_error_computed_to_analytical( + D_spe_1, permeability_spe_1, flux_values_spe_1, L, times, P_up ) - error_spe_1 = relative_error_spe_1.mean() # ##### compute analyical solution for species 2 ##### # - summation_spe_2 = np.sum( - (-1) ** n_array - * np.exp(-((np.pi * n_array) ** 2) * float(D_spe_2) / L**2 * times), - axis=0, - ) - analytical_flux_spe_2 = ( - P_up**0.5 * permeability_spe_2 / L * (2 * summation_spe_2 + 1) - ) - - # post processing - analytical_flux_spe_2 = np.abs(analytical_flux_spe_2) - flux_values_spe_2 = np.array(np.abs(flux_values[1])) - indices_spe_2 = np.where( - analytical_flux_spe_2 > 0.1 * np.max(analytical_flux_spe_2) - ) - - analytical_flux_spe_2 = analytical_flux_spe_2[indices_spe_2] - flux_values_spe_2 = flux_values_spe_2[indices_spe_2] - - # evaluate relative error compared to analytical solution - relative_error_spe_2 = np.abs( - (flux_values_spe_2 - analytical_flux_spe_2) / analytical_flux_spe_2 + D_spe_2 = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature, spe_2) + S_0_spe_2 = float(my_model.boundary_conditions[-1].S_0) + E_S_spe_2 = float(my_model.boundary_conditions[-1].E_S) + S_spe_2 = S_0_spe_2 * exp(-E_S_spe_2 / F.k_B / float(temperature)) + permeability_spe_2 = float(D_spe_2) * S_spe_2 + flux_values_spe_2 = np.array(np.abs(flux_values_spe_2)) + + error_spe_2 = relative_error_computed_to_analytical( + D_spe_2, permeability_spe_2, flux_values_spe_2, L, times, P_up ) - error_spe_2 = relative_error_spe_2.mean() for err in [error_spe_1, error_spe_2]: assert err < 0.01 diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index 3ec522846..8e2eb014a 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -6,10 +6,33 @@ import os +def relative_error_computed_to_analytical( + D, permeability, computed_flux, L, times, P_up +): + n_array = np.arange(1, 10000)[:, np.newaxis] + summation = np.sum( + (-1) ** n_array * np.exp(-((np.pi * n_array) ** 2) * float(D) / L**2 * times), + axis=0, + ) + analytical_flux = P_up**0.5 * permeability / L * (2 * summation + 1) + + # post processing + analytical_flux = np.abs(analytical_flux) + indices = np.where(analytical_flux > 0.1 * np.max(analytical_flux)) + analytical_flux = analytical_flux[indices] + computed_flux = computed_flux[indices] + + # evaulate relative error compared to analytical solution + relative_error = np.abs((computed_flux - analytical_flux) / analytical_flux) + error = relative_error.mean() + + return error + + def test_permeation_problem(mesh_size=1001): - """Test running a problem with a mobile species permeating through a 1D 0.3mm domain - asserting that the resulting concentration field is less than 1% different from a - respecitive 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 @@ -71,25 +94,11 @@ def test_permeation_problem(mesh_size=1001): S = S_0 * exp(-E_S / F.k_B / float(my_model.temperature)) permeability = float(D) * S times = np.array(times) - - n_array = np.arange(1, 10000)[:, np.newaxis] - summation = np.sum( - (-1) ** n_array * np.exp(-((np.pi * n_array) ** 2) * float(D) / L**2 * times), - axis=0, - ) - analytical_flux = P_up**0.5 * permeability / L * (2 * summation + 1) - - # post processing - analytical_flux = np.abs(analytical_flux) flux_values = np.array(np.abs(flux_values)) - indices = np.where(analytical_flux > 0.01 * np.max(analytical_flux)) - analytical_flux = analytical_flux[indices] - flux_values = flux_values[indices] - - # evaulate relative error compared to analytical solution - relative_error = np.abs((flux_values - analytical_flux) / analytical_flux) - error = relative_error.mean() + error = relative_error_computed_to_analytical( + D, permeability, flux_values, L, times, P_up + ) assert error < 0.01 @@ -171,23 +180,10 @@ def test_permeation_problem_multi_volume(tmp_path): S = S_0 * exp(-E_S / F.k_B / float(temperature)) permeability = float(D) * S times = np.array(times) - - n_array = np.arange(1, 10000)[:, np.newaxis] - summation = np.sum( - (-1) ** n_array * np.exp(-((np.pi * n_array) ** 2) * float(D) / L**2 * times), - axis=0, - ) - analytical_flux = P_up**0.5 * permeability / L * (2 * summation + 1) - - analytical_flux = np.abs(analytical_flux) flux_values = np.array(np.abs(flux_values)) - indices = np.where(analytical_flux > 0.01 * np.max(analytical_flux)) - analytical_flux = analytical_flux[indices] - flux_values = flux_values[indices] - - # evaulate relative error compared to analytical solution - relative_error = np.abs((flux_values - analytical_flux) / analytical_flux) - error = relative_error.mean() + error = relative_error_computed_to_analytical( + D, permeability, flux_values, L, times, P_up + ) assert error < 0.01