Hi am solving a nonlinear time dependent system of pde and am following the cahnHilliard demo. I got this error in the newton solver which I dont quiet understand. Below is my complete code
from dolfin import *
# Class representing initial conditions
class InitialConditions(Expression):
def eval(self, values, x):
values[0] = x[1]
values[1] = x[1]
values[2] = x[1]
values[3] = x[1]
values[4] = 0.01
values[5] = 0
def value_shape(self):
return (6,)
#constants
c = Constant(2/3)
xi0= Constant(1)
n = Constant(2)
alpha = Constant(25)
# Class interfacing with newton solver
class MyEquations(NonlinearProblem):
def __init__(self, a, L, bcs):
NonlinearProblem.__init__(self)
self.L = L
self.a = a
self.bcs = bcs
self.reset_sparsity = True
def F(self, b, x):
assemble(self.L, tensor=b, bcs = self.bcs)
def J(self, A, x):
assemble(self.a, tensor=A, reset_sparsity=self.reset_sparsity, bcs=self.bcs)
self.reset_sparsity = False
class PeriodicBoundary(SubDomain):
# left boundary is "target domain"
def inside(self, x, on_boundary):
return bool(x[0]< DOLFIN_EPS and x[0]> - DOLFIN_EPS and on_boundary)
# Map right boundary to left boundary
def map(self, x, y):
y[0] = x[0] - 1.0
y[1] = x[1]
pbc = PeriodicBoundary()
mesh = RectangleMesh(Point(-1,-1),Point(1,1),20,20) #Domain[-1,1]*[-1,1]
V1= VectorFunctionSpace(mesh, 'Lagrange',degree=2, constrained_domain = pbc)
V2=VectorFunctionSpace(mesh, 'Lagrange', degree=2, constrained_domain= pbc)
Q1 = FunctionSpace(mesh, 'Lagrange', degree=1, constrained_domain = pbc )
Q2 = FunctionSpace(mesh, 'Lagrange', degree=1, constrained_domain = pbc)
VQ = MixedFunctionSpace([V1,V2, Q1, Q2])
# Define trial and test functions
m,h,q,r = TestFunctions(VQ)
duu = TrialFunction(VQ)
H= Function(VQ) # Current solution
H0 = Function(VQ) # solution from previous converged steps
#dUU, du, dphi, dp = split(duu)
U, u, phi, p = split(H)
U0, u0, phi0, p0 = split(H0)
#Create initial conditions and interpolate
u_init = InitialConditions()
H.interpolate(u_init)
H0.interpolate(u_init)
# Sub domain for Dirichlet boundary condition
class Top(SubDomain):
def inside(self, x, on_boundary):
return (x[1]>1 - DOLFIN_EPS and on_boundary)
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return (x[1]< -1 + DOLFIN_EPS and on_boundary)
UB= Constant((-1,0))
UT= Constant((1,0))
ub= Constant((-1,0))
ut= Constant((1,0))
top= Top()
bott = Bottom()
bcUB= DirichletBC(VQ.sub(0), UB ,bott)
bcUT= DirichletBC(VQ.sub(0), UB ,top)
bcub= DirichletBC(VQ.sub(1), ub ,bott )
bcut= DirichletBC(VQ.sub(1), ub ,top )
bcs = [bcUB,bcUT, bcub,bcut]
k_phi = (phi/phi0)**n
eta = exp(-alpha*(phi-phi0))
xi = xi0* eta
dt=0.0001
f1 = phi*q*dx - phi0*q*dx - dt*div((1-phi)*U)*q*dx
f2 = inner(phi*(u-U), m)*dx + k_phi * p*div(m)*dx
f3 = p*div(h)*dx + eta*inner(grad(U) + grad(U).T, grad(h))*dx - (xi*div(U)) *div(h)*dx + c*(eta*div(U))*div(h)*dx
f4 = div(phi*u + (1-phi)*U)*r*dx
L = f1 + f2 + f3 + f4
a = derivative(L, H, duu )
problem = MyEquations(a, L, bcs)
solver = NewtonSolver()
solver.parameters['linear_solver']='lu'
solver.parameters['convergence_criterion'] = 'incremental'
solver.parameters['relative_tolerance'] = 1e-6
T_total = 50*dt
t=0.0
while (t<T_total):
t+=dt
H0 = H
solver.solve(problem, H)
The error is below as
solver.solve(problem, H)
TypeError: in method 'NewtonSolver_solve', argument 3 of type 'dolfin::GenericVector &'
I dont quiet understand this error. Could anyone help please. Thanks.