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

Ghost particles for DMSWARM (or similar)
by MIGUEL MOLINOS PEREZ 19 Nov '24

19 Nov '24

3 24

Spring 2024 PETSc meeting and release announcement 3.20
by Barry Smith 09 Nov '24

09 Nov '24

The May 23-24, 2024 PETSc annual meeting in Cologne Germany website is now available: https://cds.uni-koeln.de/en/workshops/petsc-2024/home Please register and start making your travel plans. We are pleased to announce the release of PETSc version 3.20.0 at https://petsc.org/release/download/ A list of the major changes and updates can be found at https://petsc.org/release/changes/320/ The release tarball location is changed from https://ftp.mcs.anl.gov/pub/petsc to https://web.cels.anl.gov/projects/petsc/download The final update to petsc-3.19 i.e petsc-3.19.6 is also available We recommend upgrading to PETSc 3.20.0 soon. As always, please report problems to petsc-maint(a)mcs.anl.gov and ask questions at petsc-users(a)mcs.anl.gov This release includes contributions from David Ham Iglesia Dolci Jack Betteridge Koki Sagiyama Marcin Rogowski Nathan Collier Aidan Hamilton Albert Cowie Barry Smith Blaise Bourdin Connor Ward David Andrs David Wells Duncan Campbell Hansol Suh Hong Zhang Jacob Faibussowitsch James Wright Jed Brown jeremy theler Joseph Pusztay Jose Roman Junchao Zhang Lisandro Dalcin Mark Adams Matthew Knepley Mr. Hong Zhang Pablo Brubeck Pierre Jolivet Pieter Ghysels Richard Tran Mills Satish Balay Sebastian Grimberg Stefano Zampini Suyash Tandon Toby Isaac YANG Zongze Zach Atkins and bug reports/proposed improvements received from Carl-Johan Thore Chris Douglas Daniel@OLCF Duan Junming Edoardo alinovi Philip Fackler Frederic Vi Getnet Betrie Aron Graham Hongrui Yu James Wright Jin Chen Jonas Lundgren Karsten Lettmann Boris Kaus LEONARDO MUTTI Maciej Lewkowicz maitri ksh Milan Pelletier Mr. Hong Zhang Ngoc Mai Monica Huynh Peng Wang Sébastien Gilles Daniel Taller Vibhas Mishra Wei Liu Li Xiaoye As always, thanks for your support, Barry

1 2

PetscSortInt() in PETSc 3.22
by Adrian Croucher 30 Sep '24

30 Sep '24

hi I just pulled PETSc 3.22 and tried to build my (Fortran) code with it. I'm getting compilation errors when it gets to calls to PetscSortInt(): 567 | call PetscSortInt(na, asort, ierr); CHKERRQ(ierr) | 1 Error: Type mismatch in argument ‘a’ at (1); passed INTEGER(4) to INTEGER(8) Here na is PetscInt, asort an array of PetscInt and ierr is PetscErrorCode. Do I need to make changes to get this to work again? - Adrian -- Dr Adrian Croucher Senior Research Fellow Department of Engineering Science Waipapa Taumata Rau / University of Auckland, New Zealand email: a.croucher(a)auckland.ac.nz tel: +64 (0)9 923 4611

2 1

Questions DMPlex
by Karthikeyan Chockalingam - STFC UKRI 30 Sep '24

30 Sep '24

Hi, We have been using PETSc’s block version of the AIJ matrix format to implement fully coupled multiphysics problems using finite elements. Are there specific advantages to moving toward DMPlex for finite element-based coupled multi-physics implementation? Can you still access Hypre? Does the format support running on GPUs? Thank you. 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

Re: [petsc-users] Ask for help in tracking down bugs - setting #threads for pardiso block LU solver
by Barry Smith 30 Sep '24

30 Sep '24

I should have read your initial email more carefully: ---- Loop over block PCs { PCSetType(subPc, PCLU); set subPC type PCFactorSetMatSolverType(subPc, MATSOLVERMKL_PARDISO); set solver for LU SetPardisoParams("pardiso"); set mkl related options including -mat_mkl_pardiso_65 equal to the desired #threads; return code after PetscOptionsSetValue is ok, so I presume the option is set correctly. PCSetFromOptions(subPc); } Then what I see under the debugger PCSetUp after a series of calls gets to MatLUFactorSymbolic_AIJMKL_PARDISO which in turn calls MatSetFromOptions_MKL_PARDISO(F, A) and there in these lines it seems to ignore what I’ve set for -mat_mkl_pardiso_65 and sets #threads to the default 1 PetscCall(PetscOptionsInt("-mat_mkl_pardiso_65", "Suggested number of threads to use within PARDISO", "None", threads, &threads, &flg)); if (flg) PetscSetMKL_PARDISOThreads((int)threads); ----- PCSetFromOptions() actually does not call MatSetFromOptions() on the factored matrix. There is no call to MatSetFromOptions() from PCSetFromOptions() or PetscSetFromOptions_Factor() because the factored matrix F does not yet exist (it is created in PCSetUp()). This is why your call to PCSetFromOptions(subPc); does not process any of your mat mkl options. Since MatSetFromOptions() cannot be called from PCSetFromOptions(), MatSetFromOptions_MKL_PARDISO(F, A) is called in MatFactorSymbolic_AIJMKL_PARDISO_Private() as you noted when the PC is being set up. You say "here in these lines it seems to ignore what I’ve set for -mat_mkl_pardiso_65 and sets #threads to the default 1" This is the part a do not understand, if your values have been inserted in the database earlier in the code and remain in the database they should get processed inside the MatSetFromOptions_MKL_PARDISO(F, A) called by MatFactorSymbolic_AIJMKL_PARDISO_Private(). Can you please put a break point at MatSetFromOptions_MKL_PARDISO(F, A) and then when it breaks put call PetscOptionsView(0,0) it should print out the contents of the options database so you can verify if the options you added are there. You can then step through the code in MatSetFromOptions_MKL_PARDISO to see what is happening as it tries to acess the database to see why it skips PetscCall(PetscOptionsInt("-mat_mkl_pardiso_65", "Suggested number of threads to use within MKL PARDISO", "None", threads, &threads, &flg)); if (flg) PetscSetMKL_PARDISOThreads((int)threads); > On Sep 27, 2024, at 11:58 AM, Voinov, Boris <boris.voinov(a)intel.com> wrote: > > Thanks a lot Barry, > > It looks like something wrong with my own settings though. It doesn’t make an effect. I mean this call > SetPardisoParams("pardiso"); > followed by > PCSetFromOptions(subPc); > doesn’t end up with a call of PetscSetMKL_PARDISOThreads() > This function is called later when PCSetUp is invoked. > > Eventually what I see in the end is > Option left: name:-mat_mkl_pardiso_65 value: 8 source: code listed among unused database options > That is what I’ve set with my own settings, meaning it’s in the database, I guess. > > Could it be that it hits a “not proper” section in the options database and is not being used by PCSetUp? > > Best regards, > Boris > > From: Barry Smith <bsmith(a)petsc.dev <mailto:[email protected]>> > Sent: Friday, September 27, 2024 4:21 PM > To: Voinov, Boris <boris.voinov(a)intel.com <mailto:[email protected]>> > Cc: petsc-users(a)mcs.anl.gov <mailto:[email protected]> > Subject: Re: [petsc-users] Ask for help in tracking down bugs - setting #threads for pardiso block LU solver > > > Boris, > > There is nothing obviously wrong with what you outline, so let's dive down into the code with the debugger and see what is happening. > > PetscSetMKL_PARDISOThreads() set from options is called when the command line option is provided, and this calls > > PETSC_EXTERN void PetscSetMKL_PARDISOThreads(int threads) > { > mkl_domain_set_num_threads(threads, MKL_DOMAIN_PARDISO); > } > > so in the debugger can you check that mkl_domain_set_num_threads() gets called with your requested number of threads? > > BTW: Since you are hardwiring the use of Pardiso in the code with function calls you couldsimply call PetscSetMKL_PARDISOThreads() directly > after your line of code > > PCFactorSetMatSolverType(subPc, MATSOLVERMKL_PARDISO); set solver for LU > > instead of feeding it in through the options database. > > Let us know how this turns out > > Barry > > > > > > > On Sep 27, 2024, at 7:24 AM, Voinov, Boris <boris.voinov(a)intel.com <mailto:[email protected]>> wrote: > > Hello, > > Could you please help me to figure out what’s wrong with the way I’m trying to make mkl pardiso use the number of threads more than one while setting it as a sub PC solver for ASM preconditioner. > Here is how I do this > Loop over block PCs { > PCSetType(subPc, PCLU); set subPC type > PCFactorSetMatSolverType(subPc, MATSOLVERMKL_PARDISO); set solver for LU > SetPardisoParams("pardiso"); set mkl related options including -mat_mkl_pardiso_65 equal to the desired #threads; return code after PetscOptionsSetValue is ok, so I presume the option is set correctly. > PCSetFromOptions(subPc); > } > > Then what I see under the debugger > PCSetUp after a series of calls gets to MatLUFactorSymbolic_AIJMKL_PARDISO which in turn calls MatSetFromOptions_MKL_PARDISO(F, A) and there > in these lines it seems to ignore what I’ve set for -mat_mkl_pardiso_65 and sets #threads to the default 1 > PetscCall(PetscOptionsInt("-mat_mkl_pardiso_65", "Suggested number of threads to use within PARDISO", "None", threads, &threads, &flg)); > if (flg) PetscSetMKL_PARDISOThreads((int)threads); > > This is about petsc-3.20. > PS when I do this in a standalone code which reads a matrix and solves the linear system and set -mat_mkl_pardiso_65 in the command line it’s all right but I need the number of threads to be set in other way in my code. > > Thank you and best regards, > Boris > > ------------------------------------------------------------- > Intel Ireland Limited (Branch) > Collinstown Industrial Park, Leixlip, County Kildare, Ireland > Registered Number: E902934 > > This e-mail and any attachments may contain confidential material for > the sole use of the intended recipient(s). Any review or distribution > by others is strictly prohibited. If you are not the intended > recipient, please contact the sender and delete all copies. > > > ------------------------------------------------------------- > Intel Ireland Limited (Branch) > Collinstown Industrial Park, Leixlip, County Kildare, Ireland > Registered Number: E902934 > > This e-mail and any attachments may contain confidential material for > the sole use of the intended recipient(s). Any review or distribution > by others is strictly prohibited. If you are not the intended > recipient, please contact the sender and delete all copies. >

2 3

Ask for help in tracking down bugs - setting #threads for pardiso block LU solver
by Voinov, Boris 27 Sep '24

27 Sep '24

Hello, Could you please help me to figure out what's wrong with the way I'm trying to make mkl pardiso use the number of threads more than one while setting it as a sub PC solver for ASM preconditioner. Here is how I do this Loop over block PCs { PCSetType(subPc, PCLU); set subPC type PCFactorSetMatSolverType(subPc, MATSOLVERMKL_PARDISO); set solver for LU SetPardisoParams("pardiso"); set mkl related options including -mat_mkl_pardiso_65 equal to the desired #threads; return code after PetscOptionsSetValue is ok, so I presume the option is set correctly. PCSetFromOptions(subPc); } Then what I see under the debugger PCSetUp after a series of calls gets to MatLUFactorSymbolic_AIJMKL_PARDISO which in turn calls MatSetFromOptions_MKL_PARDISO(F, A) and there in these lines it seems to ignore what I've set for -mat_mkl_pardiso_65 and sets #threads to the default 1 PetscCall(PetscOptionsInt("-mat_mkl_pardiso_65", "Suggested number of threads to use within PARDISO", "None", threads, &threads, &flg)); if (flg) PetscSetMKL_PARDISOThreads((int)threads); This is about petsc-3.20. PS when I do this in a standalone code which reads a matrix and solves the linear system and set -mat_mkl_pardiso_65 in the command line it's all right but I need the number of threads to be set in other way in my code. Thank you and best regards, Boris ------------------------------------------------------------- Intel Ireland Limited (Branch) Collinstown Industrial Park, Leixlip, County Kildare, Ireland Registered Number: E902934 This e-mail and any attachments may contain confidential material for the sole use of the intended recipient(s). Any review or distribution by others is strictly prohibited. If you are not the intended recipient, please contact the sender and delete all copies.

2 1

Re: [petsc-users] MatCopy question
by Junchao Zhang 26 Sep '24

26 Sep '24

No, Setting a bigger matrix size won't have a negative impact. Only the actual number of nonzeros matters. --Junchao Zhang On Thu, Sep 26, 2024 at 1:33 PM Frank Bramkamp <bramkamp(a)nsc.liu.se> wrote: > Thanks for the reply, > > when setting the size larger than needed, would that have any negative > impact somewhere else ?! > We typically pre-determine how many non-zero elements are present, when > using MatCreate. > > Then, the smaller size needs the same memory as the larger matrix ?! > > Thanks, Frank > > > > > On 26 Sep 2024, at 20:27, Junchao Zhang <junchao.zhang(a)gmail.com> wrote: > > Is it possible to set the size of the COMPACT matrix as large as the WIDE > matrix, then do MatAYPX(WIDE, 0, COMPACT, SUBSET_NONZERO_PATTERN), i.e, > WIDE = 0*WIDE + COMPACT? You need to build WIDE upfront once, but if you > repeat, the other MatAYPX will be fast. > > --Junchao Zhang > > > On Thu, Sep 26, 2024 at 1:13 PM Barry Smith <bsmith(a)petsc.dev> wrote: > >> >> Frank, >> >> Sorry for the confusion. MatCopy() and MatAXPY() only work for >> identically sized matrices. >> >> Even if the two matrices are the same size but one, say A, has a >> subset of nonzeros than the other, say B, what you want to do is >> non-trivial. This is because to do a copy you must MatAssembly the matrix >> first and this squeezes out the zero locations that you would like to fill >> in B. >> >> A simple solution is in the code that computes the subset of entries, >> call MatSetValues() twice, once for A and once for B. Not entirely >> satisfying but it will work and be reasonably efficient. >> >> >> Barry >> >> >> > On Sep 26, 2024, at 1:50 PM, Frank Bramkamp <bramkamp(a)nsc.liu.se> >> wrote: >> > >> > Dear PETSc team, >> > >> > >> > I would like to create two matrices. The first matrix is based on a >> compact stencil (COMPACT), and the second one has an extended stencil (WIDE) >> > So the non-zero elements of the COMPACT matrix is contained in the >> WIDE matrix as well. But the WIDE matrix shall contain some additional >> elements. >> > I try to create the COMPACT matrix first, then copy it to the WIDE >> matrix and then add some more terms to the WIDE matrix. >> > We want to create a matrix that is based on a lower order approximation >> of our scheme which serves as basis matrix for preconditioning, >> > where we want to have a more accurate jacobian which has a larger >> stencil, that is based on the lower approximation plus some additional >> non-zero elements. >> > We want to avoid that we have to call the same routines twice to setup >> the more accurate matrix, but copy the elements from the lower order >> (compact) matrix. >> > >> > >> > I try to use MatCopy to copy the COMPACT matrix into the WIDE matrix. >> Somehow I cannot figure out to do this, >> > since there is always a problem wirh non conforming sizes >> > >> > [0]PETSC ERROR: Nonconforming object sizes >> > [0]PETSC ERROR: Mat A,Mat B: global dim (4,8) (4,8) >> > >> > >> > Here is a basic example that I try. >> > Adding additional elements to the WIDE matrix is disabled right now, >> since MatCopy already gives a problem. >> > I thought that in MatCopy(COMPACT, WIDE, DIFFERENT_NONZERO_PATTERN, >> ierr) >> > the option DIFFERENT_NONZERO_PATTERN should enable to copy different >> matrix sizes into each other. >> > Or is that not possible to do ?! >> > >> > Here the COMPACT matrix is a 4x4 matrix and the WIDE matrix 8x8, for a >> basic test. >> > >> > Thanks, Frank Bramkamp >> > >> > >> > >> > >> > program petsc_matrix_example >> > #include <petsc/finclude/petsc.h> >> > use petsc >> > implicit none >> > >> > PetscErrorCode :: ierr >> > Mat :: COMPACT, WIDE >> > PetscInt :: m, n, MM, NN, d_nz >> > PetscInt :: i, j >> > PetscScalar :: additional_value >> > >> > ! Initialize PETSc >> > call PetscInitialize(PETSC_NULL_CHARACTER, ierr) >> > >> > ! Set dimensions for COMPACT matrix >> > m = 4 ! Local rows >> > n = 4 ! Local columns >> > MM = 4 ! Global rows >> > NN = 4 ! Global columns >> > d_nz = 4 ! Estimated non-zeros per row (diagonal) o_nz = 0 ! Estimated >> non-zeros per row (off-diagonal) >> > >> > ! Create COMPACT matrix >> > call MatCreateSeqAIJ(PETSC_COMM_SELF, MM, NN, PETSC_DECIDE, >> PETSC_NULL_INTEGER, COMPACT, ierr) >> > >> > ! Set some values in COMPACT matrix (example) do i = 0, M-1 >> > do i = 0, NN-1 >> > !do j = max(0, i-1), min(NN-1, i+1) >> > do j = 0,NN-1 >> > call MatSetValue(COMPACT, i, j, 1.0_PETSC_REAL_KIND, >> INSERT_VALUES, ierr) >> > end do >> > end do >> > >> > ! Assemble COMPACT matrix >> > call MatAssemblyBegin(COMPACT, MAT_FINAL_ASSEMBLY, ierr) >> > call MatAssemblyEnd(COMPACT, MAT_FINAL_ASSEMBLY, ierr) >> > >> > ! Set dimensions for WIDE matrix >> > n = 8 ! Increase local columns >> > NN = 8 ! Increase global columns >> > m = 8 >> > MM = 8 >> > d_nz = 8 >> > >> > ! Create WIDE matrix >> > call MatCreateSeqAIJ(PETSC_COMM_SELF, MM, NN, PETSC_DECIDE, >> PETSC_NULL_INTEGER, WIDE, ierr) >> > >> > ! Assemble WIDE matrix: MUST WE ASSMBLE THE MATRIX BEFORE MatCOPY ?! >> > call MatAssemblyBegin(WIDE, MAT_FINAL_ASSEMBLY, ierr) >> > call MatAssemblyEnd(WIDE, MAT_FINAL_ASSEMBLY, ierr) >> > >> > ! Copy elements from COMPACT to WIDE >> > call MatCopy(COMPACT, WIDE, DIFFERENT_NONZERO_PATTERN, ierr) >> > >> > >> > >> > !call MatAssemblyBegin(WIDE, MAT_FINAL_ASSEMBLY, ierr) >> > >> > ! Add additional elements to WIDE matrix >> > !!$ additional_value = 10.0 >> > !!$ >> > !!$ ! original >> > !!$ !do i = 0, MM-1 >> > !!$ ! do j = NN/2, NN-1 >> > !!$ ! call MatSetValue(WIDE, i, j, additional_value, INSERT_VALUES, >> ierr) >> > !!$ ! end do >> > !!$ !end do >> > !!$ >> > !!$ do i = 0, NN-1 >> > !!$ do j = max(0, i-1), min(NN-1, i+1) >> > !!$ call MatSetValue(WIDE, i, j, 10.0_PETSC_REAL_KIND, ADD_VALUES, >> ierr) >> > !!$ end do >> > !!$ end do >> > >> > >> > ! Assemble WIDE matrix again >> > call MatAssemblyBegin(WIDE, MAT_FINAL_ASSEMBLY, ierr) >> > call MatAssemblyEnd(WIDE, MAT_FINAL_ASSEMBLY, ierr) >> > >> > ! View matrices (optional) >> > call PetscPrintf(PETSC_COMM_WORLD, "COMPACT matrix:\n", ierr) >> > call MatView(COMPACT, PETSC_VIEWER_STDOUT_WORLD, ierr) >> > call PetscPrintf(PETSC_COMM_WORLD, "WIDE matrix:\n", ierr) >> > call MatView(WIDE, PETSC_VIEWER_STDOUT_WORLD, ierr) >> > >> > ! Clean up >> > call MatDestroy(COMPACT, ierr) >> > call MatDestroy(WIDE, ierr) >> > >> > ! Finalize PETSc >> > call PetscFinalize(ierr) >> > >> > end program petsc_matrix_example >> > >> > >> >> >

1 0

MatCopy question
by Frank Bramkamp 26 Sep '24

26 Sep '24

Dear PETSc team, I would like to create two matrices. The first matrix is based on a compact stencil (COMPACT), and the second one has an extended stencil (WIDE) So the non-zero elements of the COMPACT matrix is contained in the WIDE matrix as well. But the WIDE matrix shall contain some additional elements. I try to create the COMPACT matrix first, then copy it to the WIDE matrix and then add some more terms to the WIDE matrix. We want to create a matrix that is based on a lower order approximation of our scheme which serves as basis matrix for preconditioning, where we want to have a more accurate jacobian which has a larger stencil, that is based on the lower approximation plus some additional non-zero elements. We want to avoid that we have to call the same routines twice to setup the more accurate matrix, but copy the elements from the lower order (compact) matrix. I try to use MatCopy to copy the COMPACT matrix into the WIDE matrix. Somehow I cannot figure out to do this, since there is always a problem wirh non conforming sizes [0]PETSC ERROR: Nonconforming object sizes [0]PETSC ERROR: Mat A,Mat B: global dim (4,8) (4,8) Here is a basic example that I try. Adding additional elements to the WIDE matrix is disabled right now, since MatCopy already gives a problem. I thought that in MatCopy(COMPACT, WIDE, DIFFERENT_NONZERO_PATTERN, ierr) the option DIFFERENT_NONZERO_PATTERN should enable to copy different matrix sizes into each other. Or is that not possible to do ?! Here the COMPACT matrix is a 4x4 matrix and the WIDE matrix 8x8, for a basic test. Thanks, Frank Bramkamp program petsc_matrix_example #include <petsc/finclude/petsc.h> use petsc implicit none PetscErrorCode :: ierr Mat :: COMPACT, WIDE PetscInt :: m, n, MM, NN, d_nz PetscInt :: i, j PetscScalar :: additional_value ! Initialize PETSc call PetscInitialize(PETSC_NULL_CHARACTER, ierr) ! Set dimensions for COMPACT matrix m = 4 ! Local rows n = 4 ! Local columns MM = 4 ! Global rows NN = 4 ! Global columns d_nz = 4 ! Estimated non-zeros per row (diagonal) o_nz = 0 ! Estimated non-zeros per row (off-diagonal) ! Create COMPACT matrix call MatCreateSeqAIJ(PETSC_COMM_SELF, MM, NN, PETSC_DECIDE, PETSC_NULL_INTEGER, COMPACT, ierr) ! Set some values in COMPACT matrix (example) do i = 0, M-1 do i = 0, NN-1 !do j = max(0, i-1), min(NN-1, i+1) do j = 0,NN-1 call MatSetValue(COMPACT, i, j, 1.0_PETSC_REAL_KIND, INSERT_VALUES, ierr) end do end do ! Assemble COMPACT matrix call MatAssemblyBegin(COMPACT, MAT_FINAL_ASSEMBLY, ierr) call MatAssemblyEnd(COMPACT, MAT_FINAL_ASSEMBLY, ierr) ! Set dimensions for WIDE matrix n = 8 ! Increase local columns NN = 8 ! Increase global columns m = 8 MM = 8 d_nz = 8 ! Create WIDE matrix call MatCreateSeqAIJ(PETSC_COMM_SELF, MM, NN, PETSC_DECIDE, PETSC_NULL_INTEGER, WIDE, ierr) ! Assemble WIDE matrix: MUST WE ASSMBLE THE MATRIX BEFORE MatCOPY ?! call MatAssemblyBegin(WIDE, MAT_FINAL_ASSEMBLY, ierr) call MatAssemblyEnd(WIDE, MAT_FINAL_ASSEMBLY, ierr) ! Copy elements from COMPACT to WIDE call MatCopy(COMPACT, WIDE, DIFFERENT_NONZERO_PATTERN, ierr) !call MatAssemblyBegin(WIDE, MAT_FINAL_ASSEMBLY, ierr) ! Add additional elements to WIDE matrix !!$ additional_value = 10.0 !!$ !!$ ! original !!$ !do i = 0, MM-1 !!$ ! do j = NN/2, NN-1 !!$ ! call MatSetValue(WIDE, i, j, additional_value, INSERT_VALUES, ierr) !!$ ! end do !!$ !end do !!$ !!$ do i = 0, NN-1 !!$ do j = max(0, i-1), min(NN-1, i+1) !!$ call MatSetValue(WIDE, i, j, 10.0_PETSC_REAL_KIND, ADD_VALUES, ierr) !!$ end do !!$ end do ! Assemble WIDE matrix again call MatAssemblyBegin(WIDE, MAT_FINAL_ASSEMBLY, ierr) call MatAssemblyEnd(WIDE, MAT_FINAL_ASSEMBLY, ierr) ! View matrices (optional) call PetscPrintf(PETSC_COMM_WORLD, "COMPACT matrix:\n", ierr) call MatView(COMPACT, PETSC_VIEWER_STDOUT_WORLD, ierr) call PetscPrintf(PETSC_COMM_WORLD, "WIDE matrix:\n", ierr) call MatView(WIDE, PETSC_VIEWER_STDOUT_WORLD, ierr) ! Clean up call MatDestroy(COMPACT, ierr) call MatDestroy(WIDE, ierr) ! Finalize PETSc call PetscFinalize(ierr) end program petsc_matrix_example

3 2

DMPlex boundary condition from array
by Pierre Seize 25 Sep '24

25 Sep '24

Hi I would like to solve div(u) = 0 with a FE method in a spherical mesh made of tetrahedron. My problem concern the boundary condition : a Dirichlet boundary condition on the surface of the sphere where the value is given by another program, in the form of a list of indices corresponding to the vertices on the boundary and a list of values to be applied on those vertices. I am able to use DMPlexCreateFromCellListPetsc to create the DM on process 0, but the moment I distribute the DM the vertices are renumbered and I loose the correspondence with the boundary describing arrays. I have two questions: - is there a way to get the mapping between old indices and new global indices after the distribute operation? - how can I apply a boundary condition such as mine, i.e. set given values at given vertices? Thanks in advance Pierre Seize

1 0

Compiling 3.21.5 version on a Grace Hopper node failed
by Qiyue Lu 25 Sep '24

25 Sep '24

Hello, I am trying to build the PETSc version 3.21.5 (the latest, cloned from git) version on a Grace Hopper node with this configuration options: *--with-cc=cc --with-cxx=CC --with-fc=ftn --download-f2cblaslapack=1 --with-cudac=nvcc --with-cuda=1 --with-threadsafety=1* with module cudatoolkit/23.11_12.3. And I got this error: =============== /u/qiyuelu1/petsc/src/vec/vec/impls/seq/cupm/cuda/../vecseqcupm_impl.hpp: In instantiation of ‘Petsc::vec::cupm::impl::VecSeq_CUPM<<anonymous> >::SetAsync(Vec, PetscScalar, PetscDeviceContext)::<lambda()> [with Petsc::device::cupm::DeviceType T = (Petsc::device::cupm::DeviceType)0]’: /u/qiyuelu1/petsc/src/vec/vec/impls/seq/cupm/cuda/../vecseqcupm_impl.hpp:1313:157: required from ‘struct Petsc::vec::cupm::impl::VecSeq_CUPM<<anonymous> >::SetAsync(Vec, PetscScalar, PetscDeviceContext) [with Petsc::device::cupm::DeviceType T = (Petsc::device::cupm::DeviceType)0; PetscErrorCode = int; Vec = _p_Vec*; PetscScalar = double; PetscDeviceContext = _p_PetscDeviceContext*]::<lambda()>’ /u/qiyuelu1/petsc/src/vec/vec/impls/seq/cupm/cuda/../vecseqcupm_impl.hpp:1313:167: required from ‘static PetscErrorCode Petsc::vec::cupm::impl::VecSeq_CUPM<<anonymous> >::SetAsync(Vec, PetscScalar, PetscDeviceContext) [with Petsc::device::cupm::DeviceType T = (Petsc::device::cupm::DeviceType)0; PetscErrorCode = int; Vec = _p_Vec*; PetscScalar = double; PetscDeviceContext = _p_PetscDeviceContext*]’ /u/qiyuelu1/petsc/src/vec/vec/impls/seq/cupm/cuda/vecseqcupm.cu:7:40: required from here /u/qiyuelu1/petsc/src/vec/vec/impls/seq/cupm/cuda/../vecseqcupm_impl.hpp:1313:150: internal compiler error: in maybe_undo_parenthesized_ref, at cp/semantics.c:1740 PetscCallThrust(THRUST_CALL(thrust::fill, stream, dptr, dptr + n, alpha)); =============== Could you please show some hints on this issue? Thanks, Qiyue Lu

3 8