Time-dependent source

User Support
1

Hi,

I guess this is a simple question, but couldn’t find a solution fast. If I have a source:

my_model.sources = [F.Source(value=1e20, volume=1, field=0)]

then I believe it should be possible for value to be time-dependent, like 1e20*festim.t. But what if I want to have non-smooth time dependence, e.g. source is on for one hour followed by another hour off and so on. Is it possible to introduce this kind of time dependence?

Best wishes,

Mikhail

2

Hi @miklavrent!

sympy.Piecewise can be used to set up a piece-wise source in FESTIM1. Consider the following example from the FESTIM V&V book:

ion_flux = sp.Piecewise((incident_flux, F.t <= 1), (0, True))

source_term = F.ImplantationFlux(
    flux=ion_flux, imp_depth=4.5e-9, width=2.5e-9, volume=1
)

If you use FESTIM2, then you can do it with ufl.conditional as shown in FESTIM2 V&V.

Please, let me know if it helps.

3

Thank you Vladimir,

Yes this works, like:

h_source = sp.Piecewise((1e20, F.t <= 3600), (0, F.t <= 7200), (1e20, F.t <= 10800), (0, True))

my_model.sources = [F.Source(value=h_source, volume=1, field=0)]

Just if I want to do this switch on and off for say a year, do you think there is a way to avoid writing a line for h_source with about 8760 conditions? :slight_smile:

Best wishes,

Mikhail

4

I remember having to write pulsed sources like the one you describe @miklavrent but I’m afraid I can’t remember how (it’s still early here in CA so it may come back later!!)

@VVKulagin do you know if Mikhail could use a modulo here?

5

@miklavrent @remidm

Yeah, how about modulo with some periodicity. Here is an example with ImplantationFlux:

def h_source(t):
    t1 = t % 7200  # period
    return sp.Piecewise((1e20, t1 < 3600), (0, True))


source_term = F.ImplantationFlux(
    flux=h_source(F.t), imp_depth=4.5e-9, width=2.5e-9, volume=1
)
Example from FESTIM simulation

image
image614×451 17.1 KB

Draft script 
import festim as F
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt

model = F.Simulation()

sample_depth = 5e-4

vertices = np.concatenate(
    [
        np.linspace(0, 30e-9, num=200),
        np.linspace(30e-9, 3e-6, num=300),
        np.linspace(3e-6, 20e-6, num=200),
        np.linspace(20e-6, sample_depth, num=100),
    ]
)

model.mesh = F.MeshFromVertices(vertices)
# Material Setup, only W
tungsten = F.Material(
    id=1,
    D_0=4.1e-07,  # m2/s
    E_D=0.39,  # eV
)

model.materials = tungsten


def h_source(t):
    t1 = t % 7200  # period
    return sp.Piecewise((1e20, t1 < 3600), (0, True))


source_term = F.ImplantationFlux(
    flux=h_source(F.t), imp_depth=4.5e-9, width=2.5e-9, volume=1
)

model.sources = [source_term]

# boundary conditions
model.boundary_conditions = [F.DirichletBC(surfaces=[1, 2], value=0, field=0)]


model.T = 293

model.dt = F.Stepsize(
    initial_value=0.1,
    stepsize_change_ratio=1.1,
    max_stepsize=100,
    dt_min=1e-05,
)

model.settings = F.Settings(
    absolute_tolerance=1e10,
    relative_tolerance=1e-09,
    final_time=18000,  # time to reach max temp
)

derived_quantities = F.DerivedQuantities(
    [
        F.HydrogenFlux(surface=1),
    ],
)

model.exports = [derived_quantities]

model.initialise()
model.run()

plt.plot(derived_quantities.t, -np.array(derived_quantities[0].data))
plt.ylabel("Hydrogen flux")
plt.xlabel("Time")
plt.yscale("log")
plt.show()
1 Like
6

Thank you both, yes this should work, also if on and off times are different, like

def h_source(t):
t1 = t % 7200 # period
return sp.Piecewise((1e20, t1 < 5400), (0, True))

Sorry, this is probably more a Python than a FESTIM question.

Best wishes,

Mikhail

7

Yes I think this should work all the same