Skip to content

Commit

Permalink
refactoring for clarity
Browse files Browse the repository at this point in the history
  • Loading branch information
@jhdark
jhdark committed Oct 25, 2023
1 parent 6de35aa commit d254c87
Showing 1 changed file with 71 additions and 53 deletions.
124 changes: 71 additions & 53 deletions test/test_multispecies_problem.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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(
Expand All @@ -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"
Expand All @@ -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

0 comments on commit d254c87

Please sign in to comment.