|
[Sponsors] |
Reading Velocity field every time step |
![]() |
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
![]() |
![]() |
#1 |
Member
Ali
Join Date: Aug 2011
Location: Milwaukee
Posts: 34
Rep Power: 15 ![]() |
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 |
|
![]() |
![]() |
![]() |
![]() |
#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 |
|
![]() |
![]() |
![]() |
![]() |
#3 |
Member
Ali
Join Date: Aug 2011
Location: Milwaukee
Posts: 34
Rep Power: 15 ![]() |
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<nCorr; corr++) { solve ( fvm::ddt(T) + fvm::div(phi, T) - fvm::laplacian(nu, T) ); #include "continuityErrs.H" } runTime.write(); Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s" << nl << endl; } Info<< "End\n" << endl; return 0; } Code:
Info<< "Reading transportProperties\n" << endl; IOdictionary transportProperties ( IOobject ( "transportProperties", runTime.constant(), mesh, IOobject::MUST_READ_IF_MODIFIED, IOobject::NO_WRITE ) ); dimensionedScalar nu ( transportProperties.lookup("nu") ); Info<< "Reading field U\n" << endl; volVectorField U ( IOobject ( "U", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::NO_WRITE ), mesh ); Info<< "Reading field T\n" << endl; volScalarField T ( IOobject ( "T", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ); # include "createPhi.H" label TRefCell = 0; scalar TRefValue = 0.0; setRefCell(T, mesh.solutionDict().subDict("PISO"), TRefCell, TRefValue); Code:
cannot find file file: .../elbowScalarTransport/10/T at line 0. Thanks again Ali |
|
![]() |
![]() |
![]() |
![]() |
#4 |
Member
Ali
Join Date: Aug 2011
Location: Milwaukee
Posts: 34
Rep Power: 15 ![]() |
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 |
|
![]() |
![]() |
![]() |
![]() |
#5 |
Senior Member
|
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? |
|
![]() |
![]() |
![]() |
![]() |
#6 | |
Member
Ali
Join Date: Aug 2011
Location: Milwaukee
Posts: 34
Rep Power: 15 ![]() |
Quote:
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::MUST_READ, 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; } Code:
cannot find file file: /home/cssllab/OpenFOAM/cssllab-2.3.0/run/elbowScalarTransport/10/T at line 0. From function regIOobject::readStream() in file db/regIOobject/regIOobjectRead.C at line 73. |
||
![]() |
![]() |
![]() |
![]() |
#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.
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 <typename innerProduct<vector, Type>::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
__________________
When asking a question, prepare a SSCCE. |
|
![]() |
![]() |
![]() |
![]() |
#8 |
Member
Ali
Join Date: Aug 2011
Location: Milwaukee
Posts: 34
Rep Power: 15 ![]() |
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 |
|
![]() |
![]() |
![]() |
![]() |
#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; div(U) Gauss linear; I'll get back to you later on the rest of your question.
__________________
When asking a question, prepare a SSCCE. |
|
![]() |
![]() |
![]() |
![]() |
#10 |
Senior Member
Tomislav Maric
Join Date: Mar 2009
Location: Darmstadt, Germany
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) 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::MUST_READ, IOobject::NO_WRITE ), mesh ); phi = Uf & mesh.Sf(); 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
__________________
When asking a question, prepare a SSCCE. |
|
![]() |
![]() |
![]() |
![]() |
#11 |
Member
Ali
Join Date: Aug 2011
Location: Milwaukee
Posts: 34
Rep Power: 15 ![]() |
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 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 ); Ali Last edited by alib022; March 25, 2015 at 12:57. |
|
![]() |
![]() |
![]() |
![]() |
#12 |
Senior Member
|
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? |
|
![]() |
![]() |
![]() |
![]() |
#13 |
Member
Ali
Join Date: Aug 2011
Location: Milwaukee
Posts: 34
Rep Power: 15 ![]() |
Thanks,
So, is that the only way? I tried to understand the Time class but it's way to complicated! |
|
![]() |
![]() |
![]() |
![]() |
#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 |
|
![]() |
![]() |
![]() |
![]() |
#15 |
Member
Ali
Join Date: Aug 2011
Location: Milwaukee
Posts: 34
Rep Power: 15 ![]() |
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 |
|
![]() |
![]() |
![]() |
![]() |
#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 |
|
![]() |
![]() |
![]() |
![]() |
#17 |
Member
Ali
Join Date: Aug 2011
Location: Milwaukee
Posts: 34
Rep Power: 15 ![]() |
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 |
|
![]() |
![]() |
![]() |
![]() |
#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. ![]() |
|
![]() |
![]() |
![]() |
![]() |
#19 |
Member
Ali
Join Date: Aug 2011
Location: Milwaukee
Posts: 34
Rep Power: 15 ![]() |
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 |
|
![]() |
![]() |
![]() |
![]() |
#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!
|
|
![]() |
![]() |
![]() |
|
|
![]() |
||||
Thread | Thread Starter | Forum | Replies | Last Post |
Multiple floating objects | CKH | OpenFOAM Running, Solving & CFD | 14 | February 20, 2019 09:08 |
Micro Scale Pore, icoFoam | gooya_kabir | OpenFOAM Running, Solving & CFD | 2 | November 2, 2013 13:58 |
How to write k and epsilon before the abnormal end | xiuying | OpenFOAM Running, Solving & CFD | 8 | August 27, 2013 15:33 |
plot over time | fferroni | OpenFOAM Post-Processing | 7 | June 8, 2012 07:56 |
calculation diverge after continue to run | zhajingjing | OpenFOAM | 0 | April 28, 2010 04:35 |