Hi,
I’m solving a phase-field fracture problem with a governing equation similar to the Cahn-Hilliard equation in the FEniCS tutorial demo. Specifically, I need to solve a minimization problem with a constraint: 0<phi<1, where phi is a phase-field variable.
The code works when using the nonlinear solver “senes” with “vinewtonssls” method. But I got the following error when using the nonlinear solver “senes” with “vinewtonrsls” method:
*** Error: Unable to successfully call PETSc function 'SNESSetFromOptions'.
*** Reason: PETSc error code is: 60 (Nonconforming object sizes).
*** Where: This error was encountered inside /build/dolfin-pnfrym/dolfin-2019.2.0~git20201207.b495043/dolfin/nls/PETScSNESSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: unknown
***
The following is a minimal running code:
from dolfin import *
import numpy as np
import math
mesh = RectangleMesh(Point(0,0), Point(1, 2),20, 40,'crossed')
n_dimension = 2
Gc =0.05
E = 3.0e4; nu = 0.0;
l0 = 0.04;
lmbda = nu*E/(1-nu**2)
mu = E/2/(1+nu)
K=lmbda+2*mu/n_dimension
def epsilon(u):
return sym(grad(u))
def sigma(u,phi):
return g(phi)*(2.0*mu*epsilon(u)+lmbda*tr(epsilon(u))*Identity(len(u)))
def psi(u):
return 0.5*K*(tr(epsilon(u)))**2+\
mu*inner(dev(epsilon(u)),dev(epsilon(u)))
def q(u,phi):
return Gc/l0*Df(phi)+Dg(phi)*psi(u)
def zeta(phi):
return 2*Gc*l0*grad(phi)
def f(phi):
return phi**2
def Df(phi):
phi = variable(phi)
return diff(f(phi),phi)
def g(phi):
return (1-phi)**2
def Dg(phi):
return 2*(phi-1)
V = FunctionSpace(mesh, 'CG', 1)
W = VectorFunctionSpace(mesh, 'CG', 1)
bot = CompiledSubDomain(" near(x[1], 0.0) ")
top = CompiledSubDomain(" near(x[1], 2.0)")
corner= CompiledSubDomain(" near(x[0], 0.0) && x[1]<1")
zero_phi = CompiledSubDomain(" (near(x[1], 2.0,0.05) || near(x[1], 0.0,0.05)) ")
load_y = Expression("u_y", u_y = 0.0, degree=1) # Displacement loading
bc_u_bot = DirichletBC(W, Constant((0.0,0.0)), bot)
bc_u_load_y = DirichletBC(W.sub(1), load_y, top)
bc_u = [bc_u_bot, bc_u_load_y]
bc_phi =DirichletBC(V, Constant(0.0), zero_phi)
u = TrialFunction(W)
du, v = TrialFunction(W), TestFunction(W)
dp, w = TrialFunction(V), TestFunction(V)
pnew, pold = Function(V), Function(V)
unew, uold = Function(W), Function(W)
F_phi = q(uold,pnew)*w*dx + inner(zeta(pnew),grad(w))*dx
F_du = inner(sigma(unew,pold),epsilon(v))*dx
J_phi = derivative(F_phi,pnew,dp)
J = derivative(F_du,unew,du)
P_phi = NonlinearVariationalProblem(F_phi, pnew, bc_phi, J_phi)
P = NonlinearVariationalProblem(F_du, unew, bc_u, J)
solver_phi = NonlinearVariationalSolver(P_phi)
solver = NonlinearVariationalSolver(P)
"Initial constrain on phi"
pmin = project(Constant(0.0), V)
pmax = project(Constant(1.0), V)
snes_solver_parameters = {"nonlinear_solver": "snes",
"symmetric": True,
"snes_solver": {
"method": "vinewtonrsls", # "vinewtonrsls" or "vinewtonssls"
"maximum_iterations":20,
"absolute_tolerance":1e-10,
"relative_tolerance":1e-10,
"report": True,
"error_on_nonconvergence": False}}
solver_phi.parameters.update(snes_solver_parameters)
prm_phi = solver_phi.parameters
prm_phi['snes_solver']['line_search'] = 'basic'
u_r = 1
dt = 0.01
t = 0
t_end = 0.15/u_r
tol = 1e-3
conc_f = File ("./ResultsDir/phi.pvd",'compressed')
iter_max = 5000
# Staggered scheme
while t <= t_end :
load_y.u_y += dt*u_r
iter = 0
err = 1.0
P_phi.set_bounds(pmin,pmax)
while err > tol and iter <= iter_max:
iter += 1
solver.solve()
err_u = errornorm(unew,uold,norm_type = 'l2', degree_rise=0, mesh = None)/(norm(uold)+1e-20)
uold.assign(unew)
solver_phi.solve()
err_phi = errornorm(pnew,pold,norm_type = 'l2', degree_rise=0, mesh = None)/(norm(pold)+1e-20)
pold.assign(pnew)
if round(iter) % 500 == 0:
print(t,iter, err)
err = max(err_phi,err_u)
if err < tol or iter == iter_max:
uold.assign(unew)
pold.assign(pnew)
pmin.assign(pold)
conc_f << pnew
t += dt
Anyone can help? Thanks a lot!