Dear Roy, Thank you for your work and your patch.
I have looked at it and the following comments:(1) In the ghosted version of PetscVector::init(), you have placed an assert that checks whether the mesh is disjoint. Why is this necessary?
(2) In PetscVector::close(), I currently wonder why I used SCATTER_REVERSE in the call to VecGhostUpdateBegin/End(). This seems to be wrong, it should be SCATTER_FORWARD instead.
(3) It so happened that I looked at PetscVector<T>::operator=(const std::vector<T>& v) and noticed that in the ghosted case, it might be appropiate to set _is_closed to false. This is because for ghosted vectors, close() will also update the ghost values automatically. What do you think?
I attached a patch that fixes (2), but again, you might find it more convenient to perform that change manually. (The patch doesn't do anything else.)
On Wed, 4 Feb 2009, Roy Stogner wrote:
There's still something wrong with the localize() called from System::update(), but now it's not a bug so much as my failure to understand how PETSc is supposed to work: in a simple failing test case (the first update() in ex14 run on 2 procs), the index sets appear to be built correctly, but the scatter doesn't update current_local_solution the way I expect it to. The bug's either in localize() or in operator()(), both of which are fairly small simple functions but which are based on PETSc APIs I'm not familiar enough with.
Might it be that this was just due to my item (2) above? Actually, for ghosted vectors, I think localize() isn't required to be called at all, is it? I mean, close() will update the ghost values anyway, and there is nothing more to be updated.
Best Regards, Tim -- Dr. Tim Kroeger tim.kroe...@mevis.fraunhofer.de Phone +49-421-218-7710 tim.kroe...@cevis.uni-bremen.de Fax +49-421-218-4236 Fraunhofer MEVIS, Institute for Medical Image Computing Universitaetsallee 29, 28359 Bremen, Germany
Index: include/numerics/petsc_vector.h =================================================================== --- include/numerics/petsc_vector.h (Revision 3265) +++ include/numerics/petsc_vector.h (Arbeitskopie) @@ -775,9 +775,9 @@ if(this->type() == GHOSTED) { - ierr = VecGhostUpdateBegin(_vec,ADD_VALUES,SCATTER_REVERSE); + ierr = VecGhostUpdateBegin(_vec,ADD_VALUES,SCATTER_FORWARD); CHKERRABORT(libMesh::COMM_WORLD,ierr); - ierr = VecGhostUpdateEnd(_vec,ADD_VALUES,SCATTER_REVERSE); + ierr = VecGhostUpdateEnd(_vec,ADD_VALUES,SCATTER_FORWARD); CHKERRABORT(libMesh::COMM_WORLD,ierr); } Index: include/base/dof_map.h =================================================================== --- include/base/dof_map.h (Revision 3265) +++ include/base/dof_map.h (Arbeitskopie) @@ -582,7 +582,7 @@ * when projecting the solution onto the refined mesh, hence the parent's * DOF indices need to be included in the \p _send_list. */ - void augment_send_list_for_projection(const MeshBase &); + // void augment_send_list_for_projection(const MeshBase &); /** * Fills the vector di with the global degree of freedom indices @@ -681,12 +681,9 @@ * variable in the system. * Starts at index next_free_dof, and increments it to * the post-final index. - * If build_send_list is true, builds the send list. If - * false, clears and reserves the send list */ void distribute_local_dofs_var_major (unsigned int& next_free_dof, - MeshBase& mesh, - const bool build_send_list); + MeshBase& mesh); /** * Distributes the global degrees of freedom, for dofs on @@ -701,8 +698,7 @@ * false, clears and reserves the send list */ void distribute_local_dofs_node_major (unsigned int& next_free_dof, - MeshBase& mesh, - const bool build_send_list); + MeshBase& mesh); /** * Adds entries to the \p _send_list vector corresponding to DoFs Index: src/numerics/petsc_vector.C =================================================================== --- src/numerics/petsc_vector.C (Revision 3265) +++ src/numerics/petsc_vector.C (Arbeitskopie) @@ -443,18 +443,13 @@ libmesh_assert (this->_type == v._type); libmesh_assert (this->size() == v.size()); libmesh_assert (this->local_size() == v.local_size()); + libmesh_assert (this->_global_to_local_map == + v._global_to_local_map); -// this->init (v.size(), v.local_size()); -// this->_global_to_local_map = v._global_to_local_map; -// this->_is_closed = v._is_closed; -// this->_is_initialized = v._is_initialized; - if (v.size() != 0) { - int ierr = 0; - - ierr = VecCopy (v._vec, this->_vec); - CHKERRABORT(libMesh::COMM_WORLD,ierr); + int ierr = VecCopy (v._vec, this->_vec); + CHKERRABORT(libMesh::COMM_WORLD,ierr); } return *this; @@ -517,7 +512,7 @@ PetscVector<T>* v_local = libmesh_cast_ptr<PetscVector<T>*>(&v_local_in); libmesh_assert (v_local != NULL); - libmesh_assert (v_local->local_size() == this->size()); + libmesh_assert (v_local->size() == this->size()); int ierr = 0; const int n = this->size(); @@ -576,25 +571,29 @@ PetscVector<T>* v_local = libmesh_cast_ptr<PetscVector<T>*>(&v_local_in); libmesh_assert (v_local != NULL); - libmesh_assert (v_local->local_size() == this->size()); + libmesh_assert (v_local->size() == this->size()); libmesh_assert (send_list.size() <= v_local->size()); int ierr=0; - const int n_sl = send_list.size(); + const unsigned int n_sl = send_list.size(); IS is; VecScatter scatter; - std::vector<int> idx(n_sl); + std::vector<int> idx(n_sl + this->local_size()); - for (int i=0; i<n_sl; i++) + for (unsigned int i=0; i<n_sl; i++) idx[i] = static_cast<int>(send_list[i]); + for (unsigned int i = 0; i != this->local_size(); ++i) + idx[n_sl+i] = i + this->first_local_index(); // Create the index set & scatter object if (idx.empty()) - ierr = ISCreateGeneral(libMesh::COMM_WORLD, n_sl, PETSC_NULL, &is); + ierr = ISCreateGeneral(libMesh::COMM_WORLD, + n_sl+this->local_size(), PETSC_NULL, &is); else - ierr = ISCreateGeneral(libMesh::COMM_WORLD, n_sl, &idx[0], &is); + ierr = ISCreateGeneral(libMesh::COMM_WORLD, + n_sl+this->local_size(), &idx[0], &is); CHKERRABORT(libMesh::COMM_WORLD,ierr); ierr = VecScatterCreate(_vec, is, @@ -643,7 +642,7 @@ const std::vector<unsigned int>& send_list) { // Only good for serial vectors. - libmesh_assert (this->size() == this->local_size()); + // libmesh_assert (this->size() == this->local_size()); libmesh_assert (last_local_idx > first_local_idx); libmesh_assert (send_list.size() <= this->size()); libmesh_assert (last_local_idx < this->size()); Index: src/numerics/laspack_vector.C =================================================================== --- src/numerics/laspack_vector.C (Revision 3265) +++ src/numerics/laspack_vector.C (Arbeitskopie) @@ -358,7 +358,7 @@ libmesh_cast_ptr<LaspackVector<T>*>(&v_local_in); libmesh_assert (v_local != NULL); - libmesh_assert (send_list.size() == v_local->size()); + libmesh_assert (send_list.size() <= v_local->size()); *v_local = *this; } @@ -373,7 +373,7 @@ libmesh_assert (first_local_idx == 0); libmesh_assert (last_local_idx+1 == this->size()); - libmesh_assert (send_list.size() == this->size()); + libmesh_assert (send_list.size() <= this->size()); } Index: src/base/dof_map.C =================================================================== --- src/base/dof_map.C (Revision 3265) +++ src/base/dof_map.C (Arbeitskopie) @@ -668,9 +668,9 @@ // Set temporary DOF indices on this processor if (node_major_dofs) - this->distribute_local_dofs_node_major (next_free_dof, mesh, false); + this->distribute_local_dofs_node_major (next_free_dof, mesh); else - this->distribute_local_dofs_var_major (next_free_dof, mesh, false); + this->distribute_local_dofs_var_major (next_free_dof, mesh); // Get DOF counts on all processors std::vector<unsigned int> dofs_on_proc(n_proc, 0); @@ -694,9 +694,9 @@ // Set permanent DOF indices on this processor if (node_major_dofs) - this->distribute_local_dofs_node_major (next_free_dof, mesh, true); + this->distribute_local_dofs_node_major (next_free_dof, mesh); else - this->distribute_local_dofs_var_major (next_free_dof, mesh, true); + this->distribute_local_dofs_var_major (next_free_dof, mesh); libmesh_assert(next_free_dof == _end_df[proc_id]); @@ -725,6 +725,8 @@ STOP_LOG("distribute_dofs()", "DofMap"); + _send_list.clear(); + // Note that in the add_neighbors_to_send_list nodes on processor // boundaries that are shared by multiple elements are added for // each element. @@ -736,14 +738,13 @@ void DofMap::distribute_local_dofs_node_major(unsigned int &next_free_dof, - MeshBase& mesh, - const bool build_send_list) + MeshBase& mesh) { const unsigned int sys_num = this->sys_number(); const unsigned int n_vars = this->n_variables(); - // Number of elements to add to send_list - unsigned int send_list_size = 0; + // We now only add remote dofs to the _send_list + // unsigned int send_list_size = 0; // _var_first_local_df does not work with node_major dofs _var_first_local_df.resize(n_vars+1); @@ -783,17 +784,6 @@ 0, next_free_dof); next_free_dof += node->n_comp(sys_num,var); - - // If we've been requested to, add local dofs to the - // _send_list - if (build_send_list) - for (unsigned int index=0; index<node->n_comp(sys_num,var); - index++) - _send_list.push_back(node->dof_number(sys_num, - var, - index)); - else - send_list_size += node->n_comp(sys_num,var); } } } @@ -812,17 +802,6 @@ next_free_dof); next_free_dof += elem->n_comp(sys_num,var); - - // If we've been requested to, add local dofs to the - // _send_list - if (build_send_list) - for (unsigned int index=0; index<elem->n_comp(sys_num,var); - index++) - _send_list.push_back(elem->dof_number(sys_num, - var, - index)); - else - send_list_size += elem->n_comp(sys_num,var); } } // done looping over elements @@ -847,25 +826,18 @@ } } #endif // DEBUG - - if (!build_send_list) - { - _send_list.clear(); - _send_list.reserve(send_list_size); - } } void DofMap::distribute_local_dofs_var_major(unsigned int &next_free_dof, - MeshBase& mesh, - const bool build_send_list) + MeshBase& mesh) { const unsigned int sys_num = this->sys_number(); const unsigned int n_vars = this->n_variables(); - // Number of elements to add to send_list - unsigned int send_list_size = 0; + // We now only add remote dofs to the _send_list + // unsigned int send_list_size = 0; // We will cache the first local index for each variable _var_first_local_df.clear(); @@ -911,17 +883,6 @@ 0, next_free_dof); next_free_dof += node->n_comp(sys_num,var); - - // If we've been requested to, add local dofs to the - // _send_list - if (build_send_list) - for (unsigned int index=0; index<node->n_comp(sys_num,var); - index++) - _send_list.push_back(node->dof_number(sys_num, - var, - index)); - else - send_list_size += node->n_comp(sys_num,var); } } @@ -937,17 +898,6 @@ next_free_dof); next_free_dof += elem->n_comp(sys_num,var); - - // If we've been requested to, add local dofs to the - // _send_list - if (build_send_list) - for (unsigned int index=0; index<elem->n_comp(sys_num,var); - index++) - _send_list.push_back(elem->dof_number(sys_num, - var, - index)); - else - send_list_size += elem->n_comp(sys_num,var); } } // end loop on elements _var_first_local_df.push_back(next_free_dof); @@ -974,12 +924,6 @@ } } #endif // DEBUG - - if (!build_send_list) - { - _send_list.clear(); - _send_list.reserve(send_list_size); - } } @@ -989,7 +933,6 @@ START_LOG("add_neighbors_to_send_list()", "DofMap"); //------------------------------------------------------------------------- - // The _send_list now only contains entries from the local processor. // We need to add the DOFs from elements that live on neighboring processors // that are neighbors of the elements on the local processor //------------------------------------------------------------------------- @@ -1034,10 +977,12 @@ // Get the DOF indices for this neighboring element this->dof_indices (family[i], di); - // Insert the DOF indices into the send list - _send_list.insert (_send_list.end(), - di.begin(), di.end()); - } + // Insert the remote DOF indices into the send list + for (unsigned int j=0; j != di.size(); ++j) + if (di[j] < this->first_dof() || + di[j] >= this->end_dof()) + _send_list.push_back(di[j]); + } } } @@ -1072,9 +1017,11 @@ // Get the DOF indices for this neighboring element this->dof_indices (elem, di); - // Insert the DOF indices into the send list - _send_list.insert (_send_list.end(), - di.begin(), di.end()); + // Insert the remote DOF indices into the send list + for (unsigned int j=0; j != di.size(); ++j) + if (di[j] < this->first_dof() || + di[j] >= this->end_dof()) + _send_list.push_back(di[j]); } } @@ -1561,6 +1508,7 @@ +/* void DofMap::augment_send_list_for_projection(const MeshBase& mesh) { // Loop over the active local elements in the mesh. @@ -1612,6 +1560,7 @@ // The send-list might need to be sorted again. if (needs_sorting) this->sort_send_list (); } +*/ #endif // LIBMESH_ENABLE_AMR Index: src/solvers/transient_system.C =================================================================== --- src/solvers/transient_system.C (Revision 3265) +++ src/solvers/transient_system.C (Arbeitskopie) @@ -98,8 +98,12 @@ Base::init_data(); // Initialize the old & older solutions - old_local_solution->init (this->n_dofs(), false, SERIAL); - older_local_solution->init (this->n_dofs(), false, SERIAL); + old_local_solution->init (this->n_dofs(), this->n_local_dofs(), + this->get_dof_map().get_send_list(), + GHOSTED); + older_local_solution->init (this->n_dofs(), this->n_local_dofs(), + this->get_dof_map().get_send_list(), + GHOSTED); } Index: src/solvers/unsteady_solver.C =================================================================== --- src/solvers/unsteady_solver.C (Revision 3265) +++ src/solvers/unsteady_solver.C (Arbeitskopie) @@ -58,7 +58,9 @@ first_solve = false; } - old_local_nonlinear_solution->init (_system.n_dofs(), false, SERIAL); + old_local_nonlinear_solution->init (_system.n_dofs(), _system.n_local_dofs(), + _system.get_dof_map().get_send_list(), + GHOSTED); _system.get_vector("_old_nonlinear_solution").localize (*old_local_nonlinear_solution, Index: src/solvers/system.C =================================================================== --- src/solvers/system.C (Revision 3265) +++ src/solvers/system.C (Arbeitskopie) @@ -185,7 +185,9 @@ solution->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL); // Resize the current_local_solution for the current mesh - current_local_solution->init (this->n_dofs(), false, SERIAL); + current_local_solution->init (this->n_dofs(), this->n_local_dofs(), + _dof_map->get_send_list(), + GHOSTED); // from now on, no chance to add additional vectors _can_add_vectors = false; @@ -211,19 +213,24 @@ v->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL); } + const std::vector<unsigned int>& send_list = _dof_map->get_send_list (); + // Restrict the solution on the coarsened cells if (_solution_projection) { this->project_vector (*solution); current_local_solution->clear(); - current_local_solution->init(this->n_dofs()); - const std::vector<unsigned int>& send_list = _dof_map->get_send_list (); + current_local_solution->init(this->n_dofs(), + this->n_local_dofs(), send_list, + true, GHOSTED); solution->localize (*current_local_solution, send_list); } else { current_local_solution->clear(); - current_local_solution->init(this->n_dofs()); + current_local_solution->init(this->n_dofs(), + this->n_local_dofs(), send_list, + false, GHOSTED); } #endif } @@ -257,7 +264,8 @@ solution->init (this->n_dofs(), this->n_local_dofs(), true, PARALLEL); libmesh_assert (solution->size() == current_local_solution->size()); - libmesh_assert (solution->size() == current_local_solution->local_size()); + // Not true with ghosted vectors + // libmesh_assert (solution->size() == current_local_solution->local_size()); const unsigned int first_local_dof = solution->first_local_index(); const unsigned int local_size = solution->local_size(); @@ -275,7 +283,7 @@ const std::vector<unsigned int>& send_list = _dof_map->get_send_list (); // Check sizes - libmesh_assert (current_local_solution->local_size() == solution->size()); + libmesh_assert (current_local_solution->size() == solution->size()); // More processors than elements => empty send_list // libmesh_assert (!send_list.empty()); libmesh_assert (send_list.size() <= solution->size()); @@ -298,8 +306,9 @@ Utility::iota (send_list.begin(), send_list.end(), 0); // Check sizes - libmesh_assert (current_local_solution->local_size() == solution->size()); libmesh_assert (current_local_solution->size() == solution->size()); + // Not true with ghosted vectors + // libmesh_assert (current_local_solution->local_size() == solution->size()); libmesh_assert (!send_list.empty()); libmesh_assert (send_list.size() <= solution->size());
------------------------------------------------------------------------------ Create and Deploy Rich Internet Apps outside the browser with Adobe(R)AIR(TM) software. With Adobe AIR, Ajax developers can use existing skills and code to build responsive, highly engaging applications that combine the power of local resources and data with the reach of the web. Download the Adobe AIR SDK and Ajax docs to start building applications today-http://p.sf.net/sfu/adobe-com
_______________________________________________ Libmesh-devel mailing list Libmesh-devel@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/libmesh-devel