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

HDF5 corruption
by Adrian Croucher 13 May '24

13 May '24

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.croucher(a)auckland.ac.nz tel: +64 (0)9 923 4611

4 17

How to find the map between the high order coordinates of DMPlex and vertex numbering?
by Zongze Yang 15 May '23

15 May '23

Hi, if I load a `gmsh` file with second-order elements, the coordinates will be stored in a DG-P2 space. After obtaining the coordinates of a cell, how can I map the coordinates to vertex and edge? Below is some code load the gmsh file, I want to know the relation between `cl` and `cell_coords`. ``` import firedrake as fd import numpy as np # Load gmsh file (2rd) plex = fd.mesh._from_gmsh('test-fd-load-p2-rect.msh') cs, ce = plex.getHeightStratum(0) cdm = plex.getCoordinateDM() csec = dm.getCoordinateSection() coords_gvec = dm.getCoordinates() for i in range(cs, ce): cell_coords = cdm.getVecClosure(csec, coords_gvec, i) print(f'coordinates for cell {i} :\n{cell_coords.reshape([-1, 3])}') cl = dm.getTransitiveClosure(i) print('closure:', cl) break ``` Best wishes, Zongze

3 20

DMPlex PETSCViewer for Surface Component
by Mike Michell 30 Mar '23

30 Mar '23

Hi, I am a user of DMPlex object in 3D grid topology. Currently the solution field is printed out using viewer PETSCVIEWERVTK in .vtu format. By doing that the entire volume components are written in the file. I was wondering if there is an option that I can tell PETSc viewer to print out only surface component, instead of the entire volume. Best, Mike

3 4

Kokkos backend for Mat and Vec diverging when running on CUDA device.
by Fackler, Philip 27 Mar '23

27 Mar '23

This is an issue I've brought up before (and discussed in-person with Richard). I wanted to bring it up again because I'm hitting the limits of what I know to do, and I need help figuring this out. The problem can be reproduced using Xolotl's "develop" branch built against a petsc build with kokkos and kokkos-kernels enabled. Then, either add the relevant kokkos options to the "petscArgs=" line in the system test parameter file(s), or just replace the system test parameter files with the ones from the "feature-petsc-kokkos" branch. See here<https://github.com/ORNL-Fusion/xolotl/tree/feature-petsc-kokkos/benchmarks> the files that begin with "params_system_". Note that those files use the "kokkos" options, but the problem is similar using the corresponding cuda/cusparse options. I've already tried building kokkos-kernels with no TPLs and got slightly different results, but the same problem. Any help would be appreciated. Thanks, Philip Fackler Research Software Engineer, Application Engineering Group Advanced Computing Systems Research Section Computer Science and Mathematics Division Oak Ridge National Laboratory

3 16

DMPlex Reordering
by Nicholas Arnold-Medabalimi 08 Feb '23

08 Feb '23

Hi Petsc users, I have been working on how to read in meshes into a DMPlex object. The process of building the cones has been relatively straightforward. The mesh files in question have sets of faces that I use to derive the cell vertex cones. The method is basically identical to that used in DMPlexCreateFluent. After I setup the DMPlex cones and call Symmetrize and Stratify, I then load in all the coordinates corresponding to the vertices and then use DMInterpolate to generate the intermediate edges The issue that I am running into is that because I am deriving the cell-vertex relationship from independent sets of face-vertex relationships, I can end up with cells that have improper mesh ordering. For example a cell with the coordinates: point 0: 0.000000 0.000000 point 1: 0.000000 2.500000 point 2: 0.100000 0.000000 point 3: 0.100000 2.500000 As you can see instead of going around the perimeter of the cell the path from 1 to 2 instead bisects the cell. I can manually reorder these after I read in the coordinates by manually checking right-handedness but I was wondering if there is an easier way to reorder the cones? If there isn't once I do reorder the cones manually is there anything I need to do as far as station keeping on the DM? I apologize if I missed any resources on this. Thanks Nicholas -- Nicholas Arnold-Medabalimi Ph.D. Candidate Computational Aeroscience Lab University of Michigan

2 3

locate DMSwarm particles with respect to a background DMDA mesh
by Semplice Matteo 12 Jan '23

12 Jan '23

Dear Petsc developers, I am trying to use a DMSwarm to locate a cloud of points with respect to a background mesh. In the real application the points will be loaded from disk, but I have created a small demo in which * each processor creates Npart particles, all within the domain covered by the mesh, but not all in the local portion of the mesh * migrate the particles After migration most particles are not any more in the DMSwarm (how many and which ones seems to depend on the number of cpus, but it never happens that all particle survive the migration process). I am clearly missing some step, since I'd expect that a DMDA would be able to locate particles without the need to go through a DMShell as it is done in src/dm/tutorials/swarm_ex3.c.html<https://petsc.org/main/src/dm/tutorials/swarm_ex3.c.html> I attach my demo code. Could someone give me a hint? Best Matteo

5 18

PetscSF Fortran interface
by Nicholas Arnold-Medabalimi 09 Jan '23

09 Jan '23

Hi Petsc Users I'm in the process of adding some Petsc for mesh management into an existing Fortran Solver. It has been relatively straightforward so far but I am running into an issue with using PetscSF routines. Some like the PetscSFGetGraph work no problem but a few of my routines require the use of PetscSFGetLeafRanks and PetscSFGetRootRanks and those don't seem to be in the fortran interface and I just get a linking error. I also don't seem to see a PetscSF file in the finclude. Any clarification or assistance would be appreciated. Sincerely Nicholas -- Nicholas Arnold-Medabalimi Ph.D. Candidate Computational Aeroscience Lab University of Michigan

2 6

Petsc Section in DMPlex
by Nicholas Arnold-Medabalimi 08 Dec '22

08 Dec '22

Hi Petsc Users I have a question about properly using PetscSection to assign state variables to a DM. I have an existing DMPlex mesh distributed on 2 processors. My goal is to have state variables set to the cell centers. I then want to call DMPlexDistribute, which I hope will balance the mesh elements and hopefully transport the state variables to the hosting processors as the cells are distributed to a different processor count or simply just redistributing after doing mesh adaption. Looking at the DMPlex User guide, I should be able to achieve this with a single field section using SetDof and assigning the DOF to the points corresponding to cells. call DMPlexGetHeightStratum(dm,0,c0,c1,ierr) call DMPlexGetChart(dm,p0,p1,ierr) call PetscSectionCreate(PETSC_COMM_WORLD,section,ierr) call PetscSectionSetNumFields(section,1,ierr) call PetscSectionSetChart(section,p0,p1,ierr) do i = c0, (c1-1) call PetscSectionSetDof(section,i,nvar,ierr) end do call PetscSectionSetup(section,ierr) call DMSetLocalSection(dm,section,ierr) >From here, it looks like I can access and set the state vars using call DMGetGlobalVector(dmplex,state,ierr) call DMGetGlobalSection(dmplex,section,ierr) call VecGetArrayF90(state,stateVec,ierr) do i = c0, (c1-1) call PetscSectionGetOffset(section,i,offset,ierr) stateVec(offset:(offset+nvar))=state_i(:) !simplified assignment end do call VecRestoreArrayF90(state,stateVec,ierr) call DMRestoreGlobalVector(dmplex,state,ierr) To my understanding, I should be using Global vector since this is a pure assignment operation and I don't need the ghost cells. But the behavior I am seeing isn't exactly what I'd expect. To be honest, I'm somewhat unclear on a few things 1) Should be using nvar fields with 1 DOF each or 1 field with nvar DOFs or what the distinction between the two methods are? 2) Adding a print statement after the offset assignment I get (on rank 0 of 2) cell 1 offset 0 cell 2 offset 18 cell 3 offset 36 which is expected and works but on rank 1 I get cell 1 offset 9000 cell 2 offset 9018 cell 3 offset 9036 which isn't exactly what I would expect. Shouldn't the offsets reset at 0 for the next rank? 3) Does calling DMPlexDistribute also distribute the section data associated with the DOF, based on the description in DMPlexDistribute it looks like it should? I'd appreciate any insight into the specifics of this usage. I expect I have a misconception on the local vs global section. Thank you. Sincerely Nicholas -- Nicholas Arnold-Medabalimi Ph.D. Candidate Computational Aeroscience Lab University of Michigan

2 12

Using multiple MPI ranks with COO interface crashes in some cases
by Fackler, Philip 01 Dec '22

01 Dec '22

In Xolotl's "feature-petsc-kokkos" branch, I have moved our code to use the COO interface for preallocating and setting values in the Jacobian matrix. I have found that with some of our test cases, using more than one MPI rank results in a crash. Way down in the preconditioner code in petsc a Mat gets computed that has "null" for the "productsymbolic" member of its "ops". It's pretty far removed from where we compute the Jacobian entries, so I haven't been able (so far) to track it back to an error in my code. I'd appreciate some help with this from someone who is more familiar with the petsc guts so we can figure out what I'm doing wrong. (I'm assuming it's a bug in Xolotl.) Note that this is using the kokkos backend for Mat and Vec in petsc, but with a serial-only build of kokkos and kokkos-kernels. So, it's a CPU-only multiple MPI rank run. Here's a paste of the error output showing the relevant parts of the call stack: [ERROR] [0]PETSC ERROR: [ERROR] --------------------- Error Message -------------------------------------------------------------- [ERROR] [1]PETSC ERROR: [ERROR] --------------------- Error Message -------------------------------------------------------------- [ERROR] [1]PETSC ERROR: [ERROR] [0]PETSC ERROR: [ERROR] No support for this operation for this object type [ERROR] [1]PETSC ERROR: [ERROR] No support for this operation for this object type [ERROR] [0]PETSC ERROR: [ERROR] No method productsymbolic for Mat of type (null) [ERROR] No method productsymbolic for Mat of type (null) [ERROR] [0]PETSC ERROR: [ERROR] [1]PETSC ERROR: [ERROR] See https://petsc.org/release/faq/ for trouble shooting. [ERROR] See https://petsc.org/release/faq/ for trouble shooting. [ERROR] [0]PETSC ERROR: [ERROR] [1]PETSC ERROR: [ERROR] Petsc Development GIT revision: v3.18.1-115-gdca010e0e9a GIT Date: 2022-10-28 14:39:41 +0000 [ERROR] Petsc Development GIT revision: v3.18.1-115-gdca010e0e9a GIT Date: 2022-10-28 14:39:41 +0000 [ERROR] [1]PETSC ERROR: [ERROR] [0]PETSC ERROR: [ERROR] Unknown Name on a named PC0115427 by 4pf Mon Nov 14 13:22:01 2022 [ERROR] Unknown Name on a named PC0115427 by 4pf Mon Nov 14 13:22:01 2022 [ERROR] [1]PETSC ERROR: [ERROR] [0]PETSC ERROR: [ERROR] Configure options PETSC_DIR=/home/4pf/repos/petsc PETSC_ARCH=arch-kokkos-serial-debug --with-debugging=1 --with-cc=mpicc --with-cxx=mpicxx --with-fc=0 --with-cudac=0 --prefix=/home/4pf/build/petsc/serial-debug/install --with-64-bit-indices --with-shared-libraries --with-kokkos-dir=/home/4pf/build/kokkos/serial/install --with-kokkos-kernels-dir=/home/4pf/build/kokkos-kernels/serial/install [ERROR] Configure options PETSC_DIR=/home/4pf/repos/petsc PETSC_ARCH=arch-kokkos-serial-debug --with-debugging=1 --with-cc=mpicc --with-cxx=mpicxx --with-fc=0 --with-cudac=0 --prefix=/home/4pf/build/petsc/serial-debug/install --with-64-bit-indices --with-shared-libraries --with-kokkos-dir=/home/4pf/build/kokkos/serial/install --with-kokkos-kernels-dir=/home/4pf/build/kokkos-kernels/serial/install [ERROR] [1]PETSC ERROR: [ERROR] [0]PETSC ERROR: [ERROR] #1 MatProductSymbolic_MPIAIJKokkos_AB() at /home/4pf/repos/petsc/src/mat/impls/aij/mpi/kokkos/mpiaijkok.kokkos.cxx:918 [ERROR] #1 MatProductSymbolic_MPIAIJKokkos_AB() at /home/4pf/repos/petsc/src/mat/impls/aij/mpi/kokkos/mpiaijkok.kokkos.cxx:918 [ERROR] [1]PETSC ERROR: [ERROR] [0]PETSC ERROR: [ERROR] #2 MatProductSymbolic_MPIAIJKokkos() at /home/4pf/repos/petsc/src/mat/impls/aij/mpi/kokkos/mpiaijkok.kokkos.cxx:1138 [ERROR] #2 MatProductSymbolic_MPIAIJKokkos() at /home/4pf/repos/petsc/src/mat/impls/aij/mpi/kokkos/mpiaijkok.kokkos.cxx:1138 [ERROR] [1]PETSC ERROR: [ERROR] [0]PETSC ERROR: [ERROR] #3 MatProductSymbolic() at /home/4pf/repos/petsc/src/mat/interface/matproduct.c:793 [ERROR] #3 MatProductSymbolic() at /home/4pf/repos/petsc/src/mat/interface/matproduct.c:793 [ERROR] [1]PETSC ERROR: [ERROR] [0]PETSC ERROR: [ERROR] #4 MatProduct_Private() at /home/4pf/repos/petsc/src/mat/interface/matrix.c:9820 [ERROR] #4 MatProduct_Private() at /home/4pf/repos/petsc/src/mat/interface/matrix.c:9820 [ERROR] [0]PETSC ERROR: [ERROR] [1]PETSC ERROR: [ERROR] #5 MatMatMult() at /home/4pf/repos/petsc/src/mat/interface/matrix.c:9897 [ERROR] #5 MatMatMult() at /home/4pf/repos/petsc/src/mat/interface/matrix.c:9897 [ERROR] [0]PETSC ERROR: [ERROR] [1]PETSC ERROR: [ERROR] #6 PCGAMGOptProlongator_AGG() at /home/4pf/repos/petsc/src/ksp/pc/impls/gamg/agg.c:769 [ERROR] #6 PCGAMGOptProlongator_AGG() at /home/4pf/repos/petsc/src/ksp/pc/impls/gamg/agg.c:769 [ERROR] [0]PETSC ERROR: [ERROR] [1]PETSC ERROR: [ERROR] #7 PCSetUp_GAMG() at /home/4pf/repos/petsc/src/ksp/pc/impls/gamg/gamg.c:639 [ERROR] #7 PCSetUp_GAMG() at /home/4pf/repos/petsc/src/ksp/pc/impls/gamg/gamg.c:639 [ERROR] [1]PETSC ERROR: [ERROR] [0]PETSC ERROR: [ERROR] #8 PCSetUp() at /home/4pf/repos/petsc/src/ksp/pc/interface/precon.c:994 [ERROR] #8 PCSetUp() at /home/4pf/repos/petsc/src/ksp/pc/interface/precon.c:994 [ERROR] [1]PETSC ERROR: [ERROR] [0]PETSC ERROR: [ERROR] #9 KSPSetUp() at /home/4pf/repos/petsc/src/ksp/ksp/interface/itfunc.c:406 [ERROR] #9 KSPSetUp() at /home/4pf/repos/petsc/src/ksp/ksp/interface/itfunc.c:406 [ERROR] [1]PETSC ERROR: [ERROR] [0]PETSC ERROR: [ERROR] #10 KSPSolve_Private() at /home/4pf/repos/petsc/src/ksp/ksp/interface/itfunc.c:825 [ERROR] #10 KSPSolve_Private() at /home/4pf/repos/petsc/src/ksp/ksp/interface/itfunc.c:825 [ERROR] [0]PETSC ERROR: [ERROR] [1]PETSC ERROR: [ERROR] #11 KSPSolve() at /home/4pf/repos/petsc/src/ksp/ksp/interface/itfunc.c:1071 [ERROR] #11 KSPSolve() at /home/4pf/repos/petsc/src/ksp/ksp/interface/itfunc.c:1071 [ERROR] [1]PETSC ERROR: [ERROR] [0]PETSC ERROR: [ERROR] #12 PCApply_FieldSplit() at /home/4pf/repos/petsc/src/ksp/pc/impls/fieldsplit/fieldsplit.c:1246 [ERROR] #12 PCApply_FieldSplit() at /home/4pf/repos/petsc/src/ksp/pc/impls/fieldsplit/fieldsplit.c:1246 [ERROR] [1]PETSC ERROR: [ERROR] [0]PETSC ERROR: [ERROR] #13 PCApply() at /home/4pf/repos/petsc/src/ksp/pc/interface/precon.c:441 [ERROR] #13 PCApply() at /home/4pf/repos/petsc/src/ksp/pc/interface/precon.c:441 [ERROR] [1]PETSC ERROR: [ERROR] [0]PETSC ERROR: [ERROR] #14 KSP_PCApply() at /home/4pf/repos/petsc/include/petsc/private/kspimpl.h:380 [ERROR] #14 KSP_PCApply() at /home/4pf/repos/petsc/include/petsc/private/kspimpl.h:380 [ERROR] [1]PETSC ERROR: [ERROR] [0]PETSC ERROR: [ERROR] #15 KSPFGMRESCycle() at /home/4pf/repos/petsc/src/ksp/ksp/impls/gmres/fgmres/fgmres.c:152 [ERROR] #15 KSPFGMRESCycle() at /home/4pf/repos/petsc/src/ksp/ksp/impls/gmres/fgmres/fgmres.c:152 [ERROR] [1]PETSC ERROR: [ERROR] [0]PETSC ERROR: [ERROR] #16 KSPSolve_FGMRES() at /home/4pf/repos/petsc/src/ksp/ksp/impls/gmres/fgmres/fgmres.c:273 [ERROR] #16 KSPSolve_FGMRES() at /home/4pf/repos/petsc/src/ksp/ksp/impls/gmres/fgmres/fgmres.c:273 [ERROR] [1]PETSC ERROR: [ERROR] [0]PETSC ERROR: [ERROR] #17 KSPSolve_Private() at /home/4pf/repos/petsc/src/ksp/ksp/interface/itfunc.c:899 [ERROR] #17 KSPSolve_Private() at /home/4pf/repos/petsc/src/ksp/ksp/interface/itfunc.c:899 [ERROR] [0]PETSC ERROR: [ERROR] [1]PETSC ERROR: [ERROR] #18 KSPSolve() at /home/4pf/repos/petsc/src/ksp/ksp/interface/itfunc.c:1071 [ERROR] #18 KSPSolve() at /home/4pf/repos/petsc/src/ksp/ksp/interface/itfunc.c:1071 [ERROR] [0]PETSC ERROR: [ERROR] [1]PETSC ERROR: [ERROR] #19 SNESSolve_NEWTONLS() at /home/4pf/repos/petsc/src/snes/impls/ls/ls.c:210 [ERROR] #19 SNESSolve_NEWTONLS() at /home/4pf/repos/petsc/src/snes/impls/ls/ls.c:210 [ERROR] [1]PETSC ERROR: [ERROR] [0]PETSC ERROR: [ERROR] #20 SNESSolve() at /home/4pf/repos/petsc/src/snes/interface/snes.c:4689 [ERROR] #20 SNESSolve() at /home/4pf/repos/petsc/src/snes/interface/snes.c:4689 [ERROR] [1]PETSC ERROR: [ERROR] [0]PETSC ERROR: [ERROR] #21 TSStep_ARKIMEX() at /home/4pf/repos/petsc/src/ts/impls/arkimex/arkimex.c:791 [ERROR] #21 TSStep_ARKIMEX() at /home/4pf/repos/petsc/src/ts/impls/arkimex/arkimex.c:791 [ERROR] [1]PETSC ERROR: [ERROR] [0]PETSC ERROR: [ERROR] #22 TSStep() at /home/4pf/repos/petsc/src/ts/interface/ts.c:3445 [ERROR] #22 TSStep() at /home/4pf/repos/petsc/src/ts/interface/ts.c:3445 [ERROR] [1]PETSC ERROR: [ERROR] [0]PETSC ERROR: [ERROR] #23 TSSolve() at /home/4pf/repos/petsc/src/ts/interface/ts.c:3836 [ERROR] #23 TSSolve() at /home/4pf/repos/petsc/src/ts/interface/ts.c:3836 [ERROR] PetscSolver::solve: TSSolve failed. [ERROR] PetscSolver::solve: TSSolve failed. Aborting. Aborting. Thanks for the help, Philip Fackler Research Software Engineer, Application Engineering Group Advanced Computing Systems Research Section Computer Science and Mathematics Division Oak Ridge National Laboratory

4 16