You can subscribe to this list here.
2003 
_{Jan}

_{Feb}

_{Mar}

_{Apr}

_{May}

_{Jun}

_{Jul}

_{Aug}

_{Sep}
(2) 
_{Oct}
(2) 
_{Nov}
(27) 
_{Dec}
(31) 

2004 
_{Jan}
(6) 
_{Feb}
(15) 
_{Mar}
(33) 
_{Apr}
(10) 
_{May}
(46) 
_{Jun}
(11) 
_{Jul}
(21) 
_{Aug}
(15) 
_{Sep}
(13) 
_{Oct}
(23) 
_{Nov}
(1) 
_{Dec}
(8) 
2005 
_{Jan}
(27) 
_{Feb}
(57) 
_{Mar}
(86) 
_{Apr}
(23) 
_{May}
(37) 
_{Jun}
(34) 
_{Jul}
(24) 
_{Aug}
(17) 
_{Sep}
(50) 
_{Oct}
(24) 
_{Nov}
(10) 
_{Dec}
(60) 
2006 
_{Jan}
(47) 
_{Feb}
(46) 
_{Mar}
(127) 
_{Apr}
(19) 
_{May}
(26) 
_{Jun}
(62) 
_{Jul}
(47) 
_{Aug}
(51) 
_{Sep}
(61) 
_{Oct}
(42) 
_{Nov}
(50) 
_{Dec}
(33) 
2007 
_{Jan}
(60) 
_{Feb}
(55) 
_{Mar}
(77) 
_{Apr}
(102) 
_{May}
(82) 
_{Jun}
(102) 
_{Jul}
(169) 
_{Aug}
(117) 
_{Sep}
(80) 
_{Oct}
(37) 
_{Nov}
(51) 
_{Dec}
(43) 
2008 
_{Jan}
(71) 
_{Feb}
(94) 
_{Mar}
(98) 
_{Apr}
(125) 
_{May}
(54) 
_{Jun}
(119) 
_{Jul}
(60) 
_{Aug}
(111) 
_{Sep}
(118) 
_{Oct}
(125) 
_{Nov}
(119) 
_{Dec}
(94) 
2009 
_{Jan}
(109) 
_{Feb}
(38) 
_{Mar}
(93) 
_{Apr}
(88) 
_{May}
(29) 
_{Jun}
(57) 
_{Jul}
(53) 
_{Aug}
(48) 
_{Sep}
(68) 
_{Oct}
(151) 
_{Nov}
(23) 
_{Dec}
(35) 
2010 
_{Jan}
(84) 
_{Feb}
(60) 
_{Mar}
(184) 
_{Apr}
(112) 
_{May}
(60) 
_{Jun}
(90) 
_{Jul}
(23) 
_{Aug}
(70) 
_{Sep}
(119) 
_{Oct}
(27) 
_{Nov}
(47) 
_{Dec}
(54) 
2011 
_{Jan}
(22) 
_{Feb}
(19) 
_{Mar}
(92) 
_{Apr}
(93) 
_{May}
(35) 
_{Jun}
(91) 
_{Jul}
(32) 
_{Aug}
(61) 
_{Sep}
(7) 
_{Oct}
(69) 
_{Nov}
(81) 
_{Dec}
(23) 
2012 
_{Jan}
(64) 
_{Feb}
(95) 
_{Mar}
(35) 
_{Apr}
(36) 
_{May}
(63) 
_{Jun}
(98) 
_{Jul}
(70) 
_{Aug}
(171) 
_{Sep}
(149) 
_{Oct}
(64) 
_{Nov}
(67) 
_{Dec}
(126) 
2013 
_{Jan}
(108) 
_{Feb}
(104) 
_{Mar}
(171) 
_{Apr}
(133) 
_{May}
(108) 
_{Jun}
(100) 
_{Jul}
(93) 
_{Aug}
(126) 
_{Sep}
(74) 
_{Oct}
(59) 
_{Nov}
(145) 
_{Dec}
(93) 
2014 
_{Jan}
(38) 
_{Feb}
(45) 
_{Mar}
(26) 
_{Apr}
(41) 
_{May}
(125) 
_{Jun}
(70) 
_{Jul}
(61) 
_{Aug}
(66) 
_{Sep}
(60) 
_{Oct}
(110) 
_{Nov}
(27) 
_{Dec}
(30) 
2015 
_{Jan}
(43) 
_{Feb}
(67) 
_{Mar}
(71) 
_{Apr}
(92) 
_{May}
(39) 
_{Jun}
(15) 
_{Jul}
(46) 
_{Aug}
(63) 
_{Sep}
(84) 
_{Oct}
(82) 
_{Nov}
(69) 
_{Dec}
(45) 
2016 
_{Jan}
(92) 
_{Feb}
(91) 
_{Mar}
(148) 
_{Apr}
(43) 
_{May}
(58) 
_{Jun}
(117) 
_{Jul}
(92) 
_{Aug}
(140) 
_{Sep}
(49) 
_{Oct}
(33) 
_{Nov}
(85) 
_{Dec}
(40) 
2017 
_{Jan}
(41) 
_{Feb}
(36) 
_{Mar}
(49) 
_{Apr}
(41) 
_{May}
(68) 
_{Jun}

_{Jul}

_{Aug}

_{Sep}

_{Oct}

_{Nov}

_{Dec}

S  M  T  W  T  F  S 




1

2

3
(3) 
4
(3) 
5
(5) 
6
(5) 
7
(1) 
8

9
(14) 
10
(10) 
11

12
(1) 
13
(13) 
14
(2) 
15
(22) 
16
(4) 
17
(5) 
18
(1) 
19
(1) 
20
(5) 
21
(3) 
22
(5) 
23
(2) 
24

25

26

27
(2) 
28
(2) 
29
(5) 
30
(5) 


From: Jed Brown <jed@59...>  20100915 15:45:22

On Wed, 15 Sep 2010 17:34:14 +0200 (CEST), Tim Kroeger <tim.kroeger@...> wrote: > But now, let's assume > > all = [ 0 .. 9  10 .. 19  20 .. 29 ] > > and > > sub = [ 8 9  10 .. 19  20 21 22 ] > > Assume further that the matrix is tridiagonal. If I then do what you > suggested, four matrix entries are lost (that is, (9,10), (10,9), > (19,20), (20,19)), aren't they? That's why I thought the IS for the > columns should be larger than that for the rows. If the user provided > dof indices are only those already owned by the processor, it seems to > be necessary to communicate those across the processors in order to do > the correct column selection. Is this right? No, the local part of the column IS denotes the column indices that will belong to that "diagonal block". To put it differently, it defines the ownership of a vector that the matrix can be multiplied by. Suppose we want the equation y = A x to make sense. The "row IS" specifies the distribution of y, the "column IS" specifies the distribution of x. The diagonal block of A is the tensor product of the local part of the row and column ISs. > I don't find this weird. Well, we could let the user decide between > several modes here, that is KEEP_SOLUTION, CLEAR, TAKE_RHS (or > whatever names we agree to). Roy, what do you think? Usually a (tightly converged) solve is a linear operation from the righthand side to a solution, the values in the solution vector only specify the starting point for the iteration. If tight tolerances are used, then the starting vector only affects the time to solution, but not the actual value. Solving on a subdomain without modifying the solution vector for "inactive" dofs is a deeply different semantic. Jed 
From: Tim Kroeger <tim.kroeger@ce...>  20100915 15:34:34

On Wed, 15 Sep 2010, Jed Brown wrote: > On Wed, 15 Sep 2010 16:49:53 +0200 (CEST), Tim Kroeger <tim.kroeger@...> wrote: > >>> VecScatter moves entries around, and isn't the issue here. >>> Redistributing indices to satisfy some apriori size bound is somewhat >>> different (though not hard). >> >> So do I have to do something additional or not? I would think that I >> don't have to, because VecSetSizes() should be able to divide the >> index range {0,...,N1} into k (almost) equally wide intervals itself. > > The user provides a *different* std::vector on each process (the global > indices that will be owned by that process in the subproblem). Ah, that's the point. I was thinking of the situation where *the same* std::vector is used on all processors, and it hence consists of *all* global dofs in the subset. When I now think about this, I seems as if in fact your way is more sensible, in particular also for the reason that setting up the vector becomes much cheaper by this. > You only duplicate the whole matrix everywhere if every process asks > for the global problem (and the resulting matrix is weird, and > almost certainly not what you want). In the easy (static) case, > each process only asks for global indices that it already owns. If > you want to repartition these indices, then you need to do that (it > won't happen automatically). > > That is, suppose the full problem of dimension 30 is partitioned regularly > > all = [0 .. 9  10 .. 19  20 .. 29] > > Now there is some subdomain > > sub_static = [2 3  12 13 14 15  25] > > If use this IS, the subproblem will not be very well balanced, so you > might use > > sub_balanced = [2 3 12  13 14  15 25] > > But you are responsible for creating sub_balanced, it doesn't just > magically happen. Yes, I understand this now, thank you very much. But now, let's assume all = [ 0 .. 9  10 .. 19  20 .. 29 ] and sub = [ 8 9  10 .. 19  20 21 22 ] Assume further that the matrix is tridiagonal. If I then do what you suggested, four matrix entries are lost (that is, (9,10), (10,9), (19,20), (20,19)), aren't they? That's why I thought the IS for the columns should be larger than that for the rows. If the user provided dof indices are only those already owned by the processor, it seems to be necessary to communicate those across the processors in order to do the correct column selection. Is this right? >>> You probably want to set the rest of solution>vec() to the right side, >>> or zero it out, otherwise a solve is not a linear operation. The reason >>> for the first is that it makes the preconditioner look like the >>> submatrix was extended to operate on the full domain using the identity >>> (instead of zero which would make it very singular). >> >> You're right in that my current implementation makes solve not be a >> linear operation if considered on the complete vector. But what I >> want is some operation that only works on a subset of the dofs, where >> the same subset applies to the rhs and the solution. On this subset, >> solve is a linear operation, and the dofs outside the subset just >> don't take part in the operation at all, that is, they are kept >> constant. At least, in my application that's what I need. > > That makes this thing a weird beast. Are you sure it belongs in the > library? I don't find this weird. Well, we could let the user decide between several modes here, that is KEEP_SOLUTION, CLEAR, TAKE_RHS (or whatever names we agree to). Roy, what do you think? Best Regards, Tim  Dr. Tim Kroeger CeVis  Center of Complex Systems and Visualization University of Bremen tim.kroeger@... Universitaetsallee 29 tim.kroeger@... D28359 Bremen Phone +494212187710 Germany Fax +494212184236 
From: Jed Brown <jed@59...>  20100915 15:28:08

On Wed, 15 Sep 2010 09:50:25 0500 (CDT), Roy Stogner <roystgnr@...> wrote: > In the topological context I thought a domain had to be open and > connected It is always open in the subset topology, regarless of whether the subset is open and/or closed in the parent space. No need for connectedness either, x > sqrt(x^2  1), working with real values, has a perfectly nice domain (\infty,1] \cup [1,\infty). The term is also used in category theory where there is no topology either. I don't care in the least about this naming. Jed 
From: Jed Brown <jed@59...>  20100915 15:09:24

On Wed, 15 Sep 2010 16:49:53 +0200 (CEST), Tim Kroeger <tim.kroeger@...> wrote: > > VecScatter moves entries around, and isn't the issue here. > > Redistributing indices to satisfy some apriori size bound is somewhat > > different (though not hard). > > So do I have to do something additional or not? I would think that I > don't have to, because VecSetSizes() should be able to divide the > index range {0,...,N1} into k (almost) equally wide intervals itself. The user provides a *different* std::vector on each process (the global indices that will be owned by that process in the subproblem). You only duplicate the whole matrix everywhere if every process asks for the global problem (and the resulting matrix is weird, and almost certainly not what you want). In the easy (static) case, each process only asks for global indices that it already owns. If you want to repartition these indices, then you need to do that (it won't happen automatically). That is, suppose the full problem of dimension 30 is partitioned regularly all = [0 .. 9  10 .. 19  20 .. 29] Now there is some subdomain sub_static = [2 3  12 13 14 15  25] If use this IS, the subproblem will not be very well balanced, so you might use sub_balanced = [2 3 12  13 14  15 25] But you are responsible for creating sub_balanced, it doesn't just magically happen. > > You probably want to set the rest of solution>vec() to the right side, > > or zero it out, otherwise a solve is not a linear operation. The reason > > for the first is that it makes the preconditioner look like the > > submatrix was extended to operate on the full domain using the identity > > (instead of zero which would make it very singular). > > You're right in that my current implementation makes solve not be a > linear operation if considered on the complete vector. But what I > want is some operation that only works on a subset of the dofs, where > the same subset applies to the rhs and the solution. On this subset, > solve is a linear operation, and the dofs outside the subset just > don't take part in the operation at all, that is, they are kept > constant. At least, in my application that's what I need. That makes this thing a weird beast. Are you sure it belongs in the library? Jed 
From: Tim Kroeger <tim.kroeger@ce...>  20100915 14:55:09

On Wed, 15 Sep 2010, Roy Stogner wrote: > On Wed, 15 Sep 2010, Jed Brown wrote: > >> On Wed, 15 Sep 2010 09:21:30 0500 (CDT), Roy Stogner >> <roystgnr@...> wrote: >>> but "domain" is sort of PDE specific whereas LinearSolver is not. >> >> FWIW, "domain" is topological (more general than linear or even metric >> spaces). > > In the topological context I thought a domain had to be open and > connected, and I'm not sure how that applies here. E.g. {1,2,3} isn't > a connected subset of {1,2,3,4,5} under the discrete topology, and R^3 > isn't an open subset of R^5 under the norm topology. The subsets that I will be working in in my application will neither be open (because they include their boundary) nor necessarily connected (may consist of several connected components). Anyway, I will rename the methods to restrict_solve_to(), which, I think, was one of Roy's suggestions. Best Regards, Tim  Dr. Tim Kroeger CeVis  Center of Complex Systems and Visualization University of Bremen tim.kroeger@... Universitaetsallee 29 tim.kroeger@... D28359 Bremen Phone +494212187710 Germany Fax +494212184236 
From: Roy Stogner <roystgnr@ic...>  20100915 14:50:34

On Wed, 15 Sep 2010, Jed Brown wrote: > On Wed, 15 Sep 2010 09:21:30 0500 (CDT), Roy Stogner <roystgnr@...> wrote: >> but "domain" is sort of PDE specific whereas LinearSolver is not. > > FWIW, "domain" is topological (more general than linear or even metric > spaces). In the topological context I thought a domain had to be open and connected, and I'm not sure how that applies here. E.g. {1,2,3} isn't a connected subset of {1,2,3,4,5} under the discrete topology, and R^3 isn't an open subset of R^5 under the norm topology.  Roy 
From: Tim Kroeger <tim.kroeger@ce...>  20100915 14:50:07

On Wed, 15 Sep 2010, Jed Brown wrote: > On Wed, 15 Sep 2010 14:53:18 +0200 (CEST), Tim Kroeger <tim.kroeger@...> wrote: > >>> This chooses the sizes of each part, but you will still have to move >>> indices around. >> >> This is done by VecScatter, isn't it? I think that the overhead of >> doing this once before and after the solve will be small against the >> gain in performance by having balanced vectors. Well, of course this >> depends on a lot of things, but in the mean... > > VecScatter moves entries around, and isn't the issue here. > Redistributing indices to satisfy some apriori size bound is somewhat > different (though not hard). So do I have to do something additional or not? I would think that I don't have to, because VecSetSizes() should be able to divide the index range {0,...,N1} into k (almost) equally wide intervals itself. >> + size_t _solve_only_on_is_size(void)const; > > Note that you only call this once. Right for now, but I will later call it once in every overloaded solve() method. >> + ierr = MatGetSubMatrix(matrix>mat(), >> + _solve_only_on_is_local,_solve_only_on_is, >> + PETSC_DECIDE,MAT_INITIAL_MATRIX,&submat); >> + CHKERRABORT(libMesh::COMM_WORLD,ierr); > > Why do you want this nonsquare matrix? I'm pretty sure you want to > extract the symmetric block (_solve_only_on_is,_solve_only_on_is). That's what I had here before, but as far as I understood you before, this would result in duplicating the complete matrix over all processors. >> + ierr = VecScatterBegin(scatter,subsolution,solution>vec(),INSERT_VALUES,SCATTER_REVERSE); >> + CHKERRABORT(libMesh::COMM_WORLD,ierr); >> + ierr = VecScatterEnd(scatter,subsolution,solution>vec(),INSERT_VALUES,SCATTER_REVERSE); >> + CHKERRABORT(libMesh::COMM_WORLD,ierr); > > You probably want to set the rest of solution>vec() to the right side, > or zero it out, otherwise a solve is not a linear operation. The reason > for the first is that it makes the preconditioner look like the > submatrix was extended to operate on the full domain using the identity > (instead of zero which would make it very singular). You're right in that my current implementation makes solve not be a linear operation if considered on the complete vector. But what I want is some operation that only works on a subset of the dofs, where the same subset applies to the rhs and the solution. On this subset, solve is a linear operation, and the dofs outside the subset just don't take part in the operation at all, that is, they are kept constant. At least, in my application that's what I need. Best Regards, Tim  Dr. Tim Kroeger CeVis  Center of Complex Systems and Visualization University of Bremen tim.kroeger@... Universitaetsallee 29 tim.kroeger@... D28359 Bremen Phone +494212187710 Germany Fax +494212184236 
From: Jed Brown <jed@59...>  20100915 14:34:09

On Wed, 15 Sep 2010 09:21:30 0500 (CDT), Roy Stogner <roystgnr@...> wrote: > but "domain" is sort of PDE specific whereas LinearSolver is not. FWIW, "domain" is topological (more general than linear or even metric spaces). Jed 
From: Roy Stogner <roystgnr@ic...>  20100915 14:21:38

On Wed, 15 Sep 2010, Tim Kroeger wrote: > * (mainly to Roy:) Do you have some suggestion for a better name for > solve_only_on()? (Jed suggested set_active_domain().) What did I just suggest  restrict_to_subset()? I like "set_active_", and "set_active_domain" would have been fine for the System:: method, but "domain" is sort of PDE specific whereas LinearSolver is not. "set_active_indices", "set_active_dofs", "restrict_to_dofs"... I'm not especially attached to any one name; your pick.  Roy 
From: Roy Stogner <roystgnr@ic...>  20100915 14:16:32

On Wed, 15 Sep 2010, Tim Kroeger wrote: > Just to keep this clear: NumericVector::solve_only_on() will only accept a > list of dofs, nothing else, because in particular a NumericVector does not > know anything about subdomain ids. Absolutely. > On the other hand, System::solve_only_on() will only accept (a > reference to) an object of the (new) SystemSubset class. A derived > class SystemSubsetBySubdomain will exist, and users can write their > own derived classes if they need. That sounds great. > No, we don't have such a function. That is, at least in the way I thought of > this, solve_only_on() will make *all* successive calls to *any* of the > solve() methods restrict on the given subdomain. Hence, there must be a > possibility to reset this functionality. That's why I used a pointer. We > can, if you like, also use a reference and instead use a different method > solve_everywhere(void) to reset this, but I find this less intuitive. Ah, I see. I don't find it intuitive to change the internal state with a method named solve_only_on(); that sounded like it would have no side effects. Could we instead do: restrict_solves_to(); solve(); It's one extra line of code for users, but this way it's clear that the restriction call is making some state change and future solves must be affected.  Roy 
From: Tim Kroeger <tim.kroeger@ce...>  20100915 14:15:18

On Wed, 15 Sep 2010, Roy Stogner wrote: > On Wed, 15 Sep 2010, Tim Kroeger wrote: > >> On Tue, 14 Sep 2010, Roy Stogner wrote: >> >>> Are you looking at just all elements which have a shape function >>> associated with this DoF, or are you also including every element >>> which supports an associated basis function? In the former case it's >>> enough to find elements which point to this DoF's node; in the latter >>> case in an adapted mesh you've also got to check for other DoFs which >>> are constrained by this one. >> >> Good point, I didn't think of this. I'll try to sort things in my brain >> and see which of these is what I want. I see that in the latter case, I >> can use DofMap::constrain_nothing() to get the correct dof list, right? > > I don't think so. Suppose there's one hanging node where DoF 10 is > constrained in terms of DoFs 3 and 4, and another where DoF 15 is > constrained in terms of DoFs 5 and 3. DoFMap::constrain_nothing({10}) > will give you {10,3,4}, but what you'd want to be able to do is go from > {3} to {3,10,15}. I don't think there's anything in the library that > does that right now. Yes, but I meant when creating the cache. That is, I can then do it the other war round. If DoF 10 is contained in the given element of subdomain 42, I get {10,3,4}, so that in particular the minmax interval of DoF 3 gets enlarged to contain the value 42. Best Regards, Tim  Dr. Tim Kroeger CeVis  Center of Complex Systems and Visualization University of Bremen tim.kroeger@... Universitaetsallee 29 tim.kroeger@... D28359 Bremen Phone +494212187710 Germany Fax +494212184236 
From: Jed Brown <jed@59...>  20100915 14:15:14

On Wed, 15 Sep 2010 14:53:18 +0200 (CEST), Tim Kroeger <tim.kroeger@...> wrote: > > This chooses the sizes of each part, but you will still have to move > > indices around. > > This is done by VecScatter, isn't it? I think that the overhead of > doing this once before and after the solve will be small against the > gain in performance by having balanced vectors. Well, of course this > depends on a lot of things, but in the mean... VecScatter moves entries around, and isn't the issue here. Redistributing indices to satisfy some apriori size bound is somewhat different (though not hard). > Anyway, I have now attached the current state of my extension. > Questions are: > > * (mainly to Jed:) Are the things in petsc_linear_solver.{h,C} now > correct, sensible, and reasonably efficient? If yes, I will do the > same things to the remaining overloaded solve() methods. > Index: include/solvers/petsc_linear_solver.h > =================================================================== >  include/solvers/petsc_linear_solver.h (Revision 3957) > +++ include/solvers/petsc_linear_solver.h (Arbeitskopie) > @@ 122,7 +122,15 @@ > * Initialize data structures if not done so already plus much more > */ > void init (PetscMatrix<T>* matrix); > + > /** > + * After calling this method, all successive solves will be limited > + * to the given set of dofs. This mode can be disabled by calling > + * this method with \p dofs being a \p NULL pointer. > + */ > + virtual void solve_only_on (const std::vector<unsigned int>* const dofs); > + > + /** > * Call the Petsc solver. It calls the method below, using the > * same matrix for the system and preconditioner matrices. > */ > @@ 253,13 +261,35 @@ > * Krylov subspace context > */ > KSP _ksp; > + > + /** > + * PETSc index set containing the dofs on which to solve (\p NULL > + * means solve on all dofs). > + */ > + IS _solve_only_on_is; > + > + /** > + * PETSc index set containing the <em>local</em> dofs on which to > + * solve. This index set is determined when a system is solved on a > + * subset the first time. Before that point (or if no subset is > + * selected), it will be a \p NULL pointer. > + */ > + IS _solve_only_on_is_local; I don't think you want or need this, see below. > + > + /** > + * Internal method that returns the size of \p _solve_only_on_is. > + */ > + size_t _solve_only_on_is_size(void)const; Note that you only call this once. > + > }; > > > /* functions */ > template <typename T> > inline > PetscLinearSolver<T>::PetscLinearSolver () > + PetscLinearSolver<T>::PetscLinearSolver (): > + _solve_only_on_is(NULL), > + _solve_only_on_is_local(NULL) > { > if (libMesh::n_processors() == 1) > this>_preconditioner_type = ILU_PRECOND; > @@ 277,6 +307,23 @@ > } > > > + > +template <typename T> > +inline size_t > +PetscLinearSolver<T>:: > +_solve_only_on_is_size(void)const > +{ > + libmesh_assert(_solve_only_on_is!=NULL); > + > + PetscInt s; > + int ierr = ISGetSize(_solve_only_on_is,&s); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); > + > + return static_cast<size_t>(s); > +} > + > + > + > } // namespace libMesh > > > Index: src/solvers/petsc_linear_solver.C > =================================================================== >  src/solvers/petsc_linear_solver.C (Revision 3957) > +++ src/solvers/petsc_linear_solver.C (Arbeitskopie) > @@ 98,6 +98,8 @@ > { > if (this>initialized()) > { > + this>solve_only_on(NULL); > + > this>_is_initialized = false; > > int ierr=0; > @@ 294,7 +296,7 @@ > CHKERRABORT(libMesh::COMM_WORLD,ierr); > > // Set operators. The input matrix works as the preconditioning matrix >  ierr = KSPSetOperators(_ksp, matrix>mat(), matrix>mat(),SAME_NONZERO_PATTERN); > + ierr = KSPSetOperators(_ksp, matrix>mat(), matrix>mat(),DIFFERENT_NONZERO_PATTERN); > CHKERRABORT(libMesh::COMM_WORLD,ierr); > > // Set userspecified solver and preconditioner types > @@ 360,9 +362,46 @@ > > > > +template <typename T> > +void > +PetscLinearSolver<T>::solve_only_on (const std::vector<unsigned int>* const dofs) > +{ > + libmesh_assert(this>initialized()); > > + int ierr=0; > > + if(_solve_only_on_is!=NULL) > + { > + ierr = ISDestroy(_solve_only_on_is); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); > + _solve_only_on_is = NULL; > + } > > + if(_solve_only_on_is_local!=NULL) > + { > + ierr = ISDestroy(_solve_only_on_is_local); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); > + _solve_only_on_is_local = NULL; > + } > + > + if(dofs!=NULL) > + { > + PetscInt* petsc_dofs; > + ierr = PetscMalloc(dofs>size()*sizeof(PetscInt),&petsc_dofs); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); > + > + for(size_t i=0; i<dofs>size(); i++) > + { > + petsc_dofs[i] = (*dofs)[i]; > + } > + > + ierr = ISCreateGeneralNC(libMesh::COMM_WORLD,dofs>size(),petsc_dofs,&_solve_only_on_is); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); > + } > +} > + > + > + > template <typename T> > std::pair<unsigned int, Real> > PetscLinearSolver<T>::solve (SparseMatrix<T>& matrix_in, > @@ 406,10 +445,15 @@ > > // 2.1.x & earlier style > #if PETSC_VERSION_LESS_THAN(2,2,0) > + > + if(_solve_only_on_is!=NULL) > + { > + libmesh_not_implemented(); > + } > > // Set operators. The input matrix works as the preconditioning matrix > ierr = SLESSetOperators(_sles, matrix>mat(), precond>mat(), >  SAME_NONZERO_PATTERN); > + DIFFERENT_NONZERO_PATTERN); > CHKERRABORT(libMesh::COMM_WORLD,ierr); > > // Set the tolerances for the iterative solver. Use the usersupplied > @@ 431,6 +475,11 @@ > // 2.2.0 > #elif PETSC_VERSION_LESS_THAN(2,2,1) > > + if(_solve_only_on_is!=NULL) > + { > + libmesh_not_implemented(); > + } > + > // Set operators. The input matrix works as the preconditioning matrix > //ierr = KSPSetOperators(_ksp, matrix>mat(), precond>mat(), > // SAME_NONZERO_PATTERN); > @@ 475,19 +524,77 @@ > // 2.2.1 & newer style > #else > >  // Set operators. The input matrix works as the preconditioning matrix >  if(!this>same_preconditioner) >  { >  ierr = KSPSetOperators(_ksp, matrix>mat(), precond>mat(), >  SAME_NONZERO_PATTERN); >  CHKERRABORT(libMesh::COMM_WORLD,ierr); >  } > + Mat submat = NULL; > + Mat subprecond = NULL; > + Vec subrhs; > + Vec subsolution; > + VecScatter scatter; > + > + // Set operators. Also restrict rhs and solution vector to > + // subdomain if neccessary. > + if(_solve_only_on_is!=NULL) > + { > + size_t is_size = this>_solve_only_on_is_size(); > + > + ierr = VecCreate(libMesh::COMM_WORLD,&subrhs); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); > + ierr = VecSetSizes(subrhs,PETSC_DECIDE,is_size); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); > + ierr = VecSetFromOptions(subrhs); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); > + > + ierr = VecCreate(libMesh::COMM_WORLD,&subsolution); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); > + ierr = VecSetSizes(subsolution,PETSC_DECIDE,is_size); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); > + ierr = VecSetFromOptions(subsolution); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); > + > + ierr = VecScatterCreate(rhs>vec(),_solve_only_on_is, subrhs,NULL, &scatter); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); > + > + ierr = VecScatterBegin(scatter,rhs>vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); > + ierr = VecScatterEnd(scatter,rhs>vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); > + > + ierr = VecScatterBegin(scatter,solution>vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); > + ierr = VecScatterEnd(scatter,solution>vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); > + > + /* Make the local index set if it doesn't exist yet. */ > + if(_solve_only_on_is_local==NULL) > + { > + PetscInt firstIndex = 0; > + PetscInt lastIndex = 0; > + ierr = VecGetOwnershipRange(subrhs,&firstIndex,&lastIndex); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); > + > + ierr = ISCreateStride(libMesh::COMM_WORLD,lastIndexfirstIndex,firstIndex,1,&_solve_only_on_is_local); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); > + } > + > + ierr = MatGetSubMatrix(matrix>mat(), > + _solve_only_on_is_local,_solve_only_on_is, > + PETSC_DECIDE,MAT_INITIAL_MATRIX,&submat); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); Why do you want this nonsquare matrix? I'm pretty sure you want to extract the symmetric block (_solve_only_on_is,_solve_only_on_is). > + > + ierr = MatGetSubMatrix(precond>mat(), > + _solve_only_on_is_local,_solve_only_on_is, > + PETSC_DECIDE,MAT_INITIAL_MATRIX,&subprecond); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); > + > + ierr = KSPSetOperators(_ksp, submat, subprecond, > + this>same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); > + } > else >  { >  ierr = KSPSetOperators(_ksp, matrix>mat(), precond>mat(), >  SAME_PRECONDITIONER); >  CHKERRABORT(libMesh::COMM_WORLD,ierr); >  } > + { > + ierr = KSPSetOperators(_ksp, matrix>mat(), precond>mat(), > + this>same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); > + } > > // Set the tolerances for the iterative solver. Use the usersupplied > // tolerance for the relative residual & leave the others at default values. > @@ 496,8 +603,16 @@ > CHKERRABORT(libMesh::COMM_WORLD,ierr); > > // Solve the linear system >  ierr = KSPSolve (_ksp, rhs>vec(), solution>vec()); >  CHKERRABORT(libMesh::COMM_WORLD,ierr); > + if(_solve_only_on_is!=NULL) > + { > + ierr = KSPSolve (_ksp, subrhs, subsolution); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); > + } > + else > + { > + ierr = KSPSolve (_ksp, rhs>vec(), solution>vec()); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); > + } > > // Get the number of iterations required for convergence > ierr = KSPGetIterationNumber (_ksp, &its); > @@ 507,6 +622,26 @@ > ierr = KSPGetResidualNorm (_ksp, &final_resid); > CHKERRABORT(libMesh::COMM_WORLD,ierr); > > + if(_solve_only_on_is!=NULL) > + { > + ierr = VecScatterBegin(scatter,subsolution,solution>vec(),INSERT_VALUES,SCATTER_REVERSE); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); > + ierr = VecScatterEnd(scatter,subsolution,solution>vec(),INSERT_VALUES,SCATTER_REVERSE); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); You probably want to set the rest of solution>vec() to the right side, or zero it out, otherwise a solve is not a linear operation. The reason for the first is that it makes the preconditioner look like the submatrix was extended to operate on the full domain using the identity (instead of zero which would make it very singular). > + > + ierr = VecScatterDestroy(scatter); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); > + > + ierr = VecDestroy(subsolution); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); > + ierr = VecDestroy(subrhs); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); > + ierr = MatDestroy(submat); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); > + ierr = MatDestroy(subprecond); > + CHKERRABORT(libMesh::COMM_WORLD,ierr); > + } > + > #endif > > STOP_LOG("solve()", "PetscLinearSolver"); > @@ 574,7 +709,7 @@ > > // Set operators. The input matrix works as the preconditioning matrix > ierr = KSPSetOperators(_ksp, mat, mat, >  SAME_NONZERO_PATTERN); > + DIFFERENT_NONZERO_PATTERN); > CHKERRABORT(libMesh::COMM_WORLD,ierr); > > // Set the tolerances for the iterative solver. Use the usersupplied > Index: src/solvers/linear_solver.C > =================================================================== >  src/solvers/linear_solver.C (Revision 3957) > +++ src/solvers/linear_solver.C (Arbeitskopie) > @@ 115,7 +115,18 @@ > _preconditioner = preconditioner; > } > > +template <typename T> > +void > +LinearSolver<T>::solve_only_on(const std::vector<unsigned int>* const dofs) > +{ > + if(dofs!=NULL) > + { > + libmesh_not_implemented(); > + } > +} > > + > + > // > // Explicit instantiations > template class LinearSolver<Number>; 
From: Roy Stogner <roystgnr@ic...>  20100915 14:10:26

On Wed, 15 Sep 2010, Tim Kroeger wrote: > On Tue, 14 Sep 2010, Roy Stogner wrote: > >> Are you looking at just all elements which have a shape function >> associated with this DoF, or are you also including every element >> which supports an associated basis function? In the former case it's >> enough to find elements which point to this DoF's node; in the latter >> case in an adapted mesh you've also got to check for other DoFs which >> are constrained by this one. > > Good point, I didn't think of this. I'll try to sort things in my brain and > see which of these is what I want. I see that in the latter case, I can use > DofMap::constrain_nothing() to get the correct dof list, right? I don't think so. Suppose there's one hanging node where DoF 10 is constrained in terms of DoFs 3 and 4, and another where DoF 15 is constrained in terms of DoFs 5 and 3. DoFMap::constrain_nothing({10}) will give you {10,3,4}, but what you'd want to be able to do is go from {3} to {3,10,15}. I don't think there's anything in the library that does that right now.  Roy 
From: Tim Kroeger <tim.kroeger@ce...>  20100915 12:53:39

On Wed, 15 Sep 2010, Jed Brown wrote: > On Wed, 15 Sep 2010 08:42:16 +0200 (CEST), Tim Kroeger <tim.kroeger@...> wrote: > >> (One conjecture why SAME_NONZERO_PATTERN has been used here before is >> that somebody thought it means that the system matrix and the >> preconditioner matrix have the same nonzero pattern, but as far as I >> understand now, this is just not true.) > > It means that the nonzero pattern of the preconditioning matrix has not > changed since the last solve. Yes, that's also what I understood now. Since at no point in the code information about the nonzero pattern of the last solve is stored, there is no reason for using SAME_NONZERO_PATTERN at all. >> Another idea is to use VecCreateMPI(), giving only the global length >> and using PETSC_DECIDE for the local length. This way, approximately >> equal amounts of dofs are assigned to each processor, but will they >> correspond to patches of the computational domain, or will they be >> scattered all around? Well, since we are using a vector, this is of >> course in the responsiblity of the user. > > This chooses the sizes of each part, but you will still have to move > indices around. This is done by VecScatter, isn't it? I think that the overhead of doing this once before and after the solve will be small against the gain in performance by having balanced vectors. Well, of course this depends on a lot of things, but in the mean... > For a highquality partition, you could extract the > submatrix with the naive (current) partition, then run a partitioner on > the submatrix, and finally extract indices using the resulting "good" > partition. Of course this may not be worthwhile if you will only solve > with it once. I leave this for somebody else. (: Anyway, I have now attached the current state of my extension. Questions are: * (mainly to Jed:) Are the things in petsc_linear_solver.{h,C} now correct, sensible, and reasonably efficient? If yes, I will do the same things to the remaining overloaded solve() methods. * (mainly ro Roy:) Do you like the API? * (mainly to Roy:) Do you have some suggestion for a better name for solve_only_on()? (Jed suggested set_active_domain().) Best Regards, Tim  Dr. Tim Kroeger CeVis  Center of Complex Systems and Visualization University of Bremen tim.kroeger@... Universitaetsallee 29 tim.kroeger@... D28359 Bremen Phone +494212187710 Germany Fax +494212184236 
From: Jed Brown <jed@59...>  20100915 10:52:12

On Wed, 15 Sep 2010 08:42:16 +0200 (CEST), Tim Kroeger <tim.kroeger@...> wrote: > (One conjecture why SAME_NONZERO_PATTERN has been used here before is > that somebody thought it means that the system matrix and the > preconditioner matrix have the same nonzero pattern, but as far as I > understand now, this is just not true.) It means that the nonzero pattern of the preconditioning matrix has not changed since the last solve. The purpose of this is to avoid needing to recompute a symbolic factorization or some setup steps for AMG. Usually the symbolic factorization is much cheaper than numeric factorization so it makes little difference, but this eventually ceases to be true as the number of processes increase. > Another idea is to use VecCreateMPI(), giving only the global length > and using PETSC_DECIDE for the local length. This way, approximately > equal amounts of dofs are assigned to each processor, but will they > correspond to patches of the computational domain, or will they be > scattered all around? Well, since we are using a vector, this is of > course in the responsiblity of the user. This chooses the sizes of each part, but you will still have to move indices around. For a highquality partition, you could extract the submatrix with the naive (current) partition, then run a partitioner on the submatrix, and finally extract indices using the resulting "good" partition. Of course this may not be worthwhile if you will only solve with it once. Jed 
From: Tim Kroeger <tim.kroeger@ce...>  20100915 10:44:18

On Wed, 15 Sep 2010, Jed Brown wrote: > On Tue, 14 Sep 2010 09:00:39 +0200 (CEST), Tim Kroeger <tim.kroeger@...> wrote: >> Well, I could make it a vector instead, but then, on the other hand, I >> guess that in most situations the user does not have any idea of which >> order is best to use. > > That's what ordering packages are for, but defining the interface to use > std::set prohibits a smart user from providing an ordering. Yes, you're right. >> Well, what I would like to do, but don't know how, is to let subrhs be >> a vector with the same layout as rhs and with any dof being owned by >> that processor which owns the corresponding dof in rhs. Or, perhaps >> better, do an automatic load balancing for subrhs. (Likewise for >> subsolution of course.) > > Each process provides a list of global indices that it wants to own in > the subproblem. The length of this list determines the local size. > You probably want > > VecCreate(comm,X); > VecSetSizes(X,local_size,PETSC_DETERMINE); > VecSetFromOptions(X); > > which will create a serial vector in serial and a parallel vector in > parallel. The user can use a partitioning algorithm if they want to > redistribute dofs (or you could add a runtime option to use an arbitrary > algorithm, PETSc offers MatPartitioning for this). Ok, I understand. However, using VecSetSizes(X,PETSC_DECIDE,global_size) seems even more useful to me. Or am I wrong somehow? >> Also, I am unsure whether or not I have to take care about the matrix >> rows being owned by the same processors that own the corresponding >> vector dofs. > > It's all determined by the layout of the IS. If you don't want rows to > "move", then the IS should only consist of global indices that are > presently owned. Ah, I understand. That means that the current implementation of my new method will make a matrix where each entry is stored on all processors. This is certainly not efficient. What does PETSc do if the matrix rows are owned by the wrong processors? Does is crash, or will it just compensate by communicating all the information around as necessary? Anyway, if I do VecSetSizes(X,PETSC_DECIDE,global_size) as described above, the best thing would then perhaps to be to ask the vector about which dofs are local to the current processor and then tell the matrix to own the rows that correspond to these dofs (but still all columns since otherwise some components are lost). I guess that VecGetOwnershipRange() is the thing to do then, right? However, if the dofs are *not* consecutive, what will VecGetOwnershipRange() do then? >>> If dofs == NULL, this function is a noop? >> >> Yes. You see, this function is for all solver packages except PETSc, >> and they should just give an error message when somebody attempts to >> use solve_only_on(). But, of course, if a NULL pointer is used, which >> means solve on the entire domain, there is no need to give an error >> message, because this is what is done anyway, isn't it? (: > > I wrote that comment first, before realizing that solve_only_on was just > restricting the domain and didn't solve anything. Perhaps something > like set_active_domain would be a better name? Yes, Roy also seemed to be confused by this. I'll wait whether he suggests some name, and if not, I'll use you suggestion. Best Regards, Tim  Dr. Tim Kroeger CeVis  Center of Complex Systems and Visualization University of Bremen tim.kroeger@... Universitaetsallee 29 tim.kroeger@... D28359 Bremen Phone +494212187710 Germany Fax +494212184236 
From: Jed Brown <jed@59...>  20100915 10:13:24

On Tue, 14 Sep 2010 09:00:39 +0200 (CEST), Tim Kroeger <tim.kroeger@...> wrote: > Well, I could make it a vector instead, but then, on the other hand, I > guess that in most situations the user does not have any idea of which > order is best to use. That's what ordering packages are for, but defining the interface to use std::set prohibits a smart user from providing an ordering. > > So the LinearSolver doesn't take ownership of this set? It seems > > horrible to make the user responsible for that. > > Well, users are supposed to do things like > > std::set<unsigned int> a; > // Add values here > linear_solver.solve_only_on(&a); I would expect the user to have a function that selects the subdomain and registers it. When the stack frame goes away, it's destructor would be called, which is why I would normally want the solver object to reference it. But whatever. > Well, what I would like to do, but don't know how, is to let subrhs be > a vector with the same layout as rhs and with any dof being owned by > that processor which owns the corresponding dof in rhs. Or, perhaps > better, do an automatic load balancing for subrhs. (Likewise for > subsolution of course.) Each process provides a list of global indices that it wants to own in the subproblem. The length of this list determines the local size. You probably want VecCreate(comm,X); VecSetSizes(X,local_size,PETSC_DETERMINE); VecSetFromOptions(X); which will create a serial vector in serial and a parallel vector in parallel. The user can use a partitioning algorithm if they want to redistribute dofs (or you could add a runtime option to use an arbitrary algorithm, PETSc offers MatPartitioning for this). > Also, I am unsure whether or not I have to take care about the matrix > rows being owned by the same processors that own the corresponding > vector dofs. It's all determined by the layout of the IS. If you don't want rows to "move", then the IS should only consist of global indices that are presently owned. > > If dofs == NULL, this function is a noop? > > Yes. You see, this function is for all solver packages except PETSc, > and they should just give an error message when somebody attempts to > use solve_only_on(). But, of course, if a NULL pointer is used, which > means solve on the entire domain, there is no need to give an error > message, because this is what is done anyway, isn't it? (: I wrote that comment first, before realizing that solve_only_on was just restricting the domain and didn't solve anything. Perhaps something like set_active_domain would be a better name? Jed 
From: Tim Kroeger <tim.kroeger@ce...>  20100915 06:50:33

On Wed, 15 Sep 2010, Tim Kroeger wrote: > Just to keep this clear: NumericVector::solve_only_on() will only > accept a list of dofs, nothing else, because in particular a > NumericVector does not know anything about subdomain ids. On the > other hand, System::solve_only_on() will only accept (a reference to) > an object of the (new) SystemSubset class. A derived class > SystemSubsetBySubdomain will exist, and users can write their own > derived classes if they need. /NumericVector/LinearSolver/  Dr. Tim Kroeger CeVis  Center of Complex Systems and Visualization University of Bremen tim.kroeger@... Universitaetsallee 29 tim.kroeger@... D28359 Bremen Phone +494212187710 Germany Fax +494212184236 
From: Tim Kroeger <tim.kroeger@ce...>  20100915 06:42:30

On Tue, 14 Sep 2010, Roy Stogner wrote: > On Tue, 14 Sep 2010, Tim Kroeger wrote: > >> I suggest to make two methods, one of which takes a vector and the other >> takes a set. Then, the settaking method will by default just convert the >> set to a vector and call the vectortaking method. Roy, what do you think? > > I think this is one of those cases like our send_list accumulation > where the most efficient way to build the (conceptual) set is to build > a vector of nonunique indices and then sortuniq it after we're done. > If the user wants to play with the indexing order for efficiency, then > having a way to convey the order is a bonus, but in any case there'd > be no reason to pass a set at all. Even my initial idea of using one > container to signify dof ids and a different container to signify > subdomain ids was a lousy one; if we're not going to make some modular > Subdomain class, then it would be much less confusing to just use two > different method names, e.g. solve_on_dofs(vector&) and > solve_on_subdomains(vector&). Okay, so I will switch this to std::vector. Just to keep this clear: NumericVector::solve_only_on() will only accept a list of dofs, nothing else, because in particular a NumericVector does not know anything about subdomain ids. On the other hand, System::solve_only_on() will only accept (a reference to) an object of the (new) SystemSubset class. A derived class SystemSubsetBySubdomain will exist, and users can write their own derived classes if they need. >> Yes, that's a good point. Also, if the same solver object is used to solve >> different systems, there seems to be no reason to use SAME_NONZERO_PATTERN >> here. An further more, if some grid refinement/coarsening has taken place >> in the meantime, the nonzero pattern is definitely not the same as before. >> I really don't know why SAME_NONZERO_PATTERN has been used here originally. >> Roy, can you comment on this? > > Hmm... according to svn blame, David is using SAME_NONZERO_PATTERN > there because it's what I did nearby, and I'm using it because it's > what Ben did first... > > And if I had to guess I'd say that Ben did it because the way he'd > construct a separate preconditioning matrix is to make a copy > of the system matrix with the same sparsity pattern but different > data. We'd have to add a bit of cached information if we wanted to > quickly test whether two arbitrary PetscMatrix objects corresponded to > the same sparsity pattern, so the best fix is probably just to switch > those calls to DIFFERENT_NONZERO_PATTERN for now. Okay, I'll do this along the way. (One conjecture why SAME_NONZERO_PATTERN has been used here before is that somebody thought it means that the system matrix and the preconditioner matrix have the same nonzero pattern, but as far as I understand now, this is just not true.) >> I know that libMesh uses some external libraries for load balancing in >> general, but I have never dealt directly with that part of code, so I don't >> know whether and, if yes, how to involve it here. > > Even though you'd need to get a new partition the subdomain for > optimal scaling, in most cases you'd probably do better using the > existing partitioning and accepting a significant imbalance than > running in serial and being completely imbalanced. Yes, but I still don't know how to do this. Another idea is to use VecCreateMPI(), giving only the global length and using PETSC_DECIDE for the local length. This way, approximately equal amounts of dofs are assigned to each processor, but will they correspond to patches of the computational domain, or will they be scattered all around? Well, since we are using a vector, this is of course in the responsiblity of the user. Also, I still don't know whether I will have to do something with the matrix to make sure that the corresponding rows are owned by the corresponding processores. >> But, of course, if a NULL pointer is used, which means solve on the entire >> domain > > Do we need that semantic? I'd just use a reference instead of a > pointer; we've already got a function for solving on the entire > domain. No, we don't have such a function. That is, at least in the way I thought of this, solve_only_on() will make *all* successive calls to *any* of the solve() methods restrict on the given subdomain. Hence, there must be a possibility to reset this functionality. That's why I used a pointer. We can, if you like, also use a reference and instead use a different method solve_everywhere(void) to reset this, but I find this less intuitive. Best Regards, Tim  Dr. Tim Kroeger CeVis  Center of Complex Systems and Visualization University of Bremen tim.kroeger@... Universitaetsallee 29 tim.kroeger@... D28359 Bremen Phone +494212187710 Germany Fax +494212184236 
From: Tim Kroeger <tim.kroeger@ce...>  20100915 06:21:47

On Tue, 14 Sep 2010, Roy Stogner wrote: > On Tue, 14 Sep 2010, Tim Kroeger wrote: > >> For a dof of a vector the corresponds to a System, I need to know the >> maximal subdomain_id of all elements which contain this dof. > >> While it is principally possible to determine this on the fly (via >> Elem::find_point_neighbor() and then a loop over the resulting >> elements and a loop over their respective dofs), this seems to be >> very slow (I guess too slow). > > It also might be incorrect, depending on your definition of "contain". > Are you looking at just all elements which have a shape function > associated with this DoF, or are you also including every element > which supports an associated basis function? In the former case it's > enough to find elements which point to this DoF's node; in the latter > case in an adapted mesh you've also got to check for other DoFs which > are constrained by this one. Good point, I didn't think of this. I'll try to sort things in my brain and see which of these is what I want. I see that in the latter case, I can use DofMap::constrain_nothing() to get the correct dof list, right? > But when you want to cache it: > >> std::map<unsigned int, unsigned char>: Seems to be the best solution >> that I can currently think of. Can also be extended to >> std::map<unsigned int,std::pair<unsigned char,unsigned char> > so that >> minimal and maximal subdomain_id are stored in one structure. >> >> Any suggestions? > > If efficiency is important, put the local dofs in a vector and just > the ghost dofs in a map. > > If efficiency is really important, try both map and hash_map. Grep for > HASH_MAP and you'll find two places where we're using them for > multimap containers (in one of 7 incarnations depending on how up to > date your compiler is), and you can use the same #if defined() tests > to find a supported hash_map type or fallback. > > If you want to be awesome, add this as a GHOSTED option to the > DistributedVector class. PETSc may not want to play with multiple > scalar types but it'd be no problem for us to add additional > instantiations like DistributedVector<unsigned int>. Thank you for your suggestions. I will try with std::map first and see whether it takes some significant computation time in my application. If it does, I'll extend the DistributedVector class appropriately using one of the more efficient suggestions that you made. Best Regards, Tim  Dr. Tim Kroeger CeVis  Center of Complex Systems and Visualization University of Bremen tim.kroeger@... Universitaetsallee 29 tim.kroeger@... D28359 Bremen Phone +494212187710 Germany Fax +494212184236 
From: Roy Stogner <roystgnr@ic...>  20100915 04:01:08

On Tue, 14 Sep 2010, Tim Kroeger wrote: > For a dof of a vector the corresponds to a System, I need to know the > maximal subdomain_id of all elements which contain this dof. > While it is principally possible to determine this on the fly (via > Elem::find_point_neighbor() and then a loop over the resulting > elements and a loop over their respective dofs), this seems to be > very slow (I guess too slow). It also might be incorrect, depending on your definition of "contain". Are you looking at just all elements which have a shape function associated with this DoF, or are you also including every element which supports an associated basis function? In the former case it's enough to find elements which point to this DoF's node; in the latter case in an adapted mesh you've also got to check for other DoFs which are constrained by this one. But when you want to cache it: > std::map<unsigned int, unsigned char>: Seems to be the best solution > that I can currently think of. Can also be extended to > std::map<unsigned int,std::pair<unsigned char,unsigned char> > so that > minimal and maximal subdomain_id are stored in one structure. > > Any suggestions? If efficiency is important, put the local dofs in a vector and just the ghost dofs in a map. If efficiency is really important, try both map and hash_map. Grep for HASH_MAP and you'll find two places where we're using them for multimap containers (in one of 7 incarnations depending on how up to date your compiler is), and you can use the same #if defined() tests to find a supported hash_map type or fallback. If you want to be awesome, add this as a GHOSTED option to the DistributedVector class. PETSc may not want to play with multiple scalar types but it'd be no problem for us to add additional instantiations like DistributedVector<unsigned int>.  Roy 
From: Roy Stogner <roystgnr@ic...>  20100915 03:41:46

On Tue, 14 Sep 2010, Tim Kroeger wrote: > I suggest to make two methods, one of which takes a vector and the other > takes a set. Then, the settaking method will by default just convert the > set to a vector and call the vectortaking method. Roy, what do you think? I think this is one of those cases like our send_list accumulation where the most efficient way to build the (conceptual) set is to build a vector of nonunique indices and then sortuniq it after we're done. If the user wants to play with the indexing order for efficiency, then having a way to convey the order is a bonus, but in any case there'd be no reason to pass a set at all. Even my initial idea of using one container to signify dof ids and a different container to signify subdomain ids was a lousy one; if we're not going to make some modular Subdomain class, then it would be much less confusing to just use two different method names, e.g. solve_on_dofs(vector&) and solve_on_subdomains(vector&). > Yes, that's a good point. Also, if the same solver object is used to solve > different systems, there seems to be no reason to use SAME_NONZERO_PATTERN > here. An further more, if some grid refinement/coarsening has taken place in > the meantime, the nonzero pattern is definitely not the same as before. I > really don't know why SAME_NONZERO_PATTERN has been used here originally. > Roy, can you comment on this? Hmm... according to svn blame, David is using SAME_NONZERO_PATTERN there because it's what I did nearby, and I'm using it because it's what Ben did first... And if I had to guess I'd say that Ben did it because the way he'd construct a separate preconditioning matrix is to make a copy of the system matrix with the same sparsity pattern but different data. We'd have to add a bit of cached information if we wanted to quickly test whether two arbitrary PetscMatrix objects corresponded to the same sparsity pattern, so the best fix is probably just to switch those calls to DIFFERENT_NONZERO_PATTERN for now. > I know that libMesh uses some external libraries for load balancing in > general, but I have never dealt directly with that part of code, so I don't > know whether and, if yes, how to involve it here. Even though you'd need to get a new partition the subdomain for optimal scaling, in most cases you'd probably do better using the existing partitioning and accepting a significant imbalance than running in serial and being completely imbalanced. > But, of course, if a NULL pointer is used, which means > solve on the entire domain Do we need that semantic? I'd just use a reference instead of a pointer; we've already got a function for solving on the entire domain.  Roy 