This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

How can we make this Newton solver compile?

0 votes

I am trying to solve a nonlinear variational problem using the Newton solver, but I keep getting the following error:

RuntimeError:

*** Error: Unable to perform just-in-time compilation of form.

*** Reason: ffc.jit failed with message:
ArityMismatch: Multiplying expressions with overlapping form argument number 0, argument is v_0.

Here is the relevant code:

from fenics import *
import numpy as np

# Parameters
c = 1.25
l = 5.0 
P_surf = 0.0 
sigma_0 = 3.0 
gamma = 0.65 
P0 = sigma_0*gamma 
Kv = 10.0 
cm = 5.0 
# Derivative of pressure at the bottom surface.
der_p_bottom = 0.0
# Derivative of displacement at the top surface.
der_w_top = cm*gamma*sigma_0

# The initial time.
t = 0.0
# The final time
T = 2.0
# The number of time steps to take.
num_steps = 40
# The step size
dt = (T-t)/num_steps
# Width of the porous medium (y-axis length).
width = 3.0
# Thickness of the porous medium (z-axis length).
thickness = 2.0

# Create mesh and define function space.
# Mesh for the 2-D case.
mesh = RectangleMesh(Point(0.0,0.0),Point(l,width),100,10)
# Use a scalar-valued element for pressure.
P1 = FiniteElement('P', triangle, 1)
# Use a vector-valued element for displacement.
P2 = VectorElement('P', triangle, 1)
# Create a mixed element to solve the system using the fully-implicit method.
element = MixedElement(P1, P2)
# Create the function space corresponding to this mixed element.
V = FunctionSpace(mesh,element)

# Define boundary conditions for pressure.
upper_boundary = 'on_boundary && near(x[0],0)'

# An expression for the derivative of P at the bottom surface.
der_p_D = Expression('near(x[0],l) ? der_p_bottom : 0.0',degree = 1, \
    der_p_bottom=der_p_bottom, l=l)

# Set the pressure at the porous medium's top to 0 (first subsystem of V).
bc_p = DirichletBC(V.sub(0),Constant(P_surf),upper_boundary)

# Define boundary conditions for displacements.
bottom_boundary = 'on_boundary && near(x[0],%s)'%(l)

der_w_D = Expression(('near(x[0],0.0) ? der_w_top : 0.0','0.0'),degree=0, \
    der_w_top=der_w_top)

# Set the displacement at the porous medium's bottom to 0 (first dimension of
# the second subsystem of the FunctionSpace V).
bc_w = DirichletBC(V.sub(1).sub(0),Constant(0.0),bottom_boundary)

# Combine the Dirichlet boundary conditions.
bc = [bc_p, bc_w]

# Create the test functions.
v = TestFunction(V)
(v_1, v_2) = split(v)

# Define the initial conditions and previous time step values.
# The previous step's solution.
u_n = Function(V)
# Split the mixed function.
(p_n, w_n) = split(u_n)
# Create and assign the initial conditions.
u_init = Constant((P0,0.0,0.0))
u_n.interpolate(u_init)
# Create the current solution.
u_ = Function(V)
(p_, w_) = split(u_)

# Define the variational problem.
# Fluid Diffusion Equation
Eq1 = dot(p_,v_1)*dx  - dot(p_n,v_1)*dx \
    +c*dt*inner(grad(p_),grad(v_1))*dx \
    -c*dt*der_p_D*v_1*ds

# Force Equilibrium Equation
Eq2 = inner(grad(w_),grad(v_2))*dx \
    -dot(der_w_D,v_2)*ds \
    +cm*dot(grad(p_),v_2)*dx

# Add the equations to make a monolithic scheme.
F = Eq1 + Eq2

# Create the progress bar.
progress = Progress('Time-stepping')
set_log_level(PROGRESS)

# Provide the Jacobian form.
J = derivative(F,u_,v)
# Create the Newton solver.
problem = NonlinearVariationalProblem(F, u_, bc, J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters.newton_solver
prm.error_on_nonconvergence = False

# Form and solve the linear system.
for n in range(num_steps):

    # Update current time.
    t += dt

    # Solve the variational problem.
    solver.solve() # for use with the NonlinearVariationalProblem structure.

    # Update the previous solution.
    u_n.assign(u_)

    # Update the progress bar.
    progress.update(t/T)

The error occurs at:

problem = NonlinearVariationalProblem(F, u_, bc, J)

Any thoughts on why this happens? I have searched through the documentation, and I saw that UFL supports automatic functional differentiation, but the code is not working.

asked Jan 14, 2017 by SM FEniCS Novice (630 points)

1 Answer

0 votes

Found the issue. The code to compute the Jacobian was the following:

J = derivative(F,u_)

Then the code worked.

answered Jan 15, 2017 by SM FEniCS Novice (630 points)
...