Hi everyone, with the code I am working on, I have observed that when the multigrid is turned on, the residual reaches very lower values (compared to the same case run without multigrid) . I am a bit surprised about this result, as I expected an influence on the convergence speed ONLY, and not on the final solution. Any explanation? thanks

The convergence rate on the single fine mesh is not linear in the count of Work Units (fine-mesh equivalent iterations). Initially, there may be some bouncing around if starting from initial conditions that are far from the final steady-state. Then there is a relatively rapid drop in residuals as the small-wavelength errors are smoothed out. Then the (slow) asymptotic convergence rate on that mesh is reached. If there are no hindrances, the residuals will eventually be driven to machine zero (about 10^-11) in 64-bit precision. So, one of two things can give you a higher level of terminal residual norm than that. Either your code stops at a prescribed total number of iterations, before the norm can reach machine zero, or there is a subtle convergence hindrance consisting of a low-level long-range interaction between nonlinear flow features or boundary conditions that are a good distance apart. The relaxation on the fine mesh will have little to no damping effect on such interactions, unless you analyze and modify the relaxation scheme. The multigrid iteration on the other hand either reaches a lower residual norm level before any prescribed total number of WUs, or else if there is no such prescribed limit, the coarse mesh relaxations are more effective at damping the long-range interactions. Hence, the multigrid iteration will drive the residual norm lower than will the single grid, before the convergence flattens out. In some cases, the multigrid iteration will reach a lower final residual level than the single grid one, but will still not be able to reach machine zero, because of said hindrances.

Hi, thank you very much for the detailed answer, and sorry for the delay of this reply. As far as I understand, you say that the behaviour I indicated in my first message may be normal, if interactions with BC or nonlinear flow features introduce some "subtle convergence hidrance". This error should are easier damped on the coarse mesh....my question is then why?

- The nonlinear interactions tend to pruduce higher and higher wavenumbers,

on which the the coarse grid has less impact

- When the bc introduce an error, in most of the cases they consist of

spurious localised wave packets. Again, they are high wavenumbers.

In addition, don't the multigrid introduce any error during interpolation and prolungation ? Many thanks

When I used the term nonlinear, I was thinking of flow features such as shockwave location. Perhaps it applies to flow features both linear and nonlinear. You are right that local nonlinear interactions and also local boundary condition errors produce high-frequency errors (local => small-wavelength), which the fine-mesh relaxation should be better at smoothing, with a good choice of smoothing scheme. I was rather pointing at interactions between flow features and/or boundary conditions that are well separated spatially, so that the errors involved are large-wavelength (smooth and global). Very often, the relaxation on the fine mesh just propagates these low-frequency long-wavelength errors, and they bounce around from corner to corner of the domain without getting mutually reconciled. For example, an upstream shockwave location and a downstream reattachment point of a separated boundary layer may influence each other and both may oscillate without ever settling down. This can of course reflect some physical instability, but in some cases it can be the result of the weakness of the relaxation scheme in reconciling long-distance interactions because of the local nature of the relaxation scheme. If instead, you solved the entire domain by implicit time-marching with a direct solver, you would not see such bouncing around. Normally, Dirichlet boundary conditions should act as absorbers of these long-wavlength errors, because Dirichlet bcs enforce zero error, but I think the interior relaxation scheme is still capable of reflecting such long-wavelength error waves at the layer of cells next to the Dirichlet boundary. I have observed behavior such as you described, and I had arrived at this intuitive explanation of what was involved. However, one would need to do some matrix stability/eigenvalue analysis on the relaxation scheme linearized about the approximate solution to back up such an explanation.

Your explanation reminds me of an experiment I did some time ago. I was actually doing stability tests perturbing a given parallel flow, with localised disturbances. In some cases for which I was expecting the full expulsion of the perturbation, I noticed some global modes with low wavenumbers to develop. The impression was that they originated at the boundary, or, better, they were not physical but arose because of the interaction with the boundary. A famous paper review on stability, Huerre & Monkewits, reports this possibility as well, although, if I remember well, they provide an explanation based on the elliptic character of the Poisson equation in a pressure-correction model (whereas my code was a standard compressible time-marching). As I said, I am not sure 100% and have not time ;-( to look for it again, but assuming what I remember is right, there is no reason to believe that "if instead you solved the entire domain by an implicit time-marching with a direct solver you would not see such bouncing around". If you have experienced this before, well, that's very interesting, but ... are you sure of a fully consistent comparison withe the explicit? For example, was it the same time step? Back to my original question, what I still have difficulty to understand is why, even considering such long-range low wavenumbers errors, the relaxation scheme on a coarse grid should provide a better convegence rate. Its damping properties should be regulated by the CFL condition, that is, even when you change the grid, the magnitude of the amplification factor (the complex Fourier symbol \dot{z}) shouldn't change if the CFL remains constant. Are my considerations correct? Can you also tell me if the interpolation and prolungation introduce errors? Thanks much D

I did not rerun those cases with a fully implicit method, sorry if I implied otherwise. When I observed that phenomenon, I was using a line relaxation algorithm on the steady-state equations with and without multigrid. No time terms, though the relaxation algorithm could be recast as effectively introducing a time derivative of the pressure in the streamwise momentum equation.

You are right, it is possible that the low-level unsteadiness might be observed even with an implicit time-marching scheme. I do not know what case you are solving, and I should not be making sweeping statements. I assume from your posts that you are seeking a steady-state solution to a nonlinear boundary value problem (BVP) by recasting it as an initial-boundary value problem (IBVP) by introducing time derivatives arising from conservation/balance laws and using a time-marching solver with and without multigrid. This effectively uses time-marching as a relaxation path to the solution of the BVP.

When we examine the formulation and solution process for errors, the first questions that arise are existence and uniqueness of a solution to the BVP. Next is the question of stability of the solution, so that we can expect a relaxation method or solution of the IBVP to converge to the solution of the BVP without residual unsteadiness or divergence to some other stable, periodic or chaotic solution. I assumed that all these were not a concern.

Next, there is the possibility of errors introduced by the discretization. The discretization may allow spurious solution components to persist in the time marching and just bounce around. The Fourier analysis of the interior scheme would reveal any such weakness of the interior scheme for constant-coefficient linear BVPs. A matrix stability analysis would account for both interior scheme and boundary conditions. Again, such a stability analysis is only practical for linear BVPs. I am assuming that since your coarse mesh relaxation is able to drive the residuals down to obtain a steady-state, this sort of discretization error is also not a concern.

Next, we must examine the relaxation process for the inversion of the discretized nonlinear equations. Now, I am going to use the terms "implicit" and "explicit" in their original sense to refer to the degree of coupling of unknowns, whether or not time-marching is involved. Generally speaking, there is a rough correlation between degree-of-implicitness of the inversion method and the "robustness" of the same. The robustness translates to rapidity of convergence per iteration and size of solution space set (i.e., family of BVPs) for which the iteration process will converge to the desired solution. The switch from direct inversion (of linear equations) to decreasingly implicit relaxation methods brings increasing possibilities for errors to bounce around the domain and persist. I have not read the Huerre & Monkewits paper, but note that the pressure-correction method already loses some amount of implicitness by solving the equations in a segregated manner, though proponents of these methods would claim that implicit treatment is reserved for the pressure poisson equation which supposedly gathers all the elliptic influence.

My guess was that it is in this area that you are seeing some weakness of the relaxation. That the relaxation on the fine mesh alone, allows for some smooth errors to bounce around the domain, because of the local less-implicit nature of the fine-grid smoothing. When you couple a less-implicit relaxation method with multigrid acceleration, the coarse grid relaxation effectively acts to make the overall relaxation more implicit in character per fine mesh work unit. This is why you can expect better convergence rate and better damping of long-wavelength errors using multigrid. I am not sure if your assumptions are right regarding the relation between damping properties and the Courant number. The degree of damping is better for shorter wavelengths RELATIVE to the mesh spacing, and if you maintain the same Courant number on the coarse mesh as you do on the fine mesh, the coarse mesh relaxation should provide damping that is as effective on double the absolute wavelength (as measured on the fine mesh) as the fine mesh relaxation is on a given wavelength. Note that the flow features and boundaries that may be interacting to allow the errors to bounce back and forth, are separated by fewer coarse mesh intervals than they are by fine mesh intervals. This leads to better coupling (more implicitness) between them on the coarse mesh than on the fine mesh. Yes, restriction and prolongation introduce errors, too, but I expect that these are local short-wavelength errors that are quickly smoothed out.