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

Questions on QR methods in PETSc
by Leonardo De Novellis 27 Mar '26

27 Mar '26

Dear Users Support Team, I have a couple of questions regarding QR decompositions in PETSc. 1) I would like to find the least-squares solution of a rectangular system of the form Ax = b, where A is a dense and tall-skinny matrix (size around 1000000 x 10). Since A has very bad conditioning, I want to avoid iterative methods (such as KSPLSQR), since the number of iterations can get very large, and would like to use a direct QR solving method. Currently, I am running my code on only 1 core, and A is of type seqaij. With this setup, the following code works fine: call KSPSetType(ksp, KSPPREONLY, ierr) call KSPGetPC(ksp, pc, ierr) call PCSetType(pc, PCQR, ierr) call KSPSetOperators(ksp, A, A, ierr) t1 = MPI_Wtime() call KSPSolve(ksp, b, x, ierr) t2 = MPI_Wtime() I eventually want to run this in parallel on multiple cores. Will PCQR work for an mpiaij / mpidense matrix? If not, what would you suggest as a direct solving approach for this system? 2) I would also like to compute an explicit QR decomposition of A, and want to do so in parallel. Is there any way to do so in PETSc? If not, as a possible alternative, would you recommend using SLEPc BVOrthogonalize function? Kind regards, Leonardo

3 2

Different precision from MatAssembly/MatView
by Noam T. 25 Mar '26

25 Mar '26

Hello, Looking at the assembly of the Jacobian using MatView, I noticed values were somewhat different from the computed ones in the jacobian function (DMSNESSetJacobian). For the test, I used a single Q1 element, with 4 quadrature points, so the exact Jacobian matrix can be computed analytically. In the jacobian function, after looping through all quadrature points, printing the matrix to stdout shows: ---- 0.531360000000 0.066670000000 0.197333333327 0.066670000000 0.535360000000 -0.003333333327 0.197333333327 -0.003333333327 0.527413333333 ---- (showing only the 3 free DoFs). The 2x2 submatrix in [0:1] has exactly four/five nonzero digits, as shown above. The rest of the elements have a periodic 3 ending (a couple of digits at the end are off, but that's fine). This same matrix is the one added to the global jacobian during the quadrature loop, with "DMPlexMatSetClosure". After looping, calling "MatAssemblyBegin/End", then MatView: ---- Mat Object: 1 MPI process type: seqaij row 0: (0, 0.531387) (1, 0.0666733) (2, 0.197338) row 1: (0, 0.0666733) (1, 0.535387) (2, -0.0033385) row 2: (0, 0.197338) (1, -0.0033385) (2, 0.52744) ---- Values are close, but definitely not the same as computed. Are these values from "MatView" the ones actually stored in the matrix? It seems there is some precision loss in the process. Interestingly, computing the jacobian through finite differences, shows the same result for the "hand-coded" jacobian, whereas the finite differences ones is the "exact" one (same as first one above): ---- ||J - Jfd||_F/||J||_F = 4.90736e-05, ||J - Jfd||_F = 4.74265e-05 Hand-coded Jacobian ---------- Mat Object: 1 MPI process type: seqaij row 0: (0, 0.531387) (1, 0.0666733) (2, 0.197338) row 1: (0, 0.0666733) (1, 0.535387) (2, -0.0033385) row 2: (0, 0.197338) (1, -0.0033385) (2, 0.52744) Finite difference Jacobian ---------- Mat Object: 1 MPI process type: seqaij row 0: (0, 0.53136) (1, 0.06667) (2, 0.197333) row 1: (0, 0.06667) (1, 0.53536) (2, -0.00333333) row 2: (0, 0.197333) (1, -0.00333333) (2, 0.527413) Hand-coded minus finite-difference Jacobian with tolerance 1e-05 ---------- Mat Object: 1 MPI process type: seqaij row 0: (0, 2.66554e-05) row 1: (1, 2.66554e-05) row 2: (2, 2.6558e-05) ---- The main code is in Fortran, using double precision variables, in case it is relevant. Thank you. Noam

2 11

20-node hexes
by Noam T. 25 Mar '26

25 Mar '26

Hello, Is it possible to build/use a mesh consisting of 20-node hexes in PETSc? I tried reading (-dm_plex_filename) an existing mesh, with a 20-node cube, from gmsh (using the option "SecondOrderIncomplete"); statistics show 20 nodes, however the "Entities" section still only contains 8 vertices, so that may not work here. Using -petscdualspace_lagrange_XYZ I could not achieve this either. I searched some old threads in the mailing list; there are a couple, but they don't specify how they build it; they just mention they are using such mesh. Searching for "serendipity" in the documentation did not return anything either. Thanks,

2 1

Issue: Building PLC from STL for DMPlexGenerate_Tetgen and Geometry-Aware workflows
by 陈日刀 23 Mar '26

23 Mar '26

Dear PETSc Team, My name is Zhaoyang Chen, a senior undergraduate at Hunan University. I am working on a meshing workflow using DMPlex and am particularly interested in utilizing the geometry-aware features such as *DMPlexSnapToGeomModel()* and refinement snapping. My input data consists of *high-resolution STL triangulations*. I have encountered a major bottleneck at the very beginning of the pipeline: *constructing the Piecewise Linear Complex (PLC)* from the STL to be used by *DMPlexGenerate_Tetgen().* Since the downstream geometry-aware operations (like *DMPlexSnapToGeomModel()* ) rely on a valid underlying geometry representation, I need to resolve the PLC construction issue first. Specifically: 1. *STL to PLC Ingestion*: Is there a native or recommended pathway in PETSc to ingest an STL file and define it as the PLC for volume mesh generation? 2. *Compatibility with Snapping*: If the PLC is successfully built from an STL, can the resulting* DMPlex* still utilize *DMPlexSnapToGeomModel()*? Or is the "Geometry-Aware" snapping strictly reserved for analytical B-Rep models (e.g., via EGADS)? I am trying to understand if PETSc supports a workflow where a discrete surface (STL) acts as the "GeomModel" for both initial volumetric meshing and subsequent snapping/refinement. Any guidance on how to bridge STL data with the DMPlex geometry-aware pipeline would be greatly appreciated. Best regards, Zhaoyang Chen Hunan University

1 0

Request for Assistance with Parallel Matrix Solving Using petsc4py in EasyFEA
by Matthieu Noel 21 Mar '26

21 Mar '26

Dear PETSc users, My name is Matthieu Noel, and I am a research software engineer at INRIA. During my PhD thesis, I developed *EasyFEA*, an open-source finite element analysis tool documented here: https://urldefense.us/v3/__https://easyfea.readthedocs.io/en/stable/__;!!G_… . I am currently working on *Issue #26* (https://urldefense.us/v3/__https://github.com/matnoel/EasyFEA/issues/26__;!… ), where I aim to use *petsc4py* to solve a simple stationary thermal problem in parallel. The main challenge is *constructing and solving the parallel matrix system* (|Ax = b|). Current Status: * *Mesh partitioning* has been successfully performed using *GMSH*. * EasyFEA provides the *assembled matrix system* (|Ax = b|) with applied boundary conditions. * My current implementation using *petsc4py* runs without errors, but the results obtained are incorrect. I suspect the problem lies in how the *parallel matrix and vectors are assembled or solved*, particularly in handling *ghost DOFs* and ensuring proper communication between MPI ranks. I would be extremely grateful for any guidance, resources, or examples you could share to help me address this issue. My current development branch is available here: https://urldefense.us/v3/__https://github.com/matnoel/EasyFEA/tree/26_mpi__… . The PETSc function can be found at this link: https://urldefense.us/v3/__https://github.com/matnoel/EasyFEA/blob/26_mpi/E… I am currently running the script in the attached files. Thank you in advance for your time and expertise. Best regards, Matthieu Noel

2 2

Question regarding
by Aldo Bonfiglioli 20 Mar '26

20 Mar '26

Hi there, I am (successfully) using the following piece of code to insert elemental matrices into the "global" matrix "jac": > ! > if( matrix_assembly )then > ! > i = 0 > PetscCall(DMPlexGetTransitiveClosure(dm, cell, useCone, > PETSC_NULL_INTEGER, p_points, ierr)) > ! > ! loop over the vertices of the cell > ! > do j = vtxbgn, vtxend, 2 ! this is ndim-dependent > i = i + 1 > point = p_points(j) > PetscCall(PetscSectionGetDof(section, point, ndofs, ierr)) > PetscCall(PetscSectionGetOffset(section, point, offset, > ierr)) > row(i) = offset/ndofs > col(i) = offset/ndofs > end do > ! > do i = 1, nofvert > do j = 1, nofvert > dummy = transpose(eltmat(:, :, i, j)) > ir = row(i) > ic = col(j) > ! > !p_dummyisa1dviewofdummy > ! > PetscCall(MatSetValuesBlockedLocal(jac, 2, [ir], 1, [ic], p_dummy, > ADD_VALUES, ierr)) > end do > end do > ! > PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, useCone, > PETSC_NULL_INTEGER, p_points, ierr)) > ! > end if ! matrix_assembly > I would however be happy to replace all of the above with a single call to: > PetscCall(DMPlexMatSetClosure(dm, PETSC_NULL_SECTION, > PETSC_NULL_SECTION, jac, cell, values, ADD_VALUES, ierr )) which I actually succeeded to do as long as ndofs = 1. The question is: how (in which order) should I "transfer" the entries of eltmat(1:ndofs,1:ndofs,1:nofvert,1:nofvert) into the values array required by DMPlexMatSetClosure ? 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?m=002423__;!!G_uCfscf7eWS!ZusjOf_nCidORTOchtNsX_xGJZNjiAXMKKmn4rDhnR-3Y6mFzMF9Z7f7iJFeGj6tgcBzT2L5gtfxO56AMvPr3pgXHaCPziR-QL4$

2 1

VecConcatenate vs VecCreateNest for coupled problems
by Miguel Molinos 06 Mar '26

06 Mar '26

Dear all, I am solving a nonlinear system with SNES that comes from the stationary conditions of a scalar functional (L): \frac{\partial L}{\partial F} = 0, \frac{\partial L}{\partial \sigma_q} = 0 where (F) has size 3 and (\sigma_q) contains (N) particle variables, so the global unknown vector has size (3*N_rank + N). And, locally, each rank has a local vector (3+N_local). I am currently using SNES with GMRES to solve the problem since I can not compute the Jacobian analytically. I have a couple of questions: I originally constructed the solution vector as a single contiguous vector (effectively [F | sigma]) with VecConcatenate, but I recently switched to VecCreateNest to represent the block structure. However, the solver seems to perform better with VecConcatenate. Is this expected when using GMRES, if no block preconditioner is used? The system has a known null space associated with the (sigma) variables. Since I am not assembling the Jacobian explicitly, what is the recommended way to provide this null space to the solver? Should it still be attached to the operator with MatSetNullSpace when using -snes_mf? Any advice or best practices would be very helpful. Best regards, Miguel Miguel Molinos Escuela de Ingenieros de Caminos, Canales y Puertos Dep. Mecánica Medios Continuos y Teoría de Estructuras Prof. Aranguren, 3 - Ciudad Universitaria. Planta 9. 28040, Madrid ✉ m.molinos(a)upm.es <mailto:[email protected]> ✆ +34 910 67 44 28

1 0

Preconditioners with block matrix
by Praveen C 03 Mar '26

03 Mar '26

Hello I am using PETSc with deal.II. My matrix is defined as PETScWrappers::BlockSparseMatrix which uses MATNEST. Is it possible to use PCILU (in serial) or PCHYPRE with such a matrix ? thanks praveen

2 2

how is snes_test_jacobian_view computing the matrix difference?
by Matteo Semplice 03 Mar '26

03 Mar '26

Dear all, the following output from -snes_test_jacobian_view lokks clear enough to pinpoint my mistake (I must be putting jacobian elements in the wrong place), but it nevertheless left me suprised. 0 SNES Function norm 2.945516399303e-01 ---------- Testing Jacobian ------------- Testing hand-coded Jacobian, if (for double precision runs) ||J - Jfd||_F/||J||_F is O(1.e-8), the hand-coded Jacobian is probably correct. ||J - Jfd||_F/||J||_F = 0.000121529, ||J - Jfd||_F = 0.395752 Hand-coded Jacobian ---------- Mat Object: 4 MPI processes row 84: (84, 60.) (85, 0.) (3446, -80.) (3447, 0.) (6808, 20.) (6809, 0.) row 85: (84, 0.00914434) (85, 1.) (3446, 0.) (3447, 0.) (6808, 0.) (6809, 0.) Finite difference Jacobian ---------- Mat Object: 4 MPI processes type: mpiaij row 84: (84, 60.) (1806, -80.) (3528, 20.) row 85: (85, 0.999971) Hand-coded minus finite-difference Jacobian with tolerance 1e-05 ---------- Mat Object: 4 MPI processes type: mpiaij row 84: row 85: (84, 0.00914434) (85, 2.9293e-05) The difference of the two matrices shoudl be nonzero in row 84 for all elements 1806, 3528, 3446, 6808. Is it intended that the difference is shown only for matrix elements that are present in both the hand-coded and the FD jacobian? I would have expected an output like row 84: (1806,80.) (3446,-80.), etc (This was on PETSc Release Version 3.24.2) Thanks Matteo -- Prof. Matteo Semplice Università degli Studi dell’Insubria Dipartimento di Scienza e Alta Tecnologia – DiSAT Professore Associato Via Valleggio, 11 – 22100 Como (CO) – Italia tel.: +39 031 2386316

3 6

Fw: DG methods in PETSc
by Matteo Leone 20 Feb '26

20 Feb '26

Inviato da Outlook per Android<https://urldefense.us/v3/__https://aka.ms/AAb9ysg__;!!G_uCfscf7eWS!esCSra9p… > ________________________________ From: Matteo Leone <matteo4.leone(a)mail.polimi.it> Sent: Thursday, February 19, 2026 3:36:44 PM To: Mark Adams <mfadams(a)lbl.gov> Subject: R: [petsc-users] DG methods in PETSc Hello, I now have tested quite a bit the code, and I have some doubts. Note I am quite new to PETSc. I have modified a bit your code to test it and I tried to print from the functions to check if they were actually called with negative results. The Riemann solver is never called and also the solution is just static in time. I share the code. It's a bit modified and no more c but c++ (I do not know to write the csv file fom c, in case is a dealbreaker I'll try to give a look on how to do it). (For reference I use PETSc 3.24.4 by means of nix, the nix flake is in the shared docs if you are used to it lmk). There is also a small .py to handle the visualization (It was made by Claude code, but it just plots the results) Hopefully we can go deep into this PETSc DG stuff. Thanks in advace. Matteo ________________________________ Da: Mark Adams <mfadams(a)lbl.gov> Inviato: mercoledì 11 febbraio 2026 17:25 A: Matteo Leone <matteo4.leone(a)mail.polimi.it>; PETSc users list <petsc-users(a)mcs.anl.gov> Oggetto: Re: [petsc-users] DG methods in PETSc Great, and keep it on the list. Lots of people here to help! On Wed, Feb 11, 2026 at 9:46 AM Matteo Leone <matteo4.leone(a)mail.polimi.it<mailto:[email protected]>> wrote: Wow thank you so much! I was almost hopeless, I'll deep dive into it and I'll give you a feedback. Matteo ________________________________ From: Mark Adams <mfadams(a)lbl.gov<mailto:[email protected]>> Sent: Wednesday, February 11, 2026 3:38:30 PM To: Matteo Leone <matteo4.leone(a)mail.polimi.it<mailto:[email protected]>> Cc: petsc-users(a)mcs.anl.gov<mailto:[email protected]> <petsc-users(a)mcs.anl.gov<mailto:[email protected]>> Subject: Re: [petsc-users] DG methods in PETSc DG (discontinuous Galerkin) is done with a "Broken" FE. Yikes I do not see a test. Here is a test but it is not well verified. Mark On Tue, Feb 10, 2026 at 10:25 AM Matteo Leone via petsc-users <petsc-users(a)mcs.anl.gov<mailto:[email protected]>> wrote: Hello, I already posted on Reddit but just to be sure I write even here. First thanks for the job you do for PETSc, I have used it for several projects and is always nice. I am writing cause I am getting mad trying to implement DG solver in PETSc, the target is the Euler equations, however I am failing even with just the simplest transport equation (u/t + u/x = 0). I was wondering if I am missing somenthing. I tried with the DSSetReimannSolver and DualSpaces, and more, but I keep failing, I tried also with LLMs, but seems like there is no DG code with PETSc on the web, however I see many papers that do it. I was wondering if I am maybe missing something out or what. For reference I use PETSc 3.24.3 by means of nix. Thanks in advance, cheers. Matteo

2 2