March 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

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

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

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

DMPlex overlap
by Nicolas Barral 16 Apr '21

16 Apr '21

Hi all, First, I'm not sure I understand what the overlap parameter in DMPlexDistributeOverlap does. I tried the following: generate a small mesh on 1 rank with DMPlexCreateBoxMesh, then distribute it with DMPlexDistribute. At this point I have two nice partitions, with shared vertices and no overlapping cells. Then I call DMPlexDistributeOverlap with the overlap parameter set to 0 or 1, and get the same resulting plex in both cases. Why is that ? Second, I'm wondering what would be a good way to handle two overlaps and associated local vectors. In my adaptation code, the remeshing library requires a non-overlapping mesh, while the refinement criterion computation is based on hessian computations, which require a layer of overlap. What I can do is clone the dm before distributing the overlap, then manage two independent plex objects with their own local sections etc. and copy/trim local vectors manually. Is there a more automatic way to do this ? Thanks -- Nicolas

3 19

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

Newbie question: Something is wrong with this Slepc simple example for generalised eigenvalue problem
by dazza simplythebest 29 Mar '21

29 Mar '21

Dear All, I am seeking to use slepc/petsc to solve a generalised eigenvalue problem that arises in a hydrodynamic stability problem. The code is a parallelisation of an existing serial Fortran code. Before I get to grips with this target problem, I would of course like to get some relevant examples working. I have installed petsc/ slepc seemingly without any problem, and the provided slepc example fortran program ex1f.F, which solves a regular eigenvalue problem Ax = lambda x, seemed to compile and run correctly. I have now written a short program to instead solve the complex generalised problem Ax = lambda B x (see below) . This code compiles and runs w/out errors but for some reason hangs when calling EPSSolve - we enter EPSSolve but never leave. The matrices appear to be correctly assembled -all the values are correct in the Matview printout, so I am not quite sure where I have gone wrong, can anyone spot my mistake? ( Note that for the actual problem I wish to solve I have already written the code to construct the matrix, which distributes the rows across the processes and it is fully tested and working. Hence I want to specify the distribution of rows and not leave it up to a PETS_DECIDE .) I would be very grateful if someone can point out what is wrong with this small example code (see below), Many thanks, Dan. ! this program MUST be run with NGLOBAL = 10, MY_NO_ROWS = 5 ! and two MPI processes since row distribution is hard-baked into code ! program main #include <slepc/finclude/slepceps.h> use slepceps implicit none ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Declarations ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! ! Variables: ! A , B double complex operator matrices ! eps eigenproblem solver context Mat A,B EPS eps EPSType tname PetscReal tol, error PetscScalar kr, ki Vec xr, xi PetscInt NGLOBAL , MY_NO_ROWS, NL3, owner PetscInt nev, maxit, its, nconv PetscInt i,j,ic,jc PetscReal reala, imaga, realb, imagb, di, dj PetscScalar a_entry, b_entry PetscMPIInt rank PetscErrorCode ierr PetscInt,parameter :: zero = 0, one = 1, two = 2, three = 3 PetscInt M1, N1, mm1, nn1 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Beginning of program ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - call SlepcInitialize(PETSC_NULL_CHARACTER,ierr) call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr) ! make sure you set NGLOBAL = 10, MY_NO_ROWS = 5 and run with two processes NGLOBAL = 10 MY_NO_ROWS = 5 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Compute the operator matrices that define the eigensystem, Ax=kBx ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !!!!!!!!! Setup A matrix call MatCreate(PETSC_COMM_WORLD,A,ierr) call MatSetsizes(A,MY_NO_ROWS, MY_NO_ROWS ,NGLOBAL,NGLOBAL,ierr) call MatSetFromOptions(A,ierr) call MatGetSize(A,M1,N1,ierr) write(*,*)'Rank [',rank,']: global size of A is ',M1, N1 call MatGetLocalSize(A,mm1,nn1,ierr) write(*,*)'Rank [',rank,']: my local size of A is ',mm1, nn1 call MatMPIAIJSetPreallocation(A,three, PETSC_NULL_INTEGER,one, & & PETSC_NULL_INTEGER,ierr) !parallel (MPI) allocation !!!!!!!!! Setup B matrix call MatCreate(PETSC_COMM_WORLD,B,ierr) call MatSetsizes(B,MY_NO_ROWS, MY_NO_ROWS ,NGLOBAL,NGLOBAL,ierr) call MatSetFromOptions(B,ierr) call MatGetSize(B,M1,N1,ierr) write(*,*)'Rank [',rank,']: global size of B is ',M1, N1 call MatGetLocalSize(B,mm1,nn1,ierr) write(*,*)'Rank [',rank,']: my local size of B is ',mm1, nn1 call MatMPIAIJSetPreallocation(B,three, PETSC_NULL_INTEGER,one, & & PETSC_NULL_INTEGER,ierr) !parallel (MPI) allocation ! initalise call MatZeroEntries(A,ierr) call MatZeroEntries(B,ierr) ! Fill in values of A, B and assemble matrices ! Both matrices are tridiagonal with ! Aij = cmplx( (i-j)**2, (i+j)**2) ! Bij = cmplx( ij/i + j, (i/j)**2) ! (numbering from 1 ) do i = 1, NGLOBAL ! a rather crude way to distribute rows if (i < 6) owner = 0 if (i >= 6) owner = 1 if (rank /= owner) cycle do j = 1, NGLOBAL if ( abs(i-j) < 2 ) then write(*,*)rank,' : Setting ',i,j di = dble(i) ; dj = dble(j) reala = (di - dj)**2 ; imaga = (di + dj)**2 a_entry = dcmplx(reala, imaga) realb = (di*dj)/(di + dj) ; imagb = di**2/dj**2 b_entry = dcmplx(realb, imagb) ic = i -1 ; jc = j-1 ! convert to C indexing call MatSetValue(A, ic, jc, a_entry, ADD_VALUES,ierr) call MatSetValue(B, ic, jc, b_entry, ADD_VALUES,ierr) endif enddo enddo call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr) call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr) call MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr) call MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr) ! Check matrices write(*,*)'A matrix ... ' call MatView(A,PETSC_VIEWER_STDOUT_WORLD,ierr) write(*,*)'B matrix ... ' call MatView(B,PETSC_VIEWER_STDOUT_WORLD,ierr) call MatCreateVecs(A,PETSC_NULL_VEC,xr) call MatCreateVecs(A, PETSC_NULL_VEC,xi) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Create the eigensolver and display info ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! ** Create eigensolver context call EPSCreate(PETSC_COMM_WORLD,eps,ierr) ! ** Set operators.for general problem Ax = lambda B x call EPSSetOperators(eps,A, B, ierr) call EPSSetProblemType(eps,EPS_GNHEP,ierr) ! ** Set solver parameters at runtime call EPSSetFromOptions(eps,ierr) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Solve the eigensystem ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - write(*,*)'rank',rank, 'entering solver ...' call EPSSolve(eps,ierr) ! ** Free work space call EPSDestroy(eps,ierr) call MatDestroy(A,ierr) call MatDestroy(B,ierr) call VecDestroy(xr,ierr) call VecDestroy(xi,ierr) call SlepcFinalize(ierr) end

2 3

Code speedup after upgrading
by Mohammad Gohardoust 28 Mar '21

28 Mar '21

Hi, I am using a code which is based on petsc (and also parmetis). Recently I made the following changes and now the code is running about two times faster than before: - Upgraded Ubuntu 18.04 to 20.04 - Upgraded petsc 3.13.4 to 3.14.5 - This time I installed parmetis and metis directly via petsc by --download-parmetis --download-metis flags instead of installing them separately and using --with-parmetis-include=... and --with-parmetis-lib=... (the version of installed parmetis was 4.0.3 before) I was wondering what can possibly explain this speedup? Does anyone have any suggestions? Thanks, Mohammad

7 15

unsubscribe
by C B 28 Mar '21

28 Mar '21

unsubscribe

1 0