From 8bf6b041fcc1d138cfa2b480ab4ba6759e89957b Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 15 Jun 2020 20:59:26 -0600 Subject: [PATCH] Add source files for the experiment. --- 9.3/taskflow-vs-tbb/2d/step-9.cc | 2142 ++++++++++++++++++++++++++++++ 9.3/taskflow-vs-tbb/3d/step-9.cc | 2142 ++++++++++++++++++++++++++++++ 2 files changed, 4284 insertions(+) create mode 100644 9.3/taskflow-vs-tbb/2d/step-9.cc create mode 100644 9.3/taskflow-vs-tbb/3d/step-9.cc diff --git a/9.3/taskflow-vs-tbb/2d/step-9.cc b/9.3/taskflow-vs-tbb/2d/step-9.cc new file mode 100644 index 0000000..64086b9 --- /dev/null +++ b/9.3/taskflow-vs-tbb/2d/step-9.cc @@ -0,0 +1,2142 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2000 - 2020 by the deal.II authors + * + * This file is part of the deal.II library. + * + * The deal.II library is free software; you can use it, redistribute + * it, and/or modify it under the terms of the GNU Lesser General + * Public License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * The full text of the license can be found in the file LICENSE.md at + * the top level directory of deal.II. + * + * --------------------------------------------------------------------- + + * + * Author: Wolfgang Bangerth, University of Heidelberg, 2000 + */ + + +// Just as in previous examples, we have to include several files of which the +// meaning has already been discussed: +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// The following two files provide classes and information for multithreaded +// programs. In the first one, the classes and functions are declared which we +// need to do assembly in parallel (i.e. the +// WorkStream namespace). The +// second file has a class MultithreadInfo which can be used to query the +// number of processors in your system, which is often useful when deciding +// how many threads to start in parallel. +#include +#include + +#include + +// The next new include file declares a base class TensorFunction +// not unlike the Function class, but with the difference that +// TensorFunction::value returns a Tensor instead of a scalar. +#include + +#include + +// This is C++, as we want to write some output to disk: +#include +#include + +#include + +// The last step is as in previous programs: +namespace Step9 +{ + using namespace dealii; + + // @sect3{Equation data declaration} + + // Next we declare a class that describes the advection field. This, of + // course, is a vector field with as many components as there are space + // dimensions. One could now use a class derived from the + // Function base class, as we have done for boundary values and + // coefficients in previous examples, but there is another possibility in + // the library, namely a base class that describes tensor valued + // functions. This is more convenient than overriding Function::value() with + // a method that knows about multiple function components: at the end of the + // day we need a Tensor, so we may as well just use a class that returns a + // Tensor. + template + class AdvectionField : public TensorFunction<1, dim> + { + public: + virtual Tensor<1, dim> value(const Point &p) const override; + + // In previous examples, we have used assertions that throw exceptions in + // several places. However, we have never seen how such exceptions are + // declared. This can be done as follows: + DeclException2(ExcDimensionMismatch, + unsigned int, + unsigned int, + << "The vector has size " << arg1 << " but should have " + << arg2 << " elements."); + // The syntax may look a little strange, but is reasonable. The format is + // basically as follows: use the name of one of the macros + // DeclExceptionN, where N denotes the number of + // additional parameters which the exception object shall take. In this + // case, as we want to throw the exception when the sizes of two vectors + // differ, we need two arguments, so we use + // DeclException2. The first parameter then describes the + // name of the exception, while the following declare the data types of + // the parameters. The last argument is a sequence of output directives + // that will be piped into the std::cerr object, thus the + // strange format with the leading @<@< operator and the + // like. Note that we can access the parameters which are passed to the + // exception upon construction (i.e. within the Assert call) + // by using the names arg1 through argN, where + // N is the number of arguments as defined by the use of the + // respective macro DeclExceptionN. + // + // To learn how the preprocessor expands this macro into actual code, + // please refer to the documentation of the exception classes. In brief, + // this macro call declares and defines a class + // ExcDimensionMismatch inheriting from ExceptionBase which + // implements all necessary error output functions. + }; + + // The following two functions implement the interface described above. The + // first simply implements the function as described in the introduction, + // while the second uses the same trick to avoid calling a virtual function + // as has already been introduced in the previous example program. Note the + // check for the right sizes of the arguments in the second function, which + // should always be present in such functions; it is our experience that + // many if not most programming errors result from incorrectly initialized + // arrays, incompatible parameters to functions and the like; using + // assertion as in this case can eliminate many of these problems. + template + Tensor<1, dim> AdvectionField::value(const Point &p) const + { + Point value; + value[0] = 2; + for (unsigned int i = 1; i < dim; ++i) + value[i] = 1 + 0.8 * std::sin(8. * numbers::PI * p[0]); + + return value; + } + + // Besides the advection field, we need two functions describing the source + // terms (right hand side) and the boundary values. As + // described in the introduction, the source is a constant function in the + // vicinity of a source point, which we denote by the constant static + // variable center_point. We set the values of this center + // using the same template tricks as we have shown in the step-7 example + // program. The rest is simple and has been shown previously. + template + class RightHandSide : public Function + { + public: + virtual double value(const Point & p, + const unsigned int component = 0) const override; + + private: + static const Point center_point; + }; + + + template <> + const Point<1> RightHandSide<1>::center_point = Point<1>(-0.75); + + template <> + const Point<2> RightHandSide<2>::center_point = Point<2>(-0.75, -0.75); + + template <> + const Point<3> RightHandSide<3>::center_point = Point<3>(-0.75, -0.75, -0.75); + + + + // The only new thing here is that we check for the value of the + // component parameter. As this is a scalar function, it is + // obvious that it only makes sense if the desired component has the index + // zero, so we assert that this is indeed the + // case. ExcIndexRange is a global predefined exception + // (probably the one most often used, we therefore made it global instead of + // local to some class), that takes three parameters: the index that is + // outside the allowed range, the first element of the valid range and the + // one past the last (i.e. again the half-open interval so often used in the + // C++ standard library): + template + double RightHandSide::value(const Point & p, + const unsigned int component) const + { + (void)component; + Assert(component == 0, ExcIndexRange(component, 0, 1)); + const double diameter = 0.1; + return ((p - center_point).norm_square() < diameter * diameter ? + 0.1 / std::pow(diameter, dim) : + 0.0); + } + + + + // Finally for the boundary values, which is just another class derived from + // the Function base class: + template + class BoundaryValues : public Function + { + public: + virtual double value(const Point & p, + const unsigned int component = 0) const override; + }; + + + + template + double BoundaryValues::value(const Point & p, + const unsigned int component) const + { + (void)component; + Assert(component == 0, ExcIndexRange(component, 0, 1)); + + const double sine_term = std::sin(16. * numbers::PI * p.norm_square()); + const double weight = std::exp(5. * (1. - p.norm_square())); + return weight * sine_term; + } + + // @sect3{AdvectionProblem class declaration} + + // Here comes the main class of this program. It is very much like the main + // classes of previous examples, so we again only comment on the + // differences. + template + class AdvectionProblem + { + public: + AdvectionProblem(); + void run(); + + private: + void setup_system(); + + // The next set of functions will be used to assemble the + // matrix. However, unlike in the previous examples, the + // assemble_system() function will not do the work + // itself, but rather will delegate the actual assembly to helper + // functions assemble_local_system() and + // copy_local_to_global(). The rationale is that + // matrix assembly can be parallelized quite well, as the + // computation of the local contributions on each cell is entirely + // independent of other cells, and we only have to synchronize + // when we add the contribution of a cell to the global + // matrix. + // + // The strategy for parallelization we choose here is one of the + // possibilities mentioned in detail in the @ref threads module in + // the documentation. Specifically, we will use the WorkStream + // approach discussed there. Since there is so much documentation + // in this module, we will not repeat the rationale for the design + // choices here (for example, if you read through the module + // mentioned above, you will understand what the purpose of the + // AssemblyScratchData and + // AssemblyCopyData structures is). Rather, we will + // only discuss the specific implementation. + // + // If you read the page mentioned above, you will find that in + // order to parallelize assembly, we need two data structures -- + // one that corresponds to data that we need during local + // integration ("scratch data", i.e., things we only need as + // temporary storage), and one that carries information from the + // local integration to the function that then adds the local + // contributions to the corresponding elements of the global + // matrix. The former of these typically contains the FEValues and + // FEFaceValues objects, whereas the latter has the local matrix, + // local right hand side, and information about which degrees of + // freedom live on the cell for which we are assembling a local + // contribution. With this information, the following should be + // relatively self-explanatory: + struct AssemblyScratchData + { + AssemblyScratchData(const FiniteElement &fe); + AssemblyScratchData(const AssemblyScratchData &scratch_data); + + // FEValues and FEFaceValues are expensive objects to set up, so we + // include them in the scratch object so that as much data is reused + // between cells as possible. + FEValues fe_values; + FEFaceValues fe_face_values; + + // We also store a few vectors that we will populate with values on each + // cell. Setting these objects up is, in the usual case, cheap; however, + // they require memory allocations, which can be expensive in + // multithreaded applications. Hence we keep them here so that + // computations on a cell do not require new allocations. + std::vector rhs_values; + std::vector> advection_directions; + std::vector face_boundary_values; + std::vector> face_advection_directions; + + // Finally, we need objects that describe the problem's data: + AdvectionField advection_field; + RightHandSide right_hand_side; + BoundaryValues boundary_values; + }; + + struct AssemblyCopyData + { + FullMatrix cell_matrix; + Vector cell_rhs; + std::vector local_dof_indices; + }; + + void assemble_system(); + void local_assemble_system( + const typename DoFHandler::active_cell_iterator &cell, + AssemblyScratchData & scratch, + AssemblyCopyData & copy_data); + void copy_local_to_global(const AssemblyCopyData ©_data); + + + // The following functions again are the same as they were in previous + // examples, as are the subsequent variables: + void solve(); + void refine_grid(); + void output_results(const unsigned int cycle) const; + + Triangulation triangulation; + DoFHandler dof_handler; + + FE_Q fe; + + AffineConstraints hanging_node_constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; + }; + + + + // @sect3{GradientEstimation class declaration} + + // Now, finally, here comes the class that will compute the difference + // approximation of the gradient on each cell and weighs that with a power + // of the mesh size, as described in the introduction. This class is a + // simple version of the DerivativeApproximation class in the + // library, that uses similar techniques to obtain finite difference + // approximations of the gradient of a finite element field, or of higher + // derivatives. + // + // The class has one public static function estimate that is + // called to compute a vector of error indicators, and a few private functions + // that do the actual work on all active cells. As in other parts of the + // library, we follow an informal convention to use vectors of floats for + // error indicators rather than the common vectors of doubles, as the + // additional accuracy is not necessary for estimated values. + // + // In addition to these two functions, the class declares two exceptions + // which are raised when a cell has no neighbors in each of the space + // directions (in which case the matrix described in the introduction would + // be singular and can't be inverted), while the other one is used in the + // more common case of invalid parameters to a function, namely a vector of + // wrong size. + // + // Two other comments: first, the class has no non-static member functions + // or variables, so this is not really a class, but rather serves the + // purpose of a namespace in C++. The reason that we chose a + // class over a namespace is that this way we can declare functions that are + // private. This can be done with namespaces as well, if one declares some + // functions in header files in the namespace and implements these and other + // functions in the implementation file. The functions not declared in the + // header file are still in the namespace but are not callable from + // outside. However, as we have only one file here, it is not possible to + // hide functions in the present case. + // + // The second comment is that the dimension template parameter is attached + // to the function rather than to the class itself. This way, you don't have + // to specify the template parameter yourself as in most other cases, but + // the compiler can figure its value out itself from the dimension of the + // DoFHandler object that one passes as first argument. + // + // Before jumping into the fray with the implementation, let us also comment + // on the parallelization strategy. We have already introduced the necessary + // framework for using the WorkStream concept in the declaration of the main + // class of this program above. We will use it again here. In the current + // context, this means that we have to define + //
    + //
  1. classes for scratch and copy objects,
    2. + //
  2. a function that does the local computation on one cell, and
    3. + //
  3. a function that copies the local result into a global object.
    4. + //
+ // Given this general framework, we will, however, deviate from it a + // bit. In particular, WorkStream was generally invented for cases where + // each local computation on a cell adds to a global object -- for + // example, when assembling linear systems where we add local contributions + // into a global matrix and right hand side. WorkStream is designed to handle + // the potential conflict of multiple threads trying to do this addition at + // the same time, and consequently has to provide for some way to ensure that + // only one thread gets to do this at a time. Here, however, the situation is + // slightly different: we compute contributions from every cell + // individually, but then all we need to do is put them into an element of + // an output vector that is unique to each cell. Consequently, there is no + // risk that the write operations from two cells might conflict, and the + // elaborate machinery of WorkStream to avoid conflicting writes is not + // necessary. Consequently, what we will do is this: We still need a scratch + // object that holds, for example, the FEValues object. However, we only + // create a fake, empty copy data structure. Likewise, we do need the + // function that computes local contributions, but since it can already put + // the result into its final location, we do not need a copy-local-to-global + // function and will instead give the WorkStream::run() function an empty + // function object -- the equivalent to a NULL function pointer. + class GradientEstimation + { + public: + template + static void estimate(const DoFHandler &dof, + const Vector & solution, + Vector & error_per_cell); + + DeclException2(ExcInvalidVectorLength, + int, + int, + << "Vector has length " << arg1 << ", but should have " + << arg2); + DeclException0(ExcInsufficientDirections); + + private: + template + struct EstimateScratchData + { + EstimateScratchData(const FiniteElement &fe, + const Vector & solution, + Vector & error_per_cell); + EstimateScratchData(const EstimateScratchData &data); + + FEValues fe_midpoint_value; + std::vector::active_cell_iterator> + active_neighbors; + + const Vector &solution; + Vector & error_per_cell; + + std::vector cell_midpoint_value; + std::vector neighbor_midpoint_value; + }; + + struct EstimateCopyData + {}; + + template + static void + estimate_cell(const typename DoFHandler::active_cell_iterator &cell, + EstimateScratchData &scratch_data, + const EstimateCopyData & copy_data); + }; + + + + // @sect3{AdvectionProblem class implementation} + + + // Now for the implementation of the main class. Constructor, destructor and + // the function setup_system follow the same pattern that was + // used previously, so we need not comment on these three function: + template + AdvectionProblem::AdvectionProblem() + : dof_handler(triangulation) + , fe(5) + {} + + + + template + void AdvectionProblem::setup_system() + { + dof_handler.distribute_dofs(fe); + hanging_node_constraints.clear(); + DoFTools::make_hanging_node_constraints(dof_handler, + hanging_node_constraints); + hanging_node_constraints.close(); + + DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, + dsp, + hanging_node_constraints, + /*keep_constrained_dofs =*/false); + sparsity_pattern.copy_from(dsp); + + system_matrix.reinit(sparsity_pattern); + + solution.reinit(dof_handler.n_dofs()); + system_rhs.reinit(dof_handler.n_dofs()); + } + + + // one task per cell, create scratch data and copydata as needed without reuse + // + namespace taskflow_v1 + { + template + void run(const Iterator & begin, + const typename identity::type &end, + Worker worker, + Copier copier, + const ScratchData & sample_scratch_data, + const CopyData & sample_copy_data, + const unsigned int queue_length = 2 * MultithreadInfo::n_threads(), + const unsigned int chunk_size = 8) + { + if (MultithreadInfo::n_threads() == 1) + { + // need to copy the sample since it is marked const + ScratchData scratch_data = sample_scratch_data; + CopyData copy_data = sample_copy_data; // NOLINT + + for (Iterator i = begin; i != end; ++i) + { + // need to check if the function is not the zero function. To + // check zero-ness, create a C++ function out of it and check that + if (static_cast &>(worker)) + worker(i, scratch_data, copy_data); + if (static_cast &>( + copier)) + copier(copy_data); + } + + return; + } + + tf::Executor &executor = MultithreadInfo::get_taskflow_executor(); + tf::Taskflow taskflow; + + ScratchData scratch_data = sample_scratch_data; + CopyData copy_data = sample_copy_data; // NOLINT + + tf::Task last_copier; + + std::vector> copy_datas; + + unsigned int idx = 0; + for (Iterator i = begin; i != end; ++i, ++idx) + { + copy_datas.emplace_back(); + + auto worker_task = taskflow + .emplace([it = i, + idx, + &sample_scratch_data, + ©_datas, + &sample_copy_data, + &worker]() { + // std::cout << "worker " << idx << std::endl; + ScratchData scratch = sample_scratch_data; + auto & copy = copy_datas[idx]; + copy = + std::make_unique(sample_copy_data); + + worker(it, scratch, *copy.get()); + }) + .name("worker"); + + tf::Task copier_task = taskflow + .emplace([idx, ©_datas, &copier]() { + copier(*copy_datas[idx].get()); + copy_datas[idx].reset(); + }) + .name("copy"); + + worker_task.precede(copier_task); + + if (!last_copier.empty()) + last_copier.precede(copier_task); + last_copier = copier_task; + } + + + // debugging: + + executor.run(taskflow).wait(); + +#ifdef DEBUG + std::cout << "done" << std::endl; + std::ofstream f("graph.dia"); + taskflow.dump(f); + f.close(); +#endif + } + + + template + void + run(const Iterator & begin, + const typename identity::type &end, + MainClass & main_object, + void (MainClass::*worker)(const Iterator &, ScratchData &, CopyData &), + void (MainClass::*copier)(const CopyData &), + const ScratchData &sample_scratch_data, + const CopyData & sample_copy_data, + const unsigned int queue_length = 2 * MultithreadInfo::n_threads(), + const unsigned int chunk_size = 8) + { + // forward to the other function + run(begin, + end, + [&main_object, worker](const Iterator &iterator, + ScratchData & scratch_data, + CopyData & copy_data) { + (main_object.*worker)(iterator, scratch_data, copy_data); + }, + [&main_object, copier](const CopyData ©_data) { + (main_object.*copier)(copy_data); + }, + sample_scratch_data, + sample_copy_data, + queue_length, + chunk_size); + } + } // namespace taskflow_v1 + + + // Like v1, except that we do not create a big task graph, but we submit 512 + // jobs at a time. + namespace taskflow_v2 + { + template + void run(const Iterator & begin, + const typename identity::type &end, + Worker worker, + Copier copier, + const ScratchData & sample_scratch_data, + const CopyData & sample_copy_data, + const unsigned int queue_length = 10 * + MultithreadInfo::n_threads(), + const unsigned int chunk_size = 8) + { + if (MultithreadInfo::n_threads() == 1) + { + // need to copy the sample since it is marked const + ScratchData scratch_data = sample_scratch_data; + CopyData copy_data = sample_copy_data; // NOLINT + + for (Iterator i = begin; i != end; ++i) + { + // need to check if the function is not the zero function. To + // check zero-ness, create a C++ function out of it and check that + if (static_cast &>(worker)) + worker(i, scratch_data, copy_data); + if (static_cast &>( + copier)) + copier(copy_data); + } + + return; + } + + const unsigned int max_work_tasks = 512; + + tf::Executor &executor = MultithreadInfo::get_taskflow_executor(); + tf::Taskflow taskflow; + + ScratchData scratch_data = sample_scratch_data; + CopyData copy_data = sample_copy_data; // NOLINT + + tf::Task last_copier; + + // Threads::ThreadLocalStorage> + // thread_local_scratch; + + std::vector> copy_datas; + + unsigned int idx = 0; + + for (Iterator i = begin; i != end; ++i, ++idx) + { + if (idx == max_work_tasks) + { + executor.run(taskflow).wait(); + taskflow.clear(); + copy_datas.clear(); + last_copier.reset(); + idx = 0; + } + + copy_datas.emplace_back(); + + auto worker_task = taskflow + .emplace([it = i, + idx, + &sample_scratch_data, + ©_datas, + &sample_copy_data, + &worker]() { + ScratchData scratch = sample_scratch_data; + auto & copy = copy_datas[idx]; + copy = + std::make_unique(sample_copy_data); + + worker(it, scratch, *copy.get()); + }) + .name("worker"); + + tf::Task copier_task = taskflow + .emplace([idx, ©_datas, &copier]() { + copier(*copy_datas[idx].get()); + copy_datas[idx].reset(); + }) + .name("copy"); + + worker_task.precede(copier_task); + + if (!last_copier.empty()) + last_copier.precede(copier_task); + last_copier = copier_task; + } + + + // debugging: + + executor.run(taskflow).wait(); + +#ifdef DEBUG + std::cout << "done" << std::endl; + std::ofstream f("graph.dia"); + taskflow.dump(f); + f.close(); +#endif + } + + + template + void + run(const Iterator & begin, + const typename identity::type &end, + MainClass & main_object, + void (MainClass::*worker)(const Iterator &, ScratchData &, CopyData &), + void (MainClass::*copier)(const CopyData &), + const ScratchData &sample_scratch_data, + const CopyData & sample_copy_data, + const unsigned int queue_length = 2 * MultithreadInfo::n_threads(), + const unsigned int chunk_size = 8) + { + // forward to the other function + run(begin, + end, + [&main_object, worker](const Iterator &iterator, + ScratchData & scratch_data, + CopyData & copy_data) { + (main_object.*worker)(iterator, scratch_data, copy_data); + }, + [&main_object, copier](const CopyData ©_data) { + (main_object.*copier)(copy_data); + }, + sample_scratch_data, + sample_copy_data, + queue_length, + chunk_size); + } + } // namespace taskflow_v2 + + // Like v1 but work in chunks of size 8 + namespace taskflow_v3 + { + template + class Chunk + { + public: + Chunk(const unsigned int count, const CopyData ©_data) + : copy_datas(count, copy_data) + {} + std::vector copy_datas; + }; + + + template + void run(const Iterator & begin, + const typename identity::type &end, + Worker worker, + Copier copier, + const ScratchData & sample_scratch_data, + const CopyData & sample_copy_data, + const unsigned int queue_length = 10 * + MultithreadInfo::n_threads(), + const unsigned int chunk_size = 8) + { + if (MultithreadInfo::n_threads() == 1) + { + // need to copy the sample since it is marked const + ScratchData scratch_data = sample_scratch_data; + CopyData copy_data = sample_copy_data; // NOLINT + + for (Iterator i = begin; i != end; ++i) + { + // need to check if the function is not the zero function. To + // check zero-ness, create a C++ function out of it and check that + if (static_cast &>(worker)) + worker(i, scratch_data, copy_data); + if (static_cast &>( + copier)) + copier(copy_data); + } + + return; + } + + tf::Executor &executor = MultithreadInfo::get_taskflow_executor(); + tf::Taskflow taskflow; + + ScratchData scratch_data = sample_scratch_data; + CopyData copy_data = sample_copy_data; // NOLINT + + tf::Task last_copier = taskflow.placeholder(); + + std::vector>> chunks; + + unsigned int idx = 0; + unsigned int remaining_items = std::distance(begin, end); + + Iterator it = begin; + while (it != end) + { + unsigned int count = std::min(remaining_items, chunk_size); + Iterator middle = it; + std::advance(middle, count); + + chunks.emplace_back(); + + // this chunk works on [it,middle) + auto worker_task = + taskflow + .emplace([it_begin = it, + it_end = middle, + idx, + count, + &sample_scratch_data, + &chunks, + &sample_copy_data, + &worker]() { + ScratchData scratch = sample_scratch_data; + chunks[idx] = + std::make_unique>(count, sample_copy_data); + + unsigned int counter = 0; + for (Iterator it = it_begin; it != it_end; ++it, ++counter) + { + worker(it, scratch, chunks[idx].get()->copy_datas[counter]); + } + }) + .name("work"); + + tf::Task copier_task = taskflow + .emplace([idx, &chunks, &copier]() mutable { + auto chunk = chunks[idx].get(); + for (auto &cd : chunk->copy_datas) + copier(cd); + + chunks[idx].reset(); + }) + .name("copy"); + + worker_task.precede(copier_task); + + last_copier.precede(copier_task); + last_copier = copier_task; + + it = middle; + ++idx; + } + // copy_datas.emplace_back(); + + // auto worker_task = taskflow + // .emplace([it = i, + // idx, + // &thread_local_scratch, + // &sample_scratch_data, + // ©_datas, + // &sample_copy_data, + // &worker]() { + // // std::cout << "worker " << idx << + // std::endl; auto &scratch_ptr = + // thread_local_scratch.get(); if + // (!scratch_ptr.get()) + // { + // thread_local_scratch = + // std::make_unique(sample_scratch_data); + // scratch_ptr = + // thread_local_scratch.get(); + // } + // auto &scratch_ptr = + // thread_local_scratch.get(); + + // auto & copy = copy_datas[idx]; + // copy = + // std::make_unique(sample_copy_data); + + // worker(it, *scratch_ptr.get(), + // *copy.get()); + // }) + // .name("worker"); + + // tf::Task copier_task = taskflow + // .emplace([idx, ©_datas, + // &copier]() { + // copier(*copy_datas[idx].get()); + // copy_datas[idx].reset(); + // }) + // .name("copy"); + + // worker_task.precede(copier_task); + + // if (!last_copier.empty()) + // last_copier.precede(copier_task); + // last_copier = copier_task; + // } + + + // debugging: + + executor.run(taskflow).wait(); + +#ifdef DEBUG + std::cout << "done" << std::endl; + std::ofstream f("graph.dia"); + taskflow.dump(f); + f.close(); +#endif + } + + + template + void + run(const Iterator & begin, + const typename identity::type &end, + MainClass & main_object, + void (MainClass::*worker)(const Iterator &, ScratchData &, CopyData &), + void (MainClass::*copier)(const CopyData &), + const ScratchData &sample_scratch_data, + const CopyData & sample_copy_data, + const unsigned int queue_length = 2 * MultithreadInfo::n_threads(), + const unsigned int chunk_size = 8) + { + // forward to the other function + run(begin, + end, + [&main_object, worker](const Iterator &iterator, + ScratchData & scratch_data, + CopyData & copy_data) { + (main_object.*worker)(iterator, scratch_data, copy_data); + }, + [&main_object, copier](const CopyData ©_data) { + (main_object.*copier)(copy_data); + }, + sample_scratch_data, + sample_copy_data, + queue_length, + chunk_size); + } + } // namespace taskflow_v3 + + // Work in chunks of size 8 (only if we have enough items, otherwise 1) and + // use a thread-local scratch object + namespace taskflow_v4 + { + template + class Chunk + { + public: + Chunk(const unsigned int count, const CopyData ©_data) + : copy_datas(count, copy_data) + {} + std::vector copy_datas; + }; + + + template + void run(const Iterator & begin, + const typename identity::type &end, + Worker worker, + Copier copier, + const ScratchData & sample_scratch_data, + const CopyData & sample_copy_data, + const unsigned int queue_length = 10 * + MultithreadInfo::n_threads(), + const unsigned int chunk_size = 8) + { + if (MultithreadInfo::n_threads() == 1) + { + // need to copy the sample since it is marked const + ScratchData scratch_data = sample_scratch_data; + CopyData copy_data = sample_copy_data; // NOLINT + + for (Iterator i = begin; i != end; ++i) + { + // need to check if the function is not the zero function. To + // check zero-ness, create a C++ function out of it and check that + if (static_cast &>(worker)) + worker(i, scratch_data, copy_data); + if (static_cast &>( + copier)) + copier(copy_data); + } + + return; + } + + tf::Executor &executor = MultithreadInfo::get_taskflow_executor(); + tf::Taskflow taskflow; + + tf::Task last_copier = taskflow.placeholder(); + + Threads::ThreadLocalStorage> + thread_local_scratch; + + std::vector>> chunks; + + unsigned int idx = 0; + unsigned int remaining_items = std::distance(begin, end); + + const unsigned int real_chunk_size = + (remaining_items / chunk_size < 3 * MultithreadInfo::n_threads()) ? + 1 : + chunk_size; + + Iterator it = begin; + while (it != end) + { + unsigned int count = std::min(remaining_items, real_chunk_size); + Iterator middle = it; + std::advance(middle, count); + + chunks.emplace_back(); + + // this chunk works on [it,middle) + auto worker_task = + taskflow + .emplace([it_begin = it, + it_end = middle, + idx, + count, + &sample_scratch_data, + &thread_local_scratch, + &chunks, + &sample_copy_data, + &worker]() { + auto &scratch_ptr = thread_local_scratch.get(); + if (!scratch_ptr.get()) + scratch_ptr = + std::make_unique(sample_scratch_data); + + ScratchData &scratch = *scratch_ptr.get(); + chunks[idx] = + std::make_unique>(count, sample_copy_data); + + unsigned int counter = 0; + for (Iterator it = it_begin; it != it_end; ++it, ++counter) + { + worker(it, scratch, chunks[idx].get()->copy_datas[counter]); + } + }) + .name("work"); + + tf::Task copier_task = taskflow + .emplace([idx, &chunks, &copier]() mutable { + auto chunk = chunks[idx].get(); + for (auto &cd : chunk->copy_datas) + copier(cd); + + chunks[idx].reset(); + }) + .name("copy"); + + worker_task.precede(copier_task); + + last_copier.precede(copier_task); + last_copier = copier_task; + + it = middle; + ++idx; + } + + // debugging: + + executor.run(taskflow).wait(); + +#ifdef DEBUG + std::cout << "done" << std::endl; + std::ofstream f("graph.dia"); + taskflow.dump(f); + f.close(); +#endif + } + + + template + void + run(const Iterator & begin, + const typename identity::type &end, + MainClass & main_object, + void (MainClass::*worker)(const Iterator &, ScratchData &, CopyData &), + void (MainClass::*copier)(const CopyData &), + const ScratchData &sample_scratch_data, + const CopyData & sample_copy_data, + const unsigned int queue_length = 2 * MultithreadInfo::n_threads(), + const unsigned int chunk_size = 8) + { + // forward to the other function + run(begin, + end, + [&main_object, worker](const Iterator &iterator, + ScratchData & scratch_data, + CopyData & copy_data) { + (main_object.*worker)(iterator, scratch_data, copy_data); + }, + [&main_object, copier](const CopyData ©_data) { + (main_object.*copier)(copy_data); + }, + sample_scratch_data, + sample_copy_data, + queue_length, + chunk_size); + } + } // namespace taskflow_v4 + + // In the following function, the matrix and right hand side are + // assembled. As stated in the documentation of the main class above, it + // does not do this itself, but rather delegates to the function following + // next, utilizing the WorkStream concept discussed in @ref threads . + // + // If you have looked through the @ref threads module, you will have + // seen that assembling in parallel does not take an incredible + // amount of extra code as long as you diligently describe what the + // scratch and copy data objects are, and if you define suitable + // functions for the local assembly and the copy operation from local + // contributions to global objects. This done, the following will do + // all the heavy lifting to get these operations done on multiple + // threads on as many cores as you have in your system: + template + void AdvectionProblem::assemble_system() + { + if (true) + taskflow_v1::run(dof_handler.begin_active(), + dof_handler.end(), + *this, + &AdvectionProblem::local_assemble_system, + &AdvectionProblem::copy_local_to_global, + AssemblyScratchData(fe), + AssemblyCopyData()); + else + WorkStream::run(dof_handler.begin_active(), + dof_handler.end(), + *this, + &AdvectionProblem::local_assemble_system, + &AdvectionProblem::copy_local_to_global, + AssemblyScratchData(fe), + AssemblyCopyData()); + } + + + + // As already mentioned above, we need to have scratch objects for + // the parallel computation of local contributions. These objects + // contain FEValues and FEFaceValues objects (as well as some arrays), and so + // we will need to have constructors and copy constructors that allow us to + // create them. For the cell terms we need the values + // and gradients of the shape functions, the quadrature points in + // order to determine the source density and the advection field at + // a given point, and the weights of the quadrature points times the + // determinant of the Jacobian at these points. In contrast, for the + // boundary integrals, we don't need the gradients, but rather the + // normal vectors to the cells. This determines which update flags + // we will have to pass to the constructors of the members of the + // class: + template + AdvectionProblem::AssemblyScratchData::AssemblyScratchData( + const FiniteElement &fe) + : fe_values(fe, + QGauss(fe.degree + 1), + update_values | update_gradients | update_quadrature_points | + update_JxW_values) + , fe_face_values(fe, + QGauss(fe.degree + 1), + update_values | update_quadrature_points | + update_JxW_values | update_normal_vectors) + , rhs_values(fe_values.get_quadrature().size()) + , advection_directions(fe_values.get_quadrature().size()) + , face_boundary_values(fe_face_values.get_quadrature().size()) + , face_advection_directions(fe_face_values.get_quadrature().size()) + {} + + + + template + AdvectionProblem::AssemblyScratchData::AssemblyScratchData( + const AssemblyScratchData &scratch_data) + : fe_values(scratch_data.fe_values.get_fe(), + scratch_data.fe_values.get_quadrature(), + update_values | update_gradients | update_quadrature_points | + update_JxW_values) + , fe_face_values(scratch_data.fe_face_values.get_fe(), + scratch_data.fe_face_values.get_quadrature(), + update_values | update_quadrature_points | + update_JxW_values | update_normal_vectors) + , rhs_values(scratch_data.rhs_values.size()) + , advection_directions(scratch_data.advection_directions.size()) + , face_boundary_values(scratch_data.face_boundary_values.size()) + , face_advection_directions(scratch_data.face_advection_directions.size()) + {} + + + + // Now, this is the function that does the actual work. It is not very + // different from the assemble_system functions of previous + // example programs, so we will again only comment on the differences. The + // mathematical stuff closely follows what we have said in the introduction. + // + // There are a number of points worth mentioning here, though. The + // first one is that we have moved the FEValues and FEFaceValues + // objects into the ScratchData object. We have done so because the + // alternative would have been to simply create one every time we + // get into this function -- i.e., on every cell. It now turns out + // that the FEValues classes were written with the explicit goal of + // moving everything that remains the same from cell to cell into + // the construction of the object, and only do as little work as + // possible in FEValues::reinit() whenever we move to a new + // cell. What this means is that it would be very expensive to + // create a new object of this kind in this function as we would + // have to do it for every cell -- exactly the thing we wanted to + // avoid with the FEValues class. Instead, what we do is create it + // only once (or a small number of times) in the scratch objects and + // then re-use it as often as we can. + // + // This begs the question of whether there are other objects we + // create in this function whose creation is expensive compared to + // its use. Indeed, at the top of the function, we declare all sorts + // of objects. The AdvectionField, + // RightHandSide and BoundaryValues do not + // cost much to create, so there is no harm here. However, + // allocating memory in creating the rhs_values and + // similar variables below typically costs a significant amount of + // time, compared to just accessing the (temporary) values we store + // in them. Consequently, these would be candidates for moving into + // the AssemblyScratchData class. We will leave this as + // an exercise. + template + void AdvectionProblem::local_assemble_system( + const typename DoFHandler::active_cell_iterator &cell, + AssemblyScratchData & scratch_data, + AssemblyCopyData & copy_data) + { + // We define some abbreviations to avoid unnecessarily long lines: + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = + scratch_data.fe_values.get_quadrature().size(); + const unsigned int n_face_q_points = + scratch_data.fe_face_values.get_quadrature().size(); + + // We declare cell matrix and cell right hand side... + copy_data.cell_matrix.reinit(dofs_per_cell, dofs_per_cell); + copy_data.cell_rhs.reinit(dofs_per_cell); + + // ... an array to hold the global indices of the degrees of freedom of + // the cell on which we are presently working... + copy_data.local_dof_indices.resize(dofs_per_cell); + + // ... then initialize the FEValues object... + scratch_data.fe_values.reinit(cell); + + // ... obtain the values of right hand side and advection directions + // at the quadrature points... + scratch_data.advection_field.value_list( + scratch_data.fe_values.get_quadrature_points(), + scratch_data.advection_directions); + scratch_data.right_hand_side.value_list( + scratch_data.fe_values.get_quadrature_points(), scratch_data.rhs_values); + + // ... set the value of the streamline diffusion parameter as + // described in the introduction... + const double delta = 0.1 * cell->diameter(); + + // ... and assemble the local contributions to the system matrix and + // right hand side as also discussed above: + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + // Alias the AssemblyScratchData object to keep the lines from + // getting too long: + const auto &sd = scratch_data; + for (unsigned int j = 0; j < dofs_per_cell; ++j) + copy_data.cell_matrix(i, j) += + ((sd.fe_values.shape_value(i, q_point) + // (phi_i + + delta * (sd.advection_directions[q_point] * // delta beta + sd.fe_values.shape_grad(i, q_point))) * // grad phi_i) + sd.advection_directions[q_point] * // beta + sd.fe_values.shape_grad(j, q_point)) * // grad phi_j + sd.fe_values.JxW(q_point); // dx + + copy_data.cell_rhs(i) += + (sd.fe_values.shape_value(i, q_point) + // (phi_i + + delta * (sd.advection_directions[q_point] * // delta beta + sd.fe_values.shape_grad(i, q_point))) * // grad phi_i) + sd.rhs_values[q_point] * // f + sd.fe_values.JxW(q_point); // dx + } + + // Besides the cell terms which we have built up now, the bilinear + // form of the present problem also contains terms on the boundary of + // the domain. Therefore, we have to check whether any of the faces of + // this cell are on the boundary of the domain, and if so assemble the + // contributions of this face as well. Of course, the bilinear form + // only contains contributions from the inflow part of + // the boundary, but to find out whether a certain part of a face of + // the present cell is part of the inflow boundary, we have to have + // information on the exact location of the quadrature points and on + // the direction of flow at this point; we obtain this information + // using the FEFaceValues object and only decide within the main loop + // whether a quadrature point is on the inflow boundary. + for (const auto &face : cell->face_iterators()) + if (face->at_boundary()) + { + // Ok, this face of the present cell is on the boundary of the + // domain. Just as for the usual FEValues object which we have + // used in previous examples and also above, we have to + // reinitialize the FEFaceValues object for the present face: + scratch_data.fe_face_values.reinit(cell, face); + + // For the quadrature points at hand, we ask for the values of + // the inflow function and for the direction of flow: + scratch_data.boundary_values.value_list( + scratch_data.fe_face_values.get_quadrature_points(), + scratch_data.face_boundary_values); + scratch_data.advection_field.value_list( + scratch_data.fe_face_values.get_quadrature_points(), + scratch_data.face_advection_directions); + + // Now loop over all quadrature points and see whether this face is on + // the inflow or outflow part of the boundary. The normal + // vector points out of the cell: since the face is at + // the boundary, the normal vector points out of the domain, + // so if the advection direction points into the domain, its + // scalar product with the normal vector must be negative (to see why + // this is true, consider the scalar product definition that uses a + // cosine): + for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point) + if (scratch_data.fe_face_values.normal_vector(q_point) * + scratch_data.face_advection_directions[q_point] < + 0.) + // If the face is part of the inflow boundary, then compute the + // contributions of this face to the global matrix and right + // hand side, using the values obtained from the + // FEFaceValues object and the formulae discussed in the + // introduction: + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + copy_data.cell_matrix(i, j) -= + (scratch_data.face_advection_directions[q_point] * + scratch_data.fe_face_values.normal_vector(q_point) * + scratch_data.fe_face_values.shape_value(i, q_point) * + scratch_data.fe_face_values.shape_value(j, q_point) * + scratch_data.fe_face_values.JxW(q_point)); + + copy_data.cell_rhs(i) -= + (scratch_data.face_advection_directions[q_point] * + scratch_data.fe_face_values.normal_vector(q_point) * + scratch_data.face_boundary_values[q_point] * + scratch_data.fe_face_values.shape_value(i, q_point) * + scratch_data.fe_face_values.JxW(q_point)); + } + } + + // The final piece of information the copy routine needs is the global + // indices of the degrees of freedom on this cell, so we end by writing + // them to the local array: + cell->get_dof_indices(copy_data.local_dof_indices); + } + + + + // The second function we needed to write was the one that copies + // the local contributions the previous function computed (and + // put into the AssemblyCopyData object) into the global matrix and right + // hand side vector objects. This is essentially what we always had + // as the last block of code when assembling something on every + // cell. The following should therefore be pretty obvious: + template + void + AdvectionProblem::copy_local_to_global(const AssemblyCopyData ©_data) + { + hanging_node_constraints.distribute_local_to_global( + copy_data.cell_matrix, + copy_data.cell_rhs, + copy_data.local_dof_indices, + system_matrix, + system_rhs); + } + + // Here comes the linear solver routine. As the system is no longer + // symmetric positive definite as in all the previous examples, we cannot + // use the Conjugate Gradient method anymore. Rather, we use a solver that + // is more general and does not rely on any special properties of the + // matrix: the GMRES method. GMRES, like the conjugate gradient method, + // requires a decent preconditioner: we use a Jacobi preconditioner here, + // which works well enough for this problem. + template + void AdvectionProblem::solve() + { + SolverControl solver_control(std::max(1000, + system_rhs.size() / 10), + 1e-10 * system_rhs.l2_norm()); + SolverGMRES> solver(solver_control); + PreconditionJacobi> preconditioner; + preconditioner.initialize(system_matrix, 1.0); + solver.solve(system_matrix, solution, system_rhs, preconditioner); + + Vector residual(dof_handler.n_dofs()); + + system_matrix.vmult(residual, solution); + residual -= system_rhs; + std::cout << " Iterations required for convergence: " + << solver_control.last_step() << '\n' + << " Max norm of residual: " + << residual.linfty_norm() << '\n'; + + hanging_node_constraints.distribute(solution); + } + + // The following function refines the grid according to the quantity + // described in the introduction. The respective computations are made in + // the class GradientEstimation. + template + void AdvectionProblem::refine_grid() + { + Vector estimated_error_per_cell(triangulation.n_active_cells()); + + GradientEstimation::estimate(dof_handler, + solution, + estimated_error_per_cell); + + GridRefinement::refine_and_coarsen_fixed_number(triangulation, + estimated_error_per_cell, + 0.3, + 0.03); + + triangulation.execute_coarsening_and_refinement(); + } + + // This function is similar to the one in step 6, but since we use a higher + // degree finite element we save the solution in a different + // way. Visualization programs like VisIt and Paraview typically only + // understand data that is associated with nodes: they cannot plot + // fifth-degree basis functions, which results in a very inaccurate picture + // of the solution we computed. To get around this we save multiple + // patches per cell: in 2D we save 64 bilinear `cells' to the VTU + // file for each cell, and in 3D we save 512. The end result is that the + // visualization program will use a piecewise linear interpolation of the + // cubic basis functions: this captures the solution detail and, with most + // screen resolutions, looks smooth. We save the grid in a separate step + // with no extra patches so that we have a visual representation of the cell + // faces. + // + // Version 9.1 of deal.II gained the ability to write higher degree + // polynomials (i.e., write piecewise bicubic visualization data for our + // piecewise bicubic solution) VTK and VTU output: however, not all recent + // versions of ParaView and VisIt (as of 2018) can read this format, so we + // use the older, more general (but less efficient) approach here. + template + void AdvectionProblem::output_results(const unsigned int cycle) const + { + { + GridOut grid_out; + std::ofstream output("grid-" + std::to_string(cycle) + ".vtu"); + grid_out.write_vtu(triangulation, output); + } + + { + DataOut data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(solution, "solution"); + data_out.build_patches(8); + + // VTU output can be expensive, both to compute and to write to + // disk. Here we ask ZLib, a compression library, to compress the data + // in a way that maximizes throughput. + DataOutBase::VtkFlags vtk_flags; + vtk_flags.compression_level = + DataOutBase::VtkFlags::ZlibCompressionLevel::best_speed; + data_out.set_flags(vtk_flags); + + std::ofstream output("solution-" + std::to_string(cycle) + ".vtu"); + data_out.write_vtu(output); + } + } + + + // ... as is the main loop (setup -- solve -- refine), aside from the number + // of cycles and the initial grid: + template + void AdvectionProblem::run() + { + for (unsigned int cycle = 0; cycle < 2; ++cycle) + { + std::cout << "Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + { + GridGenerator::hyper_cube(triangulation, -1, 1); + triangulation.refine_global(5); // 7 + } + else + { + triangulation.refine_global(2); + } + + + std::cout << " Number of active cells: " + << triangulation.n_active_cells() << std::endl; + + setup_system(); + + std::cout << " Number of degrees of freedom: " + << dof_handler.n_dofs() << std::endl; + + Timer timer; + + const unsigned int n_phys_cores = MultithreadInfo::n_cores(); + std::cout << "MultithreadInfo::n_cores()=" << n_phys_cores << std::endl; + std::cout << "MultithreadInfo::n_threads()=" + << MultithreadInfo::n_threads() << std::endl; + + const unsigned int runs = 5; + +#ifdef DEAL_II_WITH_CPP_TASKFLOW + std::cout << "** TASKFLOW v1 **" << std::endl; + + + MultithreadInfo::set_thread_limit(); + + for (unsigned int n_cores = 2 * n_phys_cores; n_cores > 0; n_cores /= 2) + { + if (n_cores <= n_phys_cores) + MultithreadInfo::set_thread_limit(n_cores); + + std::cout << "n_cores "; + if (n_cores <= n_phys_cores) + std::cout << n_cores; + else + std::cout << "auto"; + std::cout << ' ' << std::flush; + double avg = 0.; + + for (unsigned int c = 0; c < runs; ++c) + { + timer.reset(); + timer.start(); + + taskflow_v1::run(dof_handler.begin_active(), + dof_handler.end(), + *this, + &AdvectionProblem::local_assemble_system, + &AdvectionProblem::copy_local_to_global, + AssemblyScratchData(fe), + AssemblyCopyData()); + + + timer.stop(); + const double time = timer.last_wall_time(); + avg += time; + std::cout << time << " " << std::flush; + } + avg /= runs; + std::cout << " avg: " << avg << std::endl; + } + + std::cout << "** TASKFLOW v2 **" << std::endl; + + + MultithreadInfo::set_thread_limit(); + + for (unsigned int n_cores = 2 * n_phys_cores; n_cores > 0; n_cores /= 2) + { + if (n_cores <= n_phys_cores) + MultithreadInfo::set_thread_limit(n_cores); + + std::cout << "n_cores "; + if (n_cores <= n_phys_cores) + std::cout << n_cores; + else + std::cout << "auto"; + std::cout << ' ' << std::flush; + double avg = 0.; + + for (unsigned int c = 0; c < runs; ++c) + { + timer.reset(); + timer.start(); + + taskflow_v2::run(dof_handler.begin_active(), + dof_handler.end(), + *this, + &AdvectionProblem::local_assemble_system, + &AdvectionProblem::copy_local_to_global, + AssemblyScratchData(fe), + AssemblyCopyData()); + + + timer.stop(); + const double time = timer.last_wall_time(); + avg += time; + std::cout << time << " " << std::flush; + } + avg /= runs; + std::cout << " avg: " << avg << std::endl; + } + + + std::cout << "** TASKFLOW v3 **" << std::endl; + + + MultithreadInfo::set_thread_limit(); + + for (unsigned int n_cores = 2 * n_phys_cores; n_cores > 0; n_cores /= 2) + { + if (n_cores <= n_phys_cores) + MultithreadInfo::set_thread_limit(n_cores); + + std::cout << "n_cores "; + if (n_cores <= n_phys_cores) + std::cout << n_cores; + else + std::cout << "auto"; + std::cout << ' ' << std::flush; + double avg = 0.; + + for (unsigned int c = 0; c < runs; ++c) + { + timer.reset(); + timer.start(); + + taskflow_v3::run(dof_handler.begin_active(), + dof_handler.end(), + *this, + &AdvectionProblem::local_assemble_system, + &AdvectionProblem::copy_local_to_global, + AssemblyScratchData(fe), + AssemblyCopyData()); + + + timer.stop(); + const double time = timer.last_wall_time(); + avg += time; + std::cout << time << " " << std::flush; + } + avg /= runs; + std::cout << " avg: " << avg << std::endl; + } + + + + std::cout << "** TASKFLOW v4 **" << std::endl; + + + MultithreadInfo::set_thread_limit(); + + for (unsigned int n_cores = 2 * n_phys_cores; n_cores > 0; n_cores /= 2) + { + if (n_cores <= n_phys_cores) + MultithreadInfo::set_thread_limit(n_cores); + + std::cout << "n_cores "; + if (n_cores <= n_phys_cores) + std::cout << n_cores; + else + std::cout << "auto"; + std::cout << ' ' << std::flush; + double avg = 0.; + + for (unsigned int c = 0; c < runs; ++c) + { + timer.reset(); + timer.start(); + + taskflow_v4::run(dof_handler.begin_active(), + dof_handler.end(), + *this, + &AdvectionProblem::local_assemble_system, + &AdvectionProblem::copy_local_to_global, + AssemblyScratchData(fe), + AssemblyCopyData()); + + + timer.stop(); + const double time = timer.last_wall_time(); + avg += time; + std::cout << time << " " << std::flush; + } + avg /= runs; + std::cout << " avg: " << avg << std::endl; + } +#endif +#ifdef DEAL_II_WITH_TBB + std::cout << "** TBB **" << std::endl; + + + MultithreadInfo::set_thread_limit(); + + for (unsigned int n_cores = 2 * n_phys_cores; n_cores > 0; n_cores /= 2) + { + if (n_cores <= n_phys_cores) + MultithreadInfo::set_thread_limit(n_cores); + + std::cout << "n_cores "; + if (n_cores <= n_phys_cores) + std::cout << n_cores; + else + std::cout << "auto"; + std::cout << ' ' << std::flush; + double avg = 0.; + + for (unsigned int c = 0; c < runs; ++c) + { + timer.reset(); + timer.start(); + WorkStream::run(dof_handler.begin_active(), + dof_handler.end(), + *this, + &AdvectionProblem::local_assemble_system, + &AdvectionProblem::copy_local_to_global, + AssemblyScratchData(fe), + AssemblyCopyData()); + + + timer.stop(); + const double time = timer.last_wall_time(); + avg += time; + std::cout << time << " " << std::flush; + } + avg /= runs; + std::cout << " avg: " << avg << std::endl; + } +#endif + + // assemble_system(); + + + // solve(); + // output_results(cycle); + } + } + + + + // @sect3{GradientEstimation class implementation} + + // Now for the implementation of the GradientEstimation class. + // Let us start by defining constructors for the + // EstimateScratchData class used by the + // estimate_cell() function: + template + GradientEstimation::EstimateScratchData::EstimateScratchData( + const FiniteElement &fe, + const Vector & solution, + Vector & error_per_cell) + : fe_midpoint_value(fe, + QMidpoint(), + update_values | update_quadrature_points) + , solution(solution) + , error_per_cell(error_per_cell) + , cell_midpoint_value(1) + , neighbor_midpoint_value(1) + { + // We allocate a vector to hold iterators to all active neighbors of + // a cell. We reserve the maximal number of active neighbors in order to + // avoid later reallocations. Note how this maximal number of active + // neighbors is computed here. + active_neighbors.reserve(GeometryInfo::faces_per_cell * + GeometryInfo::max_children_per_face); + } + + + template + GradientEstimation::EstimateScratchData::EstimateScratchData( + const EstimateScratchData &scratch_data) + : fe_midpoint_value(scratch_data.fe_midpoint_value.get_fe(), + scratch_data.fe_midpoint_value.get_quadrature(), + update_values | update_quadrature_points) + , solution(scratch_data.solution) + , error_per_cell(scratch_data.error_per_cell) + , cell_midpoint_value(1) + , neighbor_midpoint_value(1) + {} + + + // Next comes the implementation of the GradientEstimation + // class. The first function does not much except for delegating work to the + // other function, but there is a bit of setup at the top. + // + // Before starting with the work, we check that the vector into + // which the results are written has the right size. Programming + // mistakes in which one forgets to size arguments correctly at the + // calling site are quite common. Because the resulting damage from + // not catching such errors is often subtle (e.g., corruption of + // data somewhere in memory, or non-reproducible results), it is + // well worth the effort to check for such things. + template + void GradientEstimation::estimate(const DoFHandler &dof_handler, + const Vector & solution, + Vector & error_per_cell) + { + Assert( + error_per_cell.size() == dof_handler.get_triangulation().n_active_cells(), + ExcInvalidVectorLength(error_per_cell.size(), + dof_handler.get_triangulation().n_active_cells())); + + WorkStream::run(dof_handler.begin_active(), + dof_handler.end(), + &GradientEstimation::template estimate_cell, + std::function(), + EstimateScratchData(dof_handler.get_fe(), + solution, + error_per_cell), + EstimateCopyData()); + } + + + // Here comes the function that estimates the local error by computing the + // finite difference approximation of the gradient. The function first + // computes the list of active neighbors of the present cell and then + // computes the quantities described in the introduction for each of + // the neighbors. The reason for this order is that it is not a one-liner + // to find a given neighbor with locally refined meshes. In principle, an + // optimized implementation would find neighbors and the quantities + // depending on them in one step, rather than first building a list of + // neighbors and in a second step their contributions but we will gladly + // leave this as an exercise. As discussed before, the worker function + // passed to WorkStream::run works on "scratch" objects that keep all + // temporary objects. This way, we do not need to create and initialize + // objects that are expensive to initialize within the function that does + // the work every time it is called for a given cell. Such an argument is + // passed as the second argument. The third argument would be a "copy-data" + // object (see @ref threads for more information) but we do not actually use + // any of these here. Since WorkStream::run() insists on passing three + // arguments, we declare this function with three arguments, but simply + // ignore the last one. + // + // (This is unsatisfactory from an aesthetic perspective. It can be avoided + // by using an anonymous (lambda) function. If you allow, let us here show + // how. First, assume that we had declared this function to only take two + // arguments by omitting the unused last one. Now, WorkStream::run still + // wants to call this function with three arguments, so we need to find a + // way to "forget" the third argument in the call. Simply passing + // WorkStream::run the pointer to the function as we do above will not do + // this -- the compiler will complain that a function declared to have two + // arguments is called with three arguments. However, we can do this by + // passing the following as the third argument to WorkStream::run(): + // @code + // [](const typename DoFHandler::active_cell_iterator &cell, + // EstimateScratchData & scratch_data, + // EstimateCopyData &) + // { + // GradientEstimation::estimate_cell(cell, scratch_data); + // } + // @endcode + // This is not much better than the solution implemented below: either the + // routine itself must take three arguments or it must be wrapped by + // something that takes three arguments. We don't use this since adding the + // unused argument at the beginning is simpler. + // + // Now for the details: + template + void GradientEstimation::estimate_cell( + const typename DoFHandler::active_cell_iterator &cell, + EstimateScratchData & scratch_data, + const EstimateCopyData &) + { + // We need space for the tensor Y, which is the sum of + // outer products of the y-vectors. + Tensor<2, dim> Y; + + // First initialize the FEValues object, as well as the + // Y tensor: + scratch_data.fe_midpoint_value.reinit(cell); + + // Now, before we go on, we first compute a list of all active neighbors + // of the present cell. We do so by first looping over all faces and see + // whether the neighbor there is active, which would be the case if it + // is on the same level as the present cell or one level coarser (note + // that a neighbor can only be once coarser than the present cell, as + // we only allow a maximal difference of one refinement over a face in + // deal.II). Alternatively, the neighbor could be on the same level + // and be further refined; then we have to find which of its children + // are next to the present cell and select these (note that if a child + // of a neighbor of an active cell that is next to this active cell, + // needs necessarily be active itself, due to the one-refinement rule + // cited above). + // + // Things are slightly different in one space dimension, as there the + // one-refinement rule does not exist: neighboring active cells may + // differ in as many refinement levels as they like. In this case, the + // computation becomes a little more difficult, but we will explain + // this below. + // + // Before starting the loop over all neighbors of the present cell, we + // have to clear the array storing the iterators to the active + // neighbors, of course. + scratch_data.active_neighbors.clear(); + for (unsigned int face_n : GeometryInfo::face_indices()) + if (!cell->at_boundary(face_n)) + { + // First define an abbreviation for the iterator to the face and + // the neighbor + const auto face = cell->face(face_n); + const auto neighbor = cell->neighbor(face_n); + + // Then check whether the neighbor is active. If it is, then it + // is on the same level or one level coarser (if we are not in + // 1D), and we are interested in it in any case. + if (neighbor->is_active()) + scratch_data.active_neighbors.push_back(neighbor); + else + { + // If the neighbor is not active, then check its children. + if (dim == 1) + { + // To find the child of the neighbor which bounds to the + // present cell, successively go to its right child if + // we are left of the present cell (n==0), or go to the + // left child if we are on the right (n==1), until we + // find an active cell. + auto neighbor_child = neighbor; + while (neighbor_child->has_children()) + neighbor_child = neighbor_child->child(face_n == 0 ? 1 : 0); + + // As this used some non-trivial geometrical intuition, + // we might want to check whether we did it right, + // i.e., check whether the neighbor of the cell we found + // is indeed the cell we are presently working + // on. Checks like this are often useful and have + // frequently uncovered errors both in algorithms like + // the line above (where it is simple to involuntarily + // exchange n==1 for n==0 or + // the like) and in the library (the assumptions + // underlying the algorithm above could either be wrong, + // wrongly documented, or are violated due to an error + // in the library). One could in principle remove such + // checks after the program works for some time, but it + // might be a good things to leave it in anyway to check + // for changes in the library or in the algorithm above. + // + // Note that if this check fails, then this is certainly + // an error that is irrecoverable and probably qualifies + // as an internal error. We therefore use a predefined + // exception class to throw here. + Assert(neighbor_child->neighbor(face_n == 0 ? 1 : 0) == cell, + ExcInternalError()); + + // If the check succeeded, we push the active neighbor + // we just found to the stack we keep: + scratch_data.active_neighbors.push_back(neighbor_child); + } + else + // If we are not in 1d, we collect all neighbor children + // `behind' the subfaces of the current face and move on: + for (unsigned int subface_n = 0; subface_n < face->n_children(); + ++subface_n) + scratch_data.active_neighbors.push_back( + cell->neighbor_child_on_subface(face_n, subface_n)); + } + } + + // OK, now that we have all the neighbors, lets start the computation + // on each of them. First we do some preliminaries: find out about the + // center of the present cell and the solution at this point. The + // latter is obtained as a vector of function values at the quadrature + // points, of which there are only one, of course. Likewise, the + // position of the center is the position of the first (and only) + // quadrature point in real space. + const Point this_center = + scratch_data.fe_midpoint_value.quadrature_point(0); + + scratch_data.fe_midpoint_value.get_function_values( + scratch_data.solution, scratch_data.cell_midpoint_value); + + // Now loop over all active neighbors and collect the data we + // need. + Tensor<1, dim> projected_gradient; + for (const auto &neighbor : scratch_data.active_neighbors) + { + // Then get the center of the neighbor cell and the value of the + // finite element function at that point. Note that for this + // information we have to reinitialize the FEValues + // object for the neighbor cell. + scratch_data.fe_midpoint_value.reinit(neighbor); + const Point neighbor_center = + scratch_data.fe_midpoint_value.quadrature_point(0); + + scratch_data.fe_midpoint_value.get_function_values( + scratch_data.solution, scratch_data.neighbor_midpoint_value); + + // Compute the vector y connecting the centers of the + // two cells. Note that as opposed to the introduction, we denote + // by y the normalized difference vector, as this is + // the quantity used everywhere in the computations. + Tensor<1, dim> y = neighbor_center - this_center; + const double distance = y.norm(); + y /= distance; + + // Then add up the contribution of this cell to the Y matrix... + for (unsigned int i = 0; i < dim; ++i) + for (unsigned int j = 0; j < dim; ++j) + Y[i][j] += y[i] * y[j]; + + // ... and update the sum of difference quotients: + projected_gradient += (scratch_data.neighbor_midpoint_value[0] - + scratch_data.cell_midpoint_value[0]) / + distance * y; + } + + // If now, after collecting all the information from the neighbors, we + // can determine an approximation of the gradient for the present + // cell, then we need to have passed over vectors y which + // span the whole space, otherwise we would not have all components of + // the gradient. This is indicated by the invertibility of the matrix. + // + // If the matrix is not invertible, then the present + // cell had an insufficient number of active neighbors. In contrast to + // all previous cases (where we raised exceptions) this is, however, + // not a programming error: it is a runtime error that can happen in + // optimized mode even if it ran well in debug mode, so it is + // reasonable to try to catch this error also in optimized mode. For + // this case, there is the AssertThrow macro: it checks + // the condition like the Assert macro, but not only in + // debug mode; it then outputs an error message, but instead of + // aborting the program as in the case of the Assert + // macro, the exception is thrown using the throw command + // of C++. This way, one has the possibility to catch this error and + // take reasonable counter actions. One such measure would be to + // refine the grid globally, as the case of insufficient directions + // can not occur if every cell of the initial grid has been refined at + // least once. + AssertThrow(determinant(Y) != 0, ExcInsufficientDirections()); + + // If, on the other hand, the matrix is invertible, then invert it, + // multiply the other quantity with it, and compute the estimated error + // using this quantity and the correct powers of the mesh width: + const Tensor<2, dim> Y_inverse = invert(Y); + + const Tensor<1, dim> gradient = Y_inverse * projected_gradient; + + // The last part of this function is the one where we write into + // the element of the output vector what we have just + // computed. The address of this vector has been stored in the + // scratch data object, and all we have to do is know how to get + // at the correct element inside this vector -- but we can ask the + // cell we're on the how-manyth active cell it is for this: + scratch_data.error_per_cell(cell->active_cell_index()) = + (std::pow(cell->diameter(), 1 + 1.0 * dim / 2) * gradient.norm()); + } +} // namespace Step9 + + +// @sect3{Main function} + +// The main function is similar to the previous examples. The +// primary difference is that we use MultithreadInfo to set the maximum +// number of threads (see the documentation module @ref threads +// "Parallel computing with multiple processors accessing shared memory" +// for more information). The number of threads used is the minimum of the +// environment variable DEAL_II_NUM_THREADS and the parameter of +// set_thread_limit. If no value is given to +// set_thread_limit, the default value from the Intel Threading +// Building Blocks (TBB) library is used. If the call to +// set_thread_limit is omitted, the number of threads will be +// chosen by TBB independently of DEAL_II_NUM_THREADS. +int main() +{ + using namespace dealii; + try + { + // MultithreadInfo::set_thread_limit(); + + Step9::AdvectionProblem<2> advection_problem_2d; + advection_problem_2d.run(); + } + catch (std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} diff --git a/9.3/taskflow-vs-tbb/3d/step-9.cc b/9.3/taskflow-vs-tbb/3d/step-9.cc new file mode 100644 index 0000000..0cebbae --- /dev/null +++ b/9.3/taskflow-vs-tbb/3d/step-9.cc @@ -0,0 +1,2142 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2000 - 2020 by the deal.II authors + * + * This file is part of the deal.II library. + * + * The deal.II library is free software; you can use it, redistribute + * it, and/or modify it under the terms of the GNU Lesser General + * Public License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * The full text of the license can be found in the file LICENSE.md at + * the top level directory of deal.II. + * + * --------------------------------------------------------------------- + + * + * Author: Wolfgang Bangerth, University of Heidelberg, 2000 + */ + + +// Just as in previous examples, we have to include several files of which the +// meaning has already been discussed: +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// The following two files provide classes and information for multithreaded +// programs. In the first one, the classes and functions are declared which we +// need to do assembly in parallel (i.e. the +// WorkStream namespace). The +// second file has a class MultithreadInfo which can be used to query the +// number of processors in your system, which is often useful when deciding +// how many threads to start in parallel. +#include +#include + +#include + +// The next new include file declares a base class TensorFunction +// not unlike the Function class, but with the difference that +// TensorFunction::value returns a Tensor instead of a scalar. +#include + +#include + +// This is C++, as we want to write some output to disk: +#include +#include + +#include + +// The last step is as in previous programs: +namespace Step9 +{ + using namespace dealii; + + // @sect3{Equation data declaration} + + // Next we declare a class that describes the advection field. This, of + // course, is a vector field with as many components as there are space + // dimensions. One could now use a class derived from the + // Function base class, as we have done for boundary values and + // coefficients in previous examples, but there is another possibility in + // the library, namely a base class that describes tensor valued + // functions. This is more convenient than overriding Function::value() with + // a method that knows about multiple function components: at the end of the + // day we need a Tensor, so we may as well just use a class that returns a + // Tensor. + template + class AdvectionField : public TensorFunction<1, dim> + { + public: + virtual Tensor<1, dim> value(const Point &p) const override; + + // In previous examples, we have used assertions that throw exceptions in + // several places. However, we have never seen how such exceptions are + // declared. This can be done as follows: + DeclException2(ExcDimensionMismatch, + unsigned int, + unsigned int, + << "The vector has size " << arg1 << " but should have " + << arg2 << " elements."); + // The syntax may look a little strange, but is reasonable. The format is + // basically as follows: use the name of one of the macros + // DeclExceptionN, where N denotes the number of + // additional parameters which the exception object shall take. In this + // case, as we want to throw the exception when the sizes of two vectors + // differ, we need two arguments, so we use + // DeclException2. The first parameter then describes the + // name of the exception, while the following declare the data types of + // the parameters. The last argument is a sequence of output directives + // that will be piped into the std::cerr object, thus the + // strange format with the leading @<@< operator and the + // like. Note that we can access the parameters which are passed to the + // exception upon construction (i.e. within the Assert call) + // by using the names arg1 through argN, where + // N is the number of arguments as defined by the use of the + // respective macro DeclExceptionN. + // + // To learn how the preprocessor expands this macro into actual code, + // please refer to the documentation of the exception classes. In brief, + // this macro call declares and defines a class + // ExcDimensionMismatch inheriting from ExceptionBase which + // implements all necessary error output functions. + }; + + // The following two functions implement the interface described above. The + // first simply implements the function as described in the introduction, + // while the second uses the same trick to avoid calling a virtual function + // as has already been introduced in the previous example program. Note the + // check for the right sizes of the arguments in the second function, which + // should always be present in such functions; it is our experience that + // many if not most programming errors result from incorrectly initialized + // arrays, incompatible parameters to functions and the like; using + // assertion as in this case can eliminate many of these problems. + template + Tensor<1, dim> AdvectionField::value(const Point &p) const + { + Point value; + value[0] = 2; + for (unsigned int i = 1; i < dim; ++i) + value[i] = 1 + 0.8 * std::sin(8. * numbers::PI * p[0]); + + return value; + } + + // Besides the advection field, we need two functions describing the source + // terms (right hand side) and the boundary values. As + // described in the introduction, the source is a constant function in the + // vicinity of a source point, which we denote by the constant static + // variable center_point. We set the values of this center + // using the same template tricks as we have shown in the step-7 example + // program. The rest is simple and has been shown previously. + template + class RightHandSide : public Function + { + public: + virtual double value(const Point & p, + const unsigned int component = 0) const override; + + private: + static const Point center_point; + }; + + + template <> + const Point<1> RightHandSide<1>::center_point = Point<1>(-0.75); + + template <> + const Point<2> RightHandSide<2>::center_point = Point<2>(-0.75, -0.75); + + template <> + const Point<3> RightHandSide<3>::center_point = Point<3>(-0.75, -0.75, -0.75); + + + + // The only new thing here is that we check for the value of the + // component parameter. As this is a scalar function, it is + // obvious that it only makes sense if the desired component has the index + // zero, so we assert that this is indeed the + // case. ExcIndexRange is a global predefined exception + // (probably the one most often used, we therefore made it global instead of + // local to some class), that takes three parameters: the index that is + // outside the allowed range, the first element of the valid range and the + // one past the last (i.e. again the half-open interval so often used in the + // C++ standard library): + template + double RightHandSide::value(const Point & p, + const unsigned int component) const + { + (void)component; + Assert(component == 0, ExcIndexRange(component, 0, 1)); + const double diameter = 0.1; + return ((p - center_point).norm_square() < diameter * diameter ? + 0.1 / std::pow(diameter, dim) : + 0.0); + } + + + + // Finally for the boundary values, which is just another class derived from + // the Function base class: + template + class BoundaryValues : public Function + { + public: + virtual double value(const Point & p, + const unsigned int component = 0) const override; + }; + + + + template + double BoundaryValues::value(const Point & p, + const unsigned int component) const + { + (void)component; + Assert(component == 0, ExcIndexRange(component, 0, 1)); + + const double sine_term = std::sin(16. * numbers::PI * p.norm_square()); + const double weight = std::exp(5. * (1. - p.norm_square())); + return weight * sine_term; + } + + // @sect3{AdvectionProblem class declaration} + + // Here comes the main class of this program. It is very much like the main + // classes of previous examples, so we again only comment on the + // differences. + template + class AdvectionProblem + { + public: + AdvectionProblem(); + void run(); + + private: + void setup_system(); + + // The next set of functions will be used to assemble the + // matrix. However, unlike in the previous examples, the + // assemble_system() function will not do the work + // itself, but rather will delegate the actual assembly to helper + // functions assemble_local_system() and + // copy_local_to_global(). The rationale is that + // matrix assembly can be parallelized quite well, as the + // computation of the local contributions on each cell is entirely + // independent of other cells, and we only have to synchronize + // when we add the contribution of a cell to the global + // matrix. + // + // The strategy for parallelization we choose here is one of the + // possibilities mentioned in detail in the @ref threads module in + // the documentation. Specifically, we will use the WorkStream + // approach discussed there. Since there is so much documentation + // in this module, we will not repeat the rationale for the design + // choices here (for example, if you read through the module + // mentioned above, you will understand what the purpose of the + // AssemblyScratchData and + // AssemblyCopyData structures is). Rather, we will + // only discuss the specific implementation. + // + // If you read the page mentioned above, you will find that in + // order to parallelize assembly, we need two data structures -- + // one that corresponds to data that we need during local + // integration ("scratch data", i.e., things we only need as + // temporary storage), and one that carries information from the + // local integration to the function that then adds the local + // contributions to the corresponding elements of the global + // matrix. The former of these typically contains the FEValues and + // FEFaceValues objects, whereas the latter has the local matrix, + // local right hand side, and information about which degrees of + // freedom live on the cell for which we are assembling a local + // contribution. With this information, the following should be + // relatively self-explanatory: + struct AssemblyScratchData + { + AssemblyScratchData(const FiniteElement &fe); + AssemblyScratchData(const AssemblyScratchData &scratch_data); + + // FEValues and FEFaceValues are expensive objects to set up, so we + // include them in the scratch object so that as much data is reused + // between cells as possible. + FEValues fe_values; + FEFaceValues fe_face_values; + + // We also store a few vectors that we will populate with values on each + // cell. Setting these objects up is, in the usual case, cheap; however, + // they require memory allocations, which can be expensive in + // multithreaded applications. Hence we keep them here so that + // computations on a cell do not require new allocations. + std::vector rhs_values; + std::vector> advection_directions; + std::vector face_boundary_values; + std::vector> face_advection_directions; + + // Finally, we need objects that describe the problem's data: + AdvectionField advection_field; + RightHandSide right_hand_side; + BoundaryValues boundary_values; + }; + + struct AssemblyCopyData + { + FullMatrix cell_matrix; + Vector cell_rhs; + std::vector local_dof_indices; + }; + + void assemble_system(); + void local_assemble_system( + const typename DoFHandler::active_cell_iterator &cell, + AssemblyScratchData & scratch, + AssemblyCopyData & copy_data); + void copy_local_to_global(const AssemblyCopyData ©_data); + + + // The following functions again are the same as they were in previous + // examples, as are the subsequent variables: + void solve(); + void refine_grid(); + void output_results(const unsigned int cycle) const; + + Triangulation triangulation; + DoFHandler dof_handler; + + FE_Q fe; + + AffineConstraints hanging_node_constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; + }; + + + + // @sect3{GradientEstimation class declaration} + + // Now, finally, here comes the class that will compute the difference + // approximation of the gradient on each cell and weighs that with a power + // of the mesh size, as described in the introduction. This class is a + // simple version of the DerivativeApproximation class in the + // library, that uses similar techniques to obtain finite difference + // approximations of the gradient of a finite element field, or of higher + // derivatives. + // + // The class has one public static function estimate that is + // called to compute a vector of error indicators, and a few private functions + // that do the actual work on all active cells. As in other parts of the + // library, we follow an informal convention to use vectors of floats for + // error indicators rather than the common vectors of doubles, as the + // additional accuracy is not necessary for estimated values. + // + // In addition to these two functions, the class declares two exceptions + // which are raised when a cell has no neighbors in each of the space + // directions (in which case the matrix described in the introduction would + // be singular and can't be inverted), while the other one is used in the + // more common case of invalid parameters to a function, namely a vector of + // wrong size. + // + // Two other comments: first, the class has no non-static member functions + // or variables, so this is not really a class, but rather serves the + // purpose of a namespace in C++. The reason that we chose a + // class over a namespace is that this way we can declare functions that are + // private. This can be done with namespaces as well, if one declares some + // functions in header files in the namespace and implements these and other + // functions in the implementation file. The functions not declared in the + // header file are still in the namespace but are not callable from + // outside. However, as we have only one file here, it is not possible to + // hide functions in the present case. + // + // The second comment is that the dimension template parameter is attached + // to the function rather than to the class itself. This way, you don't have + // to specify the template parameter yourself as in most other cases, but + // the compiler can figure its value out itself from the dimension of the + // DoFHandler object that one passes as first argument. + // + // Before jumping into the fray with the implementation, let us also comment + // on the parallelization strategy. We have already introduced the necessary + // framework for using the WorkStream concept in the declaration of the main + // class of this program above. We will use it again here. In the current + // context, this means that we have to define + //
    + //
  1. classes for scratch and copy objects,
    2. + //
  2. a function that does the local computation on one cell, and
    3. + //
  3. a function that copies the local result into a global object.
    4. + //
+ // Given this general framework, we will, however, deviate from it a + // bit. In particular, WorkStream was generally invented for cases where + // each local computation on a cell adds to a global object -- for + // example, when assembling linear systems where we add local contributions + // into a global matrix and right hand side. WorkStream is designed to handle + // the potential conflict of multiple threads trying to do this addition at + // the same time, and consequently has to provide for some way to ensure that + // only one thread gets to do this at a time. Here, however, the situation is + // slightly different: we compute contributions from every cell + // individually, but then all we need to do is put them into an element of + // an output vector that is unique to each cell. Consequently, there is no + // risk that the write operations from two cells might conflict, and the + // elaborate machinery of WorkStream to avoid conflicting writes is not + // necessary. Consequently, what we will do is this: We still need a scratch + // object that holds, for example, the FEValues object. However, we only + // create a fake, empty copy data structure. Likewise, we do need the + // function that computes local contributions, but since it can already put + // the result into its final location, we do not need a copy-local-to-global + // function and will instead give the WorkStream::run() function an empty + // function object -- the equivalent to a NULL function pointer. + class GradientEstimation + { + public: + template + static void estimate(const DoFHandler &dof, + const Vector & solution, + Vector & error_per_cell); + + DeclException2(ExcInvalidVectorLength, + int, + int, + << "Vector has length " << arg1 << ", but should have " + << arg2); + DeclException0(ExcInsufficientDirections); + + private: + template + struct EstimateScratchData + { + EstimateScratchData(const FiniteElement &fe, + const Vector & solution, + Vector & error_per_cell); + EstimateScratchData(const EstimateScratchData &data); + + FEValues fe_midpoint_value; + std::vector::active_cell_iterator> + active_neighbors; + + const Vector &solution; + Vector & error_per_cell; + + std::vector cell_midpoint_value; + std::vector neighbor_midpoint_value; + }; + + struct EstimateCopyData + {}; + + template + static void + estimate_cell(const typename DoFHandler::active_cell_iterator &cell, + EstimateScratchData &scratch_data, + const EstimateCopyData & copy_data); + }; + + + + // @sect3{AdvectionProblem class implementation} + + + // Now for the implementation of the main class. Constructor, destructor and + // the function setup_system follow the same pattern that was + // used previously, so we need not comment on these three function: + template + AdvectionProblem::AdvectionProblem() + : dof_handler(triangulation) + , fe(5) + {} + + + + template + void AdvectionProblem::setup_system() + { + dof_handler.distribute_dofs(fe); + hanging_node_constraints.clear(); + DoFTools::make_hanging_node_constraints(dof_handler, + hanging_node_constraints); + hanging_node_constraints.close(); + + DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, + dsp, + hanging_node_constraints, + /*keep_constrained_dofs =*/false); + sparsity_pattern.copy_from(dsp); + + system_matrix.reinit(sparsity_pattern); + + solution.reinit(dof_handler.n_dofs()); + system_rhs.reinit(dof_handler.n_dofs()); + } + + + // one task per cell, create scratch data and copydata as needed without reuse + // + namespace taskflow_v1 + { + template + void run(const Iterator & begin, + const typename identity::type &end, + Worker worker, + Copier copier, + const ScratchData & sample_scratch_data, + const CopyData & sample_copy_data, + const unsigned int queue_length = 2 * MultithreadInfo::n_threads(), + const unsigned int chunk_size = 8) + { + if (MultithreadInfo::n_threads() == 1) + { + // need to copy the sample since it is marked const + ScratchData scratch_data = sample_scratch_data; + CopyData copy_data = sample_copy_data; // NOLINT + + for (Iterator i = begin; i != end; ++i) + { + // need to check if the function is not the zero function. To + // check zero-ness, create a C++ function out of it and check that + if (static_cast &>(worker)) + worker(i, scratch_data, copy_data); + if (static_cast &>( + copier)) + copier(copy_data); + } + + return; + } + + tf::Executor &executor = MultithreadInfo::get_taskflow_executor(); + tf::Taskflow taskflow; + + ScratchData scratch_data = sample_scratch_data; + CopyData copy_data = sample_copy_data; // NOLINT + + tf::Task last_copier; + + std::vector> copy_datas; + + unsigned int idx = 0; + for (Iterator i = begin; i != end; ++i, ++idx) + { + copy_datas.emplace_back(); + + auto worker_task = taskflow + .emplace([it = i, + idx, + &sample_scratch_data, + ©_datas, + &sample_copy_data, + &worker]() { + // std::cout << "worker " << idx << std::endl; + ScratchData scratch = sample_scratch_data; + auto & copy = copy_datas[idx]; + copy = + std::make_unique(sample_copy_data); + + worker(it, scratch, *copy.get()); + }) + .name("worker"); + + tf::Task copier_task = taskflow + .emplace([idx, ©_datas, &copier]() { + copier(*copy_datas[idx].get()); + copy_datas[idx].reset(); + }) + .name("copy"); + + worker_task.precede(copier_task); + + if (!last_copier.empty()) + last_copier.precede(copier_task); + last_copier = copier_task; + } + + + // debugging: + + executor.run(taskflow).wait(); + +#ifdef DEBUG + std::cout << "done" << std::endl; + std::ofstream f("graph.dia"); + taskflow.dump(f); + f.close(); +#endif + } + + + template + void + run(const Iterator & begin, + const typename identity::type &end, + MainClass & main_object, + void (MainClass::*worker)(const Iterator &, ScratchData &, CopyData &), + void (MainClass::*copier)(const CopyData &), + const ScratchData &sample_scratch_data, + const CopyData & sample_copy_data, + const unsigned int queue_length = 2 * MultithreadInfo::n_threads(), + const unsigned int chunk_size = 8) + { + // forward to the other function + run(begin, + end, + [&main_object, worker](const Iterator &iterator, + ScratchData & scratch_data, + CopyData & copy_data) { + (main_object.*worker)(iterator, scratch_data, copy_data); + }, + [&main_object, copier](const CopyData ©_data) { + (main_object.*copier)(copy_data); + }, + sample_scratch_data, + sample_copy_data, + queue_length, + chunk_size); + } + } // namespace taskflow_v1 + + + // Like v1, except that we do not create a big task graph, but we submit 512 + // jobs at a time. + namespace taskflow_v2 + { + template + void run(const Iterator & begin, + const typename identity::type &end, + Worker worker, + Copier copier, + const ScratchData & sample_scratch_data, + const CopyData & sample_copy_data, + const unsigned int queue_length = 10 * + MultithreadInfo::n_threads(), + const unsigned int chunk_size = 8) + { + if (MultithreadInfo::n_threads() == 1) + { + // need to copy the sample since it is marked const + ScratchData scratch_data = sample_scratch_data; + CopyData copy_data = sample_copy_data; // NOLINT + + for (Iterator i = begin; i != end; ++i) + { + // need to check if the function is not the zero function. To + // check zero-ness, create a C++ function out of it and check that + if (static_cast &>(worker)) + worker(i, scratch_data, copy_data); + if (static_cast &>( + copier)) + copier(copy_data); + } + + return; + } + + const unsigned int max_work_tasks = 512; + + tf::Executor &executor = MultithreadInfo::get_taskflow_executor(); + tf::Taskflow taskflow; + + ScratchData scratch_data = sample_scratch_data; + CopyData copy_data = sample_copy_data; // NOLINT + + tf::Task last_copier; + + // Threads::ThreadLocalStorage> + // thread_local_scratch; + + std::vector> copy_datas; + + unsigned int idx = 0; + + for (Iterator i = begin; i != end; ++i, ++idx) + { + if (idx == max_work_tasks) + { + executor.run(taskflow).wait(); + taskflow.clear(); + copy_datas.clear(); + last_copier.reset(); + idx = 0; + } + + copy_datas.emplace_back(); + + auto worker_task = taskflow + .emplace([it = i, + idx, + &sample_scratch_data, + ©_datas, + &sample_copy_data, + &worker]() { + ScratchData scratch = sample_scratch_data; + auto & copy = copy_datas[idx]; + copy = + std::make_unique(sample_copy_data); + + worker(it, scratch, *copy.get()); + }) + .name("worker"); + + tf::Task copier_task = taskflow + .emplace([idx, ©_datas, &copier]() { + copier(*copy_datas[idx].get()); + copy_datas[idx].reset(); + }) + .name("copy"); + + worker_task.precede(copier_task); + + if (!last_copier.empty()) + last_copier.precede(copier_task); + last_copier = copier_task; + } + + + // debugging: + + executor.run(taskflow).wait(); + +#ifdef DEBUG + std::cout << "done" << std::endl; + std::ofstream f("graph.dia"); + taskflow.dump(f); + f.close(); +#endif + } + + + template + void + run(const Iterator & begin, + const typename identity::type &end, + MainClass & main_object, + void (MainClass::*worker)(const Iterator &, ScratchData &, CopyData &), + void (MainClass::*copier)(const CopyData &), + const ScratchData &sample_scratch_data, + const CopyData & sample_copy_data, + const unsigned int queue_length = 2 * MultithreadInfo::n_threads(), + const unsigned int chunk_size = 8) + { + // forward to the other function + run(begin, + end, + [&main_object, worker](const Iterator &iterator, + ScratchData & scratch_data, + CopyData & copy_data) { + (main_object.*worker)(iterator, scratch_data, copy_data); + }, + [&main_object, copier](const CopyData ©_data) { + (main_object.*copier)(copy_data); + }, + sample_scratch_data, + sample_copy_data, + queue_length, + chunk_size); + } + } // namespace taskflow_v2 + + // Like v1 but work in chunks of size 8 + namespace taskflow_v3 + { + template + class Chunk + { + public: + Chunk(const unsigned int count, const CopyData ©_data) + : copy_datas(count, copy_data) + {} + std::vector copy_datas; + }; + + + template + void run(const Iterator & begin, + const typename identity::type &end, + Worker worker, + Copier copier, + const ScratchData & sample_scratch_data, + const CopyData & sample_copy_data, + const unsigned int queue_length = 10 * + MultithreadInfo::n_threads(), + const unsigned int chunk_size = 8) + { + if (MultithreadInfo::n_threads() == 1) + { + // need to copy the sample since it is marked const + ScratchData scratch_data = sample_scratch_data; + CopyData copy_data = sample_copy_data; // NOLINT + + for (Iterator i = begin; i != end; ++i) + { + // need to check if the function is not the zero function. To + // check zero-ness, create a C++ function out of it and check that + if (static_cast &>(worker)) + worker(i, scratch_data, copy_data); + if (static_cast &>( + copier)) + copier(copy_data); + } + + return; + } + + tf::Executor &executor = MultithreadInfo::get_taskflow_executor(); + tf::Taskflow taskflow; + + ScratchData scratch_data = sample_scratch_data; + CopyData copy_data = sample_copy_data; // NOLINT + + tf::Task last_copier = taskflow.placeholder(); + + std::vector>> chunks; + + unsigned int idx = 0; + unsigned int remaining_items = std::distance(begin, end); + + Iterator it = begin; + while (it != end) + { + unsigned int count = std::min(remaining_items, chunk_size); + Iterator middle = it; + std::advance(middle, count); + + chunks.emplace_back(); + + // this chunk works on [it,middle) + auto worker_task = + taskflow + .emplace([it_begin = it, + it_end = middle, + idx, + count, + &sample_scratch_data, + &chunks, + &sample_copy_data, + &worker]() { + ScratchData scratch = sample_scratch_data; + chunks[idx] = + std::make_unique>(count, sample_copy_data); + + unsigned int counter = 0; + for (Iterator it = it_begin; it != it_end; ++it, ++counter) + { + worker(it, scratch, chunks[idx].get()->copy_datas[counter]); + } + }) + .name("work"); + + tf::Task copier_task = taskflow + .emplace([idx, &chunks, &copier]() mutable { + auto chunk = chunks[idx].get(); + for (auto &cd : chunk->copy_datas) + copier(cd); + + chunks[idx].reset(); + }) + .name("copy"); + + worker_task.precede(copier_task); + + last_copier.precede(copier_task); + last_copier = copier_task; + + it = middle; + ++idx; + } + // copy_datas.emplace_back(); + + // auto worker_task = taskflow + // .emplace([it = i, + // idx, + // &thread_local_scratch, + // &sample_scratch_data, + // ©_datas, + // &sample_copy_data, + // &worker]() { + // // std::cout << "worker " << idx << + // std::endl; auto &scratch_ptr = + // thread_local_scratch.get(); if + // (!scratch_ptr.get()) + // { + // thread_local_scratch = + // std::make_unique(sample_scratch_data); + // scratch_ptr = + // thread_local_scratch.get(); + // } + // auto &scratch_ptr = + // thread_local_scratch.get(); + + // auto & copy = copy_datas[idx]; + // copy = + // std::make_unique(sample_copy_data); + + // worker(it, *scratch_ptr.get(), + // *copy.get()); + // }) + // .name("worker"); + + // tf::Task copier_task = taskflow + // .emplace([idx, ©_datas, + // &copier]() { + // copier(*copy_datas[idx].get()); + // copy_datas[idx].reset(); + // }) + // .name("copy"); + + // worker_task.precede(copier_task); + + // if (!last_copier.empty()) + // last_copier.precede(copier_task); + // last_copier = copier_task; + // } + + + // debugging: + + executor.run(taskflow).wait(); + +#ifdef DEBUG + std::cout << "done" << std::endl; + std::ofstream f("graph.dia"); + taskflow.dump(f); + f.close(); +#endif + } + + + template + void + run(const Iterator & begin, + const typename identity::type &end, + MainClass & main_object, + void (MainClass::*worker)(const Iterator &, ScratchData &, CopyData &), + void (MainClass::*copier)(const CopyData &), + const ScratchData &sample_scratch_data, + const CopyData & sample_copy_data, + const unsigned int queue_length = 2 * MultithreadInfo::n_threads(), + const unsigned int chunk_size = 8) + { + // forward to the other function + run(begin, + end, + [&main_object, worker](const Iterator &iterator, + ScratchData & scratch_data, + CopyData & copy_data) { + (main_object.*worker)(iterator, scratch_data, copy_data); + }, + [&main_object, copier](const CopyData ©_data) { + (main_object.*copier)(copy_data); + }, + sample_scratch_data, + sample_copy_data, + queue_length, + chunk_size); + } + } // namespace taskflow_v3 + + // Work in chunks of size 8 (only if we have enough items, otherwise 1) and + // use a thread-local scratch object + namespace taskflow_v4 + { + template + class Chunk + { + public: + Chunk(const unsigned int count, const CopyData ©_data) + : copy_datas(count, copy_data) + {} + std::vector copy_datas; + }; + + + template + void run(const Iterator & begin, + const typename identity::type &end, + Worker worker, + Copier copier, + const ScratchData & sample_scratch_data, + const CopyData & sample_copy_data, + const unsigned int queue_length = 10 * + MultithreadInfo::n_threads(), + const unsigned int chunk_size = 8) + { + if (MultithreadInfo::n_threads() == 1) + { + // need to copy the sample since it is marked const + ScratchData scratch_data = sample_scratch_data; + CopyData copy_data = sample_copy_data; // NOLINT + + for (Iterator i = begin; i != end; ++i) + { + // need to check if the function is not the zero function. To + // check zero-ness, create a C++ function out of it and check that + if (static_cast &>(worker)) + worker(i, scratch_data, copy_data); + if (static_cast &>( + copier)) + copier(copy_data); + } + + return; + } + + tf::Executor &executor = MultithreadInfo::get_taskflow_executor(); + tf::Taskflow taskflow; + + tf::Task last_copier = taskflow.placeholder(); + + Threads::ThreadLocalStorage> + thread_local_scratch; + + std::vector>> chunks; + + unsigned int idx = 0; + unsigned int remaining_items = std::distance(begin, end); + + const unsigned int real_chunk_size = + (remaining_items / chunk_size < 3 * MultithreadInfo::n_threads()) ? + 1 : + chunk_size; + + Iterator it = begin; + while (it != end) + { + unsigned int count = std::min(remaining_items, real_chunk_size); + Iterator middle = it; + std::advance(middle, count); + + chunks.emplace_back(); + + // this chunk works on [it,middle) + auto worker_task = + taskflow + .emplace([it_begin = it, + it_end = middle, + idx, + count, + &sample_scratch_data, + &thread_local_scratch, + &chunks, + &sample_copy_data, + &worker]() { + auto &scratch_ptr = thread_local_scratch.get(); + if (!scratch_ptr.get()) + scratch_ptr = + std::make_unique(sample_scratch_data); + + ScratchData &scratch = *scratch_ptr.get(); + chunks[idx] = + std::make_unique>(count, sample_copy_data); + + unsigned int counter = 0; + for (Iterator it = it_begin; it != it_end; ++it, ++counter) + { + worker(it, scratch, chunks[idx].get()->copy_datas[counter]); + } + }) + .name("work"); + + tf::Task copier_task = taskflow + .emplace([idx, &chunks, &copier]() mutable { + auto chunk = chunks[idx].get(); + for (auto &cd : chunk->copy_datas) + copier(cd); + + chunks[idx].reset(); + }) + .name("copy"); + + worker_task.precede(copier_task); + + last_copier.precede(copier_task); + last_copier = copier_task; + + it = middle; + ++idx; + } + + // debugging: + + executor.run(taskflow).wait(); + +#ifdef DEBUG + std::cout << "done" << std::endl; + std::ofstream f("graph.dia"); + taskflow.dump(f); + f.close(); +#endif + } + + + template + void + run(const Iterator & begin, + const typename identity::type &end, + MainClass & main_object, + void (MainClass::*worker)(const Iterator &, ScratchData &, CopyData &), + void (MainClass::*copier)(const CopyData &), + const ScratchData &sample_scratch_data, + const CopyData & sample_copy_data, + const unsigned int queue_length = 2 * MultithreadInfo::n_threads(), + const unsigned int chunk_size = 8) + { + // forward to the other function + run(begin, + end, + [&main_object, worker](const Iterator &iterator, + ScratchData & scratch_data, + CopyData & copy_data) { + (main_object.*worker)(iterator, scratch_data, copy_data); + }, + [&main_object, copier](const CopyData ©_data) { + (main_object.*copier)(copy_data); + }, + sample_scratch_data, + sample_copy_data, + queue_length, + chunk_size); + } + } // namespace taskflow_v4 + + // In the following function, the matrix and right hand side are + // assembled. As stated in the documentation of the main class above, it + // does not do this itself, but rather delegates to the function following + // next, utilizing the WorkStream concept discussed in @ref threads . + // + // If you have looked through the @ref threads module, you will have + // seen that assembling in parallel does not take an incredible + // amount of extra code as long as you diligently describe what the + // scratch and copy data objects are, and if you define suitable + // functions for the local assembly and the copy operation from local + // contributions to global objects. This done, the following will do + // all the heavy lifting to get these operations done on multiple + // threads on as many cores as you have in your system: + template + void AdvectionProblem::assemble_system() + { + if (true) + taskflow_v1::run(dof_handler.begin_active(), + dof_handler.end(), + *this, + &AdvectionProblem::local_assemble_system, + &AdvectionProblem::copy_local_to_global, + AssemblyScratchData(fe), + AssemblyCopyData()); + else + WorkStream::run(dof_handler.begin_active(), + dof_handler.end(), + *this, + &AdvectionProblem::local_assemble_system, + &AdvectionProblem::copy_local_to_global, + AssemblyScratchData(fe), + AssemblyCopyData()); + } + + + + // As already mentioned above, we need to have scratch objects for + // the parallel computation of local contributions. These objects + // contain FEValues and FEFaceValues objects (as well as some arrays), and so + // we will need to have constructors and copy constructors that allow us to + // create them. For the cell terms we need the values + // and gradients of the shape functions, the quadrature points in + // order to determine the source density and the advection field at + // a given point, and the weights of the quadrature points times the + // determinant of the Jacobian at these points. In contrast, for the + // boundary integrals, we don't need the gradients, but rather the + // normal vectors to the cells. This determines which update flags + // we will have to pass to the constructors of the members of the + // class: + template + AdvectionProblem::AssemblyScratchData::AssemblyScratchData( + const FiniteElement &fe) + : fe_values(fe, + QGauss(fe.degree + 1), + update_values | update_gradients | update_quadrature_points | + update_JxW_values) + , fe_face_values(fe, + QGauss(fe.degree + 1), + update_values | update_quadrature_points | + update_JxW_values | update_normal_vectors) + , rhs_values(fe_values.get_quadrature().size()) + , advection_directions(fe_values.get_quadrature().size()) + , face_boundary_values(fe_face_values.get_quadrature().size()) + , face_advection_directions(fe_face_values.get_quadrature().size()) + {} + + + + template + AdvectionProblem::AssemblyScratchData::AssemblyScratchData( + const AssemblyScratchData &scratch_data) + : fe_values(scratch_data.fe_values.get_fe(), + scratch_data.fe_values.get_quadrature(), + update_values | update_gradients | update_quadrature_points | + update_JxW_values) + , fe_face_values(scratch_data.fe_face_values.get_fe(), + scratch_data.fe_face_values.get_quadrature(), + update_values | update_quadrature_points | + update_JxW_values | update_normal_vectors) + , rhs_values(scratch_data.rhs_values.size()) + , advection_directions(scratch_data.advection_directions.size()) + , face_boundary_values(scratch_data.face_boundary_values.size()) + , face_advection_directions(scratch_data.face_advection_directions.size()) + {} + + + + // Now, this is the function that does the actual work. It is not very + // different from the assemble_system functions of previous + // example programs, so we will again only comment on the differences. The + // mathematical stuff closely follows what we have said in the introduction. + // + // There are a number of points worth mentioning here, though. The + // first one is that we have moved the FEValues and FEFaceValues + // objects into the ScratchData object. We have done so because the + // alternative would have been to simply create one every time we + // get into this function -- i.e., on every cell. It now turns out + // that the FEValues classes were written with the explicit goal of + // moving everything that remains the same from cell to cell into + // the construction of the object, and only do as little work as + // possible in FEValues::reinit() whenever we move to a new + // cell. What this means is that it would be very expensive to + // create a new object of this kind in this function as we would + // have to do it for every cell -- exactly the thing we wanted to + // avoid with the FEValues class. Instead, what we do is create it + // only once (or a small number of times) in the scratch objects and + // then re-use it as often as we can. + // + // This begs the question of whether there are other objects we + // create in this function whose creation is expensive compared to + // its use. Indeed, at the top of the function, we declare all sorts + // of objects. The AdvectionField, + // RightHandSide and BoundaryValues do not + // cost much to create, so there is no harm here. However, + // allocating memory in creating the rhs_values and + // similar variables below typically costs a significant amount of + // time, compared to just accessing the (temporary) values we store + // in them. Consequently, these would be candidates for moving into + // the AssemblyScratchData class. We will leave this as + // an exercise. + template + void AdvectionProblem::local_assemble_system( + const typename DoFHandler::active_cell_iterator &cell, + AssemblyScratchData & scratch_data, + AssemblyCopyData & copy_data) + { + // We define some abbreviations to avoid unnecessarily long lines: + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = + scratch_data.fe_values.get_quadrature().size(); + const unsigned int n_face_q_points = + scratch_data.fe_face_values.get_quadrature().size(); + + // We declare cell matrix and cell right hand side... + copy_data.cell_matrix.reinit(dofs_per_cell, dofs_per_cell); + copy_data.cell_rhs.reinit(dofs_per_cell); + + // ... an array to hold the global indices of the degrees of freedom of + // the cell on which we are presently working... + copy_data.local_dof_indices.resize(dofs_per_cell); + + // ... then initialize the FEValues object... + scratch_data.fe_values.reinit(cell); + + // ... obtain the values of right hand side and advection directions + // at the quadrature points... + scratch_data.advection_field.value_list( + scratch_data.fe_values.get_quadrature_points(), + scratch_data.advection_directions); + scratch_data.right_hand_side.value_list( + scratch_data.fe_values.get_quadrature_points(), scratch_data.rhs_values); + + // ... set the value of the streamline diffusion parameter as + // described in the introduction... + const double delta = 0.1 * cell->diameter(); + + // ... and assemble the local contributions to the system matrix and + // right hand side as also discussed above: + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + // Alias the AssemblyScratchData object to keep the lines from + // getting too long: + const auto &sd = scratch_data; + for (unsigned int j = 0; j < dofs_per_cell; ++j) + copy_data.cell_matrix(i, j) += + ((sd.fe_values.shape_value(i, q_point) + // (phi_i + + delta * (sd.advection_directions[q_point] * // delta beta + sd.fe_values.shape_grad(i, q_point))) * // grad phi_i) + sd.advection_directions[q_point] * // beta + sd.fe_values.shape_grad(j, q_point)) * // grad phi_j + sd.fe_values.JxW(q_point); // dx + + copy_data.cell_rhs(i) += + (sd.fe_values.shape_value(i, q_point) + // (phi_i + + delta * (sd.advection_directions[q_point] * // delta beta + sd.fe_values.shape_grad(i, q_point))) * // grad phi_i) + sd.rhs_values[q_point] * // f + sd.fe_values.JxW(q_point); // dx + } + + // Besides the cell terms which we have built up now, the bilinear + // form of the present problem also contains terms on the boundary of + // the domain. Therefore, we have to check whether any of the faces of + // this cell are on the boundary of the domain, and if so assemble the + // contributions of this face as well. Of course, the bilinear form + // only contains contributions from the inflow part of + // the boundary, but to find out whether a certain part of a face of + // the present cell is part of the inflow boundary, we have to have + // information on the exact location of the quadrature points and on + // the direction of flow at this point; we obtain this information + // using the FEFaceValues object and only decide within the main loop + // whether a quadrature point is on the inflow boundary. + for (const auto &face : cell->face_iterators()) + if (face->at_boundary()) + { + // Ok, this face of the present cell is on the boundary of the + // domain. Just as for the usual FEValues object which we have + // used in previous examples and also above, we have to + // reinitialize the FEFaceValues object for the present face: + scratch_data.fe_face_values.reinit(cell, face); + + // For the quadrature points at hand, we ask for the values of + // the inflow function and for the direction of flow: + scratch_data.boundary_values.value_list( + scratch_data.fe_face_values.get_quadrature_points(), + scratch_data.face_boundary_values); + scratch_data.advection_field.value_list( + scratch_data.fe_face_values.get_quadrature_points(), + scratch_data.face_advection_directions); + + // Now loop over all quadrature points and see whether this face is on + // the inflow or outflow part of the boundary. The normal + // vector points out of the cell: since the face is at + // the boundary, the normal vector points out of the domain, + // so if the advection direction points into the domain, its + // scalar product with the normal vector must be negative (to see why + // this is true, consider the scalar product definition that uses a + // cosine): + for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point) + if (scratch_data.fe_face_values.normal_vector(q_point) * + scratch_data.face_advection_directions[q_point] < + 0.) + // If the face is part of the inflow boundary, then compute the + // contributions of this face to the global matrix and right + // hand side, using the values obtained from the + // FEFaceValues object and the formulae discussed in the + // introduction: + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + copy_data.cell_matrix(i, j) -= + (scratch_data.face_advection_directions[q_point] * + scratch_data.fe_face_values.normal_vector(q_point) * + scratch_data.fe_face_values.shape_value(i, q_point) * + scratch_data.fe_face_values.shape_value(j, q_point) * + scratch_data.fe_face_values.JxW(q_point)); + + copy_data.cell_rhs(i) -= + (scratch_data.face_advection_directions[q_point] * + scratch_data.fe_face_values.normal_vector(q_point) * + scratch_data.face_boundary_values[q_point] * + scratch_data.fe_face_values.shape_value(i, q_point) * + scratch_data.fe_face_values.JxW(q_point)); + } + } + + // The final piece of information the copy routine needs is the global + // indices of the degrees of freedom on this cell, so we end by writing + // them to the local array: + cell->get_dof_indices(copy_data.local_dof_indices); + } + + + + // The second function we needed to write was the one that copies + // the local contributions the previous function computed (and + // put into the AssemblyCopyData object) into the global matrix and right + // hand side vector objects. This is essentially what we always had + // as the last block of code when assembling something on every + // cell. The following should therefore be pretty obvious: + template + void + AdvectionProblem::copy_local_to_global(const AssemblyCopyData ©_data) + { + hanging_node_constraints.distribute_local_to_global( + copy_data.cell_matrix, + copy_data.cell_rhs, + copy_data.local_dof_indices, + system_matrix, + system_rhs); + } + + // Here comes the linear solver routine. As the system is no longer + // symmetric positive definite as in all the previous examples, we cannot + // use the Conjugate Gradient method anymore. Rather, we use a solver that + // is more general and does not rely on any special properties of the + // matrix: the GMRES method. GMRES, like the conjugate gradient method, + // requires a decent preconditioner: we use a Jacobi preconditioner here, + // which works well enough for this problem. + template + void AdvectionProblem::solve() + { + SolverControl solver_control(std::max(1000, + system_rhs.size() / 10), + 1e-10 * system_rhs.l2_norm()); + SolverGMRES> solver(solver_control); + PreconditionJacobi> preconditioner; + preconditioner.initialize(system_matrix, 1.0); + solver.solve(system_matrix, solution, system_rhs, preconditioner); + + Vector residual(dof_handler.n_dofs()); + + system_matrix.vmult(residual, solution); + residual -= system_rhs; + std::cout << " Iterations required for convergence: " + << solver_control.last_step() << '\n' + << " Max norm of residual: " + << residual.linfty_norm() << '\n'; + + hanging_node_constraints.distribute(solution); + } + + // The following function refines the grid according to the quantity + // described in the introduction. The respective computations are made in + // the class GradientEstimation. + template + void AdvectionProblem::refine_grid() + { + Vector estimated_error_per_cell(triangulation.n_active_cells()); + + GradientEstimation::estimate(dof_handler, + solution, + estimated_error_per_cell); + + GridRefinement::refine_and_coarsen_fixed_number(triangulation, + estimated_error_per_cell, + 0.3, + 0.03); + + triangulation.execute_coarsening_and_refinement(); + } + + // This function is similar to the one in step 6, but since we use a higher + // degree finite element we save the solution in a different + // way. Visualization programs like VisIt and Paraview typically only + // understand data that is associated with nodes: they cannot plot + // fifth-degree basis functions, which results in a very inaccurate picture + // of the solution we computed. To get around this we save multiple + // patches per cell: in 2D we save 64 bilinear `cells' to the VTU + // file for each cell, and in 3D we save 512. The end result is that the + // visualization program will use a piecewise linear interpolation of the + // cubic basis functions: this captures the solution detail and, with most + // screen resolutions, looks smooth. We save the grid in a separate step + // with no extra patches so that we have a visual representation of the cell + // faces. + // + // Version 9.1 of deal.II gained the ability to write higher degree + // polynomials (i.e., write piecewise bicubic visualization data for our + // piecewise bicubic solution) VTK and VTU output: however, not all recent + // versions of ParaView and VisIt (as of 2018) can read this format, so we + // use the older, more general (but less efficient) approach here. + template + void AdvectionProblem::output_results(const unsigned int cycle) const + { + { + GridOut grid_out; + std::ofstream output("grid-" + std::to_string(cycle) + ".vtu"); + grid_out.write_vtu(triangulation, output); + } + + { + DataOut data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(solution, "solution"); + data_out.build_patches(8); + + // VTU output can be expensive, both to compute and to write to + // disk. Here we ask ZLib, a compression library, to compress the data + // in a way that maximizes throughput. + DataOutBase::VtkFlags vtk_flags; + vtk_flags.compression_level = + DataOutBase::VtkFlags::ZlibCompressionLevel::best_speed; + data_out.set_flags(vtk_flags); + + std::ofstream output("solution-" + std::to_string(cycle) + ".vtu"); + data_out.write_vtu(output); + } + } + + + // ... as is the main loop (setup -- solve -- refine), aside from the number + // of cycles and the initial grid: + template + void AdvectionProblem::run() + { + for (unsigned int cycle = 0; cycle < 1; ++cycle) + { + std::cout << "Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + { + GridGenerator::hyper_cube(triangulation, -1, 1); + triangulation.refine_global(2); // 7 + } + else + { + triangulation.refine_global(2); + } + + + std::cout << " Number of active cells: " + << triangulation.n_active_cells() << std::endl; + + setup_system(); + + std::cout << " Number of degrees of freedom: " + << dof_handler.n_dofs() << std::endl; + + Timer timer; + + const unsigned int n_phys_cores = MultithreadInfo::n_cores(); + std::cout << "MultithreadInfo::n_cores()=" << n_phys_cores << std::endl; + std::cout << "MultithreadInfo::n_threads()=" + << MultithreadInfo::n_threads() << std::endl; + + const unsigned int runs = 5; + +#ifdef DEAL_II_WITH_CPP_TASKFLOW + std::cout << "** TASKFLOW v1 **" << std::endl; + + + MultithreadInfo::set_thread_limit(); + + for (unsigned int n_cores = 2 * n_phys_cores; n_cores > 0; n_cores /= 2) + { + if (n_cores <= n_phys_cores) + MultithreadInfo::set_thread_limit(n_cores); + + std::cout << "n_cores "; + if (n_cores <= n_phys_cores) + std::cout << n_cores; + else + std::cout << "auto"; + std::cout << ' ' << std::flush; + double avg = 0.; + + for (unsigned int c = 0; c < runs; ++c) + { + timer.reset(); + timer.start(); + + taskflow_v1::run(dof_handler.begin_active(), + dof_handler.end(), + *this, + &AdvectionProblem::local_assemble_system, + &AdvectionProblem::copy_local_to_global, + AssemblyScratchData(fe), + AssemblyCopyData()); + + + timer.stop(); + const double time = timer.last_wall_time(); + avg += time; + std::cout << time << " " << std::flush; + } + avg /= runs; + std::cout << " avg: " << avg << std::endl; + } + + std::cout << "** TASKFLOW v2 **" << std::endl; + + + MultithreadInfo::set_thread_limit(); + + for (unsigned int n_cores = 2 * n_phys_cores; n_cores > 0; n_cores /= 2) + { + if (n_cores <= n_phys_cores) + MultithreadInfo::set_thread_limit(n_cores); + + std::cout << "n_cores "; + if (n_cores <= n_phys_cores) + std::cout << n_cores; + else + std::cout << "auto"; + std::cout << ' ' << std::flush; + double avg = 0.; + + for (unsigned int c = 0; c < runs; ++c) + { + timer.reset(); + timer.start(); + + taskflow_v2::run(dof_handler.begin_active(), + dof_handler.end(), + *this, + &AdvectionProblem::local_assemble_system, + &AdvectionProblem::copy_local_to_global, + AssemblyScratchData(fe), + AssemblyCopyData()); + + + timer.stop(); + const double time = timer.last_wall_time(); + avg += time; + std::cout << time << " " << std::flush; + } + avg /= runs; + std::cout << " avg: " << avg << std::endl; + } + + + std::cout << "** TASKFLOW v3 **" << std::endl; + + + MultithreadInfo::set_thread_limit(); + + for (unsigned int n_cores = 2 * n_phys_cores; n_cores > 0; n_cores /= 2) + { + if (n_cores <= n_phys_cores) + MultithreadInfo::set_thread_limit(n_cores); + + std::cout << "n_cores "; + if (n_cores <= n_phys_cores) + std::cout << n_cores; + else + std::cout << "auto"; + std::cout << ' ' << std::flush; + double avg = 0.; + + for (unsigned int c = 0; c < runs; ++c) + { + timer.reset(); + timer.start(); + + taskflow_v3::run(dof_handler.begin_active(), + dof_handler.end(), + *this, + &AdvectionProblem::local_assemble_system, + &AdvectionProblem::copy_local_to_global, + AssemblyScratchData(fe), + AssemblyCopyData()); + + + timer.stop(); + const double time = timer.last_wall_time(); + avg += time; + std::cout << time << " " << std::flush; + } + avg /= runs; + std::cout << " avg: " << avg << std::endl; + } + + + + std::cout << "** TASKFLOW v4 **" << std::endl; + + + MultithreadInfo::set_thread_limit(); + + for (unsigned int n_cores = 2 * n_phys_cores; n_cores > 0; n_cores /= 2) + { + if (n_cores <= n_phys_cores) + MultithreadInfo::set_thread_limit(n_cores); + + std::cout << "n_cores "; + if (n_cores <= n_phys_cores) + std::cout << n_cores; + else + std::cout << "auto"; + std::cout << ' ' << std::flush; + double avg = 0.; + + for (unsigned int c = 0; c < runs; ++c) + { + timer.reset(); + timer.start(); + + taskflow_v4::run(dof_handler.begin_active(), + dof_handler.end(), + *this, + &AdvectionProblem::local_assemble_system, + &AdvectionProblem::copy_local_to_global, + AssemblyScratchData(fe), + AssemblyCopyData()); + + + timer.stop(); + const double time = timer.last_wall_time(); + avg += time; + std::cout << time << " " << std::flush; + } + avg /= runs; + std::cout << " avg: " << avg << std::endl; + } +#endif +#ifdef DEAL_II_WITH_TBB + std::cout << "** TBB **" << std::endl; + + + MultithreadInfo::set_thread_limit(); + + for (unsigned int n_cores = 2 * n_phys_cores; n_cores > 0; n_cores /= 2) + { + if (n_cores <= n_phys_cores) + MultithreadInfo::set_thread_limit(n_cores); + + std::cout << "n_cores "; + if (n_cores <= n_phys_cores) + std::cout << n_cores; + else + std::cout << "auto"; + std::cout << ' ' << std::flush; + double avg = 0.; + + for (unsigned int c = 0; c < runs; ++c) + { + timer.reset(); + timer.start(); + WorkStream::run(dof_handler.begin_active(), + dof_handler.end(), + *this, + &AdvectionProblem::local_assemble_system, + &AdvectionProblem::copy_local_to_global, + AssemblyScratchData(fe), + AssemblyCopyData()); + + + timer.stop(); + const double time = timer.last_wall_time(); + avg += time; + std::cout << time << " " << std::flush; + } + avg /= runs; + std::cout << " avg: " << avg << std::endl; + } +#endif + + // assemble_system(); + + + // solve(); + // output_results(cycle); + } + } + + + + // @sect3{GradientEstimation class implementation} + + // Now for the implementation of the GradientEstimation class. + // Let us start by defining constructors for the + // EstimateScratchData class used by the + // estimate_cell() function: + template + GradientEstimation::EstimateScratchData::EstimateScratchData( + const FiniteElement &fe, + const Vector & solution, + Vector & error_per_cell) + : fe_midpoint_value(fe, + QMidpoint(), + update_values | update_quadrature_points) + , solution(solution) + , error_per_cell(error_per_cell) + , cell_midpoint_value(1) + , neighbor_midpoint_value(1) + { + // We allocate a vector to hold iterators to all active neighbors of + // a cell. We reserve the maximal number of active neighbors in order to + // avoid later reallocations. Note how this maximal number of active + // neighbors is computed here. + active_neighbors.reserve(GeometryInfo::faces_per_cell * + GeometryInfo::max_children_per_face); + } + + + template + GradientEstimation::EstimateScratchData::EstimateScratchData( + const EstimateScratchData &scratch_data) + : fe_midpoint_value(scratch_data.fe_midpoint_value.get_fe(), + scratch_data.fe_midpoint_value.get_quadrature(), + update_values | update_quadrature_points) + , solution(scratch_data.solution) + , error_per_cell(scratch_data.error_per_cell) + , cell_midpoint_value(1) + , neighbor_midpoint_value(1) + {} + + + // Next comes the implementation of the GradientEstimation + // class. The first function does not much except for delegating work to the + // other function, but there is a bit of setup at the top. + // + // Before starting with the work, we check that the vector into + // which the results are written has the right size. Programming + // mistakes in which one forgets to size arguments correctly at the + // calling site are quite common. Because the resulting damage from + // not catching such errors is often subtle (e.g., corruption of + // data somewhere in memory, or non-reproducible results), it is + // well worth the effort to check for such things. + template + void GradientEstimation::estimate(const DoFHandler &dof_handler, + const Vector & solution, + Vector & error_per_cell) + { + Assert( + error_per_cell.size() == dof_handler.get_triangulation().n_active_cells(), + ExcInvalidVectorLength(error_per_cell.size(), + dof_handler.get_triangulation().n_active_cells())); + + WorkStream::run(dof_handler.begin_active(), + dof_handler.end(), + &GradientEstimation::template estimate_cell, + std::function(), + EstimateScratchData(dof_handler.get_fe(), + solution, + error_per_cell), + EstimateCopyData()); + } + + + // Here comes the function that estimates the local error by computing the + // finite difference approximation of the gradient. The function first + // computes the list of active neighbors of the present cell and then + // computes the quantities described in the introduction for each of + // the neighbors. The reason for this order is that it is not a one-liner + // to find a given neighbor with locally refined meshes. In principle, an + // optimized implementation would find neighbors and the quantities + // depending on them in one step, rather than first building a list of + // neighbors and in a second step their contributions but we will gladly + // leave this as an exercise. As discussed before, the worker function + // passed to WorkStream::run works on "scratch" objects that keep all + // temporary objects. This way, we do not need to create and initialize + // objects that are expensive to initialize within the function that does + // the work every time it is called for a given cell. Such an argument is + // passed as the second argument. The third argument would be a "copy-data" + // object (see @ref threads for more information) but we do not actually use + // any of these here. Since WorkStream::run() insists on passing three + // arguments, we declare this function with three arguments, but simply + // ignore the last one. + // + // (This is unsatisfactory from an aesthetic perspective. It can be avoided + // by using an anonymous (lambda) function. If you allow, let us here show + // how. First, assume that we had declared this function to only take two + // arguments by omitting the unused last one. Now, WorkStream::run still + // wants to call this function with three arguments, so we need to find a + // way to "forget" the third argument in the call. Simply passing + // WorkStream::run the pointer to the function as we do above will not do + // this -- the compiler will complain that a function declared to have two + // arguments is called with three arguments. However, we can do this by + // passing the following as the third argument to WorkStream::run(): + // @code + // [](const typename DoFHandler::active_cell_iterator &cell, + // EstimateScratchData & scratch_data, + // EstimateCopyData &) + // { + // GradientEstimation::estimate_cell(cell, scratch_data); + // } + // @endcode + // This is not much better than the solution implemented below: either the + // routine itself must take three arguments or it must be wrapped by + // something that takes three arguments. We don't use this since adding the + // unused argument at the beginning is simpler. + // + // Now for the details: + template + void GradientEstimation::estimate_cell( + const typename DoFHandler::active_cell_iterator &cell, + EstimateScratchData & scratch_data, + const EstimateCopyData &) + { + // We need space for the tensor Y, which is the sum of + // outer products of the y-vectors. + Tensor<2, dim> Y; + + // First initialize the FEValues object, as well as the + // Y tensor: + scratch_data.fe_midpoint_value.reinit(cell); + + // Now, before we go on, we first compute a list of all active neighbors + // of the present cell. We do so by first looping over all faces and see + // whether the neighbor there is active, which would be the case if it + // is on the same level as the present cell or one level coarser (note + // that a neighbor can only be once coarser than the present cell, as + // we only allow a maximal difference of one refinement over a face in + // deal.II). Alternatively, the neighbor could be on the same level + // and be further refined; then we have to find which of its children + // are next to the present cell and select these (note that if a child + // of a neighbor of an active cell that is next to this active cell, + // needs necessarily be active itself, due to the one-refinement rule + // cited above). + // + // Things are slightly different in one space dimension, as there the + // one-refinement rule does not exist: neighboring active cells may + // differ in as many refinement levels as they like. In this case, the + // computation becomes a little more difficult, but we will explain + // this below. + // + // Before starting the loop over all neighbors of the present cell, we + // have to clear the array storing the iterators to the active + // neighbors, of course. + scratch_data.active_neighbors.clear(); + for (unsigned int face_n : GeometryInfo::face_indices()) + if (!cell->at_boundary(face_n)) + { + // First define an abbreviation for the iterator to the face and + // the neighbor + const auto face = cell->face(face_n); + const auto neighbor = cell->neighbor(face_n); + + // Then check whether the neighbor is active. If it is, then it + // is on the same level or one level coarser (if we are not in + // 1D), and we are interested in it in any case. + if (neighbor->is_active()) + scratch_data.active_neighbors.push_back(neighbor); + else + { + // If the neighbor is not active, then check its children. + if (dim == 1) + { + // To find the child of the neighbor which bounds to the + // present cell, successively go to its right child if + // we are left of the present cell (n==0), or go to the + // left child if we are on the right (n==1), until we + // find an active cell. + auto neighbor_child = neighbor; + while (neighbor_child->has_children()) + neighbor_child = neighbor_child->child(face_n == 0 ? 1 : 0); + + // As this used some non-trivial geometrical intuition, + // we might want to check whether we did it right, + // i.e., check whether the neighbor of the cell we found + // is indeed the cell we are presently working + // on. Checks like this are often useful and have + // frequently uncovered errors both in algorithms like + // the line above (where it is simple to involuntarily + // exchange n==1 for n==0 or + // the like) and in the library (the assumptions + // underlying the algorithm above could either be wrong, + // wrongly documented, or are violated due to an error + // in the library). One could in principle remove such + // checks after the program works for some time, but it + // might be a good things to leave it in anyway to check + // for changes in the library or in the algorithm above. + // + // Note that if this check fails, then this is certainly + // an error that is irrecoverable and probably qualifies + // as an internal error. We therefore use a predefined + // exception class to throw here. + Assert(neighbor_child->neighbor(face_n == 0 ? 1 : 0) == cell, + ExcInternalError()); + + // If the check succeeded, we push the active neighbor + // we just found to the stack we keep: + scratch_data.active_neighbors.push_back(neighbor_child); + } + else + // If we are not in 1d, we collect all neighbor children + // `behind' the subfaces of the current face and move on: + for (unsigned int subface_n = 0; subface_n < face->n_children(); + ++subface_n) + scratch_data.active_neighbors.push_back( + cell->neighbor_child_on_subface(face_n, subface_n)); + } + } + + // OK, now that we have all the neighbors, lets start the computation + // on each of them. First we do some preliminaries: find out about the + // center of the present cell and the solution at this point. The + // latter is obtained as a vector of function values at the quadrature + // points, of which there are only one, of course. Likewise, the + // position of the center is the position of the first (and only) + // quadrature point in real space. + const Point this_center = + scratch_data.fe_midpoint_value.quadrature_point(0); + + scratch_data.fe_midpoint_value.get_function_values( + scratch_data.solution, scratch_data.cell_midpoint_value); + + // Now loop over all active neighbors and collect the data we + // need. + Tensor<1, dim> projected_gradient; + for (const auto &neighbor : scratch_data.active_neighbors) + { + // Then get the center of the neighbor cell and the value of the + // finite element function at that point. Note that for this + // information we have to reinitialize the FEValues + // object for the neighbor cell. + scratch_data.fe_midpoint_value.reinit(neighbor); + const Point neighbor_center = + scratch_data.fe_midpoint_value.quadrature_point(0); + + scratch_data.fe_midpoint_value.get_function_values( + scratch_data.solution, scratch_data.neighbor_midpoint_value); + + // Compute the vector y connecting the centers of the + // two cells. Note that as opposed to the introduction, we denote + // by y the normalized difference vector, as this is + // the quantity used everywhere in the computations. + Tensor<1, dim> y = neighbor_center - this_center; + const double distance = y.norm(); + y /= distance; + + // Then add up the contribution of this cell to the Y matrix... + for (unsigned int i = 0; i < dim; ++i) + for (unsigned int j = 0; j < dim; ++j) + Y[i][j] += y[i] * y[j]; + + // ... and update the sum of difference quotients: + projected_gradient += (scratch_data.neighbor_midpoint_value[0] - + scratch_data.cell_midpoint_value[0]) / + distance * y; + } + + // If now, after collecting all the information from the neighbors, we + // can determine an approximation of the gradient for the present + // cell, then we need to have passed over vectors y which + // span the whole space, otherwise we would not have all components of + // the gradient. This is indicated by the invertibility of the matrix. + // + // If the matrix is not invertible, then the present + // cell had an insufficient number of active neighbors. In contrast to + // all previous cases (where we raised exceptions) this is, however, + // not a programming error: it is a runtime error that can happen in + // optimized mode even if it ran well in debug mode, so it is + // reasonable to try to catch this error also in optimized mode. For + // this case, there is the AssertThrow macro: it checks + // the condition like the Assert macro, but not only in + // debug mode; it then outputs an error message, but instead of + // aborting the program as in the case of the Assert + // macro, the exception is thrown using the throw command + // of C++. This way, one has the possibility to catch this error and + // take reasonable counter actions. One such measure would be to + // refine the grid globally, as the case of insufficient directions + // can not occur if every cell of the initial grid has been refined at + // least once. + AssertThrow(determinant(Y) != 0, ExcInsufficientDirections()); + + // If, on the other hand, the matrix is invertible, then invert it, + // multiply the other quantity with it, and compute the estimated error + // using this quantity and the correct powers of the mesh width: + const Tensor<2, dim> Y_inverse = invert(Y); + + const Tensor<1, dim> gradient = Y_inverse * projected_gradient; + + // The last part of this function is the one where we write into + // the element of the output vector what we have just + // computed. The address of this vector has been stored in the + // scratch data object, and all we have to do is know how to get + // at the correct element inside this vector -- but we can ask the + // cell we're on the how-manyth active cell it is for this: + scratch_data.error_per_cell(cell->active_cell_index()) = + (std::pow(cell->diameter(), 1 + 1.0 * dim / 2) * gradient.norm()); + } +} // namespace Step9 + + +// @sect3{Main function} + +// The main function is similar to the previous examples. The +// primary difference is that we use MultithreadInfo to set the maximum +// number of threads (see the documentation module @ref threads +// "Parallel computing with multiple processors accessing shared memory" +// for more information). The number of threads used is the minimum of the +// environment variable DEAL_II_NUM_THREADS and the parameter of +// set_thread_limit. If no value is given to +// set_thread_limit, the default value from the Intel Threading +// Building Blocks (TBB) library is used. If the call to +// set_thread_limit is omitted, the number of threads will be +// chosen by TBB independently of DEAL_II_NUM_THREADS. +int main() +{ + using namespace dealii; + try + { + // MultithreadInfo::set_thread_limit(); + + Step9::AdvectionProblem<3> advection_problem_2d; + advection_problem_2d.run(); + } + catch (std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} -- 2.39.5