[petsc-users] DMPlex example with manual grid layout

[Matthew Knepley]

On Wed, Oct 30, 2013 at 11:13 AM, Bernhard Reinhardt <
b.reinhardt at physik.uni-muenchen.de> wrote:

> On 28.10.2013 16:43, Matthew Knepley wrote::
>
>> On Mon, Oct 28, 2013 at 6:17 AM, Bernhard Reinhardt
>> <b.reinhardt at physik.uni-**muenchen.de<b.reinhardt at physik.uni-muenchen.de>
>> <mailto:b.reinhardt at physik.**uni-muenchen.de<b.reinhardt at physik.uni-muenchen.de>>>
>> wrote:
>>
>>     Hi there,
>>
>>     I am trying to familiarize myself with unstructured grids. For this
>>     I started with the example given in the petsc-manual (mesh of 2
>>     triangles) and tried to modify it according to my needs.
>>
>>     However, the example seems to be incomplete. At least I was not able
>>     to create a matrix without applying a DMPlexSetDimension.
>>
>>     In the end I decided to set up an an example as simple as possible:
>>     1D transport of a scalar in a bar with periodic
>>
>
> aries.
>
>>
>>     Lets say I have a bar consisting of 3 volumes and 3 faces.
>>     You could write this down as a vector of length 3 containing the
>>     scalar in the volumes and a 3x3 matrix containing the transport
>>     coefficients on the faces.
>>
>>     I have already tried many different ways of setting this up with
>>     plex but always failed. Attached you find the version that makes
>>     most sense to me. It gives me vectors of length 3 and a 3x3 matrix.
>>     However, the matrix contains no writable entries and I think the
>>     returned closures are not correct.
>>
>>     Is there a simple example with good documentation of the individual
>>     steps like in the manual which is complete and working?
>>
>>
>> I have modified your example and it is attached. Here are a few points:
>>
>> 1) You ask for the closure of point 4 in your Vec, and get 1.0. This is
>> correct.
>>
>>    The closure of 4 is {4, 1, 2}, the cell + its vertices. There are
>> only values on cells,
>>    so you get the value on 4, which is 1.0.
>>
>> 2) The default for preallocation is FEM-style topology, namely that
>> unknowns are connected to other
>>       unknowns with overlapping support. The support of your unknowns is
>> only the cell itself, so there
>>       are no connections by default. I put in a call to
>> DMPlexSetPreallocationCenterDi**mension(), which
>>       is the dimension of intermediate points. In FEM these are cells,
>> but in FVM they are vertices (or
>>       faces). I set your topology to FVM, and you should get the
>> preallocation you want.
>>
>
> Dear Matt,
>
> thank you for your advice. Now I get the right allocation for my exemplary
> problem (although I have to admit I´m not entirely convinced for the right
> reasons ;-). I expanded the bar from 3 to 5 volumes, so that the allocation
> becomes clearer. See attached file.
>
> However, I still don't get what DMPlexMatSetClosure does.


The words used are precise. Closure means the transitive closure of a point
in the graph. For example, the
closure of point 7 is

  closure(7) = {7, 2, 3}

MatSet() takes a set of indices and a set of values and updates a Mat. Now
MatSetClosure(), first calculates the indices

  for q in closure(p):
    PetscSectionGetOffset(s, q, &off)
    PetscSectionGetDof(s, q, &dof)
    indices += [ind, ind+dof)
  MatSetValues(indices, values)


> What I would like to do, is set the operator into the matrix which
> computes the change of the value at point p. So let's say in my example the
> Laplacian operator (1, -2, 1).
> E.g. Delta_p7 = A_21 p6 + A_22 p7 + A_23 p8; with A_21=1, A_22 = -2, A_23
> = 1. However, in my example MatClosure(p7) consits only of A_22.
>

What you seem to be suggesting is FD, which you can do on a structured grid
using DMDA. On an unstructured grid, what I am familiar
with is FV, which means looping over faces instead of cells, and doing
something like MatSetStar() instead of MatSetClosure().

If you are using some unstructured method which sets whole rows of the
Jacobian at once (like FD), then you would want something like
MatSetStarClosure(), and then you would need to add the flag for ghost
cells in DMPlexDistribute(). I have never done this since all
methods I have implemented are more local.

It would help to understand what you really want to do.

   Matt


> Probably my sieve-layout is not appropriate. I recognize that the
> Laplacian consists of 3 elements, but my volumes have only a cone size of
> 2. Nevertheless I hoped to be able to set at least A_21 and A_23.
>
> I tried all combinations of dof=1 or dof=0 at the faces or volumes and
> PreallocationCenterDimension 0 or 1 without elucidation. I also tried
> different mesh layouts without much succes. That is
>
> - - - - - (without any faces)
> or
>
> |-||-||-||-|-| (with double faces)
> or
>  _  _  _  _ _
> |-||-||-||-|-| (with double faces and an additional "internal" face)
>
> In the best case I get the foreseen allocation, but am stil am not able to
> set the right values in the matrix.
>
> This leads me to the question: Is the sieve-layout appropriate
> (|-|-|-|-|-)? If yes, what is the envisaged way to set the operator for
> point p in the matrix?
>
> Best regards!
>
> Bernhard
>



-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20131030/a38cb788/attachment-0001.html>