[petsc-users] Solving multiple linear systems with similar matrix sequentially
Hi, Contingency analysis of power systems requires to solve a sequence a (non-)linear systems in which the Jacobian/matrix changes by a rank-2/rank-1 update. Does the PETSc interface to GAMG allow to freeze the AMG hierarchy for a given matrix (base case) and to reuse this hierarchy in solving a list of slightly perturbed linear problems (list of contingencies)? Thanks, Domenico.
Re: [petsc-users] Solving multiple linear systems with similar matrix sequentially
On Feb 23, 2015, at 5:18 AM, domenico lahaye domenico_lah...@yahoo.com wrote: Hi, Contingency analysis of power systems requires to solve a sequence a (non-)linear systems in which the Jacobian/matrix changes by a rank-2/rank-1 update. Does the PETSc interface to GAMG allow to freeze the AMG hierarchy for a given matrix (base case) and to reuse this hierarchy in solving a list of slightly perturbed linear problems (list of contingencies)? You can do two things. Freeze the current preconditioner and solve a sequence of (new with slightly different numerical values) linear systems using that same precondition. Use KSPSetReusePreconditioner() in PETSc 3.5.x to freeze/unfreeze the preconditioner Freeze the hierarchy and coarse grid interpolations of GAMG but compute the new coarse grid operators RAP for each new linear system (this is a much cheaper operation). Use PCGAMGSetReuseInterpolation() to freeze/unfreeze the hierarchy. Barry Thanks, Domenico.
Re: [petsc-users] DMNetworkGetEdgeRange() in parallel
On Sun, Feb 22, 2015 at 3:59 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Thanks. Once I obtain that Index Set with the routine DMPlexCreateCellNumbering() (I assume that the edges in DMNetwork correspond to cells in DMPlex) can I use it to partition a vector with as many components as edges I have in my network? I do not completely understand the question. If you want a partition of the edges, you can use DMPlexCreatePartition() and its friend DMPlexDistribute(). What are you trying to do? Matt Thanks Miguel On Sun, Feb 22, 2015 at 12:15 PM, Matthew Knepley knep...@gmail.com wrote: On Sun, Feb 22, 2015 at 11:01 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Hi I noticed that the routine DMNetworkGetEdgeRange() returns the local indices for the edge range. Is there any way to obtain the global indices? So if my network has 10 edges, the processor 1 has the 0-4 edges and the processor 2, the 5-9 edges, how can I obtain this information? One of the points of DMPlex is we do not require a global numbering. Everything is numbered locally, and the PetscSF maps local numbers to local numbers in order to determine ownership. If you want to create a global numbering for some reason, you can using DMPlexCreatePointNumbering(). There are also cell and vertex versions that we use for output, so you could do it just for edges as well. Thanks, Matt Thanks Miguel -- *Miguel Angel Salazar de Troya* Graduate Research Assistant Department of Mechanical Science and Engineering University of Illinois at Urbana-Champaign (217) 550-2360 salaz...@illinois.edu -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener -- *Miguel Angel Salazar de Troya* Graduate Research Assistant Department of Mechanical Science and Engineering University of Illinois at Urbana-Champaign (217) 550-2360 salaz...@illinois.edu -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
Re: [petsc-users] DMNetworkGetEdgeRange() in parallel
Thanks, that will help me. Now what I would like to have is the following: if I have two processors and ten edges, the partitioning results in the first processor having the edges 0-4 and the second processor, the edges 5-9. I also have a global vector with as many components as edges, 10. How can I partition it so the first processor also has the 0-4 components and the second, the 5-9 components of the vector? Miguel On Feb 23, 2015 8:08 AM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: Miguel, One possible way is to store the global numbering of any edge/vertex in the component attached to it. Once the mesh gets partitioned, the components are also distributed so you can easily retrieve the global number of any edge/vertex by accessing its component. This is what is done in the DMNetwork example pf.c although the global numbering is not used for anything. Shri From: Matthew Knepley knep...@gmail.com Date: Mon, 23 Feb 2015 07:54:34 -0600 To: Miguel Angel Salazar de Troya salazardetr...@gmail.com Cc: petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov Subject: Re: [petsc-users] DMNetworkGetEdgeRange() in parallel On Sun, Feb 22, 2015 at 3:59 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Thanks. Once I obtain that Index Set with the routine DMPlexCreateCellNumbering() (I assume that the edges in DMNetwork correspond to cells in DMPlex) can I use it to partition a vector with as many components as edges I have in my network? I do not completely understand the question. If you want a partition of the edges, you can use DMPlexCreatePartition() and its friend DMPlexDistribute(). What are you trying to do? Matt Thanks Miguel On Sun, Feb 22, 2015 at 12:15 PM, Matthew Knepley knep...@gmail.com wrote: On Sun, Feb 22, 2015 at 11:01 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Hi I noticed that the routine DMNetworkGetEdgeRange() returns the local indices for the edge range. Is there any way to obtain the global indices? So if my network has 10 edges, the processor 1 has the 0-4 edges and the processor 2, the 5-9 edges, how can I obtain this information? One of the points of DMPlex is we do not require a global numbering. Everything is numbered locally, and the PetscSF maps local numbers to local numbers in order to determine ownership. If you want to create a global numbering for some reason, you can using DMPlexCreatePointNumbering(). There are also cell and vertex versions that we use for output, so you could do it just for edges as well. Thanks, Matt Thanks Miguel -- *Miguel Angel Salazar de Troya* Graduate Research Assistant Department of Mechanical Science and Engineering University of Illinois at Urbana-Champaign (217) 550-2360 salaz...@illinois.edu -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener -- *Miguel Angel Salazar de Troya* Graduate Research Assistant Department of Mechanical Science and Engineering University of Illinois at Urbana-Champaign (217) 550-2360 salaz...@illinois.edu -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
Re: [petsc-users] MatMPIAIJSetPreallocation: nnz cannot be greater than row length
Hi Barry, hi Dmitry, I set the matrix to BAIJ and back to AIJ, and the code got a bit further. But I now run into the error pasted below (Note that I'm now using --with-debugging=1): PETSC ERROR: - Error Message -- PETSC ERROR: Petsc has generated inconsistent data PETSC ERROR: MPIAIJ Matrix was assembled but is missing garray PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting. PETSC ERROR: Petsc Release Version 3.5.2, Sep, 08, 2014 PETSC ERROR: ./example-dbg on a arch-linux2-c-debug named david-Lenovo by dknez Mon Feb 23 08:05:44 2015 PETSC ERROR: Configure options --with-shared-libraries=1 --with-debugging=1 --download-suitesparse=1 --download-parmetis=1 --download-blacs=1 --download-scalapack=1 --download-mumps=1 --download-metis --download-superlu_dist --prefix=/home/dknez/software/libmesh_install/dbg_real/petsc --download-hypre PETSC ERROR: #1 MatCreateColmap_MPIAIJ_Private() line 361 in /home/dknez/software/petsc-3.5.2/src/mat/impls/aij/mpi/mpiaij.c PETSC ERROR: #2 MatSetValues_MPIAIJ() line 538 in /home/dknez/software/petsc-3.5.2/src/mat/impls/aij/mpi/mpiaij.c PETSC ERROR: #3 MatSetValues() line 1136 in /home/dknez/software/petsc-3.5.2/src/mat/interface/matrix.c PETSC ERROR: #4 add_matrix() line 765 in /home/dknez/software/libmesh-src/src/numerics/petsc_matrix.C -- This occurs when I try to set some entries of the matrix. Do you have any suggestions on how I can resolve this? Thanks! David On Sun, Feb 22, 2015 at 10:22 PM, Dmitry Karpeyev dkarp...@gmail.com wrote: On Sun Feb 22 2015 at 9:15:22 PM Barry Smith bsm...@mcs.anl.gov wrote: On Feb 22, 2015, at 9:09 PM, David Knezevic david.kneze...@akselos.com wrote: Hi Dmitry, Thanks for the suggestion. I tried MatSetType(mat,MATMPIAIJ) followed by MatXAIJSetPreallocation(...), but unfortunately this still gives me the same error as before: nnz cannot be greater than row length: local row 168 value 24 rowlength 0. I gather that the idea here is that MatSetType builds a new matrix object, and then I should be able to pre-allocate for that new matrix however I like, right? Was I supposed to clear the matrix object somehow before calling MatSetType? (I didn't do any sort of clear operation.) If the type doesn't change then MatSetType() won't do anything. You can try setting the type to BAIJ and then setting the type back to AIJ. This may/should clear out the matrix. Ah, yes. If the type is the same as before it does quit early, but changing the type and then back will clear out and rebuild the matrix. We need something like MatReset() to do the equivalent thing. As I said earlier, I'll make a dbg PETSc build, so hopefully that will help shed some light on what's going wrong for me. I think it's always a good idea to have a dbg build of PETSc when you doing things like these. Dmitry. Don't bother, what I suggested won't work. Barry Thanks, David On Sun, Feb 22, 2015 at 6:02 PM, Dmitry Karpeyev dkarp...@gmail.com wrote: David, It might be easier to just rebuild the whole matrix from scratch: you would in effect be doing all that with disassembling and resetting the preallocation. MatSetType(mat,MATMPIAIJ) or PetscObjectGetType((PetscObject)mat,type); MatSetType(mat,type); followed by MatXAIJSetPreallocation(...); should do. Dmitry. On Sun Feb 22 2015 at 4:45:46 PM Barry Smith bsm...@mcs.anl.gov wrote: Do not call for SeqAIJ matrix. Do not call before the first time you have preallocated and put entries in the matrix and done the MatAssemblyBegin/End() If it still crashes you'll need to try the debugger Barry On Feb 22, 2015, at 4:09 PM, David Knezevic david.kneze...@akselos.com wrote: Hi Barry, Thanks for your help, much appreciated. I added a prototype for MatDisAssemble_MPIAIJ: PETSC_INTERN PetscErrorCode MatDisAssemble_MPIAIJ(Mat); and I added a call to MatDisAssemble_MPIAIJ before MatMPIAIJSetPreallocation. However, I get a segfault on the call to MatDisAssemble_MPIAIJ. The segfault occurs in both serial and parallel. FYI, I'm using Petsc 3.5.2, and I'm not using a non-debug build (though I could rebuild PETSc in debug mode if you think that would help figure out what's happening here). Thanks, David On Sun, Feb 22, 2015 at 1:13 PM, Barry Smith bsm...@mcs.anl.gov wrote: David, This is an obscure little feature of MatMPIAIJ, each time you change the sparsity pattern before you call the MatMPIAIJSetPreallocation you need to call MatDisAssemble_MPIAIJ(Mat mat).This is a private PETSc function so you need to provide your own prototype for it above the function you use it in. Let us know if this resolves the problem. Barry
Re: [petsc-users] FE discretization in DMPlex
On Sun, Feb 22, 2015 at 7:05 PM, Justin Chang jychan...@gmail.com wrote:
Hi Matt,
Bringing this thread back from the dead.
1) Have you had the chance to implement things like RT and DG in DMPlex?
No, there has not been much call. Its on my stack of things to do.
2) Are examples/tests that illustrate how to do dualspaces?
I have commented the dual space routines now. Basically, you just create a
set of functionals,
and in this world a functional is just a quadrature rule. If you have a
hard time understanding
something, feel free to mail.
3) Or quantities like cell size h, jump, average?
The first thing is to declare the adjacency correctly, so that you get
neighboring cells. Once you
have that all these are simple local calculations.
Why don't we start with RT0 since it is the simplest to think about.
Thanks,
Matt
I was originally trying to implement DG and RT0 in FEniCS but I am having
lots of trouble getting the FEniCS code to scale on our university's
clusters, so that's why I want to attempt going back to PETSc's DMPlex to
do strong scaling studies.
Thanks,
Justin
On Sat, Sep 6, 2014 at 3:58 AM, Matthew Knepley knep...@gmail.com wrote:
On Fri, Sep 5, 2014 at 10:55 PM, Justin Chang jychan...@gmail.com
wrote:
Hi all,
So I understand how the FEM code works in the DMPlex examples (ex12 and
62). Pardon me if this is a silly question.
1) If I wanted to solve either the poisson or stokes using the
discontinuous Galerkin method, is there a way to do this with the built-in
DMPlex/FEM functions? Basically each cell/element has its own set of
degrees of freedom, and jump/average operations would be needed to
connect the dofs across element interfaces.
2) Or how about using something like Raviart-Thomas spaces (we'll say
lowest order for simplicity). Where the velocity dofs are not nodal
quantities, instead they are denoted by edge fluxes (or face fluxes for
tetrahedrals). Pressure would be piecewise constant.
Intuitively these should be doable if I were to write my own
DMPlex/PetscSection code, but I was wondering if the above two
discretizations are achievable in the way ex12 and ex62 are.
Lets do RT first since its easier. The primal space is
P_K = Poly_{q--1}(K) + x Poly_{q-1}(K)
so at lowest order its just Poly_1. The dual space is moments of the
normal component
of velocity on the edges. So you would write a dual space where the
functionals integrated
the normal component. This is the tricky part:
http://www.math.chalmers.se/~logg/pub/papers/KirbyLoggEtAl2010a.pdf
DG is just a generalization of this kind of thing where you need to a)
have some geometric
quantities available to the pointwise functions (like h), and also some
field quantities (like
the jump and average).
I understand exactly how I want to do the RT, BDM, BDMF, and NED
elements, and those
will be in soon. I think DG is fairly messy and am not completely sure
what I want here.
Matt
Thanks,
Justin
--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
Re: [petsc-users] DMNetworkGetEdgeRange() in parallel
Miguel, One possible way is to store the global numbering of any edge/vertex in the component attached to it. Once the mesh gets partitioned, the components are also distributed so you can easily retrieve the global number of any edge/vertex by accessing its component. This is what is done in the DMNetwork example pf.c although the global numbering is not used for anything. Shri From: Matthew Knepley knep...@gmail.commailto:knep...@gmail.com Date: Mon, 23 Feb 2015 07:54:34 -0600 To: Miguel Angel Salazar de Troya salazardetr...@gmail.commailto:salazardetr...@gmail.com Cc: petsc-users@mcs.anl.govmailto:petsc-users@mcs.anl.gov petsc-users@mcs.anl.govmailto:petsc-users@mcs.anl.gov Subject: Re: [petsc-users] DMNetworkGetEdgeRange() in parallel On Sun, Feb 22, 2015 at 3:59 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.commailto:salazardetr...@gmail.com wrote: Thanks. Once I obtain that Index Set with the routine DMPlexCreateCellNumbering() (I assume that the edges in DMNetwork correspond to cells in DMPlex) can I use it to partition a vector with as many components as edges I have in my network? I do not completely understand the question. If you want a partition of the edges, you can use DMPlexCreatePartition() and its friend DMPlexDistribute(). What are you trying to do? Matt Thanks Miguel On Sun, Feb 22, 2015 at 12:15 PM, Matthew Knepley knep...@gmail.commailto:knep...@gmail.com wrote: On Sun, Feb 22, 2015 at 11:01 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.commailto:salazardetr...@gmail.com wrote: Hi I noticed that the routine DMNetworkGetEdgeRange() returns the local indices for the edge range. Is there any way to obtain the global indices? So if my network has 10 edges, the processor 1 has the 0-4 edges and the processor 2, the 5-9 edges, how can I obtain this information? One of the points of DMPlex is we do not require a global numbering. Everything is numbered locally, and the PetscSF maps local numbers to local numbers in order to determine ownership. If you want to create a global numbering for some reason, you can using DMPlexCreatePointNumbering(). There are also cell and vertex versions that we use for output, so you could do it just for edges as well. Thanks, Matt Thanks Miguel -- Miguel Angel Salazar de Troya Graduate Research Assistant Department of Mechanical Science and Engineering University of Illinois at Urbana-Champaign (217) 550-2360tel:%28217%29%20550-2360 salaz...@illinois.edumailto:salaz...@illinois.edu -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener -- Miguel Angel Salazar de Troya Graduate Research Assistant Department of Mechanical Science and Engineering University of Illinois at Urbana-Champaign (217) 550-2360tel:%28217%29%20550-2360 salaz...@illinois.edumailto:salaz...@illinois.edu -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
Re: [petsc-users] solving multiple linear systems with same matrix (sequentially, not simultaneously)
Thanks Jed. So the nullspace would be initialized with an array of 3 petsc vectors: (u,v)= (0,1), (1,0), and (-y,x), correct? And also to be sure -- this is usefuil only for multigrid preconditioners, yes? Thanks Dan On Mon, Feb 23, 2015 at 2:11 AM, Jed Brown j...@jedbrown.org wrote: Barry Smith bsm...@mcs.anl.gov writes: I will try gamg as I know it can reduce the number of CG iterations required. (I'm guessing you mean algebraic, not geometric?) By default GAMG is an algebraic multigrid preconditioner. Look at its documentation at http://www.mcs.anl.gov/petsc/petsc-master/docs/index.html it will be a bit better than the older documentation. The documentation for GAMG is still pretty thin so feel free to ask questions. You might want to call MatSetBlockSize and MatSetNearNullSpace to define both translation and rotation modes. This is most relevant for large ice shelves. -- Daniel Goldberg, PhD Lecturer in Glaciology School of Geosciences, University of Edinburgh Geography Building, Drummond Street, Edinburgh EH8 9XP em: D dgold...@mit.eduan.goldb...@ed.ac.uk web: http://ocean.mit.edu/~dgoldberg
Re: [petsc-users] DMNetworkGetEdgeRange() in parallel
On Mon, Feb 23, 2015 at 8:42 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Thanks, that will help me. Now what I would like to have is the following: if I have two processors and ten edges, the partitioning results in the first processor having the edges 0-4 and the second processor, the edges 5-9. I also have a global vector with as many components as edges, 10. How can I partition it so the first processor also has the 0-4 components and the second, the 5-9 components of the vector? I think it would help to know what you want to accomplish. This is how you are proposing to do it.' If you just want to put data on edges, DMNetwork has a facility for that already. Thanks, Matt Miguel On Feb 23, 2015 8:08 AM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: Miguel, One possible way is to store the global numbering of any edge/vertex in the component attached to it. Once the mesh gets partitioned, the components are also distributed so you can easily retrieve the global number of any edge/vertex by accessing its component. This is what is done in the DMNetwork example pf.c although the global numbering is not used for anything. Shri From: Matthew Knepley knep...@gmail.com Date: Mon, 23 Feb 2015 07:54:34 -0600 To: Miguel Angel Salazar de Troya salazardetr...@gmail.com Cc: petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov Subject: Re: [petsc-users] DMNetworkGetEdgeRange() in parallel On Sun, Feb 22, 2015 at 3:59 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Thanks. Once I obtain that Index Set with the routine DMPlexCreateCellNumbering() (I assume that the edges in DMNetwork correspond to cells in DMPlex) can I use it to partition a vector with as many components as edges I have in my network? I do not completely understand the question. If you want a partition of the edges, you can use DMPlexCreatePartition() and its friend DMPlexDistribute(). What are you trying to do? Matt Thanks Miguel On Sun, Feb 22, 2015 at 12:15 PM, Matthew Knepley knep...@gmail.com wrote: On Sun, Feb 22, 2015 at 11:01 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Hi I noticed that the routine DMNetworkGetEdgeRange() returns the local indices for the edge range. Is there any way to obtain the global indices? So if my network has 10 edges, the processor 1 has the 0-4 edges and the processor 2, the 5-9 edges, how can I obtain this information? One of the points of DMPlex is we do not require a global numbering. Everything is numbered locally, and the PetscSF maps local numbers to local numbers in order to determine ownership. If you want to create a global numbering for some reason, you can using DMPlexCreatePointNumbering(). There are also cell and vertex versions that we use for output, so you could do it just for edges as well. Thanks, Matt Thanks Miguel -- *Miguel Angel Salazar de Troya* Graduate Research Assistant Department of Mechanical Science and Engineering University of Illinois at Urbana-Champaign (217) 550-2360 salaz...@illinois.edu -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener -- *Miguel Angel Salazar de Troya* Graduate Research Assistant Department of Mechanical Science and Engineering University of Illinois at Urbana-Champaign (217) 550-2360 salaz...@illinois.edu -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
Re: [petsc-users] MatMPIAIJSetPreallocation: nnz cannot be greater than row length
Hi Dmitry, Thanks very much for testing out the example. examples/systems_of_equations/ex8 works fine for me in serial, but it fails for me if I run with more than 1 MPI process. Can you try it with, say, 2 or 4 MPI processes? If we need access to MatReset in libMesh to get this to work, I'll be happy to work on a libMesh pull request for that. David -- David J. Knezevic | CTO Akselos | 17 Bay State Road | Boston, MA | 02215 Phone (office): +1-857-265-2238 Phone (mobile): +1-617-599-4755 Web: http://www.akselos.com On Mon, Feb 23, 2015 at 10:08 AM, Dmitry Karpeyev karp...@mcs.anl.gov wrote: David, What code are you running when you encounter this error? I'm trying to reproduce it and I tried examples/systems_of_equations/ex8, but it ran for me, ostensibly to completion. I have a small PETSc pull request that implements MatReset(), which passes a small PETSc test, but libMesh needs some work to be able to build against petsc/master because of some recent changes to PETSc. Dmitry. On Mon Feb 23 2015 at 7:17:06 AM David Knezevic david.kneze...@akselos.com wrote: Hi Barry, hi Dmitry, I set the matrix to BAIJ and back to AIJ, and the code got a bit further. But I now run into the error pasted below (Note that I'm now using --with-debugging=1): PETSC ERROR: - Error Message -- PETSC ERROR: Petsc has generated inconsistent data PETSC ERROR: MPIAIJ Matrix was assembled but is missing garray PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting. PETSC ERROR: Petsc Release Version 3.5.2, Sep, 08, 2014 PETSC ERROR: ./example-dbg on a arch-linux2-c-debug named david-Lenovo by dknez Mon Feb 23 08:05:44 2015 PETSC ERROR: Configure options --with-shared-libraries=1 --with-debugging=1 --download-suitesparse=1 --download-parmetis=1 --download-blacs=1 --download-scalapack=1 --download-mumps=1 --download-metis --download-superlu_dist --prefix=/home/dknez/software/libmesh_install/dbg_real/petsc --download-hypre PETSC ERROR: #1 MatCreateColmap_MPIAIJ_Private() line 361 in /home/dknez/software/petsc-3.5.2/src/mat/impls/aij/mpi/mpiaij.c PETSC ERROR: #2 MatSetValues_MPIAIJ() line 538 in /home/dknez/software/petsc-3.5.2/src/mat/impls/aij/mpi/mpiaij.c PETSC ERROR: #3 MatSetValues() line 1136 in /home/dknez/software/petsc-3.5.2/src/mat/interface/matrix.c PETSC ERROR: #4 add_matrix() line 765 in /home/dknez/software/libmesh-src/src/numerics/petsc_matrix.C -- This occurs when I try to set some entries of the matrix. Do you have any suggestions on how I can resolve this? Thanks! David On Sun, Feb 22, 2015 at 10:22 PM, Dmitry Karpeyev dkarp...@gmail.com wrote: On Sun Feb 22 2015 at 9:15:22 PM Barry Smith bsm...@mcs.anl.gov wrote: On Feb 22, 2015, at 9:09 PM, David Knezevic david.kneze...@akselos.com wrote: Hi Dmitry, Thanks for the suggestion. I tried MatSetType(mat,MATMPIAIJ) followed by MatXAIJSetPreallocation(...), but unfortunately this still gives me the same error as before: nnz cannot be greater than row length: local row 168 value 24 rowlength 0. I gather that the idea here is that MatSetType builds a new matrix object, and then I should be able to pre-allocate for that new matrix however I like, right? Was I supposed to clear the matrix object somehow before calling MatSetType? (I didn't do any sort of clear operation.) If the type doesn't change then MatSetType() won't do anything. You can try setting the type to BAIJ and then setting the type back to AIJ. This may/should clear out the matrix. Ah, yes. If the type is the same as before it does quit early, but changing the type and then back will clear out and rebuild the matrix. We need something like MatReset() to do the equivalent thing. As I said earlier, I'll make a dbg PETSc build, so hopefully that will help shed some light on what's going wrong for me. I think it's always a good idea to have a dbg build of PETSc when you doing things like these. Dmitry. Don't bother, what I suggested won't work. Barry Thanks, David On Sun, Feb 22, 2015 at 6:02 PM, Dmitry Karpeyev dkarp...@gmail.com wrote: David, It might be easier to just rebuild the whole matrix from scratch: you would in effect be doing all that with disassembling and resetting the preallocation. MatSetType(mat,MATMPIAIJ) or PetscObjectGetType((PetscObject)mat,type); MatSetType(mat,type); followed by MatXAIJSetPreallocation(...); should do. Dmitry. On Sun Feb 22 2015 at 4:45:46 PM Barry Smith bsm...@mcs.anl.gov wrote: Do not call for SeqAIJ matrix. Do not call before the first time you have preallocated and put entries in the matrix and done the MatAssemblyBegin/End() If it still crashes you'll need to try the debugger Barry On
Re: [petsc-users] MatMPIAIJSetPreallocation: nnz cannot be greater than row length
Hi David,
Thanks -- I do see a failure with two mpi ranks.
libMesh needs to change the way configure extracts PETSc information --
configuration data were moved:
conf -- lib/petsc-conf
${PETSC_ARCH}/conf -- ${PETSC_ARCH}/lib/petsc-conf
At one point I started looking at m4/petsc.m4, but that got put on the back
burner. For now making the relevant symlinks by hand lets you configure
and build libMesh with petsc/master.
Dmitry.
On Mon Feb 23 2015 at 9:15:44 AM David Knezevic david.kneze...@akselos.com
wrote:
Hi Dmitry,
Thanks very much for testing out the example.
examples/systems_of_equations/ex8 works fine for me in serial, but it
fails for me if I run with more than 1 MPI process. Can you try it with,
say, 2 or 4 MPI processes?
If we need access to MatReset in libMesh to get this to work, I'll be
happy to work on a libMesh pull request for that.
David
--
David J. Knezevic | CTO
Akselos | 17 Bay State Road | Boston, MA | 02215
Phone (office): +1-857-265-2238
Phone (mobile): +1-617-599-4755
Web: http://www.akselos.com
On Mon, Feb 23, 2015 at 10:08 AM, Dmitry Karpeyev karp...@mcs.anl.gov
wrote:
David,
What code are you running when you encounter this error? I'm trying to
reproduce it and
I tried examples/systems_of_equations/ex8, but it ran for me, ostensibly
to completion.
I have a small PETSc pull request that implements MatReset(), which
passes a small PETSc test,
but libMesh needs some work to be able to build against petsc/master
because of some recent
changes to PETSc.
Dmitry.
On Mon Feb 23 2015 at 7:17:06 AM David Knezevic
david.kneze...@akselos.com wrote:
Hi Barry, hi Dmitry,
I set the matrix to BAIJ and back to AIJ, and the code got a bit
further. But I now run into the error pasted below (Note that I'm now using
--with-debugging=1):
PETSC ERROR: - Error Message
--
PETSC ERROR: Petsc has generated inconsistent data
PETSC ERROR: MPIAIJ Matrix was assembled but is missing garray
PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html
for trouble shooting.
PETSC ERROR: Petsc Release Version 3.5.2, Sep, 08, 2014
PETSC ERROR: ./example-dbg on a arch-linux2-c-debug named david-Lenovo
by dknez Mon Feb 23 08:05:44 2015
PETSC ERROR: Configure options --with-shared-libraries=1
--with-debugging=1 --download-suitesparse=1 --download-parmetis=1
--download-blacs=1 --download-scalapack=1 --download-mumps=1
--download-metis --download-superlu_dist
--prefix=/home/dknez/software/libmesh_install/dbg_real/petsc
--download-hypre
PETSC ERROR: #1 MatCreateColmap_MPIAIJ_Private() line 361 in
/home/dknez/software/petsc-3.5.2/src/mat/impls/aij/mpi/mpiaij.c
PETSC ERROR: #2 MatSetValues_MPIAIJ() line 538 in
/home/dknez/software/petsc-3.5.2/src/mat/impls/aij/mpi/mpiaij.c
PETSC ERROR: #3 MatSetValues() line 1136 in /home/dknez/software/petsc-3.
5.2/src/mat/interface/matrix.c
PETSC ERROR: #4 add_matrix() line 765 in /home/dknez/software/libmesh-
src/src/numerics/petsc_matrix.C
--
This occurs when I try to set some entries of the matrix. Do you have
any suggestions on how I can resolve this?
Thanks!
David
On Sun, Feb 22, 2015 at 10:22 PM, Dmitry Karpeyev dkarp...@gmail.com
wrote:
On Sun Feb 22 2015 at 9:15:22 PM Barry Smith bsm...@mcs.anl.gov
wrote:
On Feb 22, 2015, at 9:09 PM, David Knezevic
david.kneze...@akselos.com wrote:
Hi Dmitry,
Thanks for the suggestion. I tried MatSetType(mat,MATMPIAIJ)
followed by MatXAIJSetPreallocation(...), but unfortunately this still
gives me the same error as before: nnz cannot be greater than row length:
local row 168 value 24 rowlength 0.
I gather that the idea here is that MatSetType builds a new matrix
object, and then I should be able to pre-allocate for that new matrix
however I like, right? Was I supposed to clear the matrix object somehow
before calling MatSetType? (I didn't do any sort of clear operation.)
If the type doesn't change then MatSetType() won't do anything. You
can try setting the type to BAIJ and then setting the type back to AIJ.
This may/should clear out the matrix.
Ah, yes. If the type is the same as before it does quit early, but
changing the type and then back will clear out and rebuild the matrix. We
need
something like MatReset() to do the equivalent thing.
As I said earlier, I'll make a dbg PETSc build, so hopefully that
will help shed some light on what's going wrong for me.
I think it's always a good idea to have a dbg build of PETSc when you
doing things like these.
Dmitry.
Don't bother, what I suggested won't work.
Barry
Thanks,
David
On Sun, Feb 22, 2015 at 6:02 PM, Dmitry Karpeyev dkarp...@gmail.com
wrote:
David,
It might be easier to just rebuild the whole matrix from scratch:
you would in effect be doing all that
Re: [petsc-users] DMNetworkGetEdgeRange() in parallel
I'm iterating through local edges given in DMNetworkGetEdgeRange(). For each edge, I extract or modify its corresponding value in a global petsc vector. Therefore that vector must have as many components as edges there are in the network. To extract the value in the vector, I use VecGetArray() and a variable counter that is incremented in each iteration. The array that I obtain in VecGetArray() has to be the same size than the edge range. That variable counter starts as 0, so if the array that I obtained in VecGetArray() is x_array, x_array[0] must be the component in the global vector that corresponds with the start edge given in DMNetworkGetEdgeRange() I need that global petsc vector because I will use it in other operations, it's not just data. Sorry for the confusion. Thanks in advance. Miguel On Mon, Feb 23, 2015 at 9:09 AM, Matthew Knepley knep...@gmail.com wrote: On Mon, Feb 23, 2015 at 8:42 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Thanks, that will help me. Now what I would like to have is the following: if I have two processors and ten edges, the partitioning results in the first processor having the edges 0-4 and the second processor, the edges 5-9. I also have a global vector with as many components as edges, 10. How can I partition it so the first processor also has the 0-4 components and the second, the 5-9 components of the vector? I think it would help to know what you want to accomplish. This is how you are proposing to do it.' If you just want to put data on edges, DMNetwork has a facility for that already. Thanks, Matt Miguel On Feb 23, 2015 8:08 AM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: Miguel, One possible way is to store the global numbering of any edge/vertex in the component attached to it. Once the mesh gets partitioned, the components are also distributed so you can easily retrieve the global number of any edge/vertex by accessing its component. This is what is done in the DMNetwork example pf.c although the global numbering is not used for anything. Shri From: Matthew Knepley knep...@gmail.com Date: Mon, 23 Feb 2015 07:54:34 -0600 To: Miguel Angel Salazar de Troya salazardetr...@gmail.com Cc: petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov Subject: Re: [petsc-users] DMNetworkGetEdgeRange() in parallel On Sun, Feb 22, 2015 at 3:59 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Thanks. Once I obtain that Index Set with the routine DMPlexCreateCellNumbering() (I assume that the edges in DMNetwork correspond to cells in DMPlex) can I use it to partition a vector with as many components as edges I have in my network? I do not completely understand the question. If you want a partition of the edges, you can use DMPlexCreatePartition() and its friend DMPlexDistribute(). What are you trying to do? Matt Thanks Miguel On Sun, Feb 22, 2015 at 12:15 PM, Matthew Knepley knep...@gmail.com wrote: On Sun, Feb 22, 2015 at 11:01 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Hi I noticed that the routine DMNetworkGetEdgeRange() returns the local indices for the edge range. Is there any way to obtain the global indices? So if my network has 10 edges, the processor 1 has the 0-4 edges and the processor 2, the 5-9 edges, how can I obtain this information? One of the points of DMPlex is we do not require a global numbering. Everything is numbered locally, and the PetscSF maps local numbers to local numbers in order to determine ownership. If you want to create a global numbering for some reason, you can using DMPlexCreatePointNumbering(). There are also cell and vertex versions that we use for output, so you could do it just for edges as well. Thanks, Matt Thanks Miguel -- *Miguel Angel Salazar de Troya* Graduate Research Assistant Department of Mechanical Science and Engineering University of Illinois at Urbana-Champaign (217) 550-2360 salaz...@illinois.edu -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener -- *Miguel Angel Salazar de Troya* Graduate Research Assistant Department of Mechanical Science and Engineering University of Illinois at Urbana-Champaign (217) 550-2360 salaz...@illinois.edu -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener -- *Miguel Angel Salazar de Troya* Graduate Research Assistant Department of Mechanical Science and Engineering University of Illinois at Urbana-Champaign (217)
Re: [petsc-users] MatMPIAIJSetPreallocation: nnz cannot be greater than row length
David, What code are you running when you encounter this error? I'm trying to reproduce it and I tried examples/systems_of_equations/ex8, but it ran for me, ostensibly to completion. I have a small PETSc pull request that implements MatReset(), which passes a small PETSc test, but libMesh needs some work to be able to build against petsc/master because of some recent changes to PETSc. Dmitry. On Mon Feb 23 2015 at 7:17:06 AM David Knezevic david.kneze...@akselos.com wrote: Hi Barry, hi Dmitry, I set the matrix to BAIJ and back to AIJ, and the code got a bit further. But I now run into the error pasted below (Note that I'm now using --with-debugging=1): PETSC ERROR: - Error Message -- PETSC ERROR: Petsc has generated inconsistent data PETSC ERROR: MPIAIJ Matrix was assembled but is missing garray PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting. PETSC ERROR: Petsc Release Version 3.5.2, Sep, 08, 2014 PETSC ERROR: ./example-dbg on a arch-linux2-c-debug named david-Lenovo by dknez Mon Feb 23 08:05:44 2015 PETSC ERROR: Configure options --with-shared-libraries=1 --with-debugging=1 --download-suitesparse=1 --download-parmetis=1 --download-blacs=1 --download-scalapack=1 --download-mumps=1 --download-metis --download-superlu_dist --prefix=/home/dknez/software/libmesh_install/dbg_real/petsc --download-hypre PETSC ERROR: #1 MatCreateColmap_MPIAIJ_Private() line 361 in /home/dknez/software/petsc-3.5.2/src/mat/impls/aij/mpi/mpiaij.c PETSC ERROR: #2 MatSetValues_MPIAIJ() line 538 in /home/dknez/software/petsc-3.5.2/src/mat/impls/aij/mpi/mpiaij.c PETSC ERROR: #3 MatSetValues() line 1136 in /home/dknez/software/petsc-3.5.2/src/mat/interface/matrix.c PETSC ERROR: #4 add_matrix() line 765 in /home/dknez/software/libmesh-src/src/numerics/petsc_matrix.C -- This occurs when I try to set some entries of the matrix. Do you have any suggestions on how I can resolve this? Thanks! David On Sun, Feb 22, 2015 at 10:22 PM, Dmitry Karpeyev dkarp...@gmail.com wrote: On Sun Feb 22 2015 at 9:15:22 PM Barry Smith bsm...@mcs.anl.gov wrote: On Feb 22, 2015, at 9:09 PM, David Knezevic david.kneze...@akselos.com wrote: Hi Dmitry, Thanks for the suggestion. I tried MatSetType(mat,MATMPIAIJ) followed by MatXAIJSetPreallocation(...), but unfortunately this still gives me the same error as before: nnz cannot be greater than row length: local row 168 value 24 rowlength 0. I gather that the idea here is that MatSetType builds a new matrix object, and then I should be able to pre-allocate for that new matrix however I like, right? Was I supposed to clear the matrix object somehow before calling MatSetType? (I didn't do any sort of clear operation.) If the type doesn't change then MatSetType() won't do anything. You can try setting the type to BAIJ and then setting the type back to AIJ. This may/should clear out the matrix. Ah, yes. If the type is the same as before it does quit early, but changing the type and then back will clear out and rebuild the matrix. We need something like MatReset() to do the equivalent thing. As I said earlier, I'll make a dbg PETSc build, so hopefully that will help shed some light on what's going wrong for me. I think it's always a good idea to have a dbg build of PETSc when you doing things like these. Dmitry. Don't bother, what I suggested won't work. Barry Thanks, David On Sun, Feb 22, 2015 at 6:02 PM, Dmitry Karpeyev dkarp...@gmail.com wrote: David, It might be easier to just rebuild the whole matrix from scratch: you would in effect be doing all that with disassembling and resetting the preallocation. MatSetType(mat,MATMPIAIJ) or PetscObjectGetType((PetscObject)mat,type); MatSetType(mat,type); followed by MatXAIJSetPreallocation(...); should do. Dmitry. On Sun Feb 22 2015 at 4:45:46 PM Barry Smith bsm...@mcs.anl.gov wrote: Do not call for SeqAIJ matrix. Do not call before the first time you have preallocated and put entries in the matrix and done the MatAssemblyBegin/End() If it still crashes you'll need to try the debugger Barry On Feb 22, 2015, at 4:09 PM, David Knezevic david.kneze...@akselos.com wrote: Hi Barry, Thanks for your help, much appreciated. I added a prototype for MatDisAssemble_MPIAIJ: PETSC_INTERN PetscErrorCode MatDisAssemble_MPIAIJ(Mat); and I added a call to MatDisAssemble_MPIAIJ before MatMPIAIJSetPreallocation. However, I get a segfault on the call to MatDisAssemble_MPIAIJ. The segfault occurs in both serial and parallel. FYI, I'm using Petsc 3.5.2, and I'm not using a non-debug build (though I could rebuild PETSc in debug mode if you think that would help figure out what's happening
Re: [petsc-users] MatMPIAIJSetPreallocation: nnz cannot be greater than row length
OK, sounds good. Let me know if I can help with the digging.
David
On Mon, Feb 23, 2015 at 11:10 AM, Dmitry Karpeyev karp...@mcs.anl.gov
wrote:
I just tried building against petsc/master, but there needs to be more
work on libMesh before it can work with petsc/master:
the new VecLockPush()/Pop() stuff isn't respected by vector manipulation
in libMesh.
I put a hack equivalent to MatReset() into your branch (patch attached, in
case you want to take a look at it),
but it generates the same error in MatCreateColmap that you reported
earlier. It's odd that it occurs
on the second nonlinear iteration. I'll have to dig a bit deeper to
see what's going on.
Dmitry.
On Mon Feb 23 2015 at 10:03:33 AM David Knezevic
david.kneze...@akselos.com wrote:
Hi Dmitry,
OK, good to hear we're seeing the same behavior for the example.
Regarding this comment:
libMesh needs to change the way configure extracts PETSc information --
configuration data were moved:
conf -- lib/petsc-conf
${PETSC_ARCH}/conf -- ${PETSC_ARCH}/lib/petsc-conf
At one point I started looking at m4/petsc.m4, but that got put on the
back burner. For now making the relevant symlinks by hand lets you
configure and build libMesh with petsc/master.
So do you suggest that the next step here is to build libmesh against
petsc/master so that we can try your PETSc pull request that implements
MatReset() to see if that gets this example working?
David
On Mon Feb 23 2015 at 9:15:44 AM David Knezevic
david.kneze...@akselos.com wrote:
Hi Dmitry,
Thanks very much for testing out the example.
examples/systems_of_equations/ex8 works fine for me in serial, but it
fails for me if I run with more than 1 MPI process. Can you try it with,
say, 2 or 4 MPI processes?
If we need access to MatReset in libMesh to get this to work, I'll be
happy to work on a libMesh pull request for that.
David
--
David J. Knezevic | CTO
Akselos | 17 Bay State Road | Boston, MA | 02215
Phone (office): +1-857-265-2238
Phone (mobile): +1-617-599-4755
Web: http://www.akselos.com
On Mon, Feb 23, 2015 at 10:08 AM, Dmitry Karpeyev karp...@mcs.anl.gov
wrote:
David,
What code are you running when you encounter this error? I'm trying
to reproduce it and
I tried examples/systems_of_equations/ex8, but it ran for me,
ostensibly to completion.
I have a small PETSc pull request that implements MatReset(), which
passes a small PETSc test,
but libMesh needs some work to be able to build against petsc/master
because of some recent
changes to PETSc.
Dmitry.
On Mon Feb 23 2015 at 7:17:06 AM David Knezevic
david.kneze...@akselos.com wrote:
Hi Barry, hi Dmitry,
I set the matrix to BAIJ and back to AIJ, and the code got a bit
further. But I now run into the error pasted below (Note that I'm now
using
--with-debugging=1):
PETSC ERROR: - Error Message
--
PETSC ERROR: Petsc has generated inconsistent data
PETSC ERROR: MPIAIJ Matrix was assembled but is missing garray
PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html
for trouble shooting.
PETSC ERROR: Petsc Release Version 3.5.2, Sep, 08, 2014
PETSC ERROR: ./example-dbg on a arch-linux2-c-debug named
david-Lenovo by dknez Mon Feb 23 08:05:44 2015
PETSC ERROR: Configure options --with-shared-libraries=1
--with-debugging=1 --download-suitesparse=1 --download-parmetis=1
--download-blacs=1 --download-scalapack=1 --download-mumps=1
--download-metis --download-superlu_dist
--prefix=/home/dknez/software/libmesh_install/dbg_real/petsc
--download-hypre
PETSC ERROR: #1 MatCreateColmap_MPIAIJ_Private() line 361 in
/home/dknez/software/petsc-3.5.2/src/mat/impls/aij/mpi/mpiaij.c
PETSC ERROR: #2 MatSetValues_MPIAIJ() line 538 in
/home/dknez/software/petsc-3.5.2/src/mat/impls/aij/mpi/mpiaij.c
PETSC ERROR: #3 MatSetValues() line 1136 in
/home/dknez/software/petsc-3.5.2/src/mat/interface/matrix.c
PETSC ERROR: #4 add_matrix() line 765 in /home/dknez/software/libmesh-
src/src/numerics/petsc_matrix.C
--
This occurs when I try to set some entries of the matrix. Do you have
any suggestions on how I can resolve this?
Thanks!
David
On Sun, Feb 22, 2015 at 10:22 PM, Dmitry Karpeyev dkarp...@gmail.com
wrote:
On Sun Feb 22 2015 at 9:15:22 PM Barry Smith bsm...@mcs.anl.gov
wrote:
On Feb 22, 2015, at 9:09 PM, David Knezevic
david.kneze...@akselos.com wrote:
Hi Dmitry,
Thanks for the suggestion. I tried MatSetType(mat,MATMPIAIJ)
followed by MatXAIJSetPreallocation(...), but unfortunately this still
gives me the same error as before: nnz cannot be greater than row
length:
local row 168 value 24 rowlength 0.
I gather that the idea here is that MatSetType builds a new
matrix object, and then I should be able to pre-allocate for that new
matrix however I like,
Re: [petsc-users] DMNetworkGetEdgeRange() in parallel
On Mon, Feb 23, 2015 at 9:27 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: I'm iterating through local edges given in DMNetworkGetEdgeRange(). For each edge, I extract or modify its corresponding value in a global petsc vector. Therefore that vector must have as many components as edges there are in the network. To extract the value in the vector, I use VecGetArray() and a variable counter that is incremented in each iteration. The array that I obtain in VecGetArray() has to be the same size than the edge range. That variable counter starts as 0, so if the array that I obtained in VecGetArray() is x_array, x_array[0] must be the component in the global vector that corresponds with the start edge given in DMNetworkGetEdgeRange() I need that global petsc vector because I will use it in other operations, it's not just data. Sorry for the confusion. Thanks in advance. This sounds like an assembly operation. The usual paradigm is to compute in the local space, and then communicate to get to the global space. So you would make a PetscSection that had 1 (or some) unknowns on each cell (edge) and then you can use DMCreateGlobal/LocalVector() and DMLocalToGlobal() to do this. Does that make sense? Thanks, Matt Miguel On Mon, Feb 23, 2015 at 9:09 AM, Matthew Knepley knep...@gmail.com wrote: On Mon, Feb 23, 2015 at 8:42 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Thanks, that will help me. Now what I would like to have is the following: if I have two processors and ten edges, the partitioning results in the first processor having the edges 0-4 and the second processor, the edges 5-9. I also have a global vector with as many components as edges, 10. How can I partition it so the first processor also has the 0-4 components and the second, the 5-9 components of the vector? I think it would help to know what you want to accomplish. This is how you are proposing to do it.' If you just want to put data on edges, DMNetwork has a facility for that already. Thanks, Matt Miguel On Feb 23, 2015 8:08 AM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov wrote: Miguel, One possible way is to store the global numbering of any edge/vertex in the component attached to it. Once the mesh gets partitioned, the components are also distributed so you can easily retrieve the global number of any edge/vertex by accessing its component. This is what is done in the DMNetwork example pf.c although the global numbering is not used for anything. Shri From: Matthew Knepley knep...@gmail.com Date: Mon, 23 Feb 2015 07:54:34 -0600 To: Miguel Angel Salazar de Troya salazardetr...@gmail.com Cc: petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov Subject: Re: [petsc-users] DMNetworkGetEdgeRange() in parallel On Sun, Feb 22, 2015 at 3:59 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Thanks. Once I obtain that Index Set with the routine DMPlexCreateCellNumbering() (I assume that the edges in DMNetwork correspond to cells in DMPlex) can I use it to partition a vector with as many components as edges I have in my network? I do not completely understand the question. If you want a partition of the edges, you can use DMPlexCreatePartition() and its friend DMPlexDistribute(). What are you trying to do? Matt Thanks Miguel On Sun, Feb 22, 2015 at 12:15 PM, Matthew Knepley knep...@gmail.com wrote: On Sun, Feb 22, 2015 at 11:01 AM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Hi I noticed that the routine DMNetworkGetEdgeRange() returns the local indices for the edge range. Is there any way to obtain the global indices? So if my network has 10 edges, the processor 1 has the 0-4 edges and the processor 2, the 5-9 edges, how can I obtain this information? One of the points of DMPlex is we do not require a global numbering. Everything is numbered locally, and the PetscSF maps local numbers to local numbers in order to determine ownership. If you want to create a global numbering for some reason, you can using DMPlexCreatePointNumbering(). There are also cell and vertex versions that we use for output, so you could do it just for edges as well. Thanks, Matt Thanks Miguel -- *Miguel Angel Salazar de Troya* Graduate Research Assistant Department of Mechanical Science and Engineering University of Illinois at Urbana-Champaign (217) 550-2360 salaz...@illinois.edu -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener -- *Miguel Angel Salazar de Troya* Graduate Research Assistant Department of Mechanical Science and Engineering University of Illinois at Urbana-Champaign (217) 550-2360 salaz...@illinois.edu -- What most experimenters take for granted before they begin their
Re: [petsc-users] MatMPIAIJSetPreallocation: nnz cannot be greater than row length
Hi Dmitry,
OK, good to hear we're seeing the same behavior for the example.
Regarding this comment:
libMesh needs to change the way configure extracts PETSc information --
configuration data were moved:
conf -- lib/petsc-conf
${PETSC_ARCH}/conf -- ${PETSC_ARCH}/lib/petsc-conf
At one point I started looking at m4/petsc.m4, but that got put on the
back burner. For now making the relevant symlinks by hand lets you
configure and build libMesh with petsc/master.
So do you suggest that the next step here is to build libmesh against
petsc/master so that we can try your PETSc pull request that implements
MatReset() to see if that gets this example working?
David
On Mon Feb 23 2015 at 9:15:44 AM David Knezevic
david.kneze...@akselos.com wrote:
Hi Dmitry,
Thanks very much for testing out the example.
examples/systems_of_equations/ex8 works fine for me in serial, but it
fails for me if I run with more than 1 MPI process. Can you try it with,
say, 2 or 4 MPI processes?
If we need access to MatReset in libMesh to get this to work, I'll be
happy to work on a libMesh pull request for that.
David
--
David J. Knezevic | CTO
Akselos | 17 Bay State Road | Boston, MA | 02215
Phone (office): +1-857-265-2238
Phone (mobile): +1-617-599-4755
Web: http://www.akselos.com
On Mon, Feb 23, 2015 at 10:08 AM, Dmitry Karpeyev karp...@mcs.anl.gov
wrote:
David,
What code are you running when you encounter this error? I'm trying to
reproduce it and
I tried examples/systems_of_equations/ex8, but it ran for me,
ostensibly to completion.
I have a small PETSc pull request that implements MatReset(), which
passes a small PETSc test,
but libMesh needs some work to be able to build against petsc/master
because of some recent
changes to PETSc.
Dmitry.
On Mon Feb 23 2015 at 7:17:06 AM David Knezevic
david.kneze...@akselos.com wrote:
Hi Barry, hi Dmitry,
I set the matrix to BAIJ and back to AIJ, and the code got a bit
further. But I now run into the error pasted below (Note that I'm now using
--with-debugging=1):
PETSC ERROR: - Error Message
--
PETSC ERROR: Petsc has generated inconsistent data
PETSC ERROR: MPIAIJ Matrix was assembled but is missing garray
PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html
for trouble shooting.
PETSC ERROR: Petsc Release Version 3.5.2, Sep, 08, 2014
PETSC ERROR: ./example-dbg on a arch-linux2-c-debug named david-Lenovo
by dknez Mon Feb 23 08:05:44 2015
PETSC ERROR: Configure options --with-shared-libraries=1
--with-debugging=1 --download-suitesparse=1 --download-parmetis=1
--download-blacs=1 --download-scalapack=1 --download-mumps=1
--download-metis --download-superlu_dist
--prefix=/home/dknez/software/libmesh_install/dbg_real/petsc
--download-hypre
PETSC ERROR: #1 MatCreateColmap_MPIAIJ_Private() line 361 in
/home/dknez/software/petsc-3.5.2/src/mat/impls/aij/mpi/mpiaij.c
PETSC ERROR: #2 MatSetValues_MPIAIJ() line 538 in
/home/dknez/software/petsc-3.5.2/src/mat/impls/aij/mpi/mpiaij.c
PETSC ERROR: #3 MatSetValues() line 1136 in
/home/dknez/software/petsc-3.5.2/src/mat/interface/matrix.c
PETSC ERROR: #4 add_matrix() line 765 in /home/dknez/software/libmesh-
src/src/numerics/petsc_matrix.C
--
This occurs when I try to set some entries of the matrix. Do you have
any suggestions on how I can resolve this?
Thanks!
David
On Sun, Feb 22, 2015 at 10:22 PM, Dmitry Karpeyev dkarp...@gmail.com
wrote:
On Sun Feb 22 2015 at 9:15:22 PM Barry Smith bsm...@mcs.anl.gov
wrote:
On Feb 22, 2015, at 9:09 PM, David Knezevic
david.kneze...@akselos.com wrote:
Hi Dmitry,
Thanks for the suggestion. I tried MatSetType(mat,MATMPIAIJ)
followed by MatXAIJSetPreallocation(...), but unfortunately this still
gives me the same error as before: nnz cannot be greater than row
length:
local row 168 value 24 rowlength 0.
I gather that the idea here is that MatSetType builds a new matrix
object, and then I should be able to pre-allocate for that new matrix
however I like, right? Was I supposed to clear the matrix object somehow
before calling MatSetType? (I didn't do any sort of clear operation.)
If the type doesn't change then MatSetType() won't do anything. You
can try setting the type to BAIJ and then setting the type back to AIJ.
This may/should clear out the matrix.
Ah, yes. If the type is the same as before it does quit early, but
changing the type and then back will clear out and rebuild the matrix.
We
need
something like MatReset() to do the equivalent thing.
As I said earlier, I'll make a dbg PETSc build, so hopefully that
will help shed some light on what's going wrong for me.
I think it's always a good idea to have a dbg build of PETSc when you
doing things like these.
Dmitry.
Don't bother, what I
Re: [petsc-users] MatMPIAIJSetPreallocation: nnz cannot be greater than row length
I just tried building against petsc/master, but there needs to be more work
on libMesh before it can work with petsc/master:
the new VecLockPush()/Pop() stuff isn't respected by vector manipulation in
libMesh.
I put a hack equivalent to MatReset() into your branch (patch attached, in
case you want to take a look at it),
but it generates the same error in MatCreateColmap that you reported
earlier. It's odd that it occurs
on the second nonlinear iteration. I'll have to dig a bit deeper to see
what's going on.
Dmitry.
On Mon Feb 23 2015 at 10:03:33 AM David Knezevic david.kneze...@akselos.com
wrote:
Hi Dmitry,
OK, good to hear we're seeing the same behavior for the example.
Regarding this comment:
libMesh needs to change the way configure extracts PETSc information --
configuration data were moved:
conf -- lib/petsc-conf
${PETSC_ARCH}/conf -- ${PETSC_ARCH}/lib/petsc-conf
At one point I started looking at m4/petsc.m4, but that got put on the
back burner. For now making the relevant symlinks by hand lets you
configure and build libMesh with petsc/master.
So do you suggest that the next step here is to build libmesh against
petsc/master so that we can try your PETSc pull request that implements
MatReset() to see if that gets this example working?
David
On Mon Feb 23 2015 at 9:15:44 AM David Knezevic
david.kneze...@akselos.com wrote:
Hi Dmitry,
Thanks very much for testing out the example.
examples/systems_of_equations/ex8 works fine for me in serial, but it
fails for me if I run with more than 1 MPI process. Can you try it with,
say, 2 or 4 MPI processes?
If we need access to MatReset in libMesh to get this to work, I'll be
happy to work on a libMesh pull request for that.
David
--
David J. Knezevic | CTO
Akselos | 17 Bay State Road | Boston, MA | 02215
Phone (office): +1-857-265-2238
Phone (mobile): +1-617-599-4755
Web: http://www.akselos.com
On Mon, Feb 23, 2015 at 10:08 AM, Dmitry Karpeyev karp...@mcs.anl.gov
wrote:
David,
What code are you running when you encounter this error? I'm trying to
reproduce it and
I tried examples/systems_of_equations/ex8, but it ran for me,
ostensibly to completion.
I have a small PETSc pull request that implements MatReset(), which
passes a small PETSc test,
but libMesh needs some work to be able to build against petsc/master
because of some recent
changes to PETSc.
Dmitry.
On Mon Feb 23 2015 at 7:17:06 AM David Knezevic
david.kneze...@akselos.com wrote:
Hi Barry, hi Dmitry,
I set the matrix to BAIJ and back to AIJ, and the code got a bit
further. But I now run into the error pasted below (Note that I'm now
using
--with-debugging=1):
PETSC ERROR: - Error Message
--
PETSC ERROR: Petsc has generated inconsistent data
PETSC ERROR: MPIAIJ Matrix was assembled but is missing garray
PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html
for trouble shooting.
PETSC ERROR: Petsc Release Version 3.5.2, Sep, 08, 2014
PETSC ERROR: ./example-dbg on a arch-linux2-c-debug named david-Lenovo
by dknez Mon Feb 23 08:05:44 2015
PETSC ERROR: Configure options --with-shared-libraries=1
--with-debugging=1 --download-suitesparse=1 --download-parmetis=1
--download-blacs=1 --download-scalapack=1 --download-mumps=1
--download-metis --download-superlu_dist
--prefix=/home/dknez/software/libmesh_install/dbg_real/petsc
--download-hypre
PETSC ERROR: #1 MatCreateColmap_MPIAIJ_Private() line 361 in
/home/dknez/software/petsc-3.5.2/src/mat/impls/aij/mpi/mpiaij.c
PETSC ERROR: #2 MatSetValues_MPIAIJ() line 538 in
/home/dknez/software/petsc-3.5.2/src/mat/impls/aij/mpi/mpiaij.c
PETSC ERROR: #3 MatSetValues() line 1136 in
/home/dknez/software/petsc-3.5.2/src/mat/interface/matrix.c
PETSC ERROR: #4 add_matrix() line 765 in /home/dknez/software/libmesh-
src/src/numerics/petsc_matrix.C
--
This occurs when I try to set some entries of the matrix. Do you have
any suggestions on how I can resolve this?
Thanks!
David
On Sun, Feb 22, 2015 at 10:22 PM, Dmitry Karpeyev dkarp...@gmail.com
wrote:
On Sun Feb 22 2015 at 9:15:22 PM Barry Smith bsm...@mcs.anl.gov
wrote:
On Feb 22, 2015, at 9:09 PM, David Knezevic
david.kneze...@akselos.com wrote:
Hi Dmitry,
Thanks for the suggestion. I tried MatSetType(mat,MATMPIAIJ)
followed by MatXAIJSetPreallocation(...), but unfortunately this still
gives me the same error as before: nnz cannot be greater than row
length:
local row 168 value 24 rowlength 0.
I gather that the idea here is that MatSetType builds a new matrix
object, and then I should be able to pre-allocate for that new matrix
however I like, right? Was I supposed to clear the matrix object somehow
before calling MatSetType? (I didn't do any sort of clear operation.)
If the type doesn't change then
Re: [petsc-users] DMNetworkGetEdgeRange() in parallel
I think you should call DMClone after partitioning (DMDistribute).
Shri
From: Miguel Angel Salazar de Troya
salazardetr...@gmail.commailto:salazardetr...@gmail.com
Date: Mon, 23 Feb 2015 14:15:00 -0600
To: Matthew Knepley knep...@gmail.commailto:knep...@gmail.com
Cc: Shri abhy...@mcs.anl.govmailto:abhy...@mcs.anl.gov,
petsc-users@mcs.anl.govmailto:petsc-users@mcs.anl.gov
petsc-users@mcs.anl.govmailto:petsc-users@mcs.anl.gov
Subject: Re: [petsc-users] DMNetworkGetEdgeRange() in parallel
Thanks a lot, the partition should be done before setting up the section, right?
Miguel
On Mon, Feb 23, 2015 at 2:05 PM, Matthew Knepley
knep...@gmail.commailto:knep...@gmail.com wrote:
On Mon, Feb 23, 2015 at 1:40 PM, Miguel Angel Salazar de Troya
salazardetr...@gmail.commailto:salazardetr...@gmail.com wrote:
Wouldn't including the edge variables in the global vector make the code
slower? I'm using the global vector in a TS, using one of the explicit RK
schemes. The edge variables would not be updated in the RHSFunction evaluation.
I only change the edge variables in the TSUpdate. If the global vector had the
edge variables, it would be a much larger vector, and all the vector operations
performed by the TS would be slower. Although the vector F returned by the
RHSFunction would be zero in the edge variable components. I guess that being
the vector sparse that would not be a problem.
I think I'm more interested in the PetscSection approach because it might
require less modifications in my code. However, I don't know how I could do
this. Maybe something like this?
PetscSectionCreate(PETSC_COMM_WORLD, s);
PetscSectionSetNumFields(s, 1);
PetscSectionSetFieldComponents(s, 0, 1);
// Now to set the chart, I pick the edge range
DMNetworkGetEdgeRange(dm, eStart, eEnd
PetscSectionSetChart(s, eStart, eEnd);
for(PetscInt e = eStart; c eEnd; ++e) {
PetscSectionSetDof(s, e, 1);
PetscSectionSetFieldDof(s, e, 1, 1);
It should be PetscSectionSetFieldDof(s, e, 0, 1);
}
PetscSectionSetUp(s);
Now in the manual I see this:
First you want to do:
DMClone(dm, dmEdge);
and then use dmEdge below.
DMSetDefaultSection(dm, s);
DMGetLocalVector(dm, localVec);
DMGetGlobalVector(dm, globalVec);
Setting up the default section in the DM would interfere with the section
already set up with the variables in the vertices?
Yep, thats why you would use a clone.
Thanks,
Matt
Thanks a lot for your responses.
On Mon, Feb 23, 2015 at 11:37 AM, Matthew Knepley
knep...@gmail.commailto:knep...@gmail.com wrote:
On Mon, Feb 23, 2015 at 9:27 AM, Miguel Angel Salazar de Troya
salazardetr...@gmail.commailto:salazardetr...@gmail.com wrote:
I'm iterating through local edges given in DMNetworkGetEdgeRange(). For each
edge, I extract or modify its corresponding value in a global petsc vector.
Therefore that vector must have as many components as edges there are in the
network. To extract the value in the vector, I use VecGetArray() and a variable
counter that is incremented in each iteration. The array that I obtain in
VecGetArray() has to be the same size than the edge range. That variable
counter starts as 0, so if the array that I obtained in VecGetArray() is
x_array, x_array[0] must be the component in the global vector that corresponds
with the start edge given in DMNetworkGetEdgeRange()
I need that global petsc vector because I will use it in other operations, it's
not just data. Sorry for the confusion. Thanks in advance.
This sounds like an assembly operation. The usual paradigm is to compute in the
local space, and then communicate to get to the global space. So you would make
a PetscSection that had 1 (or some) unknowns on each cell (edge) and then you
can use DMCreateGlobal/LocalVector() and DMLocalToGlobal() to do this.
Does that make sense?
Thanks,
Matt
Miguel
On Mon, Feb 23, 2015 at 9:09 AM, Matthew Knepley
knep...@gmail.commailto:knep...@gmail.com wrote:
On Mon, Feb 23, 2015 at 8:42 AM, Miguel Angel Salazar de Troya
salazardetr...@gmail.commailto:salazardetr...@gmail.com wrote:
Thanks, that will help me. Now what I would like to have is the following: if I
have two processors and ten edges, the partitioning results in the first
processor having the edges 0-4 and the second processor, the edges 5-9. I also
have a global vector with as many components as edges, 10. How can I partition
it so the first processor also has the 0-4 components and the second, the 5-9
components of the vector?
I think it would help to know what you want to accomplish. This is how you are
proposing to do it.'
If you just want to put data on edges, DMNetwork has a facility for that
already.
Thanks,
Matt
Miguel
On Feb 23, 2015 8:08 AM, Abhyankar, Shrirang G.
abhy...@mcs.anl.govmailto:abhy...@mcs.anl.gov wrote:
Miguel,
One possible way is to store the global numbering of any edge/vertex in the
component attached to it. Once the mesh gets partitioned, the components are
Re: [petsc-users] DMNetworkGetEdgeRange() in parallel
On Mon, Feb 23, 2015 at 2:15 PM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
Thanks a lot, the partition should be done before setting up the section,
right?
The partition will be automatic. All you have to do is make the local
section. The DM is already partitioned,
and the Section will inherit that.
Matt
Miguel
On Mon, Feb 23, 2015 at 2:05 PM, Matthew Knepley knep...@gmail.com
wrote:
On Mon, Feb 23, 2015 at 1:40 PM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
Wouldn't including the edge variables in the global vector make the code
slower? I'm using the global vector in a TS, using one of the explicit RK
schemes. The edge variables would not be updated in the RHSFunction
evaluation. I only change the edge variables in the TSUpdate. If the global
vector had the edge variables, it would be a much larger vector, and all
the vector operations performed by the TS would be slower. Although the
vector F returned by the RHSFunction would be zero in the edge variable
components. I guess that being the vector sparse that would not be a
problem.
I think I'm more interested in the PetscSection approach because it
might require less modifications in my code. However, I don't know how I
could do this. Maybe something like this?
PetscSectionCreate(PETSC_COMM_WORLD, s);
PetscSectionSetNumFields(s, 1);
PetscSectionSetFieldComponents(s, 0, 1);
// Now to set the chart, I pick the edge range
DMNetworkGetEdgeRange(dm, eStart, eEnd
PetscSectionSetChart(s, eStart, eEnd);
for(PetscInt e = eStart; c eEnd; ++e) {
PetscSectionSetDof(s, e, 1);
PetscSectionSetFieldDof(s, e, 1, 1);
It should be PetscSectionSetFieldDof(s, e, 0, 1);
}
PetscSectionSetUp(s);
Now in the manual I see this:
First you want to do:
DMClone(dm, dmEdge);
and then use dmEdge below.
DMSetDefaultSection(dm, s);
DMGetLocalVector(dm, localVec);
DMGetGlobalVector(dm, globalVec);
Setting up the default section in the DM would interfere with the
section already set up with the variables in the vertices?
Yep, thats why you would use a clone.
Thanks,
Matt
Thanks a lot for your responses.
On Mon, Feb 23, 2015 at 11:37 AM, Matthew Knepley knep...@gmail.com
wrote:
On Mon, Feb 23, 2015 at 9:27 AM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
I'm iterating through local edges given in DMNetworkGetEdgeRange().
For each edge, I extract or modify its corresponding value in a global
petsc vector. Therefore that vector must have as many components as edges
there are in the network. To extract the value in the vector, I use
VecGetArray() and a variable counter that is incremented in each
iteration.
The array that I obtain in VecGetArray() has to be the same size than
the edge range. That variable counter starts as 0, so if the array that I
obtained in VecGetArray() is x_array, x_array[0] must be the
component in the global vector that corresponds with the start edge given
in DMNetworkGetEdgeRange()
I need that global petsc vector because I will use it in other
operations, it's not just data. Sorry for the confusion. Thanks in
advance.
This sounds like an assembly operation. The usual paradigm is to
compute in the local space, and then communicate to get to the global
space. So you would make a PetscSection that had 1 (or some) unknowns on
each cell (edge) and then you can use DMCreateGlobal/LocalVector() and
DMLocalToGlobal() to do this.
Does that make sense?
Thanks,
Matt
Miguel
On Mon, Feb 23, 2015 at 9:09 AM, Matthew Knepley knep...@gmail.com
wrote:
On Mon, Feb 23, 2015 at 8:42 AM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
Thanks, that will help me. Now what I would like to have is the
following: if I have two processors and ten edges, the partitioning
results
in the first processor having the edges 0-4 and the second processor,
the
edges 5-9. I also have a global vector with as many components as edges,
10. How can I partition it so the first processor also has the 0-4
components and the second, the 5-9 components of the vector?
I think it would help to know what you want to accomplish. This is
how you are proposing to do it.'
If you just want to put data on edges, DMNetwork has a facility for
that already.
Thanks,
Matt
Miguel
On Feb 23, 2015 8:08 AM, Abhyankar, Shrirang G.
abhy...@mcs.anl.gov wrote:
Miguel,
One possible way is to store the global numbering of any
edge/vertex in the component attached to it. Once the mesh gets
partitioned, the components are also distributed so you can easily
retrieve
the global number of any edge/vertex by accessing its component. This
is
what is done in the DMNetwork example pf.c although the global
numbering is
not used for anything.
Shri
From: Matthew Knepley knep...@gmail.com
Date: Mon, 23 Feb 2015 07:54:34 -0600
To: Miguel Angel Salazar de Troya
Re: [petsc-users] DMNetworkGetEdgeRange() in parallel
On Mon, Feb 23, 2015 at 1:40 PM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
Wouldn't including the edge variables in the global vector make the code
slower? I'm using the global vector in a TS, using one of the explicit RK
schemes. The edge variables would not be updated in the RHSFunction
evaluation. I only change the edge variables in the TSUpdate. If the global
vector had the edge variables, it would be a much larger vector, and all
the vector operations performed by the TS would be slower. Although the
vector F returned by the RHSFunction would be zero in the edge variable
components. I guess that being the vector sparse that would not be a
problem.
I think I'm more interested in the PetscSection approach because it might
require less modifications in my code. However, I don't know how I could do
this. Maybe something like this?
PetscSectionCreate(PETSC_COMM_WORLD, s);
PetscSectionSetNumFields(s, 1);
PetscSectionSetFieldComponents(s, 0, 1);
// Now to set the chart, I pick the edge range
DMNetworkGetEdgeRange(dm, eStart, eEnd
PetscSectionSetChart(s, eStart, eEnd);
for(PetscInt e = eStart; c eEnd; ++e) {
PetscSectionSetDof(s, e, 1);
PetscSectionSetFieldDof(s, e, 1, 1);
It should be PetscSectionSetFieldDof(s, e, 0, 1);
}
PetscSectionSetUp(s);
Now in the manual I see this:
First you want to do:
DMClone(dm, dmEdge);
and then use dmEdge below.
DMSetDefaultSection(dm, s);
DMGetLocalVector(dm, localVec);
DMGetGlobalVector(dm, globalVec);
Setting up the default section in the DM would interfere with the section
already set up with the variables in the vertices?
Yep, thats why you would use a clone.
Thanks,
Matt
Thanks a lot for your responses.
On Mon, Feb 23, 2015 at 11:37 AM, Matthew Knepley knep...@gmail.com
wrote:
On Mon, Feb 23, 2015 at 9:27 AM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
I'm iterating through local edges given in DMNetworkGetEdgeRange(). For
each edge, I extract or modify its corresponding value in a global petsc
vector. Therefore that vector must have as many components as edges there
are in the network. To extract the value in the vector, I use VecGetArray()
and a variable counter that is incremented in each iteration. The array
that I obtain in VecGetArray() has to be the same size than the edge
range. That variable counter starts as 0, so if the array that I obtained
in VecGetArray() is x_array, x_array[0] must be the component in the
global vector that corresponds with the start edge given in
DMNetworkGetEdgeRange()
I need that global petsc vector because I will use it in other
operations, it's not just data. Sorry for the confusion. Thanks in advance.
This sounds like an assembly operation. The usual paradigm is to compute
in the local space, and then communicate to get to the global space. So you
would make a PetscSection that had 1 (or some) unknowns on each cell (edge)
and then you can use DMCreateGlobal/LocalVector() and DMLocalToGlobal() to
do this.
Does that make sense?
Thanks,
Matt
Miguel
On Mon, Feb 23, 2015 at 9:09 AM, Matthew Knepley knep...@gmail.com
wrote:
On Mon, Feb 23, 2015 at 8:42 AM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
Thanks, that will help me. Now what I would like to have is the
following: if I have two processors and ten edges, the partitioning
results
in the first processor having the edges 0-4 and the second processor, the
edges 5-9. I also have a global vector with as many components as edges,
10. How can I partition it so the first processor also has the 0-4
components and the second, the 5-9 components of the vector?
I think it would help to know what you want to accomplish. This is how
you are proposing to do it.'
If you just want to put data on edges, DMNetwork has a facility for
that already.
Thanks,
Matt
Miguel
On Feb 23, 2015 8:08 AM, Abhyankar, Shrirang G. abhy...@mcs.anl.gov
wrote:
Miguel,
One possible way is to store the global numbering of any
edge/vertex in the component attached to it. Once the mesh gets
partitioned, the components are also distributed so you can easily
retrieve
the global number of any edge/vertex by accessing its component. This is
what is done in the DMNetwork example pf.c although the global numbering
is
not used for anything.
Shri
From: Matthew Knepley knep...@gmail.com
Date: Mon, 23 Feb 2015 07:54:34 -0600
To: Miguel Angel Salazar de Troya salazardetr...@gmail.com
Cc: petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov
Subject: Re: [petsc-users] DMNetworkGetEdgeRange() in parallel
On Sun, Feb 22, 2015 at 3:59 PM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
Thanks. Once I obtain that Index Set with the routine
DMPlexCreateCellNumbering()
(I assume that the edges in DMNetwork correspond to cells in DMPlex)
can I
use it to partition a vector with as
Re: [petsc-users] DMNetworkGetEdgeRange() in parallel
Thanks a lot, the partition should be done before setting up the section,
right?
Miguel
On Mon, Feb 23, 2015 at 2:05 PM, Matthew Knepley knep...@gmail.com wrote:
On Mon, Feb 23, 2015 at 1:40 PM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
Wouldn't including the edge variables in the global vector make the code
slower? I'm using the global vector in a TS, using one of the explicit RK
schemes. The edge variables would not be updated in the RHSFunction
evaluation. I only change the edge variables in the TSUpdate. If the global
vector had the edge variables, it would be a much larger vector, and all
the vector operations performed by the TS would be slower. Although the
vector F returned by the RHSFunction would be zero in the edge variable
components. I guess that being the vector sparse that would not be a
problem.
I think I'm more interested in the PetscSection approach because it might
require less modifications in my code. However, I don't know how I could do
this. Maybe something like this?
PetscSectionCreate(PETSC_COMM_WORLD, s);
PetscSectionSetNumFields(s, 1);
PetscSectionSetFieldComponents(s, 0, 1);
// Now to set the chart, I pick the edge range
DMNetworkGetEdgeRange(dm, eStart, eEnd
PetscSectionSetChart(s, eStart, eEnd);
for(PetscInt e = eStart; c eEnd; ++e) {
PetscSectionSetDof(s, e, 1);
PetscSectionSetFieldDof(s, e, 1, 1);
It should be PetscSectionSetFieldDof(s, e, 0, 1);
}
PetscSectionSetUp(s);
Now in the manual I see this:
First you want to do:
DMClone(dm, dmEdge);
and then use dmEdge below.
DMSetDefaultSection(dm, s);
DMGetLocalVector(dm, localVec);
DMGetGlobalVector(dm, globalVec);
Setting up the default section in the DM would interfere with the section
already set up with the variables in the vertices?
Yep, thats why you would use a clone.
Thanks,
Matt
Thanks a lot for your responses.
On Mon, Feb 23, 2015 at 11:37 AM, Matthew Knepley knep...@gmail.com
wrote:
On Mon, Feb 23, 2015 at 9:27 AM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
I'm iterating through local edges given in DMNetworkGetEdgeRange().
For each edge, I extract or modify its corresponding value in a global
petsc vector. Therefore that vector must have as many components as edges
there are in the network. To extract the value in the vector, I use
VecGetArray() and a variable counter that is incremented in each iteration.
The array that I obtain in VecGetArray() has to be the same size than
the edge range. That variable counter starts as 0, so if the array that I
obtained in VecGetArray() is x_array, x_array[0] must be the component
in the global vector that corresponds with the start edge given in
DMNetworkGetEdgeRange()
I need that global petsc vector because I will use it in other
operations, it's not just data. Sorry for the confusion. Thanks in advance.
This sounds like an assembly operation. The usual paradigm is to compute
in the local space, and then communicate to get to the global space. So you
would make a PetscSection that had 1 (or some) unknowns on each cell (edge)
and then you can use DMCreateGlobal/LocalVector() and DMLocalToGlobal() to
do this.
Does that make sense?
Thanks,
Matt
Miguel
On Mon, Feb 23, 2015 at 9:09 AM, Matthew Knepley knep...@gmail.com
wrote:
On Mon, Feb 23, 2015 at 8:42 AM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
Thanks, that will help me. Now what I would like to have is the
following: if I have two processors and ten edges, the partitioning
results
in the first processor having the edges 0-4 and the second processor, the
edges 5-9. I also have a global vector with as many components as edges,
10. How can I partition it so the first processor also has the 0-4
components and the second, the 5-9 components of the vector?
I think it would help to know what you want to accomplish. This is how
you are proposing to do it.'
If you just want to put data on edges, DMNetwork has a facility for
that already.
Thanks,
Matt
Miguel
On Feb 23, 2015 8:08 AM, Abhyankar, Shrirang G.
abhy...@mcs.anl.gov wrote:
Miguel,
One possible way is to store the global numbering of any
edge/vertex in the component attached to it. Once the mesh gets
partitioned, the components are also distributed so you can easily
retrieve
the global number of any edge/vertex by accessing its component. This is
what is done in the DMNetwork example pf.c although the global
numbering is
not used for anything.
Shri
From: Matthew Knepley knep...@gmail.com
Date: Mon, 23 Feb 2015 07:54:34 -0600
To: Miguel Angel Salazar de Troya salazardetr...@gmail.com
Cc: petsc-users@mcs.anl.gov petsc-users@mcs.anl.gov
Subject: Re: [petsc-users] DMNetworkGetEdgeRange() in parallel
On Sun, Feb 22, 2015 at 3:59 PM, Miguel Angel Salazar de Troya
salazardetr...@gmail.com wrote:
Thanks. Once I obtain that
Re: [petsc-users] HermitianTranspose version of MatCreateTranspose.
On Feb 23, 2015, at 3:45 PM, Andrew Spott ansp6...@colorado.edu wrote: It looks like this function doesn’t exist, but it should be pretty easy to write. A few questions: - Does a new MatType need to be created? or will MATTRANSPOSEMAT work? - Is there any way to write this from outside of the PETSc library proper? Initial examination makes it seem like it isn’t possible, because _p_Mat is an incomplete type after the library has been built. The directory include/petsc-private/matimp.h is always provided with a PETSc installation so though _p_Mat is private to PETSc you can extend PETSc out side of the library proper. Barry - Is this a good idea? Or should a MatShell just be made? -Andrew
Re: [petsc-users] HermitianTranspose version of MatCreateTranspose.
Barry Smith bsm...@mcs.anl.gov writes: The directory include/petsc-private/matimp.h is always provided with a PETSc installation so though _p_Mat is private to PETSc you can extend PETSc out side of the library proper. What you don't get when accessing the private headers is a guarantee of backward compatibility (even across subminor releases) or documentation of changes in the release notes. If people are using private headers frequently, we should strive to make a public interface to fulfill their needs. signature.asc Description: PGP signature
Re: [petsc-users] HermitianTranspose version of MatCreateTranspose.
On Feb 23, 2015, at 4:05 PM, Jed Brown j...@jedbrown.org wrote: Barry Smith bsm...@mcs.anl.gov writes: The directory include/petsc-private/matimp.h is always provided with a PETSc installation so though _p_Mat is private to PETSc you can extend PETSc out side of the library proper. What you don't get when accessing the private headers is a guarantee of backward compatibility (even across subminor releases) or documentation of changes in the release notes. If people are using private headers frequently, we should strive to make a public interface to fulfill their needs. If they are writing code that they intend to share with the PETSc community (i.e. eventually it becomes a pull request into petsc so the user will not maintain it themselves for a long time)*; like I assumed in this case, then it is fine to use the private headers. Barry * For example providing a piece of functionality that we frankly just forgot about.
[petsc-users] HermitianTranspose version of MatCreateTranspose.
It looks like this function doesn’t exist, but it should be pretty easy to write. A few questions: - Does a new MatType need to be created? or will MATTRANSPOSEMAT work? - Is there any way to write this from outside of the PETSc library proper? Initial examination makes it seem like it isn’t possible, because _p_Mat is an incomplete type after the library has been built. - Is this a good idea? Or should a MatShell just be made? -Andrew
Re: [petsc-users] HermitianTranspose version of MatCreateTranspose.
I’m definitely willing to submit it as a pull request. Also, while I’m at it, I’m going to write a “duplicate” function for transpose and hermitian_transpose. Just because this seems 1) easy ( MatHermitianTranspose can return a new copy, as well as MatTranspose), and 2) necessary to use these for EPS. Also, is “transpose” a good enough MatType? Or does a new one need to be written? -Andrew On Mon, Feb 23, 2015 at 3:12 PM, Jed Brown j...@jedbrown.org wrote: Barry Smith bsm...@mcs.anl.gov writes: The directory include/petsc-private/matimp.h is always provided with a PETSc installation so though _p_Mat is private to PETSc you can extend PETSc out side of the library proper. What you don't get when accessing the private headers is a guarantee of backward compatibility (even across subminor releases) or documentation of changes in the release notes. If people are using private headers frequently, we should strive to make a public interface to fulfill their needs.
[petsc-users] On user defined PC of schur complement
I have a block matrix [A, B; C, D], and I want to use Schur complement to solve the linear system [A, B; C, D] * [x; y] = [c; d]. I want to apply a preconditioner for solving ( D - C*A^(-1)*B ) y = rhs. And it is easy for me to find an A' which is approximate to A, and A'^(-1) is easy to get. So for this part, I can easily define a preconditioner matrix P = ( D - C*A'^(-1)*B ) using PCFieldSplitSchurPrecondition. The next step is that I want to do incomplete LU factorization for P as my left and right preconditioner for solving ( D - C*A^(-1)*B ) y = rhs using gmres or bcgs. My question is since part of the PC is user defined, and part of it is PETSC defined, how can I combine them? Can I do something like: ierr = PCFieldSplitSchurPrecondition(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, P);CHKERRQ(ierr); ierr = PCSetType(pc, PCILU);CHKERRQ(ierr); Hui
Re: [petsc-users] Solving multiple linear systems with similar matrix sequentially
ML and Hypre both use RAP, as I recall. PtAP uses outer product which cannot be as fast as RAP, based on an analysis of memory access. Sequential RARt might be fast if ARt has small number of colors. We do not have parallel RARt. Hong On Mon, Feb 23, 2015 at 9:47 PM, Barry Smith bsm...@mcs.anl.gov wrote: On Feb 23, 2015, at 9:41 PM, Jed Brown j...@jedbrown.org wrote: Barry Smith bsm...@mcs.anl.gov writes: Freeze the hierarchy and coarse grid interpolations of GAMG but compute the new coarse grid operators RAP for each new linear system (this is a much cheaper operation). Use PCGAMGSetReuseInterpolation() to freeze/unfreeze the hierarchy. The RAP is often more than 50% of PCSetUp, so this might not save much. Hee,hee in some of my recent runs of GAMG the RAP was less than 25% of the time. Thus skipping the other portions could really pay off.* Barry * this could just be because the RAP portion has been optimized much more than the other parts but ...
Re: [petsc-users] HermitianTranspose version of MatCreateTranspose.
We've had a small amount of debate over the years on how to handle the Hermitian transpose and non-Hermitian transpose that never got fully resolved. Approach 1) Each (complex) matrix has a full set of transpose and Hermitian transpose operations (MatTranspose(), MatHermitianTranspose(), MatMultTranspose()), MatMultHermitianTranspose(), MatSolveTranspose(), MatSolveHermitianTranspose(), MatMatMultTranspose(), MatMatMultHermitianTranspose(), MatTranposeMatMult(), MatHermitianTransposeMatMult()...) plus there are two vector inner products; VecDot() and VecTDot(). Approach 2) Consider a (complex) vector (and hence the associated matrix operators on it) to live in the usual Hermitian inner product space or the non-Hermitian inner product space. Then one only needs a single VecDot() and MatTranspose(), MatMultTranspose() ... that just does the right thing based on what space the user has declared the vectors/matrices to be in. Approach 2) seems nicer since it only requires 1/2 the functions :-) and so long as the two vector spaces never interact directly (for example what would be the meaning of the inner product of a vector in the usual Hermitian inner product space with a vector from the non-Hermitian inner product space?) certain seems simpler. Approach 1) might be simpler for some people who like to always see exactly what they are doing. I personally wish I had started with Approach 2 (but I did not), but there could be some flaw with it I am not seeing. Barry On Feb 23, 2015, at 6:50 PM, Andrew Spott ansp6...@colorado.edu wrote: I’m definitely willing to submit it as a pull request. Also, while I’m at it, I’m going to write a “duplicate” function for transpose and hermitian_transpose. Just because this seems 1) easy ( MatHermitianTranspose can return a new copy, as well as MatTranspose), and 2) necessary to use these for EPS. Also, is “transpose” a good enough MatType? Or does a new one need to be written? -Andrew On Mon, Feb 23, 2015 at 3:12 PM, Jed Brown j...@jedbrown.org wrote: signature.asc
Re: [petsc-users] Solving multiple linear systems with similar matrix sequentially
On Feb 23, 2015, at 9:41 PM, Jed Brown j...@jedbrown.org wrote: Barry Smith bsm...@mcs.anl.gov writes: Freeze the hierarchy and coarse grid interpolations of GAMG but compute the new coarse grid operators RAP for each new linear system (this is a much cheaper operation). Use PCGAMGSetReuseInterpolation() to freeze/unfreeze the hierarchy. The RAP is often more than 50% of PCSetUp, so this might not save much. Hee,hee in some of my recent runs of GAMG the RAP was less than 25% of the time. Thus skipping the other portions could really pay off.* Barry * this could just be because the RAP portion has been optimized much more than the other parts but ...
Re: [petsc-users] Solving multiple linear systems with similar matrix sequentially
Barry Smith bsm...@mcs.anl.gov writes: Freeze the hierarchy and coarse grid interpolations of GAMG but compute the new coarse grid operators RAP for each new linear system (this is a much cheaper operation). Use PCGAMGSetReuseInterpolation() to freeze/unfreeze the hierarchy. The RAP is often more than 50% of PCSetUp, so this might not save much. signature.asc Description: PGP signature
Re: [petsc-users] solving multiple linear systems with same matrix (sequentially, not simultaneously)
Daniel Goldberg dngoldb...@gmail.com writes: Thanks Jed. So the nullspace would be initialized with an array of 3 petsc vectors: (u,v)= (0,1), (1,0), and (-y,x), correct? And also to be sure -- this is usefuil only for multigrid preconditioners, yes? Yes, and perhaps also for some exotic domain decomposition methods. signature.asc Description: PGP signature
