[petsc-users] Writing a domain decomposition code with PETSc

[Åsmund Ervik]

Dear all,

We have a two-phase incompressible Navier-Stokes solver written in
Fortran where we use PETSc for solving the pressure Poisson equation.
Since both PETSc and parallelism was an afterthought to this code, it
doesn't scale well at all, so I am tasked with re-writing the whole
thing now. Before I commit any fresh mistakes in the design of this new
code, I will ask for input on my "design decisions" so far.

I want to do domain decomposition on a structured 3D grid. I've been
trying to wrap my head around the DM and DMDA parts of PETSc, and as far
as I understand, these will help me solve the pressure Poisson equation
on a decomposed domain (and with geometric multigrid via Galerkin)
fairly easily.

The tricky part, then; it seems that I must handle "the rest" of the
domain decomposition myself. Omitting some detail, this means my code will:

* set up parameters, initial conditions, etc.
* decompose my array for the velocity field into several parts,
* time loop:
	* communicate e.g. the velocity field on the boundaries
	* each mpi worker will calculate on the local domain the
	  intermediate velocity field, the rhs to the Poisson equation
	  and set up the correct sparse matrix
	* PETSc will solve the Poisson equation to give me the pressure
	* each mpi worker will then calculate the updated
	  divergence-free velocity field
	* each mpi worker will calculate the time step (CFL condition),
	  and we choose the lowest dt among all nodes
* end time loop

Have I misunderstood anything here? At first I thought the DMDA would
give me the framework for decomposing the velocity field, handling
communication of the ghost values at the boundaries etc, but it seems
this is not the case?

One further question: is it a good idea to set up the DMDA letting PETSc
decide the number of processors in each direction, and then using this
same partition for the rest of my code?

If there are any unclear details, please ask. If it matters, I am using
the level-set and ghost-fluid methods, so the matrix for my Poisson
equation must be recomputed each time step. I believe this is the same
situation as Michele Rosso who posted on this list recently.

Best regards,
Åsmund Ervik