Pulsed source

[RemDelaporteMathurin]

@kaelyndunnell this is how you can have a pulsed particle source in FESTIM.

I made a simple example with 1 material, one species, no traps...

import festim as F
import numpy as np

from dolfinx.fem.function import Constant

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

tungsten = F.Material(D_0=1, E_D=0)

w_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=tungsten)
left = F.SurfaceSubdomain1D(id=1, x=0)
right = F.SurfaceSubdomain1D(id=2, x=1)

my_model.subdomains = [w_subdomain, left, right]

mobile_D = F.Species("D")
my_model.species = [mobile_D]

my_model.temperature = 1000


def gaussian_distribution(x):
    # depth = 3e-9
    # width = 1e-9
    # return ufl.exp(-((x[0] - depth) ** 2) / (2 * width**2))
    return 1


import ufl

# all in seconds
flat_top_duration = 10
ramp_up_duration = 2
ramp_down_duration = 2
dwelling_time = 10

total_time_pulse = flat_top_duration + ramp_up_duration + ramp_down_duration
total_time_cycle = total_time_pulse + dwelling_time


def deuterium_flux(t: Constant):
    flat_top_value = 1
    resting_value = 0
    return (
        flat_top_value
        if float(t) % total_time_cycle < total_time_pulse
        else resting_value
    )

my_model.sources = [
    PulsedSource(
        flux=deuterium_flux,
        distribution=gaussian_distribution,
        species=mobile_D,
        volume=w_subdomain,
    )
]

my_model.boundary_conditions = [
    F.DirichletBC(subdomain=left, value=0, species=mobile_D),
    F.DirichletBC(subdomain=right, value=0, species=mobile_D),
]

folder = "pulse_example"

my_model.exports = [
    F.VTXExport(f"{folder}/mobile_concentration_d.bp", field=mobile_D),
    F.TotalVolume(field=mobile_D, volume=w_subdomain),
]

my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=100)

my_model.settings.stepsize = F.Stepsize(initial_value=1)

my_model.initialise()
my_model.run()

import matplotlib.pyplot as plt

inv = my_model.exports[1]

plt.plot(inv.t, inv.data)
plt.show()

We just have to define our own class PulsedSource that takes two inputs flux and distribution and both can be callables.

class PulsedSource(F.ParticleSource):
    def __init__(self, flux, distribution, volume, species):
        self.flux = flux
        self.distribution = distribution
        super().__init__(None, volume, species)

    @property
    def time_dependent(self):
        return True

    def create_value_fenics(self, mesh, temperature, t: Constant):
        self.flux_fenics = Constant(mesh, float(self.flux(t)))
        x = ufl.SpatialCoordinate(mesh)
        self.distribution_fenics = self.distribution(x)

        self.value_fenics = self.flux_fenics * self.distribution_fenics

    def update(self, t: float):
        self.flux_fenics.value = self.flux(t)

Produces: