Issue with a Newton method implemented manually


89
views
0
8 weeks ago by
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).
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_ = Function(W)
(u_, p_, theta_) = TrialFunctions(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 :

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.
Community: FEniCS Project

1 Answer


0
8 weeks ago by
This is not an issue with your manual Newton method. The issue is found in your viscosity mu which nonlinearly depends on the form argument theta_ through the hyporbolic tangent. If you use lhs and rhs and work with TrialFunctions you need to provide a form F that can be split into bilinear and linear forms.
You can now either define a nonlinear form and compute its derivative in UFL or linearize manually and implement the resulting bilinear and linear forms in a structure similar to the one you provided.
Thank you for your answer. However, the function mu(phi(theta_)) is essential for the physical meaning of my problem. I don't see how I can linearize it. What do you mean when you talk about the jacobian ? Do you have an example ?
written 8 weeks ago by Matthieu Diaz  
Sooner or later you will have to linearize (or let the derivative function handle that for you) the problem and perform some kind of iteration, e.g. Newton. A direct solution to the non-linear problem cannot be found in general.

Basically you have a non-linear form \( F(w,\delta w)=0 \). The linearization at the variational level is performed by doing a linear expansion of the functional about a known solution \(w_k\) with increment \(\Delta w\) such that \(w_{k+1}=w_k-\Delta w\) is the solution for the next increment. (The minus sign is arbitrary but that's how it's coded in the FEniCS internal Newton solver).

Basically you then approximate
\[ F(w_{k+1},\delta w)\approx F(w_k,\delta w)+\underbrace{\frac{\text d}{\text d\varepsilon}F(w_k-\varepsilon\Delta w,\delta w)\bigg|_{\varepsilon=0}}_{=J(\Delta w,\delta w)}=0 \]
where \(J\) is the Gâteaux derivative of the functional \(F\) which is bilinear in the unknown increment \(\Delta w\) and test function \(\delta w\). This is what FEniCS can do for you with the derivative function. Since \(J\) is a bilinear form and \(F(w_k, \delta w)\) is a linear form because \(w_k\) is known, you now have the basic structure for a linear system.


A good overview is given here:
http://home.simula.no/~hpl/homepage/fenics-tutorial/release-1.0-nonabla/webm/nonlinear.html

And the non-linear Poisson demo could also help:
https://fenicsproject.org/docs/dolfin/dev/python/demos/nonlinear-poisson/demo_nonlinear-poisson.py.html
written 8 weeks ago by klunkean  
I am doing iterations for each time step that calculates the increment and I add it progressively up until it becomes small enough when I am close to the exact solution.
My purpose is to use the Cutfem library to solve that type of problems. The resolution of FE problems having unfitted boundaries in FEniCS requires the establishment of a lhs and of a rhs (I am forced to process like this because of specific custom quadrature rules).
How could I define a "a" and and "L" that I can assemble after with a tanh function ?
written 8 weeks ago by Matthieu Diaz  
What is the mathematical basis for your formulation, or to put it in other words: How is your iteration supposed to converge to some solution? A Newton procedure always uses information about the direction to find a zero.

Your variational form \(F\) must contain only bilinear and linear forms in order to extract those with the commands lhs and rhs. Only those bilinear and linear forms can be assembled into a matrix and a vector.
For a nonlinear variational form a linearization must be performed. I'll just post the most important lines. Instead of TrialFunctions, just use a Function object for your unknown:
w_ = Function(W)
(u_, p_, theta_)= split(w_)​
Then just define your form F as you did and extract the derivative in direction dw:
F = Eq1 + Eq2 + Eq3 + Stab_pressure + StabCIP

dw = TrialFunction(W)
J = derivative(F, w_, dw)

solve(F==0, w_, J=J, bcs=bcs, solver_parameters={"newton_solver":{"relative_tolerance":1e-5}})​
This would use the internal Newton solver utilizing the Jacobian J.
written 8 weeks ago by klunkean  
I already succeeded in solving the problem using the internal newton solver. I am doing a first approach with classic mixed finite elements because with unfitted finite elements I cannot use the internal non linear solver.
By several approximations, I avoided the previous problem but I have another issue. Below are my code and the issue I am facing :
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 = Constant(0.0)
    r = Constant(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_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]

# 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 

for n in range(num_steps):

    # Update current time
    t += dt
    
 
    # RESOLUTION OF THE VARIATIONAL PROBLEM
    
    w_s = Function(W)

    # 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_s.vector()[:] += lmbda*w_inc.vector()
    
    u_, p_, theta_ = w_s.split()
    w_n.vector()[:] = w_s.vector()
    
    # Update the progress bar
    progress_bar.update(t/tf)
    ​

*** -------------------------------------------------------------------------
*** 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)
​
written 8 weeks ago by Matthieu Diaz  
Well now you first cheat your way around the non-linearity in your viscosity by using the (known) temperature solution from the last time step. I don't think that's the approximation you want to make.

But your form is still non-linear because of the convective term in the linear momentum (NS) equation dot(dot(grad(u_), u_), v)*dx which is why the error is raised.

The way I see it, your code does not reflect what a Newton method is supposed to do. You mix up time increments with Newton iteration increments.

I really recommend reading into how the actual Newton method at the variational level works and what the associated linear and bilinear forms are. For a first look check out the first link in my first reply. Maybe also look into a book on the topic, preferably not one written by structural mechanicians as they will drown you in unnecessary matrix algebra.
written 8 weeks ago by klunkean  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »