[Pyrokrates]

Hey, I would like to implement a inner loop for simpleFoam but if I do, the results differ from each other:


I added just a while loop to compute UEqn and pEqn:

while (simple.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        int i = 0; // new line
        while(i<5) // new line
        {
            // Pressure-velocity SIMPLE corrector
            {
                #include "UEqn.H"
                #include "pEqn.H"
            }

            laminarTransport.correct();
            turbulence->correct();

            i++; // new line
        }

        runTime.write();

        runTime.printExecutionTime(Info);
    }

What I would expect is that for this code (dt = 1), the time step 1 should be equal to time step 5 without the inner loop. Normally I would do 5 simple loops and in my case I would do just one time loop but 5 times... Does anybody know why?


Another question:
How can I print/set relaxationFactors (which are defined in fvSolution) inside the code to change it in my loop? I know that relaxation is applied by e.g. UEqn.relax() but how can I reset it?


Thanks in advance for any help.


Pyro

[Tobermory]

Chris

remember that the SIMPLE algorithm is an iterative algorithm - it does not get to a converged solution in one iteration. Indeed, it can't do this if you have under-relaxation set. What you are seeing in your successive inner loops is just convergence of the solution. I guess you are trying to implement an alternative to the PIMPLE scheme? Have you looked at the PIMPLE scheme?

As for changing the RFs during a run - these are read from the fvSolution dictionary .. you can change these manually "on the fly" during a run, by editing the fvSolution file. I dare say that you could probably code up something using a coded function object to update the contents of the dictionary file ... I don't really recommend it though. My approach is to control the RFs from outside the run. Again, you could write a script to monitor the residuals and adjust the RFs ... you may spend more time getting the logic right, though, than just doing it manually!

Good luck!

[Pyrokrates]

Hey,


thank you for the answer.


What I´m trying to do is changing the relaxation on the fly... I´m using porous media where I get oscillations inside my velocity right next to the porous media. If I lower the relaxation, I can reduce this oscillations... Now I have the problem, that my velocity changes locally over time so I need "high relaxation" for letting U change and "low relaxation" for the correct caluclation of my U field. I know that SIMPLE loop is interative... thats why I would like to make one iterative loop to lets say more iterative loops where the relaxation changes after convergence.


If I do it manually e.g. for my initial solution I run normal simpleFoam with Urelax = 0.1... after 500 time steps I change my Urelax to 0.001 and do 1000 more steps ... until I get the correct Ufield... That is exactly what I would like to do for my time simulation (I don´t want to simulate 3600 seconds with dt=0.1 by manually change of relaxation for each step).


Thanks in advance


Pyro

[Tobi]

In my opinion your for-loop does not make sense at all. The only difference is, that you will get a faster steady-state solution during the calculation == speedup factor of 6 in your case, as you do 6 more iterations (but its not a speed-up in terms of run-time; simply you do 6 times less simple-loops; the simple-loop is simply a for-loop).

The simple-loop itself, is just an iterative algorithm around UEqn.H and pEqn.H.

So if you do 6 simple loops, you should get the same result compared to 1 simple loop and your additional 6 loops. Both do make 6 loops around UEqn.H and pEqn.H. As we don´t have any time-information here, its simply as follows:

1. Construct U-Matrix and update U (if momentumPredictor = true)
2. Construct p-Matrix and solve it, correct U based on new fluxes

With the new p-U fields, you start again on the first point and construct the momentum matrix. Here, and this is the main point, it does not matter if you go one SIMPLE iteration or one of your iterations. Hence, each iteration in your for-loop will differ except if you reached the steady state solution. Your for-loop will not gain anything in your case.
Furthermore, the relaxation factors can be set by using the <Function1> Class. So you have more freedom to set a time-dependent or what ever-dependent relaxation factor. Despite of that you can use Doxygen to search for the function you are searching:

https://www.openfoam.com/documentati...1fvMatrix.html
https://www.openfoam.com/documentati...ebf646011c212e

As you can see, the relax() function is defined in the fvMatrix class. We also have the <fvMatrix>.relax(const scalar) function. You simply can add any number here which is used as relaxation factor rather than the one provided and stored in the class (read from fvSolutions).

[Tobermory]

I have to confess that I agree with Tobi - I am confused by your discussion of "time" when you are using simpleFoam ... it's a steady state solver. Perhaps you need to tell us a little more about the problem you are trying to solve.

[Tobi]

Private message from Chris

Hey Tobi,


thanks for your answers and the help.


To go a little bit more into detail: What I would like to do is calculating concentrations inside inhomogeneous material. Therefore I programmed another equation system (lets call it conEqn.H). The new equation needs a stable velocity right next to my porous media to calculate my results correctly. Therefore I need low relaxationFactors of my simple loop as is mentioned in my post. But as the concentration spread over my porous media, the velocity changes locally (areas with low velocity increase speed others vice versa). Thus I need higher relaxation first.


I know that simpleLoop is not about time cause it iterates until the solution is converged for the specific time step... My problem is, my simulation is more about the concentration calculation meaning if I run my simulation with 1h = 3600s and a time step of 1s e.g., than I need to get my new correct flow field inside one time step. That´s why I would like to do some inner loop which is (I know) like more time steps in normal simpleFoam... But for my specific case I don´t want t=1,2,3,4,... for one convergence... I want t=1 for one convergence, t= 2 for the next and so on...


Hope you can help me with this


Thanks


Chris

I agree totally with Tobermory (sorry you don´t have your name here and probably I should now it). You do not provide enough information. However, in between the lines and based on your private message, for me it seems that you want to calculate a concentration distribution (time-based) on a frozen (steady-state) velocity field, right?

If so, simply use simpleFoam and converge your solution (only p/U). Use the U or phi field of the converged solution for a second run using e.g., scalarTransportFoam or use simpleFoam -postProcess xyFO to calculate the concentration based on a passive transport equation on the fly (FO); here you could also add source terms if necessary. If you are interested in time-based data, you need to change the ddt scheme for the passive scalar.

If you want to have an in-time velocity field and an in-time concentration field, you cannot use simpleFoam (maybe SIMPLEC + almost no relaxation + time-derivative activated; but a quick look into the source code of simpleFoam gives us the hint, that there is no time derivative implemented, so its not working and not worth to investigate). Therefore, you need to use pimpleFoam. Here, you can use the PIMPLE algorthm to converge between each time step the p-U correlation much more (but this is only needed if the Courent number is way much larger than 1; see my book for the usage of pimple algorithm for example)...


If you want to do anything else, you need to clarify your statements and what you want to do. Out of the box, I have no other idea what you would like to do.

[Pyrokrates]

Hey,


thank you for the interesting and helpful input.

Originally Posted by Tobi

In my opinion your for-loop does not make sense at all. The only difference is, that you will get a faster steady-state solution during the calculation == speedup factor of 6 in your case, as you do 6 more iterations (but its not a speed-up in terms of run-time; simply you do 6 times less simple-loops; the simple-loop is simply a for-loop).

The simple-loop itself, is just an iterative algorithm around UEqn.H and pEqn.H.

So if you do 6 simple loops, you should get the same result compared to 1 simple loop and your additional 6 loops. Both do make 6 loops around UEqn.H and pEqn.H. As we don´t have any time-information here, its simply as follows:

1. Construct U-Matrix and update U (if momentumPredictor = true)
2. Construct p-Matrix and solve it, correct U based on new fluxes

With the new p-U fields, you start again on the first point and construct the momentum matrix. Here, and this is the main point, it does not matter if you go one SIMPLE iteration or one of your iterations. Hence, each iteration in your for-loop will differ except if you reached the steady state solution. Your for-loop will not gain anything in your case.
Furthermore, the relaxation factors can be set by using the <Function1> Class. So you have more freedom to set a time-dependent or what ever-dependent relaxation factor. Despite of that you can use Doxygen to search for the function you are searching:

https://www.openfoam.com/documentati...1fvMatrix.html
https://www.openfoam.com/documentati...ebf646011c212e

As you can see, the relax() function is defined in the fvMatrix class. We also have the <fvMatrix>.relax(const scalar) function. You simply can add any number here which is used as relaxation factor rather than the one provided and stored in the class (read from fvSolutions).

Normally I would totally agree with you saying that e.g 5 simple loops are the same as 1 simple loop with an inner loop doing 5 times the same but the results differ. I have uploaded two simulations where I compare my results. What I`m doing is simulating the velocity field over a porous media (left hand side of both images). On the right hand side you can see the results (turquoise line represent the normal simpleLoop and the brown line the inner looping method).


I have done 2 simulations https://we.tl/t-N3pm5wCFzd:


1) init5000in1.png for initial flow field with 3 steps:
500 time steps with Urelax of 0.1, than
2000 time steps with Urelax of 0.001, than
2500 time steps with Urelax of 0.00001
For normal simple loop I´m doing this manually by changing Urealx after the specific number of time steps until I get the final result after 5000 steps. For the modified simple loop with inner for (or while) loop I have automated this with two arrays: One with the time steps in it rUTS(500,2000,2500) and one with the relaxation factors rU(0.1,0.001,0.00001). This way I can build two loops around UEqn and pEqn to get the 5000 steps in one single simple loop. As you can see, the results differ.


2) With the manually generated initial solution of 1) I start calculating my "final solution". Therefore I only have one time step number (5 here) and one relaxation factor (0.001 here).
As I expect for both results, with the increase of relaxation from 0.00001 (final init relaxation) to 0.001 (relaxation for final solution), the velocity shows more oscillations (which can be seen here). But you can also see the same behavior like in 1)... 5 simple loops are not the same as 1 simple loop with 5 inner for loops.


I´m not sure whats the reason for it. Maybe the algorithm uses some "backup" or "presaved" data as another inital guess or something like this. That would explain the difference for me because for 1 loop with 5 inner iterations I would not update this "guess" or something like this.



For the point of setting relaxation... I did this all (thanks for your help) with UEqn.realx(rU[rUi]). At the moment I get the rU[rUi] by a new defined dictionary in fvSolution and set it over the inner loops afterwrds. This works quite good for what I would like to do but is there also a way to get the normal defined relaxation value of fvOptions and overwrite it? Than I could leave my UEqn.realx() as default with the advantage, that I can also leave the standard relaxation factor unchanged and change it only if needed.


Here is some code to clarify the inner loop part and the way, I read my fvSolution:


simpleFoam:

while (simple.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        forAll(rU,rUi)
        {
            for(int i = 0; i < rUTS[rUi]; i++)
            {
                // Pressure-velocity SIMPLE corrector
                {
                    #include "UEqn.H"
                    #include "pEqn.H"
                }
 
                laminarTransport.correct();
                turbulence->correct();
            }

         }

        runTime.write();

        runTime.printExecutionTime(Info);
    }

createdFields:

IOdictionary fvSolution // Reading different relaxation factors
(
    IOobject
    (
        "fvSolution",
        runTime.system(),
        runTime,
        IOobject::MUST_READ
    )
);
dictionary& fvSol = fvSolution.subDict("relaxationFactors").subDict("equations");

label NRelaxU(readLabel(fvSol.lookup("NRelaxU")));

List<scalar> rU(NRelaxU); // List of relaxation parameter
List<label> rUTS(NRelaxU); // List of relaxation time steps

forAll(rU,rUi)
{
    rU[rUi] = fvSol.subDict("relaxU").lookupOrDefault<scalar>(("rU" + name(rUi)),0.1);
    rUTS[rUi] = fvSol.subDict("relaxU").lookupOrDefault<label>(("rUTS" + name(rUi)),1);
}

fvSolution:

relaxationFactors
{
    fields
    {
        p               0.1;
    }
    equations
    {
        U               0.1; //(0.1, 0.001, 0.000001) for init flow field
        k               0.1;
        epsilon         0.1;
        NRelaxU         3;
        relaxU
        {
          rU0           0.1;
          rU1           0.001;
          rU2           0.000001;
          rUTS0         500;
          rUTS1         2000;
          rUTS2         2500;
        }
    }
}

Thanks for any help


Chris

[Tobi]

No idea why you does make it so complicated:

relaxationFactors                                                               
{                                                                               
    equations                                                                   
    {                                                                           
        U                                                                       
        {                                                                       
            type    table;                                                      
            values                                                              
            (                                                                   
                (0  0.1)                                                        
                (10 0.2)                                                        
                (50 0.3)                                                        
            );                                                                  
        }

Here you have a linear blending between time 0 and 10 (from 0.1 to 0.2) and from 10 to 50 (from 0.2 to 0.3).
And, I don't expect that there are any differences between the normal SIMPLE loop and the SIMPLE loop with additional for-loop. Its the same. If you need a proof, I will give it to you.

[Pyrokrates]

Hey Tobi,

thanks you for the answer. Using this kind of linear blending result in following error:

Attempt to return dictionary entry as a primitive

I have tested openfoam-v2006, openfoam-v2106, and openfoam-8.


Problem is that I can only define the changing relaxation once per run... What I need is the change of relaxtion in each time step (I know for simpleFoam there is no real time, its more about iterations). But for the part I have programmed for my concentration I need time steps.

This is why I would like to change the relaxation every calculation step (increase Urealx to let U change local in my mesh but reduce it in the same step again to get correct pressure drop and lower fluctuations inside my velocity for the concentration calculation).




A proof for the same results of simpleFoam with and without for loop would be great because I´m confused about my results. I thought the results should be the same but as you can see in my uploaded images, it differs. Maybe its the porousMedia... I will do some testing without porous media and let you know if there is a change...


Thanks for all the time and help until now


Chris

[Pyrokrates]

Hey,


I have made some comparisons and need to correct/clarify something I wrote not correctly in my previous posts...


The first thing is, for my calculations I`m using porousSimpleFoam which I have modified with my concentration models, not the simpleFoam... I thought the two solvers should not differ much only for the part of using porous media but now it looks like they do.

I have compared four simulations doing 100 timesteps with Urelax 0.1 in a curved pipe:
(Link to the plots: https://wetransfer.com/downloads/113...1113935/8e36d6)
1) simpleFoam [red line]
2) simpleFoam (with implemented inner loop, 10 time steps with each 10 inner loops) [pink line]
3) porousSimpleFoam (without porous media) [blue line]
4) porousSimpleFoam (without porous media, with implemented inner loop, 10 time steps with each 10 inner loops) [pink line]

The results are following:
1) and 2) differ from each other.
1) differ from 3) ,which confuses me a little bit.
3) and 4) also differ.


Maybe this helps a little to solve the problem.

Thanks



Chris

[Tobi]

So here is the proof that both SIMPLE and SIMPLE + Inner-Loop are identical:



I simply made a copy of simpleFoam and changed the source file accordingly:

    Info<< "\nStarting time loop\n" << endl;                                    
                                                                                
    label iter = 0;                                                             
                                                                                
    while (simple.loop())                                                       
    {                                                                           
        Info<< "Time = " << runTime.timeName() << nl << endl;                   
                                                                                
        for (int i = 0; i<4; ++i)                                               
        {                                                                       
            // --- Pressure-velocity SIMPLE corrector                           
            {                                                                   
                #include "UEqn.H"                                               
                #include "pEqn.H"                                               
            }                                                                   
                                                                                
            laminarTransport.correct();                                         
            turbulence->correct();                                              
                                                                                
            ++iter;                                                             
        }                                                                       
                                                                                
        runTime.write();                                                        
                                                                                
        runTime.printExecutionTime(Info);                                       
                                                                                
        if (iter == 100)                                                        
        {                                                                       
            Info<< "We made 100 Iterations" << endl;                            
            runTime.writeAndEnd();                                              
        }                                                                       
    }                                                                           
                                                                                
    Info<< "End\n" << endl;                                                     
                                                                                
    return 0;

In the cases I set the end-time to 100. However, the modified case will stop, after the iter guy has 100 iterations done. Hence, both, the simpleFoam and simpleFoamInnerIteration solver give the same results (as depicted in the picture).

SimpleFoam and porousSimpleFoam might give you slightly different results. I never checked the seceond one but there is definitely some other implementation which you have to take into consideration.

You analysis seems to modify the relaxation factors and hence, it is obvious that you don't get the same guys. For me it is still unclear what you want to achieve as it does not make any sense to me. Nevertheless, as I showed and Tobermony also backed,

• Using SIMPLE or SIMPLE + inner iterations results in same quantities
• What you are doing with the relaxation - I don't know

Cheers.

[Tobi]

Hey all,


I got a private message from Chris in which he clarifies his ideas and what he want to do. I will summarize and give feedback now:

• He wants to calculate a filter which gets blocked by parcels or something else
• Therefore, he wants to calculate the velocity field in prior, re-evaluate the parcels/tracer gas (should be time-dependent) and then, based on the blocking of the filter, he wants to re-evaluate the velocity field
• That's why he is using some relaxation technique which should give him better control to re-calculate the velocity field

Okay, so far so good. The main problem is, Chris, that you don't understand SIMPLE (steady-state solver) in your topic. The SIMPLE algorithm needs a few hundred iterations (depending on the case and relaxation factors) to converge. Two-to-ten iterations might still have some bad p-U correlation (or even nonsense data), which don't make sense at all. Hence, what you should do is:

• Calculate the steady-state velocity-pressure fields (no-time dependence); so probably a few hundred iterations (the first time, afterward maybe 10 to 20, depends on your flow-field and influence of the parcels back to the p-U field)
• Calculate your parcel/tracer gas distribution (based on a time-based solution; so ddt is not zero in this case)
• Repeat the velocity-pressure field calculation; however, here, you need to add source terms to the U or p equations which determines somehow the blocking of the filter (out of the box I don't have any idea doing that without having a Lagrangian solver available)

This would be probably the way you should go. Despite of that, I honestly suggest you to use a fully coupled solver which means, parcels interact with the fluid and vice versa. As this problem is a transient one, I would not go with your approach - maybe depending on the blocking you might get acceptable results; depends mainly how often you want to recalculate the p-U field. Nevertheless, for me it looks more like a fully coupled Lagrangian-Euler problem in a transient manner, rather than something else. But of course, you can also try your approach as you should be the one who is more familiar with the physics inside your topic. But again, changing the relaxation factors after, e.g., 5 iterations in some sense will just give you a guess of the field (if not wrong at all). Hence, I would never rely on the results you will achieve with this approach; only if you split it really: p-U solving to convergence (steady-state), do some time-steps with your traces, recalculate p-U field to convergence, and so on. This is an approach similar we do it with radiation models (but no idea if this will hold for your approach, as the blocking will definitely effect the flow and depending if it is blocked fast or not, the p-U recalculation has to be performed more or less often.


The last question Chris mentioned:

• His results of simpleFoam and my derivative given above done give the same results
• Here, you probably compare not the same times. You cannot check time-step 5 of both solvers. You need to compare, e.g., time-step 100 with time-step 25 to have the same number of iterations performed. Doing so, you end up with the results I presented you.

[Tobi]

For comparison reasons and transparency:

https://holzmann-cfd.com/forum/simpl...eration.tar.gz



The solver is not attached as you simply can built it by copy simpleFoam and change the source file accordingly (as given two posts above).

[Tobermory]

Originally Posted by Tobi

As this problem is a transient one, I would not go with your approach - maybe depending on the blocking you might get acceptable results; depends mainly how often you want to recalculate the p-U field ...

if you split it really: p-U solving to convergence (steady-state), do some time-steps with your traces, recalculate p-U field to convergence, and so on. This is an approach similar we do it with radiation models (but no idea if this will hold for your approach, as the blocking will definitely effect the flow and depending if it is blocked fast or not, the p-U recalculation has to be performed more or less often.

Great answer Tobi - and the timescales are key. Another analogy to consider is transient conjugate heat transfer to an object with a big thermal inertia, where the timescale for the temperature change is very slow, whilst the fluid flow time scale is often very short. For this case, it is too expensive just to run the whole simulation as a transient. Instead, it makes sense to run the thermal analysis as a transient, with a frozen velocity field (this is cheap!). The velocity field (expensive to update) is solved for in a steady simulation, and only updated at a suitable interval.

For your case, if you are calculating a filter blockage/deposition rate in your simulation, and if the rate of deposition is slow compared to the flow timescales, then the above pseudo-steady approach for the flow field will possibly give you a reasonable estimate. Otherwise, you need a fully coupled transient as per Tobi's suggestion.

PS - I forgot to add - Tobi is absolutely right - a part converged SIMPLE solution could contain absolute junk ... and the rate of convergence of the field is totally case dependent, so there's no real way of taking advantage of a SIMPLE solution after only a handful of iterations; you probably need 100's, in which case just go to full convergence, and you are heading towards the 1st approach above.

[Tobi]

Totally agree Tobermory, thank you for bringing up another field in which we use a similar approach.


However, I want to update this talk in terms of differences, of the given implementation I presented in prior posts, if we use under-relaxation for p field.
So lets reconsider the code snippet I presented above:

while (simple.loop())                                                       
    {                                                                           
        Info<< "Time = " << runTime.timeName() << nl << endl;                   
                                                                                
        for (int i = 0; i<4; ++i)                                               
        {                                                                       
            // --- Pressure-velocity SIMPLE corrector                           
            {                                                                                                          
                #include "UEqn.H"                                                                                        
                #include "pEqn.H"                                               
            }                                                                   
                                                                                
            laminarTransport.correct();                                         
            turbulence->correct();                                              
                                                                                
            ++iter;                                                             
        }                                                                       
                                                                                
        runTime.write();                                                        
                                                                                
        runTime.printExecutionTime(Info);                                       
                                                                                
        if (iter == 100)                                                        
        {                                                                       
            Info<< "We made 100 Iterations" << endl;                            
            runTime.writeAndEnd();                                              
        }                                                                       
    }

This will give us the same result as simpleFoam, if we use no-under-relaxation for p (as given in my cases attached).
However, Chris mentioned, that using under-relaxation for the p field, we obtain different results (at least for experienced user it is obvious why; hence, I want to clarify it).

The field relaxation is simply:


Therefore, we update the new iteration (n+1) with the old iteration and the new value. Having no relaxation means , and hence, the new evaluated value is used completely. Despite of that, any value of less than 1, takes the old p field into account. And now the differences:

• simpleFoam will update the old p field each iteration (its inside the <Time> class)
• the modified simpleFoam solver will do the same and hence, if we construct a for loop around UEqn.H and pEqn.H, we need to update the old p values in each inner-loop explicitly.

Therefore, we need to have the following code:

while (simple.loop())                                                       
    {                                                                           
        Info<< "Time = " << runTime.timeName() << nl << endl;                   
                                                                                
        for (int i = 0; i<4; ++i)                                               
        {                                                                       
            // --- Pressure-velocity SIMPLE corrector                           
            {                                                                                                          
                #include "UEqn.H"
    
               // We need to store the previous iteration (overwriting the p-old)
                p.storePrevIter();                                
                                                   
                #include "pEqn.H"                                               
            }                                                                   
                                                                                
            laminarTransport.correct();                                         
            turbulence->correct();                                              
                                                                                
            ++iter;                                                             
        }                                                                       
                                                                                
        runTime.write();                                                        
                                                                                
        runTime.printExecutionTime(Info);                                       
                                                                                
        if (iter == 100)                                                        
        {                                                                       
            Info<< "We made 100 Iterations" << endl;                            
            runTime.writeAndEnd();                                              
        }                                                                       
    }

Regarding the velocity field U: We don't need to update the U field, because we are using matrix relaxation here and furthermore the U field gets updated in the pEqn.H (after we have the pressure field).

[Pyrokrates]

Hey Tobi,


thank you very much for the help. The part with the storing of p was exactly what I supposed when I said, that it looks like inner looping does not use some pre saved data correctly.


Thanks and have a nice day


Chris