April 2021 - petsc-users - mailman3.cels.anl.gov

PCFIELDSPLIT with BoomerAMG preconditioning on a Schur Complement
by Jorti, Zakariae 27 Apr '21

27 Apr '21

Hello, I have a block matrix A, [A_{00} A_{01}] [A_{10} A_{11}] , and am trying to use PCFIELDSPLIT as a preconditioner on it. The preconditioner has this shape: [inv(A_{00}) 0 ] [ 0 -ksp(S)] ; where S is the Schur Complement matrix: S = A_{11} - A_{10} * inv(A_{00}) * A_{01}, that I get by calling: MatGetSchurComplement(A, isrow0, iscol0, isrow1, iscol1, MAT_INITIAL_MATRIX, &S, MAT_SCHUR_COMPLEMENT_AINV_DIAG, MAT_IGNORE_MATRIX, NULL); I would like to set BoomerAMG preconditioner for ksp(S). So, I am calling: PCFieldSplitGetSubKSP(pc, &n, &subksp); KSPSetOperators(subksp[1], S, S); KSPGetPC(subksp[1], subpc); PCSetType(subpc,"boomeramg"); but I end up getting the following error: [0]PETSC ERROR: No support for this operation for this object type [0]PETSC ERROR: Mat type schurcomplement [0]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting. [0]PETSC ERROR: Petsc Release Version 3.14.2, Dec 03, 2020 [0]PETSC ERROR: ./mimeticcurleuler2 on a arch-linux-c-debug named nid00140 by zjorti Tue Apr 20 12:23:34 2021 [0]PETSC ERROR: Configure options CC=cc CXX=CC FC=ftn COPTFLAGS=-g CXXOPTFLAGS=-g FOPTFLAGS=-g --download-hypre --with-debugging=1 --with-cxx-dialect=C++11 [0]PETSC ERROR: #1 MatGetRow() line 561 in /global/u1/z/zjorti/Software/petsc-3.14.2/src/mat/interface/matrix.c [0]PETSC ERROR: #2 MatConvert_Basic() line 53 in /global/u1/z/zjorti/Software/petsc-3.14.2/src/mat/utils/convert.c [0]PETSC ERROR: #3 MatConvert() line 4254 in /global/u1/z/zjorti/Software/petsc-3.14.2/src/mat/interface/matrix.c [0]PETSC ERROR: #4 PCSetUp_HYPRE() line 231 in /global/u1/z/zjorti/Software/petsc-3.14.2/src/ksp/pc/impls/hypre/hypre.c [0]PETSC ERROR: #5 PCSetUp() line 1009 in /global/u1/z/zjorti/Software/petsc-3.14.2/src/ksp/pc/interface/precon.c [0]PETSC ERROR: #6 KSPSetUp() line 406 in /global/u1/z/zjorti/Software/petsc-3.14.2/src/ksp/ksp/interface/itfunc.c [0]PETSC ERROR: #7 KSPSolve_Private() line 658 in /global/u1/z/zjorti/Software/petsc-3.14.2/src/ksp/ksp/interface/itfunc.c [0]PETSC ERROR: #8 KSPSolve() line 889 in /global/u1/z/zjorti/Software/petsc-3.14.2/src/ksp/ksp/interface/itfunc.c [0]PETSC ERROR: #9 PCApply_FieldSplit_Schur() line 1133 in /global/u1/z/zjorti/Software/petsc-3.14.2/src/ksp/pc/impls/fieldsplit/fieldsplit.c [0]PETSC ERROR: #10 PCApply() line 444 in /global/u1/z/zjorti/Software/petsc-3.14.2/src/ksp/pc/interface/precon.c [0]PETSC ERROR: #11 KSP_PCApply() line 283 in /global/u1/z/zjorti/Software/petsc-3.14.2/include/petsc/private/kspimpl.h [0]PETSC ERROR: #12 KSPInitialResidual() line 65 in /global/u1/z/zjorti/Software/petsc-3.14.2/src/ksp/ksp/interface/itres.c [0]PETSC ERROR: #13 KSPSolve_GMRES() line 245 in /global/u1/z/zjorti/Software/petsc-3.14.2/src/ksp/ksp/impls/gmres/gmres.c [0]PETSC ERROR: #14 KSPSolve_Private() line 719 in /global/u1/z/zjorti/Software/petsc-3.14.2/src/ksp/ksp/interface/itfunc.c [0]PETSC ERROR: #15 KSPSolve() line 889 in /global/u1/z/zjorti/Software/petsc-3.14.2/src/ksp/ksp/interface/itfunc.c [0]PETSC ERROR: #16 SNESSolve_KSPONLY() line 51 in /global/u1/z/zjorti/Software/petsc-3.14.2/src/snes/impls/ksponly/ksponly.c [0]PETSC ERROR: #17 SNESSolve() line 4569 in /global/u1/z/zjorti/Software/petsc-3.14.2/src/snes/interface/snes.c [0]PETSC ERROR: #18 TSTheta_SNESSolve() line 185 in /global/u1/z/zjorti/Software/petsc-3.14.2/src/ts/impls/implicit/theta/theta.c [0]PETSC ERROR: #19 TSStep_Theta() line 223 in /global/u1/z/zjorti/Software/petsc-3.14.2/src/ts/impls/implicit/theta/theta.c [0]PETSC ERROR: #20 TSStep() line 3757 in /global/u1/z/zjorti/Software/petsc-3.14.2/src/ts/interface/ts.c [0]PETSC ERROR: #21 TSSolve() line 4154 in /global/u1/z/zjorti/Software/petsc-3.14.2/src/ts/interface/ts.c Am I doing something wrong here? How can I fix this? Many thanks. Best regards, Zakariae Jorti

5 5

Nesting splits based on IS
by Karin&NiKo 27 Apr '21

27 Apr '21

Dear PETSc users, I use the fieldsplit PC in an application where the splits are programmatically defined by IS using PCFieldSplitSetIS. Then the user can specify its own PC at runtime using PETSc options. My question : is it possible to define nested splits in this case as it can be done with strided splits (see snes/examples/tutorials/ex19.c with test suffix fieldsplit_4). In order to be perfectly clear : let's say I have a 3 fields problem : velocity (v split), pressure (p split) and temperature (t split). What I would like to do is something like the following but it fails : -ksp_type fgmres -pc_fieldsplit_type multiplicative -pc_type fieldsplit -pc_fieldsplit_0_fields fieldsplit_v_, fieldsplit_p_ -pc_fieldsplit_1_fields fieldsplit_t_ -prefix_push fieldsplit_0_ -ksp_type fgmres -pc_fieldsplit_schur_factorization_type upper -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_precondition a11 -prefix_pop -prefix_push fieldsplit_1_ -ksp_type fgmres -pc_type jacobi -prefix_pop -prefix_push fieldsplit_v_ -ksp_type fgmres -pc_type gamg -prefix_pop -prefix_push fieldsplit_p_ -ksp_type fgmres -pc_type jacobi -prefix_pop I thank you for your help, Nicolas

4 17

Re: [petsc-users] Can I assign different local number of row and column to the parallel matrix in different process?
by Barry Smith 25 Apr '21

25 Apr '21

> On Apr 25, 2021, at 10:04 PM, 王杰 <wangjieneu(a)hotmail.com> wrote: > > Thanks very much for your reply! > 1、 > In the first process, the global index starts from 0, and the local index also starts from 0. The MatCreateAij and MatSetValues work well. > However, in the second process, the minimum global index is 601, which doesn't coincide with the local index that starts from 0. Should I give a global_offset '601' when prepare the parameters *d_nnz and *o_nnz? No, d_nnz and *o_nnz don't have anything to do with global indices (that is determined completely by local indices set on each ran). > 2、 > In the first process，m and n are 1600 but in the second process they are 900. Is it ok to give different m and n in the two process which related to the same matrix d_A ? For example, (1600,1600)for the first process and (900,900) for the second process? Yes > 3、 > The global index in both process are not continuous. Is it necessory to use Application Orderings routines to make it continuous? PETSc global indices for matrices and vectors are always continuous; how you choose to manage your numbering of your own internal objects may or may not relate to the ordering PETSc users. When you use MatSetValue() and VecSetValues() you must use PETSc's continuous ordering. You can set MatSetLocalToGlobalMapping and VecSetLocalToGlobalMapping to allow you to use MatSetValuesLocal and VecSetValuesLocal using arbitrary local orderings. > 4、 > Is it necessory to difine the ghost nodes explicitily with the rountine DMDACreate2d(MPLComm comm, DMBoundaryType bx, DMBoundaryType by, DMDAStencilType stype, Petscint M, Petscint N, Petscint m, Petscint n, Petscint do£, Petscint s, canst Petscint lx[], canst Petscint ly[], DM '''da)? > Or I just to assign the nodes with the same global indices to both process, and the Petsc will generate ghost nodes automatically? If you use DMDACreate2d() which only makes sense for 2d structure grids the arguments determine the exact ghosting, basically the stencil width s and type of boundary conditions determine the ghosting. You cannot use another ghosting. But global indexing has nothing to do with ghosting and any process can use any global indices it would like and the Mat/Vec handle the needed communication. Barry > > 发件人: Barry Smith <bsmith(a)petsc.dev> > 发送时间: 2021年4月25日 17:22 > 收件人: 王杰 <wangjieneu(a)hotmail.com> > 抄送: petsc-users(a)mcs.anl.gov <petsc-users(a)mcs.anl.gov> > 主题: Re: [petsc-users] Can I assign different local number of row and column to the parallel matrix in different process? > > > Yes; but note that for square matrices using n == m almost always makes sense, otherwise the input vector to the matrix-vector product will have a different layout than the output vector. Code code is difficult to manage. > >> On Apr 23, 2021, at 9:29 PM, 王杰 <wangjieneu(a)hotmail.com <mailto:[email protected]>> wrote: >> >> hello, >> >> I use petsc 3.14.5 with OpenMPI 3.1.6 in Ubuntu 18.04.5. >> >> I don't well understand the routine MatCreateAIJ(MPI_Comm comm,PetscInt m, PetscInt n,PetscInt M,PetscInt N,PetscInt d_nz,PetscInt *d_nnz, PetscInt o_nz,PetscInt *o_nnz,Mat *A). >> >> I have tow process. The total number of global index is 2188. The first one deal with the global index 0~1400 with ghost node 1401~1500 and 2001~2188. The second process deal with globla index 1400~2100 with ghost node 601~700 and 701~800. This is dynamiclly assigned by load balancer. >> >> I want to assign the first 1400 rows to the first process and the remained 788 rows in the second process. Can I give different values to m and n in the routine MatCreateAIJin different process? >> >> Thanks. >> Wang Jie

1 0

Mat-Mat element wise multiplication
by Rui Silva 25 Apr '21

25 Apr '21

I am writing to you, to ask if there is any simple way to perform element-wise multiplication of 2 mpiaij matrices in Petsc? Furthermore is there any way to reassemble a matrix that has some zeros stored, to reduce the sparsity of the matrix? Thanks a lot in advance, R -- Dr. Rui Emanuel Ferreira da Silva Departamento de Física Teórica de la Materia Condensada Universidad Autónoma de Madrid, Spain Max Born Institute, Berlin, Germany https://ruiefdasilva.wixsite.com/ruiefdasilva https://mmuscles.eu/ https://www.mbi-berlin.de/index_en.html

2 1

Can I assign different local number of row and column to the parallel matrix in different process?
by 王杰 25 Apr '21

25 Apr '21

hello, I use petsc 3.14.5 with OpenMPI 3.1.6 in Ubuntu 18.04.5. I don't well understand the routine MatCreateAIJ(MPI_Comm comm,PetscInt m, PetscInt n,PetscInt M,PetscInt N,PetscInt d_nz,PetscInt *d_nnz, PetscInt o_nz,PetscInt *o_nnz,Mat *A). I have tow process. The total number of global index is 2188. The first one deal with the global index 0~1400 with ghost node 1401~1500 and 2001~2188. The second process deal with globla index 1400~2100 with ghost node 601~700 and 701~800. This is dynamiclly assigned by load balancer. I want to assign the first 1400 rows to the first process and the remained 788 rows in the second process. Can I give different values to m and n in the routine MatCreateAIJ in different process? Thanks. Wang Jie

2 1

Kokkos error at NERSC
by Mark Adams 24 Apr '21

24 Apr '21

I had Kokkos working here but I seemed to have broken it: global/homes/m/madams/petsc/arch-cori-gpu-opt-kokkos-gcc/include/setup/Kokkos_Setup_Cuda.hpp:54:2: error: #error "KOKKOS_ENABLE_CUDA defined but the compiler is not defining the __CUDACC__ macro as expected" #error \ Any ideas?

3 6

Adding a CUDA library
by Mark Adams 24 Apr '21

24 Apr '21

How would I add a CUDA library that is in the CUDA library directory? That is just add -lxxx to the link line. I tried LIBS=-lxxx and that failed in configure. Mark

3 3

VecLoad from HDF5 file
by Thibault Bridel-Bertomeu 23 Apr '21

23 Apr '21

Dear all, I have been facing a problem with VecLoad recently, even though it seems to me I did exactly like in the examples/tutorials. Basically, a program writes a vector with the HDF5 writer like this : call DMCreateGlobalVector(dm, sol, ierr); CHKERRA(ierr) call VecZeroEntries(X, ierr); CHKERRA(ierr) call PetscObjectSetName(X, "Solution", ierr); CHKERRA(ierr) < do something with X to fill it up with relevant data > call PetscViewerCreate(PETSC_COMM_WORLD, hdf5Viewer, ierr); CHKERRA(ierr) call PetscViewerSetType(hdf5Viewer, PETSCVIEWERHDF5, ierr); CHKERRA(ierr) call PetscViewerFileSetMode(hdf5Viewer, FILE_MODE_WRITE, ierr); CHKERRA(ierr); write(filename,'(A,I5.5,A)') "restart_", stepnum, ".h5" call PetscViewerFileSetName(hdf5Viewer, trim(filename), ierr); CHKERRA(ierr) call VecView(X, hdf5Viewer, ierr); CHKERRA(ierr) call PetscViewerDestroy(hdf5Viewer, ierr); CHKERRA(ierr) and the same program (but with different start-up options, say) re-reads such a file like this : call DMCreateGlobalVector(dm, sol, ierr); CHKERRA(ierr) call VecZeroEntries(sol, ierr); CHKERRA(ierr) call PetscObjectSetName(sol, "Solution", ierr); CHKERRA(ierr) call PetscViewerCreate(PETSC_COMM_WORLD, hdf5Viewer, ierr); CHKERRA(ierr) call PetscViewerSetType(hdf5Viewer, PETSCVIEWERHDF5, ierr); CHKERRA(ierr) call PetscViewerFileSetMode(hdf5Viewer, FILE_MODE_READ, ierr); CHKERRA(ierr) call PetscViewerFileSetName(hdf5Viewer, trim(restartname), ierr); CHKERRA(ierr) call VecLoad(sol, hdf5Viewer, ierr); CHKERRA(ierr) call PetscViewerDestroy(hdf5Viewer, ierr); CHKERRA(ierr) Such a dataset can be found under this link : https://drive.google.com/file/d/1owLAx5vknNhj61_5ieAwnWOR9cmkTseL/view?usp=… I would like the reader to read the /fields/Solution group basically, but I am not even sure it tries to do that. Anyhow, I got an error, saying that the size found in the file (105) does not match the expected size (25300). If I look at the shape of /fields/Solution it is given as (21, 5030, 5). First, it is weird, cause the 21 seems to be 1 + current iteration number ... but anyways we find the 5 variables and the 5030 cells. Only the reader seems to do 21 * 5 when it should be doing 5030 * 5 ... I tried adding 'PetscViewerHDF5PushGroup(hdf5Viewer, "/fields/Solution", ierr)' to force it to read that group, but it does not change anything. I would appreciate it if anyone could give me pointers on this issue ... Thank you very much in advance !! Thibault

4 20

Extracting solution at a boundary in a DMPlex/PetscFV paradigm
by Thibault Bridel-Bertomeu 22 Apr '21

22 Apr '21

Dear all, I am working with multi-D DMPlex's coupled with PetscFV and a TS to solve hyperbolic systems of equations. To add boundaries, I use the PetscDSAddBoundary as shown in some examples throughout the documentation. I do CFD, meaning I often have an object around which I want to simulate the flow. And I would like to know how some variables (pressure for instance) behave at the wall (literally _at_ the wall, and not only in the first fluid cell above the wall). After going through the documentation and the PETSc source code, I think the way to go is to use something like what ex23.c (of plex/tests) does : mark boundary faces that match the _wall_ boundary condition, create a submesh with them and somehow project the PetscFV/TS solution onto those facets based on the order of the spatial discretization (like if it's order 1 then just use the value in the first fluid cell above the boundary condition, if it's order 2 use the gradient and the value to get something at the facet, etc...). But I cannot figure out two things : 1/ how to mark the right boundary faces ? I tried using the DMPlexMarkBoundaryFaces with (what I thought was) the ID of the wall boundary condition (I mean, I DMView'd the DMPlex and I picked the right id in the Face Sets group), but then when I use the label thus created with DMPlexCreateSubmesh, it just returns ... nothing. The "submesh" is empty when I DMView it. 2/ how to project on the facets in a PetscFV paradigm ? Thank you very much in advance for any advice/pointers you could give me ! Take care and have a good day, Thibault

2 2

Reordering of entries in a DMDA
by Fabian Jakub 21 Apr '21

21 Apr '21

Dear Petsc team, I would like to play with the ordering of a DMDA application. To be exact, I would like to change the order in which a forward SOR/Gauss Seidel preconditioner iterates over the domain. e.g. change from [x,y,z,dof] to [z, yend->ystart, xend->xstart, dof] As I see it, I could a) reorder the application, i.e. how data is stored in the DMDA b) apply a mat reordering (happy to have the reordering only happening locally) c) or implement a custom PC (or using a PCShell?) I don't like a) because it would be very laborious and intrusive and I don't even know if it is worth it. Implementing a custom PC might give the most freedom but reduces solver configurability that I like. And how bad would a SOR be if it were to jump around in memory compared to the performance and memory needs of a mat reordering? I wanted to ask if you have any other suggestions and which approach you would recommend. Many thanks, Fabian

2 1