Custom boundary condition question

Hello! Is it possible to set the boundary condition as a function of time in festim? I want to modelling TDS-spectra with that condition: if time >= 100 (for example), boundary condition set to RecombinationFlux and if time < 100, boundary condition set to DirichletBC. I tried to use CustomFlux in that way:

def Kr_t(test_time):
    Kr_0 = 1e-28,
    Kr_E = -0.27,
    if test_time >= 100:
        Flux_t = -Kr_0*sp.exp(-Kr_E/8.625e-5/my_model.T)*(F.solute)^2, #recombination
    else:
        Flux_t = 0 #Dirichlet with value = 0
    return Flux_t
my_model.boundary_conditions = F.CustomFlux(
    field = 0,
    surfaces = [1],
    function = Kr_t,
    test_time = my_model.t,
)

but it gives me an error “stepsize reached minimal value”. If I set BC to only RecombinationFlux(1e-28, -0.27, 2, surfaces = [1]) or DirichletBC(surfaces=[1], value=0, field=0), it works.

Hi @cruiseraurora!

Is it possible to set the boundary condition as a function of time in festim?

Yes, it is. You can do it with the class in your example, but notice that the first 2 arguments for F.CustomFlux are T and solute fields at the boundary, as stated here. The API docs also show an example of usage:

def fun(T, solute, param1):
    return 2*T + solute - param1

my_bc = CustomFlux(
    surfaces=[1, 2],
    function=fun,
    field=0,
    param1=2*festim.x + festim.t
)

I’m not 100% confident that if-clause would work properly if you want to define a piece-wise function. Probably, you might need to use sp.Piecewise or fenics.conditional, but it requires testing.

You mean that instead of def Kr_t(test_time) I should try def Kr_t(T,solute,test_time)?

Now my code looks like this:

def Kr_t(T,solute,test_time):
    Kr_0 = 1e-28,
    Kr_E = -0.27,
    Flux_t = sp.Piecewise(
        (-Kr_0*sp.exp(-Kr_E/8.625e-5/T)*(solute)^2, test_time >= implantation_time+cold_time),
        (0, True))
    return Flux_t
my_model.boundary_conditions = F.CustomFlux(
    field = 0,
    surfaces = [1],
    function = Kr_t,
    test_time = my_model.t,
)

but I still have an error “stepsize reached minimal value”.

Yeah, I meant to use the custom flux BC as you showed, but exponentiation should be made with ** rather than the XOR operator ^.

Additionally, it’s difficult to suggest anything as you don’t provide the details about when the error appears (at the start or during the computation), what are the stepsize function, mesh, etc.

@cruiseraurora if I understand correctly you don’t want to change the value of the flux but rather change the type of the BC in the middle of the simulation, from a recombination BC (Robin type) to a fixed concentration of zero (homogeneous Dirichlet BC). Correct?

This specific case is not trivial in FESTIM. Here are your options:

  • split your problem in two simulations: the first with the RecombinationFlux, the second simulation with the Dirichlet BC
  • instead of changing the BC type, consider artificially increasing the Kr value so that the concentration is almost zero (ie. instantaneous recombination)

As of the exact reason of the stepsize reached minimal value error, this is due to your problem not converging and it hard to debug without a minimal working example (MWE) reproducing the issue

The error appears after 1.7e1 sec - console writes this
0.0 % 1.7e+01 s Elapsed time so far: 11.7
the computation slows donw and after a few seconds it prints this

Traceback (most recent call last):time so far: 13.9 s
  File "/home/wsluser/FESTIMfolder/Rusferplasma1e25.py", line 166, in <module>
    my_model.run()
  File "/home/wsluser/anaconda3/envs/festim-env/lib/python3.11/site-packages/festim/generic_simulation.py", line 422, in run
    self.run_transient()
  File "/home/wsluser/anaconda3/envs/festim-env/lib/python3.11/site-packages/festim/generic_simulation.py", line 445, in run_transient
    self.iterate()
  File "/home/wsluser/anaconda3/envs/festim-env/lib/python3.11/site-packages/festim/generic_simulation.py", line 473, in iterate
    self.h_transport_problem.update(self.t, self.dt)
  File "/home/wsluser/anaconda3/envs/festim-env/lib/python3.11/site-packages/festim/h_transport_problem.py", line 338, in update
    dt.adapt(t, nb_it, converged)
  File "/home/wsluser/anaconda3/envs/festim-env/lib/python3.11/site-packages/festim/stepsize.py", line 101, in adapt
    raise ValueError("stepsize reached minimal value")
ValueError: stepsize reached minimal value

I try to simplify CustomFlux as much as possible and make it look like DirichletBC:

def Kr_t(T,solute,test_time):
    Flux_t = 0
    return Flux_t
my_model.boundary_conditions = F.CustomFlux(
    field = 0,
    surfaces = [1],
    function = Kr_t,
)

but it gives me the same error (“stepsize reached minimal value”) at the same time (1.7e1 sec). If I change BC from CustomFlux to DirichletBC with arguments (surfaces=[1], value=0, field=0), DirichletBC causes no troubles in my computations. Why? I really apologize for my misunderstanding of this situation…

@cruiseraurora careful here, setting the flux to zero is very different from setting the concentration to zero.

Constraining the particle flux (Flux BC) enforces the -D \nabla c \cdot n value at the boundary whereas a DirichletBC enforces the c value at the boundary