Thanks for getting back to me. The matrix solves in each x-y plane are linear. The matrices depend on the z wavenumber and so are different at each x-y slice. The equations are basically Helmholtz and Poisson type. They are 3D, but when done in Fourier space, they decouple so each x-y plane can be solved independently.
I'd like to run on a few hundred processors, but if possible I'd like it to scale to more processors for higher Re. I agree that keeping the z-dimension data local is beneficial for FFTs. Thanks, Brandt On Tue, Nov 15, 2011 at 10:57 AM, Jed Brown <jedbrown at mcs.anl.gov> wrote: > 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/fa5ea81b/attachment.htm>
