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}
(73) 
_{Jun}
(29) 
_{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: Kirk, Benjamin (JSCEG311) <benjamin.kirk1@na...>  20100920 15:00:24

On 9/20/10 9:48 AM, "yunfei zhu" <yzhu2009@...> wrote: > and the sum of > all the weights are far less than 1. The weights should sum to the volume of the element in your reference coordinate system, which for us I believe is 1/6. 
From: André Luis Rossa <andre_rossa@ho...>  20100920 14:52:19

Ok Roy! Thanks for your help. Andre > Date: Sun, 19 Sep 2010 10:52:50 0500 > From: roystgnr@... > To: andre_rossa@... > CC: libmeshusers@... > Subject: Re: [Libmeshusers] Scaling the variables for error estimation > > > On Sat, 18 Sep 2010, André Luis Rossa wrote: > > > The version 0.6.4 does not support "component_scale" to scale (and > > to choose) the variables to calculate the error norm wich is used by > > the error estimator for AMR/C as the 0.6.3 version does. > > > How to reproduce the scaling error through the "SystemNorm > > error_norm" (error_estimator.h) for the kelly_error_estimator? > > For each variable i, set its scaling by > my_estimator.error_norm.set_weight(i, weight_) > > If you leave any weights unset, they should default to 1.0. >  > Roy 
From: yunfei zhu <yzhu2009@go...>  20100920 14:48:18

Hi I am trying to recover the stress base on the Tet4/10 elements. The quadrature rule for element tet4 seems quite confusing, could anyone give me an explanation? For example, If the quadrature order is SECOND, the quadrature points have only 3 Tetrahedral coordinates, and the sum of all the weights are far less than 1. I noted one of the comments in quadrature_gauss_3D.C, the tetrahedral quadrature rules are taken from pg. 222 of "The finite element method" vol.1 ed.5 by Zienkiewicz & Taylor, so I checked the details in the book. I find that for each of the quadrature point, there are 4 tetrahedral coordinates and the sum of the weights are equal to 1. Thanks in advance for help. Best wishes, yunfei 
From: Tim Kroeger <tim.kroeger@ce...>  20100920 06:39:47

On Fri, 17 Sep 2010, Roy Stogner wrote: > On Fri, 17 Sep 2010, Tim Kroeger wrote: > >> Thank you for this hint. Using this, it was easy to find the point where >> the macro is not called. It's in libMesh, however, that is in >> petsc_preconditioner.C, line 49. Actually, in that file, a lot of calls to >> CHKERRABORT() are missing. Roy, is there any reason for this? > > Derek's lazy? ;) > >> If not, I would like to add all these calls. > > Thanks! Did this now, but I wonder whether CHKERRABORT() is really the best thing to do. I guess that, at least in methods that are handed as callback functions to PETSc, CHKERRQ() would be the better choice, but that would also require passing the error code from PetscPreconditioner::apply() back to __libmesh_petsc_preconditioner_apply(). Well, essentially that isn't a problem either, other than that I already have a modified version of the file where __libmesh_petsc_preconditioner_apply() is implemented in, due to my subset stuff, and I don't want to set up a second checkedout copy of libMesh. If anybody agrees that doing this would be a good thing, I might do that later, when the subset stuff is finished. Provided that I don't forget. 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...>  20100919 15:52:58

On Sat, 18 Sep 2010, André Luis Rossa wrote: > The version 0.6.4 does not support "component_scale" to scale (and > to choose) the variables to calculate the error norm wich is used by > the error estimator for AMR/C as the 0.6.3 version does. > How to reproduce the scaling error through the "SystemNorm > error_norm" (error_estimator.h) for the kelly_error_estimator? For each variable i, set its scaling by my_estimator.error_norm.set_weight(i, weight_) If you leave any weights unset, they should default to 1.0.  Roy 
From: André Luis Rossa <andre_rossa@ho...>  20100918 00:00:48

Hi! The version 0.6.4 does not support "component_scale" to scale (and to choose) the variables to calculate the error norm wich is used by the error estimator for AMR/C as the 0.6.3 version does. How to reproduce the scaling error through the "SystemNorm error_norm" (error_estimator.h) for the kelly_error_estimator? Thanks. Andre 
From: Derek Gaston <friedmud@gm...>  20100917 15:17:08

On Fri, Sep 17, 2010 at 9:00 AM, Roy Stogner <roystgnr@...>wrote: > Derek's lazy? ;) > LOL  That usually is the case... but in this case I think it was just lack of diligence. I don't remember thinking "meh I'll just leave those out"... I actually just didn't think of them at all ;) Me being lazy is always a good guess though! I always maintain that people get into Computer Science because they are inherently lazy... why do for yourself what the computer can do for you! Derek 
From: Roy Stogner <roystgnr@ic...>  20100917 15:00:41

On Fri, 17 Sep 2010, Tim Kroeger wrote: > Thank you for this hint. Using this, it was easy to find the point where the > macro is not called. It's in libMesh, however, that is in > petsc_preconditioner.C, line 49. Actually, in that file, a lot of calls to > CHKERRABORT() are missing. Roy, is there any reason for this? Derek's lazy? ;) > If not, I would like to add all these calls. Thanks!  Roy 
From: Tim Kroeger <tim.kroeger@ce...>  20100917 14:58:14

On Fri, 17 Sep 2010, Jed Brown wrote: > On Fri, Sep 17, 2010 at 16:35, Tim Kroeger > <tim.kroeger@...> wrote: >> Can you (Jed) comment on how it can happen that KSPSolve() prints this error >> message but nicely returns 0? > > CHKERRQ is not being called somewhere in the stack, if it's in PETSc > code, please let me know where it is (or send a full trace). > > As a general rule, always set a breakpoint in PetscError. I have in > my ~/.gdbinit: > > set breakpoint pending on > break PetscError > > Note that start_in_debugger adds this automatically. Thank you for this hint. Using this, it was easy to find the point where the macro is not called. It's in libMesh, however, that is in petsc_preconditioner.C, line 49. Actually, in that file, a lot of calls to CHKERRABORT() are missing. Roy, is there any reason for this? If not, I would like to add all these calls. 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...>  20100917 14:48:59

On Fri, Sep 17, 2010 at 16:35, Tim Kroeger <tim.kroeger@...> wrote: > Can you (Jed) comment on how it can happen that KSPSolve() prints this error > message but nicely returns 0? CHKERRQ is not being called somewhere in the stack, if it's in PETSc code, please let me know where it is (or send a full trace). As a general rule, always set a breakpoint in PetscError. I have in my ~/.gdbinit: set breakpoint pending on break PetscError Note that start_in_debugger adds this automatically. Jed 
From: Tim Kroeger <tim.kroeger@ce...>  20100917 14:36:01

I'm currently testing my application with the new subset stuff, and I'm observing some PETSc behaviour that I find strange. At some place in petsc_linear_solver.C, I do: ierr = KSPSolve (_ksp, subrhs, subsolution); CHKERRABORT(libMesh::COMM_WORLD,ierr); What now happens is that the first line produces some this message: [0]PETSC ERROR:  Error Message  [0]PETSC ERROR: Nonconforming object sizes! [0]PETSC ERROR: Mat mat,Vec x: global dim 3022 2858! [0]PETSC ERROR:  [0]PETSC ERROR: Petsc Release Version 3.1.0, Patch 4, Fri Jul 30 14:42:02 CDT 2010 [0]PETSC ERROR: See docs/changes/index.html for recent updates. [0]PETSC ERROR: See docs/faq.html for hints about trouble shooting. [0]PETSC ERROR: See docs/index.html for manual pages. [0]PETSC ERROR:  [0]PETSC ERROR: 1 on a linuxgnu named tkroegerpc by tim Fri Sep 17 16:02:10 2010 [0]PETSC ERROR: Libraries linked from /home/tim/archives/packages/petsc3.1p4/linuxgnu/lib [0]PETSC ERROR: Configure run at Thu Sep 16 10:33:53 2010 [0]PETSC ERROR: Configure options withcc=mpicc withcxx=mpicxx withmpilib=/home/tim/archives/packages/mpich21.1.1p1/lib/libmpich.so withmpiinclude=/home/tim/archives/packages/mpich21.1.1p1/src/include/ withshared=1 downloadsuperlu_dist downloadparmetis withpic [0]PETSC ERROR:  [0]PETSC ERROR: MatSolve() line 2807 in src/mat/interface/matrix.c [0]PETSC ERROR: PCApply_LU() line 189 in src/ksp/pc/impls/factor/lu/lu.c [0]PETSC ERROR: PCApply() line 357 in src/ksp/pc/interface/precon.c Well, this is okay since I might have done some mistake that makes vector and matrix not match each other. The point is, however, that the value of ierr after this is still 0, so that the CHKERRABORT() macro does not do anything and the program just continues. This made it very hard for me to locate the point in the code where this error is generated. I used a debugger to trace through the program and find this point. Can you (Jed) comment on how it can happen that KSPSolve() prints this error message but nicely returns 0? Thank you in advance, 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...>  20100916 08:09:20

On Thu, 16 Sep 2010, Jed Brown wrote: > Correct about shell and MFFD matrices. If the Mat does not implement MatGetSubMatrix then this call returns a matrix > of type MatSubMatrix which works like you say (>=3.1). Great, thank you. I will update to 3.1 now (was using 3.0.0p8 before). 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...>  20100916 07:56:59

Correct about shell and MFFD matrices. If the Mat does not implement MatGetSubMatrix then this call returns a matrix of type MatSubMatrix which works like you say (>=3.1). Preconditioning that thing obviously does not come automatically. Jed On Sep 16, 2010 8:41 AM, "Tim Kroeger" <tim.kroeger@...> wrote: On Wed, 15 Sep 2010, Jed Brown wrote: > > On Wed, 15 Sep 2010 17:34:14 +0200 (CEST), Tim Kroeger < > tim.kroeger@...> wrote: > > > > >> Assume further that the matrix is tridiagonal. If I then do what you > >> suggested, four matrix e... > Ah, great! Thank you very much! Another question comes to my mind in this context: Does the whole thing also work for a shell matrix? It seems to me that for a (blackbox) shell matrix, the only method to extract some submatrix is to provide a wrapper shell matrix function which extends the input vector by zeros, calls the original shell matrix, and then removes the unwanted components of the result. The question is whether PETSc does this whenever MatGetSubMatrix() is called on a shell matrix or whether I will have to do this myself. >> I don't find this weird. Well, we could let the user decide between >> several modes here, tha... I understand what you mean. First of all, we should let the user know what happens with the "inactive" dofs, otherwise the user might assume something which is not true. Then, also, it might be sensible to give the user a choice here since, depending on the application, one or another thing might be useful. Among the possible things to do, the action of just letting the inactive dof values unchanged is, although removing some wellknown properties from the whole thing, the fastest thing to do (in particular if the number of active dofs is very small) and also the action which I need in my application. But perhaps you are right in that it should not be the default action because it is confusing. I'll add an ENUM for this and let the default be zeroing the inactive dofs. I'll also add a comment about "my" mode leading to loss of wellknown properties. Best Regards, Tim  Dr. Tim Kroeger CeVis  Center of Complex Systems and Visualization Unive... 
From: Tim Kroeger <tim.kroeger@ce...>  20100916 07:56:25

Here is the current state now. I have no more questions (other than the one which I asked in the previous mail). If Jed says that everything looks well now, I will do the thing for the missing overloaded PetscLinearSolver::solve() methods.  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...>  20100916 06:41:36

On Wed, 15 Sep 2010, Jed Brown wrote: > On Wed, 15 Sep 2010 17:34:14 +0200 (CEST), Tim Kroeger <tim.kroeger@...> wrote: > >> 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. Ah, great! Thank you very much! Another question comes to my mind in this context: Does the whole thing also work for a shell matrix? It seems to me that for a (blackbox) shell matrix, the only method to extract some submatrix is to provide a wrapper shell matrix function which extends the input vector by zeros, calls the original shell matrix, and then removes the unwanted components of the result. The question is whether PETSc does this whenever MatGetSubMatrix() is called on a shell matrix or whether I will have to do this myself. >> 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. I understand what you mean. First of all, we should let the user know what happens with the "inactive" dofs, otherwise the user might assume something which is not true. Then, also, it might be sensible to give the user a choice here since, depending on the application, one or another thing might be useful. Among the possible things to do, the action of just letting the inactive dof values unchanged is, although removing some wellknown properties from the whole thing, the fastest thing to do (in particular if the number of active dofs is very small) and also the action which I need in my application. But perhaps you are right in that it should not be the default action because it is confusing. I'll add an ENUM for this and let the default be zeroing the inactive dofs. I'll also add a comment about "my" mode leading to loss of wellknown properties. 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: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 