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

How to combine different element types into a single DMPlex?
by Eric Chamberland 01 Aug '24

01 Aug '24

Hi, while playing with DMPlexBuildFromCellListParallel, I noticed we have to specify "numCorners" which is a fixed value, then gives a fixed number of nodes for a series of elements. How can I then add, for example, triangles and quadrangles into a DMPlex? Thanks, Eric -- Eric Chamberland, ing., M. Ing Professionnel de recherche GIREF/Université Laval (418) 656-2131 poste 41 22 42

4 13

Nodes coordinates in distributed dmplex
by Alejandro D Otero 16 Jul '22

16 Jul '22

Hi, sorry if this is an obvious question, but I cannot figure out how to recover finite element nodes coordinates once I have distributed a mesh stored as a dmplex. I am using petsc4py as interface to petsc rutines. I first created a dmplex using: dm.createFromCellList() In a sequential run I got the coordinates with: Coords = dm.getCoordinates() which gave a sequential vector with the coordinates of the mesh nodes. When I distribute the mesh with: dm.distribute() each mpi process has it own dm but the indexing of the vector resulting from getCoordinates() or getCoordinatesLocal() seems not consistent with the local numbering of the cells and nodes. Which is the correct way of doing this in PETSc philosophy? Thanks in advance, Alejandro

4 9

Question about PCFieldSplit
by Jorti, Zakariae 31 Jan '22

31 Jan '22

Hello, I have a DMStag grid with one dof on each edge and face center. I want to use a PCFieldSplit preconditioner on a Jacobian matrix that I assume is already split but I am not sure how to determine the fields. In the DMStag examples (ex2.c and ex3.c), the function PCFieldSplitSetDetectSaddlePoint is used to determine those fields based on zero diagonal entries. In my case, I have a Jacobian matrix that does not have zero diagonal entries. Can I use that PCFieldSplitSetDetectSaddlePoint in this case? If not, how should I do? Should I do like this example (https://www.mcs.anl.gov/petsc/petsc-master/src/ksp/ksp/tutorials/ex43.c.html) const PetscInt Bfields[1] = {0},Efields[1] = {1}; KSPGetPC(ksp,&pc); PCFieldSplitSetBlockSize(pc,2); PCFieldSplitSetFields(pc,"B",1,Bfields,Bfields); PCFieldSplitSetFields(pc,"E",1,Efields,Efields); where my B unknowns are defined on face centers and E unknowns are defined on edge centers? One last thing, I do not know which field comes first. Is it the one defined for face dofs or edge dofs. Thank you. Best regards, Zakariae

4 14

Is it possible to keep track of original elements # after a call to DMPlexDistribute ?
by Eric Chamberland 11 Nov '21

11 Nov '21

Hi, I want to use DMPlexDistribute from PETSc for computing overlapping and play with the different partitioners supported. However, after calling DMPlexDistribute, I noticed the elements are renumbered and then the original number is lost. What would be the best way to keep track of the element renumbering? a) Adding an optional parameter to let the user retrieve a vector or "IS" giving the old number? b) Adding a DMLabel (seems a wrong good solution) c) Other idea? Of course, I don't want to loose performances with the need of this "mapping"... Thanks, Eric -- Eric Chamberland, ing., M. Ing Professionnel de recherche GIREF/Université Laval (418) 656-2131 poste 41 22 42

2 31

SLEPc: smallest eigenvalues
by Varun Hiremath 05 Oct '21

05 Oct '21

Hi All, I am trying to compute the smallest eigenvalues of a generalized system A*x= lambda*B*x. I don't explicitly know the matrix A (so I am using a shell matrix with a custom matmult function) however, the matrix B is explicitly known so I compute inv(B)*A within the shell matrix and solve inv(B)*A*x = lambda*x. To compute the smallest eigenvalues it is recommended to solve the inverted system, but since matrix A is not explicitly known I can't invert the system. Moreover, the size of the system can be really big, and with the default Krylov solver, it is extremely slow. So is there a better way for me to compute the smallest eigenvalues of this system? Thanks, Varun

2 20

PETSC installation on Cray
by Enrico 13 Sep '21

13 Sep '21

Hi, I'm having some problems installing PETSC with Cray compiler. I use this configuration: ./configure --with-cc=cc --with-cxx=CC --with-fc=0 --with-debugging=1 --with-shared-libraries=1 COPTFLAGS=-O0 CXXOPTFLAGS=-O0 and when I do make all I get the following error because of cmathcalls.h: CC-1043 craycc: ERROR File = /usr/include/bits/cmathcalls.h, Line = 55 _Complex can only be used with floating-point types. __MATHCALL (cacos, (_Mdouble_complex_ __z)); ^ Am I doing something wrong? Regards, Enrico Degregori

4 5

Mat preallocation in case of variable stencils
by Matteo Semplice 07 Sep '21

07 Sep '21

Hi. We are writing a code for a FD scheme on an irregular domain and thus the local stencil is quite variable: we have inner nodes, boundary nodes and inactive nodes, each with their own stencil type and offset with respect to the grid node. We currently create a matrix with DMCreateMatrix on a DMDA and for now have set the option MAT_NEW_NONZERO_LOCATIONS to PETSC_TRUE, but its time to render the code memory-efficient. The layout created automatically is correct for inner nodes, wrong for boundary ones (off-centered stencils) and redundant for outer nodes. After the preprocessing stage (including stencil creation) we'd be in position to set the nonzero pattern properly. Do we need to start from a Mat created by CreateMatrix? Or is it ok to call DMCreateMatrix (so that the splitting among CPUs and the block size are set by PETSc) and then call a MatSetPreallocation routine? Also, I've seen in some examples that you call the Seq and the MPI preallocation routines in a row. Does this work because MatMPIAIJSetPreallocation silently does nothing on a Seq matrix and viceversa? Thanks Matteo

4 5

PETSc 3.15.3 compiling error
by Sam Guo 01 Sep '21

01 Sep '21

Dear PETSc dev team, I am compiling petsc 3.15.3 and got following compiling error petsc/src/mat/impls/aij/mpi/mumps/mumps.c:52:31: error: missing binary operator before token "(" 52 | #if PETSC_PKG_MUMPS_VERSION_GE(5,3,0) Any idea what I did wrong? Thanks, Sam

5 32

Postdoctoral position at ETH Zurich: Geodynamics / HPC / Julia
by Patrick Sanan 31 Aug '21

31 Aug '21

The Geophysical Fluid Dynamics group at ETH Zurich (Switzerland) is seeking a postdoctoral appointee to work for about 2.5 years on an ambitious project involving developing a Julia-based library for GPU-accelerated multiphysics solvers based on pseudotransient relaxation. Of particular interest for this audience might be that a major component of the proposed work is to make these solvers available via PETSc (as a SNES implementation), thus exposing them for use within a host of existing HPC applications, including those involved in this specific project. We'll accept applications until the position is filled, but for full consideration please apply before October 1, 2021. Full information is in the ad at the following link, and please feel free to contact me directly! https://github.com/psanan/gpu4geo_postdoc_ad/ Best, Patrick

1 0

KSPSolve with MPIAIJ with non-square 'diagonal parts'
by Olivier Jamond 30 Aug '21

30 Aug '21

Hello, I am sorry because I surely miss something, but I cannot manage to solve a problem with a MPIAIJ matrix which has non-square 'diagonal parts'. I copy/paste at the bottom of this message a very simple piece of code which causes me troubles. I this code I try to do 'x=KSP(A)*b' (with gmres/jacobi), but this fails whereas a matrix multiplication 'b=A*x' seems to work. Is ksp with such a matrix supposed to work (I can't find anything in the documentation about that, so I guess that it is...) ? Many thanks, Olivier NB: this code should be launched with exactly 3 procs #include "petscsys.h" /* framework routines */ #include "petscvec.h" /* vectors */ #include "petscmat.h" /* matrices */ #include "petscksp.h" #include <vector> #include <string> #include <iostream> #include <numeric> static char help[] = "Trying to solve a linear system on a sub-block of a matrix\n\n"; int main(int argc, char **argv) { MPI_Init(NULL, NULL); PetscErrorCode ierr; ierr = PetscInitialize(&argc, &argv, NULL, help); CHKERRQ(ierr); // clang-format off std::vector<std::vector<double>> AA = { { 1, 2, 0, /**/ 0, 3, /**/ 0, 0, 4}, { 0, 5, 6, /**/ 7, 0, /**/ 0, 8, 0}, { 9, 0, 10, /**/11, 0, /**/ 0, 12, 0}, //--------------------------------------- {13, 0, 14, /**/15, 16, /**/17, 0, 0}, { 0, 18, 0, /**/19, 20, /**/21, 0, 0}, { 0, 0, 0, /**/22, 23, /**/ 1, 24, 0}, //-------------------------------------- {25, 26, 27, /**/ 0, 0, /**/28, 29, 0}, {30, 0, 0, /**/31, 32, /**/33, 0, 34}, }; std::vector<double> bb = {1., 1., 1., // 1., 1., 1., // 1., 1.}; std::vector<int> nDofsRow = {3, 3, 2}; std::vector<int> nDofsCol = {3, 2, 3}; // clang-format on int NDofs = std::accumulate(nDofsRow.begin(), nDofsRow.end(), 0); int pRank, nProc; MPI_Comm_rank(PETSC_COMM_WORLD, &pRank); MPI_Comm_size(PETSC_COMM_WORLD, &nProc); if (nProc != 3) { std::cerr << "THIS TEST MUST BE LAUNCHED WITH EXACTLY 3 PROCS\n"; abort(); } Mat A; MatCreate(PETSC_COMM_WORLD, &A); MatSetType(A, MATMPIAIJ); MatSetSizes(A, nDofsRow[pRank], nDofsCol[pRank], PETSC_DETERMINE, PETSC_DETERMINE); MatMPIAIJSetPreallocation(A, NDofs, NULL, NDofs, NULL); Vec b; VecCreate(PETSC_COMM_WORLD, &b); VecSetType(b, VECMPI); VecSetSizes(b, nDofsRow[pRank], PETSC_DECIDE); if (pRank == 0) { for (int i = 0; i < NDofs; ++i) { for (int j = 0; j < NDofs; ++j) { if (AA[i][j] != 0.) { MatSetValue(A, i, j, AA[i][j], ADD_VALUES); } } VecSetValue(b, i, bb[i], ADD_VALUES); } } MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY); VecAssemblyBegin(b); VecAssemblyEnd(b); PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_DENSE); MatView(A, PETSC_VIEWER_STDOUT_WORLD); VecView(b, PETSC_VIEWER_STDOUT_WORLD); KSP ksp; KSPCreate(PETSC_COMM_WORLD, &ksp); KSPSetOperators(ksp, A, A); KSPSetFromOptions(ksp); PC pc; KSPGetPC(ksp, &pc); PCSetFromOptions(pc); Vec x; MatCreateVecs(A, &x, NULL); ierr = KSPSolve(ksp, b, x); // this fails MatMult(A, x, b); // whereas the seems to be ok... VecView(x, PETSC_VIEWER_STDOUT_WORLD); MPI_Finalize(); return 0; }

3 2