Dual-time Stepping in 3D Compressible Solver

Crank-Shaft · November 4, 2014, 03:29

Hello CFD-Onliners,

I am currently working on the development and extension of an in-house research CFD code. This is a 3D compressible flow solver computing on structured, multiblock meshes with Euler Equations in curvilinear coordinates.

There is an urgent need to try and develop and implement a dual-time stepping algorithm due to the large domain sizes I work with and a vast range of different cell sizes. Since the explicit multistage Runge-Kutta methods rely on the global time-step, the time-step sizes are impractically short and this leads to unfeasible computational wall-clock times. Hence, I've been trying to introduce a Dual-time stepping algorithm.

Currently, this is my general understanding (referring to the notation in the attachment) -
a) I take 2 different CFL values and the larger one (physical) is used to set the global time-step however, the smaller CFL (pseudo) is used with the local time-steps for each cells and then the minima for each block is used.
b) At nth step (physical) set p solution = n solution. Then the fluxes for all the conserved variables are also taken at nth time-step.
c) A loop is initiated from the main solver program and this is checking for MaxResidual<1e-6 for each domain.
d) The subroutine for the inner iterations is called for each block in the domain. This computes the residual as shown in the attachment and the time-step ratio. Then the solution at (p+1) is updated for all the variables. I am using the residual and the flux at pth state, prior to updating for the p+1.
e) This is then returned to Main solver which enters or exits the 'while' loop depending on the MaxResidual condition mentioned.

However, the main problems I encountered are as follows -
a) The residual value as defined in the attachment doesn't seem to be decreasing at all. They tend to build up. Would this essentially require a different definition such as (e_new-e_old)/e_new or would you suggest normalising?
b) I haven't been able to use Physical CFL > 1.0 for anything significant and it simply leads to crashes and non-physical values. I've encountered significantly higher values in literature for projects using Pseudo-time, semi-implicit or dual-time stepping algorithms.
c) This still appears to explicit time-stepping even in the pseudo-time. Other than the time-step ratio being used, I don't quite understanding how the acceleration is supposed to take place. Is it purely due to the local time-stepping, which effectively iterates repeatedly for the blocks which are slow to converge while the rest essentially wait, prior to being advanced in physical time to (n+1)th.

Here are two approaches I've tried -

1) When the fluxes are taken at nth time and then the pseudo-time iterations are computed, while holding that constant, the residuals seem to fall quite rapidly. However, this only seems to represent an explicit scheme and all it's doing is simply using a larger time-step (Due to larger physical CFL user-specified value). Ultimately, if it's run for a long time, it leads to non-physical values and some cell values just deteriorate rapidly to Infinity.

2) Another attempt was to ensure that for every p-> (p+1) update, the convective and diffusive fluxes are updated. Unfortunately, with the residual definition given (IMAGE) it just tends to grow until numerical limits are reached.

It would be greatly appreciated if someone could help me understand this better, or identify any flaws in the concept and code structure described above.

Looking forward to those that have either worked with research CFD codes, with Su2 or OpenFOAM in the past or those that are familiar with the numerical methods being discussed above. Most of the content that I've found from existing literature only discuss more elaborate approaches such as multi-grid methods, semi-implicit methods or those that require conditioning matrices...

Thanks in advance!

sbaffini · November 4, 2014, 14:39

How do you expect a compressible flow solver to work for CFL larger than 1 if you don't use preconditioning or any other mean to filter out the relative pressure stiffness? Do you have a reference for this?

Crank-Shaft · November 4, 2014, 22:00

Quote:

Originally Posted by sbaffini

How do you expect a compressible flow solver to work for CFL larger than 1 if you don't use preconditioning or any other mean to filter out the relative pressure stiffness? Do you have a reference for this?

Thank you for your response. I have limited knowledge of the actual impact of trying to implement Dual Time-stepping in an existing compressible solver, especially with the stiffness effects you mentioned.

I've basically followed a combination of two different sources -
panther.umd.edu/mediawiki/images/a/a8/Dual_Time.pdf
http://congress.cimne.com/eccomas/ec...00/pdf/781.pdf

recent -
http://www.iscfdc.co.il/sites/default/files/eznss.pdf

My initial exploration led to multigrid methods, semi-implicit methods and preconditioning however, I was just searching for something which would be simpler to implement.

Would you be able to explain how the pre-conditioning would assist with the pressure-induced stiffness?

Currently, the main source of my confusion arise from the Flux term being decomposed into two terms, one requiring the solution at pth state and the other requiring the Flux Jacobian (James Baeder PDF attached). I've had difficulty understanding this and finding sources explaining this.

Another issue is the calculation of the residual itself. I have implemented it in the same way as J. Baeder above, however, this being defined as the residual is somewhat flawed since the value continuously rises uncontrollably. Hence, I might have to switch to a norm of the density or energy and then apply this as a convergence threshold?

Look forward to everyone's comments and appreciate your time.

Crank-Shaft · November 5, 2014, 01:04

I would also like to add that the compressible solver is essentially being run at low-Mach number (0.1-0.3), with appropriate corrective algorithms and for High-Reynolds number simulations.

Hence, the stiffness in the pressure field might not be a source of serious problems...

Please let us know your thoughts and if you have some suggestions.

FMDenaro · November 5, 2014, 06:35

Quote:

Originally Posted by Crank-Shaft

I would also like to add that the compressible solver is essentially being run at low-Mach number (0.1-0.3), with appropriate corrective algorithms and for High-Reynolds number simulations.

Hence, the stiffness in the pressure field might not be a source of serious problems...

Please let us know your thoughts and if you have some suggestions.

Actually, it is at low Mach that you have to face pressure stiffnes...

Crank-Shaft · November 5, 2014, 20:47

Thanks for your response Filippo. This is the first time I'm working in actual CFD development so I am unaware of some of the pitfalls of the algorithms and numerical schemes I am working with.

Do you have any suggestions on the following -

a)The actual definition of the residual term using the flux (current pth step) and the 2nd Order BDF term (with n-1, n and pth terms). This seems to keep increasing even for the simplest of domains (homogeneous cube with side-wall).

b) The flux terms in the above are being calculated at pth state, however, in the original formulation, the implicit nature called for (p+1)th state. How does one find this value? Do we need to calculate another Flux Jacobian, relating the Flux to the gradient of the conserved variable (as shown in Baeder PDF above $(F_{Q}\Delta{Q})^{p}_{x}$ or $\frac{\partial{F}}{\partial{Q}}\Delta{Q}$ )?

c) To stop the pseudo or inner-iterations, I've decided to use the above residual with the 2nd order BDF and flux terms however, I've also thought about using another definition with the L2-norm of the density or energy or $\epsilon_{resid}=\left| \frac{Q^{p}-Q^{p-1}}{Q^{p-1}}\right|$ or just after the new pseudo-time step solution is calculated, use $\epsilon_{resid}=\left| \frac{Q^{p+1}-Q^{p}}{Q^{p}}\right|$ and wait until $\epsilon_{resid}<10^{-6}$ ?

Look forward to your suggestions.

FMDenaro · November 6, 2014, 04:15

if you have no experience in developing CFD codes, your goal is quite full of troubles...
What kind of textbook are you using?

Crank-Shaft · November 6, 2014, 05:34

Quote:

Originally Posted by FMDenaro

if you have no experience in developing CFD codes, your goal is quite full of troubles...
What kind of textbook are you using?

I have used CFD for years now, but have always used commercial packages and mainly focused on applied aerodynamics problems. So this is one of the first times actually developing functions and adding features to a research CFD code.

In terms of texbooks, I've used Turbulence modelling titles from S. Pope and D.C. Wilcox. In the past I've followed "An Introduction to Computational Fluid Dynamics: A Finite Volume Approach" by Versteeg and Malalasekara.

I currently use "Implicit Large Eddy Simulation: Computing Turbulent Fluid Dynamics" by Grinstein, Margolin and Rider. Also I sometimes refer to "Riemann Solvers and Numerical Methods for Fluid Dynamics" by E. Toro.

So far the project's gone well but I am just finding this phase particularly challenging. Your suggestions are greatly appreciated.

FMDenaro · November 6, 2014, 08:00

have also a look to
http://www.cambridge.org/us/academic...cs-2nd-edition

sbaffini · November 6, 2014, 14:13

You might want to consider the following:

http://web.stanford.edu/group/uq/pdf...s/jcp_ib07.pdf

Besides the cartesian factorization and the immersed boundary method, the solver is what you are probably looking for.

Crank-Shaft · November 11, 2014, 08:19

Thank you for all your advice, everyone.

Since I was having problems with the residual definition, I've explored some of the formulation commonly used in commercial codes. I've added a new definition for the iterative residual, for both the density and energy conserved variables ( $\phi$ ) -
$\epsilon_{resid}=\left| \frac{Q^{p+1}-Q^{p}}{Q^{p}}\right|$
which is calculated immediately after the

level solution is updated in the sub-iteration loop. They seem to be working fine.

However, the problems I currently face are related to the unsteady derivative terms in the original equation, which were discretised with the 2nd order backward difference scheme. This was mentioned as being the main residual for the Newton-like iterations, in the J. Baeder article discussed earlier.

$\left| \frac{3\phi^{P}-4\phi^{n}+\phi^{n-1}}{2\Delta{t}} +F_{x} \right|$

I realised that for this value to drop by three orders of magnitude $(10^{-3})$ , it should be normalised for it to even make sense.
For anyone with past experience, is it enough to normalise this as the other energy and density conserved variables, or is something different required since it has multiple terms including the flux terms?

Thank you

sbaffini · November 11, 2014, 15:09

Why would you want to compute residuals on a time derivative? It is the time derivative itself which, for steady computations, can be interpreted as a residual. Otherwise it is just another term of the equations.

More generally, the residual of an equation is always defined as the difference (lhs-rhs). If you have d/dt in any of them you should not use them as residuals. In the DTS you want to converge within the single time step, for which the d/dt term is already discretized. Thus it doesn't go in the residual definition in the way you wrote it.

P.S. The links to the papers are all dead, thus i cannot provide any more useful comment to your procedure.

Edit: I just proved myself wrong. Now only one (the first) is dead. I'm reading the last one (the second one seems useless to me).

sbaffini · November 11, 2014, 15:44

From the code manual doesn't seem to exist any preconditioning in the code. Thus, according to my knowledge, you can't solve for incompressible flows.

From what i understand, the DTS procedure is implemented in the code in order to increase the time accuracy beyond first order, because of the way the basic equations are solved.

Crank-Shaft · November 11, 2014, 16:52

Try this -

http://panther.umd.edu/mediawiki/images/a/a8/Dual_Time.pdf

Thanks

FMDenaro · November 12, 2014, 07:09

I suppose it is a misunderstanding about the term "residual" ...

sbaffini · November 12, 2014, 14:06

Premised that i don't fully get the linearization explained in the document, i think that you only need an useful normalization for the residuals.

The density based solver in Fluent does the following:

1) Absolute residual: RMS of the physical time derivative for the given conserved variable: R=SQRT(Sum_c (dW/dt)^2_c/Ncells)

2) Globally scaled residual: the residual above is divided by its greatest value in the first 5 iterations: Rg = R/MAX(R,iter=1,..,5)

3) Locally scaled residual: R=SQRT(Sum_c (dW/dt*V)^2_c/Ncells)/(Wmax-Wmin)_domain

where V is the cell volume, W is the generic conserved variable, Ncells is the total number of cells, Sum_c is a sum over the cells and ()_c refers to a generic cell.

Your ABSOLUTE residual (as defined in the document) seems good to me. You just need to select the appropriate scaling

November 5, 2014, 01:04		#4
Crank-Shaft Senior Member Ovi Join Date: Oct 2012 Location: Sydney, Australia Posts: 166 Rep Power: 13	I would also like to add that the compressible solver is essentially being run at low-Mach number (0.1-0.3), with appropriate corrective algorithms and for High-Reynolds number simulations. Hence, the stiffness in the pressure field might not be a source of serious problems... Please let us know your thoughts and if you have some suggestions. __________________ -- Mechanical Engineering Sydney, Australia

November 5, 2014, 20:47		#6
Crank-Shaft Senior Member Ovi Join Date: Oct 2012 Location: Sydney, Australia Posts: 166 Rep Power: 13	Thanks for your response Filippo. This is the first time I'm working in actual CFD development so I am unaware of some of the pitfalls of the algorithms and numerical schemes I am working with. Do you have any suggestions on the following - a)The actual definition of the residual term using the flux (current pth step) and the 2nd Order BDF term (with n-1, n and pth terms). This seems to keep increasing even for the simplest of domains (homogeneous cube with side-wall). b) The flux terms in the above are being calculated at pth state, however, in the original formulation, the implicit nature called for (p+1)th state. How does one find this value? Do we need to calculate another Flux Jacobian, relating the Flux to the gradient of the conserved variable (as shown in Baeder PDF above $(F_{Q}\Delta{Q})^{p}_{x}$ or $\frac{\partial{F}}{\partial{Q}}\Delta{Q}$ )? c) To stop the pseudo or inner-iterations, I've decided to use the above residual with the 2nd order BDF and flux terms however, I've also thought about using another definition with the L2-norm of the density or energy or $\epsilon_{resid}=\left\| \frac{Q^{p}-Q^{p-1}}{Q^{p-1}}\right\|$ or just after the new pseudo-time step solution is calculated, use $\epsilon_{resid}=\left\| \frac{Q^{p+1}-Q^{p}}{Q^{p}}\right\|$ and wait until $\epsilon_{resid}<10^{-6}$ ? Look forward to your suggestions. __________________ -- Mechanical Engineering Sydney, Australia

November 11, 2014, 08:19		#11
Crank-Shaft Senior Member Ovi Join Date: Oct 2012 Location: Sydney, Australia Posts: 166 Rep Power: 13	Thank you for all your advice, everyone. Since I was having problems with the residual definition, I've explored some of the formulation commonly used in commercial codes. I've added a new definition for the iterative residual, for both the density and energy conserved variables ( $\phi$ ) - $\epsilon_{resid}=\left\| \frac{Q^{p+1}-Q^{p}}{Q^{p}}\right\|$ which is calculated immediately after the level solution is updated in the sub-iteration loop. They seem to be working fine. However, the problems I currently face are related to the unsteady derivative terms in the original equation, which were discretised with the 2nd order backward difference scheme. This was mentioned as being the main residual for the Newton-like iterations, in the J. Baeder article discussed earlier. $\left\| \frac{3\phi^{P}-4\phi^{n}+\phi^{n-1}}{2\Delta{t}} +F_{x} \right\|$ I realised that for this value to drop by three orders of magnitude $(10^{-3})$ , it should be normalised for it to even make sense. For anyone with past experience, is it enough to normalise this as the other energy and density conserved variables, or is something different required since it has multiple terms including the flux terms? Thank you __________________ -- Mechanical Engineering Sydney, Australia

November 11, 2014, 16:52		#14
Crank-Shaft Senior Member Ovi Join Date: Oct 2012 Location: Sydney, Australia Posts: 166 Rep Power: 13	Try this - http://panther.umd.edu/mediawiki/images/a/a8/Dual_Time.pdf Thanks __________________ -- Mechanical Engineering Sydney, Australia

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
How to export time series of variables for one point?	mary mor	OpenFOAM Post-Processing	8	July 19, 2017 10:54
plot over time	fferroni	OpenFOAM Post-Processing	7	June 8, 2012 07:56
IcoFoam parallel woes	msrinath80	OpenFOAM Running, Solving & CFD	9	July 22, 2007 02:58
Dual Time stepping	Aditya	Main CFD Forum	4	November 12, 2005 00:38
Dual Time Stepping	CFD Student	Main CFD Forum	20	October 24, 2005 15:57

November 4, 2014, 14:39		#2
sbaffini Senior Member Paolo Lampitella Join Date: Mar 2009 Location: Italy Posts: 2,151 Blog Entries: 29 Rep Power: 39	How do you expect a compressible flow solver to work for CFL larger than 1 if you don't use preconditioning or any other mean to filter out the relative pressure stiffness? Do you have a reference for this?

November 6, 2014, 04:15		#7
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,768 Rep Power: 71	if you have no experience in developing CFD codes, your goal is quite full of troubles... What kind of textbook are you using?

November 6, 2014, 08:00		#9
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,768 Rep Power: 71	have also a look to http://www.cambridge.org/us/academic...cs-2nd-edition

November 6, 2014, 14:13		#10
sbaffini Senior Member Paolo Lampitella Join Date: Mar 2009 Location: Italy Posts: 2,151 Blog Entries: 29 Rep Power: 39	You might want to consider the following: http://web.stanford.edu/group/uq/pdf...s/jcp_ib07.pdf Besides the cartesian factorization and the immersed boundary method, the solver is what you are probably looking for.

November 11, 2014, 15:09		#12
sbaffini Senior Member Paolo Lampitella Join Date: Mar 2009 Location: Italy Posts: 2,151 Blog Entries: 29 Rep Power: 39	Why would you want to compute residuals on a time derivative? It is the time derivative itself which, for steady computations, can be interpreted as a residual. Otherwise it is just another term of the equations. More generally, the residual of an equation is always defined as the difference (lhs-rhs). If you have d/dt in any of them you should not use them as residuals. In the DTS you want to converge within the single time step, for which the d/dt term is already discretized. Thus it doesn't go in the residual definition in the way you wrote it. P.S. The links to the papers are all dead, thus i cannot provide any more useful comment to your procedure. Edit: I just proved myself wrong. Now only one (the first) is dead. I'm reading the last one (the second one seems useless to me).

November 11, 2014, 15:44		#13
sbaffini Senior Member Paolo Lampitella Join Date: Mar 2009 Location: Italy Posts: 2,151 Blog Entries: 29 Rep Power: 39	From the code manual doesn't seem to exist any preconditioning in the code. Thus, according to my knowledge, you can't solve for incompressible flows. From what i understand, the DTS procedure is implemented in the code in order to increase the time accuracy beyond first order, because of the way the basic equations are solved.

November 12, 2014, 07:09		#15
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,768 Rep Power: 71	I suppose it is a misunderstanding about the term "residual" ...

November 12, 2014, 14:06		#16
sbaffini Senior Member Paolo Lampitella Join Date: Mar 2009 Location: Italy Posts: 2,151 Blog Entries: 29 Rep Power: 39	Premised that i don't fully get the linearization explained in the document, i think that you only need an useful normalization for the residuals. The density based solver in Fluent does the following: 1) Absolute residual: RMS of the physical time derivative for the given conserved variable: R=SQRT(Sum_c (dW/dt)^2_c/Ncells) 2) Globally scaled residual: the residual above is divided by its greatest value in the first 5 iterations: Rg = R/MAX(R,iter=1,..,5) 3) Locally scaled residual: R=SQRT(Sum_c (dW/dt*V)^2_c/Ncells)/(Wmax-Wmin)_domain where V is the cell volume, W is the generic conserved variable, Ncells is the total number of cells, Sum_c is a sum over the cells and ()_c refers to a generic cell. Your ABSOLUTE residual (as defined in the document) seems good to me. You just need to select the appropriate scaling