[tayo]

I want to include if-else loop in my interPhaseChangeFoam solver. This is supposed to be easy but it's taking some time now. Below is an illustration of the code.

Foam::volScalarField Foam::phaseChangeTwoPhaseMixture::Random::dotAlpha() const
{
 const volScalarField& p = alpha1_.db().lookupObject<volScalarField>("p");
 volScalarField x = this->x();
 forAll(p, celli)
 {
   if (p[celli] <= pSat_.value())
   {
     return volScalarField
     (
       x/y_*(p - pSat_)*(1 - alpha1_)
     );
   else
   {
     return volScalarField
     (
       x*(p - pSat_)*alpha1_
     );
   }
 }
}

pSat_ and y_ are dimensionedScalar.

Notice that I didn't index volScalarfields p and x in the if-else statement because I noticed that indexing makes it dimensionless afterwards. It however compiles the way it is now (ofcourse not correct) but gives the warning "control reaches end of non-void function". Please how do I fix this if-else loop? Thanks in advance.

[jherb]

Your function returns in the first iteration of your for loop. You can create a temporary field in which you store your values and return that field (or a reference to it). (Perhaps also this bug report might be of interest what can go wrong http://www.openfoam.org/mantisbt/view.php?id=912 )

[tayo]

Originally Posted by jherb

Your function returns in the first iteration of your for loop. You can create a temporary field in which you store your values and return that field (or a reference to it). (Perhaps also this bug report might be of interest what can go wrong http://www.openfoam.org/mantisbt/view.php?id=912 )

Hi Jherb, I don't really get what you mean, Please explain with an example. Thank you

[ngj]

Also, since you do not modify x, I suggest that you make it into a const reference. The fact that you are copying x will also include some unnecessary overhead.

As already stated, you loop over all cells, but it is the value of the pressure in the first cell, which completely determines the outcome of your method. You should revise the location of your return statements.

Kind regards

Niels

[tayo]

Thanks Niels & Joachim, I made the corrections to the loop and return statement as advised. The code now compiles.

Foam::volScalarField Foam::phaseChangeTwoPhaseMixture::Random::dotAlpha() const
{
 const volScalarField& p = alpha1_.db().lookupObject<volScalarField>("p");
 volScalarField x = this->x();
 forAll(p, celli)
 {
   if (p[celli] <= pSat_.value())
   {
       dotAlpha() = x/y_*(p[celli] - pSat_.value())*(1 - alpha1_);
   }
   else
   {
       dotAlpha() = x/y_*(p[celli] - pSat_.value())*alpha1_;
   }
 }
 return dotAlpha();
}

However, due to the indexing of pressure and pSat_, the dimension of the returned dotAlpha() becomes inconsistent. i.e. (p[celli] - pSat_.value()) in the equation becomes dimensionless. Because of this, it is not able to fully run. How do I ensure that dotAlpha() keeps the intended dimension after indexing p? Thanks

[ngj]

But you still change the entire field based on a single if statement on the pressure, so now it looks like the pressure in the final cell in the computational domain controls the return value.

I am, however, not certain, because there is some weird recursive nature going on with calls to dotAlpha() inside its own function, so in worst case it will scale with number of elements squared, do an immense amount of copying of memory and still return the wrong result.

As an easy starting point, you must initialise the return variable, populate the data in each element corresponding to each cell and then return the field. To initialise a vol-field, you could look in createFields.H in your favourite solver.

Kind regards

Niels

[tayo]

Hi Niels, you're right. My earlier loop (though compiled0 was recursive and thus cannot run beyond the continuity equation in the first iteration. So to initialize dotAlpha, I created a volScalarField (mAlpha) in createfield.H as shown below.

volScalarField mAlpha
(
  IOobject
  (
    "mAlpha",
    runTime.time(),
    mesh,
    IOobject::NO_READ,
    IOobject::NO_WRITE
  ),
  mesh,
  dimensionedScalar("mAlpha",dimensionSet(1,-3,-1,0,0,0,0),0)
);

mAlpha is then computed in the loop below and returned so that it will be equal to dotAlpha(). However, when I compile this way, it indicates that "mAlpha is not declared in this scope". I don't know how to call mAlpha (from createfield.H) into the loop.

Foam::volScalarField Foam::phaseChangeTwoPhaseMixture::Random::dotAlpha() const
{
 const volScalarField& p = alpha1_.db().lookupObject<volScalarField>("p");
 volScalarField x = this->x();
 forAll(p, celli)
 {
   if (p[celli] <= pSat_.value())
   {
       mAlpha = x/y_*(p[celli] - pSat_.value())*(1 - alpha1_);
   }
   else
   {
       mAlpha = x/y_*(p[celli] - pSat_.value())*alpha1_;
   }
 }
 return mAlpha;
}

Attempts to call mAlpha in the loop has only led back to the recursive iteration. Please how do correct this? Also, kindly correct if you still notice any error in the way I'm looping. Thanks

[ngj]

Well, when you define the variable in createFields.H, it is out of scope (aka unknown) for your function. I merely referred to createFields.H for your inspiration.

The following should work, no optimally, but it will do what you need. Please note that I have not tried compiling it:

Foam::volScalarField foam::phaseChangeTwoPhaseMixture::Random::dotAlpha() const
{
  volScalarField mAlpha
  (
    IOobject
    (
       "mAlpha",
       alpha_.mesh().time().timeName(),
       alpha_.mesh(),
       IOobject::NO_READ,
       IOobject::NO_WRITE
    ),
    alpha_.mesh(),
    dimensionedScalar("null", dimensionSet(1,-3,1,0,0,0,0),0),
    "zeroGradient"
  );
 
  const volScalarField& p = alpha_.db().lookupObject<volScalarField>("p");
  const volScalarField& x = this->x();

  forAll (p, celli)
  {
    if (p[celli] <= pSat_.value())
    {
        mAlpha[celli] = x[celli]/y_[celli]*(p[celli] - pSat_.value())*(1 - alpha1_[celli]);
    }
    else
    {
        mAlpha[celli] = x[celli]/y_[celli]*(p[celli] - pSat_.value())*alpha_[celli];
    }
  }
  return mAlpha;
}

[tayo]

Hi Niels, it worked after making some adjustments. Thank you .