CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Post-Processing

basic question with 'ForAll' loop

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

Like Tree1Likes
  • 1 Post By kwardle

Reply
 
LinkBack Thread Tools Display Modes
Old   November 23, 2009, 15:31
Default basic question with 'ForAll' loop
  #1
Member
 
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 8
Pascal_doran is on a distinguished road
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 re-used 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;
    }
Pascal_doran is offline   Reply With Quote

Old   November 24, 2009, 03:47
Default
  #2
Senior Member
 
Mark Olesen
Join Date: Mar 2009
Location: http://olesenm.github.io/
Posts: 777
Rep Power: 18
olesen will become famous soon enough
Quote:
Originally Posted by Pascal_doran View Post
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 re-used 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,
If you take a look at the file src/OpenFOAM/containers/Lists/UList/UList.H, you'll see that forAll is a very simple macro that cannot be having the sort of side-effects you report.
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 pseudo-code, 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;
    }
olesen is offline   Reply With Quote

Old   December 13, 2012, 16:14
Default
  #3
New Member
 
Yanxiang Shi
Join Date: Mar 2012
Location: Cambridge, MA, US
Posts: 29
Rep Power: 5
yanxiang is on a distinguished road
Did you get this problem solved? I am having exactly the same issue. The forAll loop does not update the variable at all.
yanxiang is offline   Reply With Quote

Old   December 13, 2012, 16:25
Default
  #4
Member
 
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 8
Pascal_doran is on a distinguished road
Hello Yanxiang,

That make a long time ago I think it was a simple error in my if statement:
Code:
if (diff <= 0.00002)
should have been:
Code:
if (abs(diff) <= 0.00002)
Pascal
Pascal_doran is offline   Reply With Quote

Old   December 13, 2012, 16:41
Default
  #5
New Member
 
Yanxiang Shi
Join Date: Mar 2012
Location: Cambridge, MA, US
Posts: 29
Rep Power: 5
yanxiang is on a distinguished road
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;
            }
        }
where beta is a scalar that should be bounded within 0 and 1. But after this snippet, the output of min(beta).value() still gives negative value. I was wondering what might be missing.

yanxiang
yanxiang is offline   Reply With Quote

Old   December 13, 2012, 16:55
Default
  #6
Member
 
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 8
Pascal_doran is on a distinguished road
Hi,
Your sample of code seems all right. Could you post the entire code? It will be easier for others to give help.

Pascal
Pascal_doran is offline   Reply With Quote

Old   December 13, 2012, 17:02
Default
  #7
New Member
 
Yanxiang Shi
Join Date: Mar 2012
Location: Cambridge, MA, US
Posts: 29
Rep Power: 5
yanxiang is on a distinguished road
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;
    }
}
yanxiang is offline   Reply With Quote

Old   December 13, 2012, 17:13
Default
  #8
Member
 
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 8
Pascal_doran is on a distinguished road
Hi!

I think that this line of code :
Code:
beta.correctBoundaryConditions();
is responsable for your negative betas. Try to output min(beta).value() before and after the correction of betas. To avoid this problem you can put the correction of beta before your forall loop.

Pascal
Pascal_doran is offline   Reply With Quote

Old   December 13, 2012, 17:25
Default
  #9
New Member
 
Yanxiang Shi
Join Date: Mar 2012
Location: Cambridge, MA, US
Posts: 29
Rep Power: 5
yanxiang is on a distinguished road
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
yanxiang is offline   Reply With Quote

Old   December 14, 2012, 16:04
Default
  #10
Senior Member
 
Kent Wardle
Join Date: Mar 2009
Location: Illinois, USA
Posts: 195
Rep Power: 10
kwardle is on a distinguished road
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 + ... (n-1)*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.
neeraj likes this.
kwardle is offline   Reply With Quote

Old   December 14, 2012, 18:39
Default
  #11
New Member
 
Yanxiang Shi
Join Date: Mar 2012
Location: Cambridge, MA, US
Posts: 29
Rep Power: 5
yanxiang is on a distinguished road
Quote:
beta = max(beta, 0.0);
Thanks for pointing this out. I wasn't aware of it and indeed I was looking for something like this. Anyhow, it worked like a charm. Nevertheless, I still couldn't understand why the loop doesn't give the results I wanted.

Quote:
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 + ... (n-1)*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.
You are perfectly right that alphas is in fact alpha_solid as opposed to the one defined in mpEF. And also note that I adapted the tPEF code instead of mpEF, and there the MULES is not used. But you are right. MULES ensures the boundedness. So I will try to add that to the modified tPEF to see if solution still blows up (pressure being extremely high). Will post the results once I get it working.

Thanks,
yanxiang
yanxiang is offline   Reply With Quote

Reply

Tags
forall loop

Thread Tools
Display Modes

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 Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Problem with Gmsh nishant_hull Open Source Meshers: Gmsh, Netgen, CGNS, ... 19 Today 05:10
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


All times are GMT -4. The time now is 18:26.