I wrote a simple code below to illustrate my problem. I want to solve this PDE, but with large coefficients a and b. In the example below, a = 10 and b = 10, and the solution looks symmetric, as it should be. If you change the coefficients to a = 100 and b = 100, already the solution is not symmetric. If you change to a = 1000 and b = 1000, then the solution fails to converge.
My question is, do you know which solver and preconditionner I should use to solve this problem for large a and b? Thanks.
from dolfin import *
import pdb
mesh = RectangleMesh(Point(0.0, 0.0),Point(1.0,1.0), 50,50,'crossed')
T = FunctionSpace(mesh, "CG", 1)
class Omega(SubDomain):
def inside(self, x, on_boundary):
return between(x[0], (.0, 1.0)) and between(x[1], (.0, 1.0)) and (x[1]-0.5)**2 < 0.03
class Winlet2D(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0) and (x[1]-0.5)**2 < 0.03 and on_boundary
domains = CellFunction("size_t", mesh)
omega = Omega()
domains.set_all(0)
omega.mark(domains, 1)
dx = Measure('dx')[domains]
# Boundary conditions
boundaries = FacetFunction("size_t", mesh)
winlet = Winlet2D()
winlet.mark(boundaries, 1)
inlet_temp = Constant(280.0)
bc = DirichletBC(T, inlet_temp, boundaries, 1)
a = Constant('100.0')
b = Constant('100.0')
ka = Constant('10.0')
kw = Constant('1.0')
f = Expression('10.0')
u = project(Expression('x[0]+x[1]'),T)
t = Function(T)
v = TestFunction(T)
dT = TrialFunction(T)
# Define variational form
F = (inner( ka*grad(t), grad(v) )*dx(0)\
+ inner( kw*grad(t), grad(v) )*dx(1)\
+ inner( a*dot(grad(u),grad(t)), v )*dx(0)\
+ inner( b*dot(grad(u),grad(t)), v )*dx(1)\
- inner( f, v )*dx(0)\
- inner( f, v )*dx(1))
# Solve problem
J = derivative(F,t,dT)
problem = NonlinearVariationalProblem(F,t,bc,J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['newton_solver']['linear_solver'] = 'gmres'
prm['newton_solver']['preconditioner'] = 'hypre_amg'
prm['newton_solver']['relative_tolerance'] = 1e-12
solver.solve()
filet = File("test.pvd")
filet << t
