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

Nodes coordinates in distributed dmplex
by Alejandro D Otero 16 Jul '22

16 Jul '22

Hi, sorry if this is an obvious question, but I cannot figure out how to recover finite element nodes coordinates once I have distributed a mesh stored as a dmplex. I am using petsc4py as interface to petsc rutines. I first created a dmplex using: dm.createFromCellList() In a sequential run I got the coordinates with: Coords = dm.getCoordinates() which gave a sequential vector with the coordinates of the mesh nodes. When I distribute the mesh with: dm.distribute() each mpi process has it own dm but the indexing of the vector resulting from getCoordinates() or getCoordinatesLocal() seems not consistent with the local numbering of the cells and nodes. Which is the correct way of doing this in PETSc philosophy? Thanks in advance, Alejandro

4 9

MatGetDiagonal
by Nachiket Gokhale 09 Dec '19

09 Dec '19

Any chance of making this work in serial? http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatGetDiago… Not a show stopper, I am trying to get the diagonal of some small projected, dense matrices (which come from large sparse matrices). I am running in serial because 1) Since my projected matrices are small, and 2) PETSc does not do certain matrix multiplications involving a dense matrix in parallel, -Nachiket

4 5

GAMG scaling
by Mark Adams 24 Dec '18

24 Dec '18

(Hong), what is the current state of optimizing RAP for scaling? Nate, is driving 3D elasticity problems at scaling with GAMG and we are working out performance problems. They are hitting problems at ~1.5B dof problems on a basic Cray (XC30 I think). Thanks, Mark

7 28

monitoring the convergence of fieldsplit 0 and 1
by Klaij, Christiaan 15 May '18

15 May '18

A question when using PCFieldSplit: for each linear iteration of the system, how many iterations for fielsplit 0 and 1? One way to find out is to run with -ksp_monitor, -fieldsplit_0_ksp_monitor and -fieldsplit_0_ksp_monitor. This gives the complete convergence history. Another way, suggested by Matt, is to use -ksp_monitor, -fieldsplit_0_ksp_converged_reason and -fieldsplit_1_ksp_converged_reason. This gives only the totals for fieldsplit 0 and 1 (but without saying for which one). Both ways require to somehow process the output, which is a bit inconvenient. Could KSPGetResidualHistory perhaps return (some) information on the subsystems' convergence for processing inside the code? Chris dr. ir. Christiaan Klaij CFD Researcher Research & Development E mailto:[email protected] T +31 317 49 33 44 MARIN 2, Haagsteeg, P.O. Box 28, 6700 AA Wageningen, The Netherlands T +31 317 49 39 11, F +31 317 49 32 45, I www.marin.nl

4 19

Question about memory usage in Multigrid preconditioner
by frank 22 Jan '18

22 Jan '18

Hi, I am using the CG ksp solver and Multigrid preconditioner to solve a linear system in parallel. I chose to use the 'Telescope' as the preconditioner on the coarse mesh for its good performance. The petsc options file is attached. The domain is a 3d box. It works well when the grid is 1536*128*384 and the process mesh is 96*8*24. When I double the size of grid and keep the same process mesh and petsc options, I get an "out of memory" error from the super-cluster I am using. Each process has access to at least 8G memory, which should be more than enough for my application. I am sure that all the other parts of my code( except the linear solver ) do not use much memory. So I doubt if there is something wrong with the linear solver. The error occurs before the linear system is completely solved so I don't have the info from ksp view. I am not able to re-produce the error with a smaller problem either. In addition, I tried to use the block jacobi as the preconditioner with the same grid and same decomposition. The linear solver runs extremely slow but there is no memory error. How can I diagnose what exactly cause the error? Thank you so much. Frank

7 56

Re: [petsc-users] A number of questions about DMDA with SNES and Quasi-Newton methods
by zakaryah . 03 Nov '17

03 Nov '17

OK - I've checked the Jacobian and function very thoroughly and I am confident there are no errors. I suspect that I am having problems with a bad starting point, and the SNES cannot find the global minimum from there. I know that the global minimum (with residual zero) exists in all cases but I would like the methods for finding it to be more robust to the starting value. The problem comes from the physics of finite deformations of elastic materials. In short, I have a functional of smooth maps on a 3D domain to itself. The functional contains two terms. The first term represents forces which come from external data, and in the Lagrangian this term only involves the values of the map at the point in question. The second term penalizes fluctuations in the map, and can take various forms. The simplest form is just the Dirichlet energy, but I'm also interested in the infinitesimal strain energy and the finite strain energy. The first two have terms in the Lagrangian which are first order in the second spatial derivatives of the map, while the third (finite strain energy) has terms which are up to third order in the first and second spatial derivatives of the map. It is the finite strain energy term which has been problematic. The Euler-Lagrange equations are discretized on a cubic grid, with equal interval spacing in each dimension. The map is the dependent variable, i.e. the x in F(x) = 0. I prefer Neumann boundary conditions. Because the spatial derivatives of the map are usually small, the Jacobian typically has large values in 3x3 blocks along the diagonal (which depend on the map and the external data), and up to 54 values which are functions of the spatial derivatives of the map and tend to be smaller. Do you have any advice on diagnosing and improving situations in which Newton's method finds a stationary point that is not the state with globally minimal residual? One hint is that -snes_type newtonls does not find as good a solution as -snes_type newtontr but I don't know much about these trust region methods, or how to customize and assess them. I'm grateful for any advice. On Mon, Sep 4, 2017 at 5:44 PM, zakaryah . <zakaryah(a)gmail.com> wrote: > Yes, it looks like it IS the other way around, and I think the row is > > r.c + r.i*3 + r.j*3*M + r.k*3*M*N, where r.i is in [0,M-1], r.j is in > [0,N-1], and r.k is in [0,P-1]. > > That matches the boundary conditions in the displayed Jacobian. > > On Mon, Sep 4, 2017 at 5:33 PM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > >> >> > On Sep 4, 2017, at 4:09 PM, zakaryah . <zakaryah(a)gmail.com> wrote: >> > >> > OK that is super helpful. Just to be sure - for MxNxP, the row r in >> the Jacobian is at r.i*P*N*3 + r.j*P*3 + r.k*3 + r.c? >> >> It is that way, or the other way around r.k*M*N*3 + r.j*N*3 + r.k*3 + >> r.c >> > >> > >> > On Mon, Sep 4, 2017 at 4:58 PM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: >> > >> > > On Sep 4, 2017, at 3:48 PM, zakaryah . <zakaryah(a)gmail.com> wrote: >> > > >> > > One piece of information that would be useful is what ordering PETSc >> uses for the Jacobian in the snes_test_display. Is it a natural ordering, >> or the PETSc ordering? For debugging the Jacobian manually, the natural >> ordering is much easier to work with. >> > >> > What is displayed is always the natural ordering (internally it is >> not the natural ordering). >> > >> > > For -n 1, are the orderings the same? >> > >> > yes >> > >> > >> > >> > > >> > > If I use a MatStencil r to represent a field with 3 degrees of >> freedom, and the dimensions of my 3D DMDA are MxNxP, which row of the >> Jacobian corresponds to r.i=x, r.j=y, r.k=z, r.c=f? >> > >> > Internally it is complicated but for any viewing it is just the natural >> ordering and all the degrees of freedom for a single point are next to each >> other in the vector/matrix. >> > >> > >> >> >

6 47

preconditioning matrix-free newton-krylov
by Mark W. Lohry 30 Oct '17

30 Oct '17

I'm currently using JFNK in an application where I don't have a hand-coded jacobian, and it's working well enough but as expected the scaling isn't great. What is the general process for using PC with MatMFFDComputeJacobian? Does it make sense to occasionally have petsc re-compute the jacobian via finite differences, and then recompute the preconditioner? Any that just need the sparsity structure? Are there any PCs that don't work in the matrix-free context? Are there any example codes I overlooked? Last but not least... can the Boomer-AMG preconditioner work with JFNK? To really show my ignorance of AMG, can it actually be written as a matrix P^-1(Ax-b)=0, , or is it just a linear operator? Thanks, Mark

4 4

slepc NHEP error
by Kannan, Ramakrishnan 23 Oct '17

23 Oct '17

Hello, I am running NHEP across 16 MPI processors over 16 nodes in a matrix of global size of 1,000,000x1,000,000 with approximately global 16,000,000 non-zeros. Each node has approximately 1million non-zeros. The following is my slepc code for EPS. PetscInt nev; ierr = EPSCreate(PETSC_COMM_WORLD, &eps); // CHKERRQ(ierr); ierr = EPSSetOperators(eps, A, NULL); // CHKERRQ(ierr); ierr = EPSSetProblemType(eps, EPS_NHEP); // CHKERRQ(ierr); EPSSetWhichEigenpairs(eps, EPS_LARGEST_REAL); EPSSetDimensions(eps, 100, PETSC_DEFAULT, PETSC_DEFAULT); PRINTROOT("calling epssolve"); ierr = EPSSolve(eps); // CHKERRQ(ierr); ierr = EPSGetType(eps, &type); // CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD, " Solution method: %s\n\n", type); // CHKERRQ(ierr); ierr = EPSGetDimensions(eps, &nev, NULL, NULL); // CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD, " Number of requested eigenvalues: %D\n", nev); // CHKERRQ(ierr); I am getting the following error. Attached is the entire error file for your reference. Please let me know what should I fix in this code. 2]PETSC ERROR: Argument out of range [2]PETSC ERROR: Argument 2 out of range [2]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting. [2]PETSC ERROR: Petsc Release Version 3.7.4, Oct, 02, 2016 [2]PETSC ERROR: ./miniapps on a sandybridge named nid00300 by d3s Wed Jun 14 15:32:00 2017 [13]PETSC ERROR: #1 BVDotVec_BLAS_Private() line 272 in /autofs/nccs-svm1_home1/ramki/libraries/slepc-3.7.3/src/sys/classes/bv/interface/bvblas.c [13]PETSC ERROR: #2 BVDotVec_Svec() line 150 in /autofs/nccs-svm1_home1/ramki/libraries/slepc-3.7.3/src/sys/classes/bv/impls/svec/svec.c [13]PETSC ERROR: #3 BVDotVec() line 191 in /autofs/nccs-svm1_home1/ramki/libraries/slepc-3.7.3/src/sys/classes/bv/interface/bvglobal.c [13]PETSC ERROR: #4 BVOrthogonalizeCGS1() line 81 in /autofs/nccs-svm1_home1/ramki/libraries/slepc-3.7.3/src/sys/classes/bv/interface/bvorthog.c [13]PETSC ERROR: #5 BVOrthogonalizeCGS() line 214 in /autofs/nccs-svm1_home1/ramki/libraries/slepc-3.7.3/src/sys/classes/bv/interface/bvorthog.c [13]PETSC ERROR: #6 BVOrthogonalizeColumn() line 371 in /autofs/nccs-svm1_home1/ramki/libraries/slepc-3.7.3/src/sys/classes/bv/interface/bvorthog.c [13]PETSC ERROR: #7 EPSBasicArnoldi() line 59 in /autofs/nccs-svm1_home1/ramki/libraries/slepc-3.7.3/src/eps/impls/krylov/epskrylov.c [13]PETSC ERROR: #8 EPSSolve_KrylovSchur_Default() line 203 in /autofs/nccs-svm1_home1/ramki/libraries/slepc-3.7.3/src/eps/impls/krylov/krylovschur/krylovschur.c [13]PETSC ERROR: #9 EPSSolve() line 101 in /autofs/nccs-svm1_home1/ramki/libraries/slepc-3.7.3/src/eps/interface/epssolve.c -- Regards, Ramki

3 33

Good recommendation on meshing package?
by Zou, Ling 16 Oct '17

16 Oct '17

Hi all, I know this is a bit off topic on PETSc email list. I would like to try some finite volume type of CFD algorithm with PETSc, but I found it quite troublesome to manage mesh by myself. I wonder if there is any good existing meshing package that works well with PTESc. My expectation on such a package would be: 1) I create the mesh with some tool. 2) Read this mesh with the meshing package, so I have things like node set, edge set, cell set, etc. to play with 3) discretize my PDE with the mesh 4) solve it I also understand many people here use PETSc solve their CFD problem. I would appreciate it if you could also point me to some good examples. Best, Ling

2 2

Using factored complex matrices from MUMPS as a preconditioner in PETSC
by Evan Um 11 Oct '17

11 Oct '17

Dear PETSC Users, My system matrix comes from finite element modeling and is complex and unstructured. Its typical size is a few millions-by a few millions. I wondering if I can use MUMPS parallel direct solver as a preconditioner in PETSC. For example, I want to pass factored matrices to Krylov iterative solvers such as QMR. Is there any PETSC+MUMPS example code for the purpose? Can PETSC call the latest MUMPS that supports block low rank approximation? In advance, thank you very much for your comments. Best, Evan

4 15