[petsc-users] preallocation after DMCreateMatrix?

[Matthew Knepley]

Thanks for finding this bug. I guess no one has been making matrices with
FVM. I will fix this internally, but here
is a workaround which makes the code go for now.

  Thanks,

     Matt

On Wed, Nov 29, 2017 at 6:36 AM, Matteo Semplice <matteo.semplice at unito.it>
wrote:

>
>
> 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
> > 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> 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
>
>
>


-- 
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

https://www.cse.buffalo.edu/~knepley/ <http://www.caam.rice.edu/~mk51/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20171130/7425b0de/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: testfvm.c
Type: text/x-csrc
Size: 2063 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20171130/7425b0de/attachment-0001.bin>