Solving adaptive nonlinear equation in Vector Space

variational formulation
1

Hello all,

I am trying to use adaptive meshing to solve a nonlinear equation. In my code, I am trying to minimize elastic energy (such as this example on documentation: https://fenicsproject.org/docs/dolfin/2019.1.0/python/demos/hyperelasticity/demo_hyperelasticity.py.html).

My problem is that when I use “solve(F == 0, u, bcs=bc, tol=tol, M=M)” to solve my problem, it gives me the following error:

File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py”, line 228, in solve
_solve_varproblem_adaptive(*args, **kwargs)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py”, line 321, in _solve_varproblem_adaptive
eq, u, bcs, J, tol, M, form_compiler_parameters,
ValueError: too many values to unpack (expected 8)

Please also note that since I am trying to solve a VectorFunctionSpace parameter, I have defined goal function M as below:
M = y[1] * dx(0)

where y is the deformation vector.
For solving nonlinear equation adaptively, I am using section 4 of the following document:
https://fenicsproject.org/docs/dolfin/2017.2.0/python/programmers-reference/fem/solving/solve.html

Any help is much appreciated!

2

Please supply a minimal working code example. Note that your latter link references the FEniCS 2017.2.0 API, which is now outdated.

3

Thank you very much for your reply.
Here is my code:

import numpy as np
from dolfin import *


# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"


# Create mesh and define function space
mesh = RectangleMesh(Point(-1,-1),Point(1,1), 10, 10)
V = VectorFunctionSpace(mesh, "Lagrange", 1)

dy = TrialFunction(V)            
v  = TestFunction(V)             # Test function
y  = Function(V)
u  = Function(V)                 


# Define Dirichlet boundary 

F11=1.0
F12=0.0+0.001
F21=0.0
F22=1.0-0.001

y_BC = Expression(("F11*x[0]+F12*x[1]", "F21*x[0]+F22*x[1]"),F11=F11,F12=F12,F21=F21,F22=F22,degree=2)

def x0_boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, y_BC, x0_boundary)
y=project(y_BC,V)
y_reference = Expression(("x[0]","x[1]"),degree=3)


# Elasticity model: 

E, nu = 10.0, 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))
def W(y):    
    F = grad(y)
    J  = det(F) 
    C = F.T*F  
    Ic = tr(C)

    return (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2



# Solving adaptively


M = inner(y,y)*dx()
tol = 1.e-3

W1 = W(y) * dx
dW1 = derivative(W1, y, v) # first variation (derivative at u in the direction of v)
ddW1 = derivative(dW1, y, dy) # Jacobian of dPi


solve(dW1 == 0, y, bcs=bc, tol=tol, M=M)
4

Unfortunately, I cannot reproduce your error using Docker:

 docker run -ti -v $(pwd):/home/fenics/shared -w /home/fenics/shared/ --rm quay.io/fenicsproject/stable:latest

You need to supply more information about how you have installed dolfin, and which version you are using for me to be of any more help.

5

Thanks for the follow-up.
I have installed Fenics using virtual Ubuntu on Windows 10, and this the version of dolfin I have installed:
2019.2.0.dev0

6

I can confirm that there is an issue with 2019.2.0.dev0. I have submitted a fix for it at BitBucket. So you could either add this to your source code, or go back on release to 2019.1.0, as the code works there.

7

I really appreciate your help.
Can you please tell me how I should add this fix to my source code?

8

Simply open: /usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py and change line 321/322 as indicated in my Pull request.

9

Thank you so much! The problem of my code is now resolved!

10

Hello all,

I am trying to use the “SNES” library to solve my nonlinear variational problem adaptively. However, I am getting the following error:

RuntimeError: Parameter nonlinear_variational_solver not found in Parameters object

My minimal code is as below:

import numpy as np
from dolfin import *


# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"


# Create mesh and define function space
mesh = RectangleMesh(Point(-1,-1),Point(1,1), 40, 40)
V = VectorFunctionSpace(mesh, "Lagrange", 1)

dy = TrialFunction(V)            
v  = TestFunction(V)             # Test function
y  = Function(V)
u  = Function(V)                 


# Define Dirichlet boundary 

F11=1.0
F12=0.0+0.001
F21=0.0
F22=1.0-0.001

y_BC = Expression(("F11*x[0]+F12*x[1]", "F21*x[0]+F22*x[1]"),F11=F11,F12=F12,F21=F21,F22=F22,degree=2)

def x0_boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, y_BC, x0_boundary)
y=project(y_BC,V)
y_reference = Expression(("x[0]","x[1]"),degree=3)


# Elasticity model: 

E, nu = 10.0, 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))
def W(y):    
    F = grad(y)
    J  = det(F) 
    C = F.T*F  
    Ic = tr(C)

    return (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2



# Solving adaptively

M = inner(y,y)*dx()
tol = 1.e-8

W1 = W(y) * dx
dW1 = derivative(W1, y, v) # first variation (derivative at u in the direction of v)
ddW1 = derivative(dW1, y, dy) # Jacobian of dPi

problem=NonlinearVariationalProblem(dW1,y,bcs=bc,J=ddW1)
solver = AdaptiveNonlinearVariationalSolver(problem, M)
solver.parameters['nonlinear_variational_solver'] = 'snes'
solver.solve(tol)

As in my previous problem, I am using Fenics installed on virtual Ubuntu on Windows 10, and the version of dolfin is 2019.2.0.dev0

I would be so thankful if you could help with fixing this problem.

11

Try:

solver.parameters["nonlinear_variational_solver"]["nonlinear_solver"]="snes"

A neat thing you can do is to add an embed-session:

from IPython import embed
embed()

and check what the variables of solver.parameters are, by
calling solver.parameters.keys(), solver.parameters["nonlinear_variational_solver"].keys() and so on and so forth to see what they contain

12

Hello,
I wanted to ask a follow-up question on using SNES conjugate gradient solver. Following your suggestion on adding an embed-session, I get these keys for solver.parameters["snes_solver"].keys() :

['absolute_tolerance',
 'error_on_nonconvergence',
 'krylov_solver',
 'line_search',
 'linear_solver',
 'lu_solver',
 'maximum_iterations',
 'maximum_residual_evaluations',
 'method',
 'preconditioner',
 'relative_tolerance',
 'report',
 'sign',
 'solution_tolerance']

For my problem, SNES solver converges for the lower number of meshes, however, it diverges for the higher number of meshes because of divergence reason DIVERGED_DTOL
So, I want to increase the divergence tolerance, based on this link ( SNESSetDivergenceTolerance ), but there is no key parameter for setting divergence tolerance based on the embed-session. Can you please help me with setting this parameter?

Also, for using the conjugate gradient method, what is the advantage of using fmin_cg ( scipy.optimize.fmin_cg — SciPy v1.6.1 Reference Guide ) over SNES solver?

Thank you very much!