From: heister Date: Wed, 12 Oct 2011 18:52:57 +0000 (+0000) Subject: New: parallel::distributed::Triangulation::save()/load() to store X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=00189a534efdb8eca744943f0aa259434be233ab;p=dealii-svn.git New: parallel::distributed::Triangulation::save()/load() to store the refinement information to disk. Also supports saving solution vectors using the SolutionTransfer class. git-svn-id: https://svn.dealii.org/trunk@24595 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 5fcd80ed7b..9d9bf74773 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -34,7 +34,9 @@ inconvenience this causes.

General

    -
  1. None so far. +
  2. New: parallel::distributed::Triangulation::save()/load() to store +the refinement information to disk. Also supports saving solution vectors +using the SolutionTransfer class. (Timo Heister, 2011/10/12)
diff --git a/deal.II/include/deal.II/distributed/solution_transfer.h b/deal.II/include/deal.II/distributed/solution_transfer.h index 2d1229813f..cf30c9ed44 100644 --- a/deal.II/include/deal.II/distributed/solution_transfer.h +++ b/deal.II/include/deal.II/distributed/solution_transfer.h @@ -32,30 +32,69 @@ namespace parallel * coarsening a distributed grid and * handles the necessary communication. * - *

Usage

+ * @note: It is important to note, that if you use more than one + * SolutionTransfer object at the same time, that the calls to + * prepare_*() and interpolate()/deserialize() need to be in the same order. + * + *

Note on ghost elements

+ * In a parallel computation PETSc or Trilinos vector may contain ghost + * elements or not. For reading in information with + * prepare_for_coarsening_and_refinement() or prepare_serialization() you need + * to supply vectors with ghost elements, so that all locally_active elements + * can be read. On the other hand, ghosted vectors are generally not writable, + * so for calls to interpolate() or deserialize() you need to supply + * distributed vectors without ghost elements. + * + *

Transfering a solution

+ * Here VECTOR is your favorite vector type, e.g. PETScWrappers::MPI::Vector, + * TrilinosWrappers::MPI::Vector, or corresponding blockvectors. * @verbatim - * SolutionTransfer > soltrans(dof_handler); - * // flag some cells for refinement - * // and coarsening, e.g. - * GridRefinement::refine_and_coarsen_fixed_fraction( - * *tria, error_indicators, 0.3, 0.05); - * // prepare the triangulation, - * tria->prepare_coarsening_and_refinement(); - * // prepare the SolutionTransfer object - * // for coarsening and refinement and give - * // the solution vector that we intend to - * // interpolate later, - * soltrans.prepare_for_coarsening_and_refinement(solution); - * // actually execute the refinement, - * tria->execute_coarsening_and_refinement (); - * // redistribute dofs, - * dof_handler->distribute_dofs (fe); - * // and interpolate the solution - * Vector interpolated_solution(dof_handler->n_dofs()); - * soltrans.interpolate(interpolated_solution); - * @endverbatim + SolutionTransfer soltrans(dof_handler); + // flag some cells for refinement + // and coarsening, e.g. + GridRefinement::refine_and_coarsen_fixed_fraction( + tria, error_indicators, 0.3, 0.05); + // prepare the triangulation, + tria.prepare_coarsening_and_refinement(); + // prepare the SolutionTransfer object + // for coarsening and refinement and give + // the solution vector that we intend to + // interpolate later, + soltrans.prepare_for_coarsening_and_refinement(solution); + // actually execute the refinement, + tria.execute_coarsening_and_refinement (); + // redistribute dofs, + dof_handler.distribute_dofs (fe); + // and interpolate the solution + VECTOR interpolated_solution; + //create VECTOR in the right size here + soltrans.interpolate(interpolated_solution); +@endverbatim + * + *

Use for Serialization

+ * + * This class can be used to serialize and later deserialize a distributed mesh with solution vectors + * to a file. If you use more than one DoFHandler and therefore more than one SolutionTransfer object, they need to be serialized and deserialized in the same order. + * + * If vector has the locally relevant DoFs, serialization works as follows: + *@verbatim + +parallel::distributed::SolutionTransfer sol_trans(dof_handler); +sol_trans.prepare_serialization (vector); + +triangulation.save(filename); +@endverbatim +* For deserialization the vector needs to be a distributed vector (without ghost elements): +* @verbatim +//[create coarse mesh...] +triangulation.load(filename); + +parallel::distributed::SolutionTransfer sol_trans(dof_handler); +sol_trans.deserialize (distributed_vector); +@endverbatim + * * @ingroup distributed - * @author Timo Heister, 2009 + * @author Timo Heister, 2009-2011 */ template > class SolutionTransfer @@ -113,12 +152,49 @@ namespace parallel /** - * return the size in bytes that need + * Return the size in bytes that need * to be stored per cell. */ unsigned int get_data_size() const; + /** + * Prepare the serialization of the + * given vector. The serialization is + * done by Triangulation::save(). The + * given vector needs all information + * on the locally active DoFs (it must + * be ghosted). See documentation of + * this class for more information. + */ + void prepare_serialization(const VECTOR &in); + + + /** + * Same as the function above, only + * for a list of vectors. + */ + void prepare_serialization(const std::vector &all_in); + + + /** + * Execute the deserialization of the + * given vector. This needs to be + * done after calling + * Triangulation::load(). The given + * vector must be a fully distributed + * vector without ghost elements. See + * documentation of this class for + * more information. + */ + void deserialize(VECTOR &in); + + + /** + * Same as the function above, only + * for a list of vectors. + */ + void deserialize(std::vector &all_in); private: /** @@ -173,6 +249,12 @@ namespace parallel const void* data, std::vector &all_out); + + /** + * + */ + void register_data_attach(const std::size_t size); + }; diff --git a/deal.II/include/deal.II/distributed/tria.h b/deal.II/include/deal.II/distributed/tria.h index 1d98cd845a..1208278c45 100644 --- a/deal.II/include/deal.II/distributed/tria.h +++ b/deal.II/include/deal.II/distributed/tria.h @@ -436,6 +436,21 @@ namespace parallel */ unsigned int get_checksum () const; + /** + * Save the refinement information from the coarse mesh into the given + * file. This file needs to be reachable from all nodes in the computation + * on a shared network file system. See the SolutionTransfer class + * on how to store solution vectors into this file. + */ + void save(const char* filename) const; + + /** + * Load the refinement information saved with save() back in. The mesh + * must contain the same coarse mesh that was used in save(). You need + * to load with the same number of CPUs that you saved with. + */ + void load(const char* filename); + /** * Used to inform in the callbacks of * register_data_attach() and diff --git a/deal.II/source/distributed/solution_transfer.cc b/deal.II/source/distributed/solution_transfer.cc index d1253a3fdb..0bcc2658de 100644 --- a/deal.II/source/distributed/solution_transfer.cc +++ b/deal.II/source/distributed/solution_transfer.cc @@ -93,10 +93,20 @@ namespace parallel SolutionTransfer:: prepare_for_coarsening_and_refinement (const std::vector &all_in) { - Assert(all_in.size() > 0, ExcMessage("Please transfer atleast one vector!")); + input_vectors = all_in; - SolutionTransfer *ptr = this; + register_data_attach( get_data_size() * input_vectors.size() ); + } + + + template + void + SolutionTransfer:: + register_data_attach(const std::size_t size) + { + Assert(size > 0, ExcMessage("Please transfer at least one vector!")); + //TODO: casting away constness is bad parallel::distributed::Triangulation * tria = (dynamic_cast*> @@ -105,16 +115,16 @@ namespace parallel Assert (tria != 0, ExcInternalError()); offset - = tria->register_data_attach(static_cast - (get_data_size() * input_vectors.size()), + = tria->register_data_attach(size, std_cxx1x::bind(&SolutionTransfer::pack_callback, - ptr, + this, std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::_3)); - } - + + } + template void @@ -127,6 +137,52 @@ namespace parallel + template + void + SolutionTransfer:: + prepare_serialization(const VECTOR &in) + { + std::vector all_in(1, &in); + prepare_serialization(all_in); + } + + + + template + void + SolutionTransfer:: + prepare_serialization(const std::vector &all_in) + { + prepare_for_coarsening_and_refinement (all_in); + } + + + + template + void + SolutionTransfer:: + deserialize(VECTOR &in) + { + std::vector all_in(1, &in); + deserialize(all_in); + } + + + + template + void + SolutionTransfer:: + deserialize(std::vector &all_in) + { + register_data_attach( get_data_size() * all_in.size() ); + + // this makes interpolate() happy + input_vectors.resize(all_in.size()); + + interpolate(all_in); + } + + template void SolutionTransfer:: diff --git a/deal.II/source/distributed/tria.cc b/deal.II/source/distributed/tria.cc index de7b273f7f..63db8611f3 100644 --- a/deal.II/source/distributed/tria.cc +++ b/deal.II/source/distributed/tria.cc @@ -37,6 +37,7 @@ #include #include #include +#include DEAL_II_NAMESPACE_OPEN @@ -1947,7 +1948,92 @@ namespace parallel } + template + void + Triangulation:: + save(const char* filename) const + { + int real_data_size = 0; + if (attached_data_size>0) + real_data_size = attached_data_size+sizeof(CellStatus); + + if (my_subdomain==0) + { + std::string fname=std::string(filename)+".info"; + std::ofstream f(fname.c_str()); + f << Utilities::System::get_n_mpi_processes (mpi_communicator) << " " + << real_data_size << " " + << attached_data_pack_callbacks.size() << std::endl; + } + + if (attached_data_size>0) + { + const_cast*>(this) + ->attach_mesh_data(); + } + + dealii::internal::p4est::functions::save(filename, parallel_forest, attached_data_size>0); + + dealii::parallel::distributed::Triangulation* tria + = const_cast*>(this); + + tria->n_attached_datas = 0; + tria->attached_data_size = 0; + tria->attached_data_pack_callbacks.clear(); + + // and release the data + void * userptr = parallel_forest->user_pointer; + dealii::internal::p4est::functions::reset_data (parallel_forest, 0, NULL, NULL); + parallel_forest->user_pointer = userptr; + } + + + template + void + Triangulation:: + load(const char* filename) + { + dealii::internal::p4est::functions::destroy (parallel_forest); + parallel_forest = 0; + dealii::internal::p4est::functions::connectivity_destroy (connectivity); + connectivity=0; + + unsigned int numcpus, attached_size, attached_count; + { + std::string fname=std::string(filename)+".info"; + std::ifstream f(fname.c_str()); + f >> numcpus >> attached_size >> attached_count; + if (numcpus != Utilities::System::get_n_mpi_processes (mpi_communicator)) + throw ExcInternalError(); + } + + attached_data_size = 0; + n_attached_datas = 0; + + parallel_forest = dealii::internal::p4est::functions::load ( + filename, mpi_communicator, + attached_size, attached_size>0, + this, + &connectivity); + + try + { + copy_local_forest_to_triangulation (); + } + catch (const typename Triangulation::DistortedCellList &) + { + // the underlying + // triangulation should not + // be checking for + // distorted cells + AssertThrow (false, ExcInternalError()); + } + + update_number_cache (); + } + + template unsigned int Triangulation::get_checksum () const