From: wolf Date: Tue, 5 Sep 2000 13:37:40 +0000 (+0000) Subject: Simplify multithreading scheme by putting one data object at each X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=eac962b358b860742cc0ee596e6b70447e13eaa5;p=dealii-svn.git Simplify multithreading scheme by putting one data object at each thread and other small things. git-svn-id: https://svn.dealii.org/trunk@3297 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/numerics/error_estimator.h b/deal.II/deal.II/include/numerics/error_estimator.h index b88fc3ed46..e17a9e4b8d 100644 --- a/deal.II/deal.II/include/numerics/error_estimator.h +++ b/deal.II/deal.II/include/numerics/error_estimator.h @@ -296,18 +296,20 @@ class KellyErrorEstimator /** - * All data needed by the several functions - * of the error estimator is gathered in - * this struct. It is passed as a reference - * to the seperate functions in this class. + * All data needed by the several + * functions of the error + * estimator is gathered in this + * struct. It is passed as a + * reference to the seperate + * functions in this class. * * The reason for invention of * this object is two-fold: * first, class member data is * not possible because no real * object is created (all - * functions are @p{static}), which - * is a historical + * functions are @p{static}), + * which is a historical * reason. Second, if we don't * collect the data the various * functions need somewhere at a @@ -330,33 +332,46 @@ class KellyErrorEstimator * parallel, since they are quite * often blocked by these * synchronisation points. + * + * Thus, every thread gets an + * instance of this class to work + * with and needs not allocate + * memory itself, or synchronise + * with other threads. */ struct Data { - const DoFHandler &dof; + const DoFHandler &dof_handler; const Quadrature &quadrature; const FunctionMap &neumann_bc; const Vector &solution; - vector component_mask; + const vector component_mask; const Function *coefficients; - unsigned int n_threads; - - DoFHandler::active_cell_iterator endc; - unsigned int n_components; - unsigned int n_q_points; - FaceIntegrals face_integrals; + const unsigned int n_threads; /** - * A vector to store the jump of the - * normal vectors in the quadrature - * points. - * There is one vector for every - * thread used in the estimator. - * The allocation of memory has to - * be global to enable fast use - * of multithreading + * Reference to the global + * object that stores the + * face integrals. + */ + FaceIntegrals &face_integrals; + + /** + * A vector to store the jump + * of the normal vectors in + * the quadrature points + * (i.e. a temporary + * value). This vector is not + * allocated inside the + * functions that use it, but + * rather globally, since + * memory allocation is slow, + * in particular in presence + * of multiple threads where + * synchronisation makes + * things even slower. */ - vector< vector > > phi; + vector > phi; /** * A vector for the gradients of @@ -370,18 +385,18 @@ class KellyErrorEstimator * index the number of the * quadrature point. */ - vector< vector > > > psi; + vector > > psi; /** * The same vector for a neighbor cell */ - vector< vector > > > neighbor_psi; + vector > > neighbor_psi; /** * The normal vectors of the finite * element function on one face */ - vector< vector > > normal_vectors; + vector > normal_vectors; /** * Two arrays needed for the @@ -389,8 +404,8 @@ class KellyErrorEstimator * the jumps, if they are * given. */ - vector< vector > coefficient_values1; - vector< vector > > coefficient_values; + vector coefficient_values1; + vector > coefficient_values; /** * Array for the products of @@ -398,7 +413,7 @@ class KellyErrorEstimator * weights of quadraturs * points. */ - vector< vector > JxW_values; + vector JxW_values; /** * A constructor of the @@ -409,9 +424,10 @@ class KellyErrorEstimator const Quadrature &quadrature, const FunctionMap &neumann_bc, const Vector &solution, - vector component_mask_, + const vector &component_mask, const Function *coefficients, - unsigned int n_threads); + const unsigned int n_threads, + FaceIntegrals &face_integrals); }; @@ -434,25 +450,26 @@ class KellyErrorEstimator const unsigned int this_thread); /** - * Actually do the computation on a face - * which has no hanging nodes (it is - * regular), i.e. - * either on the other side there is - * nirvana (face is at boundary), or - * the other side's refinement level - * is the same as that of this side, - * then handle the integration of + * Actually do the computation on + * a face which has no hanging + * nodes (it is regular), i.e. + * either on the other side there + * is nirvana (face is at + * boundary), or the other side's + * refinement level is the same + * as that of this side, then + * handle the integration of * these both cases together. * - * The meaning of the parameters becomes - * clear when looking at the source - * code. This function is only - * externalized from @p{estimate_error} - * to avoid ending up with a function - * of 500 lines of code. + * The meaning of the parameters + * becomes clear when looking at + * the source code. This function + * is only externalized from + * @p{estimate_error} to avoid + * ending up with a function of + * 500 lines of code. */ static void integrate_over_regular_face (Data &data, - const unsigned int this_thread, const active_cell_iterator &cell, const unsigned int face_no, FEFaceValues &fe_face_values_cell, @@ -460,15 +477,16 @@ class KellyErrorEstimator /** - * The same applies as for the function - * above, except that integration is - * over face @p{face_no} of @p{cell}, where - * the respective neighbor is refined, - * so that the integration is a bit more + * The same applies as for the + * function above, except that + * integration is over face + * @p{face_no} of @p{cell}, where + * the respective neighbor is + * refined, so that the + * integration is a bit more * complex. */ static void integrate_over_irregular_face (Data &data, - const unsigned int this_thread, const active_cell_iterator &cell, const unsigned int face_no, FEFaceValues &fe_face_values, diff --git a/deal.II/deal.II/source/numerics/error_estimator.cc b/deal.II/deal.II/source/numerics/error_estimator.cc index 0abd75a8ce..a117fc16d1 100644 --- a/deal.II/deal.II/source/numerics/error_estimator.cc +++ b/deal.II/deal.II/source/numerics/error_estimator.cc @@ -15,6 +15,7 @@ #include #include #include +#include #include #include #include @@ -30,12 +31,6 @@ #include #include - // if multithreaded include - // ThreadManager -#ifdef DEAL_II_USE_MT -# include -#endif - static @@ -55,13 +50,15 @@ KellyErrorEstimator<1>::Data::Data(const DoFHandler<1> &, const Quadrature<0> &, const FunctionMap &, const Vector &, - vector , + const vector &, const Function<1> *, - unsigned int ): - dof(* static_cast *> (0)), + const unsigned int , + FaceIntegrals &): + dof_handler(* static_cast *> (0)), quadrature(* static_cast *> (0)), neumann_bc(* static_cast (0)), - solution(* static_cast *> (0)) + solution(* static_cast *> (0)), + face_integrals (* static_cast (0)) { Assert (false, ExcInternalError()); } @@ -69,26 +66,24 @@ KellyErrorEstimator<1>::Data::Data(const DoFHandler<1> &, #else template -KellyErrorEstimator::Data::Data(const DoFHandler &dof, +KellyErrorEstimator::Data::Data(const DoFHandler &dof_handler, const Quadrature &quadrature, const FunctionMap &neumann_bc, const Vector &solution, - vector component_mask_, + const vector &component_mask, const Function *coefficients, - unsigned int n_threads): - dof(dof), - quadrature(quadrature), - neumann_bc(neumann_bc), - solution(solution), - coefficients(coefficients), - n_threads(n_threads) + const unsigned int n_threads, + FaceIntegrals &face_integrals): + dof_handler (dof_handler), + quadrature (quadrature), + neumann_bc (neumann_bc), + solution (solution), + component_mask (component_mask), + coefficients (coefficients), + n_threads (n_threads), + face_integrals (face_integrals) { - n_components = dof.get_fe().n_components(); - - // if no mask given: treat all components - component_mask = ((component_mask_.size() == 0) ? - vector(n_components, true) : - component_mask_); + const unsigned int n_components = dof_handler.get_fe().n_components(); Assert (component_mask.size() == n_components, ExcInvalidComponentMask()); Assert (count(component_mask.begin(), component_mask.end(), true) > 0, @@ -105,39 +100,25 @@ KellyErrorEstimator::Data::Data(const DoFHandler &dof, for (typename FunctionMap::const_iterator i=neumann_bc.begin(); i!=neumann_bc.end(); ++i) Assert (i->second->n_components == n_components, ExcInvalidBoundaryFunction()); - // the last cell, often needed - endc=dof.end(); - - n_q_points=quadrature.n_quadrature_points; // Init the size of a lot of vectors // needed in the calculations once // per thread. - phi.resize(n_threads); - psi.resize(n_threads); - neighbor_psi.resize(n_threads); - normal_vectors.resize(n_threads); - coefficient_values1.resize(n_threads); - coefficient_values.resize(n_threads); - JxW_values.resize(n_threads); + const unsigned int n_q_points = quadrature.n_quadrature_points; + phi.resize(n_q_points); + psi.resize(n_q_points); + neighbor_psi.resize(n_q_points); + normal_vectors.resize(n_q_points); + coefficient_values1.resize(n_q_points); + coefficient_values.resize(n_q_points); + JxW_values.resize(n_q_points); - for (unsigned int t=0;t::Data::Data(const DoFHandler &dof, template <> void KellyErrorEstimator<1>::estimate_some (Data &, const unsigned int) { + // in 1d, the @p{estimate} function + // does all the work Assert (false, ExcInternalError() ); } + template <> -void KellyErrorEstimator<1>::estimate (const DoFHandler<1> &dof, +void KellyErrorEstimator<1>::estimate (const DoFHandler<1> &dof_handler, const Quadrature<0> &, const FunctionMap &neumann_bc, const Vector &solution, @@ -164,7 +148,7 @@ void KellyErrorEstimator<1>::estimate (const DoFHandler<1> &dof, const Function<1> *coefficient, const unsigned int) { - const unsigned int n_components = dof.get_fe().n_components(); + const unsigned int n_components = dof_handler.get_fe().n_components(); // if no mask given: treat all components vector component_mask ((component_mask_.size() == 0) ? @@ -187,7 +171,7 @@ void KellyErrorEstimator<1>::estimate (const DoFHandler<1> &dof, // reserve one slot for each cell and set // it to zero - error.reinit (dof.get_tria().n_active_cells()); + error.reinit (dof_handler.get_tria().n_active_cells()); // fields to get the gradients on // the present and the neighbor cell. @@ -216,9 +200,9 @@ void KellyErrorEstimator<1>::estimate (const DoFHandler<1> &dof, // two contributions from the two // vertices of each cell. QTrapez<1> quadrature; - FEValues<1> fe_values (dof.get_fe(), quadrature, update_gradients); - DoFHandler<1>::active_cell_iterator cell = dof.begin_active(); - for (unsigned int cell_index=0; cell != dof.end(); ++cell, ++cell_index) + FEValues<1> fe_values (dof_handler.get_fe(), quadrature, update_gradients); + active_cell_iterator cell = dof_handler.begin_active(); + for (unsigned int cell_index=0; cell != dof_handler.end(); ++cell, ++cell_index) { error(cell_index) = 0; // loop over the two points bounding @@ -321,7 +305,7 @@ void KellyErrorEstimator::estimate_some (Data &data, // need not compute all values on the // neighbor cells, so using two objects // gives us a performance gain). - FEFaceValues fe_face_values_cell (data.dof.get_fe(), + FEFaceValues fe_face_values_cell (data.dof_handler.get_fe(), data.quadrature, UpdateFlags(update_gradients | update_JxW_values | @@ -329,15 +313,15 @@ void KellyErrorEstimator::estimate_some (Data &data, (data.coefficients != 0)) ? update_q_points : 0) | update_normal_vectors)); - FEFaceValues fe_face_values_neighbor (data.dof.get_fe(), + FEFaceValues fe_face_values_neighbor (data.dof_handler.get_fe(), data.quadrature, update_gradients); - FESubfaceValues fe_subface_values (data.dof.get_fe(), + FESubfaceValues fe_subface_values (data.dof_handler.get_fe(), data.quadrature, update_gradients); - DoFHandler::active_cell_iterator cell=data.dof.begin_active(); + active_cell_iterator cell=data.dof_handler.begin_active(); // calculate the start cell for // this thread. note that this way @@ -357,12 +341,13 @@ void KellyErrorEstimator::estimate_some (Data &data, // pseudorandom distribution of the // `hard' cells to the different // threads. - for (unsigned int t=0; (t::estimate_some (Data &data, // of faces with contribution zero. const unsigned char boundary_indicator = cell->face(face_no)->boundary_indicator(); - if ((boundary_indicator != 255) && + if (cell->face(face_no)->at_boundary() + && data.neumann_bc.find(boundary_indicator)==data.neumann_bc.end()) { data.face_integrals[cell->face(face_no)] = 0; @@ -409,7 +395,6 @@ void KellyErrorEstimator::estimate_some (Data &data, // then handle the integration of // these both cases together integrate_over_regular_face (data, - this_thread, cell, face_no, fe_face_values_cell, fe_face_values_neighbor); @@ -420,21 +405,25 @@ void KellyErrorEstimator::estimate_some (Data &data, // not fit into the framework of // the above function integrate_over_irregular_face (data, - this_thread,cell, face_no, + cell, face_no, fe_face_values_cell, fe_subface_values); }; - // next cell in this thread - for (unsigned int t=0;((t -void KellyErrorEstimator::estimate (const DoFHandler &dof, +void KellyErrorEstimator::estimate (const DoFHandler &dof_handler, const Quadrature &quadrature, const FunctionMap &neumann_bc, const Vector &solution, @@ -447,82 +436,96 @@ void KellyErrorEstimator::estimate (const DoFHandler &dof, #ifndef DEAL_II_USE_MT n_threads = 1; #endif + Assert (n_threads > 0, ExcInternalError()); - // all the data needed in the error- - // estimator is gathered in this stuct. - KellyErrorEstimator::Data data (dof, - quadrature, - neumann_bc, - solution, - component_mask, - coefficients, - n_threads); - // map of integrals indexed by + // Map of integrals indexed by // the corresponding face. In this map // we store the integrated jump of the - // gradient for each face. By doing so, - // we can check whether we have already - // integrated along this face by testing - // whether the respective face is already - // a key in this map. + // gradient for each face. // At the end of the function, we again // loop over the cells and collect the - // conrtibutions of the different faces + // contributions of the different faces // of the cell. - // the values for all faces are set to - // -10e20. It would cost a lot of time - // to synchronise the initialisation - // of the map in multithreaded mode. - // negative value indicates that the - // face is not calculated. - - for (active_cell_iterator cell=data.dof.begin_active(); cell!=data.endc; ++cell) + // + // The initial values for all faces + // are set to -10e20. It would cost + // a lot of time to synchronise the + // initialisation (i.e. the + // creation of new keys) of the map + // in multithreaded mode. Negative + // value indicates that the face + // has not yet been processed. + FaceIntegrals face_integrals; + for (active_cell_iterator cell=dof_handler.begin_active(); cell!=dof_handler.end(); ++cell) for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) - data.face_integrals[cell->face(face_no)]=-10e20; - + face_integrals[cell->face(face_no)]=-10e20; - // split all cells into threads - // if multithreading is used -#ifdef DEAL_II_USE_MT + // all the data needed in the error + // estimator by each of the threads + // is gathered in the following + // stuctures + // + // note that if no component mask + // was given, then treat all + // components + vector data_structures (n_threads); + for (unsigned int i=0; i(dof_handler.get_fe().n_components(), + true) : + component_mask), + coefficients, + n_threads, + face_integrals); + + // split all cells into threads if + // multithreading is used and run + // the whole thing Threads::ThreadManager thread_manager; - for (unsigned int i=0;i::estimate_some) - .collect_args (data, i)); + .collect_args (*data_structures[i], i)); thread_manager.wait(); + + // delete the structures for the + // different threads again. the + // results are in the + // face_integrals object + for (unsigned int i=0; i::faces_per_cell; ++face_no) { - Assert(data.face_integrals.find(cell->face(face_no)) != data.face_integrals.end(), + Assert(face_integrals.find(cell->face(face_no)) != face_integrals.end(), ExcInternalError()); - error(present_cell) += (data.face_integrals[cell->face(face_no)] * + error(present_cell) += (face_integrals[cell->face(face_no)] * cell->diameter() / 24); }; error(present_cell) = sqrt(error(present_cell)); @@ -536,7 +539,6 @@ void KellyErrorEstimator::estimate (const DoFHandler &dof, template <> void KellyErrorEstimator<1>::integrate_over_regular_face (Data &, - const unsigned int, const active_cell_iterator &, const unsigned int , FEFaceValues<1> &, @@ -549,7 +551,6 @@ void KellyErrorEstimator<1>::integrate_over_regular_face (Data &, template <> void KellyErrorEstimator<1>:: integrate_over_irregular_face (Data &, - const unsigned int, const active_cell_iterator &, const unsigned int , FEFaceValues<1> &, @@ -564,13 +565,15 @@ integrate_over_irregular_face (Data &, template void KellyErrorEstimator:: integrate_over_regular_face (Data &data, - const unsigned int this_thread, const active_cell_iterator &cell, const unsigned int face_no, FEFaceValues &fe_face_values_cell, FEFaceValues &fe_face_values_neighbor) { const DoFHandler::face_iterator face = cell->face(face_no); + const unsigned int n_q_points = data.quadrature.n_quadrature_points, + n_components = data.dof_handler.get_fe().n_components(); + // initialize data of the restriction // of this cell to the present face @@ -578,7 +581,7 @@ integrate_over_regular_face (Data &data, // get gradients of the finite element // function on this cell - fe_face_values_cell.get_function_grads (data.solution, data.psi[this_thread]); + fe_face_values_cell.get_function_grads (data.solution, data.psi); // now compute over the other side of // the face @@ -589,7 +592,7 @@ integrate_over_regular_face (Data &data, Assert (cell->neighbor(face_no).state() == valid, ExcInternalError()); - const DoFHandler::active_cell_iterator neighbor = cell->neighbor(face_no); + const active_cell_iterator neighbor = cell->neighbor(face_no); // find which number the current // face has relative to the neighboring @@ -604,12 +607,12 @@ integrate_over_regular_face (Data &data, // get gradients on neighbor cell fe_face_values_neighbor.get_function_grads (data.solution, - data.neighbor_psi[this_thread]); + data.neighbor_psi); // compute the jump in the gradients - for (unsigned int component=0; componentvalue_list (fe_face_values_cell.get_quadrature_points(), - data.coefficient_values1[this_thread]); - for (unsigned int component=0; componentvector_value_list (fe_face_values_cell.get_quadrature_points(), - data.coefficient_values[this_thread]); - for (unsigned int component=0; component > g(data.n_q_points, Vector(data.n_components)); + vector > g(n_q_points, Vector(n_components)); data.neumann_bc.find(boundary_indicator)->second ->vector_value_list (fe_face_values_cell.get_quadrature_points(), g); - for (unsigned int component=0; component void KellyErrorEstimator:: integrate_over_irregular_face (Data &data, - const unsigned int this_thread, const active_cell_iterator &cell, const unsigned int face_no, FEFaceValues &fe_face_values, FESubfaceValues &fe_subface_values) { const DoFHandler::cell_iterator neighbor = cell->neighbor(face_no); + const unsigned int n_q_points = data.quadrature.n_quadrature_points, + n_components = data.dof_handler.get_fe().n_components(); Assert (neighbor.state() == valid, ExcInternalError()); Assert (neighbor->has_children(), ExcInternalError()); @@ -764,20 +768,20 @@ integrate_over_irregular_face (Data &data, // store the gradient of the solution // in psi fe_subface_values.reinit (cell, face_no, subface_no); - fe_subface_values.get_function_grads (data.solution, data.psi[this_thread]); + fe_subface_values.get_function_grads (data.solution, data.psi); // restrict the finite element on the // neighbor cell to the common @p{subface}. // store the gradient in @p{neighbor_psi} fe_face_values.reinit (neighbor_child, neighbor_neighbor); - fe_face_values.get_function_grads (data.solution, data.neighbor_psi[this_thread]); + fe_face_values.get_function_grads (data.solution, data.neighbor_psi); // compute the jump in the gradients - for (unsigned int component=0; componentn_components == 1) { data.coefficients->value_list (fe_face_values.get_quadrature_points(), - data.coefficient_values1[this_thread]); - for (unsigned int component=0; componentvector_value_list (fe_face_values.get_quadrature_points(), - data.coefficient_values[this_thread]); - for (unsigned int component=0; componentface(neighbor_neighbor)] = face_integral; };