[petsc-users] preallocation after DMCreateMatrix?
Matteo Semplice
matteo.semplice at unito.it
Fri Dec 1 04:50:50 CST 2017
Thanks for the fix!
(If you need a volunteer for testing the bug-fix, drop me a line)
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 <mailto: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 <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
>
>
>
>
>
> --
> 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/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20171201/0c1bb0b4/attachment.html>
More information about the petsc-users
mailing list