# How to solve using relaxation techniques

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

 November 25, 2008, 16:17 Hi All, I would like to sol #1 Senior Member   Vishal Nandigana Join Date: Mar 2009 Location: Champaign, Illinois, U.S.A Posts: 208 Rep Power: 10 Hi All, I would like to solve the following equations using relaxation techniqes. As direct solvers are not available in OpenFoam. ddt(C1) = K1 * laplcian(C1) + K2 * laplacian (phi) ddt(C2) = K1 * laplcian(C2) + K2 * laplacian (phi) K3 * laplacian (phi) = K4 * f(C1) + K5 * f(C2) where k1,k2,k3,k4 and k5 are some constant values.. f(C1) and f(C2) means the equation is a function of C1 and C2. I would like to know which tutorial i have to follow to solve these equations. I am fairly new to OpenFoam. K Kindly do the needful help. Thanking You Vishal

 November 29, 2008, 15:59 I have no experience with PNP, #2 Senior Member   Henrik Rusche Join Date: Mar 2009 Location: Braunschweig, Niedersachsen, Germany Posts: 275 Rep Power: 10 I have no experience with PNP, but this might help. http://groups.google.com/group/sci.math.num-analysis/browse_thread/thread/ec5e96 b20472d308?pli=1 To try relaxation methods is easy with OpenFOAM - Discretise each equation, solve one after another and hope for the best. This is your first step and should yield results if the the time step and/or under-relaxation is small enough. Here are a few (uneducated) ideas to improve efficiency, but do step 1 first!: a) Substitute eqn 3 into eqn 1 and eqn 2 to yield 1a and 2a. See whether this is better. b) Substitution of eqn 2a into 1a is not possible, but you might want to create the matrix and substitute C2=H/A into 1a. Good luck, Henrik

 February 2, 2009, 08:03 There are two forms of relaxat #3 Senior Member   Henrik Rusche Join Date: Mar 2009 Location: Braunschweig, Niedersachsen, Germany Posts: 275 Rep Power: 10 There are two forms of relaxation. One is operating and on the field (Phi.relax()) to be used after solving and the other is operating on the matrix (PhiEqn.relax()). The latter is the one you want. What you coded has no effect, it only changes the starting point for the solver. Also make sure to to have a subDict relaxationFactors { Phi 0.7; } in system/fvSolution. The underrelaxation factor is 1.0 if there is no entry. Make sure underrelaxation is working by rerunning the code (after changing relaxation) and see whether the solver output changes. Henrik mm.abdollahzadeh and sharonyue like this.

 February 3, 2009, 00:22 Hi Henrik, Thanks for the r #4 Senior Member   Vishal Nandigana Join Date: Mar 2009 Location: Champaign, Illinois, U.S.A Posts: 208 Rep Power: 10 Hi Henrik, Thanks for the response... So Now the code should be something like this right.... fvScalarMatrix PhiEqn ( fvm::laplacian(Phi) == -alpha * Z1 * C1 - alpha * Z2 * C2 ); PhiEqn.relax(); I would like to know does this relax the equation or the variable Phi.. I meant does this do something like Phi = Phi * (1-relax) + Phi_old (relax) I couldn't exactly comprehend what you replied last time.. Can you throw some more light on this,,, Thanks a lot.. Regards Vishal

 February 3, 2009, 12:56 Hi Vishal, @code: Yes, if y #5 Senior Member   Henrik Rusche Join Date: Mar 2009 Location: Braunschweig, Niedersachsen, Germany Posts: 275 Rep Power: 10 Hi Vishal, @code: Yes, if you don't forget PhiEqn.solve() later. @equation: This is exactly what Phi.relax() does - Explicit relaxation. In implicit relaxation (PhiEqn.relax()) the matrix is modified before solving. The solution is still approached taking it step by step, similar to explicit relaxation, but note that even if you were relaxing only one cell, the effect would be felt throughout the whole domain, in contrast to explicit relaxation. Furthermore, the diagonal dominance of the matrix is increased which is crucial for steady-state calculations where the matrix of a standard transport equation is diagonally equal at best. @exactly comprehend: Please be more specific and will try and help. mm.abdollahzadeh and sharonyue like this.

 February 4, 2009, 17:19 Hi Henrik, Thanks for the r #6 Senior Member   Vishal Nandigana Join Date: Mar 2009 Location: Champaign, Illinois, U.S.A Posts: 208 Rep Power: 10 Hi Henrik, Thanks for the response.. Now I understood.... I just have a small query regarding implementing some boundary condition.. I would like to incorporate a fixed gradient b.c which is a function of gradient of another variable.. for example like I would like to incorporate grad(C1) = Z1*C1*grad(Phi)where Z1 - constant How this type of B.C should be incorporated... I tried this but I don't know why it is not working... The code I worte goes like this... const polyPatchList& patches = mesh.boundaryMesh(); forAll (patches, patchI){ if(C1.boundaryField()[patchI].type()=="fixedGradient"){ fixedGradientFvPatchScalarField& nameOfYourChoice = refCast(C1.bounda ryField()[patchI]); scalarField& otherNameOfYourChoice = nameOfYourChoice.gradient(); otherNameOfYourChoice = -C1*Z1*fvc::snGrad(Phi); } } Kindly let me know what is the mistake I am making... As the error I get is something like this C1 Residual C2 Residual Phi Residual Begin to solve C11#0 Foam::error::printStack(Foam:stream&) in "/home/vishal/OpenFOAM/OpenFOAM-1.5/lib/linuxGccDPOpt/libOpenFOAM.so" #1 Foam::sigFpe::sigFpeHandler(int) in "/home/vishal/OpenFOAM/OpenFOAM-1.5/lib/linuxGccDPOpt/libOpenFOAM.so" #2 Uninterpreted: [0xb7f12420] #3 Foam::divide(Foam::Field&, Foam::UList const&, Foam::UList const&) in "/home/vishal/OpenFOAM/OpenFOAM-1.5/lib/linuxGccDPOpt/libOpenFOAM.so" #4 Foam::operator/(Foam::UList const&, Foam::UList const&) in "/home/vishal/OpenFOAM/OpenFOAM-1.5/lib/linuxGccDPOpt/libOpenFOAM.so" #5 Foam::fixedGradientFvPatchField::evaluate( Foam::Pstream::commsTypes) in "/home/vishal/OpenFOAM/OpenFOAM-1.5/lib/linuxGccDPOpt/libfiniteVolume.so" #6 Foam::GeometricField::GeometricBoundaryField::evaluate() in "/home/vishal/OpenFOAM/vishal-1.5/applications/bin/linuxGccDPOpt/electroosmoticF oam" #7 Foam::fvMatrix::solve(Foam::Istream&) in "/home/vishal/OpenFOAM/OpenFOAM-1.5/lib/linuxGccDPOpt/libfiniteVolume.so" #8 main in "/home/vishal/OpenFOAM/vishal-1.5/applications/bin/linuxGccDPOpt/electroosmoticF oam" #9 __libc_start_main in "/lib/tls/i686/cmov/libc.so.6" #10 Foam::regIOobject::readIfModified() in "/home/vishal/OpenFOAM/vishal-1.5/applications/bin/linuxGccDPOpt/electroosmoticF oam" Floating point exception Kindly do the needful Thanks Regards Vishal

 February 6, 2009, 02:38 Dear Henrik, Thanks for the #7 Senior Member   Vishal Nandigana Join Date: Mar 2009 Location: Champaign, Illinois, U.S.A Posts: 208 Rep Power: 10 Dear Henrik, Thanks for the response... Here is how I have coded my equations fvScalarMatrix C1Eqn ( fvm::ddt(C1) == fvm::laplacian(D1,C1) + (D1 * Z1 * fvc::laplacian(Phi) * fvm::Sp(1.0,C1))+ (D1 * Z1 * (fvc::grad(Phi) & fvc::grad(C1))) ); C1Eqn.relax(); eqnResidualC1 = C1Eqn.solve().initialResidual(); maxResidual = max(eqnResidualC1, maxResidual); Similarly for C2 fvScalarMatrix C2Eqn ( fvm::ddt(C2) == fvm::laplacian(D2,C2) + (D2 * Z2 * fvc::laplacian(Phi) * fvm::Sp(1.0, C2))+ (D2 * Z2 * (fvc::grad(Phi) & fvc::grad(C2))) ); C2Eqn.relax(); eqnResidualC2 = C2Eqn.solve().initialResidual(); maxResidual = max(eqnResidualC2, maxResidual); and for Phi fvScalarMatrix PhiEqn ( fvm::laplacian(Phi) == -alpha * Z1 * C1 - alpha * Z2 * C2 ); PhiEqn.relax(); eqnResidualPhi = PhiEqn.solve().initialResidual(); maxResidual = max(eqnResidualPhi, maxResidual); Buy by using these equations, I don't see any relaxation being carried out even if I change my relaxation factor I don't see any change, further I wanted to solve steady state solution, but I get bizzare results with each iteration... Kindly let me know why is the relaxation not working.. I followed the same way given in the tutorial "BuoyantSimpleRadiation". Kindly give your suggestions Thanks a lot Regards Vishal Thanks Regards Vishal

 February 6, 2009, 05:10 Vishal, I think you should #8 Senior Member   Henrik Rusche Join Date: Mar 2009 Location: Braunschweig, Niedersachsen, Germany Posts: 275 Rep Power: 10 Vishal, I think you should have opened one thread for the complete problem. It also makes it easier for you to find the answers ;-). For interested third parties: BCs are dealt with in this thread http://www.cfd-online.com/cgi-bin/Op...c=1&post=22484 @code: Looks fine to me. Please post the output of the first two iterations for a run with and a run without underrelaxation as well as the relevant section in system/fvSchemes. Henrik

 February 6, 2009, 14:28 Spend an "s" on "relaxationFac #9 Senior Member   Henrik Rusche Join Date: Mar 2009 Location: Braunschweig, Niedersachsen, Germany Posts: 275 Rep Power: 10 Spend an "s" on "relaxationFactor" will make your day. Henrik

 February 6, 2009, 15:10 Dear Henrik, Thanks a lot f #10 Senior Member   Vishal Nandigana Join Date: Mar 2009 Location: Champaign, Illinois, U.S.A Posts: 208 Rep Power: 10 Dear Henrik, Thanks a lot for the suggestion, now I see the change.. But still the solution is blowing up with each iteration. I see that the residual is changing but it is not helping in giving me a converged solution. I am solving a steady state problem but still for every time step the error is picking up and I get NAN at the end... Kindly let me know where I am making an error... Regards Vishal

 February 6, 2009, 17:19 Yes, now you are at right at t #11 Senior Member   Henrik Rusche Join Date: Mar 2009 Location: Braunschweig, Niedersachsen, Germany Posts: 275 Rep Power: 10 Yes, now you are at right at the beginning where the becomes interesting: a) Why do you have ddt-terms if you are seeking a steady-state solution? b) Change all relaxation factors to 0.001 and have a look how the fields are evolving. Where is it blowing up? What's causing the trouble? c) fvc::grad(Phi) & fvc::grad(C2) is nasty. Bring this into a conservative form by using div(s*v) = s*div(v) + grad(s) . v d) Just cosmetics, but it hurts my eyes: (D2 * Z2 * fvc::laplacian(Phi) * fvm::Sp(1.0, C2)) is equal to fvm::Sp(D2 * Z2 * fvc::laplacian(Phi), C2) Henrik mm.abdollahzadeh likes this.

 February 7, 2009, 06:37 a) Okay b) Okay, I understand #13 Senior Member   Henrik Rusche Join Date: Mar 2009 Location: Braunschweig, Niedersachsen, Germany Posts: 275 Rep Power: 10 a) Okay b) Okay, I understand that the solution has bounds (C1 > 0 C2 > 0) which are violated c) One hint: s = Phi and v = grad(C2). One more thing: Do you understand what the conservative form of a transport equation is? If not, look it up. Any textbook on FVM will do. There is no chance you will solve this problem if you do not understand the basics - especially in the light of b) d) Try suggestion a) from my first quote. This gives you two conservative equations

 February 7, 2009, 16:07 Vishal, @I) I tried to sugg #14 Senior Member   Henrik Rusche Join Date: Mar 2009 Location: Braunschweig, Niedersachsen, Germany Posts: 275 Rep Power: 10 Vishal, @I) I tried to suggest the following: Replace 0 = K1 * laplacian(C1) + K2 * laplacian (phi) 0 = K1 * laplacian(C2) + K2 * laplacian (phi) K3 * laplacian (phi) = beta(C1, C2) beta (C1,C2) = K4 * f(C1) + K5 * f(C2) by 0 = K1 * laplacian(C1) + K2/K3*beta(C1,C2) 0 = K1 * laplacian(C2) + K2/K3*beta(C1,C2) This should be easier to solve because you have fewer equations and they are conservative. But you need to be careful when discretising the source to ensure boundedness at zero. What do f(C1) and f(C2) look like? Henrik

 February 7, 2009, 16:33 Dear Henrik, The equations #15 Senior Member   Vishal Nandigana Join Date: Mar 2009 Location: Champaign, Illinois, U.S.A Posts: 208 Rep Power: 10 Dear Henrik, The equations which I have quoted in the beginning of the forum is not the equations I am looking to solve... The exact equations which I would like to solve are div {(D1*grad(C1) + D1*Z1*C1*grad(Phi))} = 0 -(1) div {(D2*grad(C2) + D2*Z2*C2*grad(Phi))} = 0 -(2) laplacian (phi) = - alpha * Z1*C1 - alpha*Z2 * C2 - (3) where C1 C2 and Phi are concentrations and potential respectively and D1,D2,Z1,Z2, alpha are constants Sorry for the confusion. I am trying to solve these equations, Initially I have tried to expand each term on the LHS of equation 1/2 and wrote them as laplacian terms.. When I solve them like that then after every iteration I see a change in my C1 and C2 and eventually blowing out my solution... Hope I have made myself clear now.. I am again sorry for the confusions.. Thanks Regards Vishal

 February 8, 2009, 06:22 Shouldn't this be written as: #16 Senior Member     Dragos Join Date: Mar 2009 Posts: 648 Rep Power: 12 Shouldn't this be written as: fvScalarMatrix Eqn1( fvm::div(D1*fvc::grad(C1)+D1*Z1*C1*fvc::grad(Phi)) ); fvScalarMatrix Eqn2( fvm::div(D2*fvc::grad(C2)+D2*Z2*C2*fvc::grad(Phi)) ); fvScalarMatrix Eqn3( fvm::laplacian(Phi)+alpha*Z1*C1+alpha*Z2*C2 ); Eqn1.relax(); Eqn2.relax(); Eqn3.relax(); Eqn1.solve(); Eqn2.solve(); Eqn3.solve(); Dragos

 February 8, 2009, 11:26 Your entry div(((Z1*C1)*fvc:: #17 Senior Member     Dragos Join Date: Mar 2009 Posts: 648 Rep Power: 12 Your entry div(((Z1*C1)*fvc::grad(Phi))) Gauss upwind; in the dictionary is wrong. Write instead div(((Z1*C1)*grad(Phi))) Gauss upwind; On the other hand, I don't know why it didn't work in the previous formulation. Dragos

 February 8, 2009, 11:51 It seems there is an extra spa #18 Senior Member     Dragos Join Date: Mar 2009 Posts: 648 Rep Power: 12 It seems there is an extra space in div(((Z1*C1)*grad(Ph i)))

 February 8, 2009, 12:05 Dear Dragos, No I just chec #19 Senior Member   Vishal Nandigana Join Date: Mar 2009 Location: Champaign, Illinois, U.S.A Posts: 208 Rep Power: 10 Dear Dragos, No I just checked it again... There is no space fvm::laplacian(C1) + fvc::div(((Z1*C1)*fvc::grad(Phi))) divSchemes { default none; div(((Z1*C1)*grad(Phi))) Gauss upwind; div(((Z2*C2)*grad(Phi))) Gauss upwind; } Calculating concentration distribution Time = 0.2 C1 Residual C2 Residual Phi Residual attempt to read beyond EOF file: /home/vishal/OpenFOAM/OpenFOAM-1.5/tutorials/coupledFoam_2dSteady_majumdar_geome try/coupledFoam2dSteady_majumdar_geometry/system/fvSchemes::div(((Z1*C1)*grad(Ph i))) at line 31. From function ITstream::read(token& t) in file db/IOstreams/Tstreams/ITread.C at line 64. FOAM exiting But it is displayed like this... when I post it.. Regards Vishal

 February 9, 2009, 02:30 Hello guys, Can anyone tell m #20 Senior Member     Dragos Join Date: Mar 2009 Posts: 648 Rep Power: 12 Hello guys, Can anyone tell me why the following expression gives me a compiler syntax error? fvScalarMatrix Eqn1( fvm::div(D1*fvc::grad(C1)+D1*Z1*C1*fvc::grad(Phi)) ); Dragos

 Thread Tools Display Modes Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post samar FLUENT 1 November 9, 2007 05:15 PK FLUENT 2 February 6, 2007 17:44 Ugur FLUENT 0 September 21, 2004 10:05 Mihail Galabov Main CFD Forum 1 April 23, 2004 11:57 phil Main CFD Forum 3 November 26, 2002 08:04

All times are GMT -4. The time now is 01:59.