May 2021 - petsc-users - mailman3.cels.anl.gov

Efficient FFTShift-implementation for vectors/matrices
by Roland Richter 18 May '21

18 May '21

Dear all, I tried to implement the function fftshift from numpy (i.e. swap the half-spaces of all axis) for row vectors in a matrix by using the following code void fft_shift(Mat &fft_matrix) { PetscScalar *mat_ptr; MatDenseGetArray (fft_matrix, &mat_ptr); PetscInt r_0, r_1; MatGetOwnershipRange(fft_matrix, &r_0, &r_1); PetscInt local_row_num = r_1 - r_0; arma::cx_mat temp_mat(local_row_num, Ntime, arma::fill::zeros); for(int i = 0; i < Ntime; ++i) { const PetscInt row_shift = i * local_row_num; for(int j = 0; j < local_row_num; ++j) { const PetscInt cur_pos = j + row_shift; if(i < (int)(Ntime / 2)) temp_mat(j, i + int(Ntime / 2)) = *(mat_ptr + cur_pos); else temp_mat(j, i - int(Ntime / 2)) = *(mat_ptr + cur_pos); } } for(int i = 0; i < Ntime; ++i) { const PetscInt row_shift = i * local_row_num; for(int j = 0; j < local_row_num; ++j) { const PetscInt cur_pos = j + row_shift; *(mat_ptr + cur_pos) = temp_mat(j, i); } } MatDenseRestoreArray (fft_matrix, &mat_ptr); } but I do not like the approach of having a second matrix as temporary storage space. Are there more efficient approaches possible using PETSc-functions? Thanks! Regards, Roland Richter

2 1

reuse a real matrix for a second linear system with complex numbers
by feng wang 17 May '21

17 May '21

Dear All, I am solving a coupled system. One system is AX=B. A, X and B are all real numbers and it is solved with GMRES in petsc. Now I need to solve a second linear system, it can be represented as (A+i*w)*Z=C. i is the imaginary unit. Z and C are also complex numbers. So the Jacobian matrix of the second system is just A plus a diagonal contribution i*w. I would like solve the second system with GMRES, could petsc handle this? any comments are welcome. Thanks, Feng

2 8

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

Binary format in real vs. complex scalar type configurations
by Frederico Teixeira 12 May '21

12 May '21

Dear fellows, I hope this message finds you safe and well. I have a complex-valued matrix and its real/imaginary components in binary format . They were extracted from a solver that only works with "scalar-type=complex" configuration. I am getting weird results when I load them into a small test program that's configured with "scalar-type=real", but I believe this is expected. At the end of the day, I would like to have both real and imaginary components as real-valued matrices. Is it possible to do it? I want to test preconditioners that are tailored for this sort of problem. Regards, Frederico.

5 11

Open MPI library version error
by Daniel Stone 11 May '21

11 May '21

Hello, I have the following error: PETSc Error --- Open MPI library version FUJITSU MPI Library 4.0.0 (4.0.1fj4.0.0) does not match what PETSc was compiled with 4.0, aborting I have to configure and compile petsc in a cross compilation environment and run it in a different one. I'm trying to run src/snes/examples/tutorials/ex19 to confirm things are working. It looks like this error could be down to some confusion with string comparison ("4.0.0" =/= "4.0"), or it could be some nightmare misalignment between the two environments I have to use, and/or my efforts to configure PETSc. I'm trying to determine which. Is there some way to turn off this version checking, to see if it is just a string comparison error and maybe the rest of the program will run fine? Or could someone point me to where in the source code/configuration scripts these numbers are harvested and this error is generated, so I can try to get a better idea about what is going on? Thanks, Daniel Stone

4 4

MUMPS error
by Bhamidipati, Vikram 06 May '21

06 May '21

I configured a PetSc build to include mumps using '--download-mumps=1'. But it was not successful. I am attaching the relevant portion from configure.log. Could you advise me on what went wrong with the compilation and building of MUMPS library? Thanks, Vikram --------------------------------------------------- Vikram Bhamidipati Senior Research Engineer Computational Material Integrity Section Materials Engineering Department Mechanical Engineering Division Ph: (210) 522-2576

2 2

Error in installation of petsc4py
by arash rezaei 06 May '21

06 May '21

Hello, I want to install pets4py I use this command: python -m pip install [--user] petsc petsc4py

4 5

Equivalent function call for -mg_levels_ksp_type & -mg_levels_pc_type
by Milan Pelletier 06 May '21

06 May '21

Dear users, I would like to know what function call would allow me to achieve the same result as using the two flags -mg_levels_ksp_type richardson -mg_levels_pc_type jacobi and then PCSetFromOptions(). I would prefer not to rely on command-line options and rather do the setup "by hand". To achieve that, what should I add after these lines: KSPCreate(PETSC_COMM_WORLD, solver); KSPSetType(solver, KSPCG); KSPGetPC(solver, &pc); PCSetType(pc, PCGAMGAGG); <here I don't know what to set> Thanks for your help, Regards, Milan Pelletier

2 1

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