Hi,
I have a strange problem in which the weak form has a tanh(a+b*div(u)) term. I wish to solve the problem for a range of values of b. The problem converges fine for small b, but using the previous solution as an initial guess in a loop increasing the value of b, eventually the solver has convergence issues. Does anyone have any ideas on reasons that this might happen?
I have written a minimal `working’ example which has approximations to the constants I use in my problem and it fails to converge at a similar value of b:
from dolfin import *
from mshr import *
from ufl import tanh
L = 1.
l = 0.3
# Sub domain for Periodic boundary condition = edges of the box
class PeriodicBoundary(SubDomain):
L = 1.0
# Making the negative parts the reference parts and the positive the slave parts
# Don't think I need to worry about corners as these are rounded in my geometry but included anyway
def inside(self, x, on_boundary):
return bool((near(x[0], -L/2.) or near(x[1], -L/2.) or near(x[2], -L/2.)) and \
(not ((near(x[0], -L/2.) and near(x[1], L/2.)) or \
(near(x[0], -L/2.) and near(x[2], L/2.)) or \
(near(x[1], -L/2.) and near(x[0], L/2.)) or \
(near(x[1], -L/2.) and near(x[2], L/2.)) or \
(near(x[2], -L/2.) and near(x[0], L/2.)) or \
(near(x[2], -L/2.) and near(x[1], L/2.)))) and on_boundary)
# Map positive sides to negative sides
def map(self, x, y):
if (near(x[0], L/2.) and near(x[1], L/2.) and near(x[2], L/2.)):
y[0] = x[0] - L
y[1] = x[1] - L
y[2] = x[2] - L
elif (near(x[0], L/2.) and near(x[1], L/2.)):
y[0] = x[0] - L
y[1] = x[1] - L
y[2] = x[2]
elif (near(x[1], L/2.) and near(x[2], L/2.)):
y[0] = x[0]
y[1] = x[1] - L
y[2] = x[2] - L
elif (near(x[2], L/2.) and near(x[0], L/2.)):
y[0] = x[0] - L
y[1] = x[1]
y[2] = x[2] - L
elif near(x[0], L/2.):
y[0] = x[0] - L
y[1] = x[1]
y[2] = x[2]
elif near(x[1], L/2.):
y[0] = x[0]
y[1] = x[1] - L
y[2] = x[2]
elif near(x[2], L/2.):
y[0] = x[0]
y[1] = x[1]
y[2] = x[2] - L
else:
y[0] = 1000.
y[1] = 1000.
y[2] = 1000.
mesh = Mesh("working_example_mesh.xml")
Vh = VectorFunctionSpace(mesh, "Lagrange", 2, constrained_domain=PeriodicBoundary())
Vh1 = FunctionSpace(mesh, "Lagrange", 2, constrained_domain=PeriodicBoundary())
du = TrialFunction(Vh)
w = TestFunction(Vh)
u = Function(Vh)
bdry = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
class inner_bndry(SubDomain):
def inside(self, x, on_boundary):
return bool(not (near(x[0], -0.5*L) or near(x[0], 0.5*L) or near(x[1], -0.5*L) or near(x[1], 0.5*L) or near(x[2],-0.5*L) or near(x[2], 0.5*L)) and on_boundary)
class faces1(SubDomain):
def inside(self, x, on_boundary):
return bool((near(x[0], -0.5*L) or near(x[1], -0.5*L) or near(x[2], -0.5*L))) and on_boundary
class faces2(SubDomain):
def inside(self, x, on_boundary):
return bool((near(x[0], 0.5*L) or near(x[1], 0.5*L) or near(x[2], 0.5*L))) and on_boundary
ds = Measure("ds")(subdomain_data = bdry)
n = FacetNormal(mesh)
bcs = []
#------------------------
# Solver Settings
#------------------------
class Newton_Interface(NonlinearProblem):
def __init__(self, a, L):
NonlinearProblem.__init__(self)
self.L = L
self.a = a
#self.reset_sparsity = True
def F(self, b, x):
assemble(self.L, tensor=b)
def J(self, A, x):
assemble(self.a, tensor=A)#, reset_sparsity=self.reset_sparsity)
#self.reset_sparsity = False
parameters["krylov_solver"]["nonzero_initial_guess"] = True
solver = NewtonSolver()
solver.parameters["preconditioner"] = "hypre_amg"
solver.parameters["linear_solver"] = "gmres"
solver.parameters["convergence_criterion"] = "residual"
solver.parameters["relative_tolerance"] = 1e-8
#--------------------------------------
# Solving for u
#--------------------------------------
for c in range(0,10):
c1,c2,c3,c4,c5 = 1.38, -1.44 + 0.27656*c, 0.92, 3.0714, -0.6142*c
gu = grad(u)
L = c1*div(u)*div(w)*dx + inner(grad(w),grad(u)+gu.T)*dx \
-c4*tanh(c2 + c3*div(u))*div(w)*dx(mesh) \
-inner(w,Constant(c4 + c5)*n)*ds(0) \
+10.**(-6)*inner(u,w)*dx # workaround to set nullspace
a = derivative(L, u, du)
problem = Newton_Interface(a, L)
solver.solve(problem, u.vector())
out = tanh(c2+c3*div(u))
out = project(out,Vh1)
File('working_example_tanh.pvd') << out
I created the mesh using gmsh as it needs to be periodic and attempted to reproduce it using mshr but the problem did not even converge for small b in this case. I can send the .xml file to anyone wishing to run the minimal working code. The shape of the mesh is shown in the picture below.
Any comments would be very welcome!
