January 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

How to add a source term for PETSCFV ?
by Thibault Bridel-Bertomeu 14 Apr '21

14 Apr '21

Dear all, I have been studying ex11.c from ts/tutorials to understand how to solve an hyperbolic system of equations using PETSCFV. I first worked on the Euler equations for inviscid fluids and based on what ex11.c presents, I was able to add the right PETSc instructions in an already existing in-house code with different gas models to solve the problems in parallel (MPI) and with the AMR capabilities offered by P4EST. Now my goal is to move to Navier-Stokes equations. Theoretically the system is not completely hyperbolic and can be seen as one with an hyperbolic part (identical to the Euler equations) and a parabolic part coming from the RHS diffusion terms. I have been looking into the manual and also the sources of PETSc around the DM, DMPlex, DS and FV classes but I could not find anything that speaks to me as "adding a RHS to an hyperbolic system of equations" or "adding a source term to an hyperbolic system of equations". What's more, that source term depends on the derivatives of the context variables ... I wanted to know if anyone maybe had a suggestion regarding this issue ? Thank you very much in advance, Thibault Bridel-Bertomeu — Eng, MSc, PhD Research Engineer CEA/CESTA 33114 LE BARP Tel.: (+33)557046924 Mob.: (+33)611025322 Mail: thibault.bridelbertomeu(a)gmail.com

4 7

TSSetEventHandler and TSSetPostEventIntervalStep
by Blondel, Sophie 05 Apr '21

05 Apr '21

Hi, I am currently using TSSetEventHandler in my code to detect a random event where the solution vector gets modified during the event. Ideally, after the event happens I want the solver to use a much smaller timestep using TSSetPostEventIntervalStep. However, when I use TSSetPostEventIntervalStep the solver doesn't use the set value. I managed to reproduce the behavior by modifying ex40.c as attached. I think the issue is related to the fact that the fvalue is not technically "approaching" 0 with a random event, it is more of a step function instead. Do you have any recommendation on how to implement the behavior I'm looking for? Let me know if I can provide additional information. Best, Sophie

4 11

creation of parallel dmplex from a partitioned mesh
by Cameron Smith 01 Mar '21

01 Mar '21

Hello, We have a partitioned mesh that we want to create a DMPlex from that has the same distribution of elements (i.e., assignment of elements to processes) and vertices 'interior' to a process (i.e., vertices not on the inter-process boundary). We were trying to use DMPlexCreateFromCellListParallelPetsc() or DMPlexBuildFromCellListParallel() and found that the vertex ownership (roots in the returned Vertex SF) appears to be sequentially assigned to processes based on global vertex id. In general, this will not match our mesh distribution. As we understand, to subsequently set vertex coordinates (or other vertex data) we would have to utilize a star forest (SF) communication API to send data to the correct process. Is that correct? Alternatively, if we create a dmplex object from the elements that exist on each process using DMCreateFromCellList(), and then create a SF from mesh vertices on inter-process boundaries (using the mapping from local to global vertex ids provided by our mesh library), could we then associate the dmplex objects with the SF? Is it as simple as calling DMSetPointSF()? If manually defining the PointSF is a way forward, we would like some help understanding its definition; i.e., which entities become roots and which become leaves. In DMPlexBuildFromCellListParallel() https://gitlab.com/petsc/petsc/-/blob/753428fdb0644bc4cb7be6429ce8776c05405… the PointSF appears to contain roots for elements and vertices and leaves for owned vertices on the inter-process boundary. Is that correct? Thank-you, Cameron

3 10

How to confirm the performance of asynchronous computations
by Viet H.Q.H. 29 Jan '21

29 Jan '21

Hello Petsc developers and supporters, I would like to confirm the performance of asynchronous computations of inner product computation overlapping with matrix-vector multiplication computation by the below code. PetscLogDouble tt1,tt2; KSP ksp; //ierr = VecSet(c,one); ierr = VecSet(c,one); ierr = VecSet(u,one); ierr = VecSet(b,one); ierr = KSPCreate(PETSC_COMM_WORLD,&ksp); CHKERRQ(ierr); ierr = KSP_MatMult(ksp,A,x,Ax); CHKERRQ(ierr); ierr = PetscTime(&tt1);CHKERRQ(ierr); ierr = VecNormBegin(u,NORM_2,&norm1); ierr = PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)Ax)); ierr = KSP_MatMult(ksp,A,c,Ac); ierr = VecNormEnd(u,NORM_2,&norm1); ierr = PetscTime(&tt2);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD, "The time used for the asynchronous calculation: %f\n",tt2-tt1); CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"+ |u| = %g\n",(double) norm1); CHKERRQ(ierr); ierr = PetscTime(&tt1);CHKERRQ(ierr); ierr = VecNorm(b,NORM_2,&norm2); CHKERRQ(ierr); ierr = KSP_MatMult(ksp,A,c,Ac); ierr = PetscTime(&tt2);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD, "The time used for the synchronous calculation: %f\n",tt2-tt1); CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"+ |b| = %g\n",(double) norm2); CHKERRQ(ierr); On a cluster with two or four nodes, the asynchronous computation is always much slower than synchronous computation. The time used for the asynchronous calculation: 0.000203 + |u| = 100. The time used for the synchronous calculation: 0.000006 + |b| = 100. Are there any necessary settings on MPI or Petsc to gain performance of asynchronous computation? Thank you very much for anything you can provide. Sincerely, Viet.

7 17

Convert a 3D DMDA sub-vector to a natural 2D vector
by Sajid Ali 25 Jan '21

25 Jan '21

Hi PETSc-developers, For an application, I'd like to extract a 2D slice from a 3D DMDA vector, perform a `MatMult` on it (a discretized rotation on the 2D vector) and place the resulting vector back into the 3D vector. The approach I'd taken was to use `DMDACreatePatchIS` to create an IS that selects the slice (say this slice is [:,:,n_z] with no loss of generality). This IS is then used to extract the 2D slice via the `VecGetSubVector` utility. The issue that arises from this scheme is that the extracted 2D vector does not represent a flattened 2D vector, it is instead ordered as a DMDA vector (aka each rank having the portion of vector it owns in a column major ordering). How do I convert this DMDA ordered 2D subvector to a vector that represents a flattened 2D array? Since I can't predict how the extracted sub-vector While PETSc provides `DMDAGlobalToNatural` routines, those don't apply to extracted sub-vectors and extracting the 2D slice from a natural vector (filled by a scatter from the global vector) does not work either. On the other hand, converting the output of DMDACreatePatchIS to a natural IS via `AOPetscToApplicationIS` before extracting a subvector extracts the wrong indices from the 3D vector. Could someone tell me how I can achieve the desired goal (converting a 2D vector ordered as per a 3D DMDA grid onto a natural flattened index) ? Thanks in advance for the help! PS : Should it help, I'm attaching a slightly modified version of `src/dm/tests/ex53.c` to demonstrate the issue. The array I wish to obtain (the [:,:,2]) slice is : `[196., 194., 192., 190., 198., 197., 196., 195., 200., 200., 200.,200., 202., 203., 204., 205.]` (a row-major flattening of the 2D slice or alternatively a column-major flattening). What I obtain from the program (by creating the PatchIS and extracting the subvector) instead is : `[0 (A)] 196,198,194,197 [1 (B)] 200,202,200,203 [2 (C)] 192,196,190,195 [3 (D)] 200,204,200,205` (where each square bracket represents the rank and portion of the 2D vector being flattened via a column-major formatting with the layout shown below) |---------| | A | C | | B | D | |---------| Thank You, Sajid Ali (he/him) | PhD Candidate Applied Physics Northwestern University s-sajid-ali.github.io

2 2

Gathering distributed dense matrix only on rank 0
by Roland Richter 25 Jan '21

25 Jan '21

Hei, is there a way to gather a distributed dense matrix only on rank 0, not on the other ranks? My matrix processing/storage routines in my program currently are single-threaded and only operate on rank 0, and therefore I assume that I can ignore all other ranks. This should also save a bit time, I hope (MatCreateRedundantMatrix is currently the by far most expensive part of my program). When executing MatCreateRedundantMatrix(A, size, MPI_COMM_NULL, MAT_INITIAL_MATRIX, &B), usually all ranks get the full matrices, something I would like to avoid. Therefore, is there another way of handling this? Thanks! Regards, Roland Richter

2 1

DMPlex with dof in cells AND at nodes
by Thibault Bridel-Bertomeu 24 Jan '21

24 Jan '21

Dear all, I haven't been on this mailing list since december, so first of all, happy new year to everyone ! I am facing another challenge with the DMPlex structure : although I am doing cell-centered finite volume (using PetscFV), I also need to do some manipulations and store information that are vertex-based. The problem I have is when I run the code in parallel : when I visualize the nodal data using a vtk viewer, it shows NaN at all the vertices in close proximity to the MPI partition line. Here is how I proceed to work on the vertex-based data : First comes the creation of the DM that I'll be using for the cell-centered FVM call DMPlexCreateFromFile(comm, name, PETSC_TRUE, dm, ierr) call DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_TRUE, ierr) call DMPlexDistribute(dm, 1, PETSC_NULL_SF, dmDist, ierr) ... dm = dmDist ... call DMConstructGhostCells(dm, ..., ..., dmGhost, ierr) ... dm = dmGhost ... I set both boolean to PETSC_TRUE in the SetBasicAdjacency so even in MPI the neighbouring relationships are conserved for all strata. Is that true ? Like, does doing this actually enable to get all the transitive closure of a given vertex, even if part of it is in another MPI block ? Then I need a vertex-based DM, so I figure it's pretty much what's done in CreateMassMatrix in TS ex11.c. call DMClone(dm, dmNodal, ierr) call DMGetCoordinateSection(dm, coordSection, ierr); CHKERRA(ierr) call DMGetCoordinatesLocal(dm, coordinates, ierr); CHKERRA(ierr) call DMSetCoordinateSection(dmNodal, PETSC_DETERMINE, coordSection, ierr); call PetscSectionCreate(PETSC_COMM_WORLD, sectionNodal, ierr); CHKERRA(ierr) call DMPlexGetDepthStratum(dm, 0, vstart, vend, ierr); CHKERRA(ierr) call PetscSectionSetChart(sectionNodal, vstart, vend, ierr); CHKERRA(ierr) do inode = vstart, vend-1 ! 3 dof per node, ndim = 3 call PetscSectionSetDof(sectionNodal, inode, ndim, ierr); CHKERRA(ierr) end do call PetscSectionSetUp(sectionNodal, ierr); CHKERRA(ierr) call DMSetLocalSection(dmNodal, sectionNodal, ierr); CHKERRA(ierr) call PetscSectionDestroy(sectionNodal, ierr); CHKERRA(ierr) >From what I understood, this new dmNodal basically has all the characteristics of the ori dm, except it's primary layer is vertex-based with 3 DOF per vertex. Then I get a local vector and work on it : call DMGetLocalVector(dmNodal, nodalvec, ierr) call VecGetArrayF90(nodalvec, nodalvecarray, ierr) call DMPlexGetDepthStratum(dm, 0, vstart, vend, ierr) call DMPlexGetDepthStratum(dm, 2, fstart, fend, ierr) ! For each node ... do inode = vstart, vend-1 A = 0 ! matrix, ndim x ndim b = 0 ! vector, ndim call DMPlexGetTransitiveClosure(dmNodal, inode, PETSC_FALSE, vertexclosure, ierr); ! For each face attached to that node do iface = 1, size(vertexclosure) if ((vertexclosure(iface) < fstart) .or. (vertexclosure(iface) >= fend)) cycle A = A + ... b = b + ... end do nodalvecarray(1+(inode-vstart)*ndim:3+(inode-vstart)*ndim) = matmul(A, b) ... end do Then I want to visualize it : call VecRestoreArrayF90(nodalvec, nodalvecarray, ierr) call DMCreateGlobalVector(dmNodal, globalnodalvec, ierr) call DMLocalToGlobal(dmNodal, nodalvec, INSERT_VALUES, globalnodalvec, ierr) call PetscViewerCreate(PETSC_COMM_WORLD, vtkViewer, ierr); CHKERRA(ierr) call PetscViewerSetType(vtkViewer, PETSCVIEWERVTK, ierr); CHKERRA(ierr) call PetscViewerFileSetName(vtkViewer, "toto.vtk", ierr); CHKERRA(ierr) call VecView(globalnodalvec, vtkViewer, ierr); CHKERRA(ierr) call PetscViewerDestroy(vtkViewer, ierr); CHKERRA(ierr) In the toto.vtk file, I get vertex-based values as expected, and the whole shebang works in sequential. Now when I run it in parallel, I get vertex-based values inside the MPI blocks, but in the neighbourhood of the join between the MPI blocks, all the vertex-based values are NaN (see attachment). All in all, I think I do not understand fully the consequences of the DMPlexDistribute with useCone=true & useClosure=true ... but I cannot figure out where the NaN are coming from : is it because the vertices do not actually know all their closure through the MPI join, or is it some Global / Local error I made ? Thank you very much in advance for your support !! Thibault

2 5

MatCreateSchurComplementPmat to create a Schur complement
by Tang, Qi 21 Jan '21

21 Jan '21

Hi, When I use the MatCreateSchurComplementPmat to create a Schur complement for my own usage, I do not understand the logic when preuse=MAT_INITIAL_MATRIX. If all the matrixes are not empty, the first thing it did is 809: MatDestroy(Spmat); So that means I cannot pass in an empty matrix. It works if I first do MatCreate a dummy matrix and then pass in. Is that just some inconsistency or is there something I misunderstood? Qi

2 2

Program spends a large amount of time in PetscStrcmp and PetscCheckPointer - Why?
by Roland Richter 21 Jan '21

21 Jan '21

Hei, I noticed while benchmarking and profiling that my program spends ~2800 s in PetscStrcmp and ~2300 s in PetscCheckPointer (for a total runtime of 21000 s) when running it on four threads. I do not recall using those functions actively, especially for PetscStrcmp. The first function which is doing something useful (i.e. I know where it is in use) is zgemv_, with 574 s total runtime. Why are those functions called, and where are they necessary? Thanks! Regards, Roland Richter

3 2