I am trying to solve the following two scalar linear equations as a test case: Find u, \phi \in H_0^1(\Omega) \times H_0^1(\Omega)
(\nabla u, \nabla w_1) = (z + f, w_1)
(\nabla \phi, \nabla w_2) = (u - ud, w_2),
where f, u_d are given data.
Thus, my newton system to be solved looks like,
(\nabla \delta u, \nabla w_1) = - (\nabla u_0, \nabla w_1) + (z + f, w_1), \
-(\delta u, w_2) + (\nabla \delta\phi, \nabla w_2) = -(\nabla \phi_0, \nabla w_2) + (u_0, w_2) - (u_d, w_2),
where u = u_0 + \delta u and \phi = \phi_0 + \delta \phi.
Following is the code that I wrote to solve this system simultaneously as a Newton method:
from dolfin import *
import matplotlib.pyplot as plt
# Data
uex = Expression('sin(pi*x[0])*sin(pi*x[1])', degree = 2)
phiex = Expression('x[0]*(x[0]-1)*x[1]*(x[1]-1)', degree = 2)
za = -0.5; zb = 0.5; lam = 1
z = Expression('max(-0.5, min(0.5, (-1/1)*x[0]*(x[0]-1)*x[1]*(x[1]-1)))', degree = 2)
f = Expression('2*pow(pi, 2)*sin(pi*x[0])*sin(pi*x[1]) - max(-0.5, min(0.5, (-1/1)*x[0]*(x[0]-1)*x[1]*(x[1]-1)))', degree = 2)
ud = Expression('sin(pi*x[0])*sin(pi*x[1]) + 2*(x[0]-1)*x[0] + 2*(x[1]-1)*x[1]', degree = 2)
uD = Constant(0); tol = 1e-5; R = 1
# mesh
mesh = UnitSquareMesh(8,8)
## Solver
# Mixed function space
P1 = FiniteElement('CG', triangle, 1)
elem = MixedElement([P1, P1])
W = FunctionSpace(mesh, elem)
# Trial and test functions
w_1, w_2 = TestFunctions(W)
u, phi = TrialFunctions(W)
# Initial guesses (u0 and phi0)
uInt = interpolate(Expression('x[0]', degree=2), W.sub(0).collapse())
ass0 = FunctionAssigner(W.sub(0), uInt.function_space())
u0 = Function(W)
ass0.assign(u0.sub(0), uInt)
phiInt = interpolate(Expression('x[1]', degree=2), W.sub(1).collapse())
ass1 = FunctionAssigner(W.sub(1), phiInt.function_space())
phi0 = Function(W)
ass1.assign(phi0.sub(1), phiInt)
# boundary conditions
def boundary(x):
return abs(x[0]) < DOLFIN_EPS or abs(x[0]-1) < DOLFIN_EPS or abs(x[1]) < DOLFIN_EPS or abs(x[1]-1) < DOLFIN_EPS
bcSt = DirichletBC(W.sub(0), uD, boundary)
bcAd = DirichletBC(W.sub(1), 0, boundary)
Eq1 = dot(grad(u), grad(w_1))*dx; Eq2 = dot(grad(phi), grad(w_2))*dx
Eq3 = z*w_1*dx; Eq4 = f*w_1*dx; Eq5 = ud*w_2*dx; Eq6 = u*w_2*dx
iter = 0; R = 1
# Newton iterations
while (R > tol):
iter += 1
J = Eq1 - Eq6 + Eq2
res = - action(Eq1, u0) + Eq3 + Eq4 \
+ action(Eq6, u0) - action(Eq2, phi0) - Eq5
delta = Function(W)
solve(J == res, delta, [bcSt, bcAd])
un, phin = delta.split()
# termination criterion
error = (delta)**2*dx
R = sqrt(abs(assemble(error)))
u0.vector()[:] += un.vector()[:]
phi0.vector()[:] += phin.vector()[:]
print('|delta| = %e, iteration %d \n' % (R, iter))
# plot
plt.figure(1)
plot(u0)
plt.figure(2)
plot(phi0)
plt.show()
# Error computation
E1 = errornorm(uex, u0, norm_type='L2', degree_rise=3)
E2 = errornorm(phiex, phi0, norm_type='L2', degree_rise=3)
The newton iteration successfully terminates in 2 iterations. However, I get the following plots for u and \phi
which do not make any sense, since u and \phi are not vector fields and nowhere in the code I have mentioned this.
Secondly, when i compute the error using error norm function, I get the following error
Cannot compute error norm. Value shapes do not match.
Please, help!
Regards,
Jai