gmres convergence stronly depends on frequency

Question

Hello everybody

I try to solve a simple acoustic test case consisting of an air filled box with a point source inside. On one surface is a thin layer of porous absorbing material.

The problem is, that gmres does not converge for certain frequencies.
For example for 300, 290, 280, 270Hz the solver converges and the solution process takes about 10min each on my machine. For 230Hz the solution process takes about 2:40 for 220Hz about 1:30 and for 190Hz gmres does not converge at all.

This is counterintuitive because the solution oscillates less for lower frequencies, so I would expect that lower frequencies would converge faster.

Where can I begin to search for the problem?
Any hint would be very much appreciated.

This is the code of the problem:

#!/usr/bin/env python
# encoding: utf-8
# 2015-03-04 10:58

import math as m
import numpy as np
from dolfin import *

thickness = 0.1

L = 3.2
Lx = L
Ly = L
Lz = L
Nx = 32
Ny = 32
Nz = 32
mesh = BoxMesh(0., 0., 0., Lx, Ly, Lz, Nx, Ny, Nz)

W = FunctionSpace(mesh, 'Lagrange', 2)
W2 = W * W

########################################
# add subdomains
class Air(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] >= thickness

class Absorber(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] <= thickness

class AA_interface(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], thickness)

air = Air()
absorber = Absorber()
aai = AA_interface()

# Initialize mesh function for interior domains
domains = CellFunction("size_t", mesh)
domains.set_all(0)
absorber.mark(domains, 1)
air.mark(domains, 2)

# Initialize mesh function for internal boundary domains
internal_boundaries = FacetFunction("size_t", mesh)
internal_boundaries.set_all(0)
aai.mark(internal_boundaries, 1)

dx = Measure("dx")[domains]
dS = Measure("dS")[internal_boundaries]
########################################
frequency = 300
print frequency

# acoustic parameters
c = 343.4
rho = 1.204
omega = 2.*m.pi*frequency
agemo = Constant(1./omega)
k = omega/c
ksquared = Constant(k**2.)
k = Constant(k)

# Absorber Parameters
xi = 10000.
kappa = 1.3
sigma = 0.99
zr = Constant(kappa)
zi = Constant(-xi*sigma / (rho))
ks = Constant(kappa/sigma)
xr = Constant(xi/rho)

########################################
# define functions
(pr, pi) = TrialFunction(W2)
(w1, w2) = TestFunction(W2)

# air domain
a_air = (pr*w1*ksquared + pi*w2*ksquared -
         inner(nabla_grad(pr), nabla_grad(w1)) - 
         inner(nabla_grad(pi), nabla_grad(w2)))*dx(2)

# absorber domain
a_absorber = (pr*zr*w1*ksquared - zi*agemo*pi*w1*ksquared + 
              zi*agemo*pr*w2*ksquared + zr*pi*w2*ksquared -
              inner(nabla_grad(pr), nabla_grad(w1)) - 
              inner(nabla_grad(pi), nabla_grad(w2)))*dx(1)

a = a_air + a_absorber

# right hand side
L = Constant('0.')*(w1 - w2)*dx(1)

# interface
n = FacetNormal(mesh)

denominator = Constant(1./(kappa*kappa/(sigma*sigma) +
                       xi*xi/(rho*rho*omega*omega)))

# (w1*dot(grad(pr)('-'), n('-')) + w2*dot(grad(pi)('-'), n('-')))*dS(1) =
a_abs_interface = (- ks*      w1('-')*dot(grad(pr('-')), n('-'))
                   - xr*agemo*w1('-')*dot(grad(pi('-')), n('-'))
                   - ks*      w2('-')*dot(grad(pi('-')), n('-'))
                   + xr*agemo*w2('-')*dot(grad(pr('-')), n('-'))
                   )*dS(1)


# (w1*dot(grad(pr)('+'), n('+')) + w2*dot(grad(pi)('+'), n('+')))*dS(1) =
a_air_interface = (- denominator*ks*      w1('+')*dot(grad(pr('+')), n('+'))
                   + denominator*xr*agemo*w1('+')*dot(grad(pi('+')), n('+'))
                   - denominator*ks*      w2('+')*dot(grad(pi('+')), n('+'))
                   - denominator*xr*agemo*w2('+')*dot(grad(pr('+')), n('+'))
                   )*dS(1)


a_interface = a_abs_interface + a_air_interface
a += a_interface

A = assemble(a)
b = assemble(L)

ps0 = PointSource(W2.sub(0), Point(thickness + 1.5, Ly/2., Lz/2.), 1.)
ps1 = PointSource(W2.sub(1), Point(thickness + 1.5, Ly/2., Lz/2.), 1.)
ps0.apply(b)
ps1.apply(b)

p = Function(W2)
P = p.vector()

solver = KrylovSolver("gmres", "amg")
info(solver.default_parameters(), True)
solver.parameters['monitor_convergence'] = True
solver.parameters['absolute_tolerance'] = 1.0e-6
solver.parameters['gmres']['restart'] = 10000
solver.solve(A, P, b)

pr, pi = p.split()

plot(pr, title='pr')
plot(pi, title='pi')
interactive()

Comments

If the solution has internal resonance, then as you approach the resonant frequency, the matrix will become increasingly ill-conditioned, and it becomes more like an eigenvalue problem. This might explain why it takes longer to converge for decreasing frequencies.

Answer 1

Hi,

I once had a similar problem (with a completely different example). The problem then was that I was trying to use the 'amg' preconditioner for an hyperbolic problem. The solver then stopped to converge for high CFL numbers.

Can you try to change your preconditioner?

Comments

Thanks for your answer.

To be honest I wasn't sure which preconditioner I could use for my problem, because of the special interface conditions I have to use. So I tried the preconditioners listed here and the amg preconditioner was the only one which caused the solver to converge.

For the problem sizes I want to simulate, I need an interative solver. So I'm not sure what I could try now.