[alib022]

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

[pablodecastillo]

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

[alib022]

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:

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;
}

and for the createFields.H

    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);

and I'm getting following error:

cannot find file

file: .../elbowScalarTransport/10/T at line 0.

So that's the whole point of the solver. Reading U and calculate T. Can you please tell me what am I doing wrong?

Thanks again
Ali

[alib022]

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

[jhoepken]

Just copy the construction of U to the appropriate place in the solver and recompile it.

[alib022]

Originally Posted by jhoepken

Just copy the construction of U to the appropriate place in the solver and recompile it.

I did this:

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;
}

and I'm getting the following error:

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.

Thanks again

[tomislav_maric]

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.

1. 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....
2. 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:

    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

[alib022]

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

[tomislav_maric]

Hi,

I was writing in a hurry and forgot max:

   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.

[tomislav_maric]

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

... 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):

surfaceVectorField Uf 
( 
   IOobject
   (
       "Uf", 
       runTime.timeName(), 
       mesh, 
       IOobject::MUST_READ, 
       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

[alib022]

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.

|-- 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:

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

[jhoepken]

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.

[alib022]

Thanks,

So, is that the only way? I tried to understand the Time class but it's way to complicated!

[amin.enr]

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

[alib022]

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

[amin.enr]

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

[alib022]

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.

######################## 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

[amin.enr]

Many thanks for your code Ali.
I was able to run my problem based on that.

[alib022]

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

[crazybarney]

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!