Hi,
I would like to interface with the petsc4py.PETSc.SNES object directly to access some specific PETSc functions, but I cannot find any documentation on how to use to the SNES object in FEniCS e.g. how to convert the linear form F, Jacobian J, and the boundary conditions to PETSc objects etc. Take for example the fenics hyperelasticty demo as a starting point:
from dolfin import *
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
"eliminate_zeros": True, \
"precompute_basis_const": True, \
"precompute_ip_const": True}
# Create mesh and define function space
mesh = UnitCubeMesh(24, 16, 16)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
# Mark boundary subdomians
left = CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 1.0)
# Define Dirichlet boundary (x = 0 or x = 1)
c = Expression(("0.0", "0.0", "0.0"), degree=2)
r = Expression(("scale*0.0",
"scale*(y0 + (x[1] - y0)*cos(theta) - (x[2] - z0)*sin(theta) - x[1])",
"scale*(z0 + (x[1] - y0)*sin(theta) + (x[2] - z0)*cos(theta) - x[2])"),
scale = 0.5, y0 = 0.5, z0 = 0.5, theta = pi/3, degree=2)
bcl = DirichletBC(V, c, left)
bcr = DirichletBC(V, r, right)
bcs = [bcl, bcr]
# Define functions
du = TrialFunction(V) # Incremental displacement
v = TestFunction(V) # Test function
u = Function(V) # Displacement from previous iteration
B = Constant((0.0, -0.5, 0.0)) # Body force per unit volume
T = Constant((0.1, 0.0, 0.0)) # Traction force on the boundary
# Kinematics
d = u.geometric_dimension()
I = Identity(d) # Identity tensor
F = I + grad(u) # Deformation gradient
C = F.T*F # Right Cauchy-Green tensor
# Invariants of deformation tensors
Ic = tr(C)
J = det(F)
# Elasticity parameters
E, nu = 10.0, 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))
# Stored strain energy density (compressible neo-Hookean model)
psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2
# Total potential energy
Pi = psi*dx - dot(B, u)*dx - dot(T, u)*ds
# Compute first variation of Pi (directional derivative about u in the direction of v)
F = derivative(Pi, u, v)
# Compute Jacobian of F
J = derivative(F, u, du)
solve(F == 0, u, bcs, J=J, form_compiler_parameters=ffc_options)
Now, to replace the magic solve call with a custom SNES solver… I know how to use NonlinearProblem together with PETScSNESSolver to replace the magic solve as below, but this is NOT what I want:
class Problem(NonlinearProblem):
def __init__(self, J, F, bcs):
self.bilinear_form = J
self.linear_form = F
self.bcs = bcs
NonlinearProblem.__init__(self)
def F(self, b, x):
assemble(self.linear_form, tensor=b)
for bc in self.bcs:
bc.apply(b, x)
def J(self, A, x):
assemble(self.bilinear_form, tensor=A)
for bc in self.bcs:
bc.apply(A)
problem = Problem(J, F, bcs)
solver = PETScSNESSolver()
solver.solve(problem, u.vector())
Instead I would like to use petsc4py directly as suggested by @nate here:
Set krylov linear solver paramters in newton solver starting with:
import petsc4py, sys
petsc4py.init(sys.argv)
from petsc4py import PETSc
snes = PETSc.SNES().create(MPI.comm_world)
....
Somehow somewhere there have to be calls to snes.setFunction() and snes.setJacobian() and handling bcs etc:
snes.setFunction(???)
snes.setJacobian(???)
....
snes.solve(??)
Does anyone have any code that does something like this that they could share? It would be much appreciated! Thanks!
