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

MatZeroRows changes my sparsity pattern
by Alexander Lindsay 15 Jul '21

15 Jul '21

My interpretation of the documentation page of MatZeroRows is that if I've set MAT_KEEP_NONZERO_PATTERN to true, then my sparsity pattern shouldn't be changed by a call to it, e.g. a->imax should not change. However, at least for sequential matrices, MatAssemblyEnd is called with MAT_FINAL_ASSEMBLY at the end of MatZeroRows_SeqAIJ and that does indeed change my sparsity pattern. Is my interpretation of the documentation page wrong? Alex

3 5

output DMDA to hdf5 file?
by Matteo Semplice 15 Jul '21

15 Jul '21

Dear all, I am experimenting with hdf5+xdmf output. At https://www.xdmf.org/index.php/XDMF_Model_and_Format I read that "XDMF uses XML to store Light data and to describe the data Model. Either HDF5[3] <https://www.hdfgroup.org/HDF5> or binary files can be used to store Heavy data. The data Format is stored redundantly in both XML and HDF5." However, if I call DMView(dmda,hdf5viewer) and then I run h5ls or h5stat on the resulting h5 file, I see no "geometry" section in the file. How should I write the geometry to the HDF5 file? Here below is what I have tried. Best Matteo //Setup ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE,DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, ctx.Nx,ctx.Ny, //global dim PETSC_DECIDE,PETSC_DECIDE, //n proc on each dim 2,stWidth, //dof, stencil width NULL, NULL, //n nodes per direction on each cpu &(ctx.daAll)); ierr = DMDASetUniformCoordinates(ctx.daAll, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0); CHKERRQ(ierr); ierr = DMDASetFieldName(ctx.daAll,0,"first"); CHKERRQ(ierr); ierr = DMDASetFieldName(ctx.daAll,1,"second"); CHKERRQ(ierr); ierr = DMDAGetLocalInfo(ctx.daAll,&ctx.daInfo); CHKERRQ(ierr); ierr = DMCreateFieldDecomposition(ctx.daAll,NULL, NULL, &ctx.is, &ctx.daField); CHKERRQ(ierr); //Initial data ierr = DMCreateGlobalVector(ctx.daAll,&ctx.U0); CHKERRQ(ierr); ierr = VecISSet(ctx.U0,ctx.is[0],1.); CHKERRQ(ierr); ierr = VecISSet(ctx.U0,ctx.is[1],100.0); CHKERRQ(ierr); Write mesh and Vec to hdf5: PetscViewer viewer; ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD,"solution.h5",FILE_MODE_WRITE,&viewer); CHKERRQ(ierr); ierr = DMView(ctx.daAll , viewer); CHKERRQ(ierr); //does not output anything to solution.h5?? ierr = VecView(ctx.U0,viewer); CHKERRQ(ierr); ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); Attempt to save the two fields separately:: PetscViewer viewer; ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD,"solution.h5",FILE_MODE_WRITE,&viewer); CHKERRQ(ierr); ierr = DMView(ctx.daField[0] , viewer); CHKERRQ(ierr); //does not output anything to solution.h5?? Vec uF; ierr = VecGetSubVector(ctx.U,ctx.is[0],&uF); CHKERRQ(ierr); PetscObjectSetName((PetscObject) uF, "first"); ierr = VecView(uF,viewer); CHKERRQ(ierr); ierr = VecRestoreSubVector(ctx.U,ctx.is[0],&uF); CHKERRQ(ierr); ierr = VecGetSubVector(ctx.U,ctx.is[1],&uF); CHKERRQ(ierr); PetscObjectSetName((PetscObject) uF, "second"); ierr = VecView(uF,viewer); CHKERRQ(ierr); ierr = VecRestoreSubVector(ctx.U,ctx.is[1],&uF); CHKERRQ(ierr); ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);

2 6

Loading PETSC Grid
by Alfredo J Duarte Gomez 14 Jul '21

14 Jul '21

Good morning, We are currently developing a PETSC application on a structured grid, however, we are applying a curvilinear transformation so that the physical domain is not a rectangle (the computational domain is). As a result we are trying to load in a mesh into the PETSC dmda object to create easy visualizations. I have successfully done this through the following code: ierr=VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,6400,&co);CHKERRQ(ierr); ierr=PetscObjectSetName((PetscObject) co, "coordinates");CHKERRQ(ierr); ierr=PetscViewerHDF5Open(PETSC_COMM_WORLD,"coord2d.h5",FILE_MODE_READ,&viewer);CHKERRQ(ierr); ierr=VecLoad(co,viewer);CHKERRQ(ierr); ierr=PetscViewerDestroy(&viewer);CHKERRQ(ierr); ierr=DMSetCoordinates(da,co);CHKERRQ(ierr); Where my coordinate vector in HDF5 format is a single column vector in the order x1,y1,x2,y2,...xN,yN,etc. where 1 is the lower left corner and N. However this seems like a rather crude way of doing it, it seems like using the DMGetCoordinateDM() and DMCreateGlobalVector() is a more sophisticated way of leveraging the dmda object, however when I use: DMGetCoordinateDM(da, &cda); DMCreateGlobalVector(cda, &co); ierr=PetscObjectSetName((PetscObject) co, "coordinates");CHKERRQ(ierr); ierr=PetscViewerHDF5Open(PETSC_COMM_WORLD,"coord2d.h5",FILE_MODE_READ,&viewer);CHKERRQ(ierr); ierr=VecLoad(co,viewer);CHKERRQ(ierr); ierr=PetscViewerDestroy(&viewer);CHKERRQ(ierr); ierr=DMSetCoordinates(da,co);CHKERRQ(ierr); The loading function complains about the structure/format of my data. I haven't been able to find the format that this coordinate dmda expects from an HDF5 file in any of the documentation. Is there anywhere I could find this information? Thank you, -Alfredo -- Alfredo Duarte Graduate Research Assistant The University of Texas at Austin

2 1

error on Spock
by Mark Adams 13 Jul '21

13 Jul '21

3 6

petsc-3.15.2 now available
by Satish Balay 10 Jul '21

10 Jul '21

Dear PETSc users, The patch release petsc-3.15.2 is now available for download. https://petsc.org/release/download/ Satish

1 0

Does petsc duplicate the users communicator?
by Kozdon, Jeremy (CIV) 09 Jul '21

09 Jul '21

Sorry if this is clearly stated somewhere in the docs, I'm still getting familiar with the petsc codebase and was also unable to find the answer searching (nor could I determine where this would be done in the source). Does petsc duplicate MPI communicators? Or does the users program need to make sure that the communicator remains valid for the life of a petsc object? The attached little test code seems to suggest that there is some duplication of MPI communicators behind the scenes. This came up when working on Julia wrappers for petsc. (Julia has a garbage collector so we need to make sure that references are properly kept if needed.)

3 4

best way to output in parallel data from DMDA (levelset) finite difference simulation
by Matteo Semplice 09 Jul '21

09 Jul '21

Dear all, it seems it should be a fairly straighforward thing to do but I am struggling with the output of my finite difference simulation. I have tried adapting my XML ascii VTK output routines that work nicely for finite volumes (but I have issues at points where 4 CPU subdomains touch), with an HDF5 PetscViewer (for which I cannot get a correct xdmf file) and few other random attempts, but neither got me a fully satisfactory result. Rather than correcting my attempts, I am ready to start afresh, also since this project will end up with much larger meshes than those I am used to. In your experience, what is the best way to output data - associated to a DMDA with more than one scalar fields, so that each variable can be visualized independently (dof 0 to scalar field "A", dof 1 to field "B", etc) - compatible with paraview (or visit or any other free tool on linux, if it need be) - with decent scaling, i.e. Vec data should be written directly by each CPU - binary format so that files are not incredibly huge is a plus - (bonus point) in the (not so) long run this will become a method for an arbitrary subdomain of the Cartesian grid defined by a level-set function, so a way to distinguish the points of the DMDA that are outside the physical domain would be a bonus (for example I could easily fill the Vec with NaNs for those points as long as the output/visualization can handle this). Thanks in advance! Best Matteo

1 0

CUDA running out of memory in PtAP
by Mark Adams 09 Jul '21

09 Jul '21

I am running out of memory in GAMG. It looks like this is from the new cuSparse RAP. I was able to run Hypre with twice as much work on the GPU as this run. Are there parameters to tweek for this perhaps or can I disable it? Thanks, Mark 0 SNES Function norm 5.442539952302e-04 [2]PETSC ERROR: --------------------- Error Message -------------------------------------------------------------- [2]PETSC ERROR: GPU resources unavailable [2]PETSC ERROR: CUDA error 2 (cudaErrorMemoryAllocation) : out of memory. Reports alloc failed; this indicates the GPU has run out resources [2]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting. [2]PETSC ERROR: Petsc Development GIT revision: v3.15.1-569-g270a066c1e GIT Date: 2021-07-06 03:22:54 -0700 [2]PETSC ERROR: ../ex2 on a arch-cori-gpu-opt-gcc named cgpu11 by madams Tue Jul 6 13:37:43 2021 [2]PETSC ERROR: Configure options --with-mpi-dir=/usr/common/software/sles15_cgpu/openmpi/4.0.3/gcc --with-cuda-dir=/usr/common/software/sles15_cgpu/cuda/11.1.1 --CFLAGS=" -g -DLANDAU_DIM=2 -DLANDAU_MAX_SPECI ES=10 -DLANDAU_MAX_Q=4" --CXXFLAGS=" -g -DLANDAU_DIM=2 -DLANDAU_MAX_SPECIES=10 -DLANDAU_MAX_Q=4" --CUDAFLAGS="-g -Xcompiler -rdynamic -DLANDAU_DIM=2 -DLANDAU_MAX_SPECIES=10 -DLANDAU_MAX_Q=4" --FFLAGS=" -g " - -COPTFLAGS=" -O3" --CXXOPTFLAGS=" -O3" --FOPTFLAGS=" -O3" --download-fblaslapack=1 --with-debugging=0 --with-mpiexec="srun -G 1" --with-cuda-gencodearch=70 --with-batch=0 --with-cuda=1 --download-p4est=1 -- download-hypre=1 --with-zlib=1 PETSC_ARCH=arch-cori-gpu-opt-gcc [2]PETSC ERROR: #1 MatProductSymbolic_SeqAIJCUSPARSE_SeqAIJCUSPARSE() at /global/u2/m/madams/petsc/src/mat/impls/aij/seq/seqcusparse/ aijcusparse.cu:2622 [2]PETSC ERROR: #2 MatProductSymbolic_ABC_Basic() at /global/u2/m/madams/petsc/src/mat/interface/matproduct.c:1159 [2]PETSC ERROR: #3 MatProductSymbolic() at /global/u2/m/madams/petsc/src/mat/interface/matproduct.c:799 [2]PETSC ERROR: #4 MatPtAP() at /global/u2/m/madams/petsc/src/mat/interface/matrix.c:9626 [2]PETSC ERROR: #5 PCGAMGCreateLevel_GAMG() at /global/u2/m/madams/petsc/src/ksp/pc/impls/gamg/gamg.c:87 [2]PETSC ERROR: #6 PCSetUp_GAMG() at /global/u2/m/madams/petsc/src/ksp/pc/impls/gamg/gamg.c:663 [2]PETSC ERROR: #7 PCSetUp() at /global/u2/m/madams/petsc/src/ksp/pc/interface/precon.c:1014 [2]PETSC ERROR: #8 KSPSetUp() at /global/u2/m/madams/petsc/src/ksp/ksp/interface/itfunc.c:406 [2]PETSC ERROR: #9 KSPSolve_Private() at /global/u2/m/madams/petsc/src/ksp/ksp/interface/itfunc.c:850 [2]PETSC ERROR: #10 KSPSolve() at /global/u2/m/madams/petsc/src/ksp/ksp/interface/itfunc.c:1084 [2]PETSC ERROR: #11 SNESSolve_NEWTONLS() at /global/u2/m/madams/petsc/src/snes/impls/ls/ls.c:225 [2]PETSC ERROR: #12 SNESSolve() at /global/u2/m/madams/petsc/src/snes/interface/snes.c:4769 [2]PETSC ERROR: #13 TSTheta_SNESSolve() at /global/u2/m/madams/petsc/src/ts/impls/implicit/theta/theta.c:185 [2]PETSC ERROR: #14 TSStep_Theta() at /global/u2/m/madams/petsc/src/ts/impls/implicit/theta/theta.c:223 [2]PETSC ERROR: #15 TSStep() at /global/u2/m/madams/petsc/src/ts/interface/ts.c:3571 [2]PETSC ERROR: #16 TSSolve() at /global/u2/m/madams/petsc/src/ts/interface/ts.c:3968 [2]PETSC ERROR: #17 main() at ex2.c:699

4 9

Scaling of the Petsc Binary Viewer
by Thibault Bridel-Bertomeu 09 Jul '21

09 Jul '21

Dear all, I have been having issues with large Vec (based on DMPLex) and massive MPI I/O ... it looks like the data that is written by the Petsc Binary Viewer is gibberish for large meshes split on a high number of processes. For instance, I am using a mesh that has around 50 million cells, split on 1024 processors. The computation seems to run fine, the timestep computed from the data makes sense so I think internally everything is fine. But when I look at the solution (one example attached) it's noise - at this point it should show a bow shock developing on the left near the step. The piece of code I use is below for the output : call DMGetOutputSequenceNumber(dm, save_seqnum, save_seqval, ierr); CHKERRA(ierr) call DMSetOutputSequenceNumber(dm, -1, 0.d0, ierr); CHKERRA(ierr) write(filename,'(A,I8.8,A)') "restart_", stepnum, ".bin" call PetscViewerCreate(PETSC_COMM_WORLD, binViewer, ierr); CHKERRA(ierr) call PetscViewerSetType(binViewer, PETSCVIEWERBINARY, ierr); CHKERRA(ierr) call PetscViewerFileSetMode(binViewer, FILE_MODE_WRITE, ierr); CHKERRA(ierr); call PetscViewerBinarySetUseMPIIO(binViewer, PETSC_TRUE, ierr); CHKERRA(ierr); call PetscViewerFileSetName(binViewer, trim(filename), ierr); CHKERRA(ierr) call VecView(X, binViewer, ierr); CHKERRA(ierr) call PetscViewerDestroy(binViewer, ierr); CHKERRA(ierr) call DMSetOutputSequenceNumber(dm, save_seqnum, save_seqval, ierr); CHKERRA(ierr) I do not think there is anything wrong with it but of course I would be happy to hear your feedback. Nonetheless my question was : how far have you tested the binary mpi i/o of a Vec ? Does it make some sense that for a 50 million cell mesh split on 1024 processes, it could somehow fail ? Or is it my python drawing method that is completely incapable of handling this dataset ? (paraview displays the same thing though so I'm not sure ...) Thank you very much for your advice and help !!! Thibault

3 6

[SLEPc] Computing Smallest Eigenvalue+Eigenvector of Many Small Matrices
by Jacob Faibussowitsch 08 Jul '21

08 Jul '21

Hello PETSc/SLEPc users, Similar to a recent question I am looking for an algorithm to compute the smallest eigenvalue and eigenvector for a bunch of matrices however I have a few extra “restrictions”. All matrices have the following properties: - All matrices are the same size - All matrices are small (perhaps no larger than 12x12) - All matrices are SPD - I only need the smallest eigenpair So far my best bet seems to be Lanczos but I’m wondering if there is some wunder method I’ve overlooked. Best regards, Jacob Faibussowitsch (Jacob Fai - booss - oh - vitch)

4 4