# Reading Velocity field every time step

 Register Blogs Members List Search Today's Posts Mark Forums Read

 March 20, 2015, 11:18 Reading Velocity field every time step #1 Member   Ali Join Date: Aug 2011 Location: Milwaukee Posts: 34 Rep Power: 14 Hi everyone, I'm a newbie in OpenFoam so most likely my question might be really easy to solve, but unfortunately, I couldn't find any answer for my problem so far by searching in different resources. So the problem: I'm solving advection-diffusion equation for a field which the velocity field is known and measured in every time step. So, I need to read the velocity field each time step and then calculate the fluid concentration. So far I can calculate the concentration with edited solver based on icoFoam. But, I don't know how can I read the velocity field instead of calculating it. I really appreciate any and every comment on this problem. Thanks Ali ahmmedshakil and Kummi like this.

 March 24, 2015, 06:37 #2 Senior Member   Pablo Join Date: Mar 2009 Posts: 102 Rep Power: 17 With : forAll(U, celli) { Info << "U aire dife " << mag(U[celli]) << nl; } You can check all the velocities, so u can modify, compute or else Pablo

 March 24, 2015, 09:14 #3 Member   Ali Join Date: Aug 2011 Location: Milwaukee Posts: 34 Rep Power: 14 Hi Pablo, Thank you so much for your respond. I would be grateful if I can have more help from you. So, I'm trying to solve the following equation: φt + ∇.(uφ) = ∇.(D∇φ) where u is known for each time step. So I took icoFoam solver and edited as bellow: Code: int main(int argc, char *argv[]) { #include "setRootCase.H" #include "createTime.H" #include "createMesh.H" #include "createFields.H" #include "initContinuityErrs.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; while (runTime.loop()) { Info<< "Time = " << runTime.timeName() << nl << endl; #include "readPISOControls.H" #include "CourantNo.H" // --- PISO loop for (int corr=0; corr

 March 24, 2015, 10:52 #4 Member   Ali Join Date: Aug 2011 Location: Milwaukee Posts: 34 Rep Power: 14 Hi again, So I found out that the scalarTransportFoam does almost the thing that I'm looking for. The only problem is that, the solver just reads the velocity field as the initial condition but nothing after that. So, I'm looking to read every single step. Any thoughts? Thanks again. Ali

 March 24, 2015, 11:02 #5 Senior Member   Jens Höpken Join Date: Apr 2009 Location: Duisburg, Germany Posts: 159 Rep Power: 17 Just copy the construction of U to the appropriate place in the solver and recompile it. __________________ Blog: sourceflux.de/blog "The OpenFOAM Technology Primer": sourceflux.de/book Twitter: @sourceflux_de Interested in courses on OpenFOAM?

March 24, 2015, 11:28
#6
Member

Ali
Join Date: Aug 2011
Location: Milwaukee
Posts: 34
Rep Power: 14
Quote:
 Originally Posted by jhoepken Just copy the construction of U to the appropriate place in the solver and recompile it.
I did this:

Code:
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "createFields.H"
#include "createFvOptions.H"

simpleControl simple(mesh);

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Info<< "\nCalculating scalar transport\n" << endl;

#include "CourantNo.H"

while (simple.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;
Info<< "Reading field U\n" << endl;

volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::NO_WRITE
),
mesh
);

while (simple.correctNonOrthogonal())
{
solve
(
fvm::ddt(T)
+ fvm::div(phi, T)
- fvm::laplacian(DT, T)
==
fvOptions(T)
);
}

runTime.write();
}

Info<< "End\n" << endl;

return 0;
}
and I'm getting the following error:

Code:
cannot find file

file: /home/cssllab/OpenFOAM/cssllab-2.3.0/run/elbowScalarTransport/10/T at line 0.

in file db/regIOobject/regIOobjectRead.C at line 73.
Thanks again

 March 24, 2015, 11:59 #7 Senior Member   Tomislav Maric Join Date: Mar 2009 Location: Darmstadt, Germany Posts: 284 Blog Entries: 5 Rep Power: 21 OK, I've read the thread. From what I see, I can tell you a few hints. First of all, in OpenFOAM, volumetric fluxes are used for the transport, not the cell centered velocity field. If you have gotten the velocity field from a measurement, you should somehow use it to compute the volumetric flux at the faces. fvm::div(phi, T) Is the convective term, where phi is the volumetric flux. Point #0 Rename your concentration field from \phi to \Psi or something else, otherwise you have a 100% probability of mistaking it for OpenFOAM's \phi. fvm::ddt(Psi) + fvm::div(phi, Psi) - fvm::laplacian(lambda, Psi) == 0 Is this your equation? If yes, then we are at Point #1: Point #1 - You most likely don't need a pressure-velocity coupling algorithm, use the scalarTransportFoam solver as a start. Point #2 - Since the volumetric flux is used for the convective transport, you need to get it from the measured velocity field. Any error you did in the measurement will introduce an error in volume conservation: if you know the pressue-velocity boundary conditions of your experimental setup, at this point you *might* consider solving for p-U to correct the fluxes for volume conservation. But then, you have to measure the pressure at the boundaries with sufficient accuracy, varying in time. Don't know about that.... If you use an interpolation method to evaluate the velocity at the finite volume face to get the volumetric flux, depending on the interpolation, the related interpolation errors will mess up the fluxes, you will loose volume conservation. If you need to reconstruct the volumetric flux from a cell centered velocity field, there are some obscure methods: Radial Basis Functions with divergence free basis, or FEM/FVM hybrids, which is something out of scope I think. Can you align the mesh face centers with the measured velocity field? This way we could avoid interpolation. Is the domain a box? Your basic goal is to pre-process Psi and get phi from U. And yes, you should read U every time step since it's the input for the volumetric flux required by the convective concentration transport term. Can you please comment everything out in the time loop apart from reading U, and report what you get using this line: Info << fvc::average(fvc::div(U)).value() << endl; If U comes from a measurement, I expect something very different than 1e-15. The explicit divergence operator fvc::div, when called on the volumetric field, redirects to the divScheme. This one does the following: Code:  tmp < GeometricField ::type, fvPatchField, volMesh> > tDiv ( fvc::surfaceIntegrate ( this->mesh_.Sf() & this->tinterpScheme_().interpolate(vf) ) ); tDiv().rename("div(" + vf.name() + ')'); return tDiv; 'tinterpScheme_().interpolate(vf)' - Your cell centered velocity field from a measurement will be interpolated based on the entry in the case directory, interpolationSchemes, in the dictionary system/fvSchemes. Trust me, you need to have a hell of a fine mesh not to break volume conservation there. Any interpolation error different than machine tolerance will introduce volume conservation errors in the concentration transport. Also, paste only the header part of the file /home/cssllab/OpenFOAM/cssllab-2.3.0/run/elbowScalarTransport/10/T jhoepken, atulkjoy and faiazk like this. __________________ When asking a question, prepare a SSCCE.

 March 24, 2015, 15:38 #8 Member   Ali Join Date: Aug 2011 Location: Milwaukee Posts: 34 Rep Power: 14 Tomislav, Thanks for your great response. First off, my data is from machine measurement, means it comes in a structure mesh format. But I still have velocities at the center of each hex which I guess in my case interpolation works fine. About the test you asked me, I believe I'm doing something wrong. So, I'm trying to calculate the concentration of a fluid inside another fluid which do not react to each other. I have velocity for the main stream and now wants to calculate the concentration of the second with a source in some boundary condition. First let me ask this, how should I feed velocity fields to the solver? To be more clear, currently I am making time step folders manually and put the velocity field of each time step in the folder. Which seems OpenFOAM looks for all fields as soon as finds the excising folder. As the first step, can you pleas walk me through this? Should I some how (which I don't know how) put all available velocity fields from all time steps in folder 0 or what? As the second step, as you also mentioned, I already started to working on scalarTransformFoam. So that was the reason that I used T as the concentration. and phi is OpenFoam's phi. So, can we work through this step by step if it's possible for you? P.S. I'm getting following error for the command you asked to run Info << fvc::average(fvc::div(U)).value() << endl; has no member named ‘value’ Thanks again. Ali

 March 24, 2015, 15:44 #9 Senior Member   Tomislav Maric Join Date: Mar 2009 Location: Darmstadt, Germany Posts: 284 Blog Entries: 5 Rep Power: 21 Hi, I was writing in a hurry and forgot max: Code:  Info << max(fvc::average(fvc::div(U))).value() << endl; Also, in system/fvSchemes: div(U) Gauss linear; I'll get back to you later on the rest of your question. jhoepken likes this. __________________ When asking a question, prepare a SSCCE.

March 25, 2015, 04:27
#10
Senior Member

Tomislav Maric
Join Date: Mar 2009
Posts: 284
Blog Entries: 5
Rep Power: 21
I can't give you a step-by-step answer in source code - I would have to re-create your problem, use your data, and try this stuff.

But I can tell you approximately what you could try.

I don't know what you did manually, but be careful: OpenFOAM uses unstructured meshes, so you must be sure that the position of the U vector in a list corresponds to the proper label of the cell where U is stored. If you have a list of vectors, and you store it in a time-step directory like this

Code:
... header

internalField   nonuniform List<vector>
16384
(
(1.03544e-05 -1.06153e-05 0) // Cell 0
(3.06175e-05 -1.01041e-05 0) // Cell 1
(5.11404e-05 -9.46919e-06 0) // Cell 2
(7.20913e-05 -9.58538e-06 0) // ...
(9.31207e-05 -9.87433e-06 0)
(0.000115088 -1.02901e-05 0)
Then the velocity vector that is first in the list will map to cell 0, second velocity vector will correspond to cell 1, and so on. Now, cell 0 can be practically anywhere in a general sense, when the mesh is polyhedral or tetrahedral. With box domains meshed using blockMesh, there is a structured ordering of cells, so you can make use of that information. Open thsi image

http://postimg.org/image/jruimkh0d/

to see what I mean. The cell IDs start at zero and are increased from left to right until the end of the row.

You need to absolutely make sure that the U internal field vectors correspond to proper cell IDs.

I am stil convinced that for coarse meshes, if you use interpolation to get the flux, from a measured velocity field, you will introduce interpolation errors on top of measurement errors, and your passive scalar transport will not be volume conservative. If you can shift the measurement so that you measure the velocities directly at the face centers, you will not have any interpolation errors at least, just the measurement errors. It's kind of like the staggered grid approach.

Once you do this, you can simply initialize T concentration field, and compute phi as the dot product of the measured face centered velocity and Sf. Something like this (I am writing this directly, not debugged):

Code:
surfaceVectorField Uf
(
IOobject
(
"Uf",
runTime.timeName(),
mesh,
IOobject::NO_WRITE
),
mesh
);

phi = Uf & mesh.Sf();
Instead of reading in the U velocity, you would read Uf. Not a big difference. The issue there is then to map the read Uf to face centers, that is more problematic, but much more accurate on coarse meshes.

If you decide to stick with cell centered velocities for a first go at it, just order them based on the image I uploaded. You know the position and ordering of the measurement, so there you go. Or you can even extract measured U directly in the order as shown in the image and that's it.

Hope this helps,
Tomislav
Attached Images
 block-mesh-ids.jpg (73.9 KB, 31 views)
__________________
When asking a question, prepare a SSCCE.

 March 25, 2015, 09:25 #11 Member   Ali Join Date: Aug 2011 Location: Milwaukee Posts: 34 Rep Power: 14 Tomislav, Thank you so much for your detailed answer. #1: For the interpolation part, I'm totally agree with you. If you just simply interpolate data you will not get accurate data. To be more clear using schemes like upwind, central difference etc.. But, we solved this problem using QUICK. You can find my report on the numerical solution in the following link: https://www.overleaf.com/read/vdmryrypjzhw So, I solved that part already since I even solved my problem in Matlab and it's fully functioning. As now, I'm looking to have the solver in OpenFOAM for bigger domains. #2, About the manual part that I mentioned. I structured my U field same as you mentioned but in the following structure. So, initially, I have velocities for each time step in a folder and trying to solve for T for the same time step using available U. But, the problem is, with this structure as soon as I run the case it gives me error that it cannot find T file in the folders. Code: |-- 0 | |--T | --U |-- 1 | --U |-- 2 | --U |--constant --system I run so many test with available tutorials. Seems to me, OpenFOAM WILL NOT solve the time step that the folder is already available. Even if I put the starttime as zero. It will continue to solve from the latest available folder. So, the question for me here is how to force OpenFOAM to solve and ignore available folders and just read the U field from them. As the record, my velocity and concentration structures are as bellow: Code: volScalarField T ( IOobject ( "T", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ); Info<< "Reading field U\n" << endl; volVectorField U ( IOobject ( "U", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::NO_WRITE ), mesh ); Thanks again. Ali Last edited by alib022; March 25, 2015 at 12:57.

 March 26, 2015, 02:14 #12 Senior Member   Jens Höpken Join Date: Apr 2009 Location: Duisburg, Germany Posts: 159 Rep Power: 17 I suppose that you would have to change the Time class for that, which is something you don't want to do, as it is a fairly essential component of OpenFOAM. __________________ Blog: sourceflux.de/blog "The OpenFOAM Technology Primer": sourceflux.de/book Twitter: @sourceflux_de Interested in courses on OpenFOAM?

 March 27, 2015, 09:20 #13 Member   Ali Join Date: Aug 2011 Location: Milwaukee Posts: 34 Rep Power: 14 Thanks, So, is that the only way? I tried to understand the Time class but it's way to complicated!

 November 24, 2015, 21:12 #14 New Member   amin enr Join Date: Sep 2015 Posts: 4 Rep Power: 10 Hi Ali, My problem is exactly the same as yours, I want to read velocity field in each time step, but I haven't found any way for that, so far! What I am thinking of as a possible way, is to write (Allrun) file which does the following: 1. copies the velocity field in each time step to the current time folder 2. runs the solver for one time step 3. Shifts the start and end time in "controlDict" for one step I haven't found out how to run the code for (Allrun), but I am trying to do that. Please inform me whenever you find a solution for your problem. Amin

 November 25, 2015, 10:16 #15 Member   Ali Join Date: Aug 2011 Location: Milwaukee Posts: 34 Rep Power: 14 Hey Amin, My project changed a bit after that post. But, for another similar project I did what you are trying to do. I used octave as my Allrun controler, in that way I could copy and paste and also was able to do all sort of edits in my controlDict file. Please let me know if you need more details then maybe I can find my files and send you some part of the code. Best Ali

 December 1, 2015, 14:28 #16 New Member   amin enr Join Date: Sep 2015 Posts: 4 Rep Power: 10 Hi Ali, I haven't worked with octave and I don't have any idea about that, but I will try to find out. Anyway, it would be a great help if you can find your files and codes, because I am really confused with writing my own Allrun file. Thanks, Amin

 December 6, 2015, 11:05 #17 Member   Ali Join Date: Aug 2011 Location: Milwaukee Posts: 34 Rep Power: 14 Amin, Here is the code that I wrote in Octave for editing velocity profile as well as ControlDic file. This is not the proper way to it, but because of nature of my project I had to it in the hard way. I put here just the parts relevant to editing the controlDict file. Code: ######################## S T A R T ############################### # ---> Write back the edited controlDict #------------------------------------------------------------ fileName = "controlDict"; fidControlDict = fopen(fileName, "wt"); Line1 = ["/*--------------------------------*- C++ -*----------------------------------*\\ \n"]; Line2 = ["| ========= | | \n"]; Line3 = ["| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | \n"]; Line4 = ["| \\ / O peration | Version: 2.4.0 | \n"]; Line5 = ["| \\ / A nd | Web: www.OpenFOAM.org | \n"]; Line6 = ["| \\/ M anipulation | | \n"]; Line61 =["|-----------------------------------------------------------------------------| \n"]; Line7 = ["\\*---------------------------------------------------------------------------*/ \n"]; Line8 = ["FoamFile \n"]; Line9 = ["{ \n version 2.0; \n format ascii; \n class dictionary; \n"]; Line91 = [ 'location "system";']; Line92 = [" \n object controlDict; \n } \n"]; Line10 = ["// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // \n"]; Line11 = [" application icoFoam; \n startFrom startTime; \n"]; Line12 = [" startTime \t " num2str(folderNum) "; \n"]; Line13 = ["stopAt endTime; \n endTime \t " num2str(folderNum+0.01) "; \n" "deltaT 0.0005; \n" "writeControl timeStep; \n" "writeInterval 20; \n"]; Line14 = ["purgeWrite 0; \n" "writeFormat ascii; \n" "writePrecision 6; \n" "writeCompression off; \n" "timeFormat general; \n" "timePrecision 6; \n"]; Line15 = ["runTimeModifiable true; \n" "// ************************************************************************* // \n"]; Line = [Line1 Line2 Line3 Line4 Line5 Line6 Line61 Line62 Line63 Line7 Line8 Line9 Line91 Line92 Line10 Line11 Line12 Line13 Line14 Line15]; fputs(fidControlDict, Line); fclose(fidControlDict) #----------------------------------------------------------- # <--- End of writing back the edited controlDict ######################## E N D ############################## copyU = ["cp -a controlDict" " Run_" num2str(kLoop) "/Ensemble_" num2str(numEn) "/system/"]; unix(copyU); unix("rm controlDict"); # Runing OpenFOAM for each case OpenFOAMString = ["icoFoam -case" " Run_" num2str(kLoop) "/Ensemble_" num2str(numEn)]; unix(OpenFOAMString);` I believe with some edit you can run this code in Matlab. Regards Ali atulkjoy, Zhiheng Wang and Ramzy1990 like this.

 December 9, 2015, 17:57 #18 New Member   amin enr Join Date: Sep 2015 Posts: 4 Rep Power: 10 Many thanks for your code Ali. I was able to run my problem based on that.

 December 10, 2015, 10:41 #19 Member   Ali Join Date: Aug 2011 Location: Milwaukee Posts: 34 Rep Power: 14 Amin, Again, this is not the best way to do it. This code was a part of much bigger code that I had to write in Octave, since you dont have to do anything in Octave you can simply use bash file to do your work in much elegant way. Best Ali

 December 23, 2019, 18:42 #20 New Member   jhydome Join Date: Nov 2018 Posts: 2 Rep Power: 0 Ok, I know this thread is pretty old, but for the case that someone still has the same problem, such as me, I recommend this thread: Operation to READ_IF_PRESENT field post#6. The basic idea is that since the velocity field U has been created with creatField.h before the time loop, OpenFoam will not read from the data again during the time loop. So we need to create a new field UNew in each step, read the file to write in UNew and let U=UNew. Mine seems to work normally, good luck! Roozbeh_Sa likes this.