Combine Fenics with BDF of Scipy

1

Hello dear member,
I’m in the process of solving a broadcast reaction problem but I’m encountering a time step problem. When I have an extremely small time step, the numerical error seems to be very small, but I’m looking for a more appropriate and flexible solution. Since Fenics doesn’t give a method for the time step (BDF, Runge-Kutta), I thought of discretising the problem with fenics and solving it directly with scipy’s solve_ivp function, which implements a BDF. After doing this, I get an error that I can’t solve.
Has anyone used such a strategy before? If you have any other methods, I’d love to hear about them.
I’m using FEniCs Legacy. Here is some minimal code for the Aliev-Panfilov model.

import dolfin as df
import numpy as np
from scipy.integrate import solve_ivp

# Constants Parameters
alpha = df.Constant(0.01)
b = df.Constant(0.15)
D = df.Constant(0.45)
c = df.Constant(8)
gamma = df.Constant(0.002)
mu_1 = df.Constant(0.2)
mu_2 = df.Constant(0.3)

# Domaine spatial
mesh = df.RectangleMesh(df.Point(0, -5), df.Point(200, 5), 100, 20)

# Finite Element space 
P1 = df.FiniteElement("CG", mesh.ufl_cell(), 1)
element = df.MixedElement([P1, P1])
V = df.FunctionSpace(mesh, element)

# Conditions aux limites 
bc = df.DirichletBC(V, df.Constant((0, 0)), "on_boundary")

# Définition of the test en trial Function
v_phi, v_r = df.TestFunctions(V)
phi_r = df.Function(V)
phi, r = df.split(phi_r)
# Definition of the initials conditions 
Mh = df.FunctionSpace(mesh, P1)
phi0 = df.interpolate(df.Constant(0.0),Mh)
r0 = df.interpolate(df.Constant(0.0),Mh)


waveS0 = df.Expression("1*(x[0] < 20)", amp = 1, degree=2)
# time Discretisation (BDF)
t_span = (0, 10)  # Plage de temps
t_eval = np.linspace(t_span[0], t_span[1], num=1000)  # Évaluation temporelle
dt = t_eval[1] - t_eval[0]

# Source I(t)
I = df.Expression("t <= (dt + 1) && x[0] <= 10 ? 1 : 0", t=0, dt=dt, degree=2)
#t = df.Constant(0)
# variationnal Formulation 
F_phi = (
    (phi - phi0) * v_phi * df.dx
    + dt * (
        D * df.dot(df.grad(phi), df.grad(v_phi)) * df.dx
        + c * phi * (phi - alpha) * (1 - phi) * v_phi * df.dx
        - r * phi * v_phi * df.dx
        + I * v_phi * df.dx
    )
)
F_r = (
    (r - r0) * v_r * df.dx
    + dt
    * (
        (
            gamma
            + (r * mu_1) / (mu_2 + phi)
        )
        * (r + c * phi * (phi - b - 1))
        * v_r
        * df.dx
    )
)

# Residual for Scipy (BDF)
def residual(t, phi_r_vec):
    phi_r.vector().set_local(phi_r_vec)
    df.solve(F_phi == 0, phi_r, bcs=bc)
    df.solve(F_r == 0, phi_r, bcs=bc)
    return df.assemble(F_phi), df.assemble(F_r)

# solving
sol = solve_ivp(
    residual,
    t_span,
    phi_r.vector().get_local(),
    t_eval,
    method="BDF",
    jac=None,
    rtol=1e-6,
    atol=1e-6,
)

# Obtention of the solutions
solutions = [df.Function(V, phi_r.vector()) for phi_r in sol.y]

I have the following error.

Traceback (most recent call last):
  File "testFS3.py", line 84, in <module>
    atol=1e-6,
TypeError: solve_ivp() got multiple values for argument 'method'

Thank you in advance for your Helps.

1 Like
2

The problem is not related to FEniCS. You’re giving positional arguments to solve_ivp, which takes method as 4th positional argument while you’re also giving method=“BDF”. In essence, you’re first trying to give t_eval as method and later trying to overwrite it with method="BDF”. You should use all keyword argument or positional argument not to confuse them. scipy.integrate.solve_ivp — SciPy v1.11.2 Manual

2 Likes
3

Thanks for your reply.
Could you please advise me on another technique to apply this strategy? So discretize with FEniCs and solve with an RK45 or a BDF.