CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

Creating Helmholtz Solver: no convergence

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   October 5, 2015, 10:25
Unhappy Creating Helmholtz Solver: no convergence
  #1
New Member
 
Sergey Lesnik
Join Date: Mar 2015
Posts: 6
Rep Power: 3
serles is on a distinguished road
Hello everybody!

I tried to implement a simple solver for Helmholtz equation (wave equation in frequency domain) of the form:

Quote:
laplace(p) + k^2 * p = 0
where p is the acoustic pressure and k the complex wave number. Since OpenFOAM doesn't support complex numbers I decomposed the equations in two (introducing p = p_Re + i*p_Im and same for k) and implemented them like this:
Code:
solve
                (
                        fvm::laplacian(p_Re)
                        + fvm::Sp(k_Re_sq, p_Re)
                        == k_Im_sq * p_Im
                );

            solve
                (
                        fvm::laplacian(p_Im)
                        + fvm::Sp(k_Re_sq, p_Im)
                        == - k_Im_sq * p_Re
                );
where Im and Re stand for imaginary and real.

As one can see the problem has to be solved iteratively: first for Re(p), then for the Im(p) using the result from the previous equation and afterwards repeat from the beginning until both solutions converge.

Set-up: I tested the code on 1D benchmark with the following conditions:
  • wave number is constant for the hole domain
  • BC of the acoustic source: Dirichlet with Re(p)=const and Im(p)=const
  • BC of the opposite end: Dirichlet with Re(p)=Im(p)=0
  • other boundaries are defined as empty (since it's 1D)

Result: The solution doesn't converge for several test cases!
It converges for small frequencies and small wave numbers, but not for higher values. The pressure starts to grow with each iteration step. If it converges, I get good results which perfectly match analytical solution or results from other software.

Possible solutions: I've tried all the numerical tweaks I could find/use in OF:
  • different linear solvers (btw GAMG gives an error in almost all of the test cases)
  • different interpolation/laplacian/snGrade schemes
  • resolution of the domain from 100 cells to 10^6
  • relaxation of the equations with different relaxation factors

I've even tried the block matrix solver from the openfoam-extend which is suitable for such coupled problems - no improvement.
Next thing I want to try is the Finite Element solver from foam-extend. I suppose this will bring problems when I start to couple acoustics with fluid dynamics...

But I truly cannot believe that Finite Volume Method is not capable to solve it. I'm not a specialist in CFD, so may be I don't see some simple mistake
Does somebody have any suggestions what else I can try?

I'm open for any help or discussion.
serles is offline   Reply With Quote

Old   October 5, 2015, 13:19
Default
  #2
Senior Member
 
mkraposhin's Avatar
 
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 329
Rep Power: 11
mkraposhin is on a distinguished road
Hi, maybe you must try to present p as vector instead of 2 scalars?
mkraposhin is offline   Reply With Quote

Old   October 6, 2015, 05:23
Default
  #3
New Member
 
Sergey Lesnik
Join Date: Mar 2015
Posts: 6
Rep Power: 3
serles is on a distinguished road
Hi, thank you for the input!

But I don't see a way how to introduce such a vector in given equations.
Given the Helmholtz eqn. introduce complex p and k:
HTML Code:
Δp + k^2 ⋅ p = 0
p = p_Re + i ⋅ p_Im
k^2 = K_Re + i ⋅ K_Im
Gives two coupled eqns:
HTML Code:
Δp_Re + K_Re ⋅ p_Re = K_Im ⋅ p_Im
Δp_Im + K_Re ⋅ p_Im = - K_Im ⋅ p_Re
where the second eqn. is purely imaginary.
If I introduce a vector p_V = (p_Re; p_Im) what am I going to do with the RHS?

If I understood right, something similar, what you mean, the block matrix from the extend project is doing. It incorporates p_Re and p_Im as a vector and solves the equations implicitly. But in my case it doesn't help:
Code:
Time = 1

BiCGStab:  Solving for blockT, Initial residual = (0.018707 0.0155012), Final residual = (2.11629e-24 5.10652e-24), No Iterations 1
ExecutionTime = 0.75 s  ClockTime = 0 s

BiCGStab:  Solving for blockT, Initial residual = (0.017386 0.0209848), Final residual = (4.9824e-24 5.08883e-24), No Iterations 1
ExecutionTime = 0.77 s  ClockTime = 0 s
It converges perfectly during a single iteration but the residuals are large again when the new solution is introduced to the solver.
serles is offline   Reply With Quote

Old   October 6, 2015, 11:41
Default
  #4
Senior Member
 
mkraposhin's Avatar
 
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 329
Rep Power: 11
mkraposhin is on a distinguished road
Hi,

ok, now i understand what you are doing.

I would suggest three improvements for weakly coupled iterative solution scheme
1) Use relaxation
2) Use hybrid implicit / explicit source term.
3) Update all sources only once before solving equations

for example:

Code:
bool converged = false;
while (converged)
{

    fvScalarMatrix ReEqn
    (
         fvm::laplacian(p_Re)
          + fvm::SuSp(k_Re_sq, p_Re)
          == k_Im_sq * p_Im
    );

    fvScalarMatrix ImEqn
    (
         fvm::laplacian(p_Im)
         + fvm::SuSp(k_Re_sq, p_Im)
         == - k_Im_sq * p_Re
    );

    ReEqn.relax();
    ImEqn.relax();

    ReEqn.solve();
    ImEqn.solve();

    //converged = ....
}
Oh, I'm sorry, i read your first message again and found, that you tried relaxation already. But nevertheless, i think that relaxation will stablize convergence process.
Maybe you need to play with relaxation coefficients, from 0.5 down to 0.05

I hope, this will help

Last edited by mkraposhin; October 6, 2015 at 15:53.
mkraposhin is offline   Reply With Quote

Old   October 7, 2015, 08:07
Default
  #5
Senior Member
 
mkraposhin's Avatar
 
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 329
Rep Power: 11
mkraposhin is on a distinguished road
Hi,

I found that solution doesn't converge for cases with non-zero imaginary part of wave number,

you can see example and solver in the attachment

maybe you need different B.C. or treatment for the source due to Im(k^2)?
Attached Files
File Type: gz helmholtzFoam.tar.gz (3.1 KB, 6 views)
mkraposhin is offline   Reply With Quote

Old   October 16, 2015, 08:36
Default
  #6
New Member
 
Sergey Lesnik
Join Date: Mar 2015
Posts: 6
Rep Power: 3
serles is on a distinguished road
Hi,
thank you for the interest on the problem. I've looked at your solver and it is exactly what I was trying to do. As I found out solving the Helmholtz equation is not such a trivial task: Why it is Difficult to Solve Helmholtz Problems with Classical Iterative Methods. Not all the solvers are suggested and special preconditioners are needed.


Thus, I persued the approch with the Coupled Block Matrix solver and found out that I initialized the block matrix in a wrong way. In the first try I did it like in foam-extend-3.2/applications/solvers/coupled/blockCoupledScalarTransportFoam. They insert the implicit terms in the diagonal parts of the block matrix but not in the upper and lower parts, which I don't understand why.

After some time I found this thread: Solving Coupled Set of PDEs, where Tron has done exactly what I needed. It looks now like this:
Code:
 
// separate 1st eqn in two parts and store them in separate matrices
// 2nd matrix serves only to insert the coupled terms into the 1st one
fvScalarMatrix Ap_ReM
(
    fvm::laplacian(Ap_Re) + fvm::Sp(k_Re_sq, Ap_Re)
);

Ap_ReM.relax();

fvScalarMatrix Ap_ImM2
(
    -fvm::Sp(k_Im_sq, Ap_Im)
);

Ap_ImM2.relax();

// separate 2nd eqn in two parts and store them in separate matrices
fvScalarMatrix Ap_ImM
(
    fvm::laplacian(Ap_Im) + fvm::Sp(k_Re_sq, Ap_Im)
);

Ap_ImM.relax();

fvScalarMatrix Ap_ReM2
(
    fvm::Sp(k_Im_sq, Ap_Re)
);

Ap_ReM2.relax();

// Prepare block system
BlockLduMatrix<vector2> blockM(mesh);
BlockLduMatrix<vector2> blockM2(mesh);

// Grab block diagonal and set it to zero
Field<tensor2>& d = blockM.diag().asSquare();
d = tensor2::zero;
Field<tensor2>& d2 = blockM2.diag().asSquare();
d2 = tensor2::zero;

// Grab linear off-diagonal and set it to zero
Field<tensor2>& u = blockM.upper().asSquare();
Field<tensor2>& l = blockM.lower().asSquare();
u = tensor2::zero;
l = tensor2::zero;
Field<tensor2>& u2 = blockM2.upper().asSquare();
Field<tensor2>& l2 = blockM2.lower().asSquare();
u2 = tensor2::zero;
l2 = tensor2::zero;

vector2Field& blockX = blockT.internalField();
vector2Field& blockX2 = blockT2.internalField();
vector2Field blockB(mesh.nCells(), vector2::zero);
vector2Field blockB2(mesh.nCells(), vector2::zero);

//- Insert equations into block Matrix
blockMatrixTools::insertEquation(0, Ap_ReM, blockM, blockX, blockB);
blockMatrixTools::insertEquation(1, Ap_ImM, blockM, blockX, blockB);
blockMatrixTools::insertEquation(0, Ap_ReM2, blockM2, blockX2, blockB2);
blockMatrixTools::insertEquation(1, Ap_ImM2, blockM2, blockX2, blockB2);

//- Add off-diagonal terms and remove from block source for the diagonal
forAll(d, i)
{
    d[i](0,1) = d2[i](1,1);
    d[i](1,0) = d2[i](0,0);
    blockB[i][0] -= blockB2[i][1];
    blockB[i][1] -= blockB2[i][0];
}

//- Add off-diagonal terms and remove from block source upper and lower
forAll(u, j)
{
    u[j](0,1) = u2[j](1,1);
    u[j](1,0) = u2[j](0,0);
}

forAll(l, j)
{
    l[j](0,1) = l2[j](1,1);
    l[j](1,0) = l2[j](0,0);
}

//- Block coupled solver call for Matrix1 with inserted terms
BlockSolverPerformance<vector2> solverPerf =
    BlockLduSolver<vector2>::New
    (
        blockT.name(),
        blockM,
        mesh.solutionDict().solver(blockT.name())
    )->solve(blockX, blockB);

solverPerf.print();

// Retrieve solution
blockMatrixTools::blockRetrieve(0, Ap_Re.internalField(), blockX);
blockMatrixTools::blockRetrieve(1, Ap_Im.internalField(), blockX);

Ap_Re.correctBoundaryConditions();
Ap_Im.correctBoundaryConditions();
Anyway, my solver with this implementation works nicely for the 1D case now, but I have still some convergence problems with 2D and 3D for the desired wave numbers. There are two linear solvers I tried out: BiCGStab and GMRES. These are my issues:
  1. It doesn't matter what under-relaxation coefficients for the matrices (Ap_ReM, Ap_ImM, ) I use with both solvers. The residuals don't change at all.
  2. BiCGStab converges fast but it's unstable (what an irony). Residuals may go down to 1e-5 and then shoot to the sky.
  3. GMRES is much more stable but it converges very slow. It's better with higher value of nDirections but it needs ages to converge.
Does someone have an idea about these issues or do you have a tip to try out? The problem about under-relaxation concerns me mainly.


In the attached file you can find the solver (you need foam-extend-3.2 to compile it) and the mentioned test cases, if you want to try it out and get it to converge (faster) =))
Attached Files
File Type: gz blockCoupledHelmholtzFoam.tar.gz (9.2 KB, 10 views)
serles is offline   Reply With Quote

Old   October 16, 2015, 13:14
Default
  #7
Senior Member
 
mkraposhin's Avatar
 
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 329
Rep Power: 11
mkraposhin is on a distinguished road
Hi,
i have several colleagues, who are working in the same or close area. I shall try to interest them
__________________
Winter OpenFOAM days in Moscow - http://www.opencloudconf.ru/en/oscome.html
mkraposhin is offline   Reply With Quote

Old   October 19, 2015, 11:25
Default
  #8
New Member
 
Vladimir Platonov
Join Date: Oct 2015
Posts: 1
Rep Power: 0
through_horizon is on a distinguished road
Hi,
which preconditioner you're using with BiCGStab? You could try using GAMG (in segregated case) or AMG with F-cycle in coupled case.
through_horizon is offline   Reply With Quote

Old   October 20, 2015, 08:28
Default
  #9
New Member
 
Sergey Lesnik
Join Date: Mar 2015
Posts: 6
Rep Power: 3
serles is on a distinguished road
Quote:
Originally Posted by through_horizon View Post
Hi,
which preconditioner you're using with BiCGStab? You could try using GAMG (in segregated case) or AMG with F-cycle in coupled case.
Hi,
I've tried all the available Preconditioners with BiCGStab, Cholesky looked for as the most sufficient.
I've tried GAMG as well, but it couldn't accomplish any iteration. I also found some works telling that Multigrid-solvers are not suitable for wave equations (e.g. here chapter 2.5), which seems logical to me since the "wave-like" structure get lost on a coarser grid. I might try the Multigrid Solver from foam-extend, but from the considerations above I'm doubting it'll be any better.
serles is offline   Reply With Quote

Reply

Tags
coupled equations, frequency domain, helmholtz, openfoam

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
solver ends with code 255 after convergence meow CFX 2 October 28, 2014 20:47
convergence of density-based solver for unsteady flow zhengjg FLUENT 0 June 16, 2014 10:16
Convergence parameter CFX Solver chiragsvnit CFX 2 March 17, 2014 01:45
Naca 0012 (compressible and inviscid) flow convergence problem bipulsaha FLUENT 1 July 6, 2011 07:51
Solver stoped without creating a result file Roland CFX 1 September 4, 2006 13:09


All times are GMT -4. The time now is 12:25.