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

#1  
New Member
Sergey Lesnik
Join Date: Mar 2015
Posts: 9
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: 330
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: 9
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: 330
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: 330
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: 9
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: 330
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: 9
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. 

June 20, 2016, 05:54 

#10 
New Member
Join Date: Jun 2015
Posts: 24
Rep Power: 3 
Dear Sergey,
I was wondering if you were able to resolve your issue with solving the system. I also need to solve the Helmholtz equation, so just tired out the solver you attached for the 2D case. If I use the GMRES solver it just gets stuck at the first iteration (Time=1). Same goes for BiCGStab. Thanks Ali 

July 5, 2016, 09:13 

#11 
New Member
Sergey Lesnik
Join Date: Mar 2015
Posts: 9
Rep Power: 3 
Hi Ali,
may be it's just very slow... =) How large is your mesh or do you use one of my test cases? What value of nDirections do you use with GMRES solver? 

July 26, 2016, 05:56 

#12 
New Member
Join Date: Jun 2015
Posts: 24
Rep Power: 3 
Hi Sergey,
Thanks a lot for the reply, and sorry for the late response. Well I was looking at your 2d test case. In can't get simpler than this. Q1) Thanks for the hint regarding nDirections. Indeed as mentioned in the Chalmers University's report by Klas Jareteg on block Coupled Calculations in OpenFOAM, the AMG solver is sensitive to both nDir and maxIter. I set it to nDir=10. But now I'm facing the issue you mentioned regarding convergence. So did you figure this out? Q2) I noticed from your examples that you are or were working on sonochemical reactors which is the topic I'm going to look into. From the solver description, I was wondering if you were using the paper of R. Jamshidi et al, Chemical Engr Jr 189190 (2012) 364375? since they also discuss eventually an iterative scheme for calculation the volume fraction and hence the k_m. Q3) Regarding the boundary conditions, I was wondering why would you set the imaginary pressure to same amplitude as the real while it seems that, at first glance, it should be zero ? Thanks, Ali 

July 27, 2016, 09:37 

#13  
New Member
Sergey Lesnik
Join Date: Mar 2015
Posts: 9
Rep Power: 3 
Unfortunately not. It's sufficient for my problem at this point, since I work with stationary (in frequency domain) acoustic fields now. What I figured out is that, the problem has to do more with the preconditioner than with the solver. There are some but none has been implemented in OF yet. E.g. Analytic Incomplete LU (AILU) preconditioner for Krylov methods proposed by M. J. Gander 2000 or a Multigrid one proposed by Y. A. Erlangga 2006, whereby the latter looks much more complex. It sounds like we need to wait until somebody implements it in OF or do it ourselves =)
Exactly, I'm working with his models. Quote:
What I also learned about the BCs meanwhile is that one should use the BC of type "fixedGradient" for sonotrode/transducer/vibrating surface. This one is closer to the physics since it imposes (indirectly) a displacement of the surface, which is what actualy happens in the real world. 

August 10, 2016, 04:52 

#14 
New Member
Join Date: Jun 2015
Posts: 24
Rep Power: 3 
Hi Sergey,
Thanks for the reply. With regards to splitting the pressure amplitude at the source/BC between the real and imaginary part, it is true that the amplitude matters at the end but still slightly unsure about the end result being equivalent to case where the amplitude is entirely imposed on the real part, given the coupled nature of the problem. Perhaps since it is linear problem it wouldn't matter. My background is not acoustics either so ... The easiest way is to consider both and check whether we get the same solution. But the problem now is that we just don't have a stable algorithm to run the cases Also a question regarding the imaginary and real wave numbers values you've included in the transportProperties dict. I just would like to double check I didn't make any errors in the steps to arrive at these values from a given volume fraction of bubbles (Beta), equations 18 in in the paper of Jamshidi. I noticed you always have positive values of the wave numbers for all values of beta, which is not the case for me and can be that both are +,, or the two are of different signs. Say for f=25 kHz, and p0=1 bar=1e5 Pa, and a0=150 microns, and beta=0.1, I get kmSqr.real = 4.4e7, and kmSqr.imag =8.3e6, for which km=6.7e3 as in figure 4.a in the paper. Could you please confirm this using your code, if possible. Also, have you seen this recent ppt from openfoam workshop: http://www.openfoamworkshop.org/asse...es/OFWP053.pdf They've seem to include a new preconditioner for AMG solver, also they're working on new ones as well, so it seems your wish is coming true . Best, Ali Last edited by Ali Blues; August 10, 2016 at 11:11. 

September 28, 2016, 09:22 
Great interest in this project

#15 
Member
ms
Join Date: Mar 2009
Location: West London
Posts: 45
Rep Power: 9 
Hi everyone.
1) Is work still going on on this? 2) Did you consider using a complex solver rather than splitting real and complex? I'm sure I read somewhere that a complex solver had been implemented. I'll look around and post back if I find the solver. Best regards, Mark. 

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 