From: Wolfgang Bangerth Date: Fri, 22 Nov 2013 15:05:35 +0000 (+0000) Subject: First part of documenting the switch to WorkStream in GradientEstimation. X-Git-Tag: v8.1.0~205 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4a128c19a7926914704ac5b89f91aa302772ab49;p=dealii.git First part of documenting the switch to WorkStream in GradientEstimation. git-svn-id: https://svn.dealii.org/trunk@31762 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-9/step-9.cc b/deal.II/examples/step-9/step-9.cc index 06b54d3c9b..c10edf29eb 100644 --- a/deal.II/examples/step-9/step-9.cc +++ b/deal.II/examples/step-9/step-9.cc @@ -77,7 +77,7 @@ namespace Step9 // @sect3{AdvectionProblem class declaration} // Following we declare the main class of this program. It is very much - // alike the main classes of previous examples, so we again only comment on + // like the main classes of previous examples, so we again only comment on // the differences. template class AdvectionProblem @@ -144,8 +144,8 @@ namespace Step9 void assemble_system (); void local_assemble_system (const typename DoFHandler::active_cell_iterator &cell, - AssemblyScratchData &scratch, - AssemblyCopyData ©_data); + AssemblyScratchData &scratch, + AssemblyCopyData ©_data); void copy_local_to_global (const AssemblyCopyData ©_data); @@ -407,21 +407,15 @@ namespace Step9 // 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 if higher + // 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 one private function - // that does the actual work on an interval of all active cells. The latter - // is called by the first one in order to be able to do the computations in - // parallel if your computer has more than one processor. While the first - // function accepts as parameter a vector into which the error indicator is - // written for each cell. This vector is passed on to the second function - // that actually computes the error indicators on some cells, and the - // respective elements of the vector are written. By the way, we made it - // somewhat of a 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. + // 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 @@ -430,26 +424,66 @@ namespace Step9 // more common case of invalid parameters to a function, namely a vector of // wrong size. // - // Two annotations to this class are still in order: the first is that 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, i.e. visible to the - // outside world but not callable. 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. + // 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 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 DoF - // handler object that one passes as first argument. + // 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 + // DoF handler object that one passes as first argument. // - // Finally note that the IndexInterval typedef is introduced as - // a convenient abbreviation for an otherwise lengthy type name. + // 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 (i) classes for scratch and + // copy objects, (ii) a function that does the local computation on one + // cell, and (iii) a function that copies the local result into a global + // object. 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. 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 an 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. + // + // The second idea to make this approach work is this: If we want to write + // the result into its final destination right away, then the local worker + // function needs to already know where this destination is. Here, this is + // an element of a vector -- but which element is something that the local + // worker function (or, if we wanted to use one, a copy-local-to-global + // function) can not determine easily just knowing an iterator to a cell it + // is supposed to work on. Consequently, in addition to a cell, we need to + // pass a second piece of identifying information along: the element of the + // output vector to write into. What this means is that the work items are + // identified by two iterators: to a cell, and to an output vector + // element. Moving from one work item to the next requires incrementing both + // iterators. deal.II has a class for this, called SynchronousIterators, + // that takes a tuple of iterator types as arguments and stores an iterator + // of each type. Whenever the SynchronousIterators object is incremented, it + // increments the stored iterators in turn. Thus, this class is exactly what + // we need to do our work, and we consequently use it as the first argument + // of the worker function. We will further down below show how to create + // such an object. class GradientEstimation { public: @@ -466,7 +500,7 @@ namespace Step9 private: template - struct EstimateScratchData + struct EstimateScratchData { EstimateScratchData (const FiniteElement &fe, const Vector &solution); @@ -476,22 +510,15 @@ namespace Step9 Vector solution; }; - // There is nothing to copy but WorkStream requires a CopyData structure - template - struct EstimateCopyData - { - EstimateCopyData () {} - }; + struct EstimateCopyData + {}; template - static void estimate_cell ( - const SynchronousIterators::active_cell_iterator, - Vector::iterator> > &cell, - EstimateScratchData &scratch_data, - const EstimateCopyData ©_data); - // There is nothing to copy but WorkStream required a copy function - template - static void dummy_copy(const EstimateCopyData ©_data) {} + static + void estimate_cell (const SynchronousIterators::active_cell_iterator, + Vector::iterator> > &cell, + EstimateScratchData &scratch_data, + const EstimateCopyData ©_data); }; @@ -503,7 +530,8 @@ namespace Step9 // the function setup_system follow the same pattern that was // used previously, so we need not comment on these three function: template - AdvectionProblem::AdvectionProblem () : + AdvectionProblem::AdvectionProblem () + : dof_handler (triangulation), fe(1) {} @@ -562,12 +590,12 @@ namespace Step9 void AdvectionProblem::assemble_system () { WorkStream::run(dof_handler.begin_active(), - dof_handler.end(), - *this, - &AdvectionProblem::local_assemble_system, - &AdvectionProblem::copy_local_to_global, - AssemblyScratchData(fe), - AssemblyCopyData()); + dof_handler.end(), + *this, + &AdvectionProblem::local_assemble_system, + &AdvectionProblem::copy_local_to_global, + AssemblyScratchData(fe), + AssemblyCopyData()); // After the matrix has been assembled in parallel, we still have to @@ -601,15 +629,15 @@ namespace Step9 template AdvectionProblem::AssemblyScratchData:: AssemblyScratchData (const FiniteElement &fe) - : - fe_values (fe, - QGauss(2), - update_values | update_gradients | - update_quadrature_points | update_JxW_values), - fe_face_values (fe, - QGauss(2), - update_values | update_quadrature_points | - update_JxW_values | update_normal_vectors) + : + fe_values (fe, + QGauss(2), + update_values | update_gradients | + update_quadrature_points | update_JxW_values), + fe_face_values (fe, + QGauss(2), + update_values | update_quadrature_points | + update_JxW_values | update_normal_vectors) {} @@ -617,15 +645,15 @@ namespace Step9 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) + : + 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) {} @@ -668,8 +696,8 @@ namespace Step9 void AdvectionProblem:: local_assemble_system (const typename DoFHandler::active_cell_iterator &cell, - AssemblyScratchData &scratch_data, - AssemblyCopyData ©_data) + AssemblyScratchData &scratch_data, + AssemblyCopyData ©_data) { // First of all, we will need some objects that describe boundary values, // right hand side function and the advection field. As we will only @@ -821,12 +849,12 @@ namespace Step9 { for (unsigned int i=0; i ::EstimateScratchData (const FiniteElement &fe, const Vector &solution) - : - fe_midpoint_value(fe, - QMidpoint (), - update_values | update_quadrature_points), - solution(solution) + : + fe_midpoint_value(fe, + QMidpoint (), + update_values | update_quadrature_points), + solution(solution) {} @@ -959,11 +987,11 @@ namespace Step9 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) + : + 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) {} @@ -1002,23 +1030,20 @@ namespace Step9 // purpose. Here, in this case, Compaq's cxx compiler choked // on the code so we use this workaround with the function pointer: void (*estimate_cell_ptr) (const SynchronousIterators::active_cell_iterator,Vector::iterator> > &cell, - EstimateScratchData &scratch_data, - const EstimateCopyData ©_data) + typename DoFHandler::active_cell_iterator,Vector::iterator> > &cell, + EstimateScratchData &scratch_data, + const EstimateCopyData ©_data) = &GradientEstimation::template estimate_cell; - void (*dummy_copy) (const EstimateCopyData ©_data) - = &GradientEstimation::template dummy_copy; - typedef std_cxx1x::tuple::active_cell_iterator,Vector::iterator> - Iterators; + Iterators; SynchronousIterators begin_sync_it(Iterators(dof_handler.begin_active(), - error_per_cell.begin())); + error_per_cell.begin())); SynchronousIterators end_sync_it(Iterators(dof_handler.end(),error_per_cell.end())); WorkStream::run(begin_sync_it,end_sync_it, estimate_cell_ptr, - dummy_copy, + std_cxx1x::function (), EstimateScratchData (dof_handler.get_fe(),solution), EstimateCopyData ()); @@ -1049,10 +1074,10 @@ namespace Step9 template void GradientEstimation::estimate_cell (const SynchronousIterators::active_cell_iterator,Vector::iterator> > &cell, - EstimateScratchData &scratch_data, - const EstimateCopyData ©_data) - { + typename DoFHandler::active_cell_iterator,Vector::iterator> > &cell, + EstimateScratchData &scratch_data, + const EstimateCopyData ©_data) + { // We need space for the tensor Y, which is the sum of // outer products of the y-vectors. Tensor<2,dim> Y; @@ -1067,7 +1092,7 @@ namespace Step9 GeometryInfo::max_children_per_face); typename DoFHandler::active_cell_iterator cell_it(std_cxx1x::get<0>(cell.iterators)); - + // First initialize the FEValues object, as well as the // Y tensor: scratch_data.fe_midpoint_value.reinit (cell_it); @@ -1151,7 +1176,7 @@ namespace Step9 // as an internal error. We therefore use a predefined // exception class to throw here. Assert (neighbor_child->neighbor(face_no==0 ? 1 : 0) - ==std_cxx1x::get<0>(cell.iterators),ExcInternalError()); + ==std_cxx1x::get<0>(cell.iterators),ExcInternalError()); // If the check succeeded, we push the active neighbor // we just found to the stack we keep: @@ -1204,7 +1229,7 @@ namespace Step9 const Point neighbor_center = scratch_data.fe_midpoint_value.quadrature_point(0); scratch_data.fe_midpoint_value.get_function_values (scratch_data.solution, - neighbor_midpoint_value); + neighbor_midpoint_value); // Compute the vector y connecting the centers of the // two cells. Note that as opposed to the introduction, we denote @@ -1260,9 +1285,9 @@ namespace Step9 contract (gradient, Y_inverse, projected_gradient); *(std_cxx1x::get<1>(cell.iterators)) = (std::pow(std_cxx1x::get<0>(cell.iterators)->diameter(), - 1+1.0*dim/2) * - std::sqrt(gradient.square())); - + 1+1.0*dim/2) * + std::sqrt(gradient.square())); + } }