April 2021 - 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

Question about PCFieldSplit
by Jorti, Zakariae 31 Jan '22

31 Jan '22

Hello, I have a DMStag grid with one dof on each edge and face center. I want to use a PCFieldSplit preconditioner on a Jacobian matrix that I assume is already split but I am not sure how to determine the fields. In the DMStag examples (ex2.c and ex3.c), the function PCFieldSplitSetDetectSaddlePoint is used to determine those fields based on zero diagonal entries. In my case, I have a Jacobian matrix that does not have zero diagonal entries. Can I use that PCFieldSplitSetDetectSaddlePoint in this case? If not, how should I do? Should I do like this example (https://www.mcs.anl.gov/petsc/petsc-master/src/ksp/ksp/tutorials/ex43.c.html) const PetscInt Bfields[1] = {0},Efields[1] = {1}; KSPGetPC(ksp,&pc); PCFieldSplitSetBlockSize(pc,2); PCFieldSplitSetFields(pc,"B",1,Bfields,Bfields); PCFieldSplitSetFields(pc,"E",1,Efields,Efields); where my B unknowns are defined on face centers and E unknowns are defined on edge centers? One last thing, I do not know which field comes first. Is it the one defined for face dofs or edge dofs. Thank you. Best regards, Zakariae

4 14

PETSC installation on Cray
by Enrico 13 Sep '21

13 Sep '21

Hi, I'm having some problems installing PETSC with Cray compiler. I use this configuration: ./configure --with-cc=cc --with-cxx=CC --with-fc=0 --with-debugging=1 --with-shared-libraries=1 COPTFLAGS=-O0 CXXOPTFLAGS=-O0 and when I do make all I get the following error because of cmathcalls.h: CC-1043 craycc: ERROR File = /usr/include/bits/cmathcalls.h, Line = 55 _Complex can only be used with floating-point types. __MATHCALL (cacos, (_Mdouble_complex_ __z)); ^ Am I doing something wrong? Regards, Enrico Degregori

4 5

Rather different matrix product results on multiple processes
by Peder Jørgensgaard Olesen 03 Jun '21

03 Jun '21

Hello, When computing a matrix product of the type R = D.DT using MatMatTransposeMult() I find I get rather different results depending on the number of processes. In one example using a data set that is small compared to the application I get Frobenius norms |R| = 1.047e3 on a single process, 1.0363e3 on a single HPC node (40 cores), and 9.7307e2 on two nodes. I have ascertained that the single process result is indeed the correct one (i.e., eigenvectors of R form a proper basis for the columns of D), so naturally I'd love to be able to reproduce this result across different parallel setups. How might I achieve this? I'm attaching MWE code and the data set used for the example. Thanks in advance! Best Regards Peder Jørgensgaard Olesen PhD Student, Turbulence Research Lab Dept. of Mechanical Engineering Technical University of Denmark Niels Koppels Allé Bygning 403, Rum 105 DK-2800 Kgs. Lyngby

5 7

SLEPc: non-real singular vectors from SVD of real matrix
by Peder Jørgensgaard Olesen 02 Jun '21

02 Jun '21

Hello I've noticed that doing a singular value decomposition of a real matrix appears to result in non-real singular vectors. This should not be the case - singular vectors of a real matrix must be real-valued. In the example attached I read a matrix from a binary file (also attached), perform the SVD, and write singular vectors to an HDF5 file which I subsequently inspect using h5dump, revealing non-zero imaginary parts of vector elements as highlighted below: [pjool@svol mwes]$ h5dump -d U0_sample_svd -c "5,2" sample_svecs.h5 HDF5 "sample_svecs.h5" { DATASET "U0_sample_svd" { DATATYPE H5T_IEEE_F64LE DATASPACE SIMPLE { ( 882, 2 ) / ( 882, 2 ) } SUBSET { START ( 0, 0 ); STRIDE ( 1, 1 ); COUNT ( 5, 2 ); BLOCK ( 1, 1 ); DATA { (0,0): 0.0226108, 0.0299595, (1,0): 0.035414, 0.0469237, (2,0): 0.0276317, 0.0366122, (3,0): 0.0145344, 0.0192581, (4,0): 0.0110376, 0.0146249 } } ATTRIBUTE "complex" { DATATYPE H5T_STD_I32LE DATASPACE SCALAR DATA { (0): 1 } } } } I also extract the imaginary part of the input matrix and print its norm to ensure that the matrix is indeed real. I'm running v3.14, but I don't believe that alone should cause the issue, since it what I'm trying to do appears like a rather common and basic task. What might be the reason behind this behavior, and what can be done to resolve it? Med venlig hilsen / Best Regards Peder Jørgensgaard Olesen PhD Student, Turbulence Research Lab Dept. of Mechanical Engineering Technical University of Denmark Niels Koppels Allé Bygning 403, Rum 105 DK-2800 Kgs. Lyngby

3 4

Questions on matrix-free GMRES implementation
by feng wang 12 May '21

12 May '21

Dear All, I am new to petsc and trying to implement a matrix-free GMRES. I have assembled an approximate Jacobian matrix just for preconditioning. After reading some previous questions on this topic, my approach is: the matrix-free matrix is created as: ierr = MatCreateMFFD(*A_COMM_WORLD, iqe*blocksize, iqe*blocksize, PETSC_DETERMINE, PETSC_DETERMINE, &petsc_A_mf); CHKERRQ(ierr); ierr = MatMFFDSetFunction(petsc_A_mf, FormFunction_mf, this); CHKERRQ(ierr); KSP linear operator is set up as: ierr = KSPSetOperators(petsc_ksp, petsc_A_mf, petsc_A_pre); CHKERRQ(ierr); //petsc_A_pre is my assembled pre-conditioning matrix Before calling KSPSolve, I do: ierr = MatMFFDSetBase(petsc_A_mf, petsc_csv, petsc_baserhs); CHKERRQ(ierr); //petsc_csv is the flow states, petsc_baserhs is the pre-computed right hand side The call back function is defined as: PetscErrorCode cFdDomain::FormFunction_mf(void *ctx, Vec in_vec, Vec out_vec) { PetscErrorCode ierr; cFdDomain *user_ctx; cout << "FormFunction_mf called\n"; //in_vec: flow states //out_vec: right hand side + diagonal contributions from CFL number user_ctx = (cFdDomain*)ctx; //get perturbed conservative variables from petsc user_ctx->petsc_getcsv(in_vec); //get new right side user_ctx->petsc_fd_rhs(); //set new right hand side to the output vector user_ctx->petsc_setrhs(out_vec); ierr = 0; return ierr; } The linear system I am solving is (J+D)x=RHS. J is the Jacobian matrix. D is a diagonal matrix and it is used to stabilise the solution at the start but reduced gradually when the solution moves on to recover Newton's method. I add D*x to the true right side when non-linear function is computed to work out finite difference Jacobian, so when finite difference is used, it actually computes (J+D)*dx. The code runs but diverges in the end. If I don't do matrix-free and use my approximate Jacobian matrix, GMRES works. So something is wrong with my matrix-free implementation. Have I missed something in my implementation? Besides, is there a way to check if the finite difference Jacobian matrix is computed correctly in a matrix-free implementation? Thanks for your help in advance. Feng

3 21

MatDestroy problem with multiple matrices and SUPERLU_DIST
by Deij-van Rijswijk, Menno 06 May '21

06 May '21

Good afternoon, In our code we're using two matrices that are preconditioned using SUPERLU_dist. Upon freeing the second matrix with MatDestroy, the program segfaults with the following stacktrace (see below). This happens with PETSc versions 3.14.5 and 3.15.0, whereas version 3.11.2 does not have this problem. In the code src\mat\impls\aij\mpi\superlu_dist\superlu_dist.c I see that Petsc_Superlu_dist_keyval_Delete_Fn has been added between 3.11.2 and 3.14.5, and the comment reads that it allows reusing the communicator when only a single matrix is used at the time. I'm wondering if using multiple matrices with SUPERLU_dist is problematic here? Note: this is happening on single process and on multiple MPI processes. Best regards, Menno Deij - van Rijswijk #0 0x000015554ffc7db4 in ompi_comm_free () from /home/mdeij/install-gnu/extLibs/lib/libmpi.so.40 #1 0x0000155550021536 in PMPI_Comm_free () from /home/mdeij/install-gnu/extLibs/lib/libmpi.so.40 #2 0x00001555534262ba in superlu_gridexit (grid=0x3689da0) at /home/mdeij/install-gnu/extLibs/Linux-x86_64-Intel/superlu_dist-6.3.0/SRC/superlu_grid.c:174 #3 0x0000155553df094b in Petsc_Superlu_dist_keyval_Delete_Fn (comm=0x26a4e50, keyval=16, attr_val=0x3689d90, extra_state=0x0) at /home/mdeij/build-libs-gnu/superbuild/petsc/src/src/mat/impls/aij/mpi/superlu_dist/superlu_dist.c:97 #4 0x000015554ffc145c in ompi_attr_delete_impl () from /home/mdeij/install-gnu/extLibs/lib/libmpi.so.40 #5 0x000015554ffc3fdf in ompi_attr_delete_all () from /home/mdeij/install-gnu/extLibs/lib/libmpi.so.40 #6 0x000015554ffc7ca7 in ompi_comm_free () from /home/mdeij/install-gnu/extLibs/lib/libmpi.so.40 #7 0x0000155550021536 in PMPI_Comm_free () from /home/mdeij/install-gnu/extLibs/lib/libmpi.so.40 #8 0x00001555538f31c5 in PetscCommDestroy (comm=0x29712e0) at /home/mdeij/build-libs-gnu/superbuild/petsc/src/src/sys/objects/tagm.c:217 #9 0x00001555538f55ab in PetscHeaderDestroy_Private (h=0x29712a0) at /home/mdeij/build-libs-gnu/superbuild/petsc/src/src/sys/objects/inherit.c:121 #10 0x0000155553b44626 in MatDestroy (A=0x26e1c98) at /home/mdeij/build-libs-gnu/superbuild/petsc/src/src/mat/interface/matrix.c:1310 #11 0x0000155553b87a86 in matdestroy_ (x=0x26e1c98, ierr=0x7fffffffc9cc) dr. ir. Menno A. Deij-van Rijswijk | Researcher | Research & Development MARIN | T +31 317 49 35 06 | M.Deij(a)marin.nl<mailto:[email protected]> | www.marin.nl<http://www.marin.nl> [LinkedIn]<https://www.linkedin.com/company/marin> [YouTube] <http://www.youtube.com/marinmultimedia> [Twitter] <https://twitter.com/MARIN_nieuws> [Facebook] <https://www.facebook.com/marin.wageningen> MARIN news: IWSA Open letter release for shipping decarbonisation<https://www.marin.nl/news/iwsa-open-letter-release-for-shipping-decarbonisa…>

3 10

Parallel TS for ODE
by Francesco Brarda 05 May '21

05 May '21

Hi everyone! I am trying to solve a system of 3 ODEs (a basic SIR model) with TS. Sequentially works pretty well, but I need to switch it into a parallel version. I started working with TS not very long time ago, there are few questions I’d like to share with you and if you have any advices I’d be happy to hear. First of all, do I need to use a DM object even if the model is only time dependent? All the examples I found were using that object for the other variable when solving PDEs. When I preallocate the space for the Jacobian matrix, is it better to decide the local or global space? Best, Francesco

5 26

Detecting steady-state with TS
by Salazar De Troya, Miguel 01 May '21

01 May '21

I am solving the signed distance equation \frac{\partial \phi}{\partial t} + sign (\phi_{0})(|\nabla \phi| - 1) = 0 using a Local Discontinuous Galerkin (LDG) method as described in https://www.sciencedirect.com/science/article/pii/S0021999110005255 I am interested in solving it close to steady state. I was hoping I could measure how close to steady state the solution is by using the TSSetEventHandler infrastructure, but the handler does not have information on the time derivative. I looked at TSPSEUDO, but it forces me to use an implicit method, which I cannot provide because how the LDG method works (it calculates the fluxes solving additional equations). This makes me wonder if the LDG method is the best choice, so I am open to suggestions. Given my current progress with the LDG approach, I am wondering if there is a way to solve to steady state using explicit algorithms such as Runge-Kutta. Thanks Miguel Miguel A. Salazar de Troya Postdoctoral Researcher, Lawrence Livermore National Laboratory B141 Rm: 1085-5 Ph: 1(925) 422-6411

3 4

VecDeviceGetArray?
by Mark Adams 28 Apr '21

28 Apr '21

Have we abstracted VecCUDAGetArray[Read] ? VecDeviceGetArray? Mark

2 3