Dear Jed, Thank you for your help so far. I still have some questions:
On Mon, 13 Sep 2010, Jed Brown wrote: > On Mon, 13 Sep 2010 15:18:30 +0200 (CEST), Tim Kroeger > <[email protected]> wrote: > >> + >> + /** >> + * Set of dofs to which any solve is limited (\p NULL means solve on >> + * all dofs). >> + */ >> + const std::set<unsigned int>* _solve_only_on_dofs; > > I don't think a set is the correct data structure because the order > actually matters (it affects memory performance and the strength of > preconditioners). 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. I suggest to make two methods, one of which takes a vector and the other takes a set. Then, the set-taking method will by default just convert the set to a vector and call the vector-taking method. Roy, what do you think? > I also don't see why it needs to be kept separate from the IS, > you're duplicating information, it uses more memory and will have to > be kept consistent. Yes, you're right. >> + _solve_only_on_dofs = dofs; > > 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); Then, calling something like "delete &a;" inside the PetscLinearSolver class will certainly not be a good idea. Using a pointer as the argument is just to allow people to explicitely pass a NULL pointer. >> + { >> + ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(), >> + this->same_preconditioner ? SAME_PRECONDITIONER : >> SAME_NONZERO_PATTERN); >> + CHKERRABORT(libMesh::COMM_WORLD,ierr); >> + } > > Use of this->same_preconditioner is an unrelated change, correct? No, I just changed an if-else construction into a "?:" construction to make it less confusing (in connection with the additional if-else construction). The original code looked like this: 00477 if(!this->same_preconditioner) 00478 { 00479 ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(), 00480 SAME_NONZERO_PATTERN); 00481 CHKERRABORT(libMesh::COMM_WORLD,ierr); 00482 } 00483 else 00484 { 00485 ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(), 00486 SAME_PRECONDITIONER); 00487 CHKERRABORT(libMesh::COMM_WORLD,ierr); 00488 } The diff algorithm matched the "else" with a new, completely different "else", that's why it's difficult to see that I actually didn't change the functionality here. > Note that there are more than two values in that enum. Also, you > really have a DIFFERENT_NONZERO_PATTERN if you just switched > matrices to solve on a different index set. 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? >> + ierr = >> VecCreateSeq(libMesh::COMM_WORLD,_solve_only_on_dofs->size(),&subrhs); >> + CHKERRABORT(libMesh::COMM_WORLD,ierr); >> + >> + ierr = >> VecCreateSeq(libMesh::COMM_WORLD,_solve_only_on_dofs->size(),&subsolution); >> + CHKERRABORT(libMesh::COMM_WORLD,ierr); > > Why are you specializing this to only work in serial? I hoped that you would comment on this. (-: 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.) 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. 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. >> + ierr = VecScatterCreate(rhs->vec(),_solve_only_on_is, >> + subrhs,_solve_only_on_is_seq, >> + &scatter); >> + CHKERRABORT(libMesh::COMM_WORLD,ierr); > > Can use NULL for the second IS, it means "scatter to the whole vector". Thank you for that hint! >> +template <typename T> >> +void >> +LinearSolver<T>::solve_only_on(const std::set<unsigned int>* const dofs) >> +{ >> + if(dofs!=NULL) >> + { >> + libmesh_not_implemented(); >> + } >> +} > > If dofs == NULL, this function is a no-op? 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? (-: Best Regards, Tim -- Dr. Tim Kroeger CeVis -- Center of Complex Systems and Visualization University of Bremen [email protected] Universitaetsallee 29 [email protected] D-28359 Bremen Phone +49-421-218-7710 Germany Fax +49-421-218-4236 ------------------------------------------------------------------------------ Start uncovering the many advantages of virtual appliances and start using them to simplify application deployment and accelerate your shift to cloud computing. http://p.sf.net/sfu/novell-sfdev2dev _______________________________________________ Libmesh-users mailing list [email protected] https://lists.sourceforge.net/lists/listinfo/libmesh-users
