[petsc-users] preallocation after DMCreateMatrix?

Matthew Knepley knepley at gmail.com
Mon Dec 4 10:01:22 CST 2017 

On Fri, Dec 1, 2017 at 5:50 AM, Matteo Semplice <matteo.semplice at unito.it>
wrote:

> Thanks for the fix!
>
> (If you need a volunteer for testing the bug-fix, drop me a line)
>
Cool. Its in next, and in

  https://bitbucket.org/petsc/petsc/branch/knepley/fix-plex-fvm-adjacency

  Thanks,

    Matt


> Best,
>     Matteo
>
>
> On 30/11/2017 13:46, Matthew Knepley wrote:
>
> 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/%7Emk51/>
>
>
>


-- 
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/20171204/8d875f35/attachment.html>

More information about the petsc-users mailing list