Hi Everyone, Recently I wrote a fenics python code for solving a nonlinear
advection diffusion equation which worked very well. However, I came across
one issue which I couldn't resolve and which is illustrated by the attached
code (simpler than original project but still seems to contain the problem).
The example problem is identical to the one described here
http://fenicsproject.org/documentation/tutorial/nonlinear.html
except that the nonlinear term q(u) is (1-u)^(-2), i.e., the PDE is
grad(grad(u)/(1-u)^2) = -f
The attached code gives different results if q is written as "(1-u)**(-2)" and
"(1/(1-u))**2". Is there a reason for this? I would have thought they should
be equivalent.
Thanks,
Dave
from dolfin import *
# domain size
domainX = 100
domainY = 100
# mesh
elementSide = 1.0
meshNX = int(domainX / elementSide)
meshNY = int(domainY / elementSide)
mesh = RectangleMesh(0, 0, domainX, domainY, meshNX, meshNY)
#plot(mesh, interactive=True)
# function space
V = FunctionSpace(mesh, 'Lagrange', 1)
source = Function(V)
# create initial conditions and interpolate
s = Expression('-(2/(5*5*pi)) * exp(-(pow((x[0]-50)/5,2) + pow((x[1]-50)/5,2)))')
source.interpolate(s)
def q(u):
# appears to work
#return (1/(1-u))**2
# appears not to work
return (1-u)**(-2)
def boundary(x, on_boundary):
return on_boundary
Gamma = DirichletBC(V, Constant(0.0), boundary)
bcs = [Gamma]
u = TrialFunction(V)
v = TestFunction(V)
w = Function(V)
F = dot(q(u)*grad(u), grad(v))*dx + v*source*dx
F = action(F, w)
J = derivative(F, w, u)
problem = NonlinearVariationalProblem(F, w, bcs, J)
solver = NonlinearVariationalSolver(problem)
solver.solve()
plot(w, interactive=True)
_______________________________________________
fenics mailing list
[email protected]
http://fenicsproject.org/mailman/listinfo/fenics
