Hi everybody
No iterative solver is converging for my simulations and I don't know why.
Below is a minimal working example with multiple 'solve' lines, each with another solver.
When I use the default solver, the solution is computed almost instantly, but when I use any iterative solver (gmres, cg, minres, ...) the solution never converges.
I tested it on two different computers (Debian Jessie/fenics 1.4, Ubuntu 14.04/fenics 1.5) and it didn't work on both.
Please give me a hint,
thank you in advance!
#!/usr/bin/env python
# encoding: utf-8
import math as m
import numpy as np
from dolfin import *
Lx = 0.19
Ly = 0.01
Lz = 0.01732
Nx = 380
Ny = 1
Nz = 1
mesh = BoxMesh(0., 0., 0., Lx, Ly, Lz, Nx, Ny, Nz)
c = 343.4
rho = 1.204
########################################################################
# set periodic boundary conditions and create constrained function space
class PeriodicDomain(SubDomain):
# set periodic boundary conditions
def inside(self, x, on_boundary):
# return True if on left or bottom boundary AND NOT on one of the two
# slave edges
return bool((near(x[1], 0) or near(x[2], 0.)) and
(not ((near(x[1], Ly) and near(x[2], 0.)) or
(near(x[1], 0) and near(x[2], Lz)))) and on_boundary)
def map(self, x, y):
if near(x[1], Ly) and near(x[2], Lz):
y[0] = x[0]
y[1] = x[1] - Ly
y[2] = x[2] - Lz
elif near(x[1], Ly):
y[0] = x[0]
y[1] = x[1] - Ly
y[2] = x[2]
elif near(x[2], Lz):
y[0] = x[0]
y[1] = x[1]
y[2] = x[2] - Lz
else:
y[0] = -1000
y[1] = -1000
y[2] = -1000
periodicdomain = PeriodicDomain()
# ################################################################################
class ABC_Sommerfeld(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], Lx)
abcs = ABC_Sommerfeld()
# Initialize facet function for external boundary domain
external_boundaries = FacetFunction("size_t", mesh)
external_boundaries.set_all(0)
abcs.mark(external_boundaries, 1)
################################################################################
# Define constants and expressions
frequency = 500.
omega = 2.*m.pi*frequency
k = omega/c
p0r = 1.
p0i = 1.
f1 = Expression('2*k_*p0i_', k_=k, p0i_=p0i)
f2 = Expression('2*k_*p0r_', k_=k, p0r_=p0r)
ksquared = Constant(k**2.)
k = Constant(k)
################################################################################
# create two independent function spaces
W = FunctionSpace(mesh, 'Lagrange', 2, constrained_domain=periodicdomain)
W2 = W * W
(pr, pi) = TrialFunction(W2)
(w1, w2) = TestFunction(W2)
a_air = (pr*w1*ksquared + pi*w2*ksquared -
inner(nabla_grad(pr), nabla_grad(w1)) -
inner(nabla_grad(pi), nabla_grad(w2)))*dx()
a_right = (k*pi*w1 - k*pr*w2 )*ds(1)
L = (+f1*w1 - f2*w2)*ds(1)
a = a_air + a_right
# assemble
A = assemble(a, exterior_facet_domains=external_boundaries)
b = assemble(L, exterior_facet_domains=external_boundaries)
# Compute solution
p = Function(W2)
P = p.vector()
# works
solve(A, P, b)
# # doesn't work
# solve(A, P, b, 'gmres', 'none')
# solve(A, P, b, 'gmres', 'ilu')
# # doesn't work
# solve(A, P, b, 'minres', 'none')
# # doesn't work
# solve(A, P, b, "cg", "none")
# # doesn't work (from Fenics Q&A)
# prec = PETScPreconditioner('hypre_amg')
# PETScOptions.set('pc_hypre_boomeramg_relax_type_coarse', 'jacobi')
# solver = PETScKrylovSolver('cg', prec)
# solver.parameters['absolute_tolerance'] = 1.0e-6
# solver.parameters['relative_tolerance'] = 1.0e-10
# solver.parameters['maximum_iterations'] = 1000
# solver.parameters['monitor_convergence'] = True
# A_petsc = as_backend_type(A)
# b_petsc = as_backend_type(b)
# x_petsc = as_backend_type(P)
# solver.set_operator(A_petsc)
# solver.solve(x_petsc, b_petsc)
# solver = KrylovSolver("gmres")
# solver.parameters["monitor_convergence"] = True
# solver.parameters["error_on_nonconvergence"] = False
# solver.solve(A, P, b)
pr, pi = p.split()
plot(pr)
plot(pi)
interactive()
