Exporting fields

FESTIM has convenience classes that allow users to create XDMF and VTX exports, which can then be viewed in Paraview.

Objectives:

• Writing VTX files for species and temperature fields

• Additional options for VTX exports

• Writing XDMF files for field exports

• Exporting a custom field

• Exporting a reaction rate

Writing VTX files for species and temperature fields

Users can export concentration fields to VTX using VTXSpeciesExport, and then view their exports in ParaView. This example will discuss how to define the export, and what result should you expect to see.

Let us setup a 2D, transient problem with the following boundary conditions:

\[ \text{Right surface:} \quad c_H = 0 \]

\[ \text{Left surface:} \quad c_H = 1 \]

Show code cell source

Hide code cell source

import numpy as np
import festim as F
from dolfinx.mesh import create_unit_square
from mpi4py import MPI

mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
my_model = F.HydrogenTransportProblem()
my_model.mesh = F.Mesh(mesh)
mat = F.Material(D_0=1e-2, E_D=0)

right_surface = F.SurfaceSubdomain(id=1, locator = lambda x: np.isclose(x[0], 1.0))
left_surface = F.SurfaceSubdomain(id=2, locator = lambda x: np.isclose(x[0], 0.0))
vol = F.VolumeSubdomain(id=1, material=mat)

H = F.Species("H")

my_model.subdomains = [right_surface, left_surface, vol]
my_model.species = [H]
my_model.boundary_conditions = [
    F.FixedConcentrationBC(subdomain=right_surface, value=0, species=H),
    F.FixedConcentrationBC(subdomain=left_surface, value=1, species=H)
]
my_model.temperature = lambda x: 400 + 10*x[0] + 2*x[1]
my_model.settings = F.Settings(atol=1e-10,rtol=1e-10,stepsize=1, final_time=10)

We can export the concentration field by defining a VTXSpeciesExport object for our model. The required arguments are a filename path (which must end in .bp) and the field (species) you want to export:

H_export = F.VTXSpeciesExport(filename="H_concentration.bp", field=H)

Note

If we had several Species, we would create several VTXSpeciesExport objects, one per species.

If needed, we can also export the temperature field this way using VTXTemperatureExport

temp_export = F.VTXTemperatureExport(filename="temperature.bp")

Then, we just pass these exports to my_model.exports as a list:

my_model.exports = [H_export, temp_export]

my_model.initialise()
my_model.run()

We should expect to see two new folders called H_concentration.bp and temperature.bp. To view the results, we can use ParaView (see the ParaView section to learn more).

Exporting results for a discontinuous problem

If running a multi-material discontinuous simulation, it is necessary to also specify the subdomain during the export process.

Consider the same multi-material problem in the Materials chapter, where we use HydrogenTransportProblemDiscontinuous to solve a multi-material problem:

Show code cell source

Hide code cell source

import festim as F
import numpy as np
from dolfinx.mesh import create_unit_square
from mpi4py import MPI

fenics_mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
festim_mesh = F.Mesh(fenics_mesh)

my_model = F.HydrogenTransportProblemDiscontinuous()

material_top = F.Material(D_0=1, E_D=0, K_S_0=2, E_K_S=0)
material_bottom = F.Material(D_0=2, E_D=0, K_S_0=3, E_K_S=0)

top_volume = F.VolumeSubdomain(id=3, material=material_top, locator=lambda x: x[1] >= 0.5)
bottom_volume = F.VolumeSubdomain(id=4, material=material_bottom, locator=lambda x: x[1] <= 0.5)

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))

my_model.mesh = festim_mesh
my_model.subdomains = [top_surface, bottom_surface, top_volume, bottom_volume]

my_model.interfaces = [F.Interface(5, (bottom_volume, top_volume), penalty_term=1000)]
my_model.surface_to_volume = {
    top_surface: top_volume,
    bottom_surface: bottom_volume,
}

H = F.Species("H")
my_model.species = [H]
H.subdomains = [top_volume, bottom_volume]

my_model.temperature = 400

my_model.boundary_conditions = [
    F.FixedConcentrationBC(subdomain=top_surface, value=1.0, species=H),
    F.FixedConcentrationBC(subdomain=bottom_surface, value=0.0, species=H),
]

my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False)

We can specify separate export objects for the top and bottom domains using the subdomain argument, and should expect to see two new folders created named top.bp and bottom.bp:

top_export = F.VTXSpeciesExport(filename="top.bp", field=H, subdomain=top_volume)
bottom_export = F.VTXSpeciesExport(filename="bottom.bp", field=H, subdomain=bottom_volume)
my_model.exports = [
    top_export,
    bottom_export,
]
my_model.initialise()
my_model.run()

Finally, these .bp exports can be viewed in ParaView:

Note

For multi-material discontinuous problems, each .bp file only shows its corresponding subdomain.

Exporting fields at specific timesteps

Users can also specify which timesteps they’d like to export using the times argument (which must be a list):

export = F.VTXSpeciesExport(filename="H_concentration.bp", field=H, times=[0, 5, 10])

If no times argument is given, the export stores results for all timesteps by default.

Checkpointing

It may be helpful to store results from one simulation for later use in another (perhaps as an initial condition, see Initial condition from a checkpoint file). FESTIM includes this capability by incorporationg io4dolfinx functionality, which stores mesh information and solutions into a checkpoint.bp file. Learn more about checkpointing in DOLFINx here.

To store the species field as a checkpoint file, simply set the checkpoint argument to True:

export = F.VTXSpeciesExport(filename="H_concentration.bp", field=H, checkpoint=True)

Note

Checkpointed files cannot be viewed in ParaView.

Writing XDMF files for field exports

Warning

Exporting to VTX is preferable over XDMF, as XDMF functionality will soon be deprecated. Additionally, you cannot view transient results using XDMF.

Users can export functions to XDMF files using the XDMFExport class, which requires a filename and field:

import festim as F

H = F.Species("H")
export = F.XDMFExport(filename="my_export.xdmf", field=H)

To export this in a FESTIM simulation, add the export to your problem’s export attribute:

my_model = F.HydrogenTransportProblem()
my_model.exports = [export]

This will produce the corresponding export files (my_export.xdmf and my_export.h5).

Exporting a custom field

Added in version 2.0: CustomFieldExport was introduced in FESTIM 2.0.

Sometimes the field you want to visualise is not a species concentration or the temperature directly, but a quantity derived from them. CustomFieldExport writes an arbitrary field to a VTX file. You pass an expression (a callable) whose positional arguments can be t (time), x (spatial coordinate), T (temperature), or any species mapped through species_dependent_value.

As an example, let us export the hydrogen concentration converted from \(\mathrm{mol/m^3}\) to \(\mathrm{atoms/m^3}\) by multiplying by Avogadro’s constant:

Show code cell source

Hide code cell source

import numpy as np
import festim as F
from dolfinx.mesh import create_unit_square
from mpi4py import MPI

my_model = F.HydrogenTransportProblem()
my_model.mesh = F.Mesh(create_unit_square(MPI.COMM_WORLD, 10, 10))
mat = F.Material(D_0=1e-2, E_D=0)

right_surface = F.SurfaceSubdomain(id=1, locator=lambda x: np.isclose(x[0], 1.0))
left_surface = F.SurfaceSubdomain(id=2, locator=lambda x: np.isclose(x[0], 0.0))
vol = F.VolumeSubdomain(id=1, material=mat)

H = F.Species("H")

my_model.subdomains = [right_surface, left_surface, vol]
my_model.species = [H]
my_model.boundary_conditions = [
    F.FixedConcentrationBC(subdomain=right_surface, value=0, species=H),
    F.FixedConcentrationBC(subdomain=left_surface, value=1, species=H),
]
my_model.temperature = 400
my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, stepsize=1, final_time=5)

NA = 6.022e23  # Avogadro's constant

atoms_export = F.CustomFieldExport(
    filename="H_atoms.bp",
    expression=lambda H: H * NA,
    species_dependent_value={"H": H},
)

my_model.exports = [atoms_export]
my_model.initialise()
my_model.run()

This produces a H_atoms.bp folder that can be opened in ParaView. The computed field is also available on the export object through its function attribute:

print(f"Max concentration in atoms/m3: {atoms_export.function.x.array.max():.3e}")

Max concentration in atoms/m3: 6.022e+23

Note

The keys of species_dependent_value must match the argument names of expression. To combine several species — for example to export the total hydrogen isotope concentration \(c_\mathrm{H} + c_\mathrm{D}\) — map each argument to its Species: species_dependent_value={"cH": H, "cD": D} with expression=lambda cH, cD: cH + cD.

Exporting a reaction rate

Added in version 2.0: ReactionRateExport was introduced in FESTIM 2.0.

When a model contains reactions, it is often useful to visualise where and how fast a reaction proceeds. ReactionRateExport writes the rate of a given Reaction to a VTX file. The direction argument selects the "forward", "backward", or "both" (net) contribution.

Here we set up a 1D trapping problem where mobile hydrogen is trapped, and export the net trapping rate:

Show code cell source

Hide code cell source

my_model = F.HydrogenTransportProblem()
my_model.mesh = F.Mesh1D(np.linspace(0, 1, 100))

mobile_H = F.Species("H")
trapped_H = F.Species("H_trapped", mobile=False)
empty_traps = F.ImplicitSpecies(n=2, others=[trapped_H], name="empty_traps")
my_model.species = [mobile_H, trapped_H]

material = F.Material(D_0=1, E_D=0)
vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=material)
left_surf = F.SurfaceSubdomain1D(id=1, x=0)
right_surf = F.SurfaceSubdomain1D(id=2, x=1)
my_model.subdomains = [vol, left_surf, right_surf]

my_model.boundary_conditions = [
    F.FixedConcentrationBC(left_surf, value=10, species=mobile_H),
    F.FixedConcentrationBC(right_surf, value=0, species=mobile_H),
]
my_model.temperature = 300
my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=50)
my_model.settings.stepsize = F.Stepsize(1)

We keep a reference to the Reaction object so we can pass it to the export:

trapping_reaction = F.Reaction(
    reactant=[mobile_H, empty_traps],
    product=[trapped_H],
    k_0=0.01,
    E_k=0,
    p_0=0.1,
    E_p=0,
    volume=vol,
)
my_model.reactions = [trapping_reaction]

rate_export = F.ReactionRateExport(
    reaction=trapping_reaction,
    filename="trapping_rate.bp",
    direction="both",
)

my_model.exports = [rate_export]
my_model.initialise()
my_model.run()

The trapping rate field is written to trapping_rate.bp (viewable in ParaView) and is also accessible through the function attribute:

print(f"Max trapping rate: {rate_export.function.x.array.max():.3e}")

Max trapping rate: 1.408e-04