Skip to content

Fixes #4

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
RemDelaporteMathurin merged 3 commits into from
Oct 4, 2024
Merged

Fixes #4

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
586988b
fixes xdmf file names, reaction volumes
kaelyndunnell Oct 2, 2024
1dd662e
moved to VTX export + fixed trapping rate + black fromatting
RemDelaporteMathurin Oct 2, 2024
18df7b3
added plotting + refined mesh
RemDelaporteMathurin Oct 4, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
115 changes: 81 additions & 34 deletions simple_mb.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,24 +2,36 @@

import festim as F
import numpy as np
import h_transport_materials as htm
import h_transport_materials as htm


my_model = F.HydrogenTransportProblem()

# building 1D mesh, W mb

L = 6e-3 # m
vertices = np.linspace(0,L,num=300)
L = 6e-3 # m
vertices = np.concatenate(
[
np.linspace(0, 30e-9, num=200),
np.linspace(30e-9, 3e-6, num=300),
np.linspace(3e-6, 30e-6, num=200),
np.linspace(30e-6, L, num=200),
]
)
my_model.mesh = F.Mesh1D(vertices)


# W material parameters
w_density = 6.3382e28 # atoms/m3
w_density = 6.3382e28 # atoms/m3
tungsten_diff = htm.diffusivities.filter(material=htm.TUNGSTEN).mean()
tungsten = F.Material(D_0=tungsten_diff.pre_exp.magnitude, E_D=tungsten_diff.act_energy.magnitude, name="tungsten")
tungsten = F.Material(
D_0=tungsten_diff.pre_exp.magnitude,
E_D=tungsten_diff.act_energy.magnitude,
name="tungsten",
)

# subdomains
w_subdomain = F.VolumeSubdomain1D(id=1, borders=[0,L], material=tungsten)
w_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, L], material=tungsten)
inlet = F.SurfaceSubdomain1D(id=1, x=0)
outlet = F.SurfaceSubdomain1D(id=2, x=L)

Expand All @@ -42,20 +54,20 @@


# traps
empty_trap1 = F.ImplicitSpecies( # implicit trap 1
n=6.338e24, # 1e-4 at.fr.
empty_trap1 = F.ImplicitSpecies( # implicit trap 1
n=6.338e24, # 1e-4 at.fr.
others=[trap1_T, trap1_D],
name="empty_trap1",
)

empty_trap2 = F.ImplicitSpecies( # implicit trap 2
empty_trap2 = F.ImplicitSpecies( # implicit trap 2
n=6.338e24,
others=[trap2_T, trap2_D],
name="empty_trap2",
)

empty_trap3 = F.ImplicitSpecies( # fermi-dirac-like trap 3
n=6.338e27, # 1e-1 at.fr.
empty_trap3 = F.ImplicitSpecies( # fermi-dirac-like trap 3
n=6.338e27, # 1e-1 at.fr.
others=[trap3_T, trap3_D],
name="empty_trap3",
)
Expand All @@ -69,71 +81,74 @@
trap2_T,
trap3_D,
trap3_T,
empty_trap1,
empty_trap2,
empty_trap3
]

# hydrogen reactions - 1 per trap per species
my_model.reactions = [
F.Reaction(
k_0=1e13,
k_0=4.1e-7 / (1.1e-10**2 * 6 * w_density),
E_k=0.20,
p_0=1e13,
E_p=0.85,
volume=w_subdomain,
reactant=[mobile_D, empty_trap1],
product=trap1_D,
),
F.Reaction(
k_0=1e13,
k_0=4.1e-7 / (1.1e-10**2 * 6 * w_density),
E_k=0.2,
p_0=1e13,
E_p=0.85,
volume=w_subdomain,
reactant=[mobile_T, empty_trap1],
product=trap1_T,
),
F.Reaction(
k_0=1e13,
k_0=4.1e-7 / (1.1e-10**2 * 6 * w_density),
E_k=0.2,
p_0=1e13,
E_p=1,
volume=w_subdomain,
reactant=[mobile_D, empty_trap2],
product=trap2_D,
),
F.Reaction(
k_0=1e13,
k_0=4.1e-7 / (1.1e-10**2 * 6 * w_density),
E_k=0.2,
p_0=1e13,
E_p=1,
volume=w_subdomain,
reactant=[mobile_T, empty_trap2],
product=trap2_T,
),
F.Reaction(
k_0=1e13,
k_0=4.1e-7 / (1.1e-10**2 * 6 * w_density),
E_k=0.2,
p_0=1e13,
E_p=1.5,
volume=w_subdomain,
reactant=[mobile_D, empty_trap3],
product=trap3_D,
),
F.Reaction(
k_0=1e13,
k_0=4.1e-7 / (1.1e-10**2 * 6 * w_density),
E_k=0.2,
p_0=1e13,
E_p=1.5,
volume=w_subdomain,
reactant=[mobile_T, empty_trap3],
product=trap3_T,
)
),
]

# temperature
# temperature
# my_model.temperature = 1000 - (1000-350)/L*x
my_model.temperature = 10000
my_model.temperature = 1000

# boundary conditions

my_model.boundary_conditions = [
F.DirichletBC(subdomain=inlet, value=1e20, species=mobile_T),
F.DirichletBC(subdomain=inlet, value=1e20, species=mobile_T),
F.DirichletBC(subdomain=inlet, value=1e19, species=mobile_D),
F.DirichletBC(subdomain=outlet, value=0, species=mobile_T),
F.DirichletBC(subdomain=outlet, value=0, species=mobile_D),
Expand All @@ -147,16 +162,22 @@
folder = "multi_isotope_trapping_example"

my_model.exports = [
F.XDMFExport(f"{folder}/mobile_concentration_t.xdmf", field=mobile_T),
F.XDMFExport(f"{folder}/mobile_concentration_d.xdmf", field=mobile_D),
F.XDMFExport(f"{folder}/trapped_concentration_t1.xdmf", field=trap1_D),
F.XDMFExport(f"{folder}/trapped_concentration_t2.xdmf", field=trap1_T),
F.XDMFExport(f"{folder}/trapped_concentration_d1.xdmf", field=trap2_D),
F.XDMFExport(f"{folder}/trapped_concentration_d2.xdmf", field=trap2_T),
F.XDMFExport(f"{folder}/trapped_concentration_dt.xdmf", field=trap3_D),
F.XDMFExport(f"{folder}/trapped_concentration_dt.xdmf", field=trap3_T)
F.VTXExport(f"{folder}/mobile_concentration_t.bp", field=mobile_T),
F.VTXExport(f"{folder}/mobile_concentration_d.bp", field=mobile_D),
F.VTXExport(f"{folder}/trapped_concentration_d1.bp", field=trap1_D),
F.VTXExport(f"{folder}/trapped_concentration_t1.bp", field=trap1_T),
F.VTXExport(f"{folder}/trapped_concentration_d2.bp", field=trap2_D),
F.VTXExport(f"{folder}/trapped_concentration_t2.bp", field=trap2_T),
F.VTXExport(f"{folder}/trapped_concentration_d3.bp", field=trap3_D),
F.VTXExport(f"{folder}/trapped_concentration_t3.bp", field=trap3_T),
]

quantities = {}
for species in my_model.species:
quantity = F.TotalVolume(field=species, volume=w_subdomain)
my_model.exports.append(quantity)
quantities[species.name] = quantity

# settings

my_model.settings = F.Settings(
Expand All @@ -168,6 +189,32 @@
# run simu

my_model.initialise()
my_model.run()


import matplotlib.pyplot as plt

for name, quantity in quantities.items():
plt.plot(quantity.t, quantity.data, label=name)

plt.xlabel("Time (s)")
plt.ylabel("Total quantity (atoms/m2)")
plt.legend()
plt.yscale("log")

plt.show()

# make the same but with a stack plot

fig, ax = plt.subplots()

ax.stackplot(
quantity.t,
[quantity.data for quantity in quantities.values()],
labels=quantities.keys(),
)

print(my_model.formulation)
my_model.run()
plt.xlabel("Time (s)")
plt.ylabel("Total quantity (atoms/m2)")
plt.legend()
plt.show()