12

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 ...


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

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 ...


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 ...


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

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

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


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

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 ...


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

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

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

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

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

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 ...


4

Most people use (algebraic as well as geometric) multigrid as preconditioners these days. It's an empirical observation that that leads to faster convergence in terms of iterations, given that in a typical CG iteration, applying a multigrid preconditioner takes up the vast majority of the effort compared to all other operations. In other words, you get the ...


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

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)...


3

This question is pretty well discussed in literature. However, there are lots of questions concerning multigrid on SciComp, so I decided to compose more or less detailed answer. I. When multigrid DOES work fine as a stand–alone solver MG works fine for rather simple elliptic problems. In fact, both MG and PCG (with MG cycles as inner iterations) usually ...


3

Check Algebraic Multigrid Solvers for Complex-Valued Matrices by MacLachlan and Osterlee. On the implementation side of things, PyAMG supports complex-valued matrices. I've used it before and it works well.


3

The point important to understand when thinking about multigrid is that the lower levels of the hierarchy do not actually have to solve the problem accurately. Rather, the operators at the lower levels just need to provide good approximations of a part of the spectrum (eigenvalues) of the operator on the finest level -- specifically, they need to well ...


3

Every iterative solver -- Jacobi, SSOR, CG, etc -- starts with an initial approximation. One often just uses the zero vector, but there is nothing wrong with using the solution of the previous time step. In fact, extrapolating from previous time steps to the current one is an even better idea -- one the authors apparently missed! For some iterative solvers, ...


3

If I understand your question correctly you're solving a linear elasticity problem using conjugate gradient and it's preconditioned with a preconditioned AMG solver? It seems to me that this may be overkill for a pretty well behaved problem, and that could be why you don't see much of a speed-up. Just to elaborate a bit. I think it makes more sense to just ...


2

The regular Gauss-Seidel method is not a good smoother for higher order discretizations. From my experience with finite element multigrid, I expect that a block Gauss-Seidel with block size comparable to the width of the stencil should do the job. Example: if your stencil involves neighbors and second neighbors, use blocks of size 4 or 9. These blocks ...


2

Maybe you could make one array that's sized for the total number of points, then fill it with coarse points from the front and fine points from the back. They'll meet up somewhere in the middle (but not overlap).


2

Jacobi is usually not a great smoother, Gauss-Seidel is a better choice. You might also want to use a larger initial grid so you see whether the defect is being generated at the boundary or in the interior. Your 2x2 coarse grid is essentially all boundary condition. Mixed boundaries are tricky to implement. I'd get only Dirichlet boundaries working, and ...


2

There is currently no documented 'FEniCS way' to do this. However, since FEniCS is a pretty standard finite element code behind all the UFL and code generation magic, you can implement things like transfer operators by yourself. The only difficulty is that they have no built-in mechanism to deal with the inter-level mappings in hierarchically refined meshes ...


2

What you are thinking of is something that uses the structure of the augmented matrix to make solution of the system simpler. For example, one could be tempted to think of forming the Schur complement with regard to the bottom right $2\times 2$ block of the rewritten system. But I don't think anything like this exists. If it would, then you would have ...


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