Hi,
Thank you so much for following up on my question @nate @bhaveshshrimali
I am using SNES to solve my nonlinear variational problem, but I am getting the following error:
*** Error: Unable to solve nonlinear variational problem.
*** Reason: Unknown nonlinear solver type.
*** Where: This error was encountered inside NonlinearVariationalSolver.cpp.
Following is my minimal code:
import numpy as np
from numpy import linalg as LA
from numpy.linalg import inv
import matplotlib.pyplot as plt
from dolfin import *
from scipy import optimize
import numdifftools as nd
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), 50, 50)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
################################################################################
# Define boundary condictions
################################################################################
F11=0.5
F12=0.0 + 0.1
F21=0.0 - 0.1
F22=2.0
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)
tol = 1E-3
def x0_boundary(x, on_boundary):
return on_boundary and near(x[0], -1, tol) or near(x[0], 1, tol) or \
near(x[1], -1, tol) or near(x[1], 1, tol)
bcs = DirichletBC(V, y_BC, x0_boundary)
################################################################################
# Define functions
################################################################################
dy = TrialFunction(V) # Incremental displacement
v = TestFunction(V) # Test function
y = Function(V)
B = Constant((0, 0)) # Body force per unit area
T = Constant((0, 0)) # Traction force on the boundary
y_reference = Expression(("x[0]","x[1]"),degree=3)
################################################################################
# Elasticity model: Stored strain energy density (compressible neo-Hookean model)
################################################################################
E, nu = 10.0, 0.3
mu, blk = Constant(E/(2*(1 + nu))), Constant(E/(3*(1 - 2*nu)))
def W(y): #Stored strain energy density
F = grad(y) # Deformation gradient
C = F.T*F # Right Cauchy-Green tensor
Ic = tr(C)
J = det(F)
return (mu/2)*(Ic - 3) - mu*ln(J) + (blk/2 - mu/3)*(J-1)**2
################################################################################
# Total potential energy and solve
################################################################################
Pi = W(y)*dx -dot(B, (y-y_reference))*dx + dot(T, (y-y_reference))*ds
F = derivative(Pi, y, v) # first variation (derivative at u in the direction of v)
dF = derivative(F, y, dy) # Jacobian of F
problem=NonlinearVariationalProblem(F, y, bcs=bcs, J=dF)
solver=NonlinearVariationalSolver(problem)
solver.parameters['nonlinear_solver'] = 'snesncg'
solver.parameters["snes_solver"]["maximum_iterations"] = 50
solver.parameters["snes_solver"]["report"] = False
solver.solve()
To solve the problem, I tried implementing what Nate has suggested in this link: Using petsc4py.PETSc.SNES directly - #11 by niewiarowski
However, I am getting another error that says: TypeError: 'DirichletBC' object is not iterable
The code that I have used based on the mentioned link is as below:
import numpy as np
from numpy import linalg as LA
from numpy.linalg import inv
import matplotlib.pyplot as plt
from dolfin import *
from scipy import optimize
import numdifftools as nd
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), 50, 50)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
################################################################################
# Define boundary condictions
################################################################################
F11=0.5
F12=0.0 + 0.1
F21=0.0 - 0.1
F22=2.0
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)
tol = 1E-3
def x0_boundary(x, on_boundary):
return on_boundary and near(x[0], -1, tol) or near(x[0], 1, tol) or \
near(x[1], -1, tol) or near(x[1], 1, tol)
bcs = DirichletBC(V, y_BC, x0_boundary)
################################################################################
# Define functions
################################################################################
dy = TrialFunction(V) # Incremental displacement
v = TestFunction(V) # Test function
y = Function(V)
B = Constant((0, 0)) # Body force per unit area
T = Constant((0, 0)) # Traction force on the boundary
y_reference = Expression(("x[0]","x[1]"),degree=3)
################################################################################
# Elasticity model: Stored strain energy density (compressible neo-Hookean model)
################################################################################
E, nu = 10.0, 0.3
mu, blk = Constant(E/(2*(1 + nu))), Constant(E/(3*(1 - 2*nu)))
def W(y): #Stored strain energy density
F = grad(y) # Deformation gradient
C = F.T*F # Right Cauchy-Green tensor
Ic = tr(C)
J = det(F)
return (mu/2)*(Ic - 3) - mu*ln(J) + (blk/2 - mu/3)*(J-1)**2
################################################################################
# Total potential energy and solve
################################################################################
Pi = W(y)*dx -dot(B, (y-y_reference))*dx + dot(T, (y-y_reference))*ds
F = derivative(Pi, y, v) # first variation (derivative at u in the direction of v)
class SNESProblem():
def __init__(self, F, y, bcs):
V = y.function_space()
dy = TrialFunction(V)
self.L = F
self.a = derivative(F, y, dy)
self.bcs = bcs
self.y = y
def F(self, snes, x, F):
x = PETScVector(x)
F = PETScVector(F)
assemble(self.L, tensor=F)
for bc in self.bcs:
bc.apply(F, x)
def J(self, snes, x, J, P):
J = PETScMatrix(J)
assemble(self.a, tensor=J)
for bc in self.bcs:
bc.apply(J)
problem = SNESProblem(F, y, bcs)
import petsc4py, sys
petsc4py.init(sys.argv)
from petsc4py import PETSc
b = PETScVector() # same as b = PETSc.Vec()
J_mat = PETScMatrix()
snes = PETSc.SNES().create(MPI.comm_world)
snes.setFunction(problem.F, b.vec())
snes.setJacobian(problem.J, J_mat.mat())
snes.solve(None, problem.y.vector().vec())
I would be so thankful if you could help me to solve my nonlinear variational problem using SNES (and preferably, its conjugate gradient method (NCG)).
Thank you very much!
Mary