write cell volumes
Hi there - I am trying to modify the writeCellCentres application to get something that will write cell volumes. The heart of the code is as follows:
// Main program: int main(int argc, char *argv[]) { timeSelector::addOptions(); # include "setRootCase.H" # include "createTime.H" instantList timeDirs = timeSelector::select0(runTime, args); # include "createMesh.H" forAll(timeDirs, timeI) { runTime.setTime(timeDirs[timeI], timeI); Info<< "Time = " << runTime.timeName() << endl; // Check for new mesh mesh.readUpdate(); scalarField cv ( IOobject ( "cellVolumes", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), mesh.V() ); Info<< "Writing cellCentre positions to " << cv.name() << " in " << runTime.timeName() << endl; cv.write(); This code mainly differs from the writeCellCentres in that I replaced mesh.C() with mesh.V(). Also, where I define scalarField cv, the original code used volVectorField cc. mesh.V() is used in a variety of places in openFOAM, but apparently it is not a scalarField or a volScalarField because I get a variety of compilation errors depending on which type I specify. Does anyone know a good way to address this problem? Many Thanks. |
Hi Daniel,
you can try this: Code:
volScalarField cv |
Hi Martin - That did the trick. Thank you very much.
|
Hi,
i have the same problem. I want to know the volume of each cell and let the program prints them somewhere. With this trick where openFoam prints the volume? Thank you |
Hi Andrea,
you can loop over the cells volumes and print them with Code:
forAll(mesh.V(),celli) |
Hi Martin
thanks for the reply. Where i have to put these lines? Andrea |
1 Attachment(s)
Hi Andrea - Just in case you want it I have attached the source file that I use to get cell Volumes (Thanks to Martin:)). I put the file in a directory I created called writeCellVolumes that is in the same directory as writeCellCentres. I also copied over the the other make files (Make/files and Make/options) from writeCellCentres then update those files with the correct EXE = path names then wmake it. The application is executed within a case directory. Good Luck
|
Hi Andrea,
it depends ;-) Without knowing your solver and the operation you want to perform with cell volumes I can only guess: - you might want to create a volScalarField that will contain the results - you loop over the mesh as described above and calculate something important, which is then stored in your new volScalarField - after the loop you let your volScalarField be written out to the current time directory This action can be done after the solver is finished completely (saves time), or during each writing of time step, or even in a separated utility that can be run independently of the solver. You have to make sure that all fields you need for the computation are valid, since some solver use "{ ... }" to end field's life early and reduce memory usage. May be you can post or send your solver and your intended calculation... Martin |
Hi Martin
I just started using OpenFoam as part of my phd and so i am a very new user. At the moment i would like to simulate multiphase flows using interFoam, in simple geometry created with blockMesh. The geometry at the moment is a sphere so the cell volumes are not equal. I just want the openFoam print volumes of each cell in a file, like "points" or "owner" in the polyMesh directory, named for example "Volumes". Otherwise that openfoam prints after the simulation these volumes during post-processing, like "writeCellCentres" (something like writeCellVolumes). I dont know which is the best way to do that. Thanks |
Hi Andrea,
if you don't use a dynamic mesh (i.e. cell volumes stay constant) you can write the values immediately before "Info<< "\nStarting time loop\n" << endl;" in interFoam. To write the cell volumes to "constant/polyMesh": Code:
volScalarField cv http://openfoamwiki.net/index.php/In...IOobject_class Martin |
Thanks a lot for your help.
After adding this line i guess that you have to recompile the program...or not? andrea |
Yes, you must recompile. Simply call wmake in the interFoam directory to make a quick shot.
By the way it might be a good idea to make these kind of changes in your user application folder. Copy and rename interFoam to your favorite project name and do changes there... here is a good explanation how to do it: http://openfoamwiki.net/index.php/Ho...ure_to_icoFoam Martin |
Great...it runs...thanks a lot again!
andrea |
Hi,
i got another problem. Now I'm trying to import an external geometry using gmsh (so without using BlockMesh). I would like to do the same thing and print the volume of each cell in the constant/polymesh subdirectory. (while importing the geometry for example). Is possible? Where i have to put those line now? Thanks Andrea |
Hi Andrea,
you could append the lines from post #10 to each mesh converting tool, of course... for gmshToFoam.C it would be directly before the line "Info<< "End\n" << endl;". However it would make more sense to create a separate tool like Daniels "writeCellVolumes.C" from post #7. Replace the location where to write in the IOobject part. In your every day work you call the mesh converter for gmsh, followed by a writeCellVolumes call. Martin |
Thanks Martin
Your advice were very helpful!! Andrea |
I continue this discussion because the new problem is similar to the old.
I am using the sample Dictionary to get the value of the color function alpha1 on a particular patch (wall for my problem). Is possible (similarly to obtain the volume of each cell) to let OF prints not only the value of the color function but also the area of the surface (in m^2) for each cell on the wall patch?? Thanks Andrea |
Hi Andrea,
you might have a look at this thread: http://www.cfd-online.com/Forums/ope...ergence-2.html In post #29 there is a piece of code that calculates the wall shear stress and the surface area for each cell surface on a given patch. Replace the wall shear stress by your alpha1 and you have your code... installation hints are given in this thread, too. A parallel version is given in post #44. I recommend, that you try to implement a uniprocessor version that you can call like the writeCellVolumes tool discussed earlier. The patch name of the wall you want to examine can be passed as a command line argument or via dictionary. Martin |
Thanks a lot Martin for all your help!!
I only need the values of cell's area on the patch and i can get alpha1 on the same patch directly from the sample. So i changed your code (for parallel run, because i really need parallel) in this way. If I have not made mistakes should calculate the cell area on the patch "obstacle " and write it in area.dat. Please check if there are any errors (sure there will be...) I commented out all lines that I did not need.. #include "RASModel.H" // write matlab file { Info << "Writing Area of the cells on the patch!" << endl; word wallPatchName = "obstacle"; Info << "Searching for patch " << wallPatchName << endl; // Check, if the patch name exists: label patchID = mesh.boundaryMesh().findPatchID(wallPatchName); if (patchID < 0) { Info << "No patch found with name \"" << wallPatchName << "\"!" << endl; } // else // { // calculating tau //# include "createFields2.H" // mesh.readUpdate(); // volSymmTensorField Reff(RASModel->devReff()); // volVectorField tau // ( // IOobject // ( // "tau", // runTime.timeName(), // mesh, // IOobject::NO_READ, // IOobject::AUTO_WRITE // ), // mesh, // dimensionedVector // ( // "tau", // Reff.dimensions(), // vector::zero // ) // ); // forAll(tau.boundaryField(), patchi) // { // tau.boundaryField()[patchi] = // ( // -mesh.Sf().boundaryField()[patchi] // /mesh.magSf().boundaryField()[patchi] // ) & Reff.boundaryField()[patchi]; // } // Pointer to the patch we want to have shear stress for: const polyPatch& cPatch = mesh.boundaryMesh()[patchID]; label nFaces = cPatch.size(); // pointField faceCenterCoords(nFaces, vector(0,0,0)); pointField faceAreas(nFaces, vector(0,0,0)); // pointField wallShearStress(nFaces, vector(0,0,0)); label runFaces = 0; // collecting the interesting data for (int i=0; i<nFaces; i++) { // faceCenterCoords[runFaces] = mesh.Cf().boundaryField()[patchID][i]; faceAreas[runFaces] = mesh.Sf().boundaryField()[patchID][i]; // wallShearStress[runFaces] = tau.boundaryField()[patchID][i]; runFaces++; } // for parallel reduce(nFaces, sumOp<label>()); Info << "Number of global faces: " << nFaces << endl; // make a cute filename fileName casePath = cwd();//casePath.name(); fileName caseName = casePath.name(); fileName myFile = "area.dat"; Info << "Data is written to file " << myFile << endl; // create a new file OFstream *myStream = NULL; if (Pstream::master()) { myStream = new OFstream(myFile); } // slaves send data to master, master writes data to file if (Pstream::parRun()) { if (Pstream::myProcNo() != Pstream::masterNo()) { OPstream toMaster(Pstream::blocking, Pstream::masterNo()); // toMaster << faceCenterCoords; toMaster << faceAreas; // toMaster << wallShearStress; } else // master part { // first write own data for (int i=0; i<faceCenterCoords.size(); i++) { *myStream << i << ";" // << faceCenterCoords[i].component(0) << "; " // << faceCenterCoords[i].component(1) << "; " // << faceCenterCoords[i].component(2) << "; " // << wallShearStress[i].component(0) << "; " // << wallShearStress[i].component(1) << "; " // << wallShearStress[i].component(2) << "; " // << mag(wallShearStress[i]) << endl; << faceAreas[i] << "; " } // label run = faceCenterCoords.size(); // then receive slave data and write // pointField bufferFaceCenterCoords; pointField bufferFaceAreas; // pointField bufferWallShearStress; for (int slave=Pstream::firstSlave(); slave <= Pstream::lastSlave(); slave++) { IPstream fromSlave(Pstream::blocking, slave); // fromSlave >> bufferFaceCenterCoords; fromSlave >> bufferFaceAreas; // fromSlave >> bufferWallShearStress; for (int i=0; i<bufferFaceCenterCoords.size(); i++) { *myStream << run << ";" // << bufferFaceCenterCoords[i].component(0) << "; " // << bufferFaceCenterCoords[i].component(1) << "; " // << bufferFaceCenterCoords[i].component(2) << "; " // // << bufferWallShearStress[i].component(0) << "; " // << bufferWallShearStress[i].component(1) << "; " // << bufferWallShearStress[i].component(2) << "; " // << mag(bufferWallShearStress[i]) << endl; << bufferFaceAreas << "; " run++; } } } } else // single processor only { for (int i=0; i<nFaces; i++) { *myStream << i << "; " // coordinates on face centers: // << mesh.Cf().boundaryField()[patchID][i].component(0) << "; " // << mesh.Cf().boundaryField()[patchID][i].component(1) << "; " // << mesh.Cf().boundaryField()[patchID][i].component(2) << "; " // wall shear stress as a vector: // << tau.boundaryField()[patchID][i].component(0) << "; " // << tau.boundaryField()[patchID][i].component(1) << "; " // << tau.boundaryField()[patchID][i].component(2) << "; " // magnitude of wall shear stress: // << mag(tau.boundaryField()[patchID][i]) << endl; << mesh.Sf().boundaryField()[patchID][i] << "; " } } } } Thanks andrea |
Hi Andrea,
some minor problems removed: Code:
// write matlab file Martin |
Hi again,
I found the bug: a bracket related problem... so here comes the (hopefully) correct parallel version: Code:
#include "RASModel.H" Martin |
Hi Martin,
after wclean and wmake i got this error (i dont know if the instructions that you had given to compile simpleFoam are the same for interFoam (the solver that i use) but seems that there is some error on #include RASModel.H) Making dependency list for source file interFoam.C could not open file RASModel.H for source file interFoam.C SOURCE=interFoam.C ; g++ -m64 -Dlinux64 -DWM_DP -Wall -Wno-strict-aliasing -Wextra -Wno-unused-parameter -Wold-style-cast -Wnon-virtual-dtor -O3 -DNoRepository -ftemplate-depth-40 -I/opt/OpenFOAM/OpenFOAM-1.7.0/src/transportModels -I/opt/OpenFOAM/OpenFOAM-1.7.0/src/transportModels/incompressible/lnInclude -I/opt/OpenFOAM/OpenFOAM-1.7.0/src/transportModels/interfaceProperties/lnInclude -I/opt/OpenFOAM/OpenFOAM-1.7.0/src/turbulenceModels/incompressible/turbulenceModel -I/opt/OpenFOAM/OpenFOAM-1.7.0/src/finiteVolume/lnInclude -IlnInclude -I. -I/opt/OpenFOAM/OpenFOAM-1.7.0/src/OpenFOAM/lnInclude -I/opt/OpenFOAM/OpenFOAM-1.7.0/src/OSspecific/POSIX/lnInclude -fPIC -c $SOURCE -o Make/linux64GccDPOpt/interFoam.o In file included from interFoam.C:125: write.H:1:22: error: RASModel.H: No such file or directory /opt/OpenFOAM/OpenFOAM-1.7.0/src/finiteVolume/lnInclude/readPISOControls.H: In function â: /opt/OpenFOAM/OpenFOAM-1.7.0/src/finiteVolume/lnInclude/readPISOControls.H:11: warning: unused variable â /opt/OpenFOAM/OpenFOAM-1.7.0/src/finiteVolume/lnInclude/readPISOControls.H:14: warning: unused variable â /opt/OpenFOAM/OpenFOAM-1.7.0/src/finiteVolume/lnInclude/readPISOControls.H:3: warning: unused variable â /opt/OpenFOAM/OpenFOAM-1.7.0/src/finiteVolume/lnInclude/readPISOControls.H:8: warning: unused variable â /opt/OpenFOAM/OpenFOAM-1.7.0/src/finiteVolume/lnInclude/readPISOControls.H:11: warning: unused variable â /opt/OpenFOAM/OpenFOAM-1.7.0/src/finiteVolume/lnInclude/readPISOControls.H:14: warning: unused variable â make: *** [Make/linux64GccDPOpt/interFoam.o] Error 1 Thanks again Andrea |
Hi Andrea,
just remove the line "#include RASModel.H". It was necessary for wall shear stress calculation in the original piece of code... Martin |
I figured that was enough to remove that line. I just launched a test. I write to you in case there are errors. Thank you very much for all your help.
Andrea |
Hi Martin, everything works fine!
I am taking advantage of your patience to ask another thing. I would like to calculate the gradient of alpha1 for each cell and each time step. I tried to write this piece of code: writeGradientAlpha1.C (writeCellVolumes.C them). Please look if is ok. I do not know if i have included the right source file (.H). I guess the function to calculate gradient is fvcGrad::(field). Code:
#include "argList.H" andrea |
1 Attachment(s)
Hi Andrea,
you must read the alpha1 field from disk first, otherwise your tool doesn't know what to work on. Extract the attached tool to your user directory, for example .../OpenFOAM/andrea-1.7.1/ and call wmake in the new folder. Have fun Martin |
Hi,
can i put this code in the /home/../applications/utilities/postProcessing/miscellaneous/ like writeCellCenter? for the next help I guess I will have to pay!! Thanks andrea |
Yes, you can.
Please have a look at the file Make/options Right now the executable is written to $(FOAM_USER_APPBIN) so only you have access to it. If you want to make it available globally (for other users), you may want to change this variable to $(FOAM_APPBIN) Martin |
Yes, i know that.
Thanks a lot andrea |
Hi Martin,
I have copied all file in the right directory and then wmake. When i call writeGradientAlpha1 from my case directory i got this error: | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ Build : 1.7.0-21131bcbd876 Exec : writeGradientAlpha1 Date : Feb 14 2011 Time : 12:28:14 Host : master.cluster PID : 25375 Case : /home/aferrari/Simulation/OstacoliRandom/simm2e-2alpha30/simmetricou2e-2 nProcs : 1 SigFpe : Enabling floating point exception trapping (FOAM_SIGFPE). // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Create time Create mesh for time = 0 Time = 0 Reading field alpha1 --> FOAM FATAL IO ERROR: Unknown patchField type constantAlphaContactAngle for patch type wall Valid patchField types are : 41 ( advective buoyantPressure calculated cyclic directMapped directionMixed empty fan fixedFluxPressure fixedGradient fixedInternalValue fixedPressureCompressibleDensity fixedValue freestream freestreamPressure inletOutlet inletOutletTotalTemperature mixed oscillatingFixedValue outletInlet partialSlip processor rotatingTotalPressure sliced slip symmetryPlane syringePressure timeVaryingMappedFixedValue timeVaryingMappedTotalPressure timeVaryingTotalPressure timeVaryingUniformFixedValue timeVaryingUniformInletOutlet totalPressure totalTemperature turbulentInlet turbulentIntensityKineticEnergyInlet uniformDensityHydrostaticPressure uniformFixedValue waveTransmissive wedge zeroGradient ) file: /home/aferrari/Simulation/OstacoliRandom/simm2e-2alpha30/simmetricou2e-2/0/alpha1::boundaryField::wall.4 from line 263883 to line 263885. From function fvPatchField<Type>::New(const fvPatch&, const DimensionedField<Type, volMesh>&, const dictionary&) in file /opt/OpenFOAM/OpenFOAM-1.7.0/src/finiteVolume/lnInclude/newFvPatchField.C at line 110. I use constantAlphaContactAngle patch instead of simply wall. Thanks andrea |
Hi Andrea,
to use the writeGradientAlpha1 tool in OpenFOAM-1.7.0 replace the lines in writeGradientAlpha1/Make/options with: Code:
EXE_INC = \ wclean wmake Martin |
Thanks Martin,
the correct library for the interface properties is : -linterfaceProperties \ . I had an error with -ltwoPhaseInterfaceProperties \. Thanks again andrea |
1 Attachment(s)
Hi Martin
I tried to edit your file to calculate the curvature of the interface, according to the calculations found in interfaceProperties.C Attached my writeCurvature.C. i got this error SOURCE=writeCurvature.C ; g++ -m64 -Dlinux64 -DWM_DP -Wall -Wno-strict-aliasing -Wextra -Wno-unused-parameter -Wold-style-cast -Wnon-virtual-dtor -O3 -DNoRepository -ftemplate-depth-40 -I/opt/OpenFOAM/OpenFOAM-1.7.0/src/transportModels -I/opt/OpenFOAM/OpenFOAM-1.7.0/src/transportModels/incompressible/lnInclude -I/opt/OpenFOAM/OpenFOAM-1.7.0/src/transportModels/interfaceProperties/lnInclude -I/opt/OpenFOAM/OpenFOAM-1.7.0/src/turbulenceModels/incompressible/turbulenceModel -I/opt/OpenFOAM/OpenFOAM-1.7.0/src/finiteVolume/lnInclude -IlnInclude -I. -I/opt/OpenFOAM/OpenFOAM-1.7.0/src/OpenFOAM/lnInclude -I/opt/OpenFOAM/OpenFOAM-1.7.0/src/OSspecific/POSIX/lnInclude -fPIC -c $SOURCE -o Make/linux64GccDPOpt/writeCurvature.o In file included from /opt/OpenFOAM/OpenFOAM-1.7.0/src/OpenFOAM/lnInclude/GeometricField.C:1273, from /opt/OpenFOAM/OpenFOAM-1.7.0/src/OpenFOAM/lnInclude/GeometricField.H:586, from /opt/OpenFOAM/OpenFOAM-1.7.0/src/OpenFOAM/lnInclude/GeometricScalarField.H:38, from /opt/OpenFOAM/OpenFOAM-1.7.0/src/OpenFOAM/lnInclude/GeometricFields.H:34, from /opt/OpenFOAM/OpenFOAM-1.7.0/src/finiteVolume/lnInclude/volFields.H:37, from writeCurvature.C:35: /opt/OpenFOAM/OpenFOAM-1.7.0/src/OpenFOAM/lnInclude/GeometricFieldFunctions.C: In function â: /opt/OpenFOAM/OpenFOAM-1.7.0/src/OpenFOAM/lnInclude/GeometricFieldFunctions.C:950: instantiated from â writeCurvature.C:98: instantiated from here /opt/OpenFOAM/OpenFOAM-1.7.0/src/OpenFOAM/lnInclude/GeometricFieldFunctions.C:950: error: no matching function for call to â /opt/OpenFOAM/OpenFOAM-1.7.0/src/OpenFOAM/lnInclude/scalarField.H:77: note: candidates are: void Foam::add(Foam::Field<double>&, const Foam::scalar&, const Foam::UList<double>&) /opt/OpenFOAM/OpenFOAM-1.7.0/src/OpenFOAM/lnInclude/scalarField.H:77: note: void Foam::add(Foam::Field<double>&, const Foam::UList<double>&, const Foam::scalar&) /opt/OpenFOAM/OpenFOAM-1.7.0/src/OpenFOAM/lnInclude/labelField.H:53: note: void Foam::add(Foam::Field<int>&, const Foam::label&, const Foam::UList<int>&) /opt/OpenFOAM/OpenFOAM-1.7.0/src/OpenFOAM/lnInclude/labelField.H:53: note: void Foam::add(Foam::Field<int>&, const Foam::UList<int>&, const Foam::label&) /opt/OpenFOAM/OpenFOAM-1.7.0/src/OpenFOAM/lnInclude/GeometricFieldFunctions.C:950: instantiated from â writeCurvature.C:98: instantiated from here /opt/OpenFOAM/OpenFOAM-1.7.0/src/OpenFOAM/lnInclude/GeometricFieldFunctions.C:950: error: no matching function for call to â /opt/OpenFOAM/OpenFOAM-1.7.0/src/OpenFOAM/lnInclude/scalarField.H:77: note: candidates are: void Foam::add(Foam::Field<double>&, const Foam::scalar&, const Foam::UList<double>&) /opt/OpenFOAM/OpenFOAM-1.7.0/src/OpenFOAM/lnInclude/scalarField.H:77: note: void Foam::add(Foam::Field<double>&, const Foam::UList<double>&, const Foam::scalar&) /opt/OpenFOAM/OpenFOAM-1.7.0/src/OpenFOAM/lnInclude/labelField.H:53: note: void Foam::add(Foam::Field<int>&, const Foam::label&, const Foam::UList<int>&) /opt/OpenFOAM/OpenFOAM-1.7.0/src/OpenFOAM/lnInclude/labelField.H:53: note: void Foam::add(Foam::Field<int>&, const Foam::UList<int>&, const Foam::label&) make: *** [Make/linux64GccDPOpt/writeCurvature.o] Error 1 what the problem for you? andrea |
1 Attachment(s)
Hi Andrea,
without doing testing, just removed some typos... have a look... Martin |
It runs!
do you know what OF calculate in these lines? const surfaceVectorField& Sf = mesh.Sf(); surfaceScalarField nHatf_ = nHatfv & Sf; thanks andrea |
The face unit interface normal flux...
See the comments here, line 00113ff: http://foam.sourceforge.net/doc/Doxy...8C_source.html Martin |
Yeah,
but what does it means? If i want the geometric curvature of the interface ([1/L]) which is linked to the capillary pressure that i would like to calculate, do i need the flux? Or is OF that needs the flux to calculate the interfacial forces to put in the equation? thanks andrea |
Hi,
I still have a problem when i calculate the curvature. I find high value of curvature where there is no interface (values higher than the values at the interface!!)and that's no possible. what is wrong? thanks andrea |
Unknown patchField type constantAlphaContactAngle for patch type patch in OpenFOAM
hello, foamers
I encounter a problem, when initializes the alpha in 0 files as bellow. i don't know why. anyone can help me. thank you for your help in advance Unknown patchField type constantAlphaContactAngle for patch type patch Valid patchField types are : 62 ( advective calculated codedFixedValue codedMixed cyclic cyclicACMI cyclicAMI cyclicSlip directionMixed empty externalCoupled fan fanPressure fixedFluxPressure fixedGradient fixedInternalValue fixedJump fixedJumpAMI fixedMean fixedPressureCompressibleDensity fixedValue freestream freestreamPressure inletOutlet inletOutletTotalTemperature mapped mappedField mappedFixedInternalValue mappedFixedPushedInternalValue mixed nonuniformTransformCyclic oscillatingFixedValue outletInlet outletMappedUniformInlet partialSlip phaseHydrostaticPressure prghPressure processor processorCyclic rotatingTotalPressure sliced slip symmetry symmetryPlane syringePressure timeVaryingMappedFixedValue totalPressure totalTemperature turbulentInlet turbulentIntensityKineticEnergyInlet uniformDensityHydrostaticPressure uniformFixedGradient uniformFixedValue uniformInletOutlet uniformJump uniformJumpAMI uniformTotalPressure variableHeightFlowRate waveSurfacePressure waveTransmissive wedge zeroGradient ) |
Quote:
Code:
cv.ref() = mesh.V(); |
All times are GMT -4. The time now is 13:54. |