CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   Write Cell Volumes and Velocities for each TimeStep (https://www.cfd-online.com/Forums/openfoam-programming-development/186436-write-cell-volumes-velocities-each-timestep.html)

sisetrun April 19, 2017 04:35
Write Cell Volumes and Velocities for each TimeStep
 
Hello everybody,

I am write a little post-processing tool which writes out the cell volumes and the velocities in the cells for a specific region in the domain.

Since I am only interested in a small part of my domain, I define a distance from z_min to z_max in an so called CoordinatesDict in the constant directory. I want to write the data in a .dat file to proceed working with python.

At the moment the code looks like this...


Code:
// This file is based on OpenFOAM.
// Description
//    Write mesh cell volumes in volumes.dat file and velocity in
//    velocity.dat file in current folder.

#include "OFstream.H"
#include "argList.H"
#include "timeSelector.H"
#include "Time.H"
#include "fvMesh.H"
#include "vectorIOField.H"
#include "volFields.H"
#include "surfaceFields.H"


//#include "fvCFD.H"

using namespace Foam;

int main(int argc, char *argv[])
{
    timeSelector::addOptions();

#  include "setRootCase.H"
   
    Info << "\n***EXTRACT CELL VOLUMES AND VELOCITIES***\n" << endl;
    Info << "This applications is extacting the volumes and the velocities from a region in the domain\n" << endl;
   
#  include "createTime.H"
    instantList timeDirs = timeSelector::select0(runTime, args);

#  include "createMesh.H"
#  include "createFields.H"


    Info << "1) Reading all volumes..." << endl;
    Info << "Extracting cells between z = " << z_min.value();
    Info << " and z = " << z_max.value() << endl;
   
    OFstream volumes("volumes_" + runTime.timeName() + ".dat");   

    forAll(mesh.C(), idx)
      {

        if (mesh.C()[idx].component(2) >= z_min.value())
            if (mesh.C()[idx].component(2) <= z_max.value())
                {
                    volumes << idx << ' ' << mesh.V()[idx] << ' ' << mesh.C()[idx] << endl;
                }
      }
     
    Info << "Volumes DONE...\n" << endl;

    Info << "2) Reading velocities...\n" << endl;

    forAll(timeDirs, timeI)
    {

        runTime.setTime(timeDirs[timeI], timeI);
        Info<< "Time = " << runTime.timeName() << endl;

        #include "createMesh.H"
        #include "createFields.H"

        OFstream velocity("velocity_" + runTime.timeName() + ".dat");

        Info << "Extracting cells between z = " << z_min.value();
        Info << " and z = " << z_max.value() << endl;

        //Extract the z component of the cell centers
        Info << "Writing data to the files velocity_" + runTime.timeName() + ".dat" << "\n" << endl;

        forAll(mesh.C(), idx)
        {

          if (mesh.C()[idx].component(2) >= z_min.value())

              if (mesh.C()[idx].component(2) <= z_max.value())
              {
                  velocity << idx << ' ' << U[idx].component(0) << ' ' <<  U[idx].component(1) << ' ' << U[idx].component(2) << endl;
              }
        }

    }

    Info << "Writing done \n" << endl;

    return 0;
}
My question: is there a way to skip the

Code:
#include "createMesh.H"
#include "createFields.H"
in the forAll(timeDirs, timeI)? When I leave out this statement, the program only writes out the velocities from the first time directory. Since I work with quite large meshes, reading the mesh for every time step is very time consuming.

I tried so come around with mesh.readUpdate() but I did not work out...

Thank a lot for your input!

Best regards

WildeCat April 20, 2017 05:51
Default Write Cell Volumes and Velocities for each TimeStep
 
Hi Sisetrum

I do not know a way to get around reading

Code:
#include "createMesh.H"
#include "createFields.H"
but what you can do is rewrite to only read the last time instance by removing the forAll and modfying the runTime as shown below. I have not tested this solution, but have written post processing tools in a similar way in the past.

Code:
runTime.setTime(timeDirs[timeDirs.size()-1], timeDirs.size()-1);       
Info<< "Time = " << runTime.timeName() << endl;
This should reduce your computing time for only the latest time step, if that suffices you.

Hope I could help

sisetrun April 20, 2017 06:22
Hey,

thank you for your reply!

I had something similar to your version. As you said, this gets the data of the last time step...

What I want to get, is a tool that it works in the same way as the sample utility:
Read the mesh only one time to write out the volumes of the cells and extract the velocities from the data in the time directories.

Code:
SAMPLE:


Create time

---> CREATE THE MESH
Create mesh for time = 0

Reading set description:
    lineX1

Reading surface description:
    constantPlane

---> SAMPLE ON THE MESH FOR DIFFERENT TIME STEPS

Time = 0

Time = 0.25

Time = 0.37

End

WildeCat April 20, 2017 06:32
Links to sample
 
I recommend you look at the sample utility then.

https://cpp.openfoam.org/v3/a04479_source.html

and check how it loads the mesh in "#include "createNamedMesh.H"" before forAll.

sisetrun April 20, 2017 06:33
Thanks, I will give it a try...

WildeCat April 20, 2017 06:55
Sample code
 
Hi Sisetrun,

I have taken the liberty to modify your code to import the mesh only once. I think this should work. An idea to speed up your code is not not use "createFields.H" but write in only the variables you need specifically. I have not figured out if createFields or the equivalent variables can be written out of the time loop. To write only the variables you need look at other post processing utilities.

I Hope this helps.

Code:
// This file is based on OpenFOAM.
// Description
//    Write mesh cell volumes in volumes.dat file and velocity in
//    velocity.dat file in current folder.

#include "OFstream.H"
#include "argList.H"
#include "timeSelector.H"
#include "Time.H"
#include "fvMesh.H"
#include "vectorIOField.H"
#include "volFields.H"
#include "surfaceFields.H"


//#include "fvCFD.H"

using namespace Foam;

int main(int argc, char *argv[])
{
    timeSelector::addOptions();

#  include "setRootCase.H"
   
    Info << "\n***EXTRACT CELL VOLUMES AND VELOCITIES***\n" << endl;
    Info << "This applications is extacting the volumes and the velocities from a region in the domain\n" << endl;
   
#  include "createTime.H"
    instantList timeDirs = timeSelector::select0(runTime, args);
    // modification end
    #include "createMesh.H"

    Info << "1) Reading all volumes..." << endl;
    Info << "Extracting cells between z = " << z_min.value();
    Info << " and z = " << z_max.value() << endl;
   
    OFstream volumes("volumes_" + runTime.timeName() + ".dat");   

    forAll(mesh.C(), idx)
      {

        if (mesh.C()[idx].component(2) >= z_min.value())
            if (mesh.C()[idx].component(2) <= z_max.value())
                {
                    volumes << idx << ' ' << mesh.V()[idx] << ' ' << mesh.C()[idx] << endl;
                }
      }
     
    Info << "Volumes DONE...\n" << endl;

    Info << "2) Reading velocities...\n" << endl;

    forAll(timeDirs, timeI)
    {

        runTime.setTime(timeDirs[timeI], timeI);
        Info<< "Time = " << runTime.timeName() << endl;

    mesh.readUpdate();

        #include "createFields.H" //eventually write your own createFields.H

        OFstream velocity("velocity_" + runTime.timeName() + ".dat");

        Info << "Extracting cells between z = " << z_min.value();
        Info << " and z = " << z_max.value() << endl;

        //Extract the z component of the cell centers
        Info << "Writing data to the files velocity_" + runTime.timeName() + ".dat" << "\n" << endl;

        forAll(mesh.C(), idx)
        {

          if (mesh.C()[idx].component(2) >= z_min.value())

              if (mesh.C()[idx].component(2) <= z_max.value())
              {
                  velocity << idx << ' ' << U[idx].component(0) << ' ' <<  U[idx].component(1) << ' ' << U[idx].component(2) << endl;
              }
        }

    }

    Info << "Writing done \n" << endl;

    return 0;
}

sisetrun April 20, 2017 07:29
Awesome... Thanks a lot for your help!


All times are GMT -4. The time now is 19:30.