diff --git a/test/test_multispecies_problem.py b/test/test_multispecies_problem.py index 8f39f2520..e2c05ae3e 100644 --- a/test/test_multispecies_problem.py +++ b/test/test_multispecies_problem.py @@ -6,15 +6,20 @@ import os -def test_multispecies_permeation_problem(tmp_path): +def test_multispecies_permeation_problem(): + """Test running a permeation problem with two species permeating from the + same side of a 1D domain with different diffusion coefficients, asserting + that both species match their respecitive analytcial solutions""" + + # festim model L = 3e-04 my_mesh = F.Mesh1D(np.linspace(0, L, num=1001)) my_model = F.HydrogenTransportProblem() my_model.mesh = my_mesh my_mat = F.Material( - D_0={"D": 1.9e-7, "T": 3.8e-7}, - E_D={"D": 0.2, "T": 0.2}, + D_0={"spe_1": 1.9e-7, "spe_2": 3.8e-7}, + E_D={"spe_1": 0.2, "spe_2": 0.2}, name="my_mat", ) my_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, L], material=my_mat) @@ -26,29 +31,21 @@ def test_multispecies_permeation_problem(tmp_path): right_surface, ] - mobile_D = F.Species("D") - mobile_T = F.Species("T") - my_model.species = [mobile_D, mobile_T] + spe_1 = F.Species("spe_1") + spe_2 = F.Species("spe_2") + my_model.species = [spe_1, spe_2] temperature = Constant(my_mesh.mesh, 500.0) my_model.temperature = temperature my_model.boundary_conditions = [ - F.DirichletBC(subdomain=right_surface, value=0, species="D"), - F.DirichletBC(subdomain=right_surface, value=0, species="T"), + 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="D" + subdomain=left_surface, S_0=4.02e21, E_S=1.04, pressure=100, species="spe_1" ), F.SievertsBC( - subdomain=left_surface, S_0=4.02e21, E_S=1.04, pressure=100, species="T" - ), - ] - my_model.exports = [ - F.XDMFExport( - os.path.join(tmp_path, "mobile_concentration_D.xdmf"), field=mobile_D - ), - F.XDMFExport( - os.path.join(tmp_path, "mobile_concentration_T.xdmf"), field=mobile_T + subdomain=left_surface, S_0=4.02e21, E_S=1.04, pressure=100, species="spe_2" ), ] my_model.settings = F.Settings( @@ -59,7 +56,6 @@ def test_multispecies_permeation_problem(tmp_path): ) my_model.settings.stepsize = F.Stepsize(initial_value=1 / 20) - my_model.initialise() my_model.solver.convergence_criterion = "incremental" @@ -73,48 +69,70 @@ def test_multispecies_permeation_problem(tmp_path): times, flux_values = my_model.run() - # ---------------------- analytical solution ----------------------------- - D_1 = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature, mobile_D) - D_2 = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature, mobile_T) + # ---------------------- 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_1 = float(D_1) * S - permeability_2 = float(D_2) * S + permeability_spe_1 = float(D_spe_1) * S + permeability_spe_2 = float(D_spe_2) * S times = np.array(times) - n_array_1 = np.arange(1, 10000)[:, np.newaxis] - n_array_2 = np.arange(1, 10000)[:, np.newaxis] - summation_1 = np.sum( - (-1) ** n_array_1 - * np.exp(-((np.pi * n_array_1) ** 2) * float(D_1) / L**2 * times), + # ##### 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, ) - summation_2 = np.sum( - (-1) ** n_array_2 - * np.exp(-((np.pi * n_array_2) ** 2) * float(D_2) / L**2 * times), + 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 + ) + 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_1 = P_up**0.5 * permeability_1 / L * (2 * summation_1 + 1) - analytical_flux_2 = P_up**0.5 * permeability_2 / L * (2 * summation_2 + 1) - analytical_flux_1 = np.abs(analytical_flux_1) - analytical_flux_2 = np.abs(analytical_flux_2) - flux_values_1 = np.array(np.abs(flux_values[0])) - flux_values_2 = np.array(np.abs(flux_values[1])) - - indices_1 = np.where(analytical_flux_1 > 0.1 * np.max(analytical_flux_1)) - indices_2 = np.where(analytical_flux_2 > 0.1 * np.max(analytical_flux_2)) - analytical_flux_1 = analytical_flux_1[indices_1] - analytical_flux_2 = analytical_flux_2[indices_2] - flux_values_1 = flux_values_1[indices_1] - flux_values_2 = flux_values_2[indices_2] - - relative_error_1 = np.abs((flux_values_1 - analytical_flux_1) / analytical_flux_1) - relative_error_2 = np.abs((flux_values_2 - analytical_flux_2) / analytical_flux_2) - - error_1 = relative_error_1.mean() - error_2 = relative_error_2.mean() - - for err in [error_1, error_2]: + 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 + ) + error_spe_2 = relative_error_spe_2.mean() + + for err in [error_spe_1, error_spe_2]: assert err < 0.01