[petsc-users] 3D multi-dof DMDA and Poisson equation

[Jed Brown]

Håkon Strandenes <haakon at hakostra.net> writes:

> Hi,
>
> I am making a new, simple incompressible CFD code. I have created a grid 
> as a 3D DMDA, with ndof=4, one dof for each velocity component (u,v,w) 
> and one for the pressure (p). Currently I am creating a forward Euler 
> step, where the velocities are updated based on the velocities and 
> pressure from the previous time step. That works great. Now I need to 
> solve a Poisson equation for the pressure correction term, enforcing 
> continuity onto my solution, and this is where I have a few questions 
> for you:
>
> 1)
> How do I set up my equation system in the most efficient way? The 
> vectors from the DMDA is not contiguous in a global sense, i.e. if I map 
> the local (i,j,k) indices to a global index with
>
> globId = (info.mx*info.my)*k + (info-.mx)*j + i;

This is a "natural index", used in some cases for IO, but not for
computation.

> globId is not gontignous on any process. How do I create my coefficient 
> matrix, r.h.s. vector and solution vector? 

DMCreateGlobalVector()

it will use what is described in the Users Manual as "PETSc ordering".
Please read this section of the manual; I think it will answer your
questions.

> Is there some way I can make the local parts of the matrix and vectors
> mimic the decomposition from the DMDA? Or should I push all values
> into a different decomposition scheme, where each process holds one
> contigous slab of the equation system using the global ID's calculated
> above? Have I missed some really basic stuff here?
>
> 2)
> I know that my coefficient matrix from the Poisson equation is 
> symmetric. Currently I am still creating and inserting values on both 
> the upper and lower diagonal. Is this really necessary? Or is there a 
> "smart" way I can save half the storage by only creating the upper or 
> lower part?

-dm_mat_type sbaij

> 3)
> I solve a transient problem, but my coefficient matrix is the same for 
> all time steps and I only recompute the r.h.s. vector for each time 
> step. I have got the impression that KSPsolve "finds out" if the A 
> matrix have changed since last invocation or not, recalculating the 
> preconditioner as necessary. 

If you just call KSPSolve, it will reuse everything.  Call
KSPSetOperators (assuming you are manually doing the time
discretization) if the matrix changes.

> But still I feel that I should utilize this property more. Do you have
> any tips in how I can use this in the most efficient way possible? For
> example really, really efficient, but expensive, preconditioners?
>
> Thanks in advance.
>
> Regards,
> Håkon Strandenes
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 818 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20141008/61aed57d/attachment.pgp>