16

The main idea behind multigrid is projection. I try to think about it as follows: Suppose I want to solve a PDE on with a lot of accuracy, so I proceed to discretize the domain (let's say, using finite difference method) on a very fine grid with lots and lots of points. In the end, I setup my system of equations and I'm ready to solve it. I try using my ...


15

This site is possibly not a good place to ask for a detailed explanation with pseudocode (as stated in the FAQ, "If you can imagine an entire book that answers your question, you’re asking too much."), so you might want to start with one of the classical books on this topic (listed below) and come back with specific questions about concrete details you have ...


12

I believe comparing an iterative method (multigrid) to a direct/exact method (Thomas) in terms of exact operation count isn't really meaningful. IIRC, Thomas operation count is $8N$ for any tridiagonal system. The only time I can imagine multigrid conceivably beating that is for a trivial case of having a linear solution, and even then the cost of evaluating ...


11

Where did you come up with interpolating the "error"? (And how do you measure the error?) On the first visit to a finer grid, the entire solution $u$ must be interpolated, ideally using a higher order operator (e.g., postprocessed/reconstructed solution for FEM). This FMG interpolation is $u^h \gets \mathbb{I}_H^h u^H$. (It's okay to use the a normal ...


11

The short answer is that the Thomas algorithm will be faster than any iterative scheme for almost all cases. The exception would perhaps be applying a single iteration of a very simple iterative scheme such as Gauss-Seidel, but this is highly unlikely to give an acceptable solution. Also, this is ignoring parallel processing concerns. Multigrid is an ...


10

In some cases, (F)MG provides an algorithm with optimal properties. For instance, properly tuned FMG can solve some elliptic problems in a small number of "work units", where a work unit is defined to be the computational effort required to express the problem itself - in this case the operations to form the residual $b-Ax$ on the finest grid. This is such ...


10

Given interpolation $I_H^h$ and restriction $I_h^H$ (where restriction is typically $(I_H^h)^T$ for symmetric problems), with fine grid discretized operator $A^h$, there are two common approaches for constructing the coarse grid operator $A^H$. (Petrov-)Galerkin coarse operators This explicitly computes the matrix triple product $$ A^H = I_h^H A^h I_H^h .$$...


8

One simpler explanation - the range of the restriction operator is the coarse grid space, while the range of the interpolation operator is the fine grid space. Unless the two are equal, interpolation + restriction will not result in an identity matrix, because there will always be components of $x$ which are truncated by the restriction operator and lost.


7

I don't directly know your answer as I mainly use FAS instead of correction since I do multigrid for nonlinear problems, but some thoughts you can look into: You're applying a linear correction scheme to a linear problem, so it's not shocking that it does very well. Consider your boundary conditions: make sure you're doing them correctly, and also note that ...


7

I assume you are solving a nice elliptic problems like the Laplacian with smooth coefficients. (You should always check convergence in this friendly setting first.) Confirm that interpolation is as accurate as necessary. There usual rule for vertex-centered MG is that the sum of interpolation and restriction orders should be at least as high as the order of ...


6

If you are using a vertex-centered discretization, then state restriction should be injection rather than the full-weighted residual restriction that it appears you use. That is, replace $I_h^H$ with $\hat I_h^H$ when restricting the state. Using full-weighted restriction for state produces aliasing of high-frequency components of the state which after ...


6

This isn't a multigrid issue, it's a problem formulation issue. Consider the symmetric system $A x = b$ and suppose that $A e = 0$ for some nonzero vector $e$ (the constant when $A$ is the Laplacian with periodic or Neumann boundary conditions). This system is singular and has no solution if you choose $b$ such that $e^T b \ne 0$. Indeed, the exact periodic ...


6

Multigrid doesn't need a Cartesian (rectangular), uniform grid. What it needs is that you can define a fine and a coarse level (possibly recursively, if you want to go from a two-level to a multi-level scheme), and that you can define interpolation operators between these levels. The easiest way to explain this is if you indeed have a Cartesian grid, but you ...


6

There are "Algebraic Multigrid Methods" that can be applied to general systems $Ax=b$.


5

The smoother, restriction and prolongation can be tested independently. In essence, take a periodic domain, put in a particular Fourier mode and see if the code behaves in accordance with Fourier analysis. Also, in Briggs, Henson and McCormick's A Multigrid Tutorial you will find a section called Diagnostic tools (Chapter 4). That should help you with ...


5

Another classic: Wesseling, An Introduction to Multigrid Methods, John Wiley & Sons, 1992. Example codes can be found at MGNet


5

There are two different flavors of smoothing and stability. The spurious oscillations in convection-diffusion problems are not an artifact of the linear solver but an inevitable artifact of the discretization method. Any linear solver for non-symmetric matrices you use will give you the same wiggly answer, or it's wrong. Multi-grid is one such linear solver. ...


5

The answer depends somewhat on the discretization; for example, some boundary integral discretizations result in very well-conditioned matrices, for which a Krylov solver works just fine without having to introduce the multi-level machinery of multigrid. For ill-conditioned matrices arising from finite difference or finite element methods, Krylov solvers ...


4

It's not a question whether the condition holds or not. You choose $R$ so that it's a multiple of $I^T$ because you want the condition to hold. You want the condition to hold because that's the only way you can ensure that the multigrid operator is symmetric and amenable to, say, being a preconditioner for CG.


4

This is a pretty common problem. There are two simpler ways to what you are doing: Enforce a Dirichlet-kind condition on one pressure node by fixing it to zero. This way, your matrix no longer has a null space. The pressure may not have mean value zero, but you can fix this once (instead of once per iteration) after you have the solution of the linear ...


4

Lets say that you have the following grid composed of rectangular elements: Now if you perform your interpolation assuming a normal structured rectangular grid then you will be introducing errors associated with this inaccurate interpolation. In other words when you restrict your residual vector and when you prolong your error vector there will be errors ...


4

Definitely not. To pick one example, the book Multigrid has a plot on page 53 (Figure 2.10) that shows the decrease in the residual as a function of the number of V or W cycles. You would stop cycling when you are satisfied with the size of the residual. The source of your confusion may be because some descriptions only describe a single V-cycle. In some ...


4

Yes, it is normal to have no convergence checks in MG for a few reasons. First, if you use a different number of iterates on each pass, then the MG operator is no longer linear, and you would have to use something like FGMRES as an accelerator which can accommodate a nonlinear preconditioner. Second, FMG is an exact solver (reduces error below discretization ...


4

There are two parts of the answer. First, you don't get the identity because you throw away information during the restriction operation (if you think of a larger mesh than just the three points you consider, then you "forget" the values of ever other node during restriction), and you cannot recover this information during the transpose operation. ...


4

Post-smoothing reduces the high frequency error that is introduced by the coarse grid correction. If you visualize a correction computed on a mesh of size $2h$, then it has no kinks on half of the nodes of the fine mesh of size $h$, and strong kinks on the other half. Post-smoothing distributes this a bit and leads to a coarse grid correction that doesn't ...


3

Thank you to the posters for encouraging me to look for a bug. I found one, a subtle issue related to restriction and interpolation. I am using ghost points to treat the boundaries, and so the first interior point is at index (0,0) (boundary points are stored on negative indices). This can work, but you have to remember that the interior coarse grid points ...


3

There are many libraries that provide widely used implementations of multigrid methods. Most of the available finite element libraries do to the best of my knowledge (certainly the one I work on, deal.II, does). There are also widely used implementations of algebraic multigrid methods, most notably from the ML and hypre packages that you can most easily ...


3

Actually there's no restriction that the number of intervals in the fine grid should have to be a power of 2. This number can be either odd or even. No big deal. Take 1-dimension problem as an example, suppose the number of intervals in the fine grid is 7(Thus the number of nodes is 8), then you can use the following restriction operator Columns 1 through 7 ...


3

It is true that if the smoother is not convergent for high frequencies, that multigrid will not converge. Note that this is weaker than being globally convergent: a relaxation method that is unstable for low frequencies can be used with "compatible relaxation" to yield a stable smoother that is only convergent when combined with a coarse-grid correction. ...


3

In algebraic multigrid there are usually two steps: 1) Setup: Here we compute the $A$ matrices and interpolation matrices ($W$, $W'$) at each grid level. This is based on computing c-points and f-points which are entirely derived from information from the matrix $A$ at each level. We do not need the right hand side $b$ vector to complete this setup step. 2)...


Only top voted, non community-wiki answers of a minimum length are eligible