Solution#
Here you will find the solution for the challenge D.
import meshio
import fenics as f
import h_transport_materials as htm
import matplotlib.pyplot as plt
import numpy as np
import festim as F
def convert_med_to_xdmf(
med_file,
cell_file="mesh_domains.xdmf",
facet_file="mesh_boundaries.xdmf",
cell_type="tetra",
facet_type="triangle",
):
"""Converts a MED mesh to XDMF
Args:
med_file (str): the name of the MED file
cell_file (str, optional): the name of the file containing the
volume markers. Defaults to "mesh_domains.xdmf".
facet_file (str, optional): the name of the file containing the
surface markers.. Defaults to "mesh_boundaries.xdmf".
cell_type (str, optional): The topology of the cells. Defaults to "tetra".
facet_type (str, optional): The topology of the facets. Defaults to "triangle".
Returns:
dict, dict: the correspondance dict, the cell types
"""
msh = meshio.read(med_file)
correspondance_dict = msh.cell_tags
cell_data_types = msh.cell_data_dict["cell_tags"].keys()
for mesh_block in msh.cells:
if mesh_block.type == cell_type:
meshio.write_points_cells(
cell_file,
msh.points,
[mesh_block],
cell_data={"f": [-1 * msh.cell_data_dict["cell_tags"][cell_type]]},
)
elif mesh_block.type == facet_type:
meshio.write_points_cells(
facet_file,
msh.points,
[mesh_block],
cell_data={"f": [-1 * msh.cell_data_dict["cell_tags"][facet_type]]},
)
return correspondance_dict, cell_data_types
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[1], line 3
1 import meshio
2 import fenics as f
----> 3 import h_transport_materials as htm
4 import matplotlib.pyplot as plt
5 import numpy as np
6 import festim as F
File ~/checkouts/readthedocs.org/user_builds/festim-workshop/conda/festim1/lib/python3.11/site-packages/h_transport_materials/__init__.py:25
22 Rg = 8.314 * ureg.Pa * ureg.m**3 * ureg.mol**-1 * ureg.K**-1
23 avogadro_nb = 6.022e23 * ureg.particle * ureg.mol**-1
---> 25 from pybtex.database import parse_file
26 from pathlib import Path
28 bib_database = parse_file(str(Path(__file__).parent) + "/references.bib")
File ~/checkouts/readthedocs.org/user_builds/festim-workshop/conda/festim1/lib/python3.11/site-packages/pybtex/database/__init__.py:44
42 from pybtex.errors import report_error
43 from pybtex.py3compat import fix_unicode_literals_in_doctest, python_2_unicode_compatible
---> 44 from pybtex.plugin import find_plugin
47 # for python2 compatibility
48 def indent(text, prefix):
File ~/checkouts/readthedocs.org/user_builds/festim-workshop/conda/festim1/lib/python3.11/site-packages/pybtex/plugin/__init__.py:26
2 # Copyright (c) 2006-2021 Andrey Golovizin
3 # Copyright (c) 2014 Matthias C. M. Troffaes
4 #
(...) 21 # TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
22 # SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
25 import os.path # splitext
---> 26 import pkg_resources
28 from pybtex.exceptions import PybtexError
31 class Plugin(object):
ModuleNotFoundError: No module named 'pkg_resources'
mesh_file = "../challenge_D/challenge_mesh.med"
correspondance_dict, cell_data_types = convert_med_to_xdmf(
mesh_file,
cell_type="triangle",
facet_type="line",
cell_file="../challenge_D/mesh_domains.xdmf",
facet_file="../challenge_D/mesh_boundaries.xdmf",
)
print("This is the correspondance dict:")
print(correspondance_dict)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[2], line 3
1 mesh_file = "../challenge_D/challenge_mesh.med"
2
----> 3 correspondance_dict, cell_data_types = convert_med_to_xdmf(
4 mesh_file,
5 cell_type="triangle",
6 facet_type="line",
NameError: name 'convert_med_to_xdmf' is not defined
# IDs for volumes and surfaces (must be the same as in xdmf files)
id_solid = 6
id_inner_boundary = 7
id_outer_boundary = 8
my_model = F.Simulation(log_level=40)
# define mesh
my_model.mesh = F.MeshFromXDMF(
volume_file="../challenge_D/mesh_domains.xdmf",
boundary_file="../challenge_D/mesh_boundaries.xdmf",
)
# material properties
D_mat = htm.diffusivities.filter(material="gold")[0]
D_0 = D_mat.pre_exp.magnitude
E_D = D_mat.pre_exp.magnitude
my_mat = F.Material(id=id_solid, D_0=D_0, E_D=E_D, thermal_cond=315)
my_model.materials = F.Materials([my_mat])
# define temperature
my_model.T = F.HeatTransferProblem(transient=False)
# define boundary conditions
my_model.boundary_conditions = [
F.ImplantationDirichlet(
surfaces=id_inner_boundary,
phi=1e19,
R_p=3e-09,
D_0=D_0,
E_D=E_D,
),
F.DirichletBC(surfaces=id_outer_boundary, value=0, field="solute"),
F.DirichletBC(surfaces=id_inner_boundary, value=600, field="T"),
F.ConvectiveFlux(h_coeff=1e4, T_ext=300, surfaces=id_outer_boundary),
]
# define exports
my_derived_quantities = F.DerivedQuantities(
[
F.SurfaceFlux(field="solute", surface=id_outer_boundary),
F.TotalVolume(field="solute", volume=id_solid),
],
show_units=True,
)
my_model.exports = F.Exports([my_derived_quantities])
my_model.dt = F.Stepsize(
initial_value=1,
stepsize_change_ratio=1.05,
dt_min=1e-04,
)
my_model.settings = F.Settings(
transient=True,
maximum_iterations=30,
final_time=1.2e03,
absolute_tolerance=1e-10,
relative_tolerance=1e-10,
)
my_model.initialise()
my_model.run()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[3], line 8
4
5 id_inner_boundary = 7
6 id_outer_boundary = 8
7
----> 8 my_model = F.Simulation(log_level=40)
9
10 # define mesh
11 my_model.mesh = F.MeshFromXDMF(
NameError: name 'F' is not defined
surface_flux = np.abs(my_derived_quantities[0].data)
inv = my_derived_quantities[1].data
t = my_derived_quantities.t
plt.figure()
plt.title("Outer surface flux")
plt.plot(t, surface_flux, color="black")
plt.xlabel(r"Time (s)")
plt.ylabel(r"Surface flux (H/m2/s)")
ax = plt.gca()
ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)
plt.tight_layout
print(f"Steady surface flux = {surface_flux[-1]:.2e} H/m2/s")
plt.figure()
plt.title("Inventory")
plt.plot(t, inv, color="black")
plt.xlabel(r"Time (s)")
plt.ylabel(r"Inventory (H/m)")
ax = plt.gca()
ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)
plt.tight_layout
print(f"Steady inventory = {inv[-1]:.2e} H/m")
plt.figure()
sol = my_model.T.T
CS = f.plot(sol, cmap="inferno")
plt.colorbar(CS)
plt.title(r"Temperature (K)")
plt.figure()
sol = my_model.h_transport_problem.mobile.mobile_concentration()
CS = f.plot(sol)
plt.colorbar(CS)
plt.title(r"Mobile tritium concentration (H/m3)")
plt.show()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[4], line 1
----> 1 surface_flux = np.abs(my_derived_quantities[0].data)
2 inv = my_derived_quantities[1].data
3 t = my_derived_quantities.t
4
NameError: name 'np' is not defined