[CFDanielGER]

Hello everybody,


for a new hybrid turbulence model I want to implement I need the dimensions of each cell to calculate a new delta. I found a way of calculating the dimensions in each direction (x, y, z) in meter. See code attached:

  const faceList & ff = this->mesh().faces();
  const pointField & pp = this->mesh().points();
  volScalarField deltax = this->y_;
  volScalarField deltay = this->y_;
  volScalarField deltaz = this->y_;

  forAll(this->mesh().C(), celli)
  {
    const cell & cc = this->mesh().cells()[celli];
    labelList pLabels(cc.labels(ff));
    pointField pLocal(pLabels.size(), vector::zero);

    forAll(pLabels, pointi)
    {
      pLocal[pointi] = pp[pLabels[pointi]];
    }

    dimensionedScalar Convert ("Convert", dimensionSet(0, 1, 0, 0, 0, 0, 0), scalar(1));

    deltax = (Foam::max(pLocal & vector(1,0,0)) - Foam::min(pLocal & vector(1,0,0))) * Convert;
    deltay = (Foam::max(pLocal & vector(0,1,0)) - Foam::min(pLocal & vector(0,1,0))) * Convert;
    deltaz = (Foam::max(pLocal & vector(0,0,1)) - Foam::min(pLocal & vector(0,0,1))) * Convert;
  }

I think the code is working fine, at least it is compiling and a few iterations are running but I have to check further.


But my question is:


Is there a more effective way of computing the cell dimensions in x, y and z than using a forAll loop? I have the impression that this is rather time-consuming when using a large mesh.



Many thanks for your help!


Daniel

[olesen]

boundBox for the cell points?

BTW in your code example it appears that you are overwriting the entire fields each cell - probably not what you want.

For efficiency, it would make sense to store/access the values as a DimensionedField of vectors and dump all of it into a class as a MeshObject. This would make it demand-driven and you could handle mesh changes seamlessly.

[CFDanielGER]

Maybe I show the whole code so it becomes a bit more clear. I should also mention that I am a beginner with OpenFOAM coding and found this code in a other forum post and it seemed that this the thing I need: Cell size (x,y,z)

My whole code for calculating the variable I need with the cell dimensions looks like this. For a non-moving mesh, the dimensions of each cell should only calculated at the beginning since they are not changing in time, but I am not sure yet how to solve this issue. Maybe write another file like the cubeRootVolDelta or the maxDelta. But what changes with each iteration is the vorticity, this should therefore be calculated at each time step.

template<class BasicTurbulenceModel>
tmp<volScalarField> SAEDES<BasicTurbulenceModel>::deltaOmega
(

) const
{
  volVectorField vorticity = fvc::curl(this->U());
  volScalarField magVorticity = mag(vorticity);
  volVectorField vorticityUnit = vorticity / magVorticity;

  const faceList & ff = this->mesh().faces();
  const pointField & pp = this->mesh().points();
  volScalarField deltax = this->y_;
  volScalarField deltay = this->y_;
  volScalarField deltaz = this->y_;

  forAll(this->mesh().C(), celli)
  {
    const cell & cc = this->mesh().cells()[celli];
    labelList pLabels(cc.labels(ff));
    pointField pLocal(pLabels.size(), vector::zero);

    forAll(pLabels, pointi)
    {
      pLocal[pointi] = pp[pLabels[pointi]];
    }

    dimensionedScalar Convert ("Convert", dimensionSet(0, 1, 0, 0, 0, 0, 0), scalar(1));

    deltax = (Foam::max(pLocal & vector(1,0,0)) - Foam::min(pLocal & vector(1,0,0))) * Convert;
    deltay = (Foam::max(pLocal & vector(0,1,0)) - Foam::min(pLocal & vector(0,1,0))) * Convert;
    deltaz = (Foam::max(pLocal & vector(0,0,1)) - Foam::min(pLocal & vector(0,0,1))) * Convert;

    //Info <<deltax<< endl;
    //Info <<deltay<< endl;
    //Info <<deltaz<< endl;
  }

  return sqrt( sqr(vorticityUnit.component(0)) * deltay * deltaz + sqr(vorticityUnit.component(1)) * deltaz * deltax + sqr(vorticityUnit.component(2)) * deltax * deltay );
}

Originally Posted by olesen

boundBox for the cell points?

BTW in your code example it appears that you are overwriting the entire fields each cell - probably not what you want.

For efficiency, it would make sense to store/access the values as a DimensionedField of vectors and dump all of it into a class as a MeshObject. This would make it demand-driven and you could handle mesh changes seamlessly.

Where do you mean the entire fields will be overwritten each cell, sorry if i don't understand yet clearly. And may you please explain or describe further how your idea of storing and accessing the values look like?

Many thanks for your help!

Best regards,
Daniel

[olesen]

Originally Posted by CFDanielGER

Where do you mean the entire fields will be overwritten each cell, sorry if i don't understand yet clearly. And may you please explain or describe further how your idea of storing and accessing the values look like?

If you examine the closing statements in your loop:

  forAll(this->mesh().C(), celli)
  {
...
// !!! HERE ! !!!

     deltax = (Foam::max(pLocal & vector(1,0,0)) - Foam::min(pLocal & vector(1,0,0))) * Convert;
    deltay = (Foam::max(pLocal & vector(0,1,0)) - Foam::min(pLocal & vector(0,1,0))) * Convert;
    deltaz = (Foam::max(pLocal & vector(0,0,1)) - Foam::min(pLocal & vector(0,0,1))) * Convert;
  }

You are assigning deltax, deltay etc, which are volume fields with a single value.
After exiting the loop, the entire field will simply have the value of the last cell, which you will not notice if the mesh is uniform.

[olesen]

You would need to verify, but my first idea was something like this (NOTE: quality not checked!!).

volScalarField deltax = this->y_;
volScalarField deltay = this->y_;
volScalarField deltaz = this->y_;

... starting here
const polyMesh& pMesh = this->mesh();
const label nCells = pMesh.nCells();


boundBox cellBb;
for (label celli = 0; celli < nCells; ++celli)
{
    const labelList& cPoints = pMesh.cellPoints(celli);
    cellBb.clear();
    cellBb.add(pMesh.points(), cPoints);

   deltax[celli] = (cellBb.max().x() - cellBb.min().x());
   deltay[celli] = (cellBb.max().y() - cellBb.min().y());
   deltaz[celli] = (cellBb.max().z() - cellBb.min().z());
}

[CFDanielGER]

Many thanks for your help, now I understand! You are completely right. I just tested your piece of code and it seems to work also it seems much faster, I just have to check if the dimensions are right. I think for this the easiest way is to create IOobjects, what do you think?

But if I now want to calculate the deltaOmega for each cell, I don't have to use a forAll loop, right? I can just calculate like this:

sqrt( sqr(vorticityUnit.component(0)) * deltay * deltaz + sqr(vorticityUnit.component(1)) * deltaz * deltax + sqr(vorticityUnit.component(2)) * deltax * deltay );

[olesen]

It's fine enough the way you now have it.
Since you have complete volfields after the loop, the rest should be OK.

You will Probably want to store deltaxy, deltayz, deltazx instead. Cheap to calculate within the loop, and saves extra multiplications later.

After that can revise how you create your fields. Probably want them with the class scope. For that, autoPtr is the usual container.

[saicharan662000@gmail.com]

Hi Daniel,

what does this->Y in your code mean??
Thanks in advance

[CFDanielGER]

Originally Posted by saicharan662000@gmail.com

Hi Daniel,

what does this->Y in your code mean??
Thanks in advance

That should be the wall distance if I remember correctly. It was just for initialization.

[saicharan662000@gmail.com]

I had written a loop similar to yours for finding deltaX and wrote return statement after the loop. But when I use deltaX in return statement it is throwing an error message that deltaX is undefined

[CFDanielGER]

Originally Posted by saicharan662000@gmail.com

I had written a loop similar to yours for finding deltaX and wrote return statement after the loop. But when I use deltaX in return statement it is throwing an error message that deltaX is undefined

Did you define deltaX before the loop?

[saicharan662000@gmail.com]

Hi Daniel,
I didn't define it before loop. Because I don't know what I should assign to deltaX as u assigned this->Y. my cose looks like this


Foam::tmp<Foam::volScalarField>
Foam:haseChangeTwoPhaseMixtures::myPhaseChange:: kappaGradT() const
{

const volScalarField limalpha1(min(max((alpha1()), scalar(0)), scalar(1)));
const volVectorField& U = db().lookupObject<volVectorField>("U");
const fvMesh & mesh = U.mesh();
const unallocLabelList& neighbour = mesh.neighbour();
const unallocLabelList& owner = mesh.owner();
const volVectorField& C = mesh.C();
const volScalarField& T = alpha1().db().lookupObject<volScalarField>("T");
forAll(owner, facei)
{
vector delta = C[neighbour[facei]] - C[owner[facei]];
volScalarField gradAlpha = fvc::grad(limalpha1) & mesh.C();
volScalarField alphaPrim = mag(gradAlpha);
const dimensionedScalar len("len",dimLength,mag(delta));
const dimensionedScalar small("small",dimensionSet(0, 0, 0, 0, 0, 0, 0),1e-5);

}
return
(
(limalpha1*lambda1_+ (1-limalpha1)*lambda2_) * ((T-Tsat_)/(len))*(gradAlpha/max(alphaPrim,small))
);
}

[CFDanielGER]

Originally Posted by saicharan662000@gmail.com

Hi Daniel,
I didn't define it before loop. Because I don't know what I should assign to deltaX as u assigned this->Y. my cose looks like this


Foam::tmp<Foam::volScalarField>
Foam:haseChangeTwoPhaseMixtures::myPhaseChange:: kappaGradT() const
{

const volScalarField limalpha1(min(max((alpha1()), scalar(0)), scalar(1)));
const volVectorField& U = db().lookupObject<volVectorField>("U");
const fvMesh & mesh = U.mesh();
const unallocLabelList& neighbour = mesh.neighbour();
const unallocLabelList& owner = mesh.owner();
const volVectorField& C = mesh.C();
const volScalarField& T = alpha1().db().lookupObject<volScalarField>("T");
forAll(owner, facei)
{
vector delta = C[neighbour[facei]] - C[owner[facei]];
volScalarField gradAlpha = fvc::grad(limalpha1) & mesh.C();
volScalarField alphaPrim = mag(gradAlpha);
const dimensionedScalar len("len",dimLength,mag(delta));
const dimensionedScalar small("small",dimensionSet(0, 0, 0, 0, 0, 0, 0),1e-5);

}
return
(
(limalpha1*lambda1_+ (1-limalpha1)*lambda2_) * ((T-Tsat_)/(len))*(gradAlpha/max(alphaPrim,small))
);
}

you could try to define a vector before the loop like:

vector delta(0,0,0);

[saicharan662000@gmail.com]

Hi Daniel,
Thanks for your suggestion,I implemented it like this now I have no eroors or warnings,is it correct way??



Foam::tmp<Foam::volScalarField>
Foam:haseChangeTwoPhaseMixtures::myPhaseChange:: kappaGradT() const
{

const volScalarField limalpha1(min(max((alpha1()), scalar(0)), scalar(1)));
const volVectorField& U = db().lookupObject<volVectorField>("U");
const fvMesh & mesh = U.mesh();
const unallocLabelList& neighbour = mesh.neighbour();
const unallocLabelList& owner = mesh.owner();
const volVectorField& C = mesh.C();
const volScalarField& T = alpha1().db().lookupObject<volScalarField>("T");
vector delta(0,0,0);
const dimensionedScalar len("len",dimLength,0);
forAll(owner, facei)
{
vector delta = C[neighbour[facei]] - C[owner[facei]];
volScalarField gradAlpha = fvc::grad(limalpha1) & mesh.C();
volScalarField alphaPrim = mag(gradAlpha);
const dimensionedScalar len("len",dimLength,mag(delta));
const dimensionedScalar small("small",dimensionSet(0, 0, 0, 0, 0, 0, 0),1e-5);

}
return
(
(limalpha1*lambda1_+ (1-limalpha1)*lambda2_) * ((T-Tsat_)/(len))*(gradAlpha/max(alphaPrim,small))
);
}


is delta in loop getting used? Why I am doubting this is, I have set both delta and len to 0 before the loop ,will the get updated

[CFDanielGER]

Originally Posted by saicharan662000@gmail.com

Hi Daniel,
Thanks for your suggestion,I implemented it like this now I have no eroors or warnings,is it correct way??



Foam::tmp<Foam::volScalarField>
Foam:haseChangeTwoPhaseMixtures::myPhaseChange:: kappaGradT() const
{

const volScalarField limalpha1(min(max((alpha1()), scalar(0)), scalar(1)));
const volVectorField& U = db().lookupObject<volVectorField>("U");
const fvMesh & mesh = U.mesh();
const unallocLabelList& neighbour = mesh.neighbour();
const unallocLabelList& owner = mesh.owner();
const volVectorField& C = mesh.C();
const volScalarField& T = alpha1().db().lookupObject<volScalarField>("T");
vector delta(0,0,0);
const dimensionedScalar len("len",dimLength,0);
forAll(owner, facei)
{
vector delta = C[neighbour[facei]] - C[owner[facei]];
volScalarField gradAlpha = fvc::grad(limalpha1) & mesh.C();
volScalarField alphaPrim = mag(gradAlpha);
const dimensionedScalar len("len",dimLength,mag(delta));
const dimensionedScalar small("small",dimensionSet(0, 0, 0, 0, 0, 0, 0),1e-5);

}
return
(
(limalpha1*lambda1_+ (1-limalpha1)*lambda2_) * ((T-Tsat_)/(len))*(gradAlpha/max(alphaPrim,small))
);
}


is delta in loop getting used? Why I am doubting this is, I have set both delta and len to 0 before the loop ,will the get updated

Yes, that should be working. You are just defining delta and len beforehand so that it will be recognized for the loop. During the loop it will then get updated. You could print len or delta during your simulation so you can see if the results are meaningful. I am not sure how to print anymore but there are some posts related to this topic, just google for it.

[saicharan662000@gmail.com]

Hi daniel,
Thanks for helping me, And I am able to print them and I am getting correct results

[CFDanielGER]

Originally Posted by saicharan662000@gmail.com

Hi daniel,
Thanks for helping me, And I am able to print them and I am getting correct results

Perfect, no worries!

[saicharan662000@gmail.com]

Hi Daniel,
Can you look this post and suggest me what's wrong in it.
Multiplication of source term with div(T) is causing error
Thanks in advance

[Wang Shang]

Hi Olesen,

If I want to write it as a post-processing program, is that possible?

Regards.