Solving Matrix-Vector Equation for Coupled Nonlinear System

General
1

Hi everyone,

I’m using FEniCS 2019.1.0 to solve the following reaction-diffusion system -

\frac{\partial u_1}{\partial t} = D_1\nabla^2 u_1 - Ku_1u_2 + f_1 \tag{1}
\frac{\partial u_2}{\partial t} = D_2\nabla^2 u_2 - Ku_1u_2 + f_2 \tag{2}
\frac{\partial u_3}{\partial t} = D_3\nabla^2 u_3 +Ku_1u_2 + f_2 - Ku_3.

Now, I know that I can solve the associated linear system in each time step by creating the variational form F(u;v) = 0 and calling ‘solve(F == 0, u)’. However, to later add a point source to (1) and (2), I would like to explicitly create the matrix and right-hand side vector A and b, respectively, and apply a point source to the latter. However, when I try this, I get an error. Can I get around this error or do I have to fundamentally change the way I’m thinking about the problem?

Find attached a MWE -

from fenics import *
T = 5.0
# final time
num_steps = 500
# number of time steps
dt = T / num_steps # time step size
eps = 0.01
# diffusion coefficient
K = 10.0
# reaction rate

# Create mesh
mesh = UnitSquareMesh(30,30)

# Define function space for system of concentrations
P1 = FiniteElement('P', triangle, 1)
element = MixedElement([P1, P1, P1])
V = FunctionSpace(mesh, element)
# Define test functions
v_1, v_2, v_3 = TestFunctions(V)
# Define functions for velocity and concentrations

u = Function(V)
u_n = Function(V)
# Split system functions to access components
u_1, u_2, u_3 = split(u)
u_n1, u_n2, u_n3 = split(u_n)
# Define source terms
f_1 = Expression('pow(x[0]-0.1,2)+pow(x[1]-0.1,2)<0.05*0.05 ? 0.1 : 0',
degree=1)
f_2 = Expression('pow(x[0]-0.1,2)+pow(x[1]-0.3,2)<0.05*0.05 ? 0.1 : 0',
degree=1)
f_3 = Constant(0)
# Define expressions used in variational forms
k = Constant(dt)
K = Constant(K)
eps = Constant(eps)
# Define variational problem
F = ((u_1 - u_n1) / k)*v_1*dx \
+ eps*dot(grad(u_1), grad(v_1))*dx + K*u_1*u_2*v_1*dx \
+ ((u_2 - u_n2) / k)*v_2*dx \
+ eps*dot(grad(u_2), grad(v_2))*dx + K*u_1*u_2*v_2*dx \
+ ((u_3 - u_n3) / k)*v_3*dx \
+ eps*dot(grad(u_3), grad(v_3))*dx - K*u_1*u_2*v_3*dx + K*u_3*v_3*dx \
- f_1*v_1*dx - f_2*v_2*dx - f_3*v_3*dx

a, L = lhs(F), rhs(F)

# Time-stepping
t = 0
for n in range(num_steps):
    # Update current time
    t += dt
    
    A = assemble(a)
    b = assemble(L)
    solve(A, u.vector(), b)
    
    # Update previous solution
    u_n.assign(u)
    # Update progress bar
2

Since your problem is non-linear, this is not the associated form you solve for your problem. See for instance: Custom Newton solvers — FEniCSx tutorial which explains how to solve non-linear problems.

3

Is there an analogous tutorial for FEniCS 2019?

4

E.g: https://fenicsproject.org/olddocs/dolfin/2019.1.0/python/demos/nonlinear-poisson/demo_nonlinear-poisson.py.html

And Set krylov linear solver paramters in newton solver - #3 by nate