January 2025 - petsc-users - mailman3.cels.anl.gov

Fortran interfaces: Google Summer of Code 2025?
by Martin Diehl 18 Mar '25

18 Mar '25

Dear PETSc team, dear Barry, applications for the Google Summer of Code will start again and I was wondering if help for the re-factoring of the Fortran interfaces is still needed. Whether this makes sense depends on the progress of https://gitlab.com/petsc/petsc/-/merge_requests/7517 In contrast to the failed attempt last year, I have a student interested in working on this topic. Martin -- KU Leuven Department of Computer Science Department of Materials Engineering Celestijnenlaan 200a 3001 Leuven, Belgium

2 5

Visualizing higher order finite element output in ParaView
by Anna Dalklint 04 Feb '25

04 Feb '25

Hello, We have created a finite element code in PETSc for unstructured meshes using DMPlex. The first order meshes are created in gmsh and loaded into PETSc. To introduce higher order elements, e.g. 10 node tetrahedral elements, we start from scratch using PetscSection and loop over the relevant points it the DM to introduce additional degrees-of-freedom (example; for 10 node tets we have 4 vertices “nodes” and 6 edge “nodes”). The coordinates of the new “nodes” are obtained by interpolation using the finite element basis functions. The simulations seem to run well, but we face issues when trying to visualize the results in ParaView. We have tried to use both CGNS and HDF5+XDMF file formats for e.g. VecView. CGNS works, but the edge degrees-of-freedom appear to not be interpolated correctly (we observe oscillations in the fields, don’t know if this is a PETSc och ParaView issue). Also, we would prefer to use another file format than CGNS since it does not appear to directly allow timeseries (at least ParaView doesn’t recognize it). We haven’t got the HDF5+XDMF file format to work at all when running on more than one core (the mesh is highly distorted when saving using VecView and DMView + running the “petsc_gen_xdmf.py” script on the .h5 output file). VTU format works but then only the vertices’ degrees-of-freedom are visualized. As far as we have understood it, this is because VTU/VTK only supports degrees-of-freedom on vertices/cell level. Does anyone have any idea of how to visualize fields generated from higher order elements in ParaView? Or understand what we might be doing wrong? Best regards, Anna

4 9

local/global DMPlex Vec output
by Matteo Semplice 30 Jan '25

30 Jan '25

Dear petsc-users, I am having issues with output of parallel data attached to a DMPlex (or maybe more fundamental ones about DMPlex...). So I currently 1. create a DMPlex (DMPlexCreateGmshFromFile or DMPlexCreateBoxMesh) 2. partition it 3. and create a section for my data layout with DMPlexCreateSection(ctx.dmMesh, NULL, numComp, numDof, numBC, NULL, NULL, NULL, NULL, &sUavg) 4. DMSetLocalSection(ctx.dmMesh, sUavg) 5. create solLoc and solGlob vectors with DMCreateGlobalVector and DMCreateLocalVector 6. solve .... 7. VecView(ctx.solGlob, vtkViewer) on a .vtu file but when I load data in ParaView I get more cells than expected and it is as if the cells in the halo are put twice in output. (I could create a MWE if the above is not clear) I guess that the culprit is point (4), but if I replace it with DMSetGlobalSection then I cannot create the local vector at point (5). How should I handle this properly? In my code I need to create both local and global vectors, to perform at least GlobalToLocal and to save the global data. (On a side note, I tried also HDF5 but then it complains about the DM not having a DS...; really, any working solution that allows data to be explored with Paraview is fine) Thanks for any advice! Matteo Semplice

3 7

Update DMDA attached to DMSWARM
by MIGUEL MOLINOS PEREZ 28 Jan '25

28 Jan '25

Dear all, It seems that we can not modify the coordinates of the background mesh attached to a DMSwarm simulation. This is a feature? Best, Miguel

3 22

Re: [petsc-users] Copy dense matrix into half part of another dense matrix
by Pierre Jolivet 28 Jan '25

28 Jan '25

Please always keep the list in copy. The way you create A is not correct, I’ve attached a fixed code. If you want to keep your own distribution for A (and not the one associated to R_part), you’ll need to first call https://urldefense.us/v3/__https://petsc.org/main/manualpages/Mat/MatCreate… to redistribute A and then do a MatCopy() of the resulting Mat into R_part Thanks, Pierre $ /Volumes/Data/repositories/petsc/arch-darwin-c-debug-real/bin/mpirun -n 4 ./ex1234 Mat Object: 4 MPI processes type: mpidense 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 Mat Object: 4 MPI processes type: mpidense 2.6219599187040323e+00 1.9661197867318445e+00 1.5218640363910978e+00 3.5202261875977947e+00 3.6311893358251384e+00 2.2279492868785069e+00 2.7505403755038014e+00 3.1546072728892756e+00 1.8416294994524489e+00 2.4676055638467314e+00 2.3185625557889602e+00 2.0401666986599833e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 > On 27 Jan 2025, at 6:53 PM, medane.tchakorom(a)univ-fcomte.fr wrote: > > Re: > > This is a small reproductible example using MatDenseGetSubMatrix > > > Command: petscmpiexec -n 4 ./example > > ========================================================== > > PetscInt nlines = 8; // lines > PetscInt ncolumns = 3; // columns > PetscInt random_size = 12; > PetscInt rank; > PetscInt size; > > // Initialize PETSc > PetscInitialize(&argc, &args, NULL, NULL); > > MPI_Comm_rank(MPI_COMM_WORLD, &rank); > MPI_Comm_size(MPI_COMM_WORLD, &size); > > // R_full with all values to zero > Mat R_full; > MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, nlines, ncolumns, NULL, &R_full); > MatZeroEntries(R_full); > MatView(R_full, PETSC_VIEWER_STDOUT_WORLD); > > // Creating and setting A and S to rand values > Mat A, S; > MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, nlines / 2, random_size, NULL, &A); > MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, random_size, ncolumns, NULL, &S); > MatSetRandom(A, NULL); > MatSetRandom(S, NULL); > > // Computing R_part > Mat R_part; > MatDenseGetSubMatrix(R_full, PETSC_DECIDE, nlines / 2, PETSC_DECIDE, PETSC_DECIDE, &R_part); > MatMatMult(A, S, MAT_REUSE_MATRIX, PETSC_DECIDE, &R_part); > > // Visualizing R_full > MatDenseRestoreSubMatrix(R_full, &R_part); > MatView(R_full, PETSC_VIEWER_STDOUT_WORLD); > > // Destroying matrices > MatDestroy(&R_part); > MatDestroy(&R_full); > > PetscFinalize(); > return 0; > > ========================================================== > > > Part of the error output contains….: > > "Cannot change/reset row sizes to 1 local 4 global after previously setting them to 2 local 4 global ….” > > > > > ========================================================== > > PetscInt nlines = 8; // lines > PetscInt ncolumns = 3; // columns > PetscInt random_size = 12; > PetscInt rank; > PetscInt size; > > // Initialize PETSc > PetscInitialize(&argc, &args, NULL, NULL); > > MPI_Comm_rank(MPI_COMM_WORLD, &rank); > MPI_Comm_size(MPI_COMM_WORLD, &size); > > // R_full with all values to zero > Mat R_full; > MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, nlines, ncolumns, NULL, &R_full); > MatZeroEntries(R_full); > MatView(R_full, PETSC_VIEWER_STDOUT_WORLD); > > // Creating and setting A and S to rand values > Mat A, S; > MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, nlines / 2, random_size, NULL, &A); > MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, random_size, ncolumns, NULL, &S); > MatSetRandom(A, NULL); > MatSetRandom(S, NULL); > > // Computing R_part > Mat R_part; > MatMatMult(A, S, MAT_INITIAL_MATRIX, PETSC_DECIDE, &R_part); > MatView(R_part, PETSC_VIEWER_STDOUT_WORLD); > > Mat R_sub; > MatDenseGetSubMatrix(R_full, PETSC_DECIDE, nlines / 2, PETSC_DECIDE, PETSC_DECIDE, &R_sub); > > PetscScalar *storage = NULL; > MatDenseGetArray(R_part, &storage); > PetscScalar *storage_sub = NULL; > MatDenseGetArray(R_sub, &storage_sub); > > PetscArraycpy(storage_sub, storage, (nlines / 2) * ncolumns); > > MatDenseRestoreArray(R_part, &storage); > MatDenseRestoreArray(R_sub, &storage_sub); > > MatDenseRestoreSubMatrix(R_full, &R_sub); > > MatView(R_full, PETSC_VIEWER_STDOUT_WORLD); > > // Destroying matrices > MatDestroy(&R_part); > MatDestroy(&R_full); > > PetscFinalize(); > return 0; > ========================================================== > > > Now Using MatDenseGetArray > > > > Please let me know if I need to clarify something. > > > > > Thanks > Medane > > >> On 27 Jan 2025, at 16:26, Pierre Jolivet <pierre(a)joliv.et> wrote: >> >> >> >>> On 27 Jan 2025, at 3:52 PM, Matthew Knepley <knepley(a)gmail.com> wrote: >>> >>> On Mon, Jan 27, 2025 at 9:23 AM medane.tchakorom(a)univ-fcomte.fr <mailto:[email protected]> <medane.tchakorom(a)univ-fcomte.fr <mailto:[email protected]>> wrote: >>>> Re: >>>> >>>> MatDenseGetSubMatrix in fact could be the best alternative (cleaner code), but as mentioned earlier, I would like to use the smaller matrix R_part to get the result of a MatMatMult operation. >>>> >>>> MatMatMult(A, S, MAT_INITIAL_MATRIX, PETSC_DECIDE, &R_part); >>>> >>>> When trying to use MAT_REUSE_MATRIX, this gave an error (expected error I think). On the other side, I mentionned on MatDenseGetSubMatrix page, "The output matrix is not redistributed by PETSc”. So will R_part be a valid output matrix for MatMatMult? >>> >>> I believe it is, but I am not the expert. >>> >>> Pierre, can the SubMatrix be used as the output of a MatMatMult()? >> >> It should be, but there may be some limitations due to leading dimensions and what not. >> By looking at just the single line of code we got from you, I can see at least one issue: it should be MAT_REUSE_MATRIX, not MAT_INITIAL_MATRIX (assuming you got R_part from MatDenseGetSubMatrix). >> Feel free to share a (minimal) reproducer. >> >> Thanks, >> Pierre >> >>> Thanks, >>> >>> Matt >>> >>>> Thanks >>>> >>>> >>>>> On 27 Jan 2025, at 14:52, Matthew Knepley <knepley(a)gmail.com <mailto:[email protected]>> wrote: >>>>> >>>>> On Mon, Jan 27, 2025 at 8:42 AM Pierre Jolivet <pierre(a)joliv.et <mailto:[email protected]>> wrote: >>>>>> > On 27 Jan 2025, at 2:23 PM, medane.tchakorom(a)univ-fcomte.fr <mailto:[email protected]> wrote: >>>>>> > >>>>>> > Dear PETSc users, >>>>>> > >>>>>> > I hope this message finds you well. I don’t know If my question is relevant, but I’am currently working with DENSE type matrix, and would like to copy one matrix R_part [ n/2 x m] (resulted from a MatMatMult operation) into another dense matrix R_full [n x m]. >>>>>> > Both matrices being on the same communicator, I would like to efficiently copy R_part in the first half of R_full. >>>>>> > I have being using MatSetValues, but for large matrices, the subsequent assembling operation is costly. >>>>>> >>>>>> Could you please share the output of -log_view as well as a single file that will be generated with -info dump (e.g., dump.0, the file associated to process #0)? >>>>>> This shouldn’t be that costly, so there may be an option missing, like MAT_NO_OFF_PROC_ENTRIES. >>>>>> Anyway, if you want to optimize this, the fastest way would be to call MatDenseGetArray[Read,Write]() and then do the necessary PetscArraycpy(). >>>>> >>>>> The other alternative (which I think makes cleaner code) is to use >>>>> >>>>> https://urldefense.us/v3/__https://petsc.org/main/manualpages/Mat/MatDenseG… >>>>> >>>>> to create your R_part matrix. Then you are directly acting on the memory you want when assemble the smaller matrix. >>>>> >>>>> THanks, >>>>> >>>>> Matt >>>>> >>>>>> Thanks, >>>>>> Pierre >>>>>> >>>>>> > Please could you suggest me some strategies or functions to do this efficiently. >>>>>> > >>>>>> > Thank you for your time and assistance. >>>>>> > >>>>>> > Best regards, >>>>>> > Tchakorom Medane >>>>>> > >>>>>> >>>>> >>>>> >>>>> >>>>> -- >>>>> 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 >>>>> >>>>> https://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/__;fg!!G_uC… <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCf… > >>>> >>> >>> >>> >>> -- >>> 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 >>> >>> https://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/__;fg!!G_uC… <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCf… >

2 1

Copy dense matrix into half part of another dense matrix
by medane.tchakorom＠univ-fcomte.fr 27 Jan '25

27 Jan '25

Dear PETSc users, I hope this message finds you well. I don’t know If my question is relevant, but I’am currently working with DENSE type matrix, and would like to copy one matrix R_part [ n/2 x m] (resulted from a MatMatMult operation) into another dense matrix R_full [n x m]. Both matrices being on the same communicator, I would like to efficiently copy R_part in the first half of R_full. I have being using MatSetValues, but for large matrices, the subsequent assembling operation is costly. Please could you suggest me some strategies or functions to do this efficiently. Thank you for your time and assistance. Best regards, Tchakorom Medane

3 2

fortran TYPE IS statement vs petsc index set IS
by Klaij, Christiaan 26 Jan '25

26 Jan '25

In fortran I'm using the following structure to check the type of an incoming variable: SELECT TYPE (myvar) TYPE IS (mytype) ... END SELECT Here IS is a fortan intrinsic, so far so good. However, when I add a petsc index set as follows #include "petsc/finclude/petscksp.h" use petscksp, only: tIS IS :: myIS the compiler gets confused and thinks that the intrinsic fortran IS is the petsc index set IS, and errors-out on the SELECT TYPE: error #8245: SELECT TYPE statement must be immediately followed by CLASS IS, TYPE IS, CLASS DEFAULT or END SELECT statement. SELECT TYPE (myvar) ----^ error #6410: This name has not been declared as an array or a function. [TYPE] TYPE type(tIS) (mytype) ---------^ compilation aborted What would be the right way to deal with this problem? dr. ir. Christiaan Klaij | Senior Researcher | Research & Development T +31 317 49 33 44 | C.Klaij(a)marin.nl | https://urldefense.us/v3/__http://www.marin.nl__;!!G_uCfscf7eWS!ci7RiI8WEqh…

2 2

Issue when installing PETSc
by Nils Tilton 24 Jan '25

24 Jan '25

Dear PETSc Team, I hope this email finds you all well. I have a question regarding installation of PETSc. I recently installed PETSc on Ubuntu 20.04.2 LTS. I followed the online instructions, whereby I first downloaded the files using git, and then ran “./configure” and then “make all check.” However, I noted two issues during this process. First, after running ./configure, I got the following text: “Language used to compile PETSc: C” This gave me the impression PETSc was compiled with C, but my code is written in C++. Is the above going to be an issue? I am attaching my configure.log and make.log files if that helps. The second issue is that when I ran “make all check,” I did get some errors in the final checking stage. These don’t appear in the make.log file, so I am attaching a screenshot. I seem to recall that screenshots are frowned upon when asking questions to the PETSc team. I apologize in advance. I just couldn’t find any other way to save the text. I will add that I was able to successfully compile my C++ code that uses PETSc, but I am getting some funny behavior that could be related to the issues above. I am holding off on including those issues here to avoid complicating my question too soon. Thank you very much for all your help! Best Wishes, Nils

2 2

User Defined KSP Method in Fortran
by Frank Bramkamp 23 Jan '25

23 Jan '25

Dear PETSc team, I was planning to program a custom KSP method, some modified GMRES. We mainly use PETSc from Fortran. Therefore I wonder it is possible to have an interface to a custom KSP solver that is written in fortran. I thought of using KSPRegister to register my own routine, but that seems only available in C. Or is it possible to have a fortran/C wrapper to do that ? Thanks, Frank Bramkamp

3 4

Bug report: using a custom ksp convergence function
by Daniel Stone 21 Jan '25

21 Jan '25

Hello PETSc Community, I think I've found a bug - Go to $PETSC_DIR/src/ksp/ksp/tutorials/ open ex5f.F90, and add the following line in MyKSPConverged(), say around line 336: call KSPConvergedDefault(ksp,n,rnorm,flag,dummy,ierr) It should make sense why someone would want to do this - within a definition of custom convergence behaviour, get the default convergence flags, and, based on certain conditions, overwrite it. Now, building and running the exercise, making sure to include the flag to use the custom convergence: boston@boston-SYS-540A-TR:~/DATA_DRIVE/PETSC_DIRS/petsc_3.22/src/ksp/ksp/tutorials$ ./ex2f -my_ksp_convergence [0]PETSC ERROR: ------------------------------------------------------------------------ [0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range[...] >From my debugging, the issue is at line 1520 of iterativ.c: KSPConvergedDefaultCtx *cctx = (KSPConvergedDefaultCtx *)ctx; Because a dummy is used for the context object, this cast fails (my debugger says "cannot access memory..." for cctx). Then the crash happens when trying to access members of cctx: if (cctx->convmaxits && n >= ksp->max_it) { I'm not sure how to fix it. Crucially, we must note that that this did work in older versions of PETSc - I can repeat this test in 3.19.1, for example, and it works fine - a debugger shows the cast succeeding and producing some defaulted version of the ctx object, the subsequent crash does not happen, etc. Does anybody have any advice for working around this for the time being? A piece of software I work with uses a custom ksp convergence test in the manner described above, and is only functional with older versions of petsc because of this. Like in ex2f, I have no need of the ctx object, and just as in ex2f, I use 0 instead. Things I'll try next: defining a proper ctx object (I use 0, like in ex2f), or looking for some PETSC_NULL_CTX definition somewhere. Many Thanks, Daniel

3 6