Error defining time-dependent sources using dolfinx0.10.0

User Support
1

Hello everyone,

I am currently working with hisp and PFC-Tritium-Transport in order to determine tritium inventory in ITER. I am running a custom version of hisp (AdriaLlealS/hisp at fix-b-bins) using the festim-fenicsx branch of festim and dolfinx 0.10.0 version. I managed to run simulations without defining sources, and fixing the surface concentration of mobile species (using the implantation approximation Theory — FESTIM Documentation ).

However, when I try to define sources and use surface reactions as BCs (or setting surface concentrations to zero) I get an error that I haven’t been able to solve.

Here is a simple script relying only on festim replicating the error:

```python
import numpy as np
import ufl
import festim as F  # fenicsx branch (DOLFINx 0.10)
from dolfinx.fem.function import Constant
import numpy.typing as npt

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 = F.as_fenics_constant(self.flux(float(t)), mesh)
        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)

my_model = F.HydrogenTransportProblem()
folder = "outputs"
exports = True

L = 5e-4

def graded_vertices(L, h0, r):
    xs = [0.0]; h = h0
    while xs[-1] + h < L:
        xs.append(xs[-1] + h); h *= r
    if xs[-1] < L:
        xs.append(L)
    return np.array(xs)

vertices_graded = graded_vertices(L=L, h0=L/12e9, r=1.01)
my_model.mesh = F.Mesh1D(vertices_graded)

k_B = 8.617333262e-5
D_0 = 4.1e-7
E_D = 0.2
tungsten = F.Material(D_0=D_0, E_D=E_D, name="tungsten")

w_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, L], material=tungsten)
inlet = F.SurfaceSubdomain1D(id=1, x=0.0)
outlet = F.SurfaceSubdomain1D(id=2, x=L)
my_model.subdomains = [w_subdomain, inlet, outlet]

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

T_surface = 600.0
def temperature(t): return T_surface
my_model.temperature = temperature

implantation_range = 3e-9
width = 1e-9
def gaussian_distribution(x: npt.NDArray, mean: float, width: float) -> ufl.core.expr.Expr:
    return ufl.exp(-((x[0] - mean)**2)/(2*width**2)) / (ufl.sqrt(2*ufl.pi*width**2))

pulse_period = 100.0
pulse_on_time = 50.0
Gamma_ion_D_peak = 1e19
Gamma_ion_T_peak = 5e18

def square_pulse(t, peak):
    try:
        tau = float(t) % pulse_period
    except Exception:
        tau = t
    return peak if (0.0 <= tau < pulse_on_time) else 0.0

F_D = lambda t: square_pulse(t, Gamma_ion_D_peak)
F_T = lambda t: square_pulse(t, Gamma_ion_T_peak)
dist_fn = lambda x: gaussian_distribution(x, implantation_range, width)

my_model.sources = [
    PulsedSource(flux=F_D, distribution=dist_fn, volume=w_subdomain, species=mobile_D),
    PulsedSource(flux=F_T, distribution=dist_fn, volume=w_subdomain, species=mobile_T),
]

my_model.boundary_conditions = [
    F.FixedConcentrationBC(subdomain=inlet, value=0.0, species="D"),
    F.FixedConcentrationBC(subdomain=inlet, value=0.0, species="T"),
]

if exports:
    my_model.exports = [
        F.VTXSpeciesExport(f"{folder}/mobile_D.bp", field=mobile_D),
        F.VTXSpeciesExport(f"{folder}/mobile_T.bp", field=mobile_T),
        F.VTXTemperatureExport(f"{folder}/temperature.bp"),
    ]

quantities = {}
for sp in my_model.species:
    total = F.TotalVolume(field=sp, volume=w_subdomain)
    my_model.exports.append(total)
    quantities[sp.name] = total

final_time = 1000.0
my_model.settings = F.Settings(
    atol=1e10,
    rtol=1e-10,
    final_time=final_time,
    max_iterations=50,
)
my_model.settings.stepsize = F.Stepsize(initial_value=1e-3)

my_model.initialise()
my_model.run()
```

I included the PulsedSource class of hisp in the script. The error message hat I get is:

Traceback (most recent call last):
File “/home/ITER/llealsa/AdriaLlealS/test.py”, line 188, in
my_model.initialise()
File “/home/ITER/llealsa/.local/lib/python3.11/site-packages/festim/hydrogen_transport_problem.py”, line 326, in initialise
self.create_formulation()
File “/home/ITER/llealsa/.local/lib/python3.11/site-packages/festim/hydrogen_transport_problem.py”, line 871, in create_formulation
source.value.fenics_object
TypeError: unsupported operand type(s) for *: ‘NoneType’ and ‘Indexed’

I’ve been trying to fix this issue but haven’t had any success.

Thank you very much in advance and best regards,
Adrià

2

Hi @AdriaLlealS, I think your code block formatting needs fixing.

I’ll look into this issue

3

It’s very hard for me to run your code because of the formatting, can you try and fix it before I run it?

4

Yes of course, apologies for that. I think it should be easier now.

5

I managed to reproduce your error with the following minimal working example (only 60 line long!)

import numpy as np
import ufl
import festim as F
from dolfinx.fem.function import Constant
import numpy.typing as npt

print(F.__version__)


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 = F.as_fenics_constant(self.flux(float(t)), mesh)
        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)


my_model = F.HydrogenTransportProblem()
L = 5e-4

my_model.mesh = F.Mesh1D(np.linspace(0, L, 50))


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

w_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, L], material=tungsten)
my_model.subdomains = [w_subdomain]

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


my_model.temperature = 600.0


my_model.sources = [
    PulsedSource(
        flux=1, distribution=lambda x: 1, volume=w_subdomain, species=mobile_D
    ),
]


final_time = 1000.0
my_model.settings = F.Settings(
    atol=1e10,
    rtol=1e-10,
    final_time=final_time,
    max_iterations=50,
)
my_model.settings.stepsize = F.Stepsize(initial_value=1e-3)

my_model.initialise()
my_model.run()

@jhdark I believe this is because ParticleSource is now expecting a Value object. I’ll dig into this more

6

@AdriaLlealS the issue is that you pass None as the value argument when initialising the particle source parent class.

Instead, you can do the following:

import numpy as np
import festim as F

print(F.__version__)


class PulsedSource(F.ParticleSource):
    def __init__(self, flux, distribution, volume, species):
        self.flux = flux
        self.distribution = distribution
        value = lambda x, t: self.flux(t) * self.distribution(x)
        super().__init__(value, volume, species)


my_model = F.HydrogenTransportProblem()
L = 1

my_model.mesh = F.Mesh1D(np.linspace(0, L, 50))

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

w_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, L], material=tungsten)
boundary_left = F.SurfaceSubdomain1D(id=1, x=0.0)
boundary_right = F.SurfaceSubdomain1D(id=2, x=L)
my_model.subdomains = [w_subdomain, boundary_left, boundary_right]

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


my_model.temperature = 600.0

import ufl

my_model.sources = [
    PulsedSource(
        flux=lambda t: ufl.conditional(ufl.lt(t, 5.0), 0, 1.0),
        distribution=lambda x: ufl.conditional(ufl.lt(x[0], L / 2), 1.0, 0.0),
        volume=w_subdomain,
        species=mobile_D,
    ),
]

my_model.boundary_conditions = [
    F.FixedConcentrationBC(
        subdomain=boundary_left,
        species=mobile_D,
        value=0.0,
    ),
    F.FixedConcentrationBC(
        subdomain=boundary_right,
        species=mobile_D,
        value=0.0,
    ),
]

final_time = 10.0
my_model.settings = F.Settings(
    atol=1e-10,
    rtol=1e-10,
    final_time=final_time,
    max_iterations=50,
)
my_model.settings.stepsize = F.Stepsize(initial_value=1)

my_model.exports = [F.VTXSpeciesExport(filename="out.bp", field=mobile_D)]

my_model.initialise()
my_model.run()

In this example, the flux value is 0 before t=5 and the flux distribution is 1 on the left half and 0 otherwise.
Concentration for t < 5

image
image717×458 4.26 KB

Concentration at t=10

image
image717×455 5.47 KB

7

You can even further simplify your model by simpy doing:

def flux(t):
    return ufl.conditional(ufl.lt(t, 5.0), 0, 1.0)


def distribution(x):
    return ufl.conditional(ufl.lt(x[0], L / 2), 1.0, 0.0)


my_model.sources = [
    F.ParticleSource(
        value=lambda x, t: flux(t) * distribution(x),
        volume=w_subdomain,
        species=mobile_D,
    ),
]
8

I see, I am trying to create the flux and distribution functions as ufl conditionals now.
Thank you very much for your help!

1 Like
9

Yes for now ufl.conditional objects are easier to work with (inside FESTIM) than python functions. Though I understand python functions are often more convenient for users…