[danny123]

Hi,

I have a problem to solve in handling time step behavior in OpenFoam. For my problem, the required number of iterations as well as initial residuals in interDymFoam increases exponatiallly if the time step is time step is too large. I want to cancel this time step and divide it, say by 2 by an adaptive method.

Changing the time step is easy, I also can save U, p_rgh, phi and alpha1 for reuse later (there is likely a whole lot more required, including the mesh movement and position). But I have trouble to understand how OpenFoam handles time progress. The code lines are:

    while (runTime.run())
    {
        #include "readControls.H"
        #include "alphaCourantNo.H"
        #include "CourantNo.H"

        #include "setDeltaT.H"

        runTime++;

So, I need to uunderstand the object "runtime". Where do I find that definition. I tried grep, but no luck.

Regards,

Daniel

[danny123]

I did a couple trials. I added some code in interDymFoam.C:

while (runTime.run())
    {
        #include "readControls.H"
        #include "alphaCourantNo.H"
        #include "CourantNo.H"

        #include "setDeltaT.H"

        Info<< "deltaT = " <<  runTime.deltaTValue() << endl;
        Info<< "1 Time = " << runTime.timeName() << nl << endl;

        runTime++;

        Info<< "2 Time = " << runTime.timeName() << nl << endl;

        scalar deltaTold = -runTime.deltaTValue();

        Info<< "deltaTold = " << deltaTold << nl << endl;

        runTime.setDeltaT
        (
            deltaTold
        );

        runTime++;

        Info<< "deltaT = " <<  runTime.deltaTValue() << endl;
        Info<< "3 Time = " << runTime.timeName() << nl << endl;

The output then is:

Interface Courant Number mean: 0 max: 0
Courant Number mean: 6.84018e-09 max: 8.1219e-08
deltaT = 9.90099e-05
deltaT = 9.90099e-05
1 Time = 0

2 Time = 9.90099e-05

deltaTold = -9.90099e-05

deltaT = -1.9802e-05
3 Time = 7.92079e-05

As can be seen, runtime++ works ok even for negative numbers. So I can go back in time to redo my calculation. However, runTime.setDeltaT sets a value 1/5 of what I want. If I use a positive value, it works ok. Does anybody know why?

Regards,

Daniel

[danny123]

Ok,

I found part of the solution on this forum:

http://www.cfd-online.com/Forums/ope...tml#post298878

The factor of 5 for negative deltaT is because of a routine check in Time.C called adjustDeltaT. This prooves that negative deltaT is not "recommended".

So, at least I can recover the last time step. In order to recover the fields, I can store them beforehand in createFields.H separatly declared fields and overwrite them with those. This should work ok.

Still, I would need also to recover the old mesh, basically return to the condition prior to

mesh.update;

Does anybody know if mesh.update uses deltaT or does it use the time or even both? I will do some tests next week, but if there is somebody who has done this previously, I would appreciate some comment.

Regards,

Daniel

[chriss85]

In case anyone is still interested in this topic (as I have been just now because the Courant number is not always a good indicator for required time step size), here's how you revert the time step properly:

//make a backup of the current state
TimeState tSCurrent(runTime);
runTime++;
//restore the previous state
runTime.TimeState::operator=(tSCurrent);
//maybe set a different deltaT
runTime.setDeltaT(runTime.deltaT().value() / 2);
//And increase again
runTime++;
//If you are inside a convergence loop in your current time step, you might want to go to the next iteration
continue;

Doing this preservers the time value, deltaT, time index and some other things from the TimeState class. I'm not sure if other things need to be backed up, but so far this has been working for me.

You should also consider that if your erroneous solution diverges too much from the real result, it might not converge with the incorrect fields and a smaller time step. In this case, you also need to back up the fields from the previously valid time step and restore them. Alternatively, you can try to limit the problematic fields to reasonable boundaries to help with the convergence.

[danny123]

Thanks Chris,

I have implemented it in a slighly different way. It seems to work, but I have to debug it within the next days to see if it does what I expect:

To safe the previous time step, I use:

time_start = runTime.timeOutputValue();

Later, I can restore by the following way:

runTime.setTime
                ( 
                    time_start,
                    time_start
                );

I wonder what the difference between your code is.

The vectors of the previous time step are stored. I think this is mandatory, since you do not want to start by a guess value that has diverged.

In order to control the time step, I run setDeltaT.H not upfront, but after the 1st try. If deltaT increases, it is divided by 2 within a loop and repeats until it does. So, basically you get an implit scheme for setting the time step instead of an explicit, as it is implimented by default. So, I can use Co number, number of iterations or whatever else to correct the time step (later is an extension within setDeltaT.H).

Another 2 things to consider: never run runtime.write(), before your time step is not definitive since then you cannot cancel your time step anymore (updates timestate) and
the mesh is updated on time not deltaT. So, I can keep the code of the mesh update as it is.

The old thread is defininatly wrong. Did I forget anything?

Thangs and regards,

Daniel

[chriss85]

I think that the stored deltaT and deltaT0 values and possibly the index are not restored with your method. There is one setTime method that tries to read these values from a stored result but this only works if the previous step was written to disk afaik. I'm not sure if they are used for anything else other than calculating the d/dt terms in the equations. However, if you use a second order time scheme the deltaT0 term might be important?

I also believe that the fields from the previous time step are always used for solving the equations, they should be inside the d/dt terms in your equations. The fields in the equations themselves will not get modified if you just work on the time object. I'm not really sure on this so please correct me if I'm wrong.

Another thing, in my case I detect too large time steps by looking at the minimum/maximum temperatures. If my time step gets too large I see negative or very high temperatures which are definitely wrong. I also wait for 3 iterations and only decrease the time step if the temperature was out of bounds for 3 iterations in a row, because this problem can also appear once after decreasing the time step.

Oh and if you use an algorithm like PIMPLE you should also break the iteration by looking at the residuals instead of a fixed number of steps. The number of steps should be high enough to reach convergence even if the time step had to be decreased once or twice. It might be possible to also reset the PIMPLE iteration index, I haven't checked that yet.

[danny123]

Hello Chris,

The index is stored as long as it is the same as the value. Since the value is a scalar, overwriting the index by a scalar is kind of a "dirty" trick. But I have not figured out how to get the index. I have plotted the result after runtime.write() and the index got updated.

I think the reason that there is an index after all is only that the time folders need to be named by a string. So, the index is propably used for that purpose.

There is a method to get the previous time step value, which is xx.oldTime(). This would save some memory using this. I have not used it for now just to test the working principle. But this can be done later.

I use residual control to automatically break the number of PIMPLE loops. I have an "optimal" number of required PIMPLE loops that determines deltaT. So, if the required PIMPLE number is bigger than the set value, deltaT is decreased. Same applies for the number of inner iterations.

Regards,

Daniel

[chriss85]

You can get the index with runTime.timeIndex().

I'm not sure if reverting the fields by using (simplified) field = field.oldTime(); is valid, because the oldTime field may be lacking the oldTime().oldTime() information which might be required for the discretization of the d/dt terms.
If you intend to do this, you should try to figure this out. For my case I get along just fine with using the wrong results as a starting point for an iteration with a smaller time step. If the equations converge than the results should be correct, it's only problematic if you can't reach convergence with the wrong results.

How did you determine your optimal number of PIMPLE iterations? Did you measure the total runtime with different time step sizes? Do you think your method is applicable to general problems? Are you getting accurate results and is your Courant number low enough?

[danny123]

Thanks for the hint. I do not know the exact d/dt scheme working principle, but I think the final (confirmed) results of the previous time step are necessary to establish the differential equations. Since this procedure is repeated with every PIMPLE loop, this information has to be stored somewhere, unless you calculate them back from the equations, which would need the solution of another system of equation using a matrix solver...

Since the current results have to be set to old at some point at the end of the PIMPLE loop, it is likely runtime.write(), which does just that. You then can set your fields however you want, preferably close to the solution to be iterated. The fields of the previous time step is a possible, but not necessarily the best starting point. I use a projection method to improve that guess fields. For that method to work, I have to return to where I was before and restart over.

I have not optimized the exact number of PIMPLE loops yet, but it seems to be anywhere around 5 to 10. Beyont this, the calculation time increases a lot and results and stability goes down too. I think it depends also what resisuals you are targeting.

A big problem is stability. I have a very high viscosity, which means that an inacurate pressure equation leads to instabilities for the next time step calculations. And the problem is transient. Optimizing Co numbers has not helped a lot so far.

Regards
Daniel

[chriss85]

Looking through the code I wasn't able to find the part where the old time fields are replaced yet. However, the d/dt schemes access the field.oldTime() members and the mesh.time().deltaT() and mesh.time().deltaT0() values. For example, the euler scheme is (field - field.oldTime()) / mesh.time().deltaT(). This means that if you use something else than a first order time scheme it would be better to backup the whole TimeState. I also saw some code based on the time index in the GeometricField/Mesh classes so this needs to be handled properly as well.

Have you seen good results with extrapolating your results from the previous times? I had been thinking about this as well sometime, but never really got around to look into it properly. I would expect that this is only useful for "scaling" fields, but not for convective fields.

I usually have some more PIMPLE loops, maybe around 10-20, but this depends on the problem. Stability is not so much an issue for me anymore now that I can dynamically adjust the time step if needed, but I find that the required relaxation factors for a converging solution can change quite a bit. If you are interested in this, I've recently written about this topic here: http://www.cfd-online.com/Forums/ope...imization.html . I welcome any meaningful comments

I don't have much experience regarding solution behavior with different viscosities but I know that for highly compressible flows density based solvers can be more suitable than pressure based solvers. I have looked into rhoCentralFoam before but I fear that it's much more difficult to adapt to my problem than sonicFoam.

[danny123]

The code seems to work by now. The only (unrelated) problem is that I stop the PIMPLE loop at one point if the solution does not converge. This means that pimple.loops() is not reinitiated, so the counter of PIMPLE continues to grow instead of returning back to 0. I could reconfigure the whole Pimple thing. Do you know how the Pimple counter can be reset?

Extrapolating previous results works ok, not much better though. I think it is good practice starting an interpolation using the best guess value you can come up with. If the behavior of the flow does not oszilate from time step to time step, an extrapolation should work. If it does, the chosen time step is too big anyway.

Incompressible solvers are much easier to implement as compressible. I will eventually go in this direction for viscoelastic flows, but first things first. The pressure equation needs to be solved properly since we need the forces since this is what we can measure on the mixer. So, it is impossible to get around the problem of precision of the pressure equation.

Regards,

Daniel

[chriss85]

Not sure about PIMPLE, I usually just set the iteration number high enough so it doesn't run out of loop iterations.

In density based solvers (like rhoCentralFoam) there is no pressure equation, instead you derive the pressure from the equation of state. These solvers provide good results for cases in which high compression occurs, such as supersonic flows. Pressure based solvers (like sonicFoam) are more suited to cases in which the compressibility is not so strong.

[danny123]

Hello,

My code meanwhile works ok. The balance of specy (alpha1eq) gets residual number that are quite lower and the residual control feature works now. However, I had to change the setDelta implementation from

runTime.setDeltaT (deltaT );

to

runTime.setDeltaT (deltaT , false);

The point is, the adjustDeltaT changes the value of delta in manner that I could not reproduce. What I found is that the deltaT you get does not match a value that fits to a time step hitting a write interval, but mostly something a couple of promille smaller value than the one, which is specified or even to the value previously in place (so no change is being made).

My own thinking is that the deltaT to be set should be the one specified, but not more than the remaining gap to the next write interval. After this time step, the solver should return to the original specified deltaT. But this is not what OpenFoam seems to do.

Regards,

Daniel

[danny123]

So, I have implemeted what I thought is correct and it works. The only remaining issue is how to read the writeInterval into the solver. It should be read from controlDict, but I did not find any object in Time.C of that kind.

This is strange since objects such as the start time can be retrieved easily:

scalar startTime_ = runTime.startTime().value();

I think I could define such an object in time.C by copying the way it is done for startTime, but is there an easier, direct way?

Thanks for the help,

Regards,

Daniel

[chriss85]

I think you can find that in one of the parent classes, TimeIO.C (?).

[danny123]

Thanks for the answer. I settled on adding it in Time.C. This is certainly the most elegant way.

Regards,

Daniel