Issue with bilinear form when implementing manually a Newton method

[Matthieu DIAZ]

Hello,
I am solving a complex phase change problem using FEniCS. I could do this using the NonLinearSolver of FEniCS but I need to implement the Newton method manually in order to solve a far complex finite element resolution (using unfitted elements). I am forced to defined a lhs and a rhs to my problem.

This is my code :

from __future__ import print_function
from fenics import *
from dolfin import *

# TIME DISCRETIZATION

t0 = 0.0
tf = 0.5
dt = 0.01
num_steps = int((tf-t0)/dt) # number of time steps
t = t0

# MESH AND RELATED PROPERTIES

N = 100 
mesh = RectangleMesh(Point(0,0),Point(1,1),N,N)
norm = FacetNormal(mesh)
h = CellDiameter(mesh)


# BOUNDARIES

# No slip walls
no_slip = 'near(x[0], 0.0) || near(x[1], 0.0) || near(x[0], 1.0)'
# Uniform ambient temperature on bottom, left and right walls
solid_walls = 'near(x[0], 0.0) || near(x[1], 0.0) || near(x[0], 1.0)'
# Heat flux on top wall
class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1.0)


# MIXED FINITE ELEMENTS AND FUNCTION SPACE

V = VectorElement('P', mesh.ufl_cell(), 1) # velocity
Q = FiniteElement('P', mesh.ufl_cell(), 1) # pressure or temperature
element = MixedElement(V, Q, Q) # velocity, pressure, temperature in this order
W = FunctionSpace(mesh, element)


# PARAMETERS OF THE PDEs

# Dimensionless numbers
Pr = Constant(1.0)  # Prandtl 
Ste = Constant(1.0) # Stefan 
Ra = Constant(10000.0)  # Rayleigh

# Define symmetric gradient
def epsilon(u):
    return sym(nabla_grad(u))

# Define the viscosity and the phase change function
def phi(temp):
    Tr = 0.0
    r = 0.1
    return(0.5+0.5*tanh((Tr-temp)/r))

def mu(temp):
    muL = 1
    muS = 100000
    return(muL + (muS - muL)*phi(temp))


# INITIAL CONDITIONS

w_n = interpolate(Expression(("0.", "0.", "0.", "-1.0"), element=element), W)
u_n, p_n, theta_n = split(w_n)


# DIRICHLET BOUNDARY CONDITIONS

u_D = Constant((0.0 , 0.0))
theta_D = Constant(-1.0) # Equatl to the initial temperature
bc_velocity = DirichletBC(W.sub(0), u_D, no_slip)
bc_temperature = DirichletBC(W.sub(2), theta_D, solid_walls)
bcs = [bc_velocity, bc_temperature]


# NEUMANN BOUNDARY CONDITION 

# Initialize sub-domain instance
top = Top() 
m = MeshFunction('size_t', mesh, mesh.geometry().dim()-1)
m.set_all(0)
top.mark(m, 1)
# Measure
ds = Measure('ds', domain=mesh, subdomain_data=m) 
# Local application of the laser beam on the top wall
width = Constant(0.30)
flux = Expression(("0.","(x[0]>(0.5-e/2) && x[0]<(0.5+e/2)) ? 1000*(x[0]-(0.5-e/2))*(x[0]-(0.5+e/2)) : 0."), degree=2, e=width)


# WEAK FORMULATION

# Test functions
w = TestFunction(W)
(v, q, theta) = split(w)
 
# Current solution
w_ = TrialFunction(W)
(u_, p_, theta_) = split(w_)

# Define the buoyancy coupling term
fB = Ra/Pr*theta_*Constant((0.0,1.0))

# Navier Stokes momentum equation
Eq1 = dot(u_-u_n, v)/dt*dx\
    + dot(dot(grad(u_), u_), v)*dx \
    + 2*mu(theta_)*inner(epsilon(u_), epsilon(v))*dx \
    - div(v)*p_*dx \
    - dot(fB, v)*dx

# Mass equation
Eq2 = q*div(u_)*dx 

# Enthalpy equation
Eq3 = (theta_-theta_n)*theta/dt*dx \
    - 1/Ste*(phi(theta_)-phi(theta_n))/dt*theta*dx \
    + 1/Pr*dot(nabla_grad(theta_), nabla_grad(theta))*dx \
    - dot(theta_*u_, nabla_grad(theta))*dx\
    + dot(flux,norm)*theta*ds(1)

# Pressure stabilization
gamma = Constant(1E-7)
Stab_pressure = gamma*p_*q*dx

# CIP Stabilization
chi = Constant(0.003)
jump_grad_p = dot(grad(p_)('+'),norm('+'))-dot(grad(p_)('-'),norm('-'))
jump_grad_q = dot(grad(q)('+'),norm('+'))-dot(grad(q)('-'),norm('-'))
StabCIP = chi*avg(1/mu(theta_))*avg(h**3)*dot(jump_grad_p,jump_grad_q)*dS 

# Monolithic System
F = Eq1 + Eq2 + Eq3 + Stab_pressure + StabCIP

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

A = assemble(a)
b = assemble(L)

[bc.apply(A, b) for bc in bcs]

# SAVE SOLUTION TO VTK FILES

vtkfile1 = File('TestN/Velocity/velocity.pvd')
vtkfile2 = File('TestN/Pressure/pressure.pvd')
vtkfile3 = File('TestN/Temperature/temperature.pvd')

# Time-progress bar
progress_bar = Progress('Time-stepping')
set_log_level(PROGRESS)


# RESOLUTION FOR ALL TIME STEPS 
w_ = Function(W)

for n in range(num_steps):

    # Update current time
    t += dt
     
    # RESOLUTION OF THE VARIATIONAL PROBLEM

    # Newton method
    tol        = 1E-7             # absolute tolerance
    lmbda      = 1.0              # relaxation parameter    
    w_inc      = Function(W)      # residual
    eps        = 1                # L2 norm of the residual
    nIter      = 0                # number of iterations
    maxIter    = 25               # max iterations
   
    while eps > tol and nIter < maxIter:
        nIter += 1
        solve(A, w_inc.vector(), b)
        eps = norm(w_inc, 'L2')

        w_.vector()[:] += lmbda*w_inc.vector()
    
    u_, p_, theta_ = w_.split()
    w_n.vector()[:] = w_.vector()
    
    # Update the progress bar
    progress_bar.update(t/tf)

And I get the following error message, I don't know what to do to manage it. It is certainly linked to the definition of a bilinear form :

Traceback (most recent call last):
  File "LaserBeam_ManualNewton.py", line 136, in <module>
    a = lhs(F)
  File "/usr/local/lib/python3.5/dist-packages/ufl/formoperators.py", line 82, in lhs
    return compute_form_lhs(form)
  File "/usr/local/lib/python3.5/dist-packages/ufl/algorithms/formtransformations.py", line 375, in compute_form_lhs
    return compute_form_with_arity(form, 2)
  File "/usr/local/lib/python3.5/dist-packages/ufl/algorithms/formtransformations.py", line 341, in compute_form_with_arity
    return map_integrands(_transform, form)
  File "/usr/local/lib/python3.5/dist-packages/ufl/algorithms/map_integrands.py", line 39, in map_integrands
    for itg in form.integrals()]
  File "/usr/local/lib/python3.5/dist-packages/ufl/algorithms/map_integrands.py", line 39, in <listcomp>
    for itg in form.integrals()]
  File "/usr/local/lib/python3.5/dist-packages/ufl/algorithms/map_integrands.py", line 46, in map_integrands
    return itg.reconstruct(function(itg.integrand()))
  File "/usr/local/lib/python3.5/dist-packages/ufl/algorithms/formtransformations.py", line 337, in _transform
    e, provides = pe.visit(e)
  File "/usr/local/lib/python3.5/dist-packages/ufl/algorithms/transformer.py", line 104, in visit
    r = h(o, *[self.visit(op) for op in o.ufl_operands])
  File "/usr/local/lib/python3.5/dist-packages/ufl/algorithms/transformer.py", line 104, in <listcomp>
    r = h(o, *[self.visit(op) for op in o.ufl_operands])
  File "/usr/local/lib/python3.5/dist-packages/ufl/algorithms/transformer.py", line 104, in visit
    r = h(o, *[self.visit(op) for op in o.ufl_operands])
  File "/usr/local/lib/python3.5/dist-packages/ufl/algorithms/transformer.py", line 104, in <listcomp>
    r = h(o, *[self.visit(op) for op in o.ufl_operands])
  File "/usr/local/lib/python3.5/dist-packages/ufl/algorithms/transformer.py", line 108, in visit
    r = h(o)
  File "/usr/local/lib/python3.5/dist-packages/ufl/algorithms/formtransformations.py", line 146, in sum
    part, term_provides = self.visit(term)
  File "/usr/local/lib/python3.5/dist-packages/ufl/algorithms/transformer.py", line 104, in visit
    r = h(o, *[self.visit(op) for op in o.ufl_operands])
  File "/usr/local/lib/python3.5/dist-packages/ufl/algorithms/transformer.py", line 104, in <listcomp>
    r = h(o, *[self.visit(op) for op in o.ufl_operands])
  File "/usr/local/lib/python3.5/dist-packages/ufl/algorithms/transformer.py", line 108, in visit
    r = h(o)
  File "/usr/local/lib/python3.5/dist-packages/ufl/algorithms/formtransformations.py", line 146, in sum
    part, term_provides = self.visit(term)
  File "/usr/local/lib/python3.5/dist-packages/ufl/algorithms/transformer.py", line 104, in visit
    r = h(o, *[self.visit(op) for op in o.ufl_operands])
  File "/usr/local/lib/python3.5/dist-packages/ufl/algorithms/transformer.py", line 104, in <listcomp>
    r = h(o, *[self.visit(op) for op in o.ufl_operands])
  File "/usr/local/lib/python3.5/dist-packages/ufl/algorithms/transformer.py", line 108, in visit
    r = h(o)
  File "/usr/local/lib/python3.5/dist-packages/ufl/algorithms/formtransformations.py", line 74, in expr
    error("Found Argument in %s, this is an invalid expression." % ufl_err_str(x))
  File "/usr/local/lib/python3.5/dist-packages/ufl/log.py", line 172, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Found Argument in <Tanh id=140199958720136>, this is an invalid expression.

Hope someone will be able to help me. Thank you in advance.

[Jan Blechta]

Your form depends on the trial function in a non-linear way - there is
the term phi(theta_). That's not possible. Trial function is a second
argument of a bilinear form, hence the form has to be linear in it.

If you want to really solve non-linear problem, then reformulate it
without trial function, see documented demos or the tutorial.

Jan

[Matthieu DIAZ]

I changed my formulation in order to have this non linear function only depending on theta_n, so that the problem remains linear in theta_. However, I have the following issue :

Traceback (most recent call last):

  File "LaserBeam_ManualNewton.py", line 138, in <module>

  File "/usr/local/lib/python3.5/dist-packages/ufl/algorithms/formtransformations.py", line 222, in division
    error("Found Argument in denominator of %s , this is an invalid expression." % ufl_err_str(x))

  File "/usr/local/lib/python3.5/dist-packages/ufl/log.py", line 172, in error
    raise self._exception_type(self._format_raw(*message))

ufl.log.UFLException: Found Argument in denominator of <Division id=140253151989040> , this is an invalid expression.
Aborted (core dumped)

[Jan Blechta]

Looks like the same problem, error message is pretty
self-explanatory "argument in denominator", maybe the term

1/mu(theta_)

?

Jan

[Matthieu DIAZ]

1 put 1/mu(theta_n) instead. However, my code is returning another error. I put my new system of equations (only thing I changed) and the error I get :

# WEAK FORMULATION

# Test functions
w = TestFunction(W)
(v, q, theta) = split(w)
 
# Current solution
w_ = TrialFunction(W)
(u_, p_, theta_) = split(w_)

# Define the buoyancy coupling term
fB = Ra/Pr*theta_*Constant((0.0,1.0))

# Navier Stokes momentum equation
Eq1 = dot(u_-u_n, v)/dt*dx\
    + dot(dot(grad(u_), u_), v)*dx \

    + 2*mu(theta_n)*inner(epsilon(u_), epsilon(v))*dx \

    - div(v)*p_*dx \
    - dot(fB, v)*dx

# Mass equation
Eq2 = q*div(u_)*dx

# Enthalpy equation

DeltaT = Constant(0.001)

Eq3 = (theta_-theta_n)*theta/dt*dx \
    - 1/Ste*(phi(theta_n+DeltaT)-phi(theta_n))/dt*theta*dx \

    + 1/Pr*dot(nabla_grad(theta_), nabla_grad(theta))*dx \
    - dot(theta_*u_, nabla_grad(theta))*dx\
    + dot(flux,norm)*theta*ds(1)

# Pressure stabilization
gamma = Constant(1E-7)
Stab_pressure = gamma*p_*q*dx

# CIP Stabilization
chi = Constant(0.003)
jump_grad_p = dot(grad(p_)('+'),norm('+'))-dot(grad(p_)('-'),norm('-'))
jump_grad_q = dot(grad(q)('+'),norm('+'))-dot(grad(q)('-'),norm('-'))

StabCIP = chi*avg(1/mu(theta_n))*avg(h**3)*dot(jump_grad_p,jump_grad_q)*dS

# Monolithic System
F = Eq1 + Eq2 + Eq3 + Stab_pressure + StabCIP

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

A = assemble(a)
b = assemble(L)

[bc.apply(A, b) for bc in bcs]

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to perform just-in-time compilation of form.
*** Reason:  ffc.jit failed with message:

Traceback (most recent call last):

  File "/usr/lib/python3/dist-packages/dolfin/compilemodules/jit.py", line 142, in jit
    result = ffc.jit(ufl_object, parameters=p)
  File "/usr/local/lib/python3.5/dist-packages/ffc/jitcompiler.py", line 218, in jit
    module = jit_build(ufl_object, module_name, parameters)
  File "/usr/local/lib/python3.5/dist-packages/ffc/jitcompiler.py", line 134, in jit_build
    generate=jit_generate)
  File "/usr/local/lib/python3.5/dist-packages/dijitso/jit.py", line 167, in jit
    header, source, dependencies = generate(jitable, name, signature, params["generator"])
  File "/usr/local/lib/python3.5/dist-packages/ffc/jitcompiler.py", line 67, in jit_generate
    prefix=module_name, parameters=parameters, jit=True)
  File "/usr/local/lib/python3.5/dist-packages/ffc/compiler.py", line 143, in compile_form
    prefix, parameters, jit)
  File "/usr/local/lib/python3.5/dist-packages/ffc/compiler.py", line 185, in compile_ufl_objects
    analysis = analyze_ufl_objects(ufl_objects, kind, parameters)
  File "/usr/local/lib/python3.5/dist-packages/ffc/analysis.py", line 94, in analyze_ufl_objects
    for form in forms)
  File "/usr/local/lib/python3.5/dist-packages/ffc/analysis.py", line 94, in <genexpr>
    for form in forms)
  File "/usr/local/lib/python3.5/dist-packages/ffc/analysis.py", line 178, in _analyze_form
    do_apply_restrictions=True)
  File "/usr/local/lib/python3.5/dist-packages/ufl/algorithms/compute_form_data.py", line 388, in compute_form_data
    check_form_arity(preprocessed_form, self.original_form.arguments())  # Currently testing how fast this is
  File "/usr/local/lib/python3.5/dist-packages/ufl/algorithms/check_arities.py", line 152, in check_form_arity
    check_integrand_arity(itg.integrand(), arguments)
  File "/usr/local/lib/python3.5/dist-packages/ufl/algorithms/check_arities.py", line 145, in check_integrand_arity
    args = map_expr_dag(rules, expr, compress=False)
  File "/usr/local/lib/python3.5/dist-packages/ufl/corealg/map_dag.py", line 37, in map_expr_dag
    result, = map_expr_dags(function, [expression], compress=compress)
  File "/usr/local/lib/python3.5/dist-packages/ufl/corealg/map_dag.py", line 86, in map_expr_dags
    r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
  File "/usr/local/lib/python3.5/dist-packages/ufl/algorithms/check_arities.py", line 57, in product
    raise ArityMismatch("Multiplying expressions with overlapping form argument number {0}, argument is {1}.".format(x.number(), x))
ufl.algorithms.check_arities.ArityMismatch: Multiplying expressions with overlapping form argument number 1, argument is v_1.
.
*** Where:   This error was encountered inside jit.py.
*** Process: 0
***
*** DOLFIN version: 2017.2.0
*** Git changeset:  4c59bbdb45b95db2f07f4e3fd8985c098615527f
*** -------------------------------------------------------------------------

Aborted (core dumped)

[Jan Blechta]

Look, it's all the time the same problem: the form has to be linear in
its arguments (test function, trial function). Here you have: theta_*u_.
That's a non-linearity.

Jan

> *** fenics-...@googlegroups.com

> ***
> *** Remember to include the error message listed below and, if
> possible, *** include a *minimal* running example to reproduce the
> error. ***
> ***
> -------------------------------------------------------------------------
> *** Error: Unable to perform just-in-time compilation of form.
> *** Reason: ffc.jit failed with message:
> Traceback (most recent call last):
> File "/usr/lib/python3/dist-packages/dolfin/compilemodules/jit.py",
> line 142, in jit
> result = ffc.jit(ufl_object, parameters=p)
> File "/usr/local/lib/python3.5/dist-packages/ffc/jitcompiler.py",

> line 218 , in jit

> module = jit_build(ufl_object, module_name, parameters)
> File "/usr/local/lib/python3.5/dist-packages/ffc/jitcompiler.py",

> line 134 , in jit_build