Re: [petsc-users] Modify matrix nonzero structure

[Adrian Croucher]

Thanks again Barry, it's working fine for me now too.

- Adrian

On 29/05/24 1:27 pm, Barry Smith wrote:


   There was a bug in my fix for parallel which I have fixed. You will 
need to


   git fetch
   git checkout main
   git branch -D barry/2024-05-27/fix-bug-baij-setvaluesblocked/release
   git checkout barry/2024-05-27/fix-bug-baij-setvaluesblocked/release

    I get the same results with your example with 1,2, and 3 ranks


On May 28, 2024, at 7:45 PM, Adrian Croucher 
 wrote:


hi Barry,

Thanks, that change has fixed the error on 2 ranks for me.

When I run on 3 ranks, there is no error, but it doesn't actually add 
the extra values in to the matrix. Do you see that behaviour too?


- Adrian

On 29/05/24 4:33 am, Barry Smith wrote:


   Adrian,

   I could reproduce with 3 MPI ranks.

   Another error I had to fix. I also added a test example

   SInce I rebased the branch you will need to do something like

   git fetch
   git checkout main
   git branch -D barry/2024-05-27/fix-bug-baij-setvaluesblocked/release
  git checkout barry/2024-05-27/fix-bug-baij-setvaluesblocked/release

  Thanks for your patience

   Barry


On May 27, 2024, at 10:42 PM, Adrian Croucher 
 wrote:


Hmm, that's a bit weird. I haven't modified the test code - I 
checked the file date to make sure. I also tried deleting my PETSc 
build dir and rebuilding it, then rebuilding the test code. I still 
get the error if I run on 2 ranks with -dm_mat_type (or -mat_type) 
baij or mpibaij, but it's fine on aij or mpiaij.


My actual (non-test) code is however working ok now, in serial or 
parallel. So I don't think this should hold up merging your bugfix.


- Adrian

On 28/05/24 2:13 pm, Barry Smith wrote:
  When I run the exact code you sent with two ranks and—mat_type 
mpibaij, it runs as expected. If you modified the code in any way 
to demonstrate the bug, please send the modified code.




On May 27, 2024, at 9:37 PM, Adrian Croucher 
 wrote:


hi Barry,

On 28/05/24 7:46 am, Barry Smith wrote:


   Thanks for reporting this. It is a bug. I have a fixed branch 
*barry/2024-05-27/fix-bug-baij-setvaluesblocked/release * and 
associated merge request 
https://urldefense.us/v3/__https://gitlab.com/petsc/petsc/-/merge_requests/7578__;!!G_uCfscf7eWS!dyTRr3EmobjcjKs0kzURm7EjPlyR-L6yBosZe_dky8HADNvWRJjBCsqt6Quq20DTsiHcW255AIkxd7mdGIssOXoFXVK_OG1L$ 

Thanks very much, that does appear to fix the bug when I run my 
test program in serial.


If I run it on 2 processes, it is OK with AIJ matrix type, but 
with BAIJ I get an error (see below). Is there another problem, 
or I am doing something else wrong?


- Adrian

[0]PETSC ERROR: - Error Message 
--

[0]PETSC ERROR: Petsc has generated inconsistent data
[0]PETSC ERROR: Incorrect colmap
[0]PETSC ERROR: See https://urldefense.us/v3/__https://petsc.org/release/faq/__;!!G_uCfscf7eWS!dyTRr3EmobjcjKs0kzURm7EjPlyR-L6yBosZe_dky8HADNvWRJjBCsqt6Quq20DTsiHcW255AIkxd7mdGIssOXoFXe8LPGgS$  for trouble 
shooting.

[0]PETSC ERROR: Petsc Release Version 3.21.1, unknown
[0]PETSC ERROR: ./matmodify on a main-debug named EN438880 by 
acro018 Tue May 28 13:33:46 2024
[0]PETSC ERROR: Configure options --with-x --download-hdf5 
--download-zlib --download-netcdf --download-pnetcdf 
--download-exodusii --download-triangle --download-ptscotch 
--download-chaco --download-hypre
[0]PETSC ERROR: #1 MatSetValuesBlocked_MPIBAIJ() at 
/home/acro018/software/PETSc/code/src/mat/impls/baij/mpi/mpibaij.c:448
[0]PETSC ERROR: #2 MatSetValuesBlocked() at 
/home/acro018/software/PETSc/code/src/mat/interface/matrix.c:2030
[0]PETSC ERROR: - Error Message 
--
[0]PETSC ERROR:   It appears a new error in the code was 
triggered after a previous error, possibly because:
[0]PETSC ERROR:   -  The first error was not properly handled via 
(for example) the use of

[0]PETSC ERROR: PetscCall(TheFunctionThatErrors()); or
[0]PETSC ERROR:   -  The second error was triggered while 
handling the first error.
[0]PETSC ERROR:   Above is the traceback for the previous 
unhandled error, below the traceback for the next error
[0]PETSC ERROR:   ALL ERRORS in the PETSc libraries are fatal, 
you should add the appropriate error checking to the code

[0]PETSC ERROR: Petsc has generated inconsistent data
[0]PETSC ERROR: Incorrect colmap
[0]PETSC ERROR: See https://urldefense.us/v3/__https://petsc.org/release/faq/__;!!G_uCfscf7eWS!dyTRr3EmobjcjKs0kzURm7EjPlyR-L6yBosZe_dky8HADNvWRJjBCsqt6Quq20DTsiHcW255AIkxd7mdGIssOXoFXe8LPGgS$  for trouble 
shooting.

[0]PETSC ERROR: Petsc Release Version 3.21.1, unknown
[0]PETSC ERROR: ./matmodify on a main-debug named EN438880 by 
acro018 Tue May 28 13:33:46 2024
[0]PETSC ERROR: Configure options --with-x --download-hdf5 
--download-zlib --download-netcdf --download-pnetcdf 
--download-exodusii --download-triangle --download-ptscotch 
--download-chaco

Re: [petsc-users] Modify matrix nonzero structure

[Adrian Croucher]

hi Barry,

Thanks, that change has fixed the error on 2 ranks for me.

When I run on 3 ranks, there is no error, but it doesn't actually add 
the extra values in to the matrix. Do you see that behaviour too?


- Adrian

On 29/05/24 4:33 am, Barry Smith wrote:


   Adrian,

   I could reproduce with 3 MPI ranks.

   Another error I had to fix. I also added a test example

   SInce I rebased the branch you will need to do something like

   git fetch
   git checkout main
   git branch -D barry/2024-05-27/fix-bug-baij-setvaluesblocked/release
  git checkout barry/2024-05-27/fix-bug-baij-setvaluesblocked/release

  Thanks for your patience

   Barry


On May 27, 2024, at 10:42 PM, Adrian Croucher 
 wrote:


Hmm, that's a bit weird. I haven't modified the test code - I checked 
the file date to make sure. I also tried deleting my PETSc build dir 
and rebuilding it, then rebuilding the test code. I still get the 
error if I run on 2 ranks with -dm_mat_type (or -mat_type) baij or 
mpibaij, but it's fine on aij or mpiaij.


My actual (non-test) code is however working ok now, in serial or 
parallel. So I don't think this should hold up merging your bugfix.


- Adrian

On 28/05/24 2:13 pm, Barry Smith wrote:
  When I run the exact code you sent with two ranks and—mat_type 
mpibaij, it runs as expected. If you modified the code in any way to 
demonstrate the bug, please send the modified code.




On May 27, 2024, at 9:37 PM, Adrian Croucher 
 wrote:


hi Barry,

On 28/05/24 7:46 am, Barry Smith wrote:


   Thanks for reporting this. It is a bug. I have a fixed branch 
*barry/2024-05-27/fix-bug-baij-setvaluesblocked/release * and 
associated merge request 
https://urldefense.us/v3/__https://gitlab.com/petsc/petsc/-/merge_requests/7578__;!!G_uCfscf7eWS!bEqLWafwqHADoOBGkMAUu-5jnxuJg95froHeBJtxmkLCGFqP1MFGgnLKOgYILALF3hYWj_cgiKOh5rEeun-2sybwJhqMoBbg$ 

Thanks very much, that does appear to fix the bug when I run my 
test program in serial.


If I run it on 2 processes, it is OK with AIJ matrix type, but with 
BAIJ I get an error (see below). Is there another problem, or I am 
doing something else wrong?


- Adrian

[0]PETSC ERROR: - Error Message 
--

[0]PETSC ERROR: Petsc has generated inconsistent data
[0]PETSC ERROR: Incorrect colmap
[0]PETSC ERROR: See https://urldefense.us/v3/__https://petsc.org/release/faq/__;!!G_uCfscf7eWS!bEqLWafwqHADoOBGkMAUu-5jnxuJg95froHeBJtxmkLCGFqP1MFGgnLKOgYILALF3hYWj_cgiKOh5rEeun-2sybwJvu9iBaG$  for trouble 
shooting.

[0]PETSC ERROR: Petsc Release Version 3.21.1, unknown
[0]PETSC ERROR: ./matmodify on a main-debug named EN438880 by 
acro018 Tue May 28 13:33:46 2024
[0]PETSC ERROR: Configure options --with-x --download-hdf5 
--download-zlib --download-netcdf --download-pnetcdf 
--download-exodusii --download-triangle --download-ptscotch 
--download-chaco --download-hypre
[0]PETSC ERROR: #1 MatSetValuesBlocked_MPIBAIJ() at 
/home/acro018/software/PETSc/code/src/mat/impls/baij/mpi/mpibaij.c:448
[0]PETSC ERROR: #2 MatSetValuesBlocked() at 
/home/acro018/software/PETSc/code/src/mat/interface/matrix.c:2030
[0]PETSC ERROR: - Error Message 
--
[0]PETSC ERROR:   It appears a new error in the code was triggered 
after a previous error, possibly because:
[0]PETSC ERROR:   -  The first error was not properly handled via 
(for example) the use of

[0]PETSC ERROR: PetscCall(TheFunctionThatErrors()); or
[0]PETSC ERROR:   -  The second error was triggered while handling 
the first error.
[0]PETSC ERROR:   Above is the traceback for the previous unhandled 
error, below the traceback for the next error
[0]PETSC ERROR:   ALL ERRORS in the PETSc libraries are fatal, you 
should add the appropriate error checking to the code

[0]PETSC ERROR: Petsc has generated inconsistent data
[0]PETSC ERROR: Incorrect colmap
[0]PETSC ERROR: See https://urldefense.us/v3/__https://petsc.org/release/faq/__;!!G_uCfscf7eWS!bEqLWafwqHADoOBGkMAUu-5jnxuJg95froHeBJtxmkLCGFqP1MFGgnLKOgYILALF3hYWj_cgiKOh5rEeun-2sybwJvu9iBaG$  for trouble 
shooting.

[0]PETSC ERROR: Petsc Release Version 3.21.1, unknown
[0]PETSC ERROR: ./matmodify on a main-debug named EN438880 by 
acro018 Tue May 28 13:33:46 2024
[0]PETSC ERROR: Configure options --with-x --download-hdf5 
--download-zlib --download-netcdf --download-pnetcdf 
--download-exodusii --download-triangle --download-ptscotch 
--download-chaco --download-hypre
[0]PETSC ERROR: #1 MatSetValuesBlocked_MPIBAIJ() at 
/home/acro018/software/PETSc/code/src/mat/impls/baij/mpi/mpibaij.c:448
[0]PETSC ERROR: #2 MatAssemblyEnd_MPIBAIJ() at 
/home/acro018/software/PETSc/code/src/mat/impls/baij/mpi/mpibaij.c:943
[0]PETSC ERROR: #3 MatAssemblyEnd() at 
/home/acro018/software/PETSc/code/src/mat/interface/matrix.c:5820

[0]PETSC ERROR: #4 matmodify.F90:40



  Barry


On May 26, 2024, at 10:45 PM, Adrian Croucher 
 wrote:


hi,

I've been trying

Re: [petsc-users] Modify matrix nonzero structure

[Adrian Croucher]

Hmm, that's a bit weird. I haven't modified the test code - I checked 
the file date to make sure. I also tried deleting my PETSc build dir and 
rebuilding it, then rebuilding the test code. I still get the error if I 
run on 2 ranks with -dm_mat_type (or -mat_type) baij or mpibaij, but 
it's fine on aij or mpiaij.


My actual (non-test) code is however working ok now, in serial or 
parallel. So I don't think this should hold up merging your bugfix.


- Adrian

On 28/05/24 2:13 pm, Barry Smith wrote:
  When I run the exact code you sent with two ranks and—mat_type 
mpibaij, it runs as expected. If you modified the code in any way to 
demonstrate the bug, please send the modified code.




On May 27, 2024, at 9:37 PM, Adrian Croucher 
 wrote:


hi Barry,

On 28/05/24 7:46 am, Barry Smith wrote:


   Thanks for reporting this. It is a bug. I have a fixed branch 
*barry/2024-05-27/fix-bug-baij-setvaluesblocked/release * and 
associated merge request 
https://urldefense.us/v3/__https://gitlab.com/petsc/petsc/-/merge_requests/7578__;!!G_uCfscf7eWS!dtDmIxQzHvnehLcuzdtJxWThxViUdu26r8fjum7PCvvKQp2PyZxpEZ0D1gvG1OsEvRs-sAisLmK-51Z7g4GONJEssSbd8xUJ$ 

Thanks very much, that does appear to fix the bug when I run my test 
program in serial.


If I run it on 2 processes, it is OK with AIJ matrix type, but with 
BAIJ I get an error (see below). Is there another problem, or I am 
doing something else wrong?


- Adrian

[0]PETSC ERROR: - Error Message 
--

[0]PETSC ERROR: Petsc has generated inconsistent data
[0]PETSC ERROR: Incorrect colmap
[0]PETSC ERROR: See 
https://urldefense.us/v3/__https://petsc.org/release/faq/__;!!G_uCfscf7eWS!dtDmIxQzHvnehLcuzdtJxWThxViUdu26r8fjum7PCvvKQp2PyZxpEZ0D1gvG1OsEvRs-sAisLmK-51Z7g4GONJEssfaywInx$
  for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.21.1, unknown
[0]PETSC ERROR: ./matmodify on a main-debug named EN438880 by acro018 
Tue May 28 13:33:46 2024
[0]PETSC ERROR: Configure options --with-x --download-hdf5 
--download-zlib --download-netcdf --download-pnetcdf 
--download-exodusii --download-triangle --download-ptscotch 
--download-chaco --download-hypre
[0]PETSC ERROR: #1 MatSetValuesBlocked_MPIBAIJ() at 
/home/acro018/software/PETSc/code/src/mat/impls/baij/mpi/mpibaij.c:448
[0]PETSC ERROR: #2 MatSetValuesBlocked() at 
/home/acro018/software/PETSc/code/src/mat/interface/matrix.c:2030
[0]PETSC ERROR: - Error Message 
--
[0]PETSC ERROR:   It appears a new error in the code was triggered 
after a previous error, possibly because:
[0]PETSC ERROR:   -  The first error was not properly handled via 
(for example) the use of

[0]PETSC ERROR: PetscCall(TheFunctionThatErrors()); or
[0]PETSC ERROR:   -  The second error was triggered while handling 
the first error.
[0]PETSC ERROR:   Above is the traceback for the previous unhandled 
error, below the traceback for the next error
[0]PETSC ERROR:   ALL ERRORS in the PETSc libraries are fatal, you 
should add the appropriate error checking to the code

[0]PETSC ERROR: Petsc has generated inconsistent data
[0]PETSC ERROR: Incorrect colmap
[0]PETSC ERROR: See 
https://urldefense.us/v3/__https://petsc.org/release/faq/__;!!G_uCfscf7eWS!dtDmIxQzHvnehLcuzdtJxWThxViUdu26r8fjum7PCvvKQp2PyZxpEZ0D1gvG1OsEvRs-sAisLmK-51Z7g4GONJEssfaywInx$
  for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.21.1, unknown
[0]PETSC ERROR: ./matmodify on a main-debug named EN438880 by acro018 
Tue May 28 13:33:46 2024
[0]PETSC ERROR: Configure options --with-x --download-hdf5 
--download-zlib --download-netcdf --download-pnetcdf 
--download-exodusii --download-triangle --download-ptscotch 
--download-chaco --download-hypre
[0]PETSC ERROR: #1 MatSetValuesBlocked_MPIBAIJ() at 
/home/acro018/software/PETSc/code/src/mat/impls/baij/mpi/mpibaij.c:448
[0]PETSC ERROR: #2 MatAssemblyEnd_MPIBAIJ() at 
/home/acro018/software/PETSc/code/src/mat/impls/baij/mpi/mpibaij.c:943
[0]PETSC ERROR: #3 MatAssemblyEnd() at 
/home/acro018/software/PETSc/code/src/mat/interface/matrix.c:5820

[0]PETSC ERROR: #4 matmodify.F90:40



  Barry


On May 26, 2024, at 10:45 PM, Adrian Croucher 
 wrote:


hi,

I've been trying creating a matrix with DMCreateMatrix() and then 
adding extra blocks of nonzeros into it using 
MatSetValuesBlocked(), but getting some unexpected results if I set 
the matrix type to BAIJ. It seems to behave as expected if I use 
matrix type AIJ.


I've attached a minimal example program. It reads in the DMPlex 
from file, sets up a section on it, creates a matrix (blocksize 2) 
and then inserts a single 2x2 block at global block indices (0,7). 
It views the matrix before and after the insertion.


If I run with "-dm_mat_type aij" it gives the expected results, but 
with "-dm_mat_type baij" it doesn't - e.g. if run in serial, it 
adds the new nonzeros in the right place but also adds a

Re: [petsc-users] Modify matrix nonzero structure

[Adrian Croucher]

hi Barry,

On 28/05/24 7:46 am, Barry Smith wrote:


   Thanks for reporting this. It is a bug. I have a fixed branch 
*barry/2024-05-27/fix-bug-baij-setvaluesblocked/release * and 
associated merge request 
https://urldefense.us/v3/__https://gitlab.com/petsc/petsc/-/merge_requests/7578__;!!G_uCfscf7eWS!e6Dvthc2K4dBRQ77JiOwCRZkPap2KZithO5HtDmQtF_1UGlC1EcRKJgoykBF54djkl1-_gepUhj95-78_S9dbu-68mw6cemp$ 

Thanks very much, that does appear to fix the bug when I run my test 
program in serial.


If I run it on 2 processes, it is OK with AIJ matrix type, but with BAIJ 
I get an error (see below). Is there another problem, or I am doing 
something else wrong?


- Adrian

[0]PETSC ERROR: - Error Message 
--

[0]PETSC ERROR: Petsc has generated inconsistent data
[0]PETSC ERROR: Incorrect colmap
[0]PETSC ERROR: See 
https://urldefense.us/v3/__https://petsc.org/release/faq/__;!!G_uCfscf7eWS!e6Dvthc2K4dBRQ77JiOwCRZkPap2KZithO5HtDmQtF_1UGlC1EcRKJgoykBF54djkl1-_gepUhj95-78_S9dbu-68p6qYDHU$
  for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.21.1, unknown
[0]PETSC ERROR: ./matmodify on a main-debug named EN438880 by acro018 
Tue May 28 13:33:46 2024
[0]PETSC ERROR: Configure options --with-x --download-hdf5 
--download-zlib --download-netcdf --download-pnetcdf --download-exodusii 
--download-triangle --download-ptscotch --download-chaco --download-hypre
[0]PETSC ERROR: #1 MatSetValuesBlocked_MPIBAIJ() at 
/home/acro018/software/PETSc/code/src/mat/impls/baij/mpi/mpibaij.c:448
[0]PETSC ERROR: #2 MatSetValuesBlocked() at 
/home/acro018/software/PETSc/code/src/mat/interface/matrix.c:2030
[0]PETSC ERROR: - Error Message 
--
[0]PETSC ERROR:   It appears a new error in the code was triggered after 
a previous error, possibly because:
[0]PETSC ERROR:   -  The first error was not properly handled via (for 
example) the use of

[0]PETSC ERROR:  PetscCall(TheFunctionThatErrors()); or
[0]PETSC ERROR:   -  The second error was triggered while handling the 
first error.
[0]PETSC ERROR:   Above is the traceback for the previous unhandled 
error, below the traceback for the next error
[0]PETSC ERROR:   ALL ERRORS in the PETSc libraries are fatal, you 
should add the appropriate error checking to the code

[0]PETSC ERROR: Petsc has generated inconsistent data
[0]PETSC ERROR: Incorrect colmap
[0]PETSC ERROR: See 
https://urldefense.us/v3/__https://petsc.org/release/faq/__;!!G_uCfscf7eWS!e6Dvthc2K4dBRQ77JiOwCRZkPap2KZithO5HtDmQtF_1UGlC1EcRKJgoykBF54djkl1-_gepUhj95-78_S9dbu-68p6qYDHU$
  for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.21.1, unknown
[0]PETSC ERROR: ./matmodify on a main-debug named EN438880 by acro018 
Tue May 28 13:33:46 2024
[0]PETSC ERROR: Configure options --with-x --download-hdf5 
--download-zlib --download-netcdf --download-pnetcdf --download-exodusii 
--download-triangle --download-ptscotch --download-chaco --download-hypre
[0]PETSC ERROR: #1 MatSetValuesBlocked_MPIBAIJ() at 
/home/acro018/software/PETSc/code/src/mat/impls/baij/mpi/mpibaij.c:448
[0]PETSC ERROR: #2 MatAssemblyEnd_MPIBAIJ() at 
/home/acro018/software/PETSc/code/src/mat/impls/baij/mpi/mpibaij.c:943
[0]PETSC ERROR: #3 MatAssemblyEnd() at 
/home/acro018/software/PETSc/code/src/mat/interface/matrix.c:5820

[0]PETSC ERROR: #4 matmodify.F90:40



  Barry


On May 26, 2024, at 10:45 PM, Adrian Croucher 
 wrote:


hi,

I've been trying creating a matrix with DMCreateMatrix() and then 
adding extra blocks of nonzeros into it using MatSetValuesBlocked(), 
but getting some unexpected results if I set the matrix type to BAIJ. 
It seems to behave as expected if I use matrix type AIJ.


I've attached a minimal example program. It reads in the DMPlex from 
file, sets up a section on it, creates a matrix (blocksize 2) and 
then inserts a single 2x2 block at global block indices (0,7). It 
views the matrix before and after the insertion.


If I run with "-dm_mat_type aij" it gives the expected results, but 
with "-dm_mat_type baij" it doesn't - e.g. if run in serial, it adds 
the new nonzeros in the right place but also adds a whole lot of 
other duplicated entries in block row 0.


Is there something I'm not understanding about BAIJ, or about 
MatSetValuesBlocked()? or possibly some other mistake?


- Adrian

On 20/05/24 12:24 pm, Barry Smith wrote:


   You can call MatSetOption(mat,MAT_NEW_NONZERO_LOCATION_ERR) then 
insert the new values. If it is just a handful of new insertions the 
extra time should be small.


    Making a copy of the matrix won't give you a new matrix that is 
any faster to insert into so best to just use the same matrix.


  Barry


On May 19, 2024, at 7:44 PM, Adrian Croucher 
 wrote:


This Message Is From an External Sender
This message came from outside your organization.
hi,

I have a Jacobian matrix created using DMCreate

Re: [petsc-users] Modify matrix nonzero structure

[Adrian Croucher]

hi,

I've been trying creating a matrix with DMCreateMatrix() and then adding 
extra blocks of nonzeros into it using MatSetValuesBlocked(), but 
getting some unexpected results if I set the matrix type to BAIJ. It 
seems to behave as expected if I use matrix type AIJ.


I've attached a minimal example program. It reads in the DMPlex from 
file, sets up a section on it, creates a matrix (blocksize 2) and then 
inserts a single 2x2 block at global block indices (0,7). It views the 
matrix before and after the insertion.


If I run with "-dm_mat_type aij" it gives the expected results, but with 
"-dm_mat_type baij" it doesn't - e.g. if run in serial, it adds the new 
nonzeros in the right place but also adds a whole lot of other 
duplicated entries in block row 0.


Is there something I'm not understanding about BAIJ, or about 
MatSetValuesBlocked()? or possibly some other mistake?


- Adrian

On 20/05/24 12:24 pm, Barry Smith wrote:


   You can call MatSetOption(mat,MAT_NEW_NONZERO_LOCATION_ERR) then 
insert the new values. If it is just a handful of new insertions the 
extra time should be small.


    Making a copy of the matrix won't give you a new matrix that is 
any faster to insert into so best to just use the same matrix.


  Barry


On May 19, 2024, at 7:44 PM, Adrian Croucher 
 wrote:


This Message Is From an External Sender
This message came from outside your organization.
hi,

I have a Jacobian matrix created using DMCreateMatrix(). What would be
the best way to add extra nonzero entries into it?

I'm guessing that DMCreateMatrix() allocates the storage so the nonzero
structure can't really be easily modified. Would it be a case of
creating a new matrix, copying the nonzero entries from the original one
and then adding the extra ones, before calling MatSetUp() or similar? If
so, how exactly would you copy the nonzero structure from the original
matrix?

Background: the flow problem I'm solving (on a DMPlex with finite volume
method) has complex source terms that depend on the solution (e.g.
pressure), and can also depend on other source terms. A simple example
is when fluid is extracted from one location, with a pressure-dependent
flow rate, and some of it is then reinjected in another location. This
can result in poor nonlinear solver convergence. I think the reason is
that there are effectively missing Jacobian entries in the row for the
reinjection cell, which should have an additional dependence on the
solution in the cell where fluid is extracted.

- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
Waipapa Taumata Rau / University of Auckland, New Zealand
email:a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611





--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
Waipapa Taumata Rau / University of Auckland, New Zealand
email:a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611


3x3grid.exo
Description: Binary data
program matmodify

#include 

  use petsc
  implicit none

  PetscInt, parameter :: overlap = 1
  DM :: dm, dist_dm
  PetscErrorCode :: ierr
  PetscSection :: section
  PetscReal, parameter :: vals(4) = [-1.d0, -2.d0, -3.d0, -4.d0]
  Mat :: J

  call PetscInitialize(PETSC_NULL_CHARACTER, ierr); CHKERRQ(ierr)

  call DMPlexCreateFromFile(PETSC_COMM_WORLD, "3x3grid.exo", "mesh", &
   PETSC_TRUE, dm, ierr); CHKERRQ(ierr)
  call DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE, ierr)
  call DMSetFromOptions(dm, ierr)

  call DMPlexDistribute(dm, overlap, PETSC_NULL_SF, dist_dm, ierr)
  CHKERRQ(ierr)
  if (dist_dm .ne. PETSC_NULL_DM) then
 call DMDestroy(dm, ierr); CHKERRQ(ierr)
 dm = dist_dm
  end if

  call dm_set_fields(dm, num_components = [1,1])
  section = dm_create_section(dm, num_components = [1,1], field_dim = [3,3])
  call DMSetSection(dm, section, ierr); CHKERRQ(ierr)

  call DMCreateMatrix(dm, J, ierr)
  call MatView(J, PETSC_VIEWER_STDOUT_WORLD, ierr)

  call MatSetOption(J, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE, ierr)
  call MatSetValuesBlocked(J, 1, [0], 1, [7], vals, INSERT_VALUES, ierr)

  call MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY, ierr); CHKERRQ(ierr)
  call MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY, ierr); CHKERRQ(ierr)

  call MatView(J, PETSC_VIEWER_STDOUT_WORLD, ierr)

  call MatDestroy(J, ierr)
  call DMDestroy(dm, ierr); CHKERRQ(ierr)
  call PetscFinalize(ierr); CHKERRQ(ierr)

contains

  subroutine dm_set_fields(dm, num_components)
!! Sets fields in the DM.

DM, intent(in out) :: dm
PetscInt, intent(in) :: num_components(:) !! Number of components in each field
! Locals:
PetscInt :: dim, f
PetscFV :: fvm
PetscErrorCode :: ierr

call DMGetDimension(dm, dim, ierr); CHKERRQ(ierr)
call DMClearFields(dm, ierr); CHKERRQ(ierr)

associate(num_fields => size(num_components))
  do f = 1, num_fields
 call PetscFVCreate(PETSC_COMM_WORLD, fvm, ierr); CHKERRQ(ierr)
 call PetscFVSetFrom

Re: [petsc-users] Modify matrix nonzero structure

[Adrian Croucher]

Great, it sounds like this might be easier than I expected. Thanks very 
much.


Did you have any thoughts on my diagnosis of the problem (the poor 
nonlinear solver convergence being caused by missing Jacobian elements 
representing interaction between the sources)?


- Adrian

On 20/05/24 12:41 pm, Matthew Knepley wrote:

On Sun, May 19, 2024 at 8:25 PM Barry Smith  wrote:

You can call MatSetOption(mat,MAT_NEW_NONZERO_LOCATION_ERR) then
insert the new values. If it is just a handful of new insertions
the extra time should be small. Making a copy of the matrix won't
give you a new matrix that is any faster to
ZjQcmQRYFpfptBannerStart
This Message Is From an External Sender
This message came from outside your organization.
ZjQcmQRYFpfptBannerEnd

   You can call MatSetOption(mat,MAT_NEW_NONZERO_LOCATION_ERR)
then insert the new values. If it is just a handful of new
insertions the extra time should be small.

    Making a copy of the matrix won't give you a new matrix that
is any faster to insert into so best to just use the same matrix.


Let me add to Barry's answer. The preallocation infrastructure is now 
not strictly necessary. It is possible to just add all your nonzeros 
in and assembly,  and the performance will be pretty good (uses 
hashing etc). So if just adding a few nonzeros does not work, we can 
go this route.


  Thanks,

     Matt

  Barry



On May 19, 2024, at 7:44 PM, Adrian Croucher
 wrote:

This Message Is From an External Sender
This message came from outside your organization.
hi,

I have a Jacobian matrix created using DMCreateMatrix(). What would be
the best way to add extra nonzero entries into it?

I'm guessing that DMCreateMatrix() allocates the storage so the nonzero
structure can't really be easily modified. Would it be a case of
creating a new matrix, copying the nonzero entries from the original one
and then adding the extra ones, before calling MatSetUp() or similar? If
so, how exactly would you copy the nonzero structure from the original
matrix?

Background: the flow problem I'm solving (on a DMPlex with finite volume
method) has complex source terms that depend on the solution (e.g.
pressure), and can also depend on other source terms. A simple example
is when fluid is extracted from one location, with a pressure-dependent
flow rate, and some of it is then reinjected in another location. This
can result in poor nonlinear solver convergence. I think the reason is
that there are effectively missing Jacobian entries in the row for the
reinjection cell, which should have an additional dependence on the
solution in the cell where fluid is extracted.

- Adrian




--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
Waipapa Taumata Rau / University of Auckland, New Zealand
email:a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

[petsc-users] Modify matrix nonzero structure

[Adrian Croucher]

 hi, I have a Jacobian matrix created using DMCreateMatrix(). What would be the best way to add extra nonzero entries into it? I'm guessing that DMCreateMatrix() allocates the storage so the nonzero structure can't really be easily modified. 




ZjQcmQRYFpfptBannerStart




  

  
	This Message Is From an External Sender
  
  
This message came from outside your organization.
  



 
  


ZjQcmQRYFpfptBannerEnd




hi,

I have a Jacobian matrix created using DMCreateMatrix(). What would be 
the best way to add extra nonzero entries into it?

I'm guessing that DMCreateMatrix() allocates the storage so the nonzero 
structure can't really be easily modified. Would it be a case of 
creating a new matrix, copying the nonzero entries from the original one 
and then adding the extra ones, before calling MatSetUp() or similar? If 
so, how exactly would you copy the nonzero structure from the original 
matrix?

Background: the flow problem I'm solving (on a DMPlex with finite volume 
method) has complex source terms that depend on the solution (e.g. 
pressure), and can also depend on other source terms. A simple example 
is when fluid is extracted from one location, with a pressure-dependent 
flow rate, and some of it is then reinjected in another location. This 
can result in poor nonlinear solver convergence. I think the reason is 
that there are effectively missing Jacobian entries in the row for the 
reinjection cell, which should have an additional dependence on the 
solution in the cell where fluid is extracted.

- Adrian

-- 
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
Waipapa Taumata Rau / University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] GMSH entities

[Adrian Croucher]

hi Matt,

On 15/05/24 6:05 am, Matthew Knepley wrote:

On Tue, May 14, 2024 at 9:07 AM Matthew Knepley  wrote:

On Mon, May 13, 2024 at 10:04 PM Adrian Croucher
 wrote:

On 14/05/24 1:44 pm, Matthew Knepley wrote:


I wish GMsh was clearer about what is optional:
https://urldefense.us/v3/__https://gmsh.info/doc/texinfo/gmsh.html*MSH-file-format__;Iw!!G_uCfscf7eWS!dbZ8R9PicaYRx4Y2IkPxvCoLlJ4lpsMbIAFBKbKZK5h5-2OTX5Ne-AwgOBHv2mJaHqg9uo5orDQYs2PaoAcpDmS65Au1g9uw$ 
<https://urldefense.us/v3/__https://gmsh.info/doc/texinfo/gmsh.html*MSH-file-format__;Iw!!G_uCfscf7eWS!dbZ8R9PicaYRx4Y2IkPxvCoLlJ4lpsMbIAFBKbKZK5h5-2OTX5Ne-AwgOBHv2mJaHqg9uo5orDQYs2PaoAcpDmS65Au1g9uw$ >


They do talk about it, but not exhaustively. GMsh always
writes and $Entities block from what I can tell.
I can make it optional, it just might take until after the
PETSc Meeting.


Looks like $Entities are optional:

https://urldefense.us/v3/__https://gitlab.onelab.info/gmsh/gmsh/-/commit/b5feba2af57181ffa946d3f0c494b014603c6efa__;!!G_uCfscf7eWS!dbZ8R9PicaYRx4Y2IkPxvCoLlJ4lpsMbIAFBKbKZK5h5-2OTX5Ne-AwgOBHv2mJaHqg9uo5orDQYs2PaoAcpDmS65ISaNx6E$ 
<https://urldefense.us/v3/__https://gitlab.onelab.info/gmsh/gmsh/-/commit/b5feba2af57181ffa946d3f0c494b014603c6efa__;!!G_uCfscf7eWS!dbZ8R9PicaYRx4Y2IkPxvCoLlJ4lpsMbIAFBKbKZK5h5-2OTX5Ne-AwgOBHv2mJaHqg9uo5orDQYs2PaoAcpDmS65ISaNx6E$ >


I can also load a GMSH 4.1 file without $Entities into GMSH
itself and it doesn't complain, suggesting that they are
indeed optional.

Yes, but they are not careful to specify when a file can be
inconsistent. For instance, omitting the $Entities, but then
specifying entity numbers in the $Nodes block. I think they also
thought this was inconsistent, but then got user complaints. The
minimal example they show does exactly this.

If the $Entities aren't strictly needed for anything in DMPlex
(which I'm guessing they aren't, as the GMSH file format 2.2
doesn't even have them) then it would be useful not to require
them.

I put in some code for this:
https://urldefense.us/v3/__https://gitlab.com/petsc/petsc/-/merge_requests/7546__;!!G_uCfscf7eWS!dbZ8R9PicaYRx4Y2IkPxvCoLlJ4lpsMbIAFBKbKZK5h5-2OTX5Ne-AwgOBHv2mJaHqg9uo5orDQYs2PaoAcpDmS65MnfH9Hr$ 
<https://urldefense.us/v3/__https://gitlab.com/petsc/petsc/-/merge_requests/7546__;!!G_uCfscf7eWS!dbZ8R9PicaYRx4Y2IkPxvCoLlJ4lpsMbIAFBKbKZK5h5-2OTX5Ne-AwgOBHv2mJaHqg9uo5orDQYs2PaoAcpDmS65MnfH9Hr$ >


It just ignores entity numbers when there is no section.


This merged, so now it should be fixed for you.


Thanks, that seems to fix the problem. Great!

Cheers, Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
Waipapa Taumata Rau / University of Auckland, New Zealand
email:a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] GMSH entities

[Adrian Croucher]

On 14/05/24 1:44 pm, Matthew Knepley wrote:

I wish GMsh was clearer about what is optional: 
https://urldefense.us/v3/__https://gmsh.info/doc/texinfo/gmsh.html*MSH-file-format__;Iw!!G_uCfscf7eWS!aCypJMAzwWHAJLXGNJmSthDjbHcU-8_MdsaXZ4d1r1RKyL0bqIv5ZuLmQtV6ve4XcTjapf38-bsdDOLEDAlhWoaaGpmjDCYs$  
<https://urldefense.us/v3/__https://gmsh.info/doc/texinfo/gmsh.html*MSH-file-format__;Iw!!G_uCfscf7eWS!aCypJMAzwWHAJLXGNJmSthDjbHcU-8_MdsaXZ4d1r1RKyL0bqIv5ZuLmQtV6ve4XcTjapf38-bsdDOLEDAlhWoaaGpmjDCYs$ > 

They do talk about it, but not exhaustively. GMsh always writes and 
$Entities block from what I can tell.

I can make it optional, it just might take until after the PETSc Meeting.


Looks like $Entities are optional:

https://urldefense.us/v3/__https://gitlab.onelab.info/gmsh/gmsh/-/commit/b5feba2af57181ffa946d3f0c494b014603c6efa__;!!G_uCfscf7eWS!aCypJMAzwWHAJLXGNJmSthDjbHcU-8_MdsaXZ4d1r1RKyL0bqIv5ZuLmQtV6ve4XcTjapf38-bsdDOLEDAlhWoaaGjNKbSSU$ 

I can also load a GMSH 4.1 file without $Entities into GMSH itself and 
it doesn't complain, suggesting that they are indeed optional.


If the $Entities aren't strictly needed for anything in DMPlex (which 
I'm guessing they aren't, as the GMSH file format 2.2 doesn't even have 
them) then it would be useful not to require them.


- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
Waipapa Taumata Rau / University of Auckland, New Zealand
email:a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

[petsc-users] GMSH entities

[Adrian Croucher]

 hi, We often create meshes in GMSH format using the meshio library. This works OK if we stick to GMSH file format 2. 2. If we use GMSH file format 4. 1, DMPlex can't read them because it expects the "Entities" section to be present: [0]PETSC




ZjQcmQRYFpfptBannerStart




  

  
	This Message Is From an External Sender
  
  
This message came from outside your organization.
  



 
  


ZjQcmQRYFpfptBannerEnd




hi,

We often create meshes in GMSH format using the meshio library. This 
works OK if we stick to GMSH file format 2.2.

If we use GMSH file format 4.1, DMPlex can't read them because it 
expects the "Entities" section to be present:

[0]PETSC ERROR: Unexpected data in file
[0]PETSC ERROR: File is not a valid Gmsh file, expecting $Entities, not 
$Nodes
[0]PETSC ERROR: See https://urldefense.us/v3/__https://petsc.org/release/faq/__;!!G_uCfscf7eWS!ZfBS1KM5EBZ7ZJIu6lBKFcclVMmXteXsW8m9HBEZ5tIf0u_3duEFt9eXKF7FcorQSAQqD5SbJrYh4C8rX676S_IMI_sp6naX$ for trouble shooting.
[0]PETSC ERROR: Petsc Development GIT revision: 
v3.21.1-124-g2d06e2faec8  GIT Date: 2024-05-08 19:31:33 +
[0]PETSC ERROR: waiwera on a main-debug named EN438880 by acro018 Tue 
May 14 11:25:54 2024
[0]PETSC ERROR: Configure options --with-x --download-hdf5 
--download-zlib --download-netcdf --download-pnetcdf --download-exodusii 
--download-triangle --download-ptscotch --download-chaco --download-hypre
[0]PETSC ERROR: #1 GmshExpect() at 
/home/acro018/software/PETSc/code/src/dm/impls/plex/plexgmsh.c:270
[0]PETSC ERROR: #2 DMPlexCreateGmsh() at 
/home/acro018/software/PETSc/code/src/dm/impls/plex/plexgmsh.c:1608
[0]PETSC ERROR: #3 DMPlexCreateGmshFromFile() at 
/home/acro018/software/PETSc/code/src/dm/impls/plex/plexgmsh.c:1469
[0]PETSC ERROR: #4 DMPlexCreateFromFile() at 
/home/acro018/software/PETSc/code/src/dm/impls/plex/plexcreate.c:5804

By default meshio doesn't seem to write the Entities section. From what 
I can gather, it is optional.

Am I right that this section is not optional in DMPlex?

- Adrian

-- 
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
Waipapa Taumata Rau / University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] HDF5 time step count

[Adrian Croucher]

hi Matt,

On 11/05/24 4:12 am, Matthew Knepley wrote:


Thanks. I tried it out and the error below was raised, looks like
it's when it tries to destroy the viewer. It still runs ok when
DMGetOutputSequenceLength() isn't called.

Sorry, it looks like I did not close the group. I pushed an update. If 
that is alright, I will merge it.


Thanks, it does appear to work now. That's going to be very useful!

- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
Waipapa Taumata Rau / University of Auckland, New Zealand
email:a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] HDF5 time step count

[Adrian Croucher]

hi Matt,

On 10/05/24 12:15 am, Matthew Knepley wrote:


I just tried to test it, but there doesn't seem to be a Fortran
interface for DMGetOutputSequenceLength().

Pushed.


Thanks. I tried it out and the error below was raised, looks like it's 
when it tries to destroy the viewer. It still runs ok when 
DMGetOutputSequenceLength() isn't called.


- Adrian

[0]PETSC ERROR: Error in external library
[0]PETSC ERROR: Error in HDF5 call H5Fclose() Status -1
[0]PETSC ERROR: See 
https://urldefense.us/v3/__https://petsc.org/release/faq/__;!!G_uCfscf7eWS!YwgL4tYSuYus7pShvp7QSDY_SmHzGbe32BN2ELvCuZw7P6ZOTBcfFsRiN3G0zVuU29FqKaNqYShrLLflVLcbn25Zbw8CAsGo$
  for trouble shooting.
[0]PETSC ERROR: Petsc Development GIT revision: 
v3.21.1-124-g2d06e2faec8  GIT Date: 2024-05-08 19:31:33 +
[0]PETSC ERROR: ./initial_test on a main-debug named EN438880 by acro018 
Fri May 10 16:55:34 2024
[0]PETSC ERROR: Configure options --with-x --download-hdf5 
--download-zlib --download-netcdf --download-pnetcdf --download-exodusii 
--download-triangle --download-ptscotch --download-chaco --download-hypre
[0]PETSC ERROR: #1 PetscViewerFileClose_HDF5() at 
/home/acro018/software/PETSc/code/src/sys/classes/viewer/impls/hdf5/hdf5v.c:107
[0]PETSC ERROR: #2 PetscViewerDestroy_HDF5() at 
/home/acro018/software/PETSc/code/src/sys/classes/viewer/impls/hdf5/hdf5v.c:126
[0]PETSC ERROR: #3 PetscViewerDestroy() at 
/home/acro018/software/PETSc/code/src/sys/classes/viewer/interface/view.c:101

[0]PETSC ERROR: #4 ../src/initial.F90:497

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
Waipapa Taumata Rau / University of Auckland, New Zealand
email:a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] HDF5 time step count

[Adrian Croucher]

hi Matt,

On 9/05/24 3:00 am, Matthew Knepley wrote:

Sorry about the delay. I had lost track of this. Can you look at branch

  knepley/feature-hdf5-seq-len

I have not made a test yet, but if this works for you, I will make a 
test and merge it in.


Thanks for looking at that.

I just tried to test it, but there doesn't seem to be a Fortran 
interface for DMGetOutputSequenceLength().


Another odd thing I found was that DMPlexConstructGhostCells() seemed to 
complain if I passed in PETSC_NULL_INTEGER for the parameter 
numGhostCells (which works ok in v3.21.1) - the compiler complained 
about a rank-1/scalar mismatch. Can you no longer pass in a null integer 
for that? It works if I declare a dummy integer and pass it in.


- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
Waipapa Taumata Rau / University of Auckland, New Zealand
email:a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] PETSc options

[Adrian Croucher]

On 7/05/24 1:11 am, Barry Smith wrote:



   Yes, just call PetscOptionsHasName() with -v and ignore the result.


Thanks Barry, that does the trick. Yes, I was using PETSc <3.19 before 
which is why it didn't show up until now.


- Adrian



On Mon, May 6, 2024 at 1:04 AM Adrian Croucher 
 wrote:


hi,

My code has some optional command line arguments -v and -h for output of
version number and usage help. These are processed using Fortran's
get_command_argument().

Since updating PETSc to version 3.21, I get some extra warnings after
the output:

acro018@EN438880:~$ waiwera -v
1.5.0b1
WARNING! There are options you set that were not used!
WARNING! could be spelling mistake, etc!
There is one unused database option. It is:
Option left: name:-v (no value) source: command line

That didn't used to happen. What should I do to make them go away?

Regards, Adrian





--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
Waipapa Taumata Rau / University of Auckland, New Zealand
email:a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

[petsc-users] PETSc options

[Adrian Croucher]

 hi, My code has some optional command line arguments -v and -h for output of version number and usage help. These are processed using Fortran's get_command_argument(). Since updating PETSc to version 3. 21, I get some extra warnings after the




ZjQcmQRYFpfptBannerStart




  

  
	This Message Is From an External Sender
  
  
This message came from outside your organization.
  



 
  


ZjQcmQRYFpfptBannerEnd




hi,

My code has some optional command line arguments -v and -h for output of 
version number and usage help. These are processed using Fortran's 
get_command_argument().

Since updating PETSc to version 3.21, I get some extra warnings after 
the output:

acro018@EN438880:~$ waiwera -v
1.5.0b1
WARNING! There are options you set that were not used!
WARNING! could be spelling mistake, etc!
There is one unused database option. It is:
Option left: name:-v (no value) source: command line

That didn't used to happen. What should I do to make them go away?

Regards, Adrian

-- 
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
Waipapa Taumata Rau / University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] HDF5 time step count

[Adrian Croucher]

hi Matt & all,

I just had a query from one of my users which prompted me to see if any 
progress had been made on the issue below - using PETSc to get the 
number of time steps in an HDF5 file.


I can't see anything new in PETSc on this - I did try using 
PetscViewerHDF5ReadSizes() to see if that would do it, but it seems it 
doesn't. If I use that on the "time" dataset it just returns 1.


Regards, Adrian

On 11/10/21 2:08 pm, Adrian Croucher wrote:

On 10/11/21 11:59 AM, Matthew Knepley wrote:
On Sun, Oct 10, 2021 at 6:51 PM Adrian Croucher 
 wrote:


hi

Is there any way to query the PETSc HDF5 viewer to find the
number of
time steps in the file?

A common use case I have is that an HDF5 file from a previous
simulation
is used to get initial conditions for a subsequent run. The most
common
thing you want to do is restart from the last set of results in the
previous output. To do that you need to know how many time steps
there
are, so you can set the output index to be the last one.

I thought maybe I could just query the size of the "time"
dataset, but I
can't even see any obvious way to do that using the viewer functions.


There is nothing in there that does it right now. Do you know how to 
do it in HDF5?

If so, I can put it in. Otherwise, I will have to learn more HDF5 :)


I haven't actually tried this myself but it looks like what you do is:

1) get the dataspace for the dataset (in our case the "time" dataset):

hid_t dspace = H5Dget_space(dset);

2) Get the dimensions of the dataspace:

const int ndims = 1;

hsize_t dims[ndims];
H5Sget_simple_extent_dims(dspace, dims, NULL);

The first element of dims should be the number of time steps. Here 
I've assumed the number of dimensions of the time dataset is 1. In 
general you can instead query the rank of the dataspace using 
H5Sget_simple_extent_ndims() to get the rank ndims.


Regards, Adrian


--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
Waipapa Taumata Rau / University of Auckland, New Zealand
email:a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] Cray perftools

[Adrian Croucher]

hi Jed,

On 2/22/22 11:45, Jed Brown wrote:
If you can share before/after output from -log_view, it would likely 
help localize.


Unfortunately there aren't any "before" logs, because -log_view wasn't 
used while everything was going well.


It looks like it might be related to I/O partly based on comparing logs 
with those from runs on other machines (not very comparable ones, 
admittedly) and also experimenting with things like turning file (HDF5) 
output from the code off.




Another unintrusive thing (if you're allowed to run Linux perf) is to
Looks like that is not available on the cluster, but will check and see 
if some module can be loaded to get it.


On 2/22/22 12:10, Barry Smith wrote:


  Have any of the modules loaded changed? Or the shared libraries that are in 
any of the modules?
Will check. We build most of the important dependencies (e.g. 
blaslapack, netcdf, exodusii, hdf5, scotch, chaco) ourselves using 
--download-xxx in the PETSc build.


- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

[petsc-users] Cray perftools

[Adrian Croucher]

hi,

We have our PETSc-based code compiled on a Cray XC-50 machine, and it 
has just recently started running about 2.5 times slower on there. 
Neither the code nor PETSc has been recompiled lately.


Turning the PETSc logging on, it appears to be spending more time on I/O 
than it used to.


The cluster admins have suggested we rebuild with the Cray "perftools" 
module loaded to get profiling info. It's a slight hassle to rebuild 
everything, so I wondered, would this actually tell us anything that we 
don't already know from the PETSc logs?


- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] HDF5 timestepping in PETSc 3.16

[Adrian Croucher]

Any response on this?

This is a bit of a showstopper for me - I can't upgrade to PETSc 3.16 if 
it does not allow my users to read their HDF5 files created using 
earlier versions of PETSc.


So far I can't see a workaround. Possibly the timestepping functions 
need some kind of optional parameter to specify what the default 
timestepping attribute should be, if it's not present in the file 
(rather than just assuming it's false)?


Regards, Adrian

On 10/14/21 4:19 PM, Adrian Croucher wrote:

hi

I am just testing out PETSc 3.16 and making the necessary changes to 
my code. Amongst other things I now have to add a 
PetscViewerHDF5PushTimestepping() call before starting to output 
time-dependent results to HDF5 using a PetscViewer.


I now also have to add this call before reading in sets of previously 
computed time-dependent results (for restarting a simulation from the 
results of a previous run).


The problem with this is that if I try to read in the results of any 
previous run, computed with an earlier version of PETSc (< 3.16), an 
error is raised because the time-dependent datasets in the file do not 
have the 'timestepping' attribute.


Is there something else I need to do to make this work?

- Adrian


--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

[petsc-users] HDF5 timestepping in PETSc 3.16

[Adrian Croucher]

hi

I am just testing out PETSc 3.16 and making the necessary changes to my 
code. Amongst other things I now have to add a 
PetscViewerHDF5PushTimestepping() call before starting to output 
time-dependent results to HDF5 using a PetscViewer.


I now also have to add this call before reading in sets of previously 
computed time-dependent results (for restarting a simulation from the 
results of a previous run).


The problem with this is that if I try to read in the results of any 
previous run, computed with an earlier version of PETSc (< 3.16), an 
error is raised because the time-dependent datasets in the file do not 
have the 'timestepping' attribute.


Is there something else I need to do to make this work?

- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] HDF5 time step count

[Adrian Croucher]

hi Matt,

On 10/11/21 11:59 AM, Matthew Knepley wrote:
On Sun, Oct 10, 2021 at 6:51 PM Adrian Croucher 
mailto:a.crouc...@auckland.ac.nz>> wrote:


hi

Is there any way to query the PETSc HDF5 viewer to find the number of
time steps in the file?

A common use case I have is that an HDF5 file from a previous
simulation
is used to get initial conditions for a subsequent run. The most
common
thing you want to do is restart from the last set of results in the
previous output. To do that you need to know how many time steps
there
are, so you can set the output index to be the last one.

I thought maybe I could just query the size of the "time" dataset,
but I
can't even see any obvious way to do that using the viewer functions.


There is nothing in there that does it right now. Do you know how to 
do it in HDF5?

If so, I can put it in. Otherwise, I will have to learn more HDF5 :)


I haven't actually tried this myself but it looks like what you do is:

1) get the dataspace for the dataset (in our case the "time" dataset):

hid_t dspace = H5Dget_space(dset);

2) Get the dimensions of the dataspace:

const int ndims = 1;

hsize_t dims[ndims];
H5Sget_simple_extent_dims(dspace, dims, NULL);

The first element of dims should be the number of time steps. Here I've 
assumed the number of dimensions of the time dataset is 1. In general 
you can instead query the rank of the dataspace using 
H5Sget_simple_extent_ndims() to get the rank ndims.


Regards, Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

[petsc-users] HDF5 time step count

[Adrian Croucher]

hi

Is there any way to query the PETSc HDF5 viewer to find the number of 
time steps in the file?


A common use case I have is that an HDF5 file from a previous simulation 
is used to get initial conditions for a subsequent run. The most common 
thing you want to do is restart from the last set of results in the 
previous output. To do that you need to know how many time steps there 
are, so you can set the output index to be the last one.


I thought maybe I could just query the size of the "time" dataset, but I 
can't even see any obvious way to do that using the viewer functions.


Regards, Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] HDF5 corruption

[Adrian Croucher]

hi Jed,

It looked to me like a call to h5f_flush() is all that is required.

Some people said there would be a performance hit (maybe ~ 10% slower), 
which would be the trade-off for increased reliability. So if this were 
made available via PetscViewerFlush(), I'd probably make it optional in 
my code so the user could decide for themselves if it was worth it for them.


Do you think flushing would be a better option than closing/opening the 
file between writes?


Regards, Adrian

On 10/8/21 9:36 AM, Jed Brown wrote:

Adrian Croucher  writes:

> hi,
>
> One of the users of my PETSc-based code has reported that HDF5 output
> files can be corrupted and unusable if e.g. the run is killed. I've 
just
> done a bit of reading about this and it appears to be a known issue 
with

> HDF5.
>
> Some people suggest flushing the HDF5 file periodically to help prevent
> data loss. I had a look at PetscViewerFlush() but it doesn't seem to be
> implemented for the HDF5 viewer- is that correct?

Correct, but I think we can and should implement it. In your research 
just now, were there subtleties beyond this call?


https://portal.hdfgroup.org/display/HDF5/H5F_FLUSH 
<https://portal.hdfgroup.org/display/HDF5/H5F_FLUSH>


--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

[petsc-users] HDF5 corruption

[Adrian Croucher]

hi,

One of the users of my PETSc-based code has reported that HDF5 output 
files can be corrupted and unusable if e.g. the run is killed. I've just 
done a bit of reading about this and it appears to be a known issue with 
HDF5.


Some people suggest flushing the HDF5 file periodically to help prevent 
data loss. I had a look at PetscViewerFlush() but it doesn't seem to be 
implemented for the HDF5 viewer- is that correct?


I am currently calling PetscViewerHDF5Open() once at the start of the 
run and closing it at the end. Some other people suggest doing this 
before and after each write instead. Is there likely to be a significant 
performance penalty in doing that?


I gather HDF5 journaling has been promised for a while to get around 
this problem but as far as I can see it hasn't materialised yet...


Regards, Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] Boundary ghost cell types in PETSc 3.15

[Adrian Croucher]

On 30/06/21 11:17 pm, Matthew Knepley wrote:


Yes, src/dm/f90-mod/petscdm.h is missing PYRAMID. I will make a
fix right now.


Here it is: https://gitlab.com/petsc/petsc/-/merge_requests/4140 
<https://gitlab.com/petsc/petsc/-/merge_requests/4140>



Aha, yes that's fixed it.

Thanks Matt!

- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] flux vector

[Adrian Croucher]

hi

On 15/06/21 12:21 pm, Matthew Knepley wrote:


I admit that it is becoming complicated, but this can also be handled.

The corner case only shows up when we distribute with overlap. Thus we 
could,


1) Distribute the mesh without overlap

2) Mark the global boundary using the former algorithm

3) Distribute overlap = 1 using DMPlexDistributeOverlap(). 
DMPlexDistribute() just does these two steps at once if you ask for 
overlap.


4) Use the former algorithm, with the boundary label we just made, to 
mark flux faces.


That sounds like it might work. Do you think this perhaps slightly 
simpler version might also work?


1) distribute mesh without overlap

2) label boundary faces using DMPlexMarkBoundaryFaces()

3) DMPlexDistributeOverlap()

4) Un-label boundary faces with support size > 1 (i.e. the shared faces 
on the non-overlapped mesh- I think) - so the boundary label is now just 
the global boundary


5) loop over faces as in my original algorithm, labelling flux faces if 
they are either on the open boundary, or are not on the global boundary


Unfortunately I can't test this at the moment - it looks like the 
Fortran interface for DMPlexDistributeOverlap() is missing?


- Adrian

--

Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] flux vector

[Adrian Croucher]

hi Matt,

On 14/06/21 9:54 pm, Matthew Knepley wrote:
Okay, I think it is not so hard to get what you want in parallel. 
There are only two kinds of faces with supportSize == 1:


  a) Faces on the global boundary

  b) Faces which are "shared"


I think there is an unfortunate corner case in which faces can be both 
on the global boundary *and* shared.


For example: consider a square 2D mesh which is partitioned down the 
middle, so one process has the left half. With overlap = 1, a row of 
ghost cells will be added along the right-hand edge of this half of the 
mesh. The two faces at top and bottom of this row of ghost cells are on 
the global boundary and also shared.


With your algorithm these would be labelled as flux faces, but they 
shouldn't be (as they are on the global boundary). I can't see a way to 
eliminate those kinds of shared faces - can you?


- Adrian



It is the second set that is somewhat confusing because PetscSF does 
not have 2-sided information by default. However, it can make it.

There is a two-step check for "shared":

  1) Is the face in the PetscSF? Here you just check for it in the 
sorted "locals" array from PetscSFGetGraph()


  2) Is the face ghosted on another process? You can get this from 
PetscSFGetRootRanks().


I just wrote a small function to check for "shared" points. After 
that, I think you can just run


  1) After distribution, loop overall faces on current process

     If face on open boundary, label face as flux face

     else:

       if face has supportSize != 1 or (supportSize == 1 && shared), 
label face as flux face


  Thanks,

     Matt


--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] flux vector

[Adrian Croucher via petsc-users]

hi, thanks for the suggestions!

On 12/06/21 12:19 am, Matthew Knepley wrote:
However, using overlap = 1 put in a bunch of new faces. We do not care 
about the ones on the process boundary. They will be
handled by the other process. We do care about faces between two ghost 
cells, since they will be a false positive. Luckily, these

are labeled by the "ghost" label.


I think we *do* have to care about the faces on the process boundary 
(assuming by that you mean the faces with support size 1 on the outside 
edge of the partition ghost cells), because they can be ghosts of flux 
faces on other processes. If we ignore them they will have no dofs on 
the current process, but they could have dofs on another one. That 
inconsistency is the problem that causes the error when you create the 
global section.


Also, I don't think we want to exclude faces between partition ghost 
cells. Those faces will definitely be ghosts of a flux face on another 
process. Again, if we exclude them the dofs will not be consistent 
across processes. We might not actually compute a flux on those faces 
locally, but there has to be a space for them anyway.


I have found a simple algorithm now which seems to work in all my test 
cases, though I still wonder if there is a better way. The algorithm is:


1) label the global boundary faces before distribution (all faces with 
support size 1), using DMPlexMarkBoundaryFaces()


2) label any open boundary faces (on which Dirichlet BCs are applied) - 
I didn't mention these before, but they need to be included as flux faces


3) after distribution, loop over all faces on current process:

    if face on open boundary: label face as a flux face

    else:

  if face not on global boundary: label face as a flux face

Here the only test of support size is in step 1), which ensures that the 
dofs are always consistent between processes. Ultimately, the support 
size on the local process is not really relevant or reliable.


The main thing I don't like about this algorithm is that it relies on 
looping over all mesh faces in serial during step 1). I would rather not 
add to the serial part of the code and would prefer if it could all be 
done in parallel, after distribution. Maybe I'm worrying unnecessarily, 
and just getting the support size of each face is cheap enough that this 
won't be a significant issue?


- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

[petsc-users] flux vector

[Adrian Croucher]

hi

I have a DMPlex in a finite volume simulation and I want to create a Vec 
containing flux data only on the mesh faces where fluxes are computed 
(e.g. not on closed boundaries). Basically these are all the faces with 
support size 2 (aside from some complications related to boundary 
conditions).


To do this I clone the DM using DMClone() and create a section on it 
with DMPlexCreateSection(), passing in a label which has value 1 on all 
the required faces. The problem then consists of labelling the 
appropriate faces.


The tricky part is that when a face exists on more than one process, you 
have to make sure it has the same label value on all processes, 
otherwise an error results when you try and create a global vector on 
the DM. Because I distribute the DMPlex with overlap = 1, a face can 
have different support sizes on different processes.


So on each process, I need to label a face if it has support size 2 on 
*any* process.


Would it be possible to use PetscSFBcast and the point SF from 
DMGetPointSF() to somehow broadcast the support sizes? Or is there a 
simpler way?


- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

[petsc-users] make check vs make test

[Adrian Croucher]

hi

In the latest PETSc version (3.13), 'make test' seems to be doing a much 
more comprehensive set of tests than it used to.


That is ok except that I have a CI pipeline which builds PETSc and my 
code, and then runs my tests. It also does a 'make test' after building 
PETSc and that is now making the CI time out.


At the end of the PETSc build process it says I should run 'make check' 
to "check if the libraries are working". I can't seem to find any 
documentation on what 'make check' actually does, but it looks like it 
runs a much reduced set of tests. Currently I am running both 'make 
check' and 'make test'.


For this use case should it be sufficient to run 'make check' and omit 
'make test'?


- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] Segfault in DMPlexDistributeFieldIS with -log_view

[Adrian Croucher via petsc-users]

It's OK, I found the problem.

I was accidentally passing in a DM which hadn't been created yet. This 
didn't matter when not using -log_view, because it looks like the DM 
isn't actually used for anything else inside DMPlexDistributeFieldIS(). 
But when you run with -log_view it tries to get its communicator.


- Adrian

On 13/09/19 2:31 PM, Adrian Croucher wrote:

hi

My code is using DMPlexDistributeFieldIS() to distribute an index set, 
and it seems to work ok, except if I run with -log_view.


In that case I get the error below.

The code (Fortran) looks like this:

    call PetscSectionCreate(PETSC_COMM_WORLD, dist_section, ierr)
    CHKERRQ(ierr)
    call ISCreate(PETSC_COMM_WORLD, dist_index_set, ierr)
    CHKERRQ(ierr)
    call DMPlexDistributeFieldIS(self%dm, sf, section, &
 index_set, dist_section, &
 dist_index_set, ierr); CHKERRQ(ierr)
    call PetscSectionDestroy(dist_section, ierr); CHKERRQ(ierr)
    call ISDestroy(index_set, ierr); CHKERRQ(ierr)
    index_set = dist_index_set



--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] Segfault in DMPlexDistributeFieldIS with -log_view

[Adrian Croucher via petsc-users]

PS if I run with -start_in_debugger and do a backtrace I get the following:

Thread 1 "waiwera" received signal SIGSEGV, Segmentation fault.
0x7fe12231b761 in PetscObjectComm (obj=0xfffe)
    at /home/acro018/software/PETSc/code/src/sys/objects/gcomm.c:33
33    return obj->comm;
(gdb) bt
#0  0x7fe12231b761 in PetscObjectComm (obj=0xfffe)
    at /home/acro018/software/PETSc/code/src/sys/objects/gcomm.c:33
#1  0x7fe12231488c in PetscLogEventBeginDefault (event=133, t=0,
    o1=0xfffe, o2=0x0, o3=0x0, o4=0x0)
    at 
/home/acro018/software/PETSc/code/src/sys/logging/utils/eventlog.c:647

#2  0x7fe1231bcbbd in DMPlexDistributeFieldIS (dm=0xfffe,
    pointSF=0x55857b89eb40, originalSection=0x55857b7bdaf0,
    originalIS=0x55857b793270, newSection=0x55857b8a26d0, 
newIS=0x7ffc52486d58)
    at 
/home/acro018/software/PETSc/code/src/dm/impls/plex/plexdistribute.c:832

#3  0x7fe1233e5c9a in dmplexdistributefieldis_ (
    dm=0x55857a4cdea8 ,
    pointSF=0x55857a4cdf20 , 
originalSection=0x7ffc52486d88,

    originalIS=0x55857a4cdf10 , newSection=0x7ffc52486d50,
    newIS=0x7ffc52486d58, __ierr=0x7ffc52486d4c)
    at 
/home/acro018/software/PETSc/code/src/dm/impls/plex/ftn-auto/plexdistributef.c:135

#4  0x55857a202cc8 in mesh_module::mesh_distribute_index_set (self=...,
    sf=..., section=..., index_set=...) at ../src/mesh.F90:3393
#5  0x55857a22b60c in mesh_module::mesh_distribute (self=...)
    at ../src/mesh.F90:152
#6  0x55857a224626 in mesh_module::mesh_configure (self=..., 
gravity=...,

    json=0x55857b574180, logfile=..., err=0) at ../src/mesh.F90:826
#7  0x55857a1baf14 in flow_simulation_module::flow_simulation_init (
---Type  to continue, or q  to quit---Quit

On 13/09/19 2:31 PM, Adrian Croucher wrote:

hi

My code is using DMPlexDistributeFieldIS() to distribute an index set, 
and it seems to work ok, except if I run with -log_view.


In that case I get the error below.

The code (Fortran) looks like this:

    call PetscSectionCreate(PETSC_COMM_WORLD, dist_section, ierr)
    CHKERRQ(ierr)
    call ISCreate(PETSC_COMM_WORLD, dist_index_set, ierr)
    CHKERRQ(ierr)
    call DMPlexDistributeFieldIS(self%dm, sf, section, &
 index_set, dist_section, &
 dist_index_set, ierr); CHKERRQ(ierr)
    call PetscSectionDestroy(dist_section, ierr); CHKERRQ(ierr)
    call ISDestroy(index_set, ierr); CHKERRQ(ierr)
    index_set = dist_index_set

I'm running the master branch.

Any clues?

- Adrian

--

[0]PETSC ERROR: 

[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, 
probably memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or 
-on_error_attach_debugger
[0]PETSC ERROR: or see 
https://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind
[0]PETSC ERROR: or try http://valgrind.org on GNU/linux and Apple Mac 
OS X to find memory corruption errors

[0]PETSC ERROR: likely location of problem given in stack below
[0]PETSC ERROR: -  Stack Frames 

[1]PETSC ERROR: 

[1]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, 
probably memory access out of range
[1]PETSC ERROR: Try option -start_in_debugger or 
-on_error_attach_debugger
[1]PETSC ERROR: or see 
https://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind
[1]PETSC ERROR: or try http://valgrind.org on GNU/linux and Apple Mac 
OS X to find memory corruption errors

[1]PETSC ERROR: likely location of problem given in stack below
[1]PETSC ERROR: -  Stack Frames 

[1]PETSC ERROR: Note: The EXACT line numbers in the stack are not 
available,
[1]PETSC ERROR:   INSTEAD the line number of the start of the 
function

[1]PETSC ERROR:   is given.
[1]PETSC ERROR: [1] PetscLogEventBeginDefault line 642 
/home/acro018/software/PETSc/code/src/sys/logging/utils/eventlog.c
[1]PETSC ERROR: [1] DMPlexDistributeFieldIS line 831 
/home/acro018/software/PETSc/code/src/dm/impls/plex/plexdistribute.c
[0]PETSC ERROR: Note: The EXACT line numbers in the stack are not 
available,
[0]PETSC ERROR:   INSTEAD the line number of the start of the 
function

[0]PETSC ERROR:   is given.
[0]PETSC ERROR: [0] PetscLogEventBeginDefault line 642 
/home/acro018/software/PETSc/code/src/sys/logging/utils/eventlog.c
[0]PETSC ERROR: [1]PETSC ERROR: - Error Message 
--

[1]PETSC ERROR: Signal received
[1]PETSC ERROR: See 
https://www.mcs.anl.gov/petsc/documentation/faq.html for trouble 
shooting.
[1]PETSC ERROR: Petsc Development GIT revision: 
v3.11.3-1582-gcb66735359  GIT Date: 2019-08-04 16:01:27 -0500
[1]PETSC ERROR: waiwera on a linux-gnu-c-debug named en-354401 by 
acro01

[petsc-users] Segfault in DMPlexDistributeFieldIS with -log_view

[Adrian Croucher via petsc-users]

hi

My code is using DMPlexDistributeFieldIS() to distribute an index set, 
and it seems to work ok, except if I run with -log_view.


In that case I get the error below.

The code (Fortran) looks like this:

    call PetscSectionCreate(PETSC_COMM_WORLD, dist_section, ierr)
    CHKERRQ(ierr)
    call ISCreate(PETSC_COMM_WORLD, dist_index_set, ierr)
    CHKERRQ(ierr)
    call DMPlexDistributeFieldIS(self%dm, sf, section, &
 index_set, dist_section, &
 dist_index_set, ierr); CHKERRQ(ierr)
    call PetscSectionDestroy(dist_section, ierr); CHKERRQ(ierr)
    call ISDestroy(index_set, ierr); CHKERRQ(ierr)
    index_set = dist_index_set

I'm running the master branch.

Any clues?

- Adrian

--

[0]PETSC ERROR: 

[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, 
probably memory access out of range

[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[0]PETSC ERROR: or see 
https://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind
[0]PETSC ERROR: or try http://valgrind.org on GNU/linux and Apple Mac OS 
X to find memory corruption errors

[0]PETSC ERROR: likely location of problem given in stack below
[0]PETSC ERROR: -  Stack Frames 

[1]PETSC ERROR: 

[1]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, 
probably memory access out of range

[1]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[1]PETSC ERROR: or see 
https://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind
[1]PETSC ERROR: or try http://valgrind.org on GNU/linux and Apple Mac OS 
X to find memory corruption errors

[1]PETSC ERROR: likely location of problem given in stack below
[1]PETSC ERROR: -  Stack Frames 


[1]PETSC ERROR: Note: The EXACT line numbers in the stack are not available,
[1]PETSC ERROR:   INSTEAD the line number of the start of the function
[1]PETSC ERROR:   is given.
[1]PETSC ERROR: [1] PetscLogEventBeginDefault line 642 
/home/acro018/software/PETSc/code/src/sys/logging/utils/eventlog.c
[1]PETSC ERROR: [1] DMPlexDistributeFieldIS line 831 
/home/acro018/software/PETSc/code/src/dm/impls/plex/plexdistribute.c

[0]PETSC ERROR: Note: The EXACT line numbers in the stack are not available,
[0]PETSC ERROR:   INSTEAD the line number of the start of the function
[0]PETSC ERROR:   is given.
[0]PETSC ERROR: [0] PetscLogEventBeginDefault line 642 
/home/acro018/software/PETSc/code/src/sys/logging/utils/eventlog.c
[0]PETSC ERROR: [1]PETSC ERROR: - Error Message 
--

[1]PETSC ERROR: Signal received
[1]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html 
for trouble shooting.
[1]PETSC ERROR: Petsc Development GIT revision: 
v3.11.3-1582-gcb66735359  GIT Date: 2019-08-04 16:01:27 -0500
[1]PETSC ERROR: waiwera on a linux-gnu-c-debug named en-354401 by 
acro018 Fri Sep 13 14:20:41 2019
[1]PETSC ERROR: Configure options --with-x --download-hdf5 
--download-zlib --download-netcdf --download-pnetcdf --download-exodusii 
--download-triangle --download-ptscotch --download-chaco --download-hypre

[1]PETSC ERROR: #1 User provided function() line 0 in  unknown file
[0] DMPlexDistributeFieldIS line 831 
/home/acro018/software/PETSc/code/src/dm/impls/plex/plexdistribute.c
[0]PETSC ERROR: - Error Message 
--

[0]PETSC ERROR: Signal received
[0]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html 
for trouble shooting.
[0]PETSC ERROR: Petsc Development GIT revision: 
v3.11.3-1582-gcb66735359  GIT Date: 2019-08-04 16:01:27 -0500
[0]PETSC ERROR: waiwera on a linux-gnu-c-debug named en-354401 by 
acro018 Fri Sep 13 14:20:41 2019
[0]PETSC ERROR: Configure options --with-x --download-hdf5 
--download-zlib --download-netcdf --download-pnetcdf --download-exodusii 
--download-triangle --download-ptscotch --download-chaco --download-hypre

[0]PETSC ERROR: #1 User provided function() line 0 in  unknown file

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

[petsc-users] DMLocalToLocal for DMPlex

[Adrian Croucher via petsc-users]

hi

If I have a local vector and I want to update all its ghost values from 
its non-ghost values, should I use DMLocalToLocalBegin()/End() ?


I have tried it and it gives me an error: "This DM does not support 
local to local maps".


The DM is a DMPlex. Is the local-to-local operation not implemented for 
DMPlex?


Or should I be using something else to do this?

- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] Copying to Vec with boundary cells

[Adrian Croucher via petsc-users]

Actually it looks like a simple VecCopy() works in this situation.

It checks the local sizes of the two vectors, and decides they are still 
the same, so is happy to do the copy.


- Adrian

On 18/07/19 10:31 AM, Adrian Croucher wrote:

hi

I have a Vec created on a distributed DMPlex (using 
DMCreateGlobalVector() or DMCreateLocalVector()).


I then add FV boundary condition cells to the DMPlex, using 
DMPlexConstructGhostCells().


I can then create a new Vec on the new DM. What is the best way to 
copy the values from the old Vec to the new one?


Regards, Adrian


--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

[petsc-users] Copying to Vec with boundary cells

[Adrian Croucher via petsc-users]

hi

I have a Vec created on a distributed DMPlex (using 
DMCreateGlobalVector() or DMCreateLocalVector()).


I then add FV boundary condition cells to the DMPlex, using 
DMPlexConstructGhostCells().


I can then create a new Vec on the new DM. What is the best way to copy 
the values from the old Vec to the new one?


Regards, Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] DMPlexDistributeField

[Adrian Croucher via petsc-users]

Probably don't want the two extra PetscSFView() calls (lines 1687, 1688) 
though- presumably they were just for temporary debugging?


- Adrian

On 11/07/19 2:18 PM, Adrian Croucher wrote:



On 10/07/19 3:40 PM, Matthew Knepley wrote:


Sorry this took so long. I put in a PR for this:

https://bitbucket.org/petsc/petsc/pull-requests/1858/knepley-fix-plex-distribute-overlap/diff

I think it fixes your problem.


Thanks very much Matt, that does fix it. Good stuff!

- Adrian


--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] DMPlexDistributeField

[Adrian Croucher via petsc-users]

On 10/07/19 3:40 PM, Matthew Knepley wrote:


Sorry this took so long. I put in a PR for this:

https://bitbucket.org/petsc/petsc/pull-requests/1858/knepley-fix-plex-distribute-overlap/diff

I think it fixes your problem.


Thanks very much Matt, that does fix it. Good stuff!

- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] DMPlexDistributeField

[Adrian Croucher via petsc-users]

hi Matt,

On 29/06/19 12:07 AM, Matthew Knepley wrote:


Okay, I ran through this. I am attaching my C version which was easier 
for me to play with in the debugger, but it should not be
hard to put the few lines into yours. The problem comes from using 
overlap in the redistribution. You can see this by running


  mpiexec -n 2 ./distgeom2 -overlap 1 -overlap2 0

which works fine. If overlap2 is 1, then the face data is shifted.



Yes, that's the same thing I found.


Something in DistributeField() does not know that we changed
the topology after we distributed by adding overlap cells. This is 
strange since DMPlexMigrate() works on the coordinates for
the overlapped mesh, however they always start at 0 in the vector, so 
it must be that some offset is missing. I am looking for it now.



It looked to me like the problem might lie not in 
DMPlexDistributeField() but in DMPlexDistribute(), because if you view 
the redistribution SF coming out of the second call to 
DMPlexDistribute() it looks wrong.


DMPlexDistributeField() appeared to be doing what it was told, but the 
SF it was given didn't look right to me.


- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] DMPlexDistributeField

[Adrian Croucher via petsc-users]

On 28/06/19 10:09 AM, Zhang, Junchao wrote:


Check how the graph is created and then whether the parameters to 
PetscSFSetGraph() are correct.



Yes, unfortunately I don't have a good enough understanding of how 
DMPlexDistribute() works to see what the problem is.


- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] DMPlexDistributeField

[Adrian Croucher via petsc-users]

hi

On 28/06/19 3:14 AM, Zhang, Junchao wrote:

You can dump relevant SFs to make sure their graph is correct.



Yes, I'm doing that, and the graphs don't look correct.

- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] DMPlexDistributeField

[Adrian Croucher via petsc-users]

hi

On 25/06/19 11:22 AM, Adrian Croucher wrote:


I know you originally recommended using overlap = 0 for the initial 
distribution and only adding overlap for the redistribution. But then 
Stefano indicated that it should work with overlap now. And it would 
simplify my code if I could use overlap for the initial distribution 
(because if dual porosity cells are not being used, then there is no 
second redistribution).


Actually I have also tried this test with overlap = 0 for the initial 
distribution, and overlap = 1 for the redistribution, and that doesn't 
seem to give correct results for the face geometry either. So that would 
suggest that, if there is a problem, it probably doesn't lie in the 
stuff that removes overlap cells from the original DM.


I am running PETSc master 3.11.2. I see 3.11.3 is just out, not sure if 
there are any changes there which might affect this?


- Adrian


Is this a bug or is there something I'm doing wrong?

- Adrian


--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] DMPlexDistributeField

[Adrian Croucher via petsc-users]

hi

Thanks Matt for the explanation about this.

I have been trying a test which does the following:

1) read in DMPlex from file

2) distribute it, with overlap = 1, using DMPlexDistribute()

3) create FVM cell and face geometry vectors using 
DMPlexComputeGeometryFVM()


4) re-distribute, again with overlap = 1, using DMPlexDistribute()

5) distribute the cell and face geometry vectors using 
DMPlexDistributeField()



Steps 4) and 5) should do essentially nothing, because the mesh has 
already been distributed (but in my actual non-test code, there is 
additional stuff between steps 3) and 4) where dual porosity cells are 
added to the DM).


So I expect the cell and face geometry vectors to be essentially 
unchanged from the redistribution. And the redistribution SF (from the 
second distribution) should be just an identity mapping on the cell and 
face points (except for the overlap ghost points).


This is true for the cells, but not the faces. I've attached the example 
code and mesh. It is a simple mesh with 10 cells in a horizontal line, 
each cell 50x50x50 m.


If I run on 2 processes, there are 5 cells (points 0 - 4) on each rank, 
with centroids at 25, 75, 125, 175 and 225 m on rank 0, and 275, 325, 
375, 425 and 475 m on rank 1. The internal faces are the points 36, 42, 
47 and 52 on rank 0, and 34, 37, 42, 47 and 52 on rank 1. On rank 0 
these should have centroids at 50, 100, 150 and 200 m respectively; on 
rank 1 they should be at 250, 300, 350 and 400 m. This is true before 
redistribution.


After redistribution, the cells centroids are still correct, and the 
face data on rank 1 are OK, but the face data on rank 0 are all wrong.


If you look at the redistribution SF the entries for the rank 0 face 
data are 36 <- (0,40), 42 <- (0,46), 47 <- (0,51), 52 <- (0,56), instead 
of the expected 36 <- (0,36), 42 <- (0,42), 47 <- (0,47), 52 <- (0,52). 
The SF for the rank 1 faces is OK.


If you change the overlap from 1 to 0, it works as expected. So it looks 
to me like something isn't quite right with the SF for faces when there 
is overlap. On rank 0 all the entries seem to be shifted up by 4.


I know you originally recommended using overlap = 0 for the initial 
distribution and only adding overlap for the redistribution. But then 
Stefano indicated that it should work with overlap now. And it would 
simplify my code if I could use overlap for the initial distribution 
(because if dual porosity cells are not being used, then there is no 
second redistribution).


Is this a bug or is there something I'm doing wrong?

- Adrian

On 23/06/19 4:39 PM, Matthew Knepley wrote:
On Fri, Jun 21, 2019 at 12:49 AM Adrian Croucher via petsc-users 
mailto:petsc-users@mcs.anl.gov>> wrote:


I have been trying to get this FVM geometry data re-distribution
to work
using DMPlexDistributeField().

It seems to be working OK for the cell geometry data (cell volumes
and
centroids). But it is making a mess of the face geometry data (face
normals and centroids).

Should I even expect DMPlexDistributeField() to work for
redistributing
a vector of data defined on mesh faces? Or is there some reason I
haven't thought of, which means that will never work?


Sorry this took a long time. The place I was in France did not have 
internet.

Here is how this stuff works:

  1) You start with a Section and local vector. The Section describes 
layout of data in
       the local vector by mapping mesh points to {# dof, offset}. For 
a small example,
       suppose I had two triangles sharing an edge on a sequential 
mesh for 2 procs.

       The mesh points would be

         [0, 1]:  Cells
         [2, 5]:  Vertices, where 3,4 are shared
         [6, 10]: Edges, where 8 is shared

       A Section for face normals would then look like

        Process 0
        [0, 5]: {0, 0}   Meaning no variables lie on cells or vertices
        6:       {2, 0}   One vector per face
        7:       {2, 2}
        8:       {2, 4}
        9:       {2, 6}
        10      {2, 8}
        Process 1
        empty

     The vector would just have the face normal values in the 
canonical order. You can use

     PetscSectionView() to check that yours looks similar.

  2) Now we add a PetscSF describing the redistribution. An SF is a 
map from a given set of
      integers (leaves) to pairs (int, rank) called roots. Many leaves 
can point to one root. To begin,

      we provide an SF mapping mesh points to the new distribution

        Process 0
        0 -> {0, 0}
        1 -> {2, 0}
        2 -> {3, 0}
        3 -> {4, 0}
        4 -> {6, 0}
        5 -> {7, 0}
        6 -> {8, 0}
        Process 1
        0 -> {1, 0}
        1 -> {4, 0}
        2 -> {3, 0}
        3 -> {5, 0}
        4 -> {8, 0}
        5 -> {9, 0}
        6 -> {10, 0}

     This tells Petsc how to move the points to the parallel 
distribution. You can use PetscS

Re: [petsc-users] DMPlexDistributeField

[Adrian Croucher via petsc-users]

I have been trying to get this FVM geometry data re-distribution to work 
using DMPlexDistributeField().


It seems to be working OK for the cell geometry data (cell volumes and 
centroids). But it is making a mess of the face geometry data (face 
normals and centroids).


Should I even expect DMPlexDistributeField() to work for redistributing 
a vector of data defined on mesh faces? Or is there some reason I 
haven't thought of, which means that will never work?


The only example I could find anywhere of DMPlexDistributeField() being 
used is in DMPlexDistributeCoordinates() so I've been basing what I'm 
doing on that.


(I think I have answered my own questions below by experiment- 1) local 
vectors should work (they do in DMPlexDistributeCoordinates); 2) 
probably doesn't matter; 3) yes.)


- Adrian

On 6/06/19 1:42 PM, Adrian Croucher wrote:

hi

I have some questions about using the DMPlexDistributeField() 
function. I have finite volume mesh geometry data stored in two local 
vectors created using DMPlexComputeGeometryFVM(), and I need to 
redistribute these after calling DMPlexDistribute() to redistribute my 
mesh. (I need the geometry data before redistribution, so I can't just 
wait until after redistribution to create them.)


So I figured DMPlexDistributeField() looks like the thing to use for 
that, using the SF that comes out of DMPlexDistribute().


1) Does DMPlexDistributeField() work on local vectors or do they have 
to be global ones?


2) It takes a 'dm' parameter and the documentation says this is "The 
DMPlex object", but is that the original DM (before redistribution) or 
the redistributed one, or does it not matter? It looks like it only 
uses the DM to get the vector type.


3) It looks like you need to manually create the newSection and newVec 
output parameters before passing them in to this routine, is that 
correct?


- Adrian


--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] parallel dual porosity

[Adrian Croucher via petsc-users]

On 30/05/19 5:26 PM, Stefano Zampini wrote:

Matt,

 redistribution with overlapped mesh is fixed in master (probably also 
in maint)



Even better. Thanks very much...

- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] parallel dual porosity

[Adrian Croucher via petsc-users]

On 30/05/19 2:45 PM, Matthew Knepley wrote:


Hmm, I had not thought about that. It will not do that at all. We have 
never rebalanced a simulation
using overlap cells. I would have to write the code that strips them 
out. Not hard, but more code.
If you only plan on redistributing once, you can wait until then to 
add the overlap cells.



So for now at least, you would suggest doing the initial distribution 
with overlap = 0, and then the redistribution with overlap = 1?


(I shouldn't need to redistribute more than once.)

- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] parallel dual porosity

[Adrian Croucher via petsc-users]

hi

On 28/05/19 11:32 AM, Matthew Knepley wrote:


I would not do that. It should be much easier, and better from a 
workflow standpoint,
to just redistribute in parallel. We now have several test examples 
that redistribute

in parallel, for example

https://bitbucket.org/petsc/petsc/src/cd762eb66180d8d1fcc3950bd19a3c1b423f4f20/src/dm/impls/plex/examples/tests/ex1.c#lines-486

Let us know if you have problems.


If you use DMPlexDistribute() a second time to redistribute the mesh, 
presumably it will first delete all the partition ghost cells created 
from the initial distribution (and create new ones after redistribution).


If I have also used DMPlexConstructGhostCells() to create boundary 
condition ghost cells on the distributed mesh, what happens to them when 
I redistribute? Do they get copied over to the redistributed mesh? or is 
it better not to add them until the mesh has been redistributed?


- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] parallel dual porosity

[Adrian Croucher via petsc-users]

Oh, so DMPlexDistribute() can redistribute an already-distributed mesh? 
I didn't know that- had assumed it was only for distributing a serial mesh.


If so, that sounds like a much better idea. I'll check it out.

Thanks!

- Adrian

On 28/05/19 11:32 AM, Matthew Knepley wrote:


I would not do that. It should be much easier, and better from a 
workflow standpoint,
to just redistribute in parallel. We now have several test examples 
that redistribute

in parallel, for example

https://bitbucket.org/petsc/petsc/src/cd762eb66180d8d1fcc3950bd19a3c1b423f4f20/src/dm/impls/plex/examples/tests/ex1.c#lines-486

Let us know if you have problems.

  Thanks,

     Matt


--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

[petsc-users] parallel dual porosity

[Adrian Croucher via petsc-users]

hi

A couple of years back I was asking questions here about implementing 
"dual porosity" finite volume methods via PETSc (in which flow in 
fractured media is represented by adding extra "matrix" cells nested 
inside the original mesh cells).


At the time I was asking about how to solve the resulting linear 
equations more efficiently (I still haven't worked on that part of it 
yet, so at present it's still just using a naive linear solve which 
doesn't take advantage of the particular sparsity pattern), and about 
how to add the extra cells into the DMPlex mesh, which I figured out how 
to do.


It is working OK except that strong scaling performance is not very 
good, if dual porosity is applied over only part of the mesh. I think 
the reason is that I read the mesh in and distribute it, then add the 
dual porosity cells in parallel on each process. So some processes can 
end up with more cells than others, in which case the load balancing is bad.


I'm considering trying to change it so that I add the dual porosity 
cells to the DMPlex in serial, before distribution, to regain decent 
load balancing.


To do that, I'd also need to compute the cell centroids in serial (as 
they are often used to identify which cells should have dual porosity 
applied), using DMPlexComputeGeometryFVM(). The geometry vectors would 
then have to be distributed later, I guess using something like 
DMPlexDistributeField().


Should I expect a significant performance hit from calling 
DMPlexComputeGeometryFVM() on the serial mesh compared with doing it (as 
now) on the distributed mesh? It will increase the serial fraction of 
the code but as it's only done once at the start I'm hoping the benefits 
will outweigh the costs.


- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email:a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] pkg-config

[Adrian Croucher]

OK thanks Jed, that's fixed it.

- Adrian

On 8/10/18 2:14 PM, Jed Brown wrote:

This was a bug, now fixed in 'maint' and will be in the next patch
release.  You can repair the PETSc.pc, update to 'maint', or apply this
2-line patch.

https://bitbucket.org/petsc/petsc/commits/e1e675de8c3f48a6c27b15412e00ded76072f735

Adrian Croucher  writes:


hi

I am experimenting with using pkg-config to help build my code, and
using the PETSc.pc file in $PETSC_DIR/$PETSC_ARCH/lib/pkgconfig.

At present it is not working, because the PETSc.pc file sets Libs:
-L${libdir} -lpetsc

where libdir=${prefix}/lib, and (in my case at least) prefix is equal to
$PETSC_DIR.

But the PETSc library is not in $PETSC_DIR/lib, it is in
$PETSC_DIR/$PETSC_ARCH/lib.

Have I messed something up in my PETSc config so that the PETSc.pc file
is not quite right?

- Adrian


--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

[petsc-users] pkg-config

[Adrian Croucher]

hi

I am experimenting with using pkg-config to help build my code, and 
using the PETSc.pc file in $PETSC_DIR/$PETSC_ARCH/lib/pkgconfig.


At present it is not working, because the PETSc.pc file sets Libs: 
-L${libdir} -lpetsc


where libdir=${prefix}/lib, and (in my case at least) prefix is equal to 
$PETSC_DIR.


But the PETSc library is not in $PETSC_DIR/lib, it is in 
$PETSC_DIR/$PETSC_ARCH/lib.


Have I messed something up in my PETSc config so that the PETSc.pc file 
is not quite right?


- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] v3.9.3 config problem

[Adrian Croucher]

Whew, that is looking much healthier now. Thanks very much for tracking 
that down...


- Adrian


On 20/08/18 15:24, Smith, Barry F. wrote:

  This is insanely bad for users (including me as a user).

  You shouldn't need to delete anything.

  You just need the additional flag of either --download-zlib or 
--with-zlib (if it is installed on your system).

  The new branch barry/fix-insanity-netcdf/maint will give you a very 
useful error message if you don't pass in the zlib. So people won't waste their 
time on this same problem in the future.

   Sorry to waste your time,

 Barry

   The netcdf code didn't build netcdf4 without access to zlib but it 
didn't bother to tell you it wasn't building it. Then when exodusii needed 
netcdf4 it assumed it was there (though it wasn't) and hence the bad compile. 
I'm also surprised this didn't happen before to other users.



--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] v3.9.3 config problem

[Adrian Croucher]

On 20/08/18 14:01, Adrian Croucher wrote:
Should I try wiping out the externalpackages directory as well as 
$PETSC_ARCH?


I did try this anyway. First it gave the following errors:

Error! missing file: 
/home/acro018/software/PETSc/code/externalpackages-bak/sowing-1.1.16d/src/sys/makefile
Error! missing file: 
/home/acro018/software/PETSc/code/externalpackages-bak/sowing-1.1.16d/src/doctext/test/makefile
Error! missing file: 
/home/acro018/software/PETSc/code/externalpackages-bak/sowing-1.1.16d/src/bfort/testing/makefile


Not sure if these matter. It continued but then got the same errors 
trying to build Exodus.


- Adrian


- Adrian

On 20/08/18 14:00, Smith, Barry F. wrote:

    Trying to reproduce it now but I'm pretty sure we've tested it.

    Do you have some old files lying around that might be causing the 
problem?


/home/acro018/software/PETSc/code/linux-gnu-c-opt/externalpackages/git.exodusii/packages/seacas/libraries/exodus/src/ex_create_par.c:375:9: 
error: 'netcdf4_mode' undeclared (first use in this function)

  if (netcdf4_mode == -1) {

Maybe externalpackages has downloads from a previous build still 
located that don't work with the latest exodusii?



    Barry




On Aug 19, 2018, at 8:51 PM, Adrian Croucher 
 wrote:


hi

I just pulled the latest PETSc master branch (v3.9.3) but am having 
trouble building it. It's been running fine on v3.8.4.


I deleted my $PETSC_ARCH directory and tried to configure. But it 
seems to be falling over when trying to build Exodus. I've attached 
the configure.log.


Any ideas?

- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611






--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] v3.9.3 config problem

[Adrian Croucher]

Should I try wiping out the externalpackages directory as well as 
$PETSC_ARCH?


- Adrian

On 20/08/18 14:00, Smith, Barry F. wrote:

Trying to reproduce it now but I'm pretty sure we've tested it.

Do you have some old files lying around that might be causing the problem?

/home/acro018/software/PETSc/code/linux-gnu-c-opt/externalpackages/git.exodusii/packages/seacas/libraries/exodus/src/ex_create_par.c:375:9:
 error: 'netcdf4_mode' undeclared (first use in this function)
  if (netcdf4_mode == -1) {

Maybe externalpackages has downloads from a previous build still located that 
don't work with the latest exodusii?


Barry





On Aug 19, 2018, at 8:51 PM, Adrian Croucher  wrote:

hi

I just pulled the latest PETSc master branch (v3.9.3) but am having trouble 
building it. It's been running fine on v3.8.4.

I deleted my $PETSC_ARCH directory and tried to configure. But it seems to be 
falling over when trying to build Exodus. I've attached the configure.log.

Any ideas?

- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611




--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] Newton methods that converge all the time

[Adrian Croucher]

hi Henrik,


On 01/12/17 22:24, Buesing, Henrik wrote:

The beauty of the pressure-enthalpy formulation is that you do not 
need to initialize the saturations with small epsilon when you go into 
two-phase. You can directly compute a correct enthalpy (and from it a 
saturation) and "jump" into the two-phase region.

Anyhow, no one says that you are not jumping in and out of the two-phase 
region. Consequently, there are examples for which either method works best.

I am trying to simulate a supercritical reservoir (T>450 °C, p>35 MPa) from 
surface down to 3.5 km depth. (It is in Italy so geothermal gradient is large).

There is a steam region, which forms and either I get osciallations in enthalpy 
(saturations) or small time-steps kill me. I want to simulate a quasi 
steady-state after 1 million years.


We regularly simulate systems like this, with steam zones near the top, 
and sometimes supercritical fluid at the bottom. We generally run steady 
states up to time step sizes of 1E15 s or more. We have done quite a bit 
of work on getting the variable-switching approach to work well under 
these conditions.


- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] Newton methods that converge all the time

[Adrian Croucher]

   Please describe in some detail how you are handling phase change. If
you have if () tests of any sort in your FormFunction() or
FormJacobian() this can kill Newton's method. If you are using "variable
switching" this WILL kill Newtons' method. Are you monkeying with phase
definitions in TSPostStep or with SNESLineSearchSetPostCheck(). This
will also kill Newton's method.


I'm doing variable switching (in a geothermal flow application) with 
Newton's method (in SNESLineSearchSetPostCheck()) and it generally works 
fine.


For pure water (no other components present) my variables are pressure 
and temperature for single-phase (liquid or vapour) and pressure and 
vapour saturation for two-phase.


You have to be pretty careful how you do the switching though.

- Adrian
--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] ISGlobalToLocalMappingApplyBlock

[Adrian Croucher]

I've debugged into the ISGlobalToLocalMappingApplyBlock() function and 
it seems to me the bounds checking in there is not correct when the 
blocksize is > 1.


It checks against the same bounds, scaled up by the blocksize, in both 
the block and non-block versions of the function. I think for the block 
version the bounds should not be scaled.


I've just created a pull request 
(acroucher/fix-IS-global-to-local-mapping-block) with a suggested fix.


- Adrian


On 16/11/17 11:52, Adrian Croucher wrote:
I actually attached the wrong test program last time- I've attached 
the right one here, which is much simpler. It test global indices 0, 
1, ... 9.


If I run on 2 processes, the local indices it returns are:

rank 0: 0, 1, 2, 3, 4, 0, 0, 0, -253701943, 0
rank 1: -1, -1, -1, -1, -1, -1, -1, -1, -1, -1

The results I expected are:

rank 0: 0, 1, 2, 3, 4, -1, -1, -1, -1, -1
rank 1: -1, -1, -1, -1, -1, 0, 1, 2, 3, 4

So the results for global indices 0, 1,... 4 are what I expected, on 
both ranks. But the results for global indices 5, 6, ... 9 are not.


I tried increasing the blocksize to 3 or 4, and the results were 
exactly the same.


It only gives the results I expected if I change the blocksize to 1.

- Adrian



--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] ISGlobalToLocalMappingApplyBlock

[Adrian Croucher]

I actually attached the wrong test program last time- I've attached the 
right one here, which is much simpler. It test global indices 0, 1, ... 9.


If I run on 2 processes, the local indices it returns are:

rank 0: 0, 1, 2, 3, 4, 0, 0, 0, -253701943, 0
rank 1: -1, -1, -1, -1, -1, -1, -1, -1, -1, -1

The results I expected are:

rank 0: 0, 1, 2, 3, 4, -1, -1, -1, -1, -1
rank 1: -1, -1, -1, -1, -1, 0, 1, 2, 3, 4

So the results for global indices 0, 1,... 4 are what I expected, on 
both ranks. But the results for global indices 5, 6, ... 9 are not.


I tried increasing the blocksize to 3 or 4, and the results were exactly 
the same.


It only gives the results I expected if I change the blocksize to 1.

- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

program testl2g

  ! Test PETSc local-to-global mapping

#include 

  use petsc

  implicit none
  PetscInt, parameter :: blocksize = 1
  PetscInt, parameter :: n = 5, nin = 10
  PetscMPIInt :: rank
  PetscInt :: i
  PetscInt, allocatable :: indices(:)
  ISLocalToGlobalMapping :: l2g
  PetscInt :: global(nin), local(nin), nout
  PetscErrorCode :: ierr

  call PetscInitialize(PETSC_NULL_CHARACTER, ierr); CHKERRQ(ierr)
  call MPI_comm_rank(PETSC_COMM_WORLD, rank, ierr)

  select case (rank)
  case (0)
 indices = [(i, i = 0, 4)]
  case (1)
 indices = [(i, i = 5, 9)]
  end select

  call ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, blocksize, n, indices, &
   PETSC_COPY_VALUES, l2g, ierr); CHKERRQ(ierr)
  call ISLocalToGlobalMappingView(l2g, PETSC_VIEWER_STDOUT_WORLD, ierr); CHKERRQ(ierr)

  global = [(i, i = 0, 9)]
  call ISGlobalToLocalMappingApplyBlock(l2g, IS_GTOLM_MASK, nin, &
   global, nout, local, ierr); CHKERRQ(ierr)
  write(*,*) 'rank:', rank, 'global:', global, 'local:', local, 'nout:', nout
  
  call ISLocalToGlobalMappingDestroy(l2g, ierr); CHKERRQ(ierr)
  call PetscFinalize(ierr); CHKERRQ(ierr)

end program testl2g

Re: [petsc-users] ISGlobalToLocalMappingApplyBlock

[Adrian Croucher]

hi Dave,


On 15/11/17 21:34, Dave May wrote:



Or am I wrong to expect this to give the same results regardless of
blocksize?



Yep.


Maybe I am not using this function correctly then.

The man page says it "Provides the local block numbering for a list of 
integers specified with a block global numbering."


So I thought if I put in global block indices it would give me the 
corresponding local block indices- which would be the same regardless of 
the size of each block.


However the large negative number being printed looks an uninitialized 
variable. This seems odd as with mode = MASK nout should equal N and 
any requested block indices not in the IS should result in -1 being 
inserted in your local_indices array.


What's the value of nout?


nout returns 1 on both ranks, as expected.

- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

[petsc-users] ISGlobalToLocalMappingApplyBlock

[Adrian Croucher]

hi

I'm trying to use ISGlobalToLocalMappingApplyBlock() and am a bit 
puzzled about the results it's giving.


I've attached a small test to illustrate. It just sets up a 
local-to-global mapping with 10 elements. Running on two processes the 
first has global indices 0 - 4 and the the second has 5 - 9. I then try 
to find the local index corresponding to global index 8.


If I set the blocksize parameter to 1, it correctly gives the results -1 
on rank 0 and 3 on rank 1.


But if I set the blocksize to 2 (or more), the results are -253701943 on 
rank 0 and -1 on rank 1. Neither of these are what I expected- I thought 
they should be the same as in the blocksize 1 case.


I'm presuming the global indices I pass in to 
ISGlobalToLocalMappingApplyBlock() should be global block indices (i.e. 
not scaled up by blocksize). If I do scale them up it doesn't give the 
answers I expect either.


Or am I wrong to expect this to give the same results regardless of 
blocksize?


Cheers, Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

program testl2g

  ! Test PETSc local-to-global mapping

#include 

  use petsc

  implicit none

  PetscInt, parameter :: dof = 1, overlap = 1, n = 30
  PetscInt :: dim, i, f, start_face, end_face, num_cells
  PetscReal, parameter :: eps = 1.e-3
  PetscReal :: area
  PetscReal, pointer :: centroid(:), normal(:)
  PetscMPIInt :: rank
  DM :: dm, dist_dm, ghost_dm
  PetscErrorCode :: ierr
  character(40) :: filename = "col30.msh"
  PetscInt, parameter :: bdy_face_index = 94
  character(len = 16) :: boundary_label_name = "boundary"
  PetscDS :: ds
  PetscFV :: fvm
  PetscInt :: global_indices(0: n-1), local_indices(0: n-1), nout
  ISLocalToGlobalMapping :: l2g

  call PetscInitialize(PETSC_NULL_CHARACTER, ierr); CHKERRQ(ierr)
  call MPI_comm_rank(PETSC_COMM_WORLD, rank, ierr)

  call DMPlexCreateFromFile(PETSC_COMM_WORLD, filename, &
   PETSC_TRUE, dm, ierr); CHKERRQ(ierr)

  call DMPlexSetAdjacencyUseCone(dm, PETSC_TRUE, ierr); CHKERRQ(ierr)
  call DMPlexSetAdjacencyUseClosure(dm, PETSC_FALSE, ierr); CHKERRQ(ierr)
  call DMGetDimension(dm, dim, ierr); CHKERRQ(ierr)

  call DMPlexDistribute(dm, overlap, PETSC_NULL_SF, dist_dm, ierr)
  CHKERRQ(ierr)
  if (dist_dm .ne. PETSC_NULL_DM) then
 call DMDestroy(dm, ierr); CHKERRQ(ierr)
 dm = dist_dm
  end if

  call DMCreateLabel(dm, boundary_label_name, ierr); CHKERRQ(ierr)
  call DMPlexGetHeightStratum(dm, 1, start_face, end_face, ierr); CHKERRQ(ierr)
  allocate(centroid(3), normal(3))
  do f = start_face, end_face - 1
 call DMPlexGetSupportSize(dm, f, num_cells, ierr); CHKERRQ(ierr)
 if (num_cells == 1) then
call DMPlexComputeCellGeometryFVM(dm, f, area, centroid, normal, ierr)
if (normal(dim) > eps) then
   call DMSetLabelValue(dm, boundary_label_name, f, 1, ierr)
end if
 end if
  end do
  deallocate(centroid, normal)
  
  call DMPlexConstructGhostCells(dm, boundary_label_name, &
   PETSC_NULL_INTEGER, ghost_dm, ierr); CHKERRQ(ierr)
  if (ghost_dm .ne. PETSC_NULL_DM) then
 call DMDestroy(dm, ierr); CHKERRQ(ierr);
 dm = ghost_dm
  end if

  call DMView(dm, PETSC_VIEWER_STDOUT_WORLD, ierr); CHKERRQ(ierr)

  call PetscFVCreate(PETSC_COMM_WORLD, fvm, ierr); CHKERRQ(ierr)
  call PetscFVSetNumComponents(fvm, dof, ierr); CHKERRQ(ierr)
  call PetscFVSetSpatialDimension(fvm, dim, ierr); CHKERRQ(ierr)
  call DMGetDS(dm, ds, ierr); CHKERRQ(ierr)
  call PetscDSAddDiscretization(ds, fvm, ierr); CHKERRQ(ierr)

  call DMGetLocalToGlobalMapping(dm, l2g, ierr); CHKERRQ(ierr)
  call ISLocalToGlobalMappingView(l2g, PETSC_VIEWER_STDOUT_WORLD, ierr); CHKERRQ(ierr)
  
  global_indices = [(i, i = 0, n-1)]
  call ISGlobalToLocalMappingApplyBlock(l2g, IS_GTOLM_MASK, n, global_indices, &
   nout, local_indices, ierr); CHKERRQ(ierr)
  write(*,*) 'rank:', rank, '(global, local):', &
   [(global_indices(i), local_indices(i), i = 0, n - 1)]

  call DMDestroy(dm, ierr); CHKERRQ(ierr)
  call PetscFinalize(ierr); CHKERRQ(ierr)

end program testl2g

Re: [petsc-users] DMPlex ghost cells

[Adrian Croucher]

Would it also work if I just didn't use DMPlexGetHybridBounds() to 
differentiate between partition ghost cells and boundary ghost cells, 
but instead checked the value of the "ghost" DMLabel?


It appears that DMPlexConstructGhostCells() labels boundary ghosts with 
the ghost label value 1, and partition ghosts with value 2 (in 
DMPlexShiftLabels_Internal()).


If that is reliable, then it might save me having to make sure the 
boundary ghosts always come last.


- Adrian


On 04/08/17 15:11, Matthew Knepley wrote:
On Thu, Aug 3, 2017 at 10:07 PM, Adrian Croucher 
<a.crouc...@auckland.ac.nz <mailto:a.crouc...@auckland.ac.nz>> wrote:


I am adding cells onto a DMPlex (for dual porosity). I notice that
normally the ghost cells in a DMPlex are at the end of the cell list.

Does it matter if the cells I add come after those ghost cells, or
should the ghosts always come last?

It seems to me it probably shoudn't matter if the partition ghosts
(created by DMPlexDistribute()) are not at the end. In my cell
loops I just check the value of the ghost DMLabel and skip them if
necessary.

But perhaps the boundary ghosts (created by
DMPlexConstructGhostCells()) need to come last? Otherwise the
output from DMPlexGetHybridBounds() (to get the bounds on
'interior' cells, i.e. non-boundary-ghosts) probably wouldn't make
much sense?


Yes, the only reason to have them last is to have GetHybridBounds() 
work for that case.


  Thanks,

 Matt


- Adrian

-- 
Dr Adrian Croucher

Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz <mailto:a.crouc...@auckland.ac.nz>
tel: +64 (0)9 923 4611 <tel:%2B64%20%280%299%20923%204611>




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

http://www.caam.rice.edu/~mk51/ <http://www.caam.rice.edu/%7Emk51/>


--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

[petsc-users] DMPlex ghost cells

[Adrian Croucher]

I am adding cells onto a DMPlex (for dual porosity). I notice that 
normally the ghost cells in a DMPlex are at the end of the cell list.


Does it matter if the cells I add come after those ghost cells, or 
should the ghosts always come last?


It seems to me it probably shoudn't matter if the partition ghosts 
(created by DMPlexDistribute()) are not at the end. In my cell loops I 
just check the value of the ghost DMLabel and skip them if necessary.


But perhaps the boundary ghosts (created by DMPlexConstructGhostCells()) 
need to come last? Otherwise the output from DMPlexGetHybridBounds() (to 
get the bounds on 'interior' cells, i.e. non-boundary-ghosts) probably 
wouldn't make much sense?


- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] specifying vertex coordinates using DMPlexCreateFromCellListParallel

[Adrian Croucher]

hi


That would great if it allows elements with mixed faces,


Yes, as I've already mentioned to Matt, I need this too for 6-node wedge 
elements (probably the same thing as your 'prisms').


- Adrian


Thank you
-Hassan

From: Matthew Knepley [mailto:knep...@gmail.com]
Sent: Thursday, June 29, 2017 10:36 AM
To: Hassan Raiesi 
Cc: petsc-users@mcs.anl.gov
Subject: Re: [petsc-users] specifying vertex coordinates using 
DMPlexCreateFromCellListParallel

On Thu, Jun 29, 2017 at 9:32 AM, Hassan Raiesi 
> 
wrote:
Ok,

The problem is that the faceSizes varies for pyramids and prisms, they have 
faces with variable number of nodes (i.e for pyramids, four triangular faces 
and one quad face, this I don?t see how to specify?
What if I report all faces as quad faces, for triangular faces, there would be 
a repeated vertex, would it make any problem?

No, you are right, that will not work. The whole interpolation code would have 
to be rewritten to allow faces with
variable size. Its not hard, but it is time-consuming. Since I am moving this 
summer, this will not be done in a timely
manner unfortunately.

Re: [petsc-users] Jacobian matrix for dual porosity model

[Adrian Croucher]

On 21/06/17 23:12, Matthew Knepley wrote:


From the above description, its not clear to me that you want this 
topology. For example, in the case that each cell gets an internal cell,
it seems like you could just handle this in the function space for 
each cell. Even for multiple internal cells, it seems like function space
parameterization is a better option. Topology is supposed to indicate 
the support of basis functions, but here that job is done. I would

just treat it as some sort of augmented approximation space.


If I understand what you mean, I considered doing something like that- 
basically just defining extra degrees of freedom in the cells where dual 
porosity is to be applied.


It seemed to me that if I then went ahead and created the Jacobian 
matrix using DMCreateMatrix(), it would give me extra nonzero entries 
that shouldn't be there - interactions between the dual porosity 
variables in neighbouring cells. Is there any way to avoid that?


- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] Jacobian matrix for dual porosity model

[Adrian Croucher]

On 16/06/17 11:37, Adrian Croucher wrote:

On 16/06/17 01:19, Matthew Knepley wrote:


Thanks for those ideas, very helpful.

If I try this approach (forming whole Jacobian matrix and using
PCFieldSplit Schur), I guess I will first need to set up a
modified DMPlex for the whole fracture + matrix mesh- so I can
use it to create vectors and the Jacobian matrix (with the right
sparsity pattern), and also to work out the coloring for finite
differencing.

Would that be straight-forward to do? Currently my DM just comes
from DMPlexCreateFromFile(). Presumably you can use
DMPlexInsertCone() or similar to add points into it?


You can certainly modify the mesh. I need to get a better idea what 
kind of modification, and then
I can suggest a way to do it. What do you start with, and what 
exactly do you want to add?



The way dual porosity is normally implemented in a finite volume 
context is to add an extra matrix rock cell 'inside' each of the 
original cells (which now represent the fractures, and have their 
volumes reduced accordingly), with a connection between the fracture 
cell and its corresponding matrix rock cell, so fluid can flow between 
them.


More generally there can be multiple matrix rock cells for each 
fracture cell, in which case further matrix rock cells are nested 
inside the first one, again with connections between them. There are 
formulae for computing the appropriate effective matrix rock cell 
volumes and connection areas, typically based on a 'fracture spacing' 
parameter which determines how fractured the rock is.


So in a DMPlex context it would mean somehow adding extra DAG points 
representing the internal cells and faces for each of the original 
cells. I'm not sure how that would be done.


I've been having a go at this- copying the cones and supports from the 
original DM into a new DM, which will represent the topology of the dual 
porosity mesh, but adding one new cell for each original cell, and a 
face connecting the two. I've based my attempt loosely on some of the 
stuff in SplitFaces() in TS ex11.c.


It isn't working though, as yet- when I construct the new DM and view 
it, the numbers of cells, faces etc. don't make sense, and the depth 
label says it's got 7 strata, where it should still just have 4 like the 
original DM.


I may have just messed it up somehow, but it occurred to me that maybe 
DMPlex won't like having bits of its DAG having different depths? With 
what I've been trying, the dual porosity parts of the mesh would just 
specify DAG points representing the new cells and the corresponding 
faces, but not for any edges or vertices of those faces. I hoped that 
would be sufficient for generating the nonzero structure of the Jacobian 
matrix. But I don't know if DMPlex will allow a DAG like that?


- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] Jacobian matrix for dual porosity model

[Adrian Croucher]

On 16/06/17 01:19, Matthew Knepley wrote:


Thanks for those ideas, very helpful.

If I try this approach (forming whole Jacobian matrix and using
PCFieldSplit Schur), I guess I will first need to set up a
modified DMPlex for the whole fracture + matrix mesh- so I can use
it to create vectors and the Jacobian matrix (with the right
sparsity pattern), and also to work out the coloring for finite
differencing.

Would that be straight-forward to do? Currently my DM just comes
from DMPlexCreateFromFile(). Presumably you can use
DMPlexInsertCone() or similar to add points into it?


You can certainly modify the mesh. I need to get a better idea what 
kind of modification, and then
I can suggest a way to do it. What do you start with, and what exactly 
do you want to add?



The way dual porosity is normally implemented in a finite volume context 
is to add an extra matrix rock cell 'inside' each of the original cells 
(which now represent the fractures, and have their volumes reduced 
accordingly), with a connection between the fracture cell and its 
corresponding matrix rock cell, so fluid can flow between them.


More generally there can be multiple matrix rock cells for each fracture 
cell, in which case further matrix rock cells are nested inside the 
first one, again with connections between them. There are formulae for 
computing the appropriate effective matrix rock cell volumes and 
connection areas, typically based on a 'fracture spacing' parameter 
which determines how fractured the rock is.


So in a DMPlex context it would mean somehow adding extra DAG points 
representing the internal cells and faces for each of the original 
cells. I'm not sure how that would be done.



Another approach, which might be easier, would be not to construct a DM 
for the whole dual-porosity mesh, but create the Jacobian matrix as a 
MatNest, with the fracture cell part created using DMCreateMatrix() as 
now, with the original DM, and the other parts of the Jacobian 
(representing fracture-matrix and matrix-matrix interactions) created 
'by hand'- because of the local one-dimensional nature of the matrix 
rock cell parts of the mesh, these would be just a bunch of block 
diagonal matrices anyway.


I've just been looking at SNES example ex70.c where something a bit like 
this is done.


I think once the Jacobian matrix is created, you can create a coloring 
for finite differencing on it just from the matrix itself, and don't 
actually need a corresponding DM. So this approach might work, without 
needing to construct a dual-porosity DM.


- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (

Re: [petsc-users] Jacobian matrix for dual porosity model

[Adrian Croucher]

On 14/06/17 07:45, Jed Brown wrote:

Barry Smith <bsm...@mcs.anl.gov> writes:


On Jun 13, 2017, at 10:06 AM, Jed Brown <j...@jedbrown.org> wrote:

Adrian Croucher <a.crouc...@auckland.ac.nz> writes:


One way might be to form the whole Jacobian but somehow use a modified
KSP solve which would implement the reduction process, do a KSP solve on
the reduced system of size n, and finally back-substitute to find the
unknowns in the matrix rock cells.

You can do this with PCFieldSplit type Schur, but it's a lot heavier
than you might like.

Is it clear that it would produce much overhead compared to doing a custom "reduction to a 
smaller problem". Perhaps he should start with this and then profiling can show if there are 
any likely benefits to "specializing more"?

Yeah, that would be reasonable.  We don't have a concept of sparsity for
preconditioners so don't have a clean way to produce the exact (sparse)
Schur complement.  Computing this matrix using coloring should be
relatively inexpensive due to the independence in each cell and its
tridiagonal structure.


Thanks for those ideas, very helpful.

If I try this approach (forming whole Jacobian matrix and using 
PCFieldSplit Schur), I guess I will first need to set up a modified 
DMPlex for the whole fracture + matrix mesh- so I can use it to create 
vectors and the Jacobian matrix (with the right sparsity pattern), and 
also to work out the coloring for finite differencing.


Would that be straight-forward to do? Currently my DM just comes from 
DMPlexCreateFromFile(). Presumably you can use DMPlexInsertCone() or 
similar to add points into it?


- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

[petsc-users] Jacobian matrix for dual porosity model

[Adrian Croucher]

 similar 
efficiency? I didn't see anything currently in PETSc specifically for 
solving block-tridiagonal systems.


Thanks, Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] DMPlex export to hdf5/vtk for triangle/prism mesh

[Adrian Croucher]

On 30/05/17 14:55, Matthew Knepley wrote:
On Mon, May 29, 2017 at 9:52 PM, Adrian Croucher 
<a.crouc...@auckland.ac.nz <mailto:a.crouc...@auckland.ac.nz>> wrote:


On 30/05/17 14:45, Matthew Knepley wrote:




What kind of basis are you expecting? A tensor product?


At present we don't even need basis functions, because we're just
doing flow simulation and it's all finite volume.


Okay good. Now for cell geometry. What kind of deformation do you 
allow in the wedge?


As in Fabian's application, these elements arise from meshes which have 
a simple layered structure in the vertical, but are unstructured in the 
horizontal (can be mixtures of quads and triangles in our case- in fact 
the triangles usually only occur where there is local refinement).


So for us these wedges are just horizontal triangles projected downwards 
in the vertical- not really deformed at all.



and what do you want
to know? For FV, we are providing the centroid and volume. If that is 
enough, we could be done quickly.


Yes, just centroid and volume would be enough.

- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] DMPlex export to hdf5/vtk for triangle/prism mesh

[Adrian Croucher]

On 30/05/17 14:45, Matthew Knepley wrote:



What kind of basis are you expecting? A tensor product?


At present we don't even need basis functions, because we're just doing 
flow simulation and it's all finite volume.


However further down the track we will also be doing rock mechanics on 
the same mesh, using finite elements. For that, tensor product basis 
would be fine.


- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] DMPlex export to hdf5/vtk for triangle/prism mesh

[Adrian Croucher]

On 30/05/17 12:25, Matthew Knepley wrote:


Sorry about not keeping up to date on that. I had not really thought 
about it working until Fabian suggested it.

So, it looks like XDMF output works. I am making a test now.

However, other stuff will not, like refinement, interpolation, cell 
geometry, and other discretization stuff.


What do you need working?


We'll definitely need interpolation and cell geometry, but that might be 
about it. We won't need refinement.


- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] DMPlex export to hdf5/vtk for triangle/prism mesh

[Adrian Croucher]

hi,

I was asking about support for exactly these 6-node wedge elements in 
DMPlex back in January.


At the time, there was no support for them. Has there been some progress 
since then?


We are going to need them before we can release our software, which 
we're aiming to do by the end of the year.


Cheers, Adrian


Message: 4
Date: Fri, 26 May 2017 22:40:40 -0500
From: Matthew Knepley 
To: "Fabian.Jakub" 
Cc: PETSc 
Subject: Re: [petsc-users] DMPlex export to hdf5/vtk for
triangle/prism mesh
Message-ID:

[petsc-users] SNES for 1D problem?

[Adrian Croucher]

hi

If part of my code needs a one-dimensional nonlinear solver (i.e. to 
solve for a single scalar variable), which for some problems will be 
called many times, does it still make sense to use SNES for that? or are 
there overheads from the vector machinery that make it less efficient 
for that case than rolling your own simple solver?


- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611

Re: [petsc-users] DMPlex hdf5 example

[Adrian Croucher]

hi


Does anybody have a simple example of HDF5 I/O for transient data with DMPlex?
I can?t figure out how to handle time steps: there does not seem to be a way to 
specify appending time steps when using PetscObjectViewFromOptions, and 
PetscViewerHDF5PushGroup and seem to not have any effect when viewing a Vec 
obtained tom a DMPlex.

I've done this, based on what's done in the Vec tutorial example ex19.c.

But to get the timesteps to work properly with DMPlex I replaced the 
PetscViewerHDF5SetTimestep() call with DMSetOutputSequenceNumber().


See thread here: 
http://lists.mcs.anl.gov/pipermail/petsc-users/2015-October/027353.html


- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 84611

Re: [petsc-users] block matrix in serial

[Adrian Croucher]

hi

On 18/05/16 02:03, Matthew Knepley wrote:


I have finally pushed what I think are the right fixes into master. 
Can you take a look and let me know

if it is fixed for you?


Yes, it looks like the block size is being set correctly now. Thanks 
very much!


Cheers, Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 84611

Re: [petsc-users] snes failures

[Adrian Croucher]

Might be unrelated, but I am also seeing a lot of SNES failures, on 
problems that used to converge fine, since I pulled the latest PETSc 
'next' branch yesterday.


This is using a FD Jacobian so it isn't a Jacobian coding problem. Same 
behaviour with LU linear solver, so it's not the linear solver either.


I'm trying to bisect and find out where it stopped working.

Cheers, Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 84611

Re: [petsc-users] block matrix in serial

[Adrian Croucher]

hi

Matt, did you (or anyone else) ever get a chance to look into these two bugs?

AFAICT the blocksize discovery code hasn't been changed since then.

I'm currently running a problem where the linear solver is struggling (running 
out of iterations),
and I'm wondering if it might perform better if the matrix blocksize were being 
set correctly.
May be nothing to do with it of course, but would at least be useful to 
eliminate that possibility.

Cheers, Adrian
//>//>/On 16/09/15 12:40, Matthew Knepley wrote: />/ > On Tue, Sep 15, 2015 at 9:05 PM, Adrian Croucher />/ > <https://lists.mcs.anl.gov/mailman/listinfo/petsc-users> 
<mailto:a.croucher at auckland.ac.nz 
<https://lists.mcs.anl.gov/mailman/listinfo/petsc-users>>> />/ > wrote: />/ > />/ > hi />/ > />/ > I have a test code (attached) that sets up a finite volume mesh />/ > using DMPlex, with 2 degrees of freedom per cell. />/ > />/ > I then create a matrix using DMCreateMatrix(), having used />/ > DMSetMatType() to set the matrix type to MATBAIJ or MATMPIBAIJ, to />/ > take advantage of the block structure. />/ > />/ > This works OK and gives me the expected matrix structure when I />/ > run on > 1 processor, but neither MATBAIJ or MATMPIBAIJ works if I />/ > run it in serial: />/ > />/ > />/ > Plex should automatically discover the block size from the Section. />/ > If not, it uses block size 1. I have to look at the example to see />/ > why the discovery is not working. Do you have any constraints? />//>/It looks like block discovery in parallel effectively always />/determines a block size of 1. Running with -mat_view ::ascii_info gives: />//>/Mat Object: 2 MPI processes />/type: mpibaij />/rows=20, cols=20 />/total: nonzeros=112, allocated nonzeros=112 />/total number of mallocs used during MatSetValues calls =0 />/block size is 1 />/^^^ />//>/The block size discovery does this: />//>/for (p = pStart; p < pEnd; ++p) { />/ierr = PetscSectionGetDof(sectionGlobal, p, );CHKERRQ(ierr); />/ierr = PetscSectionGetConstraintDof(sectionGlobal, p, />/);CHKERRQ(ierr); />/if (dof-cdof) { />/if (bs < 0) { />/bs = dof-cdof; />/} else if (bs != dof-cdof) { />//* Layout does not admit a pointwise block size */ />/bs = 1; />/break; />/} />/} />/} />//>/In parallel, some of the dofs are remote (encoded as negative). So we />/run through seeing (dof - cdof) == 2, until we hit a "remote" point at />/when we see (dof - cdof) != 2 and then we break out and set bs = 1. />//>/In serial, we correctly determine bs == 2. But then />/DMGetLocalToGlobalMapping always does />//>/ierr = ISLocalToGlobalMappingCreate(PETSC_COMM_SELF, 1,size, />/ltog, PETSC_OWN_POINTER, >ltogmap);CHKERRQ(ierr); />//>//>/i.e. is set with block size == 1. />//>/So there are two bugs here. />//>/1. In parallel, block size discovery in Mat creation has gone wrong />//

Crap. Hoist on my own petard. Okay I will fix this.



/2. (Always), the lgmap has block size of 1, irrespective of the />/discovered 
block size from the section. />//

Yep. This can also be fixed. It should work regardless, but would be better
with blocking.

  Thanks

 Matt

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 84611

Re: [petsc-users] SNES function domain error

[Adrian Croucher]

hi

On 03/12/15 18:43, Barry Smith wrote:

I absolutely do NOT recommend doing variable switching with SNES. Our 
Newton solvers are not designed with variable switching in mind so you get 
these absurd limitations like needing SNESLINESEARCHBASIC

It would be great if someone who understood variable switching would either 
1) contribute a Newton that supports variable switching to PETSc or 2) enhance 
the current SNESSolve_NewtonLS to support variable switching (I don't know 
which is easier).


I had resigned myself to not being able to use a line search, partly 
because the nonlinear solver in the TOUGH2 simulator (which we've used 
for many years to solve geothermal flow problems with phase changes and 
variable switching) doesn't use one either, and it generally works fine 
without one. Without the line search, SNES is doing much the same thing 
and doesn't seem to be upset by the variable switching either. It's only 
the error handling that's causing trouble.


But also I couldn't really see how a line search could be made to work 
with variable switching. When you switch variables, you're effectively 
doing a discontinuous leap into another part of the solution space. How 
can searching along a line between the old phase conditions and the new 
ones be made to make sense? Are there literature examples where people 
have done it?


Cheers, Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 84611

[petsc-users] SNES function domain error

[Adrian Croucher]

hi

I am trying to use SNESSetFunctionDomainError() to make SNES abort when 
my solution vector goes out of the function domain. However, it's not 
always working.


Specifically, I check the solution in a function called after the 
linesearch (set using SNESLineSearchSetPostCheck()). If it has gone out 
of the function domain I call SNESSetFunctionDomainError(). (I also 
check for phase changes at this point and do variable switching if 
necessary, which seems to work fine as long as I only use 
SNESLINESEARCHBASIC.)


This causes the linesearch reason to return SNES_LINESEARCH_FAILED_DOMAIN.

In SNESSolve_NEWTONLS() (ls.c:251) it checks the linesearch reason, but 
then also checks if (snes->stol*xnorm > ynorm). If that is true then it 
resets the snes->reason to SNES_CONVERGED_SNORM_RELATIVE and returns, 
i.e. it overrides the function domain error. So it thinks everything is 
fine, even though the 'converged' solution is out of the function domain.


Is this how it's meant to work? If so, what can I do to ensure the SNES 
always aborts if the solution goes out of the function domain? At 
present it's only working if (snes->stol*xnorm <= ynorm) as well.


Cheers, Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 84611

Re: [petsc-users] HDF5 vec output with DMPlex

[Adrian Croucher]

hi

When I read the HDF5 file back in again using PetscViewerHDF5Open() and 
VecLoad(), is there any easy way to navigate to the last time step in 
the results?


Or, equivalently, to find how many timesteps there are in there, so I 
could then use PetscViewerHDF5SetTimestep() before doing VecLoad() ?


Cheers, Adrian

On 17/10/15 03:01, Matthew Knepley wrote:


Now I remember. I did not want the output to depend on the viewer.

Does your example work if you replace PetscViewerHDF5SetTimestep() 
with DMSetOutputSequenceNumber()?


  Thanks,

Matt


--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 84611

[petsc-users] HDF5 vec output with DMPlex

[Adrian Croucher]

hi,

I am wanting to output time-stepped simulation results to an HDF5 file, 
from a PETSc Vec created on a DMPlex using DMCreateGlobalVector().


I tried doing it the way it's done in 
src/vec/vec/examples/tutorials/ex19.c, using PetscViewerHDF5Open(), 
PetscViewerHDF5PushGroup() and PetscViewerHDF5SetTimestep(), but the 
resulting HDF5 file only ever seems to have results for the last time 
step in it.


It works OK for a Vec created simply using VecCreate(), as in the ex19 
example. Is there any reason it shouldn't work with DMPlex?


Cheers, Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 84611

Re: [petsc-users] HDF5 vec output with DMPlex

[Adrian Croucher]

hi

On 16/10/15 16:54, Matthew Knepley wrote:


No, that should work. This is how we do it in PyLith. There must be 
some setup problem. When you look at the HDF5 file. Does

it have a dimension for timestep?


No, it doesn't.

Is there any extra setup I might need to do to get it to work? I've 
attached a minimal example code, mesh and output.


Cheers, Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 84611

program testview

  implicit none

#include 

  MPI_Comm :: comm
  PetscViewer :: viewer
  PetscInt :: i, dim
  PetscInt, parameter :: dof = 1
  PetscInt, parameter :: nt = 3
  PetscErrorCode :: ierr
  character(40), parameter :: filename = "column.exo"
  DM :: dm
  Vec :: X
  PetscFV :: fvm
  PetscDS :: ds

  call PetscInitialize(PETSC_NULL_CHARACTER, ierr); CHKERRQ(ierr)
  comm = PETSC_COMM_WORLD

  call DMPlexCreateExodusFromFile(comm, filename, PETSC_TRUE, dm, &
   ierr); CHKERRQ(ierr)

  call DMPlexSetAdjacencyUseCone(dm, PETSC_TRUE, ierr)
  CHKERRQ(ierr)
  call DMPlexSetAdjacencyUseClosure(dm, PETSC_FALSE, ierr)
  CHKERRQ(ierr)

  call PetscFVCreate(comm, fvm, ierr); CHKERRQ(ierr)
  call PetscFVSetNumComponents(fvm, dof, ierr); CHKERRQ(ierr)
  call DMGetDimension(dm, dim, ierr); CHKERRQ(ierr)
  call PetscFVSetSpatialDimension(fvm, dim, ierr); CHKERRQ(ierr)
  call DMGetDS(dm, ds, ierr); CHKERRQ(ierr)
  call PetscDSAddDiscretization(ds, fvm, ierr); CHKERRQ(ierr)

  call DMCreateGlobalVector(dm, X, ierr); CHKERRQ(ierr)
  call PetscObjectSetName(X, "vec", ierr); CHKERRQ(ierr)

  call PetscViewerHDF5Open(comm, "view.h5", FILE_MODE_WRITE, viewer, ierr)
  CHKERRQ(ierr)
  call PetscViewerHDF5PushGroup(viewer, "/", ierr); CHKERRQ(ierr)

  do i = 0, nt - 1
 call VecSet(X, dble(i), ierr); CHKERRQ(ierr)
 call PetscViewerHDF5SetTimestep(viewer, i, ierr); CHKERRQ(ierr)
 call VecView(X, viewer, ierr); CHKERRQ(ierr)
  end do

  call PetscViewerHDF5PopGroup(viewer, ierr); CHKERRQ(ierr)
  call PetscViewerDestroy(viewer, ierr); CHKERRQ(ierr)

  call VecDestroy(X, ierr); CHKERRQ(ierr)
  call DMDestroy(dm, ierr); CHKERRQ(ierr)
  call PetscFinalize(ierr); CHKERRQ(ierr)

end program testview


column.exo
Description: Binary data


view.h5
Description: HDF file

[petsc-users] block matrix in serial

[Adrian Croucher]

hi

I have a test code (attached) that sets up a finite volume mesh using 
DMPlex, with 2 degrees of freedom per cell.


I then create a matrix using DMCreateMatrix(), having used 
DMSetMatType() to set the matrix type to MATBAIJ or MATMPIBAIJ, to take 
advantage of the block structure.


This works OK and gives me the expected matrix structure when I run on > 
1 processor, but neither MATBAIJ or MATMPIBAIJ works if I run it in serial:


[0]PETSC ERROR: - Error Message 
--

[0]PETSC ERROR: Petsc has generated inconsistent data
[0]PETSC ERROR: Blocksize of layout 2 must match that of mapping 1
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html 
for trouble shooting.
[0]PETSC ERROR: Petsc Development GIT revision: v3.6.1-811-gaf4f2ea GIT 
Date: 2015-08-18 17:44:25 -0500
[0]PETSC ERROR: ./testmat on a linux-gnu-c-opt named des108 by acro018 
Wed Sep 16 13:51:27 2015
[0]PETSC ERROR: Configure options --download-netcdf --download-exodusii 
--with-hdf5-dir=/usr --download-triangle --download-ptscotch -

-download-chaco
[0]PETSC ERROR: #1 PetscLayoutSetISLocalToGlobalMapping() line 248 in 
/home/acro018/software/PETSc/code/src/vec/is/utils/pmap.c
[0]PETSC ERROR: #2 MatSetLocalToGlobalMapping() line 1876 in 
/home/acro018/software/PETSc/code/src/mat/interface/matrix.c
[0]PETSC ERROR: #3 DMCreateMatrix_Plex() line 853 in 
/home/acro018/software/PETSc/code/src/dm/impls/plex/plex.c
[0]PETSC ERROR: #4 DMCreateMatrix() line 961 in 
/home/acro018/software/PETSc/code/src/dm/interface/dm.c

--

It only works in serial if I use MATAIJ.

Is there a reason why the block matrix types don't appear to work in 
serial, at least in this case? (I'm running the PETSc 'next' branch.)


- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 84611

program testmat

  ! Test constructing matrix from DMPlex

  implicit none

#include 

  PetscInt, parameter :: dof = 2, overlap = 1
  MPI_Comm :: comm
  DM :: dm, dist_dm
  PetscErrorCode :: ierr
  character(40) :: filename = "col10.exo"
  PetscInt :: dim
  PetscFV :: fvm
  PetscDS :: ds
  Mat :: M

  call PetscInitialize(PETSC_NULL_CHARACTER, ierr); CHKERRQ(ierr)
  comm = PETSC_COMM_WORLD

  call DMPlexCreateExodusFromFile(comm, filename, PETSC_TRUE, dm, ierr); CHKERRQ(ierr)

  call DMPlexDistribute(dm, overlap, PETSC_NULL_OBJECT, dist_dm, ierr)
  CHKERRQ(ierr)
  if (dist_dm /= 0) then
 call DMDestroy(dm, ierr); CHKERRQ(ierr)
 dm = dist_dm
  end if

  call PetscFVCreate(comm, fvm, ierr); CHKERRQ(ierr)
  call PetscFVSetNumComponents(fvm, dof, ierr); CHKERRQ(ierr)
  call DMGetDimension(dm, dim, ierr); CHKERRQ(ierr)
  call PetscFVSetSpatialDimension(fvm, dim, ierr); CHKERRQ(ierr)
  call DMGetDS(dm, ds, ierr); CHKERRQ(ierr)
  call PetscDSAddDiscretization(ds, fvm, ierr); CHKERRQ(ierr)

  call DMSetMatType(dm, MATBAIJ, ierr); CHKERRQ(ierr)
  call DMCreateMatrix(dm, M, ierr); CHKERRQ(ierr)
  call MatView(M, PETSC_VIEWER_STDOUT_WORLD, ierr); CHKERRQ(ierr)

  call MatDestroy(M, ierr); CHKERRQ(ierr)
  call DMDestroy(dm, ierr); CHKERRQ(ierr)
  call PetscFinalize(ierr); CHKERRQ(ierr)

end program testmat



col10.exo
Description: Binary data

[petsc-users] DMPLex ghost indexing

[Adrian Croucher]

hi

I have a distributed DMPlex representing a finite volume mesh, and first 
I want to check if I understand the way the indexing of ghost cells 
works, based on my experiments.


As I understand it, there are two types of ghost cells:

1) those created by DMPlexDistribute() (which I call with overlap = 1), 
representing the cells over the edge of the processor partition boundary
2) those created by DMPlexConstructGhostCells(), which I call to add 
ghost cells on the physical boundary of the mesh (using a DMLabel), for 
applying boundary conditions


The type 1 'partition' ghost cells only appear in local vectors, not 
global, and are added on the end after the interior cells.


The type 2 'boundary' ghost cells appear in both local and global 
vectors. In a local vector they are added on the end after the partition 
ghost cells. In a global vector they appear straight after the interior 
cells.



Is this correct? If so, I have a question:

If I'm computing fluxes between my finite volume cells, I need a local 
vector to get access to both types of ghost cells.


But if I'm doing another calculation (computing fluid properties in the 
cells, in this case) that does not need access to the partition ghosts, 
but needs to be carried out for all interior cells *and* the boundary 
ghost cells, I figure it would be preferable to do this with a global 
vector (to save needless parallel communication scattering a global 
vector into a local one).


However if I use DMPlexGetHeightStratum() to get the start and end 
cells, the indices it returns apply to a local vector with the partition 
ghosts included.


Is it safe just to subtract off the number of partition ghost cells from 
the end cell index returned by DMPlexGetHeightStratum(), if I want to 
loop over all interior and boundary ghosts cells in a global vector?


Cheers, Adrian

--
Dr Adrian Croucher
Department of Engineering Science
University of Auckland
New Zealand
tel 64-9-373-7599 ext 84611

[petsc-users] DMDAVecGetArrayReadF90

[Adrian Croucher]

hi, I'm using a SNES to solve a nonlinear problem on a DMDA, using Fortran.

In my function to be minimized I need an array on the solution vector. 
If I use DMDAVecGetArrayF90(), I get an 'object is in wrong state' error 
because the solution vector has been locked read-only by the SNES.


If I'm using DMPLex I use VecGetArrayReadF90() to get a read-only array, 
but for DMDA there doesn't seem to be a corresponding 
DMDAVecGetArrayReadF90(). Is there something else I should use instead?


- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 84611

Re: [petsc-users] DMPlex and HDF5 vector ordering

[Adrian Croucher]

On 09/05/15 00:43, Matthew Knepley wrote:


Ata and Blaise have a pull request coming that creates a natural 
ordering for a Plex, similar to the
one used by DMDA, so you get output that is invariant to process 
number. It may take until the end

of the summer to get it fully integrated, but it is very close.


Great, that sounds just like what I need. I'll wait until I see it 
merged into next branch before giving it a go (not totally urgent right 
now).


Thanks!

- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 84611

[petsc-users] DMPlex and HDF5 vector ordering

[Adrian Croucher]

hi,

I create a Vec on a DMPlex using DMPlexCreateGlobalVec(), then write it 
to HDF5 using PetscViewerHDF5Open() and VecView().


I then try to read it back in later (in another program, but using the 
same DMPlex) using PetscViewerHDF5Open() and VecLoad().


It looks like the ordering of the final vector entries in the second 
program depends on how many processors I use. If they are the same in 
both programs, I get the right ordering, but if they aren't, I don't. Is 
that expected? If so, is there any way to guarantee the right ordering 
when I read the Vec back in?


- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 84611

[petsc-users] FD Dirichlet boundary conditions with ghosted DMDA

[Adrian Croucher]

hi,

I'm setting up a simple 1-D finite difference test problem with 
Dirichlet boundary conditions at each end, using DMDA.


I used DMDACreate1d() with the DM_BOUNDARY_GHOSTED option, and set the 
global array dimension to be the total number of nodes (including 
boundary nodes) minus 2- so the boundary nodes aren't included in the mesh.


When I'm evaluating the differential equation on this mesh, I use 
DMGetLocalVector() and DMGlobalToLocalBegin/End() to scatter to my 
ghosted vector, then DMDAVecGetArray() to get an array on the local vector.


When I want to apply the BCs at each end, I call DMDAGetGhostCorners() 
to find the indices of the ghost points.  I can now plug the required 
boundary values into my array.


At present, if I'm running on multiple processors, I'm just checking if 
rank=0 or rank=size-1 to make sure the BCs only get applied at the 
actual ends of the mesh, and not at the ends of the internal partitions 
as well.


This works, but is there a better way to do it? It seems potentially 
unreliable- as it assumes rank 0 is always the left hand end and rank 
size-1 always the right-hand end. Also it won't really work in 2D or 3D.


Thanks!

- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 84611

Re: [petsc-users] FD Dirichlet boundary conditions with ghosted DMDA

[Adrian Croucher]

   Use either the results for DMDAGetGhostCorners() or DMDAGetCorners() to 
determine if the first or last node is local. This extends trivially to 2 and 
3d. See for example src/ksp/ksp/examples/tutorials/ex34.c ComputeMatrix() if 
(i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1) { etc


Aha, of course. Thanks, that works great.

- Adrian

--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 84611