Initial conditions

book/_toc.yml

@@ -26,6 +26,11 @@ parts:
     - file: content/species_reactions/species
     - file: content/species_reactions/reactions
     - file: content/species_reactions/surface_reactions
+
+  - caption: Initial Conditions
+    chapters:
+    - file: content/initial_conditions/concentration
+    - file: content/initial_conditions/temperature
 
   - caption: Temperature
     chapters:

book/content/initial_conditions/concentration.md

@@ -0,0 +1,174 @@
+---
+jupytext:
+  formats: ipynb,md:myst
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.18.1
+kernelspec:
+  display_name: festim-workshop
+  language: python
+  name: python3
+---
+
+# Concentration #
+
+Initial conditions are an important part of transient simulations for both hydrogen transport and heat-transfer problems. FESTIM defaults to zero values for initial conditions, but can be set using FESTIM's `InitialConditionBase`. 
+
+This tutorial discusses how to set initial conditions for particle concentrations.
+
+Objectives:
+* Setting the initial concentration for a single species
+* Defining initial concentration conditions for multiple species
+
++++
+
+## Setting the initial concentration for a single species ##
+
+The `InitialConcentration` class can be used for defining initial conditions for the concentrations, which must defined on a volume subdomain in a FESTIM simulation:
+
+```{code-cell} ipython3
+import festim as F
+
+material = F.Material(D_0=1, E_D=0)
+vol = F.VolumeSubdomain(id=1, material=material)
+H = F.Species("H")
+
+IC = F.InitialConcentration(value=2, species=H, volume=vol)
+```
+
+Initial conditions can also be a function of space `x` and temperature `T`, such as:
+
+$$
+
+\text{IC(x, y, T)} = 2x + 3y + 5T
+
+$$
+
+```{code-cell} ipython3
+my_function = lambda x,T: 2*x[0] + 3*x[1] + 5*T
+IC = F.InitialConcentration(value=my_function, species=H, volume=vol)
+```
+
+## Multi-species initial conditions ##
+
++++
+
+The same class `InitialConcentration` can be used for multiple species, but a separate initial condition is required for each species.
+
+Consider the following 2D example, where we have three species (hydrogen, deuterium, and tritium) with initial and boundary conditions for species' concentrations:
+
+$$ \text{IC (hydrogen)} = 4 $$
+
+$$ \text{IC (deuterium)} = 5 $$
+
+$$ \text{IC (tritium)} = 6 $$
+
+$$ \text{Concentration BC (left, right, top) = 0} $$
+
+First, we can create a 2D unit square mesh and mark the subdomains:
+
+```{code-cell} ipython3
+import festim as F
+from dolfinx.mesh import create_unit_square
+from mpi4py import MPI
+import numpy as np
+
+nx, ny = 50, 50
+mesh = create_unit_square(MPI.COMM_WORLD, nx, ny)
+
+model = F.HydrogenTransportProblem()
+model.mesh = F.Mesh(mesh)
+
+material = F.Material(D_0=1e-2, E_D=0)
+
+top_surface = F.SurfaceSubdomain(id=1, locator = lambda x: np.isclose(x[1], 1.0))
+bottom_surface = F.SurfaceSubdomain(id=2, locator = lambda x: np.isclose(x[1], 0.0))
+right_surface = F.SurfaceSubdomain(id=3, locator = lambda x: np.isclose(x[0], 1.0))
+left_surface = F.SurfaceSubdomain(id=4, locator = lambda x: np.isclose(x[0], 0.0))
+
+vol = F.VolumeSubdomain(id=1, material=material)
+
+model.subdomains = [top_surface, bottom_surface, left_surface, right_surface, vol]
+```
+
+Now, let's define our species and initial conditions. We must create an `InitialConcentration` for each species, then we can pass the initial conditions into the `initial_conditions` attribute:
+
+```{code-cell} ipython3
+H = F.Species("H")
+D = F.Species("D")
+T = F.Species("T")
+
+IC_H = F.InitialConcentration(value=4, species=H, volume=vol)
+IC_D = F.InitialConcentration(value=5, species=D, volume=vol)
+IC_T = F.InitialConcentration(value=6, species=T, volume=vol)
+
+ICs = [IC_H, IC_D, IC_T]
+model.species = [H, D, T]
+
+model.initial_conditions = ICs
+```
+
+```{note}
+To pass the initial conditions into FESTIM, you must use a list!
+```
+
++++
+
+We also define the zero concentration boundary conditions for each species. Let's run the simulation for 5 seconds with a stepsize of 1 second, and visualize the final concentration profiles (hydrogen on the left, deuterium in the middle, and tritium on the right):
+
+```{code-cell} ipython3
+model.boundary_conditions = [
+    F.FixedConcentrationBC(subdomain=left_surface, value=0, species=H),
+    F.FixedConcentrationBC(subdomain=right_surface, value=0, species=H),
+    F.FixedConcentrationBC(subdomain=top_surface, value=0, species=H),
+    F.FixedConcentrationBC(subdomain=left_surface, value=0, species=D),
+    F.FixedConcentrationBC(subdomain=right_surface, value=0, species=D),
+    F.FixedConcentrationBC(subdomain=top_surface, value=0, species=D),
+    F.FixedConcentrationBC(subdomain=left_surface, value=0, species=T),
+    F.FixedConcentrationBC(subdomain=right_surface, value=0, species=T),
+    F.FixedConcentrationBC(subdomain=top_surface, value=0, species=T),
+]
+model.temperature = 400
+model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=50, stepsize=5)
+
+model.initialise()
+model.run()
+```
+
+```{code-cell} ipython3
+:tags: [hide-input]
+
+import pyvista
+from dolfinx import plot
+
+def make_ugrid(solution):
+    topology, cell_types, geometry = plot.vtk_mesh(solution.function_space)
+    u_grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
+    u_grid.point_data["c"] = solution.x.array.real
+    u_grid.set_active_scalars("c")
+    return u_grid
+
+pyvista.start_xvfb()
+pyvista.set_jupyter_backend("html")
+
+H_grid = make_ugrid(H.post_processing_solution)
+D_grid = make_ugrid(D.post_processing_solution)
+T_grid = make_ugrid(T.post_processing_solution)
+
+pl = pyvista.Plotter(shape=(1, 3))
+pl.subplot(0, 0)
+_ = pl.add_mesh(H_grid,label='h')
+pl.view_xy()
+
+pl.subplot(0, 1)
+_ = pl.add_mesh(D_grid,label='d')
+pl.view_xy()
+
+pl.subplot(0, 2)
+_ = pl.add_mesh(T_grid,label='t')
+pl.view_xy()
+
+pl.show()
+```

book/content/initial_conditions/temperature.md

@@ -0,0 +1,134 @@
+---
+jupytext:
+  formats: ipynb,md:myst
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.18.1
+kernelspec:
+  display_name: festim-workshop
+  language: python
+  name: python3
+---
+
+# Temperature #
+
+This tutorial discusses how to set initial conditions in FESTIM heat-transfer problems using the `InitialConditionBase`.
+
+Objectives:
+* Defining temperature initial conditions
+* Solving a heat-transfer problem with initial conditions
+
++++
+
+## Defining temperature initial conditions ##
+
+Similar to the __[concentration](https://festim-workshop.readthedocs.io/en/festim2/content/initial_conditions/concentration.html)__, we can define temperature initial conditions on volume subdomains using `InitialTemperature`:
+
+```{code-cell} ipython3
+import festim as F
+
+material = F.Material(D_0=1, E_D=0, thermal_conductivity=1, heat_capacity=2,density=1)
+vol = F.VolumeSubdomain(id=1, material=material)
+
+IC = F.InitialTemperature(value=400, volume=vol)
+```
+
+```{note}
+Since initial conditions are only used in transient simulations, `thermal_conductivity`, `heat_capacity`, and `density` must be defined for the material. Learn more about defining thermal properties __[here](https://festim-workshop.readthedocs.io/en/festim2/content/material/material_basics.html#defining-thermal-properties)__.
+```
+
++++
+
+Temperature initial conditions can be also be a function of space `x`, such as:
+
+$$ \text{IC(x, y ,z)} = 2\sin(x) + 4 \cos(y) - 2 \exp(-z) $$
+
+```{code-cell} ipython3
+import numpy as np
+
+my_temperature = lambda x: 2*np.sin(x[0]) + 4*np.cos(x[1]) - 2*np.exp(-x[2])
+
+IC = F.InitialTemperature(value=my_temperature, volume=vol)
+```
+
+## Solving a heat-transfer problem with initial conditions ##
+
+Consider the following 2D heat transfer problem with boundary and initial conditions:
+
+$$ \text{IC} = 300 \text{K} $$
+
+$$ \text{Temperature BC (left)} = 400 \text{K} $$
+
+$$ \text{Temperature BC (right)} = 350 \text{K} $$
+
+Let us first create our mesh and subdomains:
+
+```{code-cell} ipython3
+import festim as F
+from dolfinx.mesh import create_unit_square
+from mpi4py import MPI
+import numpy as np
+
+nx, ny = 50, 50
+mesh = create_unit_square(MPI.COMM_WORLD, nx, ny)
+
+model = F.HeatTransferProblem()
+model.mesh = F.Mesh(mesh)
+
+material = F.Material(D_0=1, E_D=0, thermal_conductivity=1, heat_capacity=2,density=1)
+
+top_surface = F.SurfaceSubdomain(id=1, locator = lambda x: np.isclose(x[1], 1.0))
+bottom_surface = F.SurfaceSubdomain(id=2, locator = lambda x: np.isclose(x[1], 0.0))
+right_surface = F.SurfaceSubdomain(id=3, locator = lambda x: np.isclose(x[0], 1.0))
+left_surface = F.SurfaceSubdomain(id=4, locator = lambda x: np.isclose(x[0], 0.0))
+
+vol = F.VolumeSubdomain(id=1, material=material)
+
+model.subdomains = [top_surface, bottom_surface, left_surface, right_surface, vol]
+```
+
+We can pass our initial conditions using the `initial_conditions` attribute, which requires a list. Let's run the simulation for 5 seconds and see the results:
+
+```{code-cell} ipython3
+model.initial_conditions = [
+    F.InitialTemperature(value=300, volume=vol)
+]
+
+model.boundary_conditions = [
+    F.FixedTemperatureBC(subdomain=left_surface, value=400),
+    F.FixedTemperatureBC(subdomain=right_surface, value=350),
+]
+
+model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=5, stepsize=1)
+
+model.initialise()
+model.run()
+```
+
+```{code-cell} ipython3
+:tags: [hide-input]
+
+import pyvista
+from dolfinx import plot
+
+pyvista.start_xvfb()
+pyvista.set_jupyter_backend("html")
+
+T = model.u
+
+topology, cell_types, geometry = plot.vtk_mesh(T.function_space)
+u_grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
+u_grid.point_data["T"] = T.x.array.real
+u_grid.set_active_scalars("T")
+u_plotter = pyvista.Plotter()
+u_plotter.add_mesh(u_grid, cmap="inferno", show_edges=False)
+
+u_plotter.view_xy()
+
+if not pyvista.OFF_SCREEN:
+    u_plotter.show()
+else:
+    figure = u_plotter.screenshot("temperature.png")
+```

environment.yml

@@ -12,7 +12,7 @@ dependencies:
   - ipykernel
   - nbconvert
   - fenics-dolfinx
-  - festim=2.0a7
+  - festim=2.0a8
   - python-gmsh
   - shapely
   - scipy