Hi Bhavesh, two things (see below).
On 25 Jun 2021, at 14:12, Bhavesh Shrimali <bhavesh.shrimali@gmail.com> wrote:
Also, for definiteness, the complete code:
from firedrake import * parameters["form_compiler"]["quadrature_degree"]=5 mesh = UnitCubeMesh(5,5,5) P2 = FiniteElement("CG", "tetrahedron", 2) B4 = FiniteElement("B", "tetrahedron", 4) FB3 = FiniteElement("FacetBubble", "tetrahedron", 3) Vh = VectorFunctionSpace(mesh, P2 + B4 + FB3) Qh = FunctionSpace(mesh, "DG", 1) Wh = Vh * Qh w = Function(Wh) dw = TrialFunction(Wh) v = TestFunction(Wh) u, p = split(w) mu, lmbda = Constant(1.), Constant(1.e3) Js = (lmbda + p + sqrt((lmbda + p)**2. + 4.*lmbda*mu ))/(2.*lmbda) F = Identity(mesh.geometric_dimension()) + grad(u) J = det(F) I1 = inner(F, F) Pi = p * (J - Js) + mu/2 * (I1 - 3) - mu * ln(Js) + lmbda/2 * (Js - 1)**2 Pi *= dx a = derivative(Pi, w, v) Jac = derivative(a, w, dw) bcs = [ DirichletBC(Wh.sub(0).sub(0), Constant(0.), 1), DirichletBC(Wh.sub(0).sub(1), Constant(0.), 3), DirichletBC(Wh.sub(0).sub(2), Constant(0.), 5), DirichletBC(Wh.sub(0).sub(0), Constant(0.5), 2) ] problem = NonlinearVariationalProblem(a, w, bcs, J=Jac) solver = NonlinearVariationalSolver(problem) solver_parameters={'snes_monitor': None, 'snes_view': None, 'ksp_monitor_true_residual': None, 'snes_converged_reason': None, 'ksp_converged_reason': None, 'snes_type': 'test', 'mat_type': 'aij'}
solver.parameters.update(solver_parameters)
You need to provide this dictionary of solver parameters when you construct the NonlinearVariationalSolver, rather than updating parameters after the fact. so: solver = NLVS(problem, solver_parameters=...) Second, it looks like the problem is that the nonlinear is not converging, even though the linear solve is fine. I suspect that you need to start with better initial guesses. Can you do continuation in a parameter (such as lambda) to ramp up the nonlinearity? You may also need to add a linesearch (e.g. 'snes_linesearch_type': 'bt') If I set backtracking linesearch and do: for val in [1, 10, 100, 1000]: lmbda.assign(val) solver.solve() I can converge a problem on a 4x4x4 mesh. Lawrence