Mixed Time Dependant Equation : "No Jacobian form Specified"

Question

Hi,

I have a complex equation to solve, in order to solve it I tried to mix the examples : "A system of advection–diffusion–reaction equations ", "A nonlinear Poisson equation" and "The heat equation" because my problem use different methods detailed in these examples.

I have the following error :

No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 0.000e+00 (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton solver finished in 0 iterations and 0 linear solver iterations.

I tried lots of things that I found here but nothing change.

Here is my code :

#coding=utf8

##############################

from fenics import *

##############################

mesh = RectangleMesh(Point(0.0, 0.0), Point(15.0, 5.0), 90, 30)

T = 5.0            # final time
num_steps = 10    # number of time steps
dt = T / num_steps # time step size

# Parameters

Dn       = 1.
chi_0    = 1.
k_1      = 1.
rho_0    = 1.
omega    = 1.
mu       = 1.
lambda_0 = 1.

# 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 pressure and concentrations

#w = Function(W)
u = Function(V)
u_n = Function(V)

# Split systeme functions to access components
u_1, u_2, u_3 = split(u)
u_n1, u_n2, u_n3 = split(u_n)

# Define expressions used in variational forms

k = Constant(dt)
Dn = Constant(Dn)
chi_0 = Constant(chi_0)
k_1 = Constant(k_1)
rho_0 = Constant(rho_0)
omega = Constant(omega)
mu = Constant(mu)
lambda_0 = Constant(lambda_0)

# Define variational problem

F = ((u_1 - u_n1) / k)*v_1*dx \
  - Dn*dot(grad(u_1), grad(v_1))*dx \
  + ((chi_0*k_1*u_1) / (k_1 + u_3))*dot(grad(u_3), grad(v_1))*dx \
  + rho_0*u_1*dot(grad(u_2), grad(v_1))*dx \
  + ((u_2 - u_n2) / k)*v_2*dx \
  - omega*u_1*v_2*dx + mu*u_1*u_2*v_2*dx \
  + ((u_3 - u_n3) / k)*v_3*dx + lambda_0*u_1*u_3*v_3*dx

# Create VTK files for visualization output

vtkfile_u_1 = File('endothecial_cells_migration/u_1.pvd')
vtkfile_u_2 = File('endothecial_cells_migration/u_2.pvd')
vtkfile_u_3 = File('endothecial_cells_migration/u_3.pvd')

# Create progress bar

progress = Progress('Time-stepping')
set_log_level(PROGRESS)

# Time-stepping
t = 0
for n in range(num_steps):

    # Update current time
    t += dt

    # Solve variationnal problem for time step
    solve(F == 0, u)

    # Save solution to file (VTK)
    _u_1, _u_2, _u_3 = u.split()
    vtkfile_u_1 << (_u_1, t)
    vtkfile_u_2 << (_u_2, t)
    vtkfile_u_3 << (_u_3, t)

    # Update previous solution
    u_n.assign(u)

    # UPdate progress bar
    progress.update(t / T)

# Hold plot
interactive ()

I think that the problem come from the expression of F. Or maybe I need to change the arguments in the solver. But I didn't manage to correct the error.

Thank you for your help.

Comments

I can't see what your initial system state is. I.e. is $u(x, 0) = 0$? Have you tried perturbing the system?

---

I did not create initial value : I added this and it worked (no error) but I need to work on my initial condition (and parameters) because the results are weird

# Define initial condition
P = FunctionSpace(mesh, P1)
u_D = Expression('exp(-a*pow(x[0], 2))',
                degree=2, a=0.03)
u_n1 = interpolate(u_D, P)

Thank you

Answer 1

Whenever I encounter this error(nan residual), I change the initial condition, i.e.

 from numpy import ones
 u.vector().set_local(ones(u.vector().size()))

Comments

It actually worked but considering this initial condition the system doesn't evolve because there is no source term for u_1. I have to change the initial condition in order to let the system evolve.

Thank you.

---

Maybe one has to add

 u.vector().apply("")

---

I don't understand where i should add this.

I just add initial condition but it return me an error : Its tell me taht the Newton solver did not converge.

# Define initial condition
u_D1 = Expression('pow(x[0]-0.1,2)+pow(x[1]-5,2)<3 ? 1 : 0',
                 degree=1)
u_n1 = interpolate(u_D1, P)
u_D2 = Expression('(-atan(x[0]/3-2.5)/2.6)+1',
                 degree=1)
u_n2 = interpolate(u_D2, P)
u_D3 = Expression('(atan(x[0]/3-2.5)/2.6)+1',
                 degree=1)
u_n3 = interpolate(u_D3, P)

---

Your problem is that you want to set a specific non-zero initial condition for a mixed problem. You need to use the FunctionAssigner for this

checkout
https://fenicsproject.org/qa/10496/how-to-make-deep-copies-with-split-u?show=10507#a10507
https://fenicsproject.org/qa/11331/existing-solution-as-one-element-of-initial-condition?show=11331#q11331