December 2020 - 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

Error occurred when version change
by 朱昭霓 30 Dec '20

30 Dec '20

Dear Sir or Madam:  Hello, I am a new PETSc learner. There are some confusions when I using different versions of PETSc. For some reasons, I have to change my codes from the version 3.6.3 to 3.4.5. Newton methods in 3.6.3 version converged successfully with the time step size dt=100. While with the same dt=100 and almost the same codes, Newton methods in version 3.4.5 failed to converge. The debug test indicated that there may be something wrong during the line search process. I have tried to search the changes in different versions but failed. Could you please give me some suggestions? Thank you so much.  Best regards, Zhaoni Zhu

4 5

Files saved with PetscViewerBinaryOpen - Binary file format ?
by C B 29 Dec '20

29 Dec '20

Hello PETSc community ! I would like to read in another code a matrix / rhs / sol saved by PETSc. I am dealing with large matrices and for precision I would like to use the binary format. To save a matrix I am using the following lines (please let me know if I should make any changes), I am using an old version PETSc 2.1.5 for backward compatibility. ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,PETSC_BINARY_CREATE,&petsc_view);CHKERRQ(ierr); ierr = VecView(solut,petsc_view);CHKERRQ(ierr); ierr = PetscViewerDestroy(petsc_view);CHKERRQ(ierr); And to save a vector: ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, PETSC_BINARY_CREATE, &petsc_view);CHKERRQ(ierr); ierr = VecView(solut, petsc_view);CHKERRQ(ierr); ierr = PetscViewerDestroy(petsc_view);CHKERRQ(ierr); The files are saved without any problems, and now I need to write the code to read the binary files. Please, would anyone point me to the format of these binary files and/or a code snippet to take as an example ? Thank you very much in advance ! Cheers

2 2

Calculating adjoint of more than one cost function separately
by Salazar De Troya, Miguel 29 Dec '20

29 Dec '20

Hello, I am interested in calculating the gradients of an optimization problem with one goal and one constraint functions which need TSAdjoint for their adjoints. I’d like to call each of their adjoints in different calls, but it does not seem to be possible without making compromises. For instance, one could set TSCreateQuadratureTS() and TSSetCostGradients() with different quadratures (and their gradients) for each adjoint call (one at a time). This would evaluate the cost functions in the backwards run though, whereas one typically computes the cost functions in a different routine than the adjoint call (like in line searches evaluations) One could also set TSCreateQuadratureTS() with the goal and the constraint functions to be evaluated at the forward run (as typically done when computing the cost function). The problem would be that the adjoint call now requires two sets of gradients for TSSetCostGradients() and their adjoint are calculated together, costing twice if your routines for the cost and the constraint gradients are separated. The only solution I can think of is to set TSCreateQuadratureTS() with both the goal and constraint functions in the forward run. Then, in the adjoint calls, reset TSCreateQuadratureTS() with just the cost function I am interested in (either the goal or the constraint) and set just a single TSSetCostGradients(). Will this work? Are there better alternatives? Even if successful, there is the problem that the trajectory goes back to the beginning when we perform a TSAdjointSolve() call. Subsequent calls to TSAdjointSolve() (for instance for another cost function) are invalid because the trajectory is not set at the end of the simulation. One needs to call the forward problem to bring it back to the end. Is there a quick way to set the trajectory state to the last time step without having to run the forward problem? I am attaching an example to illustrate this issue. One can uncomment lines 120-122 to obtain the right value of the derivative. Thanks Miguel Miguel A. Salazar de Troya Postdoctoral Researcher, Lawrence Livermore National Laboratory B141 Rm: 1085-5 Ph: 1(925) 422-6411

2 4

TS tutorial ex11 in Fortran
by Thibault Bridel-Bertomeu 29 Dec '20

29 Dec '20

Dear all, Is there somewhere a version of the TS tutorial ex11.c in Fortran ? I am looking into building in F90 (let's say that it is an unavoidable constraint) an unstructured 3D solver of the Euler equations using the "new" features of PETSc - mostly DMPlex & PetscFV - but I think there are some interfaces missing and I find it hard to find workarounds in Fortran. I would be grateful if anyone could please give me some pointers ... 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

5 31

mapping from the original ordering to the PETSc/global ordering for the imported unstructured mesh
by Hui LIU 23 Dec '20

23 Dec '20

Dear the maintainer of the PETSc: I am implementing the CPU parallel finite element calculation of the unstructured grids by using the great PETSc platform. The unstructured grids are imported from other open access software, for example, Gmsh. Herein we call the node numbering for the mesh generated by the Gmsh as the original ordering, before distributing the mesh to each processor. After distributing this unstructured mesh to each process, is it possible to find the mapping between the original ordering (before distributing) and the PETSc/global ordering (after distributing)? If yes, how to extract or find it? If no, could you give me some suggestions for better applying the displacement constraints and load boundary conditions? Another question is how to distinguish the nodes subject to displacement boundary constraints and the nodes subject to the external force loads. Thank you very much. Looking forward to your reply. Best regards, Hui Liu Wuhan University

2 1

Support for full jacobianP in TSSetIJacobianP
by Salazar De Troya, Miguel 22 Dec '20

22 Dec '20

Hello, I am working on hooking up TSAdjoint with pyadjoint through the firedrake-ts interface (https://github.com/IvanYashchuk/firedrake-ts) I have done most of the implementation and now I am just testing for corner cases. One of them is when the design variable is multiplying the first derivative term. It would be the case of F(Udot,U,P,t) = G(U,P,t) in https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/Sensitivity/TSS… . I imagine that one could think of refactoring the “P” in the left hand side to the right hand side, but this is not trivial when “P” is a discontinuous field over the domain. I think it would be better to include the case of F(Udot,U,P,t) = G(U,P,t) in TSSetIJacobianP and I am volunteering to do it. Given the current implementation of TSAdjoint, is this something feasible? Thanks Miguel Miguel A. Salazar de Troya Postdoctoral Researcher, Lawrence Livermore National Laboratory B141 Rm: 1085-5 Ph: 1(925) 422-6411

2 10