
[Sponsors] 
October 5, 2015, 10:25 
Creating Helmholtz Solver: no convergence

#1  
New Member
Sergey Lesnik
Join Date: Mar 2015
Posts: 6
Rep Power: 3 
Hello everybody!
I tried to implement a simple solver for Helmholtz equation (wave equation in frequency domain) of the form: Quote:
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 ); 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. Setup: I tested the code on 1D benchmark with the following conditions:
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:
I've even tried the block matrix solver from the openfoamextend which is suitable for such coupled problems  no improvement. Next thing I want to try is the Finite Element solver from foamextend. 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. 

October 5, 2015, 13:19 

#2 
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 329
Rep Power: 11 
Hi, maybe you must try to present p as vector instead of 2 scalars?


October 6, 2015, 05:23 

#3 
New Member
Sergey Lesnik
Join Date: Mar 2015
Posts: 6
Rep Power: 3 
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 HTML Code:
Δp_Re + K_Re ⋅ p_Re = K_Im ⋅ p_Im Δp_Im + K_Re ⋅ p_Im =  K_Im ⋅ p_Re 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.11629e24 5.10652e24), No Iterations 1 ExecutionTime = 0.75 s ClockTime = 0 s BiCGStab: Solving for blockT, Initial residual = (0.017386 0.0209848), Final residual = (4.9824e24 5.08883e24), No Iterations 1 ExecutionTime = 0.77 s ClockTime = 0 s 

October 6, 2015, 11:41 

#4 
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 329
Rep Power: 11 
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 = .... } 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. 

October 7, 2015, 08:07 

#5 
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 329
Rep Power: 11 
Hi,
I found that solution doesn't converge for cases with nonzero 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)? 

October 16, 2015, 08:36 

#6 
New Member
Sergey Lesnik
Join Date: Mar 2015
Posts: 6
Rep Power: 3 
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 foamextend3.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 offdiagonal 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 offdiagonal 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 offdiagonal 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();
In the attached file you can find the solver (you need foamextend3.2 to compile it) and the mentioned test cases, if you want to try it out and get it to converge (faster) =)) 

October 16, 2015, 13:14 

#7 
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 329
Rep Power: 11 
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 

October 19, 2015, 11:25 

#8 
New Member
Vladimir Platonov
Join Date: Oct 2015
Posts: 1
Rep Power: 0 
Hi,
which preconditioner you're using with BiCGStab? You could try using GAMG (in segregated case) or AMG with Fcycle in coupled case. 

October 20, 2015, 08:28 

#9  
New Member
Sergey Lesnik
Join Date: Mar 2015
Posts: 6
Rep Power: 3 
Quote:
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 Multigridsolvers are not suitable for wave equations (e.g. here chapter 2.5), which seems logical to me since the "wavelike" structure get lost on a coarser grid. I might try the Multigrid Solver from foamextend, but from the considerations above I'm doubting it'll be any better. 

Tags 
coupled equations, frequency domain, helmholtz, openfoam 
Thread Tools  
Display Modes  


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 densitybased 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 