From: Derek G. <der...@in...> - 2010-08-31 14:52:48
|
So... up until now I've only ever needed one "ghosted" vector... for the solution vector. And, just as we do in petsc_nonlinear_solver.C I've just been bastardizing System::update() and current_local_solution to "ghostize" a vector. But now that won't work because I need two such vectors simultaneously. What is the current state of the ghosted vector situation? What is the right way to take a Vec from a Petsc callback and form a ghosted PetscVector out of it that I can pass to my own (assembly like) routines? Thanks! Derek |
From: Roy S. <roy...@ic...> - 2010-08-31 16:18:31
|
On Tue, 31 Aug 2010, Derek Gaston wrote: > So... up until now I've only ever needed one "ghosted" vector... for > the solution vector. And, just as we do in petsc_nonlinear_solver.C > I've just been bastardizing System::update() and > current_local_solution to "ghostize" a vector. > > But now that won't work because I need two such vectors > simultaneously. What is the current state of the ghosted vector > situation? What is the right way to take a Vec from a Petsc > callback and form a ghosted PetscVector out of it that I can pass to > my own (assembly like) routines? In theory, assuming your Vec v is already ghosted, all you need is: PetscVector<Number> wrapper(v); and then use wrapper. The constructor PetscVector::PetscVector(Vec v) checks to see whether v is ghosted, serial, or parallel (and how big v is) and updates our metadata cache accordingly. We do assume that v has already been initialized and isn't in mid-assembly, which should be true in your callback function. But that's in theory. In practice: 1. We also call this->_vec = v; this works in current PETSc since a Vec is really just a pointer to the real structure, but I'm not sure if it's official standard behavior that won't break someday. 2. How well tested is this code path? You might be the second or even first guinea pig. --- Roy |
From: Derek G. <der...@in...> - 2010-08-31 16:26:59
|
On Aug 31, 2010, at 10:18 AM, Roy Stogner wrote: > In theory, assuming your Vec v is already ghosted, all you need is: > PetscVector<Number> wrapper(v); > and then use wrapper. The constructor PetscVector::PetscVector(Vec v) > checks to see whether v is ghosted, serial, or parallel (and how big v > is) and updates our metadata cache accordingly. We do assume that v > has already been initialized and isn't in mid-assembly, which should > be true in your callback function. This is what I tried first.... and it failed. I think that one of these vectors is not ghosted already. Just to be clear... this is in a function that I've set as a callback using SNESLineSearchSetPostCheck... and the two vectors I'm working with are the updated solution and the newton increment. Based on your explanation, I'm going to see if one of the vectors is already ghosted... and if it is then it will work fine. The _other_ one I'll use the current_local_solution trick on. This is a bit of a hack though.... > 2. How well tested is this code path? You might be the second or even > first guinea pig. Someone has to be! Derek |
From: Derek G. <der...@in...> - 2010-08-31 18:50:38
|
On Aug 31, 2010, at 10:26 AM, Derek Gaston wrote: > Based on your explanation, I'm going to see if one of the vectors is already ghosted... and if it is then it will work fine. The _other_ one I'll use the current_local_solution trick on. This is a bit of a hack though.... Unfortunately.... this is not the case... they are both coming through with their type set to PARALLEL.... so that plan is out. Any other ideas quickly? Is there an easy way to create a ghosted vector out of a parallel one? If not... I'll just localize them for now. Derek |
From: Jed B. <je...@59...> - 2010-08-31 19:08:43
|
On Tue, 31 Aug 2010 11:18:17 -0500 (CDT), Roy Stogner <roy...@ic...> wrote: > > > On Tue, 31 Aug 2010, Derek Gaston wrote: > > > So... up until now I've only ever needed one "ghosted" vector... for > > the solution vector. And, just as we do in petsc_nonlinear_solver.C > > I've just been bastardizing System::update() and > > current_local_solution to "ghostize" a vector. > > > > But now that won't work because I need two such vectors > > simultaneously. What is the current state of the ghosted vector > > situation? What is the right way to take a Vec from a Petsc > > callback and form a ghosted PetscVector out of it that I can pass to > > my own (assembly like) routines? If a VecGhost was passed in to SNES, then all the Krylov vectors will also be VecGhost. If not, create a new VecGhost (VecCreateGhost() or VecDuplicate() an existing VecGhost) and VecCopy from the plain Vec to the ghost Vec. > But that's in theory. In practice: > > 1. We also call this->_vec = v; this works in current PETSc since a > Vec is really just a pointer to the real structure, but I'm not sure > if it's official standard behavior that won't break someday. I'm not sure what you're after here, but typedef struct _p_Vec *Vec; is pretty fundamental to the object model, it's not subject to change. You might want to PetscObjectReference((PetscObject)v) depending on the owneship semantics you are after. Jed |
From: Derek G. <der...@in...> - 2010-08-31 19:25:14
|
On Aug 31, 2010, at 1:11 PM, Jed Brown wrote: > If a VecGhost was passed in to SNES, then all the Krylov vectors will > also be VecGhost. If not, create a new VecGhost (VecCreateGhost() or > VecDuplicate() an existing VecGhost) and VecCopy from the plain Vec to > the ghost Vec. Well, that is odd then. Because I'm sure that libMesh is passing in ghosted vectors (hmmm... maybe the residual (rhs) isn't ghosted by default?). Maybe libMesh is failing at detecting that they are ghosted? It looks like libMesh is looking for a local to global mapping to decide if it's ghosted or not. Is that the proper way? Derek |
From: Jed B. <je...@59...> - 2010-08-31 20:02:22
|
On Tue, 31 Aug 2010 13:24:44 -0600, Derek Gaston <der...@in...> wrote: > On Aug 31, 2010, at 1:11 PM, Jed Brown wrote: > > > If a VecGhost was passed in to SNES, then all the Krylov vectors > > will also be VecGhost. If not, create a new VecGhost > > (VecCreateGhost() or VecDuplicate() an existing VecGhost) and > > VecCopy from the plain Vec to the ghost Vec. > > Well, that is odd then. Because I'm sure that libMesh is passing in > ghosted vectors (hmmm... maybe the residual (rhs) isn't ghosted by > default?). Maybe libMesh is failing at detecting that they are > ghosted? It looks like libMesh is looking for a local to global > mapping to decide if it's ghosted or not. Is that the proper way? Not really, you can set a LocalToGlobalMapping even for unghosted, in order to use VecSetValuesLocal, though libmesh itself may never use this functionality. Looking at the source, there isn't a public way to check if a Vec has a local form. This could be added, but I'm not convinced it's necessary (or would be good design to use). If you really need it, I can think of two solutions: /* Nasty because this is only available when the source tree is around. */ #include <../src/vec/vec/impls/mpi/pvecimpl.h> if (vec->localrep) { /* is ghosted */ } OR PetscPushErrorHandler(PetscReturnErrorHandler,0); ierr = VecGhostGetLocalForm(X,&loc); PetscPopErrorHandler(); if (!err) { /* is ghosted */ } Note that all serial Vecs are trivially ghosted (the updates are no-ops). > Jed: how much overhead will be incurred by having all vectors in my > SNES system be ghosted? Will it be a noticeable slowdown or > communication? No extra communication, the only price is the extra memory for the ghosted part that isn't being used. That could be significant for very small or poorly shaped subdomains. PETSc could and maybe should deal with this by creating the Krylov space with non-ghosted vectors, but we currently always use vanilla VecDuplicate() so that the vector operations will be done with the appropriate types (e.g. VecCUDA). > Yeah; the real trouble is that ghosted vectors have to fall back to > serial when ghosting isn't available, e.g. on Trilinos. I'm going to > assume that you don't want to implement ghosted support there? I assume you mean working with a separate vector for the local form. > Otherwise the only thing to do is add a non-default code path, where > instead of letting PETSc work on solution while we work to keep > current_local_solution sync'ed up we do the opposite. PETSc is always going to work with a parallel vector (VecGhost is parallel). The different code path is that with VecGhost, you would normally VecGhostUpdate then VecGhostGetLocalForm instead of VecScatter to an independent local vector. Are you talking about this or something else? Jed |
From: Roy S. <roy...@ic...> - 2010-08-31 19:44:34
|
On Tue, 31 Aug 2010, Derek Gaston wrote: > Well, that is odd then. Because I'm sure that libMesh is passing in > ghosted vectors (hmmm... maybe the residual (rhs) isn't ghosted by > default?). Actually, the only vectors which currently default to ghosted are old_local_solution, older_local_solution in TransientSystems, old_local_nonlinear_solution in unsteady DiffSystems, and current_local_solution in any System. But we pass solution and rhs to the petsc nonlinear solver, and both of those are just plain parallel. --- Roy |
From: Derek G. <der...@in...> - 2010-08-31 19:47:45
|
On Aug 31, 2010, at 1:44 PM, Roy Stogner wrote: > Actually, the only vectors which currently default to ghosted are > old_local_solution, older_local_solution in TransientSystems, > old_local_nonlinear_solution in unsteady DiffSystems, and > current_local_solution in any System. But we pass solution and rhs to > the petsc nonlinear solver, and both of those are just plain parallel. That would explain it. Is there any way to modify that behavior... if I'm the only one that needs a fully ghosted system I don't want to slow everyone else down. Jed: how much overhead will be incurred by having all vectors in my SNES system be ghosted? Will it be a noticeable slowdown or communication? Derek |
From: Roy S. <roy...@ic...> - 2010-08-31 19:51:28
|
On Tue, 31 Aug 2010, Derek Gaston wrote: > > On Aug 31, 2010, at 1:44 PM, Roy Stogner wrote: > >> Actually, the only vectors which currently default to ghosted are >> old_local_solution, older_local_solution in TransientSystems, >> old_local_nonlinear_solution in unsteady DiffSystems, and >> current_local_solution in any System. But we pass solution and rhs to >> the petsc nonlinear solver, and both of those are just plain parallel. > > That would explain it. Is there any way to modify that behavior... > if I'm the only one that needs a fully ghosted system I don't want > to slow everyone else down. Yeah; the real trouble is that ghosted vectors have to fall back to serial when ghosting isn't available, e.g. on Trilinos. I'm going to assume that you don't want to implement ghosted support there? Otherwise the only thing to do is add a non-default code path, where instead of letting PETSc work on solution while we work to keep current_local_solution sync'ed up we do the opposite. --- Roy |
From: Derek G. <der...@in...> - 2010-08-31 20:03:47
|
On Aug 31, 2010, at 1:51 PM, Roy Stogner wrote: > Yeah; the real trouble is that ghosted vectors have to fall back to > serial when ghosting isn't available, e.g. on Trilinos. I'm going to > assume that you don't want to implement ghosted support there? We do. And I actually have someone assigned to it... but he has also been assigned other things with higher priorities..... > Otherwise the only thing to do is add a non-default code path, where > instead of letting PETSc work on solution while we work to keep > current_local_solution sync'ed up we do the opposite. That's an option... but what about just making a compile-time option that makes NumericVector::build() create ghosted vectors by default? Derek |
From: Roy S. <roy...@ic...> - 2010-08-31 20:13:12
|
On Tue, 31 Aug 2010, Derek Gaston wrote: > On Aug 31, 2010, at 1:51 PM, Roy Stogner wrote: > >> Otherwise the only thing to do is add a non-default code path, where >> instead of letting PETSc work on solution while we work to keep >> current_local_solution sync'ed up we do the opposite. > > That's an option... but what about just making a compile-time option > that makes NumericVector::build() create ghosted vectors by default? We don't want to create ghosted instead of serial vectors - sometimes code is going to ask for a serial vector because it really needs data that's not on the local processor or its ghost cells. (and because it's too lazy to build a proper extended send list...) But it would be sufficient for your case to just create ghosted instead of parallel vectors. We can't do that in general in NumericVector::build() because that doesn't have access to DofMap's send_list. But we could do that as an option (run-time, even - make it a flag in the System class) in System::add_vector() and *System::init(). --- Roy |
From: Jed B. <je...@59...> - 2010-08-31 20:19:18
|
On Tue, 31 Aug 2010 15:13:02 -0500 (CDT), Roy Stogner <roy...@ic...> wrote: > We don't want to create ghosted instead of serial vectors - sometimes > code is going to ask for a serial vector because it really needs data > that's not on the local processor or its ghost cells. (and because > it's too lazy to build a proper extended send list...) That wouldn't make sense because ghosted vectors are parallel vectors, not serial vectors. > But it would be sufficient for your case to just create ghosted > instead of parallel vectors. We can't do that in general in > NumericVector::build() because that doesn't have access to DofMap's > send_list. But we could do that as an option (run-time, even - make > it a flag in the System class) in System::add_vector() and > *System::init(). I don't remember how this is done in libmesh, but it would be typical to create a solution/residual vector from the function space (mesh), and the overlap (distribution and ghost indices) is available in that context. Jed |
From: Roy S. <roy...@ic...> - 2010-08-31 20:25:04
|
On Tue, 31 Aug 2010, Jed Brown wrote: >> But it would be sufficient for your case to just create ghosted >> instead of parallel vectors. We can't do that in general in >> NumericVector::build() because that doesn't have access to DofMap's >> send_list. But we could do that as an option (run-time, even - make >> it a flag in the System class) in System::add_vector() and >> *System::init(). > > I don't remember how this is done in libmesh, but it would be typical to > create a solution/residual vector from the function space (mesh), and > the overlap (distribution and ghost indices) is available in that > context. Yes; that's the context available in the System class. However when we're just building a parallel vector, we only call NumericVector::build() with the global and local sizes, not with the entire context that exists in the caller. The best solution is to add a special case alternative that uses the ghosted version of build() for those half dozen calls, not to try and figure out how to break our own encapsulation within the non-ghosted build() call. --- Roy |
From: Jed B. <je...@59...> - 2010-08-31 20:41:10
|
On Tue, 31 Aug 2010 15:24:55 -0500 (CDT), Roy Stogner <roy...@ic...> wrote: > > On Tue, 31 Aug 2010, Jed Brown wrote: > > >> But it would be sufficient for your case to just create ghosted > >> instead of parallel vectors. We can't do that in general in > >> NumericVector::build() because that doesn't have access to DofMap's > >> send_list. But we could do that as an option (run-time, even - make > >> it a flag in the System class) in System::add_vector() and > >> *System::init(). > > > > I don't remember how this is done in libmesh, but it would be typical to > > create a solution/residual vector from the function space (mesh), and > > the overlap (distribution and ghost indices) is available in that > > context. > > Yes; that's the context available in the System class. > > However when we're just building a parallel vector, we only call > NumericVector::build() with the global and local sizes, not with the > entire context that exists in the caller. Okay, I just peeked at system.C. I think the way I would implement is to add the ghosted interface to your NumericVector class (when a send-list is provided). When running with the PETSc backend, there is only one slot (for the ghosted parallel vector), and getting the local form calls VecGhostGetLocalForm. With other backends, or if the user explicitly doesn't want to use ghosted, then NumericVector would perform the scatter and give the user access to the local vector. Then you could get rid of all ghost-related special cases in system.C. Jed |
From: Derek G. <fri...@gm...> - 2010-08-31 20:26:40
|
On Aug 31, 2010, at 2:24 PM, Roy Stogner wrote: > However when we're just building a parallel vector, we only call > NumericVector::build() with the global and local sizes, not with the > entire context that exists in the caller. The best solution is to add > a special case alternative that uses the ghosted version of build() > for those half dozen calls, not to try and figure out how to break our > own encapsulation within the non-ghosted build() call. I'm open for whatever solution. If you tell me what you want I'll take a crack at an implementation. Derek |
From: Roy S. <roy...@ic...> - 2010-08-31 20:37:04
|
On Tue, 31 Aug 2010, Derek Gaston wrote: > On Aug 31, 2010, at 2:24 PM, Roy Stogner wrote: > >> However when we're just building a parallel vector, we only call >> NumericVector::build() with the global and local sizes, not with the >> entire context that exists in the caller. The best solution is to add >> a special case alternative that uses the ghosted version of build() >> for those half dozen calls, not to try and figure out how to break our >> own encapsulation within the non-ghosted build() call. > > I'm open for whatever solution. If you tell me what you want I'll take a crack at an implementation. The effect you want should be to add a get_send_list() argument and change PARALLEL to GHOSTED in system.C lines 230, 246, 262, 310, 593. Try that first and see if it works. If it turns out to be what you need, then we could add it as an option to the library; maybe create a new System::init_parallel_vector() method that gets called on those lines and either inits a PARALLEL by default or a GHOSTED if some System::all_ghosted_vectors boolean is true? --- Roy |
From: Derek G. <fri...@gm...> - 2010-08-31 20:38:06
|
On Aug 31, 2010, at 2:36 PM, Roy Stogner wrote: > The effect you want should be to add a get_send_list() argument and > change PARALLEL to GHOSTED in system.C lines 230, 246, 262, 310, 593. > Try that first and see if it works. If it turns out to be what you > need, then we could add it as an option to the library; maybe create a > new System::init_parallel_vector() method that gets called on those > lines and either inits a PARALLEL by default or a GHOSTED if some > System::all_ghosted_vectors boolean is true? Great - I'll give it a shot. Derek |
From: Derek G. <fri...@gm...> - 2010-09-01 17:39:49
|
On Aug 31, 2010, at 2:36 PM, Roy Stogner wrote: > The effect you want should be to add a get_send_list() argument and > change PARALLEL to GHOSTED in system.C lines 230, 246, 262, 310, 593. > Try that first and see if it works. If it turns out to be what you > need, then we could add it as an option to the library; maybe create a > new System::init_parallel_vector() method that gets called on those > lines and either inits a PARALLEL by default or a GHOSTED if some > System::all_ghosted_vectors boolean is true? Ok - I made these changes... and everything works for me. I'm going to take your suggestion and implement init_parallel_vector().... with the bool you suggested unless someone has objections. Derek |
From: Roy S. <roy...@ic...> - 2010-09-01 17:58:19
|
On Wed, 1 Sep 2010, Derek Gaston wrote: > On Aug 31, 2010, at 2:36 PM, Roy Stogner wrote: > >> The effect you want should be to add a get_send_list() argument and >> change PARALLEL to GHOSTED in system.C lines 230, 246, 262, 310, 593. >> Try that first and see if it works. If it turns out to be what you >> need, then we could add it as an option to the library; maybe create a >> new System::init_parallel_vector() method that gets called on those >> lines and either inits a PARALLEL by default or a GHOSTED if some >> System::all_ghosted_vectors boolean is true? > > Ok - I made these changes... and everything works for me. I'm going > to take your suggestion and implement init_parallel_vector().... > with the bool you suggested unless someone has objections. I'm not sure I understood Jed's suggestion - eliminate parallel non-ghosted vectors entirely, then write our own ghosting support to handle the non-PETSc backends? That seems like a lot of work to get slightly more robust support for a feature that most people wouldn't need. We already support non-ghosted vectors, we already support ghosted vectors (efficiently with PETSc or inefficiently otherwise), why not just let people choose between them even for the "internal" vectors? My single-bool suggestion isn't the best way to do things. Ideally any time we're reinitializing an existing vector (lines 246, 262, 310) we want to query it and reinitialize it with the same type it already had. We could use a bool (or better, an enum ParallelType) to determine how to initialize everything - perhaps a System member for System::solution and as another optional argument to System::add_vector() for others. But you've been sounding like you're in a hurry; if you want to start with the single-bool now and get things just right after whatever deadline is past, that's fine too. --- Roy |
From: Derek G. <fri...@gm...> - 2010-09-01 19:02:07
|
On Sep 1, 2010, at 11:58 AM, Roy Stogner wrote: > My single-bool suggestion isn't the best way to do things. Ideally > any time we're reinitializing an existing vector (lines 246, 262, 310) > we want to query it and reinitialize it with the same type it already > had. We could use a bool (or better, an enum ParallelType) to > determine how to initialize everything - perhaps a System member for > System::solution and as another optional argument to > System::add_vector() for others. Essentially, you're suggesting a std::map between vector name and ParallelType? That way each vector can be reinitialized with different types? Then you can pass ParallelType into add_vector()? And like you mention you'll need member variables for the parallel type of "solution" and "rhs". Sounds reasonable to me. I'll get on it. > But you've been sounding like you're in a hurry; if you want to start > with the single-bool now and get things just right after whatever > deadline is past, that's fine too. Well... I'm always in a hurry! That's no excuse for hacking terrible code! ;-) Derek |
From: Roy S. <roy...@ic...> - 2010-09-01 19:22:01
|
On Wed, 1 Sep 2010, Derek Gaston wrote: > Essentially, you're suggesting a std::map between vector name and > ParallelType? No need - just check NumericVector::type() for each. > Then you can pass ParallelType into add_vector()? Right. > And like you mention you'll need member variables for the parallel > type of "solution" and "rhs". rhs (and adjoint_solution, and everything other than solution) currently gets created via add_vector, so if anyone wants to enforce non-default vector behavior they can just create the vector themselves before ExplicitSystem::init (or ImplicitSystem::adjoint_solve, whatever) gets to it. Solution is now kind of an ugly special-case. Originally we had several vectors that weren't created via add_vector, and I changed that, but now I don't remember why I didn't change solution along with them. Maybe now's the right time to do that? And then we don't need any member variables - the user can just do an add_vector() with whatever options they prefer. We'd also get to (have to!) simplify a lot of code which currently looks like "do X to solution, then loop over all other vectors and do X to each"; it would instead just "loop over all vectors and do X to each" --- Roy |
From: Jed B. <je...@59...> - 2010-09-01 19:51:09
|
On Wed, 1 Sep 2010 12:58:11 -0500 (CDT), Roy Stogner <roy...@ic...> wrote: > I'm not sure I understood Jed's suggestion - eliminate parallel > non-ghosted vectors entirely, then write our own ghosting support to > handle the non-PETSc backends? That seems like a lot of work to get > slightly more robust support for a feature that most people wouldn't > need. We already support non-ghosted vectors, we already support > ghosted vectors (efficiently with PETSc or inefficiently otherwise), > why not just let people choose between them even for the "internal" > vectors? Abstractly, there are only two "views" that a user needs: the parallel one where each degree of freedom appears exactly once, and the "local" one where a single process can look at it's subdomain plus unowned interface dofs. With VecGhost, these get to share the same array and you use VecGhostUpdate to update ghost values and sum into global values, otherwise they are separate arrays and you map between with VecScatter. Since you have your own Vector class, you can add an extra slot for the local form and update it with VecScatter instead of using the PETSc-managed local form that would be updated with VecGhostUpdate. If this whole thing is more than about ten lines of code, something is wrong. But maybe I've missed something unique about your use case. Jed |
From: Roy S. <roy...@ic...> - 2010-09-01 22:55:59
|
On Wed, 1 Sep 2010, Jed Brown wrote: > On Wed, 1 Sep 2010 12:58:11 -0500 (CDT), Roy Stogner <roy...@ic...> wrote: >> I'm not sure I understood Jed's suggestion - eliminate parallel >> non-ghosted vectors entirely, then write our own ghosting support to >> handle the non-PETSc backends? That seems like a lot of work to get >> slightly more robust support for a feature that most people wouldn't >> need. We already support non-ghosted vectors, we already support >> ghosted vectors (efficiently with PETSc or inefficiently otherwise), >> why not just let people choose between them even for the "internal" >> vectors? > > Abstractly, there are only two "views" that a user needs: the parallel > one where each degree of freedom appears exactly once, and the "local" > one where a single process can look at it's subdomain plus unowned > interface dofs. With VecGhost, these get to share the same array and > you use VecGhostUpdate to update ghost values and sum into global > values, otherwise they are separate arrays and you map between with > VecScatter. Since you have your own Vector class, you can add an extra > slot for the local form and update it with VecScatter instead of using > the PETSc-managed local form that would be updated with VecGhostUpdate. > If this whole thing is more than about ten lines of code, something is > wrong. > > But maybe I've missed something unique about your use case. There are indeed a lot of vectors we end up using for which user code doesn't need access to ghost data, and for which we wouldn't necessarily want default storage of ghost degrees of freedom. On the other hand: Tim's initial interface to PETSc ghost vectors was a couple hundred lines. If you can support ghost data on all NumericVector subclasses together in 10 lines, we'd love a patch! (this was actually the way Ben had originally intended to add ghosted vector support - but when Tim offered to write the PETSc-specific support, we'd hardly want to turn that down in favor of something else we might write someday) --- Roy |
From: Derek G. <fri...@gm...> - 2010-09-01 22:14:51
|
On Sep 1, 2010, at 1:21 PM, Roy Stogner wrote: > rhs (and adjoint_solution, and everything other than solution) > currently gets created via add_vector, so if anyone wants to enforce > non-default vector behavior they can just create the vector themselves > before ExplicitSystem::init (or ImplicitSystem::adjoint_solve, > whatever) gets to it. I would say that this is at best inconvenient... and at worst breaks object-oriented ideas. Users shouldn't have to know the names and exactly what vectors each type of solver creates. > Solution is now kind of an ugly special-case. Originally we had > several vectors that weren't created via add_vector, and I changed > that, but now I don't remember why I didn't change solution along with > them. Maybe now's the right time to do that? And then we don't need > any member variables - the user can just do an add_vector() with > whatever options they prefer. We'd also get to (have to!) simplify a > lot of code which currently looks like "do X to solution, then loop > over all other vectors and do X to each"; it would instead just "loop > over all vectors and do X to each" I'll look into adding Solution with add_vector... that sounds like the right way to go. So my compromise is this: add a member variable to system called "default_parallel_vector_type" that defaults to PARALLEL. add_vector will get an optional argument for the ParallelType (defaulting to PARALLEL). add_vector will also store off the ParallelType for each vector in a map (I _do_ think this is necessary because there is no other way to capture the ParallelType within add_vector). All internal vectors will be added using add_vector and passing in "default_parallel_vector_type" (or GHOSTED or PARALLEL if only those make sense... like for current_local_solution). Finally, I'm also going to create a helper function in System called init_vector that just takes the name of the vector to be inited... this is because there are a few places where vectors get explicitly inited. I think that should cover it. If you _want_ to you will be able to override the defaults by either setting default_parallel_vector_type OR by doing as you mentioned and calling add_vector yourself with a different type. Derek |