February 2012 - petsc-users - mailman3.cels.anl.gov

MatStructure flag in KSPSetOperators
by Kuhlemann, Verena 22 Feb '12

22 Feb '12

Hi, I am not sure that I understand the MatStructure flag in the call of KSPSetOperators correctly. If I use KSPSetOperators(ksp, A, P, DIFFERENT_NONZERO_PATTERN) and then solve two linear systems with the same matrix operators but different rhs is the preconditioner reassembled in every solve. I.e., if I call KSPSolve(ksp,b1,x1); KSPSolve(ksp,b2,x2); And if I use KSPSetOperators(ksp, A, P, SAME_PRECONDITIONER) the preconditioner is only setup one time and reused later. Or is the flag only important if I change the matrix operators in between. Thank you, Verena ________________________________ This e-mail message (including any attachments) is for the sole use of the intended recipient(s) and may contain confidential and privileged information. If the reader of this message is not the intended recipient, you are hereby notified that any dissemination, distribution or copying of this message (including any attachments) is strictly prohibited. If you have received this message in error, please contact the sender by reply e-mail message and destroy all copies of the original message (including attachments).

2 1

Add a matrix to MatNest twice
by Thomas Witkowski 21 Feb '12

21 Feb '12

Is it possible to have the same matrix of type MatMPIAIJ twice in a MatNest? Thomas

3 5

check for NULL pointer in VecCreateGhost
by Mohammad Mirzadeh 21 Feb '12

21 Feb '12

Hi, Is there an internal check to see if the 'ghost' pointer in VecCreateGhost() is NULL ? This can happen, for example, if the code is run in serial and the there are no ghost points (hence n_ghost = 0 and ghost = NULL(implementation dependent of course) ). AOXXXXToYYYY functions seem to crash if the corresponding pointer is NULL. Thanks, Mohammad

3 7

convergence problem in spherical coordinates
by Patrick Alken 21 Feb '12

21 Feb '12

Hello all, I am having great difficulty solving a 3D finite difference equation in spherical coordinates. I am solving the equation in a spherical shell region S(a,b), with the boundary conditions being that the function is 0 on both boundaries (r = a and r = b). I haven't imposed any boundary conditions on theta or phi which may be a reason its not converging. The phi boundary condition would be that the function is periodic in phi, but I don't know if this needs to be put into the matrix somehow? I nondimensionalized the equation before solving which helped a little bit. I've also scaled the matrix and RHS vectors by their maximum element to make all entries <= 1. I've tried both direct and iterative solvers. The direct solvers give a fairly accurate solution for small grids but seem unstable for larger grids. The PETSc iterative solvers converge for very small grids but for medium to large grids don't converge at all. When running with the command (for a small grid): *> ./main -ksp_converged_reason -ksp_monitor_true_residual -pc_type svd -pc_svd_monitor* I get the output: SVD: condition number 5.929088512946e+03, 0 of 1440 singular values are (nearly) zero SVD: smallest singular values: 2.742809162118e-04 2.807446554985e-04 1.548488288425e-03 1.852332719983e-03 2.782708934678e-03 SVD: largest singular values : 1.590835571953e+00 1.593368145758e+00 1.595771695877e+00 1.623691828398e+00 1.626235829632e+00 0 KSP preconditioned resid norm 2.154365616645e+03 true resid norm 8.365589263063e+00 ||r(i)||/||b|| 1.000000000000e+00 1 KSP preconditioned resid norm 4.832753933427e-10 true resid norm 4.587845792963e-12 ||r(i)||/||b|| 5.484187244549e-13 Linear solve converged due to CONVERGED_RTOL iterations 1 When plotting the output of this SVD solution, it looks pretty good, but svd isn't practical for larger grids. Using the command (on the same grid): *> ./main -ksp_converged_reason -ksp_monitor_true_residual -ksp_compute_eigenvalues -ksp_gmres_restart 1000 -pc_type none* The output is attached. There do not appear to be any 0 eigenvalues. The solution here is much less accurate than the SVD case since it didn't converge. I've also tried the -ksp_diagonal_scale -ksp_diagonal_scale_fix options which don't help very much. Any advice on how to trouble shoot this would be greatly appreciated. Some things I've checked already: 1) there aren't any 0 rows in the matrix 2) using direct solvers on very small grids seems to give decent solutions 3) there don't appear to be any 0 singular values or eigenvalues Perhaps the matrix has a null space, but I don't know how I would find out what the null space is? Is there a tutorial on how to do this? Thanks in advance!

2 2

partition an MPIAIJ matrix
by recrusader 20 Feb '12

20 Feb '12

Dear PETSc Developers, I have an MPIAIJ matrix, I want to use MatPartitioning functions to repartition the matrix. Before doing this, I need to do MatCreateMPIAdj to generate an Adj Mat. Could you give me some hints about how to get "i" and "j" ( http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCreateMP…) from existing MPIAIJ matrix? Thank you very much, Yujie

3 10

Help with MatTranspose
by Bisht, Gautam 20 Feb '12

20 Feb '12

Hi, I have a MPIAIJ matrix A as listed below: row 0: (36, 1) (37, 1) (42, 1) (43, 1) row 1: (37, 1) (38, 1) (43, 1) (44, 1) row 2: (38, 1) (39, 1) (44, 1) (45, 1) row 3: (39, 1) (40, 1) (45, 1) (46, 1) row 4: (40, 1) (41, 1) (46, 1) (47, 1) When I try to transpose the matrix, while running on 2 processors, the code crashes. Attached below is the F90 code. I would appreciate any help in figuring out why the code crashes. Thanks, -Gautam.

2 2

Re: [petsc-users] Starting point for Stokes fieldsplit
by Max Rudolph 20 Feb '12

20 Feb '12

Matt, Thank you for your help. > On Mon, Feb 20, 2012 at 2:05 PM, Max Rudolph <rudolph at berkeley.edu> wrote: > > > Hi Dave, > > Thanks for your help. > > > > Max > > > > Hey Max, > > > > Without knowing anything about the specific application related to > > your Stokes problem, or information about the mesh you are using, I > > have a couple of questions and suggestions which might help. > > > > > > The test case that I am working with is isoviscous convection, benchmark > > case 1a from Blankenbach 1989. > > > > 1) If A, is your stokes operator A = ( K,B ; B^T, 0 ), what is your > > precondition operator? > > Specifically, what is in the (2,2) slot in the precondioner? - i.e. > > what matrix are you you applying -stokes_fieldsplit_1_pc_type jacobi > > -stokes_fieldsplit_1_ksp_type preonly to? > > Is it the identity as in the SpeedUp notes? > > > > > > I think that this is the problem. The (2,2) slot in the LHS matrix is all > > zero (pressure does not appear in the continuity equation), so I think that > > the preconditioner is meaningless. I am still confused as to why this > > choice of preconditioner was suggested in the tutorial, and what is a > > better choice of preconditioner for this block? Should I be using one of > > the Schur complement methods instead of the additive or multiplicative > > field split? > > > > Its not suggested, it is demonstrated. Its the first logical choice, since > Jacobi gives the identity for a 0 block (see > http://www.jstor.org/pss/2158202) Its > not meaningless. All the better preconditioners involve either a Schur > complement (also shown in the tutorial), or an auxiliary operator which is > more > difficult to setup and thus not shown. Thank you for clarifying this. > > > 2) This choice > > -stokes_fieldsplit_0_pc_type ml -stokes_fieldsplit_0_ksp_type preonly > > may simply not be a very effective and degrade the performance of the > > outer solver. > > I'd make the solver for the operator in the (1,1) slot much stronger, > > for example > > -stokes_fieldsplit_0_ksp_type gmres > > -stokes_fieldsplit_0_ksp_rtol 1.0e-4 > > -stokes_fieldsplit_0_mg_levels_ksp_type gmres > > -stokes_fieldsplit_0_mg_levels_pc_type bjacobi > > -stokes_fieldsplit_0_mg_levels_ksp_max_it 4 > > > > Add a monitor on this solver (-stokes_fieldsplit_0_ksp_XXX) to see how > > ML is doing. > > > > 3) Using -stokes_pc_fieldsplit_type MULTIPLICATIVE should reduce the > > number of outer iterations by a factor of two, but it will use more > > memory. > > > > 4) You should use a flexible Krylov method on the outer most solve > > (-stokes_ksp_XXX) as the preconditioner is varying between each outer > > iteration. Use -stokes_ksp_type fgmres or -stokes_ksp_type gcr > > > > > > Thanks for pointing this out. I made that change. > > > > 5) Depending on how the physical problem is scaled > > (non-dimensionalised), the size of the residuals associated with the > > momentum and continuity equation make be quite different. You are > > currently use the entire residual from (u,p) to determine when to stop > > iterating. You might want to consider writing a monitor which examines > > the these residuals independently. > > > > > > I think that I have scaled the problem correctly. I (slowly) obtain a > > sufficiently accurate solution using as options only: > > -stokes_ksp_atol 1e-5 -stokes_ksp_rtol 1e-5 > > -stokes_ksp_monitor_true_residual -stokes_ksp_norm_type UNPRECONDITIONED > > > > How do you know the problem is scaled correctly? Have you looked at norms > of the residuals for the two systems > Thanks, > > Matt > > > > Cheers, > > Dave Yes, here are the norms computed for the P, X, and Y components, following the last residual that ksp_monitor_true_residual returned: 383 KSP unpreconditioned resid norm 1.121628211019e-03 true resid norm 1.121628224178e-03 ||r(i)||/||b|| 9.626787321554e-10 P, X, Y residual norms 5.340336e-02, 4.463404e-02, 2.509621e-02

3 2

Unhandled exception at 0x00000001408cca07 ...
by TAY wee-beng 20 Feb '12

20 Feb '12

Hi, I'm runing my CFD code in windows visual studio 2008 with ifort 64bit mpich2 64bit. I managed to build the PETSc library after doing some modifications - [petsc-maint #105754] Error : Cannot determine Fortran module include flag. There is no error building my CFD code. However, when running my code, I got the error at: call PetscInitialize(PETSC_NULL_CHARACTER,ierr) It jumps to : *ierr = PetscMemzero(name,256); if (*ierr) return; May I know what's wrong? It was working when I run it in compaq visual fortran under windows xp 32bit. -- Yours sincerely, TAY wee-beng

2 2

PETSc Matrix Partitioning and SPMVM
by Bibrak Qamar 20 Feb '12

20 Feb '12

Hello all, The way PETSc MatCreateMPIAIJ distributes an N*N square matrix is row wise. And internally every processor stores its local matrix into two sub matrices one Diagonal and other part is off-Diagonal. (more here --> http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCreateMP… ). So when PETSc MatMul (v = A.x) is called it basically tries to hide communication of vector (x) by overlapping computation. As for as I understand every process first initiates a kind of non-Blocking Bcast of vector (x) and continues to do MVM of diagonal submatrix and then waits till communication is done and finally does MVM for the off-Diagonal matrix. My question is this (since I am new) what was the historical reason PETSc opted for this Diagonal and off-Diagonal storage? And what if I want to change the overlapping strategy of MVM by lets say introducing a ring based communication of vector (x), then I have to partition the local matrix into not two sub-matrices but P sub-matrices (here P = number of processors). Does PETSc provide this facility or one has to go from scratch to implement different techniques to store local matrix? Thanks Bibrak

2 1

Starting point for Stokes fieldsplit
by Max Rudolph 20 Feb '12

20 Feb '12

I am solving a 2D Stokes flow problem using a finite volume discretization, velocity (dof 0,1)-pressure (dof 2) formulation. Until now I have mostly used MUMPS. I am trying to use the PCFieldSplit interface now, with the goal of trying out the multigrid functionality provided through ML. I looked at these (http://www.mcs.anl.gov/petsc/documentation/tutorials/Speedup10.pdf) talk slides for some guidance on command line options and tried. When assembling the linear system, I obtain the global degree-of-freedom indices using DMDAGetGlobalIndices. The solver does not appear to be converging and I was wondering if someone could explain to me why this might be happening. Thanks for your help. Max petscmpiexec -n $1 $2 $3 \ -stokes_pc_fieldsplit_0_fields 0,1 -stokes_pc_fieldsplit_1_fields 2 \ -stokes_pc_type fieldsplit -stokes_pc_fieldsplit_type additive \ -stokes_fieldsplit_0_pc_type ml -stokes_fieldsplit_0_ksp_type preonly \ -stokes_fieldsplit_1_pc_type jacobi -stokes_fieldsplit_1_ksp_type preonly \ -stokes_ksp_view \ -stokes_ksp_monitor_true_residual \ Residual norms for stokes_ solve. 0 KSP preconditioned resid norm 2.156200561011e-07 true resid norm 1.165111661413e+06 ||r(i)||/||b|| 1.000000000000e+00 1 KSP preconditioned resid norm 2.156200561011e-07 true resid norm 1.165111661413e+06 ||r(i)||/||b|| 1.000000000000e+00 2 KSP preconditioned resid norm 5.176186579848e-08 true resid norm 5.747484453195e+06 ||r(i)||/||b|| 4.932990239085e+00 3 KSP preconditioned resid norm 5.054067588022e-08 true resid norm 5.763708451168e+06 ||r(i)||/||b|| 4.946915082953e+00 4 KSP preconditioned resid norm 3.556841873413e-08 true resid norm 4.649778784249e+06 ||r(i)||/||b|| 3.990843915003e+00 5 KSP preconditioned resid norm 3.177972840516e-08 true resid norm 4.677248322326e+06 ||r(i)||/||b|| 4.014420657890e+00 6 KSP preconditioned resid norm 3.100188857346e-08 true resid norm 4.857618195959e+06 ||r(i)||/||b|| 4.169229745814e+00 7 KSP preconditioned resid norm 3.045495907075e-08 true resid norm 5.030091724356e+06 ||r(i)||/||b|| 4.317261504580e+00 8 KSP preconditioned resid norm 2.993896937859e-08 true resid norm 5.213745290794e+06 ||r(i)||/||b|| 4.474888942808e+00 9 KSP preconditioned resid norm 2.944838631679e-08 true resid norm 5.403745734510e+06 ||r(i)||/||b|| 4.637963822245e+00 10 KSP preconditioned resid norm 2.898115557667e-08 true resid norm 5.596205405394e+06 ||r(i)||/||b|| 4.803149423985e+00 11 KSP preconditioned resid norm 2.853548106431e-08 true resid norm 5.788329048801e+06 ||r(i)||/||b|| 4.968046617765e+00 12 KSP preconditioned resid norm 2.810975464206e-08 true resid norm 5.978285004145e+06 ||r(i)||/||b|| 5.131083313418e+00 13 KSP preconditioned resid norm 2.770253127208e-08 true resid norm 6.164551408716e+06 ||r(i)||/||b|| 5.290953316217e+00 14 KSP preconditioned resid norm 2.731250834421e-08 true resid norm 6.346469882864e+06 ||r(i)||/||b|| 5.447091547575e+00 15 KSP preconditioned resid norm 2.693850812065e-08 true resid norm 6.523343466838e+06 ||r(i)||/||b|| 5.598899816114e+00

2 1