Could AMReX use external algebraic multigrid solvers?

[zongy17]

    We usually don't construct an `A` matrix explicitly.  `MLLinOp` is an abstract class, the implementations are in its subclasses.

Originally posted by @WeiqunZhang in https://github.com/AMReX-Codes/amrex/issues/1747#issuecomment-790913749

① First, I currently have a linear system (matrix A, right-hand-side b, and sometimes initial guess x0) derived from an application. It is based on structured mesh and unknowns are located at cell center, so I want to use the geometric multigrid solver in AMReX to solve it effeciently.

But I could not find a way to modify or fill the values of matrix A by myself.

Do you have any suggestions about what I want?

② Second, for a given problem (such as those written in amrex-tutorials/ExampleCodes/LinearSolvers/ABecLaplacian_C or NodalPoisson or MultiComponent), does AMReX provide any ways to export the matrix A, right-hand-side b, and initial guess x0 in assembled forms (such as CSR representation or COO like matrix-market format)?

If export is not supported, I guess there might be some places to see or print the values of matrix A, vector b and x0 in distributed form ?

[WeiqunZhang]

The linear solvers in AMReX solve PDEs whose discretized forms are linear systems. So the starting point for us is a PDE not a matrix A. We do not support exporting the matrix in any formats. As mentioned already we don't have the matrix explicitly. What we have is something like

template <typename T>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void mlpoisson_adotx (int i, int j, Array4<T> const& y,
                      Array4<T const> const& x,
                      T dhx, T dhy) noexcept
{
    y(i,j,0) = dhx * (x(i-1,j,0) - T(2.)*x(i,j,0) + x(i+1,j,0))
        +      dhy * (x(i,j-1,0) - T(2.)*x(i,j,0) + x(i,j+1,0));
}

For the example above, we don't need a matrix A to do A dot x. As for vectors, we store them as MultiFabs, which is AMReX's distributed data container. And we access the vector using multi-dimensional indices, not a single row index.

Does your matrix come from a PDE? If so, what is the PDE? If it's similar to what we have, you might be able to use the solvers we have directly or modify them for your system.

[zongy17]

Thanks for your comments! Btw, now that AMReX does not maintain a matrix explicitly, the setup of geometric multigrid also only relies on its operator form. Is my understanding right? If so, what coarsening types are adopted in AMReX's GMG? It seems that Galerkin-type approach that need to do triple-matrix-product is hard to implement with matrix-free form.

[WeiqunZhang]

As you have mentioned, in the algebraic approach the coarse level matrix can be generated by A_c = R*A_f*P. So it's like PDE -> A_f -> A_c. In GMG, one can generate A_c directly from discretizing the PDE at a coarser resolution.

[zongy17]

I see. Thanks a lot for your comments!

[zongy17]

Hi, @WeiqunZhang , thanks a lot for your previous comments!

Recently, I am developing an algebraic multigrid solver for semi-structured-grid problems (including Block-Structured Mesh, Adative Mesh Refinement, ..., and more types of grids could refer to https://arxiv.org/pdf/2205.14273.pdf or https://arxiv.org/pdf/2103.11962.pdf ), and looking for suitable problems for tests. AMReX is a well-acknowledged framework, so I wonder if I could find some helps here.

Our solver could work as an alternative to HYPRE's SMG, PFMG, SysPFMG, and BoomerAMG. All of them only rely on the matrix on the finest-level grid. I guess this is different from the geometric multigrid embedded in AMReX that generates A_c directly from discretizing the PDE at a coarser resolution (Is my understanding correct?).

But I also notice that in AMReX's documentation (https://amrex-codes.github.io/amrex/docs_html/LinearSolvers_Chapter.html?highlight=multigrid), it said

The default solution technique is geometric multigrid, but AMReX includes native BiCGStab solvers for a single level as well as interfaces to the hypre library.

I am a bit confused about the meaning of "single level" in this statement. Does it mean that AMReX itself constructs a hierarchy of grids, and could use native BiCGStab or HYPRE's solvers (such as the *MG mentioned above) to solve the systems of each level independently? But it seems that only the finest level (but not other levels) needs to be solved?

[WeiqunZhang]

I guess this is different from the geometric multigrid embedded in AMReX that generates A_c directly from discretizing the PDE at a coarser resolution (Is my understanding correct?).

Yes

I am a bit confused about the meaning of "single level" in this statement.

It could use bicgstab or hypre at the bottom of a multigrid v-cycle. If the v-cycle has only one level, it effectively uses bicgstab or hypre for the entire problem without using our multigrid.

[zongy17]

Thanks a lot for your comments! @WeiqunZhang

I wonder how AMReX generates the hierarchy of grids. In the website (https://amrex-codes.github.io/amrex/overview.html), it said

Note that there is no parent-child relationship between coarser grids and finer grids as the grid structure dynamically evolves over time.

So does it start from a coarsest-level grid, and then refine it again and again in the areas of interest (i.e., the generation direction is from coarse to fine)? Or does it coarsen the areas that are not of interest in the given finest-level grid (i.e., the direction is from fine to coarse)?

[WeiqunZhang]

Because this thread is about multigrid, I want to make it clear that the sentence you quoted is not about multigrid. If your question is about AMR grids, then the answer is we have the coarsest level grid first and then fine levels are added.

[zongy17]

Yes, I was talking about AMR grids exactly. Thanks so much for your detailed explanation. It really helps a lot! Thanks!

[zongy17]

Hi, @WeiqunZhang , thanks for your comments. And I have a question based on your previous two comments.

1. AMReX could use bicgstab or hypre at the bottom of a multigrid v-cycle. If the v-cycle has only one level, it effectively uses bicgstab or hypre for the entire problem without using our multigrid.

2. For the AMR grids, we have the coarsest level grid first and then fine levels are added.

Combining these two comments, it seems that they imply that hypre's multigrids could only be used on a level without any refinement. Is my understanding correct? When the v-cycle has only one level (typically by setting max_coarsening_level = 0, right?), that level is the coarsest one without refinement, according to your 2nd comment. Therefore, there is no refinement in the problem that hypre is actually solving, right?

To make my question more clearly, I would like to cite hypre's SSAMG (https://arxiv.org/pdf/2205.14273.pdf) to illustrate here. There are two parts in the above figure, in which the blue part is a local refinement of the red part. The unknowns are cell-centered. Could AMReX call hypre to solve the refined system as a whole (i.e., blue and red parts are solved simultaneously)?

[WeiqunZhang]

I fail to understand the connection between (1) and (2). The linear solver does not care about how AMR grids are generated. Suppose there are two AMR levels. (Again I fail to see the order of how the two levels are created matters in any way whatsoever.) There are three cases for the linear solver. Case 1: single level solve for AMR level 0. Case 2: single level solve for AMR level 1. Case 3: two-level composite solve for both AMR levels together. I assume there is no question for case 1. For case 2, the AMR level 0 is not part of the unknowns, but it provides Dirichlet BC at the AMR coarse/fine boundary. Hypre can be used at the bottom of the v-cycle, where it may have a total of 1, 2, or 3, or ... multigrid levels. (This again has nothing to do with AMR levels.) For case 3, the v-cycle has at least two multigrid levels, because we coarsen the fine AMR level once to have the same cell size as the coarse AMR level. Hypre can be used at the bottom of the v-cyle, which has at least two levels. (This again has nothing to do with the order of how AMR levels are created.) So for your last question, the answer is no. AMReX cannot use hypre to do the whole composite solve. Nevertheless, you might be interested in how we use hypre to do the whole composite solve in Castro. For example, https://github.com/AMReX-Astro/Castro/blob/main/Source/radiation/HypreMultiABec.cpp

[zongy17]

Thanks for your comments! It is actually a bit complicated, and easy to confuse the levels of algebraic multigrid (linear) solver and AMR itself.

I suppose that the level 0 is the coarser one, and 1 is the finer one in the three cases you mentioned. I completely understand your Case 1.

But for Case 2, hypre can be used at the bottom of the v-cycle. What does 'the bottom of v-cycle' refer to here? If it refers to the AMR hierarchy, its bottom seems to be level 0?

For Case 3, does a composite solve mean that all the unknowns in level 0 and 1 (e.g., in the blue and red parts in my cited figure) are coupled in a linear system and get solved simultaneously?

[WeiqunZhang]

For case 2, the AMR level 0 is not part of the unknowns, but it provides Dirichlet BC at the AMR coarse/fine boundary. Hypre can be used at the bottom of the v-cycle, where it may have a total of 1, 2, or 3, or ... multigrid levels. (This again has nothing to do with AMR levels.)

This was what I already said. No, it has nothing to do with AMR levels. For case 2, the linear solver is given a problem to solve. The problem involves a number of rectangular boxes of uniform cells. Whether the grids come from AMR level 0, or 1 or 3 does not matter. It tries to solve the problem using the multigrid method. So it just does the normal multigrid stuff. It does not care about AMR and it has nothing to do with AMR. We usually try to coarsen as much as possible. So it might be coarsened 8 times. Again it has nothing to do with AMR. We call the coarsest MG level the bottom level. A special case is we do not do any coarsening at all and use hypre directly. We still call it the bottom level. Again it has nothing to do with AMR. The bottom level grids happen to come from AMR level 1 in the example. But it does not matter and it does not care.

For case 3, yes the unknown includes both blue and red cells. Yes they are coupled at the coarse/fine AMR boundary. Yes they are solved simultaneously. And no, hypre is not used for fine AMR level's data in this example (this example! not case 2! not a general statement about whether hypre is used for fine AMR levels). Yes, hypre can be used at the bottom level, which might be the coarse AMR level coarsened 0 times or 1 time or 2 times, or 3 times or more.

[WeiqunZhang]

If you have more questions. Let's focus on one question at a time.

[zongy17]

Thanks! It is very kind of you to make me learn a lot. I am a developer of linear solver (especially multigrid), and have minimal knowledge about AMR. To make discussion more clearly, I try drawing a picture in the following.

Is this precisely what happens in Case 2? Linear solvers (e.g., BoomerAMG or PFMG in hypre, or others) take a matrix of a specific AMR level (maybe 0, 1, ..., or 3) as the input, and return the solution on this AMR level. Then AMReX will continue its computation of the AMR hierarchy. Is my understanding correct?

[WeiqunZhang]

Yes

[WeiqunZhang]

Your picture is like we set max_coarsening_level to 0 and call hypre. The user could also set max_coarsening_level to 1 and use hypre at the bottom. In that case, amrex does a 2-level v-cycle itself. What you labeled as "direct solver" would be a call to hypre.

[WeiqunZhang]

We call the bottom tip of the v-cycle the bottom level. If the v is just a dot, I still call it the bottom level in our discussions. In our current setup, hypre is only used at the bottom level, whether the v-cycle is just a dot or includes a single AMR level (could be level 0 or 1 or 2) or includes multiple AMR levels.

[WeiqunZhang]

In your three level case, the linear solver might need to solve (1) level 0 only, (2) level 1 only, (3) level 2 only, (4) levels 0 & 1, (5) levels 1 & 2, (6) levels 0, 1, & 2.

[zongy17]

Is this what you mean?

In AMReX, we could make solve on different single level (e.g. in my picture, level 2 is solved by purple or yellow v-cycle, and level 1 is solved by green v-cycle), and set max_coarsening_level.

The colored v-cycles are still in AMReX. And only the bottom tip will call external solvers like hypre's BoomerAMG or PFMG. Is my understanding right?

[WeiqunZhang]

It seems that we are back to square one. For a single level solve on AMR level 2, you should think level 0 and level 1 do NOT exist. Suppose there are n2 cells on level 2 and n1 cells on level 1. If we do a single level solve on level 2 and set max_coarsening_level to 1. The coarse multigrid level has n2/8 cells in 3d. And it is usually the case n2/8 < n1. Again, For a single level solve on AMR level 2, you should think level 0 and level 1 do NOT exist, other than the level 1 data provide boundary conditions for level 2. This was reply to a deleted comment.

[WeiqunZhang]

For the current comment above, the answer is yes.

[zongy17]

Hi @WeiqunZhang , I draw another picture here. Is this what you mean?

[WeiqunZhang]

Yes. I also want to note that the depth of the v-cyle in amrex can range from 0 to as much as possible by geometric coarsening.

[WeiqunZhang]

And the picture you are drawing are for single level solves only.

[zongy17]

Thanks a lot for helping me to understand.

the depth of the v-cyle in amrex can range from 0 to as much as possible by geometric coarsening

Does "as much as possible" mean that the box width must be a number of power of two?

the picture you are drawing are for single level solves only.

Yes. I think I have understand what happens in Case 2 thanks to your help. So let's move on the next question about Case 3. How about the composite solve? I suppose that the cell-centered knowns on level 0, 1, 2 are marked as "x" in the picture here.

They are coupled at the coarse/fine AMR boundary, and solved simultaneously. You mentioned that AMReX could not use hypre in this case. What technique is used in AMReX to solve them together?

[WeiqunZhang]

The number of cells does not need to be power of 2, even though it's preferred. If we have 24 cells, we can coarsen to 12 cells, then 6 cells, and then 3 cells. At that point, we cannot coarsen it further.

For a composite solve with levels 0, 1 and 2, level 2 is the top of the v-cycle, level 1 is the level below in the v-cycle, level 0 is even further below, and AMR level 0 may have more coarse multigrid levels. In the hierarchy of the v-cycle grids, the top level covers only a fraction of the physical space of the AMR level 1. We can do smoothing on the top level, then we can restrict the residual to the region below AMR level 2. That region only covers part of the AMR level 1 grids. We compute the residual in the coarse cells next to the coarse/fine boundary using both the coarse data and the latest fine AMR level data. For AMR level 1 cells away from the coarse/fine boundary, we compute residual as usual. We now have the residual on AMR level 1. We can repeat the process like this. In the end, when the iteration converges, the stencil at the coarse/fine boundary is a composite stencil that uses both coarse and fine AMR data. In the compost solve case, we can use hypre at the multigrid bottom level just like the single level case.

[zongy17]

Thanks for your detailed comment! I am curious that which kind of method does AMReX usually prefer? It seems that the composite solve is more sophisticated, and thus maybe more efficient?

Actually, the case of composite solve is what we are seeking for. Thanks for your reminding that AMReX could not use hypre in such cases, and recommending Castro (https://github.com/AMReX-Astro/Castro/blob/main/Source/radiation/HypreMultiABec.cpp).

But unfortunately, when I ran the radiation cases of Castro, I found Castro not working as my expectations. It seemed that it failed to do a true composite solve. I raised the question with more details in https://github.com/AMReX-Astro/Castro/issues/2730. Could you help me to figure it out? @WeiqunZhang Thank you so much!