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

[Dave May]

For your pressue Poisson solve, you can use DMDAGetReducedDMDA to generate
a DMDA with only 1 dof. This reduced DMDA will have the same parallel
layout as your DMDA with 4 dofs.

See
http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMDAGetReducedDMDA.html

Cheers,
  Dave

On 15 October 2014 15:41, Håkon Strandenes <haakon at hakostra.net> wrote:

> Thanks for the reply, both Jed and Matthew. However, I am still a bit
> unsure on a small detail:
>
> I have created a 3D DMDA with ndof = 4, one for each velocity component
> (u,v,w) and one for the pressure (p). I only solve a Poisson equation to
> correct p, the velocities are solved explicitly.
>
> My question is this time:
> When I create a matrix with DMCreateMatrix(), I get four times as many
> rows as I really need, i.e. it assumes that I want to solve equations for
> all my DOF's. The same happens when I create a vector with
> DMCreateGlobalVector(), it gets one element for each dof. This is not
> desirable, since I only solve equations for one unknown per gridpoint. If
> my domain (for debugginf purposes) have 250 points, I want the local
> portion of my matrix to have 250 rows. Instead it gets 1000.
>
> Was I wrong when I created my DMDA with ndof = 4? Should I have created
> one DMDA for the velocity with ndof = 3 and one with ndof = 1 to use as a
> basis for the pressure solution? Or have I again missed something
> fundamentally?
>
> Regards,
> Håkon
>
>
>
> On 08. okt. 2014 22:07, Jed Brown wrote:
>
>> 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 --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20141015/e7649e50/attachment.html>