Skip to content

Commit

Permalink
updated tests for multispecies
Browse files Browse the repository at this point in the history
  • Loading branch information
@jhdark
jhdark committed Oct 20, 2023
1 parent f64076d commit 9bf4373
Show file tree
Hide file tree
Showing 3 changed files with 107 additions and 104 deletions.
21 changes: 11 additions & 10 deletions test/test_dirichlet_bc.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,13 +57,12 @@ def test_callable_for_value():
subdomain = F.SurfaceSubdomain1D(1, x=1)
vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat)
value = lambda x, t: 1.0 + x[0] + t
species = "test"
species = F.Species("test")

bc = F.DirichletBC(subdomain, value, species)

my_model = F.HydrogenTransportProblem(
mesh=F.Mesh(mesh),
subdomains=[subdomain, vol_subdomain],
mesh=F.Mesh(mesh), subdomains=[subdomain, vol_subdomain], species=[species]
)

my_model.define_function_space()
Expand Down Expand Up @@ -96,13 +95,12 @@ def test_value_callable_x_t_T():
subdomain = F.SurfaceSubdomain1D(1, x=1)
vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat)
value = lambda x, t, T: 1.0 + x[0] + t + T
species = "test"
species = F.Species("test")

bc = F.DirichletBC(subdomain, value, species)

my_model = F.HydrogenTransportProblem(
mesh=F.Mesh(mesh),
subdomains=[subdomain, vol_subdomain],
mesh=F.Mesh(mesh), subdomains=[subdomain, vol_subdomain], species=[species]
)

my_model.define_function_space()
Expand Down Expand Up @@ -138,13 +136,14 @@ def test_callable_t_only():
subdomain = F.SurfaceSubdomain1D(1, x=1)
vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat)
value = lambda t: 1.0 + t
species = "test"
species = F.Species("test")

bc = F.DirichletBC(subdomain, value, species)

my_model = F.HydrogenTransportProblem(
mesh=F.Mesh(mesh),
subdomains=[subdomain, vol_subdomain],
species=[species],
)

my_model.define_function_space()
Expand Down Expand Up @@ -180,13 +179,14 @@ def test_callable_x_only():
subdomain = F.SurfaceSubdomain1D(1, x=1)
vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat)
value = lambda x: 1.0 + x[0]
species = "test"
species = F.Species("test")

bc = F.DirichletBC(subdomain, value, species)

my_model = F.HydrogenTransportProblem(
mesh=F.Mesh(mesh),
subdomains=[subdomain, vol_subdomain],
species=[species],
)

my_model.define_function_space()
Expand Down Expand Up @@ -230,13 +230,14 @@ def test_create_formulation(value):
# BUILD
subdomain = F.SurfaceSubdomain1D(1, x=1)
vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat)
species = "test"
species = F.Species("test")

bc = F.DirichletBC(subdomain, value, species)

my_model = F.HydrogenTransportProblem(
mesh=F.Mesh(mesh),
subdomains=[subdomain, vol_subdomain],
species=[species],
)

my_model.define_function_space()
Expand Down Expand Up @@ -276,7 +277,7 @@ def test_integration_with_HTransportProblem(value):
subdomains=[vol_subdomain, subdomain],
)
my_model.species = [F.Species("H")]
my_bc = F.DirichletBC(subdomain, value, my_model.species[0])
my_bc = F.DirichletBC(subdomain, value, "H")
my_model.boundary_conditions = [my_bc]

my_model.temperature = fem.Constant(my_model.mesh.mesh, 550.0)
Expand Down
2 changes: 1 addition & 1 deletion test/test_material.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ def test_define_diffusion_coefficient():
T, D_0, E_D = 10, 1.2, 0.5

my_mat = F.Material(D_0=D_0, E_D=E_D)
D = my_mat.get_diffusion_coefficient(test_mesh.mesh, T)
D = my_mat.get_diffusion_coefficient(test_mesh.mesh, T, species="dummy")

D_analytical = D_0 * np.exp(-E_D / F.k_B / T)

Expand Down
188 changes: 95 additions & 93 deletions test/test_permeation_problem.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,9 @@ def test_permeation_problem(mesh_size=1001):

# -------------------------- analytical solution -------------------------------------

D = my_mat.get_diffusion_coefficient(my_mesh.mesh, my_model.temperature)
D = my_mat.get_diffusion_coefficient(
my_mesh.mesh, my_model.temperature, my_model.species[0]
)

S_0 = float(my_model.boundary_conditions[-1].S_0)
E_S = float(my_model.boundary_conditions[-1].E_S)
Expand Down Expand Up @@ -89,95 +91,95 @@ def test_permeation_problem(mesh_size=1001):
assert error < 0.01


def test_permeation_problem_multi_volume():
L = 3e-04
vertices = np.linspace(0, L, num=1001)

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_1 = F.VolumeSubdomain1D(id=1, borders=[0, L / 4], material=my_mat)
my_subdomain_2 = F.VolumeSubdomain1D(id=2, borders=[L / 4, L / 2], material=my_mat)
my_subdomain_3 = F.VolumeSubdomain1D(
id=3, borders=[L / 2, 3 * L / 4], material=my_mat
)
my_subdomain_4 = F.VolumeSubdomain1D(id=4, borders=[3 * L / 4, 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_1,
my_subdomain_2,
my_subdomain_3,
my_subdomain_4,
left_surface,
right_surface,
]

mobile_H = F.Species("H")
my_model.species = [mobile_H]

temperature = Constant(my_mesh.mesh, 500.0)
my_model.temperature = temperature

my_model.boundary_conditions = [
F.DirichletBC(subdomain=right_surface, value=0, species="H"),
F.SievertsBC(
subdomain=left_surface, S_0=4.02e21, E_S=1.04, pressure=100, species="H"
),
]
my_model.exports = [F.VTXExport("test.bp", field=mobile_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()

# -------------------------- analytical solution -------------------------------------
D = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature)

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 = 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]

relative_error = np.abs((flux_values - analytical_flux) / analytical_flux)

error = relative_error.mean()

assert error < 0.01
# def test_permeation_problem_multi_volume():
# L = 3e-04
# vertices = np.linspace(0, L, num=1001)

# 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_1 = F.VolumeSubdomain1D(id=1, borders=[0, L / 4], material=my_mat)
# my_subdomain_2 = F.VolumeSubdomain1D(id=2, borders=[L / 4, L / 2], material=my_mat)
# my_subdomain_3 = F.VolumeSubdomain1D(
# id=3, borders=[L / 2, 3 * L / 4], material=my_mat
# )
# my_subdomain_4 = F.VolumeSubdomain1D(id=4, borders=[3 * L / 4, 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_1,
# my_subdomain_2,
# my_subdomain_3,
# my_subdomain_4,
# left_surface,
# right_surface,
# ]

# mobile_H = F.Species("H")
# my_model.species = [mobile_H]

# temperature = Constant(my_mesh.mesh, 500.0)
# my_model.temperature = temperature

# my_model.boundary_conditions = [
# F.DirichletBC(subdomain=right_surface, value=0, species="H"),
# F.SievertsBC(
# subdomain=left_surface, S_0=4.02e21, E_S=1.04, pressure=100, species="H"
# ),
# ]
# my_model.exports = [F.VTXExport("test.bp", field=mobile_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()

# # -------------------------- analytical solution -------------------------------------
# D = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature)

# 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 = 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]

# relative_error = np.abs((flux_values - analytical_flux) / analytical_flux)

# error = relative_error.mean()

# assert error < 0.01

0 comments on commit 9bf4373

Please sign in to comment.