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

Calling for short user presentations in a PETSc BoF
by Junchao Zhang 19 Dec '25

19 Dec '25

Dear PETSc/TAO Community, We are soliciting PETSc users to share their usage experiences, application successes, and ongoing challenges in an online Zoom Birds-of-a-Feather (BoF) session, to be held between February 10~12, 2026. We are seeking approximately five short user presentations, each consisting of a 5-minute talk followed by 2 minutes of questions. If you are interested in presenting, please contact petsc-maint(a)mcs.anl.gov with your talk title, a brief abstract, and your preferred time slot (11:00 AM, 1:00 PM, or 3:00 PM EST). The BoF is hosted by the Consortium for the Advancement of Scientific Software (CASS <https://urldefense.us/v3/__https://cass.community/__;!!G_uCfscf7eWS!cGKRFA2… >) and led by the PESO <https://urldefense.us/v3/__https://pesoproject.org/__;!!G_uCfscf7eWS!cGKRFA… > (Partnering for Scientific Software Ecosystem Stewardship) project. The PETSc session will last 90 minutes and will take place on one day between February 10 and 12, 2026 (exact date to be finalized). Please note that the session will not be recorded. Our preferred time slot is 11:00 AM EST (5:00 PM UTC) to better accommodate European participants, although alternative options at 1:00 PM or 3:00 PM EST are also under consideration. In addition to user presentations, during the session, PETSc developers will highlight recent advances developed following the Exascale Computing Project, including the new PETSc Fortran bindings, PetscRegressor, TaoTerm, updates to PETSc GPU backends, mixed-precision support in PETSc/MUMPS, and integration with OpenFOAM, among other topics. The program will also include an open discussion of emerging PETSc research directions, such as leveraging agentic artificial intelligence to enhance and exploit the PETSc knowledge base. The BoF will provide insight into PETSc’s near-term development roadmap and offer a forum for user feedback on desired features and improvements. Active participation and questions from the audience are strongly encouraged, enabling the PETSc team to better align future development with community needs. The agenda will be posted once the program is finalized. Thank you, and we look forward to your participation. Junchao Zhang On behalf of the PETSc team

1 0

Handling inactive (zero-occupancy) equations in large SNES systems
by MIGUEL MOLINOS PEREZ 16 Dec '25

16 Dec '25

Dear all, I am working with a large nonlinear system solved with SNES, where a significant fraction of the unknowns are temporarily inactive due to a physical parameter being zero (e.g. zero occupancy / zero weight). For those DOF the corresponding equilibrium equation is physically inactive, but the unknown still appears in the global vector and in couplings of neighboring particles (Im using dmswarm). At the moment, these inactive equations contribute with a zero residual (F_i=0), which (I think) leads to poor conditioning and convergence issues for large problems. My question is about best numerical practice in this situation. For the position field, should I do something like F_i = q_i - q_(i,n)? Where q_(i,n) is the position of the particle at the previous configuration. Best regards, Miguel

3 3

Trouble when viewing a subDM in vtk format
by Aldo Bonfiglioli 15 Dec '25

15 Dec '25

Dear developers, I wrote a code that extracts subDMs corresponding to the various strata in the Face Sets. I run into troubles when I view a subDM or a Vec attached to the subDM using the VTK format. More precisely, the problem only occurs on more than one processor when the rank=0 processor has no points on a given subDM. For instance, when the attached reproducer is run on 2 procs, the u_01.vtu file (global u Vec mapped to the subDM corresponding to stratum=1) only includes the header, but no data. All other u_0?.vtu files can successfully be loaded and viewed in paraview. The problem does NOT arise when I view the same objects in HDF5 format. However, my problem in using the HDF5 lies in the fact that: while the hdf5 file obtained with DMView can be post-processed with "petsc/lib/petsc/bin/petsc_gen_xdmf.py" to create a xmf file readable by paraview I do not know how to view the field(s) associated with the DM when the hdf5 file is obtained from VecView. The reproducer compiles with the latest petsc release. Thanks, Aldo -- Dr. Aldo Bonfiglioli Associate professor of Fluid Mechanics Dipartimento di Ingegneria Universita' della Basilicata V.le dell'Ateneo Lucano, 10 85100 Potenza ITALY tel:+39.0971.205203 fax:+39.0971.205215 web: https://urldefense.us/v3/__http://docenti.unibas.it/site/home/docente.html?…

2 3

Use of flag dm_plex_high_order_view
by Noam T. 11 Dec '25

11 Dec '25

Hello, In an old question about obtaining the node connectivity of a cell with a high order approximation space, the use of the flag "-dm_plex_high_order_view" for visualization purposes was brought up, as a way to refine the grid down to linear space. Link: https://urldefense.us/v3/__https://gitlab.com/petsc/petsc/-/blob/main/src/d… I'm trying to make use of this flag to see the refinement, but I see no difference with higher order approximations. Perhaps I am misunderstanding its use? I thought that by using it, one could see a "subdivision" of each element. Say, a single triangle, FE approximation space order 2 (3 corner nodes, 3 mid-edge nodes), would be refined into e.g. 4 linear triangles. The code (Fortran) looks something along these lines: ! Create a DM from the mesh file DM :: cdm PetscInt :: K, cdim, fedim PetscBool :: simplex PetscFE :: cfe PetscViewer :: viewer PetscErrorCode :: ierr PetsCallA(DMSetFromOptions(cdm, ierr)) PetsCallA(DMGetDimension(cdm, cdim, ierr)) PetsCallA(DMPlexSimplex(cdm, simplex, ierr)) PetsCallA(PetscFECreateDefault(PETSC_COMM_WORLD, cdim, cdim, simplex, "ho_", K, cfe, ierr)) PetsCallA(PetscFEGetDimension(cfe, fedim, ierr)) PetsCallA(PetscViewerCreate(PETSC_COMM_WORLD, viewer, ierr)) PetsCallA(PetscViewerDrawOpen(PETSC_COMM_WORLD, PETSC_NULL_CHARACTER, PETSC_NULL_CHARACTER, 100, 100, 800, 800, viewer, ierr)) PetsCallA(PetscViewerSetFromOptions(viewer, ierr)) ! also using flag -draw_pause 5 PetsCallA(DMView(cdm, viewer, ierr)) [...] With the following list of flags: -ho_petscspace_degree K -ho_petscdualspace_lagrange_node_type equispaced -ho_petscdualspace_lagrange_node_endpoints 1 -dm_plex_high_order_view -options_left Using "-options_left" shows that "there are no unused options"; so "-dm_plex_high_order_view" is used somehow; it is at least required for the call to "DMPlexCreateHighOrderSurrogate_Internal" within the "draw" functions. >From the CoordinateDM I can see the additional nodes created for the higher order approximation (e.g. mid-edge, mid-face nodes), so it seems the FE space is correct. Regardless of the order "K", the mesh plotted with "DMView" is always the same, corresponding to the linear case. Thank you. Noam

2 1

Petsc veccuda device to host copy
by Wenbo Zhao 07 Dec '25

07 Dec '25

Hi, we are using petsc's veccuda and found that the data in the host array obtained via VecGetArrayRead is partially updated sometime. Vec vgpu, vcpu; iterations: // ksp solve a * vgpu=b const PetscScalar * agpu; PetscScalar * acpu; VecGetArrayRead(vgpu, &agpu); VecGetArray(vcpu, &acpu); PetscArraycpy (acpu, agpu,size); // check updating std::cout << agpu[0] << agpu [size-1]<<std::endl; // we found that agpu[0] is last iterations value, agpu[size-1] updated from device value, randomly // use acpu values to update matrix a .... Petsc 3.21.1 is used. And manual said, For vectors that may also have array data in GPU memory, for example, VECCUDA, this call ensures the CPU array has the most recent array values by copying the data from the GPU memory if needed. Best wishes, Wenbo

2 1

TAO PDIPM handling of objective evaluation failures (NaN / PDE non-convergence)
by Simon Wiesheier 05 Dec '25

05 Dec '25

Dear PETSc developers and users, I am considering using TAO’s Primal-Dual Interior-Point Method (PDIPM) for a constrained optimization problem in solid mechanics. The objective involves solving a nonlinear PDE (hyperelasticity) for each parameter vector, and for some parameter combinations the PDE solver may fail to converge or produce non-physical states. With MATLAB’s fmincon, it is possible to signal such failures by returning NaN/Inf for the objective, and the solver will then backtrack or try a different step without crashing. My questions are: 1. How does TAO’s PDIPM handle cases where the user objective or gradient callback returns NaN/Inf (e.g., due to PDE solver failure)? 2. Is there a recommended way in TAO/PETSc to gracefully signal an evaluation failure (like “bad point in parameter space”) so that the algorithm can back off and try a smaller step, instead of aborting? 3. If the recommended pattern is *not* to return NaNs, what is the best practice in TAO for such PDE-constrained problems? Any guidance on how TAO/PDIPM is intended to behave in the presence of evaluation failures would be greatly appreciated. Best regards, Simon

2 2

Questions about memory usage in Petsc
by neil liu 03 Dec '25

03 Dec '25

Dear users and developers, I am recently running a large system from Nedelec element, 14 million dofs (complex number). A little confused about the memory there. Then I tried a small system (34,000 dofs) to see the memory usage. It was solved with MUMPS with 1 rank. Then I used PetscMemoryGetCurrentUsage() to show the memory used there. The pseudocode is PetscMemoryGetCurrentUsage (*Memory 1:* 64.237M) KSPset KSPsolve (INFOG(18) (size of all MUMPS internal data allocated during factorization: value on the most memory consuming processor): *Memory 2:* 408 MB) PetscMemoryGetCurrentUsage (*Memory 3: *54.307M) [0] Maximum memory PetscMalloc()ed 49.45MB maximum size of entire process 424MB (*Memory 4:* 54.307M) The following is my understanding, please correct me if I am wrong, It seems the difference between Memory 1 and 3 is approximately the size of 30 Krylov vectors (complex). It seems Memory 4 is not the summation of Memory 2 and 3; but on the same order of magnitude. It is a little confusing here. Thanks, Xiaodong

2 2

error setting the type of the TAO solver
by Moral Sanchez, Elena 02 Dec '25

02 Dec '25

Hi, I am trying to initialize a LCL TAO solver with petsc4py: from petsc4py import PETSc solver = PETSc.TAO().create() solver.setType(PETSc.TAO.Type.LCL) The last line throws the following error: Traceback (most recent call last): File "<string>", line 3, in <module> File "petsc4py/PETSc/TAO.pyx", line 183, in petsc4py.PETSc.TAO.setType petsc4py.PETSc.Error: error code 86 [0] TaoSetType() at /petsc/src/tao/interface/taosolver.c:2164 [0] Unknown type. Check for miss-spelling or missing package: https://urldefense.us/v3/__https://petsc.org/release/install/install/*exter… [0] Unable to find requested Tao type lcl However, hasattr(solver.Type(), 'LCL') returns True. The same happens with any other PETSc.TAO.Type. What am I missing here? Cheers, Elena

1 0

Bug report: petcds fortran module missing?
by Berends, C.J. (Tijn) 27 Nov '25

27 Nov '25

Dear petsc folks, I am trying to set up a very basic example program in Fortran, using Petsc to solve a simple Poisson problem (d2u/dx2 = f). However, I am running into problems when I try to call PetscDSSetResidual to set the function pointers for the residual. According to the documentation, this function is part of the petcds module, but while the header file (petscds.h) exists in my Petsc include directotry, the Fortran module (petscds.mod) does not. Therefore, my code won't compile. I have currently got Petsc installed via Homebrew. I have tried cloning the Petsc git repository and building that, using the configuration option --with-fortran-interfaces=1, but still the petscds.mod file is not there. Am I doing something wrong, or can it be that this particular Fortran interface is just missing? If you need any additional information from me, please let me know. Kind regards, Tijn Berends dr. C. J. (Tijn) Berends Post-doc (palaeo)glaciology Institute for Marine and Atmospheric research Utrecht (IMAU), Utrecht University, The Netherlands Buys Ballot Building (BBG), Room 6.67

2 1

[WARNING: UNSCANNABLE EXTRACTION FAILED]GMRES plus BlockJacobi behave differently for seemlingy identical matrices
by Nils Schween 26 Nov '25

26 Nov '25

Dear PETSc users, Dear PETSc developers, in our software we are solving a linear system with PETSc using GMRES in conjunction with a BlockJacobi preconditioner, i.e. the default of the KSP object. We have two versions of the system matrix, say A and B. The difference between them is the non-zero pattern. The non-zero pattern of matrix B is a subset of the one of matrix A. Their values should be identical. We solve the linear system, using A yields a solution after some iterations, whereas using B does not converge. I created binary files of the two matrices, the right-hand side, and wrote a small PETSc programm, which loads them and demonstrates the issue. I attach the files to this email. We would like to understand why the solver-preconditioner combination works in case A and not in case B. Can you help us finding this out? To test if the two matrices are identical, I substracted them and computed the Frobenius norm of the result. It is zero. To give you more context, we solve a system of partial differential equations that models astrophysical plasmas. It is essentially a system of advection-reaction equations. We use a discontinuous Galerkin (dG) method. Our code relies on the finite element library library deal.ii and its PETSc interface. The system matrices A and B are the result of the (dG) discretisation. We GMRES with a BlockJaboci preconditioner, because we do not know any better. I tested the code I sent with PETSc 3.24.0 and 3.19.1 on my workstation, i.e. Linux home-desktop 6.17.2-arch1-1 #1 SMP PREEMPT_DYNAMIC Sun, 12 Oct 2025 12:45:18 +0000 x86_64 GNU/Linux I use OpenMPI 5.0.8 and I compiled with mpicc, which in my cases use gcc. In case you need more information. Please let me know. Any help is appreciated. Thank you, Nils -- Nils Schween Phone: +49 6221 516 557 Mail: nils.schween(a)mpi-hd.mpg.de PGP-Key: 4DD3DCC0532EE96DB0C1F8B5368DBFA14CB81849 Max Planck Institute for Nuclear Physics Astrophysical Plasma Theory (APT) Saupfercheckweg 1, D-69117 Heidelberg https://www.mpi-hd.mpg.de/mpi/en/research/scientific-divisions-and-groups/i…

3 5