[petsc-users] preallocation after DMCreateMatrix?

[Matteo Semplice]

On 29/11/2017 12:46, Matthew Knepley wrote:
> On Wed, Nov 29, 2017 at 2:26 AM, Matteo Semplice 
> <matteo.semplice at unito.it <mailto:matteo.semplice at unito.it>> wrote:
>
>     On 25/11/2017 02:05, Matthew Knepley wrote:
>>     On Fri, Nov 24, 2017 at 4:21 PM, Matteo Semplice
>>     <matteo.semplice at unito.it <mailto:matteo.semplice at unito.it>> wrote:
>>
>>         Hi.
>>
>>         The manual for DMCreateMatrix says "Notes: This properly
>>         preallocates the number of nonzeros in the sparse matrix so
>>         you do not need to do it yourself", so I got the impression
>>         that one does not need to call the preallocation routine for
>>         the matrix and indeed in most examples listed in the manual
>>         page for DMCreateMatrix this is not done or (KSP tutorial
>>         ex4) it is called declaring 0 entries per row.
>>
>>         However, if read in a mesh in a DMPlex ore create a DMDA and
>>         then call DMCreateMatrix, the resulting matrix errors out
>>         when I call MatSetValues. I have currently followed the
>>         suggestion of the error message and call MatSetOption(A,
>>         MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE), but I'd like to
>>         fix this properly.
>>
>>
>>     It sounds like your nonzero pattern does not obey the topology.
>>     What nonzero pattern are you trying to input?
>
>     Hi.
>
>         The problem with the DMDA was a bug in our code. Sorry for the
>     noise.
>
>     On the other hand, with the DMPLex, I still experience problems.
>     It's a FV code and I reduced it to the case of the simple case of
>     a laplacian operator: I need non-diagonal entries at (i,j) if cell
>     i and cell j have a common face.
>     Here below is my code that reads in a grid of 4 cells (unit square
>     divided by the diagonals), create a section with 1 dof per cell,
>     creates a matrix and assembles it.
>
>
>       ierr = DMPlexCreateFromFile(PETSC_COMM_WORLD, "square1.msh",
>     PETSC_TRUE, &dm);CHKERRQ(ierr);
>
>
> Can you show me the mesh?
>
>          ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr);
>
> and then run with -dm_view ::ascii_info_detail
>
>
>       ierr = DMPlexSetAdjacencyUseCone(dm, PETSC_TRUE); CHKERRQ(ierr);
>       ierr = DMPlexSetAdjacencyUseClosure(dm, PETSC_FALSE); CHKERRQ(ierr);
>
>
> This looks like the right FVM adjacency, but your matrix is diagonal 
> it appears below. TS ex11
> has an identical call, but produces the correct matrix, which is why I 
> want to look at your mesh.

I have put the DMViewFromOptions call after the SetAdjacency calls and 
this is the output:

matteo at signalkuppe:~/software/petscMplex$ ./mplexMatrix -dm_view 
::ascii_info_detail
DM Object: 1 MPI processes
   type: plex
Mesh 'DM_0x557496ad8380_0':
orientation is missing
cap --> base:
[0] Max sizes cone: 3 support: 4
[0]: 4 ----> 9
[0]: 4 ----> 11
[0]: 4 ----> 12
[0]: 5 ----> 10
[0]: 5 ----> 11
[0]: 5 ----> 15
[0]: 6 ----> 14
[0]: 6 ----> 15
[0]: 6 ----> 16
[0]: 7 ----> 12
[0]: 7 ----> 13
[0]: 7 ----> 16
[0]: 8 ----> 9
[0]: 8 ----> 10
[0]: 8 ----> 13
[0]: 8 ----> 14
[0]: 9 ----> 0
[0]: 9 ----> 1
[0]: 10 ----> 0
[0]: 10 ----> 2
[0]: 11 ----> 0
[0]: 12 ----> 1
[0]: 13 ----> 1
[0]: 13 ----> 3
[0]: 14 ----> 2
[0]: 14 ----> 3
[0]: 15 ----> 2
[0]: 16 ----> 3
base <-- cap:
[0]: 0 <---- 9 (0)
[0]: 0 <---- 10 (0)
[0]: 0 <---- 11 (0)
[0]: 1 <---- 12 (0)
[0]: 1 <---- 13 (0)
[0]: 1 <---- 9 (-2)
[0]: 2 <---- 10 (-2)
[0]: 2 <---- 14 (0)
[0]: 2 <---- 15 (0)
[0]: 3 <---- 14 (-2)
[0]: 3 <---- 13 (-2)
[0]: 3 <---- 16 (0)
[0]: 9 <---- 4 (0)
[0]: 9 <---- 8 (0)
[0]: 10 <---- 8 (0)
[0]: 10 <---- 5 (0)
[0]: 11 <---- 5 (0)
[0]: 11 <---- 4 (0)
[0]: 12 <---- 4 (0)
[0]: 12 <---- 7 (0)
[0]: 13 <---- 7 (0)
[0]: 13 <---- 8 (0)
[0]: 14 <---- 8 (0)
[0]: 14 <---- 6 (0)
[0]: 15 <---- 6 (0)
[0]: 15 <---- 5 (0)
[0]: 16 <---- 7 (0)
[0]: 16 <---- 6 (0)
coordinates with 1 fields
   field 0 with 2 components
Process 0:
   (   4) dim  2 offset   0 0. 0.
   (   5) dim  2 offset   2 0. 1.
   (   6) dim  2 offset   4 1. 1.
   (   7) dim  2 offset   6 1. 0.
   (   8) dim  2 offset   8 0.5 0.5

For the records, the mesh is loaded from the (gmsh generated) file

==== square1.msh =====
$MeshFormat
2.2 0 8
$EndMeshFormat
$Nodes
5
1 0 0 0
2 0 1 0
3 1 1 0
4 1 0 0
5 0.5 0.5 0
$EndNodes
$Elements
12
1 15 2 0 1 1
2 15 2 0 2 2
3 15 2 0 3 3
4 15 2 0 4 4
5 1 2 0 1 4 3
6 1 2 0 2 3 2
7 1 2 0 3 2 1
8 1 2 0 4 1 4
9 2 2 0 6 1 5 2
10 2 2 0 6 1 4 5
11 2 2 0 6 2 5 3
12 2 2 0 6 3 5 4
$EndElements
===================

Thanks,
     Matteo


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20171129/c60c54c5/attachment.html>