
[Sponsors] 
November 23, 2009, 15:31 
basic question with 'ForAll' loop

#1 
Member
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 10 
Hi all,
Strange things happen in my forAll loop. I have those two variables 'n' and 'somme_E' that are reinitialized at the end of each iteration of the for loop (which is confirmed by using the 'info' command), but when those same variables are reused inside the forAll loop (at the next iteration) they keep their previous value (before they were reinitialized) . I wonder why. I have no error message just the wrong output. Can anybody help me? Thank you very much, Pascal Code:
scalar coord_z = 0.0; scalar coord_z_dernier = 0.0; scalar diff = 0.0; scalar plan_xy = 1./128.; scalar somme_E = 0.0; scalar dz = 1./64.; int n = 0; for (int iter = 1; iter <= 64; iter ++) { forAll(U, cellI) { coord_z = mesh.C()[cellI].component(2); diff = coord_z  plan_xy; if (diff <= 0.00002) { somme_E = somme_E + mag(U[cellI]) * mag(U[cellI]); coord_z_dernier = mesh.C()[cellI].component(2); n = n+1; } } Info<< coord_z_dernier << " " << somme_E << " " << n << endl; plan_xy = plan_xy + dz; somme_E = 0.; n = 0; Info<< n << " " << somme_E << endl; } 

November 24, 2009, 03:47 

#2  
Senior Member
Mark Olesen
Join Date: Mar 2009
Location: http://olesenm.github.io/
Posts: 800
Rep Power: 21 
Quote:
I would suspect something else is going wrong with your logic. However, why are you reinitializing things anyhow instead of just using scoped variables? For example (as pseudocode, without ANY promise that it does what you really want  or even if it compiles), Code:
scalar plan_xy = 1./128.; const scalar dz = 1./64.; for (int iter = 1; iter <= 64; iter ++) { label n = 0; scalar somme_E = 0.0; scalar coord_z_dernier = 0.0; // debug only ??? forAll(U, cellI) { scalar coord_z = mesh.C()[cellI].component(2); scalar diff = mag(coord_z  plan_xy); if (diff <= 0.00002) { somme_E += magSqr(U[cellI])); coord_z_dernier = coord_z; n++; } } Info<< coord_z_dernier << " " << somme_E << " " << n << endl; plan_xy += dz; } 

December 13, 2012, 16:14 

#3 
New Member
Join Date: Mar 2012
Posts: 29
Rep Power: 7 
Did you get this problem solved? I am having exactly the same issue. The forAll loop does not update the variable at all.


December 13, 2012, 16:25 

#4 
Member
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 10 
Hello Yanxiang,
That make a long time ago I think it was a simple error in my if statement: Code:
if (diff <= 0.00002) Code:
if (abs(diff) <= 0.00002) 

December 13, 2012, 16:41 

#5 
New Member
Join Date: Mar 2012
Posts: 29
Rep Power: 7 
Thanks, Pascal. Now it looks like I have a different issue. So I have the following code
Code:
forAll (beta, celli) { if (neg(beta[celli])) { beta[celli] = 0; } } yanxiang 

December 13, 2012, 16:55 

#6 
Member
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 10 
Hi,
Your sample of code seems all right. Could you post the entire code? It will be easier for others to give help. Pascal 

December 13, 2012, 17:02 

#7 
New Member
Join Date: Mar 2012
Posts: 29
Rep Power: 7 
Hi Pascal,
Basically I am modifying the twoPhaseEulerFoam code, adding a couple of lines to the alphaEqn.H. Here is what I have. The above code is at the end. Thanks a lot for your help in advance Code:
{ word scheme("div(phi,alpha)"); word schemer("div(phir,alpha)"); surfaceScalarField phic("phic", phi); surfaceScalarField phir("phir", phia  phib); if (g0.value() > 0.0) { surfaceScalarField alphaf(fvc::interpolate(alpha)); surfaceScalarField phipp(ppMagf*fvc::snGrad(alpha)*mesh.magSf()); phir += phipp; phic += fvc::interpolate(alpha)*phipp; } for (int acorr=0; acorr<nAlphaCorr; acorr++) { fvScalarMatrix alphaEqn ( fvm::ddt(alpha) + fvm::div(phic, alpha, scheme) + fvm::div(fvc::flux(phir, beta, schemer), alpha, schemer) + fvm::div(fvc::flux(phia, alphas, schemer), alpha, schemer) ); if (g0.value() > 0.0) { ppMagf = rUaAf*fvc::interpolate ( (1.0/(rhoa*(alpha + scalar(0.0001)))) *g0*min(exp(preAlphaExp*(alpha  alphaMax)), expMax) ); alphaEqn = fvm::laplacian ( (fvc::interpolate(alpha) + scalar(0.0001))*ppMagf, alpha, "laplacian(alphaPpMag,alpha)" ); } alphaEqn.relax(); alphaEqn.solve(); #include "packingLimiter.H" beta = scalar(1)  alpha  alphas; // =============== forAll (beta, celli) { if (neg(beta[celli])) { beta[celli] = 0; } } // =============== beta.correctBoundaryConditions(); Info<< "Liquid phase volume fraction = " << alpha.weightedAverage(mesh.V()).value() << " Min(alpha) = " << min(alpha).value() << " Max(alpha) = " << max(alpha).value() << endl; Info<< "Gas phase volume fraction = " << beta.weightedAverage(mesh.V()).value() << " Min(beta) = " << min(beta).value() << " Max(beta) = " << max(beta).value() << endl; } } 

December 13, 2012, 17:13 

#8 
Member
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 10 
Hi!
I think that this line of code : Code:
beta.correctBoundaryConditions(); Pascal 

December 13, 2012, 17:25 

#9 
New Member
Join Date: Mar 2012
Posts: 29
Rep Power: 7 
Hi Pascal,
I tried your suggestions and I also ommented out that line. Nothing really changes. I still got negative values and I think that's the reason why the solution is not stable and blows up at some point. Thanks, yanxiang 

December 14, 2012, 16:04 

#10 
Senior Member
Kent Wardle
Join Date: Mar 2009
Location: Illinois, USA
Posts: 214
Rep Power: 14 
Yanxiang,
I don't understand why you need to do this as a loop over all cells. Can't you take advantage of the 'field operations' part of FOAM and just do a max() on beta as in: beta = max(beta, 0.0); Have you tried this instead? Should do the same thing but with one line of code. Also, have you looked at where in the domain your negative values are occurring? I guess alphas is alpha_solid (and not alphas as defined in mpEF which is 0*alpha1 + 1*alpha2 + 2*alpha3 + ... (n1)*alphan)? If you are getting negative values of beta, that means that alpha+alphas is > 1 in places. I see that you have also changed the solution scheme and are not using MULES::explicitSolve. Note that MULES has a limiter built in to keep the sum of volume fractions equal to 1. With out this it is not surprising you are getting inaccurate solution of your phase fraction fields. This was important change in the development of mpEF. As to how to fix it here...not sure I can help there. If you indeed need to keep three volume fractions, perhaps mpEF needs another look. 

December 14, 2012, 18:39 

#11  
New Member
Join Date: Mar 2012
Posts: 29
Rep Power: 7 
Quote:
Quote:
Thanks, yanxiang 

Tags 
forall loop 
Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Problem with Gmsh  nishant_hull  Open Source Meshers: Gmsh, Netgen, CGNS, ...  23  August 5, 2015 02:09 
OpenFOAM basic Question  a_yoshi  OpenFOAM  1  October 29, 2009 05:13 
Basic Question  a_yoshi  OpenFOAM  2  October 29, 2009 01:45 
Basic Poiseuille flow simulation question  Ashish  Main CFD Forum  0  October 2, 2007 13:05 
NACA0012 geometry/design software needed  Franny  Main CFD Forum  13  July 7, 2007 15:57 