Hi All,
I would like to sol
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 
I have no experience with PNP,
I have no experience with PNP, but this might help.
http://groups.google.com/group/sci.math.numanalysis/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 underrelaxation 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 
There are two forms of relaxat
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 
Hi Henrik,
Thanks for the r
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 * (1relax) + 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 
Hi Vishal,
@code: Yes, if y
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 steadystate 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. 
Hi Henrik,
Thanks for the r
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<fixedgradientfvpatchscalarfield>(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:http://www.cfdonline.com/OpenFOAM_D...part/proud.gifstream&) in "/home/vishal/OpenFOAM/OpenFOAM1.5/lib/linuxGccDPOpt/libOpenFOAM.so" #1 Foam::sigFpe::sigFpeHandler(int) in "/home/vishal/OpenFOAM/OpenFOAM1.5/lib/linuxGccDPOpt/libOpenFOAM.so" #2 Uninterpreted: [0xb7f12420] #3 Foam::divide(Foam::Field<double>&, Foam::UList<double> const&, Foam::UList<double> const&) in "/home/vishal/OpenFOAM/OpenFOAM1.5/lib/linuxGccDPOpt/libOpenFOAM.so" #4 Foam::operator/(Foam::UList<double> const&, Foam::UList<double> const&) in "/home/vishal/OpenFOAM/OpenFOAM1.5/lib/linuxGccDPOpt/libOpenFOAM.so" #5 Foam::fixedGradientFvPatchField<double>::evaluate( Foam::Pstream::commsTypes) in "/home/vishal/OpenFOAM/OpenFOAM1.5/lib/linuxGccDPOpt/libfiniteVolume.so" #6 Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>::GeometricBoundaryField::evaluate() in "/home/vishal/OpenFOAM/vishal1.5/applications/bin/linuxGccDPOpt/electroosmoticF oam" #7 Foam::fvMatrix<double>::solve(Foam::Istream&) in "/home/vishal/OpenFOAM/OpenFOAM1.5/lib/linuxGccDPOpt/libfiniteVolume.so" #8 main in "/home/vishal/OpenFOAM/vishal1.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/vishal1.5/applications/bin/linuxGccDPOpt/electroosmoticF oam" Floating point exception Kindly do the needful Thanks Regards Vishal 
Dear Henrik,
Thanks for the
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 
Vishal,
I think you should
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.cfdonline.com/cgibin/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 
Spend an "s" on "relaxationFac
Spend an "s" on "relaxationFactor" will make your day.
Henrik 
Dear Henrik,
Thanks a lot f
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 
Yes, now you are at right at t
Yes, now you are at right at the beginning where the becomes interesting:
a) Why do you have ddtterms if you are seeking a steadystate 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 
Dear Henrik,
Thanks for the
Dear Henrik,
Thanks for the reply. "a) Why do you have ddtterms if you are seeking a steadystate solution? " Though I used ddt terms, I mentioned steadyState in my fvschemes, it should have the same effect right. "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? " The reason for the blowing up of the solution is because of the changes in C1 and C2 with each iteration... The C1 and C2 initially are small (inlet and outlet = 0.1 mM).. and with each iteration they are increasing.. Initially they are around 0.5 mM for C1 and C2 is reaching negative values around 0.4 mM (I don't want this to happen), then after few iterations they are just blowing out reaching 10000 mM and all... "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 " I am not able to understand how exactly I have to write instead of fvc::grad(Phi) & fvc::grad(C2) in OpenFoam "d)(D2 * Z2 * fvc::laplacian(Phi) * fvm::Sp(1.0, C2)) is equal to fvm::Sp(D2 * Z2 * fvc::laplacian(Phi), C2) " I did change this... Basically I have my equations as tau1 = D1*grad(C1)  D1*Z1*C1*grad(Phi) tau 2 = D2*grad(C2)  D2*Z2*C2*grad(Phi) where C1 C2 and Phi are concentrations and potential respectively D1,D2,Z1,Z2 are constants This is the set of coupled equations I am solving... dC1/dt =  div (tau1) dC2/dt =  div (tau2) For Steady State, the left hand side is 0. laplacian (phi) =  alpha * Z1*C1  alpha*Z2 * C2 where alpha is a constant I have expanded the div(tau1) and similarly div(tau2) and wrote everything as laplacian and gradient terms.. So I don't have any divergence terms as such hence even if I get negative values, I am not able to rectify them. Kindly give your suggestions as to how I should go about solving these issues and whether my equations are correct for solving this kind of problem. Thanks Regards Vishal 
a) Okay
b) Okay, I understand
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 
Vishal,
@I) I tried to sugg
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 
Dear Henrik,
The equations
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 
Shouldn't this be written as:
Shouldn't this be written as:
Dragos 
Your entry div(((Z1*C1)*fvc::
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 
It seems there is an extra spa
It seems there is an extra space in div(((Z1*C1)*grad(Ph i)))

Dear Dragos,
No I just chec
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/OpenFOAM1.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 
Hello guys,
Can anyone tell m
Hello guys,
Can anyone tell me why the following expression gives me a compiler syntax error?
Dragos 
All times are GMT 4. The time now is 20:55. 