Non-cartesian meshes#
Objectives
Tips and tricks to better handle high concentration gradients or discontinuities
FESTIM can run simulations in cylindrical or spherical coordinates.
This is extremely useful to simulate axisymmetric systems. Here we show two examples with a pipe and a sphere.
For instance, on the geometry below, instead of simulating the full 3D domain, one can simply simulate a 2D cross section and assume axial symmetry. If the hollow cylinder is infinitely high, then it is even possible to only simulate a 1D domain.
The cylindrical system of coordinates can be applied to 1D and 2D domains. The spherical system of coordinates can only be applied to 1D domains.
Pipe#
The first example is a 1D axisymmetric simulation representing a pipe.
import festim as F
import numpy as np
my_model = F.Simulation()
Cylindrical coordinates can be set in the type argument of Mesh:
inner_radius = 0.1e-2
thickness = 0.2e-2
my_model.mesh = F.MeshFromVertices(
np.linspace(inner_radius, inner_radius + thickness, num=100),
type="cylindrical",
)
We’ll use HTM data for the properties of steel:
import h_transport_materials as htm
steel_D = htm.diffusivities.filter(material=htm.Steel).mean()
steel_S = htm.solubilities.filter(material=htm.Steel).mean()
steel = F.Material(
id=1,
D_0=steel_D.pre_exp.magnitude,
E_D=steel_D.act_energy.magnitude,
S_0=steel_S.pre_exp.magnitude,
E_S=steel_S.act_energy.magnitude,
)
my_model.materials = steel
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[3], line 1
----> 1 import h_transport_materials as htm
2
3 steel_D = htm.diffusivities.filter(material=htm.Steel).mean()
4 steel_S = htm.solubilities.filter(material=htm.Steel).mean()
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'
We set a Sievert’s boundary condition at the inner surface and a \(c=0\) boundary condition at the outer surface.
my_model.boundary_conditions = [
F.SievertsBC(surfaces=[1], S_0=steel.S_0, E_S=steel.E_S, pressure=1e5),
F.DirichletBC(surfaces=[2], value=0, field="solute"),
]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[4], line 2
1 my_model.boundary_conditions = [
----> 2 F.SievertsBC(surfaces=[1], S_0=steel.S_0, E_S=steel.E_S, pressure=1e5),
3 F.DirichletBC(surfaces=[2], value=0, field="solute"),
4 ]
NameError: name 'steel' is not defined
Note: at the moment, AverageVolume and TotalVolume are not implemented for non-cartesian coordinates (see Issue #723)
my_model.T = 500
my_model.settings = F.Settings(
absolute_tolerance=1e10,
relative_tolerance=1e-10,
final_time=10000,
)
my_model.dt = F.Stepsize(10, stepsize_change_ratio=1.1)
derived_quantities = F.DerivedQuantities(
[
F.SurfaceFluxCylindrical(field="solute", surface=1),
F.SurfaceFluxCylindrical(field="solute", surface=2),
],
show_units=True,
)
txt_export = F.TXTExport(field="solute", filename="task14/solute.txt")
my_model.exports = [derived_quantities, txt_export]
my_model.initialise()
my_model.run()
import matplotlib.pyplot as plt
cmap = plt.get_cmap("viridis")
data = np.genfromtxt("task14/solute.txt", skip_header=1, delimiter=",")
for i in range(1, data.shape[1], 3):
plt.plot(data[:, 0] / 1e-3, data[:, i], color=cmap(i/data.shape[1]))
plt.xlabel("x [mm]")
plt.ylabel("c [H/m³]")
plt.show()
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
Cell In[6], line 5
1 import matplotlib.pyplot as plt
2
3 cmap = plt.get_cmap("viridis")
4
----> 5 data = np.genfromtxt("task14/solute.txt", skip_header=1, delimiter=",")
6 for i in range(1, data.shape[1], 3):
7 plt.plot(data[:, 0] / 1e-3, data[:, i], color=cmap(i/data.shape[1]))
8
File ~/checkouts/readthedocs.org/user_builds/festim-workshop/conda/festim1/lib/python3.11/site-packages/numpy/lib/_npyio_impl.py:1978, in genfromtxt(fname, dtype, comments, delimiter, skip_header, skip_footer, converters, missing_values, filling_values, usecols, names, excludelist, deletechars, replace_space, autostrip, case_sensitive, defaultfmt, unpack, usemask, loose, invalid_raise, max_rows, encoding, ndmin, like)
1976 fname = os.fspath(fname)
1977 if isinstance(fname, str):
-> 1978 fid = np.lib._datasource.open(fname, 'rt', encoding=encoding)
1979 fid_ctx = contextlib.closing(fid)
1980 else:
File ~/checkouts/readthedocs.org/user_builds/festim-workshop/conda/festim1/lib/python3.11/site-packages/numpy/lib/_datasource.py:192, in open(path, mode, destpath, encoding, newline)
155 """
156 Open `path` with `mode` and return the file object.
157
(...) 188
189 """
191 ds = DataSource(destpath)
--> 192 return ds.open(path, mode, encoding=encoding, newline=newline)
File ~/checkouts/readthedocs.org/user_builds/festim-workshop/conda/festim1/lib/python3.11/site-packages/numpy/lib/_datasource.py:529, in DataSource.open(self, path, mode, encoding, newline)
526 return _file_openers[ext](found, mode=mode,
527 encoding=encoding, newline=newline)
528 else:
--> 529 raise FileNotFoundError(f"{path} not found.")
FileNotFoundError: task14/solute.txt not found.
The surface flux is computed for a full angular coverage by default. The units of the surface flux is H/m/s (per unit height).
plt.plot(
derived_quantities[0].t,
np.abs(derived_quantities[0].data),
linewidth=2,
label="Inner surface",
)
plt.plot(
derived_quantities[1].t,
np.abs(derived_quantities[1].data),
linewidth=2,
label="Outer surface",
)
plt.xlabel("Time [s]")
plt.ylabel("Absolute flux [H/m s]")
plt.legend()
plt.show()
Sphere#
The second example is a plain tungsten sphere of radius 0.1 mm. Similarly we can set the system of coordinates to spherical.
my_model = F.Simulation()
radius = 1e-4
my_model.mesh = F.MeshFromVertices(
np.linspace(0, radius, num=100),
type="spherical",
)
tungsten_D = htm.diffusivities.filter(material=htm.Tungsten).mean()
tungsten_S = htm.solubilities.filter(material=htm.Tungsten).mean()
tungsten = F.Material(
id=1,
D_0=tungsten_D.pre_exp.magnitude,
E_D=tungsten_D.act_energy.magnitude,
S_0=tungsten_S.pre_exp.magnitude,
E_S=tungsten_S.act_energy.magnitude,
)
my_model.materials = tungsten
my_model.boundary_conditions = [
F.SievertsBC(surfaces=[2], S_0=tungsten.S_0, E_S=tungsten.E_S, pressure=1e5),
]
my_model.T = 800
my_model.settings = F.Settings(
absolute_tolerance=1e10,
relative_tolerance=1e-10,
final_time=100,
)
my_model.dt = F.Stepsize(1, stepsize_change_ratio=1.1)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[8], line 9
5 np.linspace(0, radius, num=100),
6 type="spherical",
7 )
8
----> 9 tungsten_D = htm.diffusivities.filter(material=htm.Tungsten).mean()
10 tungsten_S = htm.solubilities.filter(material=htm.Tungsten).mean()
11 tungsten = F.Material(
12 id=1,
NameError: name 'htm' is not defined
The spherical equivalent class of SurfaceFlux is SurfaceFluxSpherical:
derived_quantities = F.DerivedQuantities(
[
F.SurfaceFluxSpherical(field="solute", surface=2),
],
show_units=True,
)
txt_export = F.TXTExport(field="solute", filename="task14/solute_spherical.txt")
my_model.exports = [derived_quantities, txt_export]
my_model.initialise()
my_model.run()
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[10], line 1
----> 1 my_model.initialise()
2 my_model.run()
File ~/checkouts/readthedocs.org/user_builds/festim-workshop/conda/festim1/lib/python3.11/site-packages/festim/generic_simulation.py:328, in Simulation.initialise(self)
324 set_log_level(self.log_level)
326 self.t = 0 # reinitialise t to zero
--> 328 if self.settings.chemical_pot:
329 self.mobile = festim.Theta()
330 else:
AttributeError: 'NoneType' object has no attribute 'chemical_pot'
import matplotlib.pyplot as plt
data = np.genfromtxt("task14/solute_spherical.txt", skip_header=1, delimiter=",")
cmap = plt.get_cmap("viridis")
for i in range(1, data.shape[1], 3):
plt.plot(
data[:, 0] / 1e-3, data[:, i], color=cmap(i / data.shape[1])
)
plt.xlabel("x [mm]")
plt.ylabel("c [H/m³]")
plt.show()
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
Cell In[11], line 3
1 import matplotlib.pyplot as plt
2
----> 3 data = np.genfromtxt("task14/solute_spherical.txt", skip_header=1, delimiter=",")
4 cmap = plt.get_cmap("viridis")
5 for i in range(1, data.shape[1], 3):
6 plt.plot(
File ~/checkouts/readthedocs.org/user_builds/festim-workshop/conda/festim1/lib/python3.11/site-packages/numpy/lib/_npyio_impl.py:1978, in genfromtxt(fname, dtype, comments, delimiter, skip_header, skip_footer, converters, missing_values, filling_values, usecols, names, excludelist, deletechars, replace_space, autostrip, case_sensitive, defaultfmt, unpack, usemask, loose, invalid_raise, max_rows, encoding, ndmin, like)
1976 fname = os.fspath(fname)
1977 if isinstance(fname, str):
-> 1978 fid = np.lib._datasource.open(fname, 'rt', encoding=encoding)
1979 fid_ctx = contextlib.closing(fid)
1980 else:
File ~/checkouts/readthedocs.org/user_builds/festim-workshop/conda/festim1/lib/python3.11/site-packages/numpy/lib/_datasource.py:192, in open(path, mode, destpath, encoding, newline)
155 """
156 Open `path` with `mode` and return the file object.
157
(...) 188
189 """
191 ds = DataSource(destpath)
--> 192 return ds.open(path, mode, encoding=encoding, newline=newline)
File ~/checkouts/readthedocs.org/user_builds/festim-workshop/conda/festim1/lib/python3.11/site-packages/numpy/lib/_datasource.py:529, in DataSource.open(self, path, mode, encoding, newline)
526 return _file_openers[ext](found, mode=mode,
527 encoding=encoding, newline=newline)
528 else:
--> 529 raise FileNotFoundError(f"{path} not found.")
FileNotFoundError: task14/solute_spherical.txt not found.
plt.plot(
derived_quantities[0].t,
np.abs(derived_quantities[0].data),
label="Outer surface flux",
)
plt.xlabel("Time [s]")
plt.ylabel("Absolute flux [H/s]")
plt.legend()
plt.show()
