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

Mat indices for DMPlex jacobian
by Matteo Semplice 25 Sep '24

25 Sep '24

Dear petsc, I need to hand-code a jacobian and I can't figure out how to translate the DMplex points/fields to matrix indices. My DMPlex has a section with m fields per cell (which makes for n>m dof per cell since some fields are vector). Say I want to insert an nxn block for the row corresponding to cell c and coloumn of its neighbour d. I guess that I should call either MatSetValues/MatSetValuesLocal or the blocked variants, but how do I find the row/col indices to pass in starting from the c/d dmplex points? And, while I am at it, which MatSetValues version standard/Local standard/blocked is best? I looked for a petsc example, but fails to find it: if there's one, can you just point me to it? Thanks in advance. Matteo Semplice

2 2

Inquiry about reading an exodus II file from coreform Cubit
by neil liu 22 Sep '24

22 Sep '24

Dear Petsc developers and users, I am trying to read an exodus II file from coreform cubit, but without success, Then I used petsc's built-in exodus file, /share/petsc/datafiles/meshes/sevenside.exo. This file can be read by petsc successfully. And I did a test. This file, sevenside.exo, is imported into coreform Cubit, and then is saved as a new exodus file. This new exodus file can not be read by petsc successfully. [0]PETSC ERROR: Invalid argument [0]PETSC ERROR: Unrecognized element type SHELL [0]PETSC ERROR: See https://urldefense.us/v3/__https://petsc.org/release/faq/__;!!G_uCfscf7eWS!… for trouble shooting. [0]PETSC ERROR: Petsc Release Version 3.21.1, Apr 26, 2024 [0]PETSC ERROR: ./app on a arch-linux-c-opt named localhost.localdomain by xiaodongliu Fri Sep 20 09:36:16 2024 [0]PETSC ERROR: Configure options -download-mumps -download-scalapack --with-cc=gcc --with-fc=gfortran --with-cxx=g++ --download-fblaslapack --download-mpich --with-scalar-type=complex --download-triangle --with-debugging=no --download-parmetis --download-metis -download-ptscotch --download-bison --download-hdf5 -download-mmg=/home/xiaodongliu/Documents/3rdPartyLib/d5c43d1bcefe598d51428f6a7fee10ec29478b79.tar.gz --download-ctetgen --download-pragmatic --download-eigen --download-netcdf=/home/xiaodongliu/Documents/3rdPartyLib/netcdf-c-4.9.2-p1.tar.gz --download-zlib --download-pnetcdf --download-exodusii [0]PETSC ERROR: #1 ExodusGetCellType_Internal() at /home/xiaodongliu/Documents/petsc-with-docs-3.21.1/petsc-3.21.1/src/dm/impls/plex/plexexodusii.c:1470 [0]PETSC ERROR: #2 DMPlexCreateExodus() at /home/xiaodongliu/Documents/petsc-with-docs-3.21.1/petsc-3.21.1/src/dm/impls/plex/plexexodusii.c:1551 [0]PETSC ERROR: #3 DMPlexCreateExodusFromFile() at /home/xiaodongliu/Documents/petsc-with-docs-3.21.1/petsc-3.21.1/src/dm/impls/plex/plexexodusii.c:1390 Thanks a lot. Xiaodong Compose: Inquiry about reading an exodus II file from coreform Cubit [image: Minimize][image: Pop-out][image: Close] PETSc users list

4 6

Inquiry about reading exodus file for Petsc.
by neil liu 20 Sep '24

20 Sep '24

Dear Petsc developers and users, I am trying to import an exodus II file from coreform Cubit. It is a little

1 0

Question on PCGALERKIN from command line
by Edoardo alinovi 19 Sep '24

19 Sep '24

Hi Barry, I have a quick question on PC Galerking as I am trying to do the same thing reported in this post: https://urldefense.us/v3/__https://lists.mcs.anl.gov/pipermail/petsc-users/… This my command line: mpirun -np 4 flubio_coupled* -UPeqn_pc_type composite -UPeqn_pc_composite_pcs galerkin,bjacobi -UPeqn_sub_0_galerkin_ksp_type preonly -UPeqn_sub_0_galerkin_pc_type none* And this is the error: [0]PETSC ERROR: --------------------- Error Message -------------------------------------------------------------- [0]PETSC ERROR: Object is in wrong state [0]PETSC ERROR: Must set operator of PCGALERKIN KSP with PCGalerkinGetKSP()/KSPSetOperators() Do I need to set something for PCGALERKING in my code to make it work from command line? Thank you!

2 1

Error: PCHYPRESetPoissonMatrix_HYPRE
by Karthikeyan Chockalingam - STFC UKRI 19 Sep '24

19 Sep '24

Hello, I would like to make the following hypre call HYPRE_AMSSetBetaPoissonMatrix(ams, NULL); So it does look like ams_beta_is_zero has to be true https://urldefense.us/v3/__https://gitlab.com/petsc/petsc/-/blob/main/src/k… So would I end up calling PCHYPRESetPoissonMatrix_HYPRE with isalpha being false. https://urldefense.us/v3/__https://gitlab.com/petsc/petsc/-/blob/main/src/k… But I get the following error from Libmesh: is it because of version incompatibility between PETSc and Hypre? vector_fe_ex3.C:141:3: error: 'PCHYPRESetPoissonMatrix_HYPRE' was not declared in this scope; did you mean 'PCHYPRESetBetaPoissonMatrix'? 141 | PCHYPRESetPoissonMatrix_HYPRE(pc, A, false); | ^~~~~~~~~~~~~~~~~~~~~~~~~~~~ I look forward to your response. Kind regards, Karthik -- Karthik Chockalingam, Ph.D. Senior Research Software Engineer High Performance Systems Engineering Group Hartree Centre | Science and Technology Facilities Council karthikeyan.chockalingam(a)stfc.ac.uk<mailto:[email protected]> [signature_3970890138]

2 1

[MPI GPU Aware] KSP_DIVERGED
by LEDAC Pierre 18 Sep '24

18 Sep '24

Hi all, We are using PETSc 3.20 in our code and running succesfully several solvers on Nvidia GPU with OpenMPI library which are not GPU aware (so I need to add the flag -use_gpu_aware_mpi 0). But now, when using OpenMPI GPU Aware library (OpenMPI 4.0.5 ou 4.1.5 from NVHPC), some parallel calculations failed with KSP_DIVERGED_ITS or KSP_DIVERGED_DTOL with several configurations. It may run wells on a small test case with (matrix is symmetric): -ksp_type cg -pc_type gamg -pc_gamg_type classical But suddenly with a number of devices for instance bigger than 4 or 8, it may fail. If I switch to another solver (BiCGstab), it may converge: -ksp_type bcgs -pc_type gamg -pc_gamg_type classical The more sensitive cases where it diverges are the following: -ksp_type cg -pc_type hypre -pc_hypre_type boomeramg -ksp_type cg -pc_type gamg -pc_gamg_type classical And the bcgs turnaroud doesn't work each time... It seems to work without problem with aggregation (at least 128 GPUs on my simulation): -ksp_type cg -pc_type gamg -pc_gamg_type agg So I guess there is a weird thing happening in my code during the solve in PETSc with MPI GPU Aware, as all the previous configurations works with non GPU aware MPI. Here is the -ksp_view log during one fail with the first configuration: KSP Object: () 8 MPI processes type: cg maximum iterations=10000, nonzero initial guess tolerances: relative=0., absolute=0.0001, divergence=10000. left preconditioning using UNPRECONDITIONED norm type for convergence test PC Object: () 8 MPI processes type: hypre HYPRE BoomerAMG preconditioning Cycle type V Maximum number of levels 25 Maximum number of iterations PER hypre call 1 Convergence tolerance PER hypre call 0. Threshold for strong coupling 0.7 Interpolation truncation factor 0. Interpolation: max elements per row 0 Number of levels of aggressive coarsening 0 Number of paths for aggressive coarsening 1 Maximum row sums 0.9 Sweeps down 1 Sweeps up 1 Sweeps on coarse 1 Relax down l1scaled-Jacobi Relax up l1scaled-Jacobi Relax on coarse Gaussian-elimination Relax weight (all) 1. Outer relax weight (all) 1. Maximum size of coarsest grid 9 Minimum size of coarsest grid 1 Not using CF-relaxation Not using more complex smoothers. Measure type local Coarsen type PMIS Interpolation type ext+i SpGEMM type cusparse linear system matrix = precond matrix: Mat Object: () 8 MPI processes type: mpiaijcusparse rows=64000, cols=64000 total: nonzeros=311040, allocated nonzeros=311040 total number of mallocs used during MatSetValues calls=0 not using I-node (on process 0) routines I didn't succeed for the moment creating a reproducer with ex.c examples... Did you see this kind of behaviour before? Should I update my PETSc version ? Thanks for any advice, Pierre LEDAC Commissariat à l’énergie atomique et aux énergies alternatives Centre de SACLAY DES/ISAS/DM2S/SGLS/LCAN Bâtiment 451 – point courrier n°43 F-91191 Gif-sur-Yvette +33 1 69 08 04 03 +33 6 83 42 05 79

3 9

Re: [petsc-users] How to create an array of vectors of type Vec (for use with slepc EPSGetInvariantSubspace)
by dazza simplythebest 18 Sep '24

18 Sep '24

Dear Jose, Many thanks for the guidance. I haven't yet extensively tested the program, but I will cautiously claim that it seems to be working on the small test problem described below (i.e. a small modification of your ex1f.F program). I am currently using version 3.20, but I will upgrade to 3.22 when it comes out. For the Laplace problem with n=50 , using the gd method and seeking to find the two leading eigenmodes (nev=2) using a space of dimension ncv = 30, with default tolerance, on the first solution (using a randomized initial subspace) the algorithm requires 91 iterations to converge, whereas for the second solution (using the previous solution invariant subspace as an initial subspace) it finds the solution in 2 iterations. Many thanks once again and best wishes, Dan. >The Fortran interface is not fully tested. In particular, this function is not employed in any of the Fortran examples. >The Fortran interface has undergone successive improvements in recent times. I have tried the following: > Vec, pointer :: subspace_receive(:) > ... > PetscCall(EPSGetConverged(eps, nconv, ierr)) > PetscCall(MatCreateVecs(A, PETSC_NULL_VEC, xr, ierr)) > PetscCall(VecDuplicateVecsF90(xr, nconv, subspace_receive, ierr)) > PetscCall(EPSGetInvariantSubspace(eps, subspace_receive, ierr)) > ... > PetscCall(VecDestroyVecsF90(nconv, subspace_receive, ierr)) > PetscCall(VecDestroy(xr, ierr)) >I have tested this in the 'main' branch and it works as expected (note that 'main' will become version 3.22 in two weeks). You can try with version 3.21 and it may work, > probably not in older versions. * Jose > El 16 sept 2024, a las 11:07, dazza simplythebest <sayosale(a)hotmail.com> escribió: > > Dear All and Jose, > I am now trying to code up the extraction of the invariant subspace, but I seem to be running into > a compilation problem. I am using the slepc example file ex1f.F as a base, and have added the following lines: > > !Declaration > Vec subspace_receive > . > . > . > ! Executed code (following EPSSolve) > call EPSGetConverged(eps,nconv,ierr); CHKERRA(ierr) > call VecDuplicateVecs(xr, nconv, subspace_receive, ierr); CHKERRA(ierr) > call EPSGetInvariantSubspace(eps, subspace_receive, ierr) ; CHKERRA(ierr). > . > . > > However, I get the following compilation error when I include the call EPSGetInvariantSubspace(eps, subspace_receive, ierr) line: > > > minitest.F(133): error #8284: If the actual argument is scalar, the dummy argument shall be scalar unless the actual argument is of type character or is an element of an array that is not assumed shape, pointer, or polymorphic. [B] > call EPSGetInvariantSubspace(eps, subspace_receive, ierr) ; if (ierr .ne. 0) then;call PetscErrorF(ierr,133,"minitest.F");call MPIU_Abort(PETSC_COMM_SELF,ierr);endif > -----------^ > compilation aborted for minitest.F (code 1) > > > If I comment out the call EPSGetInvariantSubspace(eps, subspace_receive, ierr) line the code compiles, links and runs successfully. > I presume that I must be making a silly mistake in the above lines, but I can't quite see what it is. (There is some discussion in the > manual of there being a different situation if the matrix is real but the eigenvalues are complex, but for this case the eigenvalues > are real). Can you offer any guidance ?! I can provide the whole code of course if that would help. > Many thanks and best wishes, > Dan. > > > > > > From: Jose E. Roman <jroman(a)dsic.upv.es> > Sent: Monday, September 2, 2024 7:06 AM > To: dazza simplythebest <sayosale(a)hotmail.com> > Cc: petsc-users(a)mcs.anl.gov <petsc-users(a)mcs.anl.gov> > Subject: Re: [petsc-users] How to create an array of vectors of type Vec (for use with slepc EPSGetInvariantSubspace) > You can use VecDuplicateVecs() to create an array of Vec. > Yes, you can use VecLoad() to read the data from a binary file. > > Jose > > > > El 2 sept 2024, a las 8:56, dazza simplythebest <sayosale(a)hotmail.com> escribió: > > > > Dear All, > > I am seeking to perform a large number of eigenvalue calculations (related to a fluid dynamics > > problem) using the slepc framework. I have written the program that calls slepc and validated its solutions against > > known results for special cases - the Jacobi-Davidson iteration seemed to perform best on this problem. > > Since in the 'production-level' calculation I will step through a parameter space which should > > mean only small changes to the eigenproblem between each calculation, it seems to make sense > > to use the previous solution as the initial guess for the next calculation. > > > > The sequence to execute this strategy would seem to be: > > To write out the current soln subspace: > > > > call EPSGetInvariantSubspace > > [ call VecView -> binary file (if a file storage is required for possible restart)] > > > > To then read in this data to use it as the initial subspace for the next run: > > > > call EPSSetInitialSpace > > > > My question is that EPSGetInvariantSubspace needs an 'array of nconv vectors', how does one create > > such an array ? My understanding is that the Vec structure is just a one-dimensional array (with contents > > possibly scattered in non-contiguous memory locations !?) , so do we just need to create a long Vec > > that stores multiple vectors ? How can we do this (the examples I have seen thus far use a MatCreateVecs to specify the size > > of the Vec)? > > > > If loading such a subspace from a binary file (created maybe by VecView), presumably VecLoad is the > > best way to load such an array of vectors from the binary file ? > > > > Best wishes and many thanks, > > Dan.

1 0

Residual at the fieldsplit level
by miguel.salazar 17 Sep '24

17 Sep '24

Hello, I have this simple example on Firedrake to illustrate my point. I am solving for a two-component poisson equation (uncoupled). Only the first component has a non-zero residual. ``` import firedrake as fd from firedrake import inner, grad, dx, sin, pi N = 10 mesh = fd.UnitSquareMesh(N, N) V = fd.FunctionSpace(mesh, "CG", 1) W = V * V u = fd.Function(W) v = fd.TestFunction(W) a = inner(grad(u[0]), grad(v[0])) * dx + inner(grad(u[1]), grad(v[1])) * dx x = fd.SpatialCoordinate(mesh) f = fd.Function(V) f.interpolate(fd.Constant(1e4) * sin(x[0] * pi) * sin(2 * x[1] * pi)) L = f * v[0] * dx F = a - L bcs = [fd.DirichletBC(W.sub(0), fd.Constant(2.0), (1,))] def snes_firedrake_residual(F, u, bcs): for bcs_ in bcs: bcs_.apply(u) residual = fd.assemble(F, bcs=bcs, zero_bc_nodes=True) with residual.dat.vec_ro as r: print("Initial residual:", r.norm()) snes_firedrake_residual(F, u, bcs) problem = fd.NonlinearVariationalProblem(F, u, bcs=bcs) solver_mumps_assembled = { "ksp_type": "preonly", "ksp_monitor": None, "pc_type": "python", "pc_python_type": "firedrake.AssembledPC", "assembled_pc_type": "lu", "assembled_pc_factor_mat_solver_type": "mumps", "assembled_mat_mumps_icntl_14": 200, "assembled_mat_mumps_icntl_24": 1, } solver_fieldsplit = { "mat_type": "matfree", "snes_type": "newtonls", "ksp_type": "fgmres", "ksp_rtol": 1e-1, "ksp_monitor": None, "pc_type": "fieldsplit", "pc_fieldsplit_type": "additive", "fieldsplit_0": solver_mumps_assembled, "fieldsplit_1": solver_mumps_assembled, } solver = fd.NonlinearVariationalSolver(problem, solver_parameters=solver_fieldsplit) solver.solve() ``` The PETSc output is as follow ``` Initial residual: 462.13689530272404 0 SNES Function norm 4.621368953027e+02 Residual norms for firedrake_0_ solve. 0 KSP Residual norm 4.621368953027e+02 Residual norms for firedrake_0_fieldsplit_0_ solve. 0 KSP Residual norm 1.000000000000e+00 1 KSP Residual norm 3.501082228626e-15 Residual norms for firedrake_0_fieldsplit_1_ solve. 0 KSP Residual norm 0.000000000000e+00 1 KSP Residual norm 0.000000000000e+00 1 KSP Residual norm 1.612167203819e-12 1 SNES Function norm 1.599350481360e-12 ``` Using the fieldsplit additive preconditioner, the problem converges in a single KSP iteration, as expected. However, I do not understand why the residual of fieldsplit_0 (1e+0) does not coincide with the outer residual (462.13689530272404). It should be the case given that only fieldsplit_0 has a non-zero residual contribution. The fact that it is just 1 is suspicious. Is there something about how the fieldsplit works that I am missing? Thanks, Miguel MIGUEL ANGEL SALAZAR DE TROYA Head of Software Engineering miguel.salazar(a)corintis.com<mailto:[email protected]> Corintis SA EPFL Innovation Park Building C 1015 Lausanne [cid:2024-08-1609_02_11-re_mailmigrationfromgoogletooffice365-sebastien.gobel@corintis.com-co_4456a1bf-dd3e-46b7-85f9-21f278e66a79.png] Here at Corintis we care for your privacy. That is why we have taken appropriate measures to ensure that the data you have provided to us is always secure.

2 9

Inquiry about metric-based adaptation using MMG in petsc
by neil liu 16 Sep '24

16 Sep '24

Dear Petsc developers and users, I am trying to explore adaptive mesh refinement (tetradedra) with Petsc. It seems case 12 ($PETSC DIR/src/snes/tutorials/ex12.c, from paper, https://urldefense.us/v3/__https://arxiv.org/pdf/2201.02806__;!!G_uCfscf7eW… ) is a good example. However when I tried, ./ex12 -run_type full -dm_plex_dim 3 -dm_distribute -dm_plex_box_faces 10,10,10 -bc_type dirichlet -petscspace_degree 1 -variable_coefficient ball -snes_converged _reason ::ascii_info_detail -ksp_type cg -pc_type sor -snes_adapt_sequence 3 -adaptor_target_num 1000 -dm_plex_metric_h_max 0.5 -dm_adaptor mmg L_2 Error: 1.55486 Nonlinear solve converged due to CONVERGED_FNORM_RELATIVE iterations 2 it shows the following error, [0]PETSC ERROR: --------------------- Error Message -------------------------------------------------------------- [0]PETSC ERROR: Argument out of range [0]PETSC ERROR: Global metric normalization factor must be in (0, inf). Is the input metric positive-definite? Do you have any suggestions here? Thanks, Xiaodong

3 4

Re: [petsc-users] How to create an array of vectors of type Vec (for use with slepc EPSGetInvariantSubspace)
by Jose E. Roman 16 Sep '24

16 Sep '24

Please respond to the list. The Fortran interface is not fully tested. In particular, this function is not employed in any of the Fortran examples. The Fortran interface has undergone successive improvements in recent times. I have tried the following: Vec, pointer :: subspace_receive(:) ... PetscCall(EPSGetConverged(eps, nconv, ierr)) PetscCall(MatCreateVecs(A, PETSC_NULL_VEC, xr, ierr)) PetscCall(VecDuplicateVecsF90(xr, nconv, subspace_receive, ierr)) PetscCall(EPSGetInvariantSubspace(eps, subspace_receive, ierr)) ... PetscCall(VecDestroyVecsF90(nconv, subspace_receive, ierr)) PetscCall(VecDestroy(xr, ierr)) I have tested this in the 'main' branch and it works as expected (note that 'main' will become version 3.22 in two weeks). You can try with version 3.21 and it may work, probably not in older versions. Jose > El 16 sept 2024, a las 11:07, dazza simplythebest <sayosale(a)hotmail.com> escribió: > > Dear All and Jose, > I am now trying to code up the extraction of the invariant subspace, but I seem to be running into > a compilation problem. I am using the slepc example file ex1f.F as a base, and have added the following lines: > > !Declaration > Vec subspace_receive > . > . > . > ! Executed code (following EPSSolve) > call EPSGetConverged(eps,nconv,ierr); CHKERRA(ierr) > call VecDuplicateVecs(xr, nconv, subspace_receive, ierr); CHKERRA(ierr) > call EPSGetInvariantSubspace(eps, subspace_receive, ierr) ; CHKERRA(ierr). > . > . > > However, I get the following compilation error when I include the call EPSGetInvariantSubspace(eps, subspace_receive, ierr) line: > > > minitest.F(133): error #8284: If the actual argument is scalar, the dummy argument shall be scalar unless the actual argument is of type character or is an element of an array that is not assumed shape, pointer, or polymorphic. [B] > call EPSGetInvariantSubspace(eps, subspace_receive, ierr) ; if (ierr .ne. 0) then;call PetscErrorF(ierr,133,"minitest.F");call MPIU_Abort(PETSC_COMM_SELF,ierr);endif > -----------^ > compilation aborted for minitest.F (code 1) > > > If I comment out the call EPSGetInvariantSubspace(eps, subspace_receive, ierr) line the code compiles, links and runs successfully. > I presume that I must be making a silly mistake in the above lines, but I can't quite see what it is. (There is some discussion in the > manual of there being a different situation if the matrix is real but the eigenvalues are complex, but for this case the eigenvalues > are real). Can you offer any guidance ?! I can provide the whole code of course if that would help. > Many thanks and best wishes, > Dan. > > > > > > From: Jose E. Roman <jroman(a)dsic.upv.es> > Sent: Monday, September 2, 2024 7:06 AM > To: dazza simplythebest <sayosale(a)hotmail.com> > Cc: petsc-users(a)mcs.anl.gov <petsc-users(a)mcs.anl.gov> > Subject: Re: [petsc-users] How to create an array of vectors of type Vec (for use with slepc EPSGetInvariantSubspace) > You can use VecDuplicateVecs() to create an array of Vec. > Yes, you can use VecLoad() to read the data from a binary file. > > Jose > > > > El 2 sept 2024, a las 8:56, dazza simplythebest <sayosale(a)hotmail.com> escribió: > > > > Dear All, > > I am seeking to perform a large number of eigenvalue calculations (related to a fluid dynamics > > problem) using the slepc framework. I have written the program that calls slepc and validated its solutions against > > known results for special cases - the Jacobi-Davidson iteration seemed to perform best on this problem. > > Since in the 'production-level' calculation I will step through a parameter space which should > > mean only small changes to the eigenproblem between each calculation, it seems to make sense > > to use the previous solution as the initial guess for the next calculation. > > > > The sequence to execute this strategy would seem to be: > > To write out the current soln subspace: > > > > call EPSGetInvariantSubspace > > [ call VecView -> binary file (if a file storage is required for possible restart)] > > > > To then read in this data to use it as the initial subspace for the next run: > > > > call EPSSetInitialSpace > > > > My question is that EPSGetInvariantSubspace needs an 'array of nconv vectors', how does one create > > such an array ? My understanding is that the Vec structure is just a one-dimensional array (with contents > > possibly scattered in non-contiguous memory locations !?) , so do we just need to create a long Vec > > that stores multiple vectors ? How can we do this (the examples I have seen thus far use a MatCreateVecs to specify the size > > of the Vec)? > > > > If loading such a subspace from a binary file (created maybe by VecView), presumably VecLoad is the > > best way to load such an array of vectors from the binary file ? > > > > Best wishes and many thanks, > > Dan.

1 0