From: peterrum Date: Sun, 2 Jun 2019 14:35:48 +0000 (+0200) Subject: Add consensus algorithm and use it compute intex owners X-Git-Tag: v9.2.0-rc1~1404^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F8300%2Fhead;p=dealii.git Add consensus algorithm and use it compute intex owners --- diff --git a/doc/news/changes/minor/20190603PeterMunch b/doc/news/changes/minor/20190603PeterMunch new file mode 100644 index 0000000000..a3c2a4bbf9 --- /dev/null +++ b/doc/news/changes/minor/20190603PeterMunch @@ -0,0 +1,4 @@ +New: The function Utilities::MPI::compute_index_owner() provides the owner of +the ghost degrees of freedom relative to a set of locally owned ones. +
+(Peter Munch, 2019/06/03) diff --git a/doc/news/changes/minor/20190603PeterMunch-1 b/doc/news/changes/minor/20190603PeterMunch-1 new file mode 100644 index 0000000000..7e53962a9a --- /dev/null +++ b/doc/news/changes/minor/20190603PeterMunch-1 @@ -0,0 +1,5 @@ +Improved: The function Utilities::MPI::compute_index_owner() is used in +Utilities::MPI::Partitioner to avoid global all-to-all communication, +speeding up computations on many MPI ranks. +
+(Peter Munch, 2019/06/03) diff --git a/include/deal.II/base/mpi.h b/include/deal.II/base/mpi.h index 7068ce92fe..d5aa725229 100644 --- a/include/deal.II/base/mpi.h +++ b/include/deal.II/base/mpi.h @@ -23,6 +23,7 @@ #include #include +#include #include #if !defined(DEAL_II_WITH_MPI) && !defined(DEAL_II_WITH_PETSC) @@ -708,6 +709,495 @@ namespace Utilities const T & object_to_send, const unsigned int root_process = 0); + /** + * An interface to be able to use the ConsensusAlgorithm classes. The main + * functionality of the implementations is to return a list of process ranks + * this process wants data from and to deal with the optional payload of the + * messages sent/received by the ConsensusAlgorithm classes. + * + * There are two kinds of messages: + * - send/request message: A message consisting of a data request + * which should be answered by another process. This message is + * considered as a request message by the receiving rank. + * - recv message: The answer to a send/request message. + * + * @tparam T1 the type of the elements of the vector to sent + * @tparam T2 the type of the elements of the vector to received + * + * @note Since the payloads of the messages are optional, users have + * to deal with buffers themselves. The ConsensusAlgorithm classes 1) + * deliver only references to empty vectors (of size 0) + * the data to be sent can be inserted to or read from, and + * 2) communicate these vectors blindly. + * + * @author Peter Munch, 2019 + */ + template + class ConsensusAlgorithmProcess + { + public: + /** + * @return A vector of ranks this process wants to send a request to. + * + * @note This is the only method which has to be implemented since the + * payloads of the messages are optional. + */ + virtual std::vector + compute_targets() = 0; + + /** + * Add to the request to the process with the specified rank a payload. + * + * @param[in] other_rank rank of the process + * @param[out] send_buffer data to be sent part of the request (optional) + * + * @note The buffer is empty. Before using it, you have to set its size. + */ + virtual void + pack_recv_buffer(const int other_rank, std::vector &send_buffer); + + /** + * Prepare the buffer where the payload of the answer of the request to + * the process with the specified rank is saved in. The most obvious task + * is to resize the buffer, since it is empty when the function is called. + * + * @param[in] other_rank rank of the process + * @param[out] recv_buffer data to be sent part of the request (optional) + */ + virtual void + prepare_recv_buffer(const int other_rank, std::vector &recv_buffer); + + /** + * Prepare the buffer where the payload of the answer of the request to + * the process with the specified rank is saved in. + * + * @param[in] other_rank rank of the process + * @param[in] buffer_recv received payload (optional) + * @param[out] request_buffer payload to be sent as part of the request + * (optional) + * + * @note The request_buffer is empty. Before using it, you have to set + * its size. + */ + virtual void + process_request(const unsigned int other_rank, + const std::vector &buffer_recv, + std::vector & request_buffer); + + /** + * Process the payload of the answer of the request to the process with + * the specified rank. + * + * @param[in] other_rank rank of the process + * @param[in] recv_buffer data to be sent part of the request (optional) + */ + virtual void + unpack_recv_buffer(const int other_rank, + const std::vector &recv_buffer); + }; + + /** + * The task of the implementations of this class is to provide the + * communication patterns to retrieve data from other processes in a + * dynamic-sparse way. + * + * Dynamic-sparse means in this context: + * - by the time this function is called, the other processes do + * not know yet that they have to answer requests + * - each process only has to communicate with a small subset of + * processes of the MPI communicator + * + * Naturally, the user has to provide: + * - a communicator + * - for each rank a list of ranks of processes this process should + * communicate to + * - a functionality to pack/unpack data to be sent/received + * + * The latter two features should be implemented in a class derived from + * ConsensusAlgorithmProcess. + * + * @tparam T1 the type of the elements of the vector to sent + * @tparam T2 the type of the elements of the vector to received + * + * @author Peter Munch, 2019 + */ + template + class ConsensusAlgorithm + { + public: + ConsensusAlgorithm(ConsensusAlgorithmProcess &process, + const MPI_Comm & comm); + + /** + * Destructor. + */ + virtual ~ConsensusAlgorithm() = default; + + virtual void + run() = 0; + + protected: + /** + * Reference to the process provided by the user. + */ + ConsensusAlgorithmProcess &process; + + /** + * MPI communicator. + */ + const MPI_Comm &comm; + + /** + * Rank of this process. + */ + const unsigned int my_rank; + + /** + * Number of processes in the communicator. + */ + const unsigned int n_procs; + }; + + /** + * This class implements ConsensusAlgorithm, using only point-to-point + * communications and a single IBarrier. + * + * @note This class closely follows the paper Hoefner et. al. "Scalable + * Communication Protocols for Dynamic Sparse Data Exchange". + * Since the algorithm shown there is not considering payloads, the + * algorithm has been modified here in such a way that synchronous + * sends (Issend) have been replaced by equivalent Isend/Irecv, where + * Irecv receives the answer to a request (with payload). + * + * @tparam T1 the type of the elements of the vector to sent + * @tparam T2 the type of the elements of the vector to received + * + * @author Peter Munch, 2019 + */ + template + class ConsensusAlgorithm_NBX : public ConsensusAlgorithm + { + public: + // Unique tags to be used during Isend and Irecv + static const unsigned int tag_request = 12; + static const unsigned int tag_delivery = 13; + + /** + * Constructor. + * + * @param process Process to be run during consensus algorithm. + * @param comm MPI Communicator + */ + ConsensusAlgorithm_NBX(ConsensusAlgorithmProcess &process, + const MPI_Comm & comm); + + /** + * Destructor. + */ + virtual ~ConsensusAlgorithm_NBX() = default; + + /** + * Run consensus algorithm. + */ + virtual void + run() override; + + private: +#ifdef DEAL_II_WITH_MPI + /** + * List of processes this process wants to send requests to. + */ + std::vector targets; + + /** + * Buffers for sending requests. + */ + std::vector> send_buffers; + + /** + * Requests for sending requests. + */ + std::vector send_requests; + + /** + * Buffers for receiving answers to requests. + */ + std::vector> recv_buffers; + + + /** + * Requests for receiving answers to requests. + */ + std::vector recv_requests; + + /** + * Buffers for sending answers to requests. + */ + std::vector> request_buffers; + + /** + * Requests for sending answers to requests. + */ + std::vector> request_requests; + + // request for barrier + MPI_Request barrier_request; +#endif + +#ifdef DEBUG + /** + * List of processes who have made a request to this process. + */ + std::set requesting_processes; +#endif + + /** + * Check if all request answers have been received by this rank. + */ + bool + check_own_state(); + + /** + * Signal to all other ranks that this rank has received all request + * answers via entering IBarrier. + */ + void + signal_finish(); + + /** + * Check if all ranks have received all their request answers, i.e. + * all ranks have reached the IBarrier. + */ + bool + check_global_state(); + + /** + * A request message from another rank has been received: process the + * request and send an answer. + */ + void + process_requests(); + + /** + * Start to send all requests via ISend and post IRecvs for the incoming + * answer messages. + */ + void + start_communication(); + + /** + * After all rank has received all answers, the MPI data structures can be + * freed and the received answers can be processed. + */ + void + clean_up_and_end_communication(); + }; + + /** + * This class implements ConsensusAlgorithm, using a two step approach. In + * the first step the source ranks are determined and in the second step + * a static sparse data exchange is performed. + * + * @note In contrast to ConsensusAlgorithm_NBX, this class splits the same + * task into two distinct steps. In the first step, all processes are + * identified who want to send a request to this process. In the + * second step, the data is exchanged. However, since - in the second + * step - now it is clear how many requests have to be answered, i.e. + * when this process can stop waiting for requests, no IBarrier is + * needed. + * + * @note The function compute_point_to_point_communication_pattern() is used + * to determine the source processes, which implements a + * PEX-algorithm from Hoefner et. al. "Scalable Communication + * Protocols for Dynamic Sparse Data Exchange" + * + * @tparam T1 the type of the elements of the vector to sent + * @tparam T2 the type of the elements of the vector to received + * + * @author Peter Munch, 2019 + */ + template + class ConsensusAlgorithm_PEX : public ConsensusAlgorithm + { + public: + // Unique tags to be used during Isend and Irecv + static const unsigned int tag_request = 14; + static const unsigned int tag_delivery = 15; + + /** + * Constructor. + * + * @param process Process to be run during consensus algorithm. + * @param comm MPI Communicator + */ + ConsensusAlgorithm_PEX(ConsensusAlgorithmProcess &process, + const MPI_Comm & comm); + + /** + * Destructor. + */ + virtual ~ConsensusAlgorithm_PEX() = default; + + /** + * Run consensus algorithm. + */ + virtual void + run() override; + + private: +#ifdef DEAL_II_WITH_MPI + /** + * List of ranks of processes this processes wants to send a request to. + */ + std::vector targets; + + /** + * List of ranks of processes wanting to send a request to this process. + */ + std::vector sources; + + // data structures to send and receive requests + + /** + * Buffers for sending requests. + */ + std::vector> send_buffers; + + /** + * Buffers for receiving answers to requests. + */ + std::vector> recv_buffers; + + /** + * Requests for sending requests and receiving answers to requests. + */ + std::vector send_and_recv_buffers; + + /** + * Buffers for sending answers to requests. + */ + std::vector> requests_buffers; + + /** + * Requests for sending answers to requests. + */ + std::vector requests_answers; +#endif + + /** + * The ith request message from another rank has been received: process + * the request and send an answer. + */ + void + process_requests(int index); + + /** + * Start to send all requests via ISend and post IRecvs for the incoming + * answer messages. + */ + unsigned int + start_communication(); + + /** + * After all answers have been exchanged, the MPI data structures can be + * freed and the received answers can be processed. + */ + void + clean_up_and_end_communication(); + }; + + /** + * A class which delegates its task to other ConsensusAlgorithm + * implementations depending on the number of processes in the + * MPI communicator. For a small number of processes it uses + * ConsensusAlgorithm_PEX and for large number of processes + * ConsensusAlgorithm_NBX. The threshold depends if the program is + * compiled in debug or release mode. + * + * @tparam T1 the type of the elements of the vector to sent + * @tparam T2 the type of the elements of the vector to received + * + * @author Peter Munch, 2019 + */ + template + class ConsensusAlgorithmSelector : public ConsensusAlgorithm + { + public: + /** + * Constructor. + * + * @param process Process to be run during consensus algorithm. + * @param comm MPI Communicator. + */ + ConsensusAlgorithmSelector(ConsensusAlgorithmProcess &process, + const MPI_Comm & comm); + + /** + * Destructor. + */ + virtual ~ConsensusAlgorithmSelector() = default; + + /** + * Run consensus algorithm. The function call is delegated to another + * ConsensusAlgorithm implementation. + */ + virtual void + run() override; + + private: + // Pointer to the actual ConsensusAlgorithm implementation. + std::shared_ptr> consensus_algo; + }; + + /** + * Given a partitioned index set space, compute the owning MPI process rank + * of each element of a second index set according to the partitioned index + * set. A natural usage of this function is to compute for each ghosted + * degree of freedom the MPI rank of the process owning that index. + * + * One might think: "But we know which rank a ghost DoF belongs to based on + * the subdomain id of the cell it is on". But this heuristic fails for DoFs + * on interfaces between ghost cells with different subdomain_ids, or + * between a ghost cell and an artificial cell. Furthermore, this class + * enables a completely abstract exchange of information without the help of + * the grid in terms of neighbors. + * + * The first argument passed to this function, @p owned_indices, must + * uniquely partition an index space between all processes. + * Otherwise, there are no limitations on this argument: In particular, + * there is no need in partitioning + * the index space into contiguous subsets. Furthermore, there are no + * limitations + * on the second index set @p indices_to_look_up as long as the size matches + * the first one. It can be chosen arbitrarily and independently on each + * process. In the case that the second index set also contains locally + * owned indices, these indices will be treated correctly and the rank of + * this process is returned for those entries. + * + * @note This is a collective operation: all processes within the given + * communicator have to call this function. Since this function does not + * use MPI_Alltoall or MPI_Allgather, but instead uses non-blocking + * point-to-point communication instead, and only a single non-blocking + * barrier, it reduces the memory consumption significantly. This function + * is suited for large-scale simulations with >100k MPI ranks. + * + * @param[in] owned_indices Index set with indices locally owned by this + * process. + * @param[in] indices_to_look_up Index set containing indices of which the + * user is interested the rank of the owning process. + * @param[in] comm MPI communicator. + * + * @return List containing the MPI process rank for each entry in the index + * set @p indices_to_look_up. The order coincides with the order + * within the ElementIterator. + * + * @author Peter Munch, 2019 + */ + std::vector + compute_index_owner(const IndexSet &owned_indices, + const IndexSet &indices_to_look_up, + const MPI_Comm &comm); + #ifndef DOXYGEN // declaration for an internal function that lives in mpi.templates.h namespace internal diff --git a/source/base/mpi.cc b/source/base/mpi.cc index 81e9e32ebf..986d4007a3 100644 --- a/source/base/mpi.cc +++ b/source/base/mpi.cc @@ -27,6 +27,8 @@ #include #include +#include +#include #ifdef DEAL_II_WITH_TRILINOS # ifdef DEAL_II_WITH_MPI @@ -812,7 +814,938 @@ namespace Utilities #endif } + template + void + ConsensusAlgorithmProcess::process_request(const unsigned int, + const std::vector &, + std::vector &) + { + // noting to do + } + + + + template + void + ConsensusAlgorithmProcess::pack_recv_buffer(const int, + std::vector &) + { + // noting to do + } + + + + template + void + ConsensusAlgorithmProcess::prepare_recv_buffer(const int, + std::vector &) + { + // noting to do + } + + + + template + void + ConsensusAlgorithmProcess::unpack_recv_buffer( + const int, + const std::vector &) + { + // noting to do + } + + + + template + ConsensusAlgorithm::ConsensusAlgorithm( + ConsensusAlgorithmProcess &process, + const MPI_Comm & comm) + : process(process) + , comm(comm) + , my_rank(this_mpi_process(comm)) + , n_procs(n_mpi_processes(comm)) + {} + + + + template + ConsensusAlgorithm_NBX::ConsensusAlgorithm_NBX( + ConsensusAlgorithmProcess &process, + const MPI_Comm & comm) + : ConsensusAlgorithm(process, comm) + {} + + + + template + void + ConsensusAlgorithm_NBX::run() + { + // 1) send requests and start receiving the answers + start_communication(); + + // 2) answer requests and check if all requests of this process have been + // answered + while (!check_own_state()) + process_requests(); + + // 3) signal to all other processes that all requests of this process have + // been answered + signal_finish(); + + // 4) nevertheless, this process has to keep on answering (potential) + // incoming requests until all processes have received the + // answer to all requests + while (!check_global_state()) + process_requests(); + + // 5) process the answer to all requests + clean_up_and_end_communication(); + } + + + + template + bool + ConsensusAlgorithm_NBX::check_own_state() + { +#ifdef DEAL_II_WITH_MPI + int all_receive_requests_are_done; + const auto ierr = MPI_Testall(recv_requests.size(), + recv_requests.data(), + &all_receive_requests_are_done, + MPI_STATUSES_IGNORE); + AssertThrowMPI(ierr); + + return all_receive_requests_are_done; +#else + return true; +#endif + } + + + + template + void + ConsensusAlgorithm_NBX::signal_finish() + { +#ifdef DEAL_II_WITH_MPI +# if DEAL_II_MPI_VERSION_GTE(3, 0) + const auto ierr = MPI_Ibarrier(this->comm, &barrier_request); + AssertThrowMPI(ierr); +# else + AssertThrow( + false, + ExcMessage( + "ConsensusAlgorithm_NBX uses MPI 3.0 features. You should compile with at least MPI 3.0.")); +# endif +#endif + } + + + + template + bool + ConsensusAlgorithm_NBX::check_global_state() + { +#ifdef DEAL_II_WITH_MPI + int all_ranks_reached_barrier; + const auto ierr = MPI_Test(&barrier_request, + &all_ranks_reached_barrier, + MPI_STATUSES_IGNORE); + AssertThrowMPI(ierr); + return all_ranks_reached_barrier; +#else + return true; +#endif + } + + + + template + void + ConsensusAlgorithm_NBX::process_requests() + { +#ifdef DEAL_II_WITH_MPI + // check if there is a request pending + MPI_Status status; + int request_is_pending; + const auto ierr = MPI_Iprobe( + MPI_ANY_SOURCE, tag_request, this->comm, &request_is_pending, &status); + AssertThrowMPI(ierr); + + if (request_is_pending) // request is pending + { + // get rank of requesting process + const auto other_rank = status.MPI_SOURCE; + +# ifdef DEBUG + Assert(requesting_processes.find(other_rank) == + requesting_processes.end(), + ExcMessage("Process is requesting a second time!")); + requesting_processes.insert(other_rank); +# endif + + std::vector buffer_recv; + // get size of of incoming message + int number_amount; + auto ierr = MPI_Get_count(&status, + internal::mpi_type_id(buffer_recv.data()), + &number_amount); + AssertThrowMPI(ierr); + + // allocate memory for incoming message + buffer_recv.resize(number_amount); + ierr = MPI_Recv(buffer_recv.data(), + number_amount, + internal::mpi_type_id(buffer_recv.data()), + other_rank, + tag_request, + this->comm, + &status); + AssertThrowMPI(ierr); + + // allocate memory for answer message + request_buffers.emplace_back(); + request_requests.emplace_back(new MPI_Request); + + // process request + auto &request_buffer = request_buffers.back(); + this->process.process_request(other_rank, + buffer_recv, + request_buffer); + + // start to send answer back + ierr = MPI_Isend(request_buffer.data(), + request_buffer.size(), + internal::mpi_type_id(request_buffer.data()), + other_rank, + tag_delivery, + this->comm, + request_requests.back().get()); + AssertThrowMPI(ierr); + } +#endif + } + + + + template + void + ConsensusAlgorithm_NBX::start_communication() + { +#ifdef DEAL_II_WITH_MPI + // 1) + targets = this->process.compute_targets(); + const auto n_targets = targets.size(); + + // 2) allocate memory + recv_buffers.resize(n_targets); + recv_requests.resize(n_targets); + send_requests.resize(n_targets); + send_buffers.resize(n_targets); + + { + // 4) send and receive + for (unsigned int i = 0; i < n_targets; i++) + { + const unsigned int rank = targets[i]; + const unsigned int index = i; + + // translate index set to a list of pairs + auto &send_buffer = send_buffers[index]; + this->process.pack_recv_buffer(rank, send_buffer); + + // start to send data + auto ierr = MPI_Isend(send_buffer.data(), + send_buffer.size(), + internal::mpi_type_id(send_buffer.data()), + rank, + tag_request, + this->comm, + &send_requests[index]); + AssertThrowMPI(ierr); + + // start to receive data + auto &recv_buffer = recv_buffers[index]; + this->process.prepare_recv_buffer(rank, recv_buffer); + ierr = MPI_Irecv(recv_buffer.data(), + recv_buffer.size(), + internal::mpi_type_id(recv_buffer.data()), + rank, + tag_delivery, + this->comm, + &recv_requests[index]); + AssertThrowMPI(ierr); + } + } +#endif + } + + + + template + void + ConsensusAlgorithm_NBX::clean_up_and_end_communication() + { +#ifdef DEAL_II_WITH_MPI + // clean up + { + auto ierr = MPI_Waitall(send_requests.size(), + send_requests.data(), + MPI_STATUSES_IGNORE); + AssertThrowMPI(ierr); + + ierr = MPI_Waitall(recv_requests.size(), + recv_requests.data(), + MPI_STATUSES_IGNORE); + AssertThrowMPI(ierr); + + ierr = MPI_Wait(&barrier_request, MPI_STATUS_IGNORE); + AssertThrowMPI(ierr); + + for (auto &i : request_requests) + { + const auto ierr = MPI_Wait(i.get(), MPI_STATUS_IGNORE); + AssertThrowMPI(ierr); + } + +# ifdef DEBUG + // note: IBarrier seems to make problem during testing, this additional + // Barrier seems to help + MPI_Barrier(this->comm); +# endif + } + + // unpack data + { + for (unsigned int i = 0; i < targets.size(); i++) + this->process.unpack_recv_buffer(targets[i], recv_buffers[i]); + } +#endif + } + + + + /** + * A re-implementation of compute_point_to_point_communication_pattern + * using the ConsensusAlgorithm. + */ + class ConsensusAlgorithmProcessTargets + : public ConsensusAlgorithmProcess + { + public: + ConsensusAlgorithmProcessTargets(std::vector &target) + : target(target) + {} + + using T1 = int; + using T2 = int; + + virtual void + process_request(const unsigned int other_rank, + const std::vector &buffer_recv, + std::vector & request_buffer) override + { + (void)buffer_recv; + (void)request_buffer; + this->sources.push_back(other_rank); + } + + /** + * Simply return the user-provided list. + * + * @return List of processes this process wants to send requests to. + */ + virtual std::vector + compute_targets() override + { + return target; + } + + /** + * The result of the consensus algorithm. + * @return Sorted list of ranks of processes wanting to send a request to + * this process. + */ + std::vector + get_result() + { + std::sort(sources.begin(), sources.end()); + return sources; + } + + private: + /** + * List of processes this process wants to send requests to. + */ + const std::vector ⌖ + + /** + * List of ranks of processes wanting to send a request to this process. + */ + std::vector sources; + }; + + + + template + ConsensusAlgorithm_PEX::ConsensusAlgorithm_PEX( + ConsensusAlgorithmProcess &process, + const MPI_Comm & comm) + : ConsensusAlgorithm(process, comm) + {} + + + + template + void + ConsensusAlgorithm_PEX::run() + { + // 1) send requests and start receiving the answers + // especially determine how many requests are expected + const unsigned int n_requests = start_communication(); + + // 2) answer requests + for (unsigned int request = 0; request < n_requests; request++) + process_requests(request); + + // 3) process answers + clean_up_and_end_communication(); + } + + + + template + void + ConsensusAlgorithm_PEX::process_requests(int index) + { +#ifdef DEAL_II_WITH_MPI + MPI_Status status; + MPI_Probe(MPI_ANY_SOURCE, tag_request, this->comm, &status); + + // get rank of incoming message + const auto other_rank = status.MPI_SOURCE; + + std::vector buffer_recv; + + // get size of incoming message + int number_amount; + auto ierr = MPI_Get_count(&status, + internal::mpi_type_id(buffer_recv.data()), + &number_amount); + AssertThrowMPI(ierr); + + // allocate memory for incoming message + buffer_recv.resize(number_amount); + ierr = MPI_Recv(buffer_recv.data(), + number_amount, + internal::mpi_type_id(buffer_recv.data()), + other_rank, + tag_request, + this->comm, + &status); + AssertThrowMPI(ierr); + + // process request + auto &request_buffer = requests_buffers[index]; + this->process.process_request(other_rank, buffer_recv, request_buffer); + + // start to send answer back + ierr = MPI_Isend(request_buffer.data(), + request_buffer.size(), + MPI_UNSIGNED, + other_rank, + tag_delivery, + this->comm, + &requests_answers[index]); + AssertThrowMPI(ierr); +#else + (void)index; +#endif + } + + + template + unsigned int + ConsensusAlgorithm_PEX::start_communication() + { +#ifdef DEAL_II_WITH_MPI + // 1) determine with which processes this process wants to communicate + targets = this->process.compute_targets(); + + // 2) determine who wants to communicate with this process + const bool use_nbx = false; + if (!use_nbx) + { + sources = + compute_point_to_point_communication_pattern(this->comm, targets); + } + else + { + ConsensusAlgorithmProcessTargets process(targets); + ConsensusAlgorithm_NBX + consensus_algorithm(process, this->comm); + consensus_algorithm.run(); + sources = process.get_result(); + } + + const auto n_targets = targets.size(); + const auto n_sources = sources.size(); + + // 2) allocate memory + recv_buffers.resize(n_targets); + send_buffers.resize(n_targets); + send_and_recv_buffers.resize(2 * n_targets); + + requests_answers.resize(n_sources); + requests_buffers.resize(n_sources); + + // 4) send and receive + for (unsigned int i = 0; i < n_targets; i++) + { + const unsigned int rank = targets[i]; + + // pack data which should be sent + auto &send_buffer = send_buffers[i]; + this->process.pack_recv_buffer(rank, send_buffer); + + // start to send data + auto ierr = MPI_Isend(send_buffer.data(), + send_buffer.size(), + internal::mpi_type_id(send_buffer.data()), + rank, + tag_request, + this->comm, + &send_and_recv_buffers[n_targets + i]); + AssertThrowMPI(ierr); + + // start to receive data + auto &recv_buffer = recv_buffers[i]; + this->process.prepare_recv_buffer(rank, recv_buffer); + ierr = MPI_Irecv(recv_buffer.data(), + recv_buffer.size(), + MPI_UNSIGNED, + rank, + tag_delivery, + this->comm, + &send_and_recv_buffers[i]); + AssertThrowMPI(ierr); + } + + return sources.size(); +#else + return 0; +#endif + } + + + + template + void + ConsensusAlgorithm_PEX::clean_up_and_end_communication() + { +#ifdef DEAL_II_WITH_MPI + // finalize all MPI_Requests + MPI_Waitall(send_and_recv_buffers.size(), + send_and_recv_buffers.data(), + MPI_STATUSES_IGNORE); + MPI_Waitall(requests_answers.size(), + requests_answers.data(), + MPI_STATUSES_IGNORE); + + // unpack received data + for (unsigned int i = 0; i < targets.size(); i++) + this->process.unpack_recv_buffer(targets[i], recv_buffers[i]); +#endif + } + + + + template + ConsensusAlgorithmSelector::ConsensusAlgorithmSelector( + ConsensusAlgorithmProcess &process, + const MPI_Comm & comm) + : ConsensusAlgorithm(process, comm) + { + // Depending on the number of processes we switch between implementations. + // We reduce the threshold for debug mode to be able to test also the + // non-blocking implementation. This feature is tested by: + // tests/multigrid/transfer_matrix_free_06.with_mpi=true.with_p4est=true.with_trilinos=true.mpirun=15.output +#ifdef DEAL_II_WITH_MPI +# if DEAL_II_MPI_VERSION_GTE(3, 0) +# ifdef DEBUG + if (Utilities::MPI::n_mpi_processes(comm) > 14) +# else + if (Utilities::MPI::n_mpi_processes(comm) > 99) +# endif + consensus_algo.reset(new ConsensusAlgorithm_NBX(process, comm)); + else +# endif +#endif + consensus_algo.reset(new ConsensusAlgorithm_PEX(process, comm)); + } + + + + template + void + ConsensusAlgorithmSelector::run() + { + consensus_algo->run(); + } + + + + namespace ComputeIndexOwner + { + struct Dictionary + { + static const unsigned int tag_setup = 11; + + std::vector actually_owning_ranks; + + types::global_dof_index dofs_per_process; + std::pair local_range; + types::global_dof_index local_size; + types::global_dof_index size; + + void + reinit(const IndexSet &owned_indices, const MPI_Comm &comm) + { + this->partition(owned_indices, comm); + +#ifdef DEAL_II_WITH_MPI + unsigned int my_rank = this_mpi_process(comm); + + types::global_dof_index dic_local_rececived = 0; + std::map relevant_procs_map; + + // 2) collect relevant processes and process local dict entries + { + std::vector relevant_procs; + for (auto i : owned_indices) + { + unsigned int other_rank = this->dof_to_dict_rank(i); + if (other_rank == my_rank) + { + this->actually_owning_ranks[i - this->local_range.first] = + my_rank; + dic_local_rececived++; + } + else if (relevant_procs.empty() || + relevant_procs.back() != other_rank) + relevant_procs.push_back(other_rank); + } + + { + unsigned int c = 0; + for (auto i : relevant_procs) + relevant_procs_map[i] = c++; + } + } + + const unsigned int n_relevant_procs = relevant_procs_map.size(); + std::vector>> + buffers(n_relevant_procs); + std::vector request(n_relevant_procs); + + // 3) send messages with local dofs to the right dict process + { + std::vector> temp( + n_relevant_procs); + + // collect dofs of each dict process + for (auto i : owned_indices) + { + unsigned int other_rank = this->dof_to_dict_rank(i); + if (other_rank != my_rank) + temp[relevant_procs_map[other_rank]].push_back(i); + } + + // send dofs to each process + for (auto rank_pair : relevant_procs_map) + { + const int rank = rank_pair.first; + const int index = rank_pair.second; + + // create index set and compress data to be sent + auto & indices_i = temp[index]; + IndexSet is(this->size); + is.add_indices(indices_i.begin(), indices_i.end()); + is.compress(); + + // translate index set to a list of pairs + auto &buffer = buffers[index]; + for (auto interval = is.begin_intervals(); + interval != is.end_intervals(); + interval++) + buffer.emplace_back(*interval->begin(), interval->last() + 1); + + // send data + const auto ierr = MPI_Isend(buffer.data(), + buffer.size() * 2, + DEAL_II_DOF_INDEX_MPI_TYPE, + rank, + tag_setup, + comm, + &request[index]); + AssertThrowMPI(ierr); + } + } + + + // 4) receive messages until all dofs in dict are processed + while (this->local_size != dic_local_rececived) + { + // wait for an incoming message + MPI_Status status; + auto ierr = MPI_Probe(MPI_ANY_SOURCE, tag_setup, comm, &status); + AssertThrowMPI(ierr); + + // retrieve size of incoming message + int number_amount; + ierr = MPI_Get_count(&status, + DEAL_II_DOF_INDEX_MPI_TYPE, + &number_amount); + AssertThrowMPI(ierr); + + const auto other_rank = status.MPI_SOURCE; + + // receive message + Assert(number_amount % 2 == 0, ExcInternalError()); + std::vector< + std::pair> + buffer(number_amount / 2); + ierr = MPI_Recv(buffer.data(), + number_amount, + DEAL_II_DOF_INDEX_MPI_TYPE, + other_rank, + tag_setup, + comm, + &status); + AssertThrowMPI(ierr); + + // process message: loop over all intervals + for (auto interval : buffer) + for (types::global_dof_index i = interval.first; + i < interval.second; + i++) + { + this->actually_owning_ranks[i - this->local_range.first] = + other_rank; + dic_local_rececived++; + } + } + + // 5) make sure that all messages have been sent + const auto ierr = + MPI_Waitall(n_relevant_procs, request.data(), MPI_STATUSES_IGNORE); + AssertThrowMPI(ierr); +#else + (void)owned_indices; + (void)comm; +#endif + } + + unsigned int + dof_to_dict_rank(const types::global_dof_index i) + { + return i / dofs_per_process; + } + + private: + void + partition(const IndexSet &owned_indices, const MPI_Comm &comm) + { +#ifdef DEAL_II_WITH_MPI + const unsigned int n_procs = n_mpi_processes(comm); + const unsigned int my_rank = this_mpi_process(comm); + + size = owned_indices.size(); + dofs_per_process = (size + n_procs - 1) / n_procs; + local_range.first = std::min(dofs_per_process * my_rank, size); + local_range.second = std::min(dofs_per_process * (my_rank + 1), size); + local_size = local_range.second - local_range.first; + + actually_owning_ranks.resize(local_size); +#else + (void)owned_indices; + (void)comm; +#endif + } + }; + + class ConsensusAlgorithmProcess + : public dealii::Utilities::MPI:: + ConsensusAlgorithmProcess + { + public: + ConsensusAlgorithmProcess(const IndexSet & owned_indices, + const IndexSet & indices_to_look_up, + const MPI_Comm & comm, + std::vector &owning_ranks) + : owned_indices(owned_indices) + , indices_to_look_up(indices_to_look_up) + , comm(comm) + , my_rank(this_mpi_process(comm)) + , n_procs(n_mpi_processes(comm)) + , owning_ranks(owning_ranks) + { + this->dict.reinit(owned_indices, comm); + } + + const IndexSet & owned_indices; + const IndexSet & indices_to_look_up; + const MPI_Comm & comm; + const unsigned int my_rank; + const unsigned int n_procs; + std::vector &owning_ranks; + + Dictionary dict; + + std::map> temp; + std::map> recv_indices; + + virtual void + process_request(const unsigned int other_rank, + const std::vector &buffer_recv, + std::vector &request_buffer) override + { + (void)other_rank; + Assert(buffer_recv.size() % 2 == 0, ExcInternalError()); + for (unsigned int j = 0; j < buffer_recv.size(); j += 2) + for (auto i = buffer_recv[j]; i < buffer_recv[j + 1]; i++) + request_buffer.push_back( + dict.actually_owning_ranks[i - dict.local_range.first]); + } + + virtual std::vector + compute_targets() override + { + std::vector targets; + + // 1) collect relevant processes and process local dict entries + { + unsigned int index = 0; + for (auto i : indices_to_look_up) + { + unsigned int other_rank = dict.dof_to_dict_rank(i); + if (other_rank == my_rank) + owning_ranks[index] = + dict.actually_owning_ranks[i - dict.local_range.first]; + else if (targets.empty() || targets.back() != other_rank) + targets.push_back(other_rank); + index++; + } + } + + + for (auto i : targets) + { + recv_indices[i] = {}; + temp[i] = {}; + } + + // 3) collect indices for each process + { + unsigned int index = 0; + for (auto i : indices_to_look_up) + { + unsigned int other_rank = dict.dof_to_dict_rank(i); + if (other_rank != my_rank) + { + recv_indices[other_rank].push_back(index); + temp[other_rank].push_back(i); + } + index++; + } + } + + Assert(targets.size() == recv_indices.size() && + targets.size() == temp.size(), + ExcMessage("Size does not match!")); + + return targets; + } + + virtual void + pack_recv_buffer( + const int other_rank, + std::vector &send_buffer) override + { + // create index set and compress data to be sent + auto & indices_i = temp[other_rank]; + IndexSet is(dict.size); + is.add_indices(indices_i.begin(), indices_i.end()); + is.compress(); + + for (auto interval = is.begin_intervals(); + interval != is.end_intervals(); + interval++) + { + send_buffer.push_back(*interval->begin()); + send_buffer.push_back(interval->last() + 1); + } + } + + virtual void + prepare_recv_buffer(const int other_rank, + std::vector &recv_buffer) override + { + recv_buffer.resize(recv_indices[other_rank].size()); + } + + virtual void + unpack_recv_buffer( + const int other_rank, + const std::vector &recv_buffer) override + { + Assert(recv_indices[other_rank].size() == recv_buffer.size(), + ExcMessage("Sizes do not match!")); + + for (unsigned int j = 0; j < recv_indices[other_rank].size(); j++) + owning_ranks[recv_indices[other_rank][j]] = recv_buffer[j]; + } + }; + + } // namespace ComputeIndexOwner + + + + std::vector + compute_index_owner(const IndexSet &owned_indices, + const IndexSet &indices_to_look_up, + const MPI_Comm &comm) + { + Assert(owned_indices.size() == indices_to_look_up.size(), + ExcMessage("IndexSets have to have the same sizes.")); + + std::vector owning_ranks(indices_to_look_up.n_elements()); + + // Step 1: setup dictionary + // The input owned_indices can be partitioned arbitrarily. In the + // dictionary, the index set is statically repartitioned among the + // processes again and extended with information with the actual owner + // of that the index. + ComputeIndexOwner::ConsensusAlgorithmProcess process(owned_indices, + indices_to_look_up, + comm, + owning_ranks); + + // Step 2: read dictionary + // Communicate with the process who owns the index in the static + // partition (i.e. in the partition). This process returns the actual + // owner of the index. + ConsensusAlgorithmSelector + consensus_algorithm(process, comm); + consensus_algorithm.run(); + + return owning_ranks; + } #include "mpi.inst" } // end of namespace MPI diff --git a/source/base/partitioner.cc b/source/base/partitioner.cc index 673c52f450..9f94b2090a 100644 --- a/source/base/partitioner.cc +++ b/source/base/partitioner.cc @@ -233,56 +233,32 @@ namespace Utilities } } - // Allocate memory for data that will be exported - std::vector expanded_ghost_indices( - n_ghost_indices_data); unsigned int n_ghost_targets = 0; - if (n_ghost_indices_data > 0) - { - // Create first a vector of ghost_targets from the list of ghost - // indices and then push back new values. When we are done, copy the - // data to that field of the partitioner. This way, the variable - // ghost_targets will have exactly the size we need, whereas the - // vector filled with emplace_back might actually be too long. - unsigned int current_proc = 0; - ghost_indices_data.fill_index_vector(expanded_ghost_indices); - types::global_dof_index current_index = expanded_ghost_indices[0]; - while (current_index >= first_index[current_proc + 1]) - current_proc++; - AssertIndexRange(current_proc, n_procs); - - // since DoFs are contiguous, populate a vector which stores - // a process rank and the number of ghosts - std::vector> ghost_targets_temp( - 1, std::pair(current_proc, 0)); - n_ghost_targets++; - - // find which process is the owner of other indices: - for (unsigned int iterator = 1; iterator < n_ghost_indices_data; - ++iterator) - { - current_index = expanded_ghost_indices[iterator]; - while (current_index >= first_index[current_proc + 1]) - current_proc++; - AssertIndexRange(current_proc, n_procs); - // if we found a new target (i.e. higher rank) then adjust the - // pair.second in the last element so that it stores the total - // number of ghosts owned by pair.first - if (ghost_targets_temp[n_ghost_targets - 1].first < current_proc) - { - ghost_targets_temp[n_ghost_targets - 1].second = - iterator - ghost_targets_temp[n_ghost_targets - 1].second; - ghost_targets_temp.emplace_back(current_proc, iterator); - n_ghost_targets++; - } - } - // adjust the last element so that pair.second stores the number of - // elements - ghost_targets_temp[n_ghost_targets - 1].second = - n_ghost_indices_data - - ghost_targets_temp[n_ghost_targets - 1].second; - ghost_targets_data = ghost_targets_temp; - } + { + const auto index_owner = + Utilities::MPI::compute_index_owner(this->locally_owned_range_data, + this->ghost_indices_data, + this->communicator); + + ghost_targets_data.clear(); + + if (index_owner.size() > 0) + { + ghost_targets_data.emplace_back(index_owner[0], 0); + for (auto i : index_owner) + { + Assert(i >= ghost_targets_data.back().first, + ExcInternalError( + "Expect result of compute_index_owner to be sorted")); + if (i == ghost_targets_data.back().first) + ghost_targets_data.back().second++; + else + ghost_targets_data.emplace_back(i, 1); + } + } + + n_ghost_targets = ghost_targets_data.size(); + } // find the processes that want to import to me { std::vector send_buffer(n_procs, 0); @@ -316,9 +292,9 @@ namespace Utilities // now that we know how many indices each process will receive from // ghosts, send and receive indices for import data. non-blocking receives // and blocking sends - std::vector expanded_import_indices( - n_import_indices_data); { + std::vector expanded_import_indices( + n_import_indices_data); unsigned int current_index_start = 0; std::vector import_requests(import_targets_data.size() + n_ghost_targets); @@ -339,6 +315,10 @@ namespace Utilities // use non-blocking send for ghost indices stored in // expanded_ghost_indices + std::vector expanded_ghost_indices; + if (n_ghost_indices_data > 0) + ghost_indices_data.fill_index_vector(expanded_ghost_indices); + current_index_start = 0; for (unsigned int i = 0; i < n_ghost_targets; i++) { diff --git a/source/multigrid/mg_transfer_internal.cc b/source/multigrid/mg_transfer_internal.cc index 71918c5e41..1edb807932 100644 --- a/source/multigrid/mg_transfer_internal.cc +++ b/source/multigrid/mg_transfer_internal.cc @@ -54,8 +54,6 @@ namespace internal {} }; - - template void fill_copy_indices( @@ -197,24 +195,66 @@ namespace internal tria->level_ghost_owners(); std::map> send_data; - // * find owners of the level dofs and insert into send_data - // accordingly - for (const auto &dofpair : send_data_temp) + std::sort(send_data_temp.begin(), + send_data_temp.end(), + [](const DoFPair &lhs, const DoFPair &rhs) { + if (lhs.level < rhs.level) + return true; + if (lhs.level > rhs.level) + return false; + + if (lhs.level_dof_index < rhs.level_dof_index) + return true; + if (lhs.level_dof_index > rhs.level_dof_index) + return false; + + if (lhs.global_dof_index < rhs.global_dof_index) + return true; + else + return false; + }); + send_data_temp.erase( + std::unique(send_data_temp.begin(), + send_data_temp.end(), + [](const DoFPair &lhs, const DoFPair &rhs) { + return (lhs.level == rhs.level) && + (lhs.level_dof_index == rhs.level_dof_index) && + (lhs.global_dof_index == rhs.global_dof_index); + }), + send_data_temp.end()); + + for (unsigned int level = 0; level < n_levels; ++level) { - std::set::iterator it; - for (it = neighbors.begin(); it != neighbors.end(); ++it) - { - if (locally_owned_mg_dofs_per_processor[dofpair.level][*it] - .is_element(dofpair.level_dof_index)) - { - send_data[*it].push_back(dofpair); - break; - } - } - // Is this level dof not owned by any of our neighbors? That would - // certainly be a bug! - Assert(it != neighbors.end(), - ExcMessage("could not find DoF owner.")); + const IndexSet &is_local = mg_dof.locally_owned_mg_dofs(level); + + std::vector level_dof_indices; + std::vector global_dof_indices; + for (const auto &dofpair : send_data_temp) + if (dofpair.level == level) + { + level_dof_indices.push_back(dofpair.level_dof_index); + global_dof_indices.push_back(dofpair.global_dof_index); + } + + IndexSet is_ghost(is_local.size()); + is_ghost.add_indices(level_dof_indices.begin(), + level_dof_indices.end()); + + AssertThrow(level_dof_indices.size() == is_ghost.n_elements(), + ExcMessage("Size does not match!")); + + const auto index_owner = + Utilities::MPI::compute_index_owner(is_local, + is_ghost, + tria->get_communicator()); + + AssertThrow(level_dof_indices.size() == index_owner.size(), + ExcMessage("Size does not match!")); + + for (unsigned int i = 0; i < index_owner.size(); i++) + send_data[index_owner[i]].emplace_back(level, + global_dof_indices[i], + level_dof_indices[i]); } // * send diff --git a/tests/multigrid/transfer_matrix_free_06.with_mpi=true.with_p4est=true.with_trilinos=true.mpirun=15.output b/tests/multigrid/transfer_matrix_free_06.with_mpi=true.with_p4est=true.with_trilinos=true.mpirun=15.output new file mode 100644 index 0000000000..084511f9af --- /dev/null +++ b/tests/multigrid/transfer_matrix_free_06.with_mpi=true.with_p4est=true.with_trilinos=true.mpirun=15.output @@ -0,0 +1,178 @@ + +DEAL::FE: FESystem<2>[FE_Q<2>(1)^2] +DEAL::no. cells: 250 +DEAL::Diff prolongate l1: 0.00000 +DEAL::Diff prolongate l2: 0.00000 +DEAL::Diff prolongate l3: 9.61481e-17 +DEAL::Diff prolongate l4: 2.54384e-16 +DEAL::Diff prolongate l5: 2.97645e-16 +DEAL::Diff prolongate l6: 3.06884e-16 +DEAL::Diff restrict l1: 0.00000 +DEAL::Diff restrict add l1: 0.00000 +DEAL::Diff restrict l2: 0.00000 +DEAL::Diff restrict add l2: 0.00000 +DEAL::Diff restrict l3: 1.04148e-15 +DEAL::Diff restrict add l3: 8.88178e-16 +DEAL::Diff restrict l4: 9.55050e-16 +DEAL::Diff restrict add l4: 9.42055e-16 +DEAL::Diff restrict l5: 1.42178e-15 +DEAL::Diff restrict add l5: 1.60119e-15 +DEAL::Diff restrict l6: 1.33227e-15 +DEAL::Diff restrict add l6: 1.50598e-15 +DEAL:: +DEAL::no. cells: 1887 +DEAL::Diff prolongate l1: 5.55112e-17 +DEAL::Diff prolongate l2: 2.71948e-16 +DEAL::Diff prolongate l3: 5.75552e-16 +DEAL::Diff prolongate l4: 6.31096e-16 +DEAL::Diff prolongate l5: 7.90669e-16 +DEAL::Diff prolongate l6: 1.02785e-15 +DEAL::Diff restrict l1: 8.30815e-16 +DEAL::Diff restrict add l1: 8.88178e-16 +DEAL::Diff restrict l2: 1.42178e-15 +DEAL::Diff restrict add l2: 1.77636e-15 +DEAL::Diff restrict l3: 3.87307e-15 +DEAL::Diff restrict add l3: 4.02752e-15 +DEAL::Diff restrict l4: 3.21773e-15 +DEAL::Diff restrict add l4: 3.74196e-15 +DEAL::Diff restrict l5: 3.55098e-15 +DEAL::Diff restrict add l5: 4.13027e-15 +DEAL::Diff restrict l6: 4.89947e-15 +DEAL::Diff restrict add l6: 5.45706e-15 +DEAL:: +DEAL::FE: FESystem<2>[FE_Q<2>(3)^2] +DEAL::no. cells: 88 +DEAL::Diff prolongate l1: 4.25485e-16 +DEAL::Diff prolongate l2: 1.29017e-15 +DEAL::Diff prolongate l3: 2.00231e-15 +DEAL::Diff prolongate l4: 2.31946e-15 +DEAL::Diff prolongate l5: 1.85168e-15 +DEAL::Diff restrict l1: 7.69185e-16 +DEAL::Diff restrict add l1: 7.69185e-16 +DEAL::Diff restrict l2: 2.67435e-15 +DEAL::Diff restrict add l2: 2.85221e-15 +DEAL::Diff restrict l3: 4.72265e-15 +DEAL::Diff restrict add l3: 5.11185e-15 +DEAL::Diff restrict l4: 5.29601e-15 +DEAL::Diff restrict add l4: 6.06877e-15 +DEAL::Diff restrict l5: 4.22786e-15 +DEAL::Diff restrict add l5: 4.66822e-15 +DEAL:: +DEAL::no. cells: 510 +DEAL::Diff prolongate l1: 1.95148e-15 +DEAL::Diff prolongate l2: 4.10380e-15 +DEAL::Diff prolongate l3: 3.44413e-15 +DEAL::Diff prolongate l4: 4.64940e-15 +DEAL::Diff prolongate l5: 5.43550e-15 +DEAL::Diff restrict l1: 5.04518e-15 +DEAL::Diff restrict add l1: 5.79447e-15 +DEAL::Diff restrict l2: 1.14867e-14 +DEAL::Diff restrict add l2: 1.25214e-14 +DEAL::Diff restrict l3: 1.02271e-14 +DEAL::Diff restrict add l3: 1.09412e-14 +DEAL::Diff restrict l4: 1.17679e-14 +DEAL::Diff restrict add l4: 1.32526e-14 +DEAL::Diff restrict l5: 1.38716e-14 +DEAL::Diff restrict add l5: 1.55288e-14 +DEAL:: +DEAL::FE: FESystem<3>[FE_Q<3>(1)^2] +DEAL::no. cells: 15 +DEAL::Diff prolongate l1: 0.00000 +DEAL::Diff prolongate l2: 0.00000 +DEAL::Diff restrict l1: 0.00000 +DEAL::Diff restrict add l1: 0.00000 +DEAL::Diff restrict l2: 0.00000 +DEAL::Diff restrict add l2: 0.00000 +DEAL:: +DEAL::no. cells: 1791 +DEAL::Diff prolongate l1: 2.41669e-16 +DEAL::Diff prolongate l2: 3.91794e-16 +DEAL::Diff prolongate l3: 1.01512e-15 +DEAL::Diff prolongate l4: 1.55516e-15 +DEAL::Diff restrict l1: 2.51215e-15 +DEAL::Diff restrict add l1: 2.17558e-15 +DEAL::Diff restrict l2: 2.94575e-15 +DEAL::Diff restrict add l2: 2.94575e-15 +DEAL::Diff restrict l3: 6.51541e-15 +DEAL::Diff restrict add l3: 7.00761e-15 +DEAL::Diff restrict l4: 9.59492e-15 +DEAL::Diff restrict add l4: 1.07043e-14 +DEAL:: +DEAL::FE: FESystem<3>[FE_Q<3>(3)^2] +DEAL::no. cells: 1 +DEAL:: +DEAL::no. cells: 223 +DEAL::Diff prolongate l1: 6.16822e-15 +DEAL::Diff prolongate l2: 9.02624e-15 +DEAL::Diff prolongate l3: 6.05143e-15 +DEAL::Diff restrict l1: 2.49384e-14 +DEAL::Diff restrict add l1: 2.64056e-14 +DEAL::Diff restrict l2: 3.99229e-14 +DEAL::Diff restrict add l2: 4.12483e-14 +DEAL::Diff restrict l3: 3.05119e-14 +DEAL::Diff restrict add l3: 3.08762e-14 +DEAL:: +DEAL::FE: FESystem<2>[FE_Q<2>(2)^2] +DEAL::no. cells: 250 +DEAL::Diff prolongate l1: 4.16500e-08 +DEAL::Diff prolongate l2: 1.14933e-07 +DEAL::Diff prolongate l3: 4.36323e-07 +DEAL::Diff prolongate l4: 4.56555e-07 +DEAL::Diff prolongate l5: 6.44943e-07 +DEAL::Diff prolongate l6: 6.88182e-07 +DEAL::Diff restrict l1: 1.30930e-07 +DEAL::Diff restrict add l1: 1.61321e-07 +DEAL::Diff restrict l2: 3.31534e-07 +DEAL::Diff restrict add l2: 3.70787e-07 +DEAL::Diff restrict l3: 9.06752e-07 +DEAL::Diff restrict add l3: 1.06748e-06 +DEAL::Diff restrict l4: 1.12259e-06 +DEAL::Diff restrict add l4: 1.27800e-06 +DEAL::Diff restrict l5: 1.32543e-06 +DEAL::Diff restrict add l5: 1.66341e-06 +DEAL::Diff restrict l6: 1.39826e-06 +DEAL::Diff restrict add l6: 1.79262e-06 +DEAL:: +DEAL::no. cells: 1887 +DEAL::Diff prolongate l1: 3.33987e-07 +DEAL::Diff prolongate l2: 7.05824e-07 +DEAL::Diff prolongate l3: 1.48105e-06 +DEAL::Diff prolongate l4: 1.34450e-06 +DEAL::Diff prolongate l5: 1.51185e-06 +DEAL::Diff prolongate l6: 1.87422e-06 +DEAL::Diff restrict l1: 6.89860e-07 +DEAL::Diff restrict add l1: 9.08311e-07 +DEAL::Diff restrict l2: 1.47143e-06 +DEAL::Diff restrict add l2: 1.75614e-06 +DEAL::Diff restrict l3: 3.26568e-06 +DEAL::Diff restrict add l3: 3.98218e-06 +DEAL::Diff restrict l4: 2.87348e-06 +DEAL::Diff restrict add l4: 3.39202e-06 +DEAL::Diff restrict l5: 3.24632e-06 +DEAL::Diff restrict add l5: 3.95530e-06 +DEAL::Diff restrict l6: 3.95620e-06 +DEAL::Diff restrict add l6: 4.98782e-06 +DEAL:: +DEAL::FE: FESystem<3>[FE_Q<3>(2)^2] +DEAL::no. cells: 15 +DEAL::Diff prolongate l1: 1.41616e-07 +DEAL::Diff prolongate l2: 2.80670e-07 +DEAL::Diff restrict l1: 6.42811e-07 +DEAL::Diff restrict add l1: 6.42811e-07 +DEAL::Diff restrict l2: 1.11616e-06 +DEAL::Diff restrict add l2: 1.15446e-06 +DEAL:: +DEAL::no. cells: 1791 +DEAL::Diff prolongate l1: 1.19894e-06 +DEAL::Diff prolongate l2: 1.53705e-06 +DEAL::Diff prolongate l3: 2.59694e-06 +DEAL::Diff prolongate l4: 3.66785e-06 +DEAL::Diff restrict l1: 3.54901e-06 +DEAL::Diff restrict add l1: 3.80668e-06 +DEAL::Diff restrict l2: 4.41015e-06 +DEAL::Diff restrict add l2: 4.66617e-06 +DEAL::Diff restrict l3: 7.02679e-06 +DEAL::Diff restrict add l3: 7.83739e-06 +DEAL::Diff restrict l4: 9.56069e-06 +DEAL::Diff restrict add l4: 1.05378e-05 +DEAL::