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

How to find the map between the high order coordinates of DMPlex and vertex numbering?
by Zongze Yang 15 May '23

15 May '23

Hi, if I load a `gmsh` file with second-order elements, the coordinates will be stored in a DG-P2 space. After obtaining the coordinates of a cell, how can I map the coordinates to vertex and edge? Below is some code load the gmsh file, I want to know the relation between `cl` and `cell_coords`. ``` import firedrake as fd import numpy as np # Load gmsh file (2rd) plex = fd.mesh._from_gmsh('test-fd-load-p2-rect.msh') cs, ce = plex.getHeightStratum(0) cdm = plex.getCoordinateDM() csec = dm.getCoordinateSection() coords_gvec = dm.getCoordinates() for i in range(cs, ce): cell_coords = cdm.getVecClosure(csec, coords_gvec, i) print(f'coordinates for cell {i} :\n{cell_coords.reshape([-1, 3])}') cl = dm.getTransitiveClosure(i) print('closure:', cl) break ``` Best wishes, Zongze

3 20

Compiling PETSC with Intel OneAPI compilers and OpenMPI
by Vanella, Marcos (Fed) 15 May '23

15 May '23

Hello, I'm trying to compile the PETSc library version 3.19.1 with OpenMPI 4.1.4 and the OneAPI 2022 Update 2 Intel Compiler suite on a Mac with OSX Ventura 13.3.1. I can compile PETSc in debug mode with this configure and make lines. I can run the PETSC tests, which seem fine. When I compile the library in optimized mode, either using -O3 or O1, for example configuring with: $ ./configure --prefix=/opt/petsc-oneapi22u3 --with-blaslapack-dir=/opt/intel/oneapi/mkl/2022.2.1 COPTFLAGS='-m64 -O1 -g -diag-disable=10441' CXXOPTFLAGS='-m64 -O1 -g -diag-disable=10441' FOPTFLAGS='-m64 -O1 -g' LDFLAGS='-m64' --with-debugging=0 --with-shared-libraries=0 --download-make and using mpicc (icc), mpif90 (ifort) from Open MPI, the static lib compiles. Yet, I see right off the bat this segfault error in the first PETSc example: $ make PETSC_DIR=/Users/mnv/Documents/Software/petsc-3.19.1 PETSC_ARCH=arch-darwin-c-opt test /Users/mnv/Documents/Software/petsc-3.19.1/arch-darwin-c-opt/bin/make --no-print-directory -f /Users/mnv/Documents/Software/petsc-3.19.1/gmakefile.test PETSC_ARCH=arch-darwin-c-opt PETSC_DIR=/Users/mnv/Documents/Software/petsc-3.19.1 test /opt/intel/oneapi/intelpython/latest/bin/python3 /Users/mnv/Documents/Software/petsc-3.19.1/config/gmakegentest.py --petsc-dir=/Users/mnv/Documents/Software/petsc-3.19.1 --petsc-arch=arch-darwin-c-opt --testdir=./arch-darwin-c-opt/tests Using MAKEFLAGS: --no-print-directory -- PETSC_ARCH=arch-darwin-c-opt PETSC_DIR=/Users/mnv/Documents/Software/petsc-3.19.1 CC arch-darwin-c-opt/tests/sys/classes/draw/tests/ex1.o In file included from /Users/mnv/Documents/Software/petsc-3.19.1/include/petscsys.h(44), from /Users/mnv/Documents/Software/petsc-3.19.1/src/sys/classes/draw/tests/ex1.c(4): /Users/mnv/Documents/Software/petsc-3.19.1/include/petscsystypes.h(68): warning #2621: attribute "warn_unused_result" does not apply here PETSC_ERROR_CODE_TYPEDEF enum PETSC_ERROR_CODE_NODISCARD { ^ CLINKER arch-darwin-c-opt/tests/sys/classes/draw/tests/ex1 TEST arch-darwin-c-opt/tests/counts/sys_classes_draw_tests-ex1_1.counts not ok sys_classes_draw_tests-ex1_1 # Error code: 139 # [excess:98681] *** Process received signal *** # [excess:98681] Signal: Segmentation fault: 11 (11) # [excess:98681] Signal code: Address not mapped (1) # [excess:98681] Failing at address: 0x7f # [excess:98681] *** End of error message *** # -------------------------------------------------------------------------- # Primary job terminated normally, but 1 process returned # a non-zero exit code. Per user-direction, the job has been aborted. # -------------------------------------------------------------------------- # -------------------------------------------------------------------------- # mpiexec noticed that process rank 0 with PID 0 on node excess exited on signal 11 (Segmentation fault: 11). # -------------------------------------------------------------------------- ok sys_classes_draw_tests-ex1_1 # SKIP Command failed so no diff I see the same segfault error in all PETSc examples. Any help is mostly appreciated, I'm starting to work with PETSc. Our plan is to use the linear solver from PETSc for the Poisson equation on our numerical scheme and test this on a GPU cluster. So also, any guideline on how to interface PETSc with a fortran code and personal experience is also most appreciated! Marcos

4 14

problem for using PCBDDC
by ziming xiong 15 May '23

15 May '23

Hello sir, I am a PhD student and am trying to use the PCBDDC method in petsc to solve the matrix, but the final result is wrong. So I would like to ask you a few questions. First I will describe the flow of my code, I first used the finite element method to build the total matrix in CSR format (total boundary conditions have been imposed), where I did not build the total matrix, but only the parameters ia, ja,value in CSR format, through which the parameters of the metis (xadj, adjncy) are derived. The matrix is successfully divided into 2 subdomains using metis. After getting the global index of the points of each subdomain by the part parameter of metis. I apply ISLocalToGlobalMappingCreate to case mapping and use ISGlobalToLocalMappingApply to convert the global index of points within each process to local index and use MatSetValueLocal to populate the corresponding subdomain matrix for each process. Here I am missing the relationship of the boundary points between subdomains, and by using ISGlobalToLocalMappingApply (I use IS_GTOLM_MASK to get the points outside the subdomains converted to -1) I can get the index of the missing relationship in the global matrix as well as the value. After creating the global MATIS use MatISSetLocalMat to synchronize the subdomain matrix to the global MATIS. After using MatSetValues to add the relationship of the boundary points between subdomains into the global MATIS. The final calculation is performed, but the final result is not correct. My question is: 1. in PetscCall(MatAssemblyBegin(matIS, MAT_FINAL_ASSEMBLY)). PetscCall(MatAssemblyEnd(matIS, MAT_FINAL_ASSEMBLY)). After that, when viewing the matrix by PetscCall(MatView(matIS,PETSC_VIEWER_STDOUT_WORLD));, each process will output the non-zero items of the matrix separately, but this index is the local index is this normal? 2. I found that after using MatSetValues to add the relationship of boundary points between subdomains into the global MATIS, the calculation result does not change. Why is this? Can I interpolate directly into the global MATIS if I know the global matrix index of the missing relations? Best regards, Ziming XIONG

2 1

Help with KSPSetConvergenceTest
by Edoardo alinovi 13 May '23

13 May '23

Hello guys, Today I am about to write a custom convergence test for KSP doing the following job: - if the number of ksp iterations is less than a given threshold, iterate until that threshold is met - if the number of ksp iterations is bigger than the threshold, use the standard convergence checks. As far as I understood from the examples, this is the function I need to use: KSPSetConvergenceTest <https://petsc.org/main/manualpages/KSP/KSPSetConvergenceTest/>(ksp,MyKSPConverged,0,PETSC_NULL_FUNCTION,ierr) With MyKSPConverge a subroutine (I am in fortran) where my custom checks are implemented. Here is how I have coded MyKSPConverged: * subroutine MyKSPConverged(ksp,n,rnorm,flag,dummy,ierr) KSP ksp PetscErrorCode ierr PetscInt n,dummy KSPConvergedReason flag PetscReal rnorm if (n>1) then call KSPConvergedDefault(ksp, n, rnorm, flag, ierr) else flag = 0 endif ierr = 0 end subroutine MyKSPConverged* However I get this error: [0]PETSC ERROR: KSPSolve has not converged, reason DIVERGED_DTOL which I believe to be triggered by *KSPConvergedDefault. *This is odd, as this simulation converge perfectly if I comment ou my dodgy convergence test. *Do you have any advice for me?* *thank you! *

4 16

DMSWARM particle coordinates per rank
by Matthew Young 12 May '23

12 May '23

Does setting up a PIC-type DMSWARM with an associated cell DM guarantee that each MPI rank will own the particles with coordinates inside the bounds of the portion of the grid it owns? --Matt ========================== Matthew Young, PhD (he/him) Research Scientist II Space Science Center University of New Hampshire Matthew.Young(a)unh.edu ==========================

2 5

DMPlex, is there an API to get a list of Boundary points?
by Ferrand, Jesus A. 11 May '23

11 May '23

Greetings. I terms of dm-plex terminology, I need a list points corresponding to the boundary (i.e., height-1 points whose support is of size 1). Sorry if this is trivial, but I've been looking at the list of APIs for DM-Plex and couldn't find/discern something that addresses this need. I guess I could loop through all height-1 points and call DMPlexGetSupportSize() and assemble said list manually. The assumption for the plex I'm working with is that it has had DMPlexInterpolate()called on it. Sincerely: J.A. Ferrand Embry-Riddle Aeronautical University - Daytona Beach - FL Ph.D. Candidate, Aerospace Engineering M.Sc. Aerospace Engineering B.Sc. Aerospace Engineering B.Sc. Computational Mathematics Phone: (386)-843-1829 Email(s): ferranj2(a)my.erau.edu jesus.ferrand(a)gmail.com

2 1

issues with VecSetValues in petsc 3.19
by Edoardo alinovi 10 May '23

10 May '23

Hello Barry, Matt, Jed, I have just installed the latest and greatest version of petsc and I am hitting a problem I did not have in previous releases. Here is the error: *[1]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------[1]PETSC ERROR: Object is in wrong state[1]PETSC ERROR: Not for unassembled vector, did you call VecAssemblyBegin()/VecAssemblyEnd()?[1]PETSC ERROR: WARNING! There are option(s) set that were not used! Could be the program crashed before they were used or a spelling mistake, etc!* the vector x I am filling is created with: - call VecDuplicate(this%rhs, x, ierr) rhs is allocated and I called VecAssemblyBegin()/VecAssemblyEnd() on, do I need to call it on duplicated vectors as well now on? Thank you!

4 17

PCMG questions
by Mark Adams 09 May '23

09 May '23

I have a MG hierarchy that I construct manually with DMRefine and DMPlexExtrude. * The solver works great with chevy/sor but with chevy/sor it converges slowly or I get indefinite PC errors from CG. And the eigen estimates in cheby are really high, like 10-15. * I tried turning galerkin=none and I got this error. Any thoughts on either of these issues? Thanks, Mark [0]PETSC ERROR: --------------------- Error Message -------------------------------------------------------------- [0]PETSC ERROR: Must call DMShellSetGlobalVector() or DMShellSetCreateGlobalVector() [0]PETSC ERROR: WARNING! There are option(s) set that were not used! Could be the program crashed before they were used or a spelling mistake, etc! [0]PETSC ERROR: Option left: name:-ksp_converged_reason (no value) source: command line [0]PETSC ERROR: Option left: name:-mg_levels_esteig_ksp_type value: cg source: command line [0]PETSC ERROR: Option left: name:-mg_levels_pc_type value: sor source: command line [0]PETSC ERROR: Option left: name:-options_left (no value) source: command line [0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting. [0]PETSC ERROR: Petsc Development GIT revision: v3.19.1-224-g9ed82936d20 GIT Date: 2023-05-07 12:33:48 -0400 [0]PETSC ERROR: ./ex96 on a arch-macosx-gnu-O named MarksMac-302.local by markadams Tue May 9 12:26:52 2023 [0]PETSC ERROR: Configure options CFLAGS="-g -Wall" CXXFLAGS="-g -Wall" COPTFLAGS=-O CXXOPTFLAGS=-O --with-cc=/usr/local/opt/llvm/bin/clang --with-cxx=/usr/local/opt/llvm/bin/clang++ --download-mpich --with-strict-petscerrorcode --download-triangle=1 --with-x=0 --with-debugging=0 --download-hdf5=1 PETSC_ARCH=arch-macosx-gnu-O [0]PETSC ERROR: #1 DMCreateGlobalVector_Shell() at /Users/markadams/Codes/petsc/src/dm/impls/shell/dmshell.c:210 [0]PETSC ERROR: #2 DMCreateGlobalVector() at /Users/markadams/Codes/petsc/src/dm/interface/dm.c:1022 [0]PETSC ERROR: #3 DMGetNamedGlobalVector() at /Users/markadams/Codes/petsc/src/dm/interface/dmget.c:377 [0]PETSC ERROR: #4 DMRestrictHook_SNESVecSol() at /Users/markadams/Codes/petsc/src/snes/interface/snes.c:649 [0]PETSC ERROR: #5 DMRestrict() at /Users/markadams/Codes/petsc/src/dm/interface/dm.c:3407 [0]PETSC ERROR: #6 PCSetUp_MG() at /Users/markadams/Codes/petsc/src/ksp/pc/impls/mg/mg.c:1074 [0]PETSC ERROR: #7 PCSetUp() at /Users/markadams/Codes/petsc/src/ksp/pc/interface/precon.c:994 [0]PETSC ERROR: #8 KSPSetUp() at /Users/markadams/Codes/petsc/src/ksp/ksp/interface/itfunc.c:406 [0]PETSC ERROR: #9 KSPSolve_Private() at /Users/markadams/Codes/petsc/src/ksp/ksp/interface/itfunc.c:824 [0]PETSC ERROR: #10 KSPSolve() at /Users/markadams/Codes/petsc/src/ksp/ksp/interface/itfunc.c:1070 [0]PETSC ERROR: #11 SNESSolve_KSPONLY() at /Users/markadams/Codes/petsc/src/snes/impls/ksponly/ksponly.c:48 [0]PETSC ERROR: #12 SNESSolve() at /Users/markadams/Codes/petsc/src/snes/interface/snes.c:4663 [0]PETSC ERROR: #13 main() at ex96.c:433

2 5

About case ksp/tutorial/ex36.cxx
by neil liu 09 May '23

09 May '23

Hello, Petsc Developers, I am trying to compile ksp/tutorial/ex36.cxx like make ex36, it shows an error " Documents/petsc-3.19.1/include/petscdmmoab.h:10:10: fatal error: moab/Core.hpp: No such file or directory #include <moab/Core.hpp> /*I "moab/Core.hpp" I*/ ^~~~~~~~~~~~~~~ compilation terminated. " Did I miss something? In addition, I tried to know how to use DMAddboundary. It seems all the examples for DMAddboundary are related to a PetscFE object. Does that mean DMAddboundary can only be used with PetscFE? Thanks, Xiaodong

2 1

Fwd: Understanding index sets for PCGASM
by LEONARDO MUTTI 09 May '23

09 May '23

---------- Forwarded message --------- Da: LEONARDO MUTTI <leonardo.mutti01(a)universitadipavia.it> Date: mar 9 mag 2023 alle ore 18:29 Subject: Re: [petsc-users] Understanding index sets for PCGASM To: Matthew Knepley <knepley(a)gmail.com> Thank you for your answer, but I am still confused, sorry. Consider https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/ksp/tests/ex71f.F90 on one processor. Let M=12 for the sake of simplicity, i.e. we deal with a 12x12 2D grid, hence, a 144x144 matrix. Let NSubx = 3, so that on the grid we do 3 vertical and 3 horizontal subdivisions. We should obtain 9 subdomains that are grids of 4x4 nodes each, thus corresponding to 9 submatrices of size 16x16. In my run I obtain NSub = 9 (great) and subdomain_IS(i), i=1,...,9, reads: *IS Object: 1 MPI process* * type: general* *Number of indices in set 16* *0 0* *1 1* *2 2* *3 3* *4 12* *5 13* *6 14* *7 15* *8 24* *9 25* *10 26* *11 27* *12 36* *13 37* *14 38* *15 39* *IS Object: 1 MPI process* * type: general* *Number of indices in set 16* *0 4* *1 5* *2 6* *3 7* *4 16* *5 17* *6 18* *7 19* *8 28* *9 29* *10 30* *11 31* *12 40* *13 41* *14 42* *15 43* *IS Object: 1 MPI process* * type: general* *Number of indices in set 16* *0 8* *1 9* *2 10* *3 11* *4 20* *5 21* *6 22* *7 23* *8 32* *9 33* *10 34* *11 35* *12 44* *13 45* *14 46* *15 47* *IS Object: 1 MPI process* * type: general* *Number of indices in set 16* *0 48* *1 49* *2 50* *3 51* *4 60* *5 61* *6 62* *7 63* *8 72* *9 73* *10 74* *11 75* *12 84* *13 85* *14 86* *15 87* *IS Object: 1 MPI process* * type: general* *Number of indices in set 16* *0 52* *1 53* *2 54* *3 55* *4 64* *5 65* *6 66* *7 67* *8 76* *9 77* *10 78* *11 79* *12 88* *13 89* *14 90* *15 91* *IS Object: 1 MPI process* * type: general* *Number of indices in set 16* *0 56* *1 57* *2 58* *3 59* *4 68* *5 69* *6 70* *7 71* *8 80* *9 81* *10 82* *11 83* *12 92* *13 93* *14 94* *15 95* *IS Object: 1 MPI process* * type: general* *Number of indices in set 16* *0 96* *1 97* *2 98* *3 99* *4 108* *5 109* *6 110* *7 111* *8 120* *9 121* *10 122* *11 123* *12 132* *13 133* *14 134* *15 135* *IS Object: 1 MPI process* * type: general* *Number of indices in set 16* *0 100* *1 101* *2 102* *3 103* *4 112* *5 113* *6 114* *7 115* *8 124* *9 125* *10 126* *11 127* *12 136* *13 137* *14 138* *15 139* *IS Object: 1 MPI process* * type: general* *Number of indices in set 16* *0 104* *1 105* *2 106* *3 107* *4 116* *5 117* *6 118* *7 119* *8 128* *9 129* *10 130* *11 131* *12 140* *13 141* *14 142* *15 143* As you said, no number here reaches 144. But the number stored in subdomain_IS are 9x16= #subdomains x 16, whereas I would expect, also given your latest reply, 9x16x16x2=#subdomains x submatrix height x submatrix width x length of a (row,column) pair. It would really help me if you could briefly explain how the output above encodes the subdivision into subdomains. Many thanks again, Leonardo Il giorno mar 9 mag 2023 alle ore 16:24 Matthew Knepley <knepley(a)gmail.com> ha scritto: > On Tue, May 9, 2023 at 10:05 AM LEONARDO MUTTI < > leonardo.mutti01(a)universitadipavia.it> wrote: > >> Great thanks! I can now successfully run >> https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/ksp/tests/ex71f.F90. >> >> Going forward with my experiments, let me post a new code snippet (very >> similar to ex71f.F90) that I cannot get to work, probably I must be >> setting up the IS objects incorrectly. >> >> I have an 8x8 matrix A=diag(1,1,2,2,...,2) and a vector b=(0.5,...,0.5). >> We have only one processor, and I want to solve Ax=b using GASM. In >> particular, KSP is set to preonly, GASM is the preconditioner and it uses >> on each submatrix an lu direct solver (sub_ksp = preonly, sub_pc = lu). >> >> For the GASM algorithm, I divide A into diag(1,1) and diag(2,2,...,2). >> For simplicity I set 0 overlap. Now I want to use GASM to solve Ax=b. The >> code follows. >> >> #include <petsc/finclude/petscmat.h> >> #include <petsc/finclude/petscksp.h> >> #include <petsc/finclude/petscpc.h> >> USE petscmat >> USE petscksp >> USE petscpc >> USE MPI >> >> Mat :: A >> Vec :: b, x >> PetscInt :: M, I, J, ISLen, NSub >> PetscMPIInt :: size >> PetscErrorCode :: ierr >> PetscScalar :: v >> KSP :: ksp >> PC :: pc >> IS :: subdomains_IS(2), inflated_IS(2) >> PetscInt,DIMENSION(4) :: indices_first_domain >> PetscInt,DIMENSION(36) :: indices_second_domain >> >> call PetscInitialize(PETSC_NULL_CHARACTER, ierr) >> call MPI_Comm_size(PETSC_COMM_WORLD, size, ierr) >> >> >> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! >> ! INTRO: create matrix and right hand side >> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! >> >> WRITE(*,*) "Assembling A,b" >> >> M = 8 >> call MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, >> & M, M, PETSC_DEFAULT_INTEGER, PETSC_NULL_INTEGER, >> & PETSC_DEFAULT_INTEGER, PETSC_NULL_INTEGER,A, ierr) >> DO I=1,M >> DO J=1,M >> IF ((I .EQ. J) .AND. (I .LE. 2 )) THEN >> v = 1 >> ELSE IF ((I .EQ. J) .AND. (I .GT. 2 )) THEN >> v = 2 >> ELSE >> v = 0 >> ENDIF >> call MatSetValue(A, I-1, J-1, v, INSERT_VALUES, ierr) >> END DO >> END DO >> >> call MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr) >> call MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr) >> >> call VecCreate(PETSC_COMM_WORLD,b,ierr) >> call VecSetSizes(b, PETSC_DECIDE, M,ierr) >> call VecSetFromOptions(b,ierr) >> >> do I=1,M >> v = 0.5 >> call VecSetValue(b,I-1,v, INSERT_VALUES,ierr) >> end do >> >> call VecAssemblyBegin(b,ierr) >> call VecAssemblyEnd(b,ierr) >> >> >> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! >> ! FIRST KSP/PC SETUP >> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! >> >> WRITE(*,*) "KSP/PC first setup" >> >> call KSPCreate(PETSC_COMM_WORLD, ksp, ierr) >> call KSPSetOperators(ksp, A, A, ierr) >> call KSPSetType(ksp, 'preonly', ierr) >> call KSPGetPC(ksp, pc, ierr) >> call KSPSetUp(ksp, ierr) >> call PCSetType(pc, PCGASM, ierr) >> >> >> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! >> ! GASM, SETTING SUBDOMAINS >> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! >> >> WRITE(*,*) "Setting GASM subdomains" >> >> ! Let's create the subdomain IS and inflated_IS >> ! They are equal if no overlap is present >> ! They are 1: 0,1,8,9 >> ! 2: 10,...,15,18,...,23,...,58,...,63 >> >> indices_first_domain = [0,1,8,9] ! corresponds to diag(1,1) >> do I=0,5 >> do J=0,5 >> indices_second_domain(I*6+1+J) = 18 + J + 8*I ! corresponds >> to diag(2,2,...,2) >> !WRITE(*,*) I*6+1+J, 18 + J + 8*I >> end do >> end do >> >> ! Convert into IS >> ISLen = 4 >> call ISCreateGeneral(PETSC_COMM_WORLD,ISLen,indices_first_domain, >> & PETSC_COPY_VALUES, subdomains_IS(1), ierr) >> call ISCreateGeneral(PETSC_COMM_WORLD,ISLen,indices_first_domain, >> & PETSC_COPY_VALUES, inflated_IS(1), ierr) >> ISLen = 36 >> call ISCreateGeneral(PETSC_COMM_WORLD,ISLen,indices_second_domain, >> & PETSC_COPY_VALUES, subdomains_IS(2), ierr) >> call ISCreateGeneral(PETSC_COMM_WORLD,ISLen,indices_second_domain, >> & PETSC_COPY_VALUES, inflated_IS(2), ierr) >> >> NSub = 2 >> call PCGASMSetSubdomains(pc,NSub, >> & subdomains_IS,inflated_IS,ierr) >> call PCGASMDestroySubdomains(NSub, >> & subdomains_IS,inflated_IS,ierr) >> >> >> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! >> ! GASM: SET SUBSOLVERS >> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! >> >> WRITE(*,*) "Setting subsolvers for GASM" >> >> call PCSetUp(pc, ierr) ! should I add this? >> >> call PetscOptionsSetValue(PETSC_NULL_OPTIONS, >> & "-sub_pc_type", "lu", ierr) >> call PetscOptionsSetValue(PETSC_NULL_OPTIONS, >> & "-sub_ksp_type", "preonly", ierr) >> >> call KSPSetFromOptions(ksp, ierr) >> call PCSetFromOptions(pc, ierr) >> >> >> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! >> ! DUMMY SOLUTION: DID IT WORK? >> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! >> >> WRITE(*,*) "Solve" >> >> call VecDuplicate(b,x,ierr) >> call KSPSolve(ksp,b,x,ierr) >> >> call MatDestroy(A, ierr) >> call KSPDestroy(ksp, ierr) >> call PetscFinalize(ierr) >> >> This code is failing in multiple points. At call PCSetUp(pc, ierr) it >> produces: >> >> *[0]PETSC ERROR: Argument out of range* >> *[0]PETSC ERROR: Scatter indices in ix are out of range* >> *...* >> *[0]PETSC ERROR: #1 VecScatterCreate() at >> ***\src\vec\is\sf\INTERF~1\vscat.c:736* >> *[0]PETSC ERROR: #2 PCSetUp_GASM() at >> ***\src\ksp\pc\impls\gasm\gasm.c:433* >> *[0]PETSC ERROR: #3 PCSetUp() at ***\src\ksp\pc\INTERF~1\precon.c:994* >> >> And at call KSPSolve(ksp,b,x,ierr) it produces: >> >> *forrtl: severe (157): Program Exception - access violation* >> >> >> The index sets are setup coherently with the outputs of e.g. >> https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/ksp/tests/output/ex71f_1…: >> in particular each element of the matrix A corresponds to a number from 0 >> to 63. >> > > This is not correct, I believe. The indices are row/col indices, not > indices into dense blocks, so for > your example, they are all in [0, 8]. > > Thanks, > > Matt > > >> Note that each submatrix does not represent some physical subdomain, the >> subdivision is just at the algebraic level. >> I thus have the following questions: >> >> - is this the correct way of creating the IS objects, given my >> objective at the beginning of the email? Is the ordering correct? >> - what am I doing wrong that is generating the above errors? >> >> Thanks for the patience and the time. >> Best, >> Leonardo >> >> Il giorno ven 5 mag 2023 alle ore 18:43 Barry Smith <bsmith(a)petsc.dev> >> ha scritto: >> >>> >>> Added in *barry/2023-05-04/add-pcgasm-set-subdomains *see also >>> https://gitlab.com/petsc/petsc/-/merge_requests/6419 >>> >>> Barry >>> >>> >>> On May 4, 2023, at 11:23 AM, LEONARDO MUTTI < >>> leonardo.mutti01(a)universitadipavia.it> wrote: >>> >>> Thank you for the help. >>> Adding to my example: >>> >>> >>> * call PCGASMSetSubdomains(pc,NSub, subdomains_IS, >>> inflated_IS,ierr) call >>> PCGASMDestroySubdomains(NSub,subdomains_IS,inflated_IS,ierr)* >>> results in: >>> >>> * Error LNK2019 unresolved external symbol PCGASMDESTROYSUBDOMAINS >>> referenced in function ... * >>> >>> * Error LNK2019 unresolved external symbol PCGASMSETSUBDOMAINS >>> referenced in function ... * >>> I'm not sure if the interfaces are missing or if I have a compilation >>> problem. >>> Thank you again. >>> Best, >>> Leonardo >>> >>> Il giorno sab 29 apr 2023 alle ore 20:30 Barry Smith <bsmith(a)petsc.dev> >>> ha scritto: >>> >>>> >>>> Thank you for the test code. I have a fix in the branch >>>> barry/2023-04-29/fix-pcasmcreatesubdomains2d >>>> <https://gitlab.com/petsc/petsc/-/commits/barry/2023-04-29/fix-pcasmcreatesu…> with >>>> merge request https://gitlab.com/petsc/petsc/-/merge_requests/6394 >>>> >>>> The functions did not have proper Fortran stubs and interfaces so I >>>> had to provide them manually in the new branch. >>>> >>>> Use >>>> >>>> git fetch >>>> git checkout barry/2023-04-29/fix-pcasmcreatesubdomains2d >>>> <https://gitlab.com/petsc/petsc/-/commits/barry/2023-04-29/fix-pcasmcreatesu…> >>>> ./configure etc >>>> >>>> Your now working test code is in src/ksp/ksp/tests/ex71f.F90 I had >>>> to change things slightly and I updated the error handling for the latest >>>> version. >>>> >>>> Please let us know if you have any later questions. >>>> >>>> Barry >>>> >>>> >>>> >>>> >>>> On Apr 28, 2023, at 12:07 PM, LEONARDO MUTTI < >>>> leonardo.mutti01(a)universitadipavia.it> wrote: >>>> >>>> Hello. I am having a hard time understanding the index sets to feed >>>> PCGASMSetSubdomains, and I am working in Fortran (as a PETSc novice). To >>>> get more intuition on how the IS objects behave I tried the following >>>> minimal (non) working example, which should tile a 16x16 matrix into 16 >>>> square, non-overlapping submatrices: >>>> >>>> #include <petsc/finclude/petscmat.h> >>>> #include <petsc/finclude/petscksp.h> >>>> #include <petsc/finclude/petscpc.h> >>>> USE petscmat >>>> USE petscksp >>>> USE petscpc >>>> >>>> Mat :: A >>>> PetscInt :: M, NSubx, dof, overlap, NSub >>>> INTEGER :: I,J >>>> PetscErrorCode :: ierr >>>> PetscScalar :: v >>>> KSP :: ksp >>>> PC :: pc >>>> IS :: subdomains_IS, inflated_IS >>>> >>>> call PetscInitialize(PETSC_NULL_CHARACTER , ierr) >>>> >>>> !-----Create a dummy matrix >>>> M = 16 >>>> call MatCreateAIJ(MPI_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, >>>> & M, M, >>>> & PETSC_DEFAULT_INTEGER, PETSC_NULL_INTEGER, >>>> & PETSC_DEFAULT_INTEGER, PETSC_NULL_INTEGER, >>>> & A, ierr) >>>> >>>> DO I=1,M >>>> DO J=1,M >>>> v = I*J >>>> CALL MatSetValue (A,I-1,J-1,v, >>>> & INSERT_VALUES , ierr) >>>> END DO >>>> END DO >>>> >>>> call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY , ierr) >>>> call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY , ierr) >>>> >>>> !-----Create KSP and PC >>>> call KSPCreate(PETSC_COMM_WORLD,ksp, ierr) >>>> call KSPSetOperators(ksp,A,A, ierr) >>>> call KSPSetType(ksp,"bcgs",ierr) >>>> call KSPGetPC(ksp,pc,ierr) >>>> call KSPSetUp(ksp, ierr) >>>> call PCSetType(pc,PCGASM, ierr) >>>> call PCSetUp(pc , ierr) >>>> >>>> !-----GASM setup >>>> NSubx = 4 >>>> dof = 1 >>>> overlap = 0 >>>> >>>> call PCGASMCreateSubdomains2D(pc, >>>> & M, M, >>>> & NSubx, NSubx, >>>> & dof, overlap, >>>> & NSub, subdomains_IS, inflated_IS, ierr) >>>> >>>> call ISView(subdomains_IS, PETSC_VIEWER_STDOUT_WORLD, ierr) >>>> >>>> call KSPDestroy(ksp, ierr) >>>> call PetscFinalize(ierr) >>>> >>>> Running this on one processor, I get NSub = 4. >>>> If PCASM and PCASMCreateSubdomains2D are used instead, I get NSub = 16 >>>> as expected. >>>> Moreover, I get in the end "forrtl: severe (157): Program Exception - >>>> access violation". So: >>>> 1) why do I get two different results with ASM, and GASM? >>>> 2) why do I get access violation and how can I solve this? >>>> In fact, in C, subdomains_IS, inflated_IS should pointers to IS >>>> objects. As I see on the Fortran interface, the arguments to >>>> PCGASMCreateSubdomains2D are IS objects: >>>> >>>> subroutine PCGASMCreateSubdomains2D(a,b,c,d,e,f,g,h,i,j,z) >>>> import tPC,tIS >>>> PC a ! PC >>>> PetscInt b ! PetscInt >>>> PetscInt c ! PetscInt >>>> PetscInt d ! PetscInt >>>> PetscInt e ! PetscInt >>>> PetscInt f ! PetscInt >>>> PetscInt g ! PetscInt >>>> PetscInt h ! PetscInt >>>> IS i ! IS >>>> IS j ! IS >>>> PetscErrorCode z >>>> end subroutine PCGASMCreateSubdomains2D >>>> Thus: >>>> 3) what should be inside e.g., subdomains_IS? I expect it to contain, >>>> for every created subdomain, the list of rows and columns defining the >>>> subblock in the matrix, am I right? >>>> >>>> Context: I have a block-tridiagonal system arising from space-time >>>> finite elements, and I want to solve it with GMRES+PCGASM preconditioner, >>>> where each overlapping submatrix is on the diagonal and of size 3x3 blocks >>>> (and spanning multiple processes). This is PETSc 3.17.1 on Windows. >>>> >>>> Thanks in advance, >>>> Leonardo >>>> >>>> >>>> >>> > > -- > 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://www.cse.buffalo.edu/~knepley/ > <http://www.cse.buffalo.edu/~knepley/> >

2 2