[petsc-users] Using petsc for banded matrices and 2D finite differences

[Jed Brown]

On Tue, Nov 15, 2011 at 09:43, Brandt Belson <bbelson at princeton.edu> wrote:

> I'm writing a 3D incompressible fluids solver for transitional and
> turbulent boundary layers, and would like to make use of petsc if possible.
> At each time-step I'll need to solve matrix equations arising from finite
> differences in two dimensions (x and y) on a structured grid. The matrix is
> block tri/penta-diagonal, depending on the stencil, and the blocks are also
> tri/penta-diagonal. Correct me if I'm wrong, but I believe these types of
> matrix equations can be solved directly and cheaply on one node.
>
> The two options I'm comparing are:
> 1. Distribute the data in z and solve x-y plane matrices with LAPACK or
> other serial or shared-memory libraries. Then do an MPI all-to-all to
> distribute the data in x and/or y, and do all computations in z. This
> method allows all calculations to be done with all necessary data available
> to one node, so serial or shared-memory algorithms can be used. The
> disadvantages are that the MPI all-to-all can be expensive and the number
> of nodes is limited by the number of points in the z direction.
>
> 2. Distribute the data in x-y only and use petsc to do matrix solves in
> the x-y plane across nodes. The data would always be contiguous in z. The
> possible disadvantage is that the x-y plane matrix solves could be slower.
> However, there is no need for an all-to-all and the number of nodes is
> roughly limited by nx*ny instead of nz.
>
> The size of the grid will be about 1000 x 100 x 500 in x, y, and z, so
> matrices would be about 100,000 x 100,000, but the grid size could vary.
>
> For anyone interested, the derivatives in x and y are done with compact
> finite differences, and in z with discrete Fourier transforms. I also hope
> to make use petsc4py and python.
>

Are the problems in the x-y planes linear or nonlinear? Are the
coefficients the same in each level? What equations are being solved in
this direction? (You can do much better than "banded" solvers, ala LAPACK,
for these plane problems, but the best methods will depend on the problem.
Fortunately, you can don't have to change the code to change the method.)

About how many processes would you like to run this problem size on? Since
the Fourier direction is densely coupled, it would be convenient to keep it
local, but it's probably worth partitioning if your desired subdomain sizes
get very small.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20111115/afd357fd/attachment.htm>