[Pascal_doran]

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

    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;
    }

[olesen]

Originally Posted by Pascal_doran

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),

    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;
    }

[yanxiang]

Did you get this problem solved? I am having exactly the same issue. The forAll loop does not update the variable at all.

[Pascal_doran]

Hello Yanxiang,

That make a long time ago I think it was a simple error in my if statement:

if (diff <= 0.00002)

should have been:

if (abs(diff) <= 0.00002)

Pascal

[yanxiang]

Thanks, Pascal. Now it looks like I have a different issue. So I have the following 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

[Pascal_doran]

Hi,
Your sample of code seems all right. Could you post the entire code? It will be easier for others to give help.

Pascal

[yanxiang]

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

{
    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;
    }
}

[Pascal_doran]

Hi!

I think that this line of 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

[yanxiang]

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

[kwardle]

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.

[yanxiang]

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.

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