From: Martin Kronbichler Date: Mon, 12 Jan 2009 16:38:50 +0000 (+0000) Subject: Provide a function for generating nice output of computing times. Use that in step-32. X-Git-Tag: v8.0.0~8143 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=ef55eb92cdc6b770f328c770d146ffed383f970b;p=dealii.git Provide a function for generating nice output of computing times. Use that in step-32. git-svn-id: https://svn.dealii.org/trunk@18185 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/include/base/timer.h b/deal.II/base/include/base/timer.h index 28cf778628..a42f49fdb0 100644 --- a/deal.II/base/include/base/timer.h +++ b/deal.II/base/include/base/timer.h @@ -14,6 +14,9 @@ #define __deal2__timer_h #include +#include +#include +#include DEAL_II_NAMESPACE_OPEN @@ -167,6 +170,80 @@ class Timer bool running; }; + + +/** + * This class can be used to generate formatted output from time + * measurements of different subsections in a program. It is possible to + * create several sections that perform certain aspects of the program. A + * section can be entered several times. By changing the options in + * OutputFrequency and OutputType, the user can choose whether output should + * be generated every time a section is joined or just in the end of the + * program. Moreover, it is possible to show CPU times, wall times or both. + * + * TODO: Write a more extensive documentation. + * + * @ingroup utilities + * @author M. Kronbichler, 2009. + */ +class TimerOutput +{ + public: + // Sets whether to generate output every + // time we exit a section, just in the + // end, or both. + enum OutputFrequency {every_call, summary, every_call_and_summary} + output_frequency; + + // Sets whether to show CPU times, wall + // times, or both CPU and wall times. + enum OutputType {cpu_times, wall_times, cpu_and_wall_times} + output_type; + + // Constructor that takes std::cout as + // output stream. + TimerOutput (std::ostream &stream, + const enum OutputFrequency output_frequency, + const enum OutputType output_type); + + // Constructor that takes a + // ConditionalOStream to write output to. + TimerOutput (ConditionalOStream &stream, + const enum OutputFrequency output_frequency, + const enum OutputType output_type); + + // Destructor. Calls print_summary() in + // case the option for writing the + // summary output is set. + ~TimerOutput(); + + // Open a section by given a string name + // of it. In case the name already + // exists, that section is done once + // again. + void enter_section (const std::string §ion_name); + + // Leave a section. If no name is given, + // the last section that was entered is + // left. + void exit_section (const std::string §ion_name = std::string()); + + // Print a formatted table that + // summarizes the time consumed in the + // various sections. + void print_summary (); + + private: + Timer timer_all; + std::vector section_timers; + std::vector section_names; + std::vector section_total_cpu_times; + std::vector section_total_wall_times; + std::vector section_n_calls; + ConditionalOStream out_stream; + std::vector active_sections; +}; + DEAL_II_NAMESPACE_CLOSE #endif diff --git a/deal.II/base/source/timer.cc b/deal.II/base/source/timer.cc index eda6649283..f20c70dfff 100644 --- a/deal.II/base/source/timer.cc +++ b/deal.II/base/source/timer.cc @@ -13,13 +13,18 @@ #include - +#include +#include // these includes should probably be properly // ./configure'd using the AC_HEADER_TIME macro: #include #include +#ifdef DEAL_II_COMPILER_SUPPORTS_MPI +#include +#endif + // on SunOS 4.x, getrusage is stated in the man pages and exists, but // is not declared in resource.h. declare it ourselves @@ -131,4 +136,347 @@ void Timer::reset () running = false; } + + +/* ---------------------------- TimerOutput -------------------------- */ + +TimerOutput::TimerOutput (std::ostream &stream, + const enum OutputFrequency output_frequency, + const enum OutputType output_type) + : + output_frequency (output_frequency), + output_type (output_type), + out_stream (stream, true) +{} + + + +TimerOutput::TimerOutput (ConditionalOStream &stream, + const enum OutputFrequency output_frequency, + const enum OutputType output_type) + : + output_frequency (output_frequency), + output_type (output_type), + out_stream (stream) +{} + + + +TimerOutput::~TimerOutput() +{ + while (active_sections.size() > 0) + exit_section(); + + if (output_frequency != every_call) + print_summary(); +} + + + +void +TimerOutput::enter_section (const std::string §ion_name) +{ + Assert (section_name.empty() == false, + ExcMessage ("Section string is empty.")); + + unsigned int this_section_number = numbers::invalid_unsigned_int; + + // check whether the requested section + // already exists + for (unsigned int i=0; i::iterator position_to_delete + = active_sections.begin() + active_index_to_delete; + active_sections.erase(position_to_delete); +} + + + +void +TimerOutput::print_summary () +{ + // in case we want to write CPU times + if (output_type != wall_times) + { + double total_cpu_time; + +#ifdef DEAL_II_COMPILER_SUPPORTS_MPI + { + double my_cpu_time = timer_all(); + int mpiInitialized; + MPI_Initialized(&mpiInitialized); + + if( mpiInitialized ) + { + MPI_Allreduce (&my_cpu_time, &total_cpu_time, 1, MPI_DOUBLE, MPI_SUM, + MPI_COMM_WORLD); + } + else + total_cpu_time = my_cpu_time; + } +#else + total_cpu_time = timer_all(); +#endif + + // check that the sum of all times is + // less or equal than the total + // time. otherwise, we might have + // generated a lot of overhead in this + // function. + double check_time = 0.; + for (unsigned int i=0; i total_cpu_time) + { + total_cpu_time = check_time; + out_stream << std::endl << "Sum of counted times is larger than total time. " + << "Timer function may have introduced too much overhead." << std::endl; + } + + // now generate a nice table + out_stream << "\n\n" + << "+---------------------------------------------+------------" + << "+------------+\n" + << "| Total CPU time elapsed since start |"; + std::cout.width(10); + std::cout.precision(3); + out_stream << total_cpu_time << "s | |\n"; + out_stream << "| | " + << "| |\n"; + out_stream << "| Section | no. calls |"; + std::cout.width(10); + std::cout.precision(3); + out_stream << " CPU time " << " | % of total |\n"; + out_stream << "+---------------------------------------------+------------" + << "+------------+"; + for (unsigned int i=0; i total_wall_time) + { + total_wall_time = check_time; + out_stream << std::endl + << "Sum of counted times is larger than total time. " + << "Timer function may have introduced too much overhead." + << std::endl; + } + + // now generate a nice table + out_stream << "\n\n" + << "+---------------------------------------------+------------" + << "+------------+\n" + << "| Total wallclock time elapsed from start |"; + std::cout.width(10); + std::cout.precision(3); + out_stream << total_wall_time << "s | |\n"; + out_stream << "| | " + << "| |\n"; + out_stream << "| Section | no. calls |"; + std::cout.width(10); + std::cout.precision(3); + out_stream << " CPU time " << " | % of total |\n"; + out_stream << "+---------------------------------------------+------------" + << "+------------+"; + for (unsigned int i=0; i + // Time measurements. +#include + #include #include #include @@ -65,7 +68,6 @@ // names into global namespace using namespace dealii; - // @sect3{Equation data} @@ -95,10 +97,11 @@ namespace EquationData template double - TemperatureInitialValues::value (const Point &p, + TemperatureInitialValues::value (const Point &, const unsigned int) const { - return (p.norm() < 0.55+0.02*std::sin(p[0]*20) ? 1 : 0); + //return (p.norm() < 0.55+0.02*std::sin(p[0]*20) ? 1 : 0); + return 0.; } @@ -130,10 +133,32 @@ namespace EquationData template double - TemperatureRightHandSide::value (const Point &, - const unsigned int /*component*/) const + TemperatureRightHandSide::value (const Point &p, + const unsigned int component) const { - return 0; + // return 0; + Assert (component == 0, + ExcMessage ("Invalid operation for a scalar function.")); + + Assert ((dim==2) || (dim==3), ExcNotImplemented()); + + static const Point source_centers[3] + = { (dim == 2 ? Point(.3,.1) : Point(.3,.5,.1)), + (dim == 2 ? Point(.45,.1) : Point(.45,.5,.1)), + (dim == 2 ? Point(.75,.1) : Point(.75,.5,.1)) }; + static const double source_radius + = (dim == 2 ? 1./32 : 1./8); + + return ((source_centers[0].distance (p) < source_radius) + || + (source_centers[1].distance (p) < source_radius) + || + (source_centers[2].distance (p) < source_radius) + ? + 1 + : + 0); + } @@ -278,7 +303,6 @@ class BoussinesqFlowProblem void output_results () const; void refine_mesh (const unsigned int max_grid_level); - static double compute_viscosity(const std::vector &old_temperature, const std::vector &old_old_temperature, @@ -286,13 +310,12 @@ class BoussinesqFlowProblem const std::vector > &old_old_temperature_grads, const std::vector &old_temperature_laplacians, const std::vector &old_old_temperature_laplacians, - const std::vector > &present_velocity_values, + const std::vector > &old_velocity_values, + const std::vector > &old_old_velocity_values, const std::vector &gamma_values, const double global_u_infty, const double global_T_variation, - const double global_Omega_diameter, - const double cell_diameter, - const double old_time_step); + const double cell_diameter); const Epetra_Comm &trilinos_communicator; @@ -300,6 +323,7 @@ class BoussinesqFlowProblem ConditionalOStream pcout; Triangulation triangulation; + double global_Omega_diameter; const unsigned int stokes_degree; FESystem stokes_fe; @@ -311,6 +335,7 @@ class BoussinesqFlowProblem TrilinosWrappers::BlockSparseMatrix stokes_preconditioner_matrix; TrilinosWrappers::MPI::BlockVector stokes_solution; + TrilinosWrappers::BlockVector old_stokes_solution; TrilinosWrappers::MPI::BlockVector stokes_rhs; @@ -334,12 +359,16 @@ class BoussinesqFlowProblem double old_time_step; unsigned int timestep_number; - boost::shared_ptr Amg_preconditioner; - boost::shared_ptr Mp_preconditioner; + boost::shared_ptr Amg_preconditioner; + boost::shared_ptr Mp_preconditioner; + boost::shared_ptr T_preconditioner; bool rebuild_stokes_matrix; - bool rebuild_temperature_matrices; bool rebuild_stokes_preconditioner; + bool rebuild_temperature_matrices; + bool rebuild_temperature_preconditioner; + + TimerOutput computing_timer; }; @@ -369,8 +398,11 @@ BoussinesqFlowProblem::BoussinesqFlowProblem () old_time_step (0), timestep_number (0), rebuild_stokes_matrix (true), + rebuild_stokes_preconditioner (true), rebuild_temperature_matrices (true), - rebuild_stokes_preconditioner (true) + rebuild_temperature_preconditioner (true), + computing_timer (pcout, TimerOutput::summary, + TimerOutput::wall_times) {} @@ -513,35 +545,34 @@ BoussinesqFlowProblem::get_extrapolated_temperature_range () const template double BoussinesqFlowProblem:: -compute_viscosity(const std::vector &old_temperature, - const std::vector &old_old_temperature, - const std::vector > &old_temperature_grads, - const std::vector > &old_old_temperature_grads, - const std::vector &old_temperature_laplacians, - const std::vector &old_old_temperature_laplacians, - const std::vector > &present_velocity_values, - const std::vector &gamma_values, - const double global_u_infty, - const double global_T_variation, - const double global_Omega_diameter, - const double cell_diameter, - const double old_time_step) +compute_viscosity (const std::vector &old_temperature, + const std::vector &old_old_temperature, + const std::vector > &old_temperature_grads, + const std::vector > &old_old_temperature_grads, + const std::vector &old_temperature_laplacians, + const std::vector &old_old_temperature_laplacians, + const std::vector > &old_velocity_values, + const std::vector > &old_old_velocity_values, + const std::vector &gamma_values, + const double global_u_infty, + const double global_T_variation, + const double cell_diameter) { - const double beta = 0.015 * dim; - const double alpha = 1; + const double beta = 0.04 * dim; + const double alpha = 2; if (global_u_infty == 0) return 5e-3 * cell_diameter; const unsigned int n_q_points = old_temperature.size(); - // Stage 1: calculate residual double max_residual = 0; double max_velocity = 0; for (unsigned int q=0; q < n_q_points; ++q) { - const Tensor<1,dim> u = present_velocity_values[q]; + const Tensor<1,dim> u = (old_velocity_values[q] + + old_old_velocity_values[q]) / 2; const double dT_dt = (old_temperature[q] - old_old_temperature[q]) / old_time_step; @@ -567,15 +598,18 @@ compute_viscosity(const std::vector &old_temperature, return (beta * max_velocity * std::min (cell_diameter, - std::pow(cell_diameter,alpha) * max_residual / global_scaling)); + std::pow(cell_diameter,alpha) * + max_residual / global_scaling)); } + // @sect4{BoussinesqFlowProblem::setup_dofs} template void BoussinesqFlowProblem::setup_dofs () { + computing_timer.enter_section("Setup dof systems"); std::vector stokes_sub_blocks (dim+1,0); stokes_sub_blocks[dim] = 1; @@ -664,6 +698,13 @@ void BoussinesqFlowProblem::setup_dofs () sp.compress(); stokes_matrix.reinit (sp); + /*std::cout << "Processor " << trilinos_communicator.MyPID() + << " stokes(0,0) rows: " + << stokes_matrix.block(0,0).matrix->NumMyRows() + << ", nnz: " + << stokes_matrix.block(0,0).matrix->NumMyNonzeros() + << std::endl;*/ + } { @@ -696,6 +737,7 @@ void BoussinesqFlowProblem::setup_dofs () 0, trilinos_communicator); { + T_preconditioner.reset (); temperature_mass_matrix.clear (); temperature_stiffness_matrix.clear (); temperature_matrix.clear (); @@ -711,12 +753,15 @@ void BoussinesqFlowProblem::setup_dofs () } stokes_solution.reinit (stokes_partitioner); + old_stokes_solution.reinit (stokes_partitioner); stokes_rhs.reinit (stokes_partitioner); temperature_solution.reinit (temperature_partitioner); old_temperature_solution.reinit (temperature_partitioner); old_old_temperature_solution.reinit (temperature_partitioner); temperature_rhs.reinit (temperature_partitioner); + + computing_timer.exit_section(); } @@ -791,6 +836,7 @@ BoussinesqFlowProblem::build_stokes_preconditioner () if (rebuild_stokes_preconditioner == false) return; + computing_timer.enter_section (" Build Stokes preconditioner"); pcout << " Rebuilding Stokes preconditioner..." << std::flush; assemble_stokes_preconditioner (); @@ -806,16 +852,17 @@ BoussinesqFlowProblem::build_stokes_preconditioner () Amg_preconditioner->initialize(stokes_preconditioner_matrix.block(0,0), TrilinosWrappers::PreconditionAMG::AdditionalData - (true, true, 5e-2, null_space, 3, 0, false)); + (true, true, 1e-2, null_space, 3, 0, false)); Mp_preconditioner = boost::shared_ptr (new TrilinosWrappers::PreconditionIC()); Mp_preconditioner->initialize (stokes_preconditioner_matrix.block(1,1)); - pcout << std::endl; - rebuild_stokes_preconditioner = false; + + pcout << std::endl; + computing_timer.exit_section(); } @@ -826,6 +873,8 @@ void BoussinesqFlowProblem::assemble_stokes_system () { pcout << " Assembling..." << std::flush; + computing_timer.enter_section (" Assemble Stokes system"); + if (rebuild_stokes_matrix == true) stokes_matrix=0; @@ -860,7 +909,7 @@ void BoussinesqFlowProblem::assemble_stokes_system () std::vector > grads_phi_u (dofs_per_cell); std::vector div_phi_u (dofs_per_cell); std::vector phi_p (dofs_per_cell); - + const FEValuesExtractors::Vector velocities (0); const FEValuesExtractors::Scalar pressure (dim); @@ -934,6 +983,7 @@ void BoussinesqFlowProblem::assemble_stokes_system () rebuild_stokes_matrix = false; pcout << std::endl; + computing_timer.exit_section(); } @@ -947,7 +997,8 @@ void BoussinesqFlowProblem::assemble_temperature_matrix () { if (rebuild_temperature_matrices == false) return; - + + computing_timer.enter_section (" Assemble temperature matrices"); temperature_mass_matrix = 0; temperature_stiffness_matrix = 0; @@ -1015,6 +1066,9 @@ void BoussinesqFlowProblem::assemble_temperature_matrix () temperature_stiffness_matrix.compress(); rebuild_temperature_matrices = false; + rebuild_temperature_preconditioner = true; + + computing_timer.exit_section(); } @@ -1038,6 +1092,15 @@ void BoussinesqFlowProblem::assemble_temperature_system () temperature_matrix.add (time_step, temperature_stiffness_matrix); } temperature_matrix.compress(); + + if (rebuild_temperature_preconditioner == true) + { + T_preconditioner = boost::shared_ptr + (new TrilinosWrappers::PreconditionIC()); + T_preconditioner->initialize (temperature_matrix); + + rebuild_temperature_preconditioner = false; + } temperature_rhs = 0; @@ -1059,8 +1122,8 @@ void BoussinesqFlowProblem::assemble_temperature_system () std::vector local_dof_indices (dofs_per_cell); - std::vector > present_velocity_values (n_q_points); - + std::vector > old_velocity_values (n_q_points); + std::vector > old_old_velocity_values (n_q_points); std::vector old_temperature_values (n_q_points); std::vector old_old_temperature_values(n_q_points); @@ -1081,7 +1144,10 @@ void BoussinesqFlowProblem::assemble_temperature_system () global_T_range = get_extrapolated_temperature_range(); const double global_Omega_diameter = GridTools::diameter (triangulation); - const TrilinosWrappers::BlockVector localized_stokes_solution (stokes_solution); + const TrilinosWrappers::BlockVector + localized_stokes_solution (stokes_solution); + const TrilinosWrappers::BlockVector + localized_old_stokes_solution (old_stokes_solution); const FEValuesExtractors::Vector velocities (0); @@ -1117,9 +1183,11 @@ void BoussinesqFlowProblem::assemble_temperature_system () temperature_right_hand_side.value_list (temperature_fe_values.get_quadrature_points(), gamma_values); - + stokes_fe_values[velocities].get_function_values (localized_stokes_solution, - present_velocity_values); + old_velocity_values); + stokes_fe_values[velocities].get_function_values (localized_old_stokes_solution, + old_old_velocity_values); const double nu = compute_viscosity (old_temperature_values, @@ -1128,12 +1196,12 @@ void BoussinesqFlowProblem::assemble_temperature_system () old_old_temperature_grads, old_temperature_laplacians, old_old_temperature_laplacians, - present_velocity_values, + old_velocity_values, + old_old_velocity_values, gamma_values, global_u_infty, global_T_range.second - global_T_range.first, - global_Omega_diameter, cell->diameter(), - old_time_step); + cell->diameter()); for (unsigned int q=0; q::assemble_temperature_system () phi_T[k] = temperature_fe_values.shape_value (k, q); } - const double old_T = old_temperature_values[q]; - const double old_old_T = old_old_temperature_values[q]; - - const Tensor<1,dim> old_grad_T = old_temperature_grads[q]; - const Tensor<1,dim> old_old_grad_T = old_old_temperature_grads[q]; - - - const Tensor<1,dim> present_u = present_velocity_values[q]; - - if (use_bdf2_scheme == true) - { - for (unsigned int i=0; i ext_grad_T + = (use_bdf2_scheme ? + (old_temperature_grads[q] * + (1+time_step/old_time_step) + - + old_old_temperature_grads[q] * + time_step / old_time_step) + : + old_temperature_grads[q]); + + const Tensor<1,dim> extrapolated_u + = (use_bdf2_scheme ? + (old_velocity_values[q] * (1+time_step/old_time_step) - + old_old_velocity_values[q] * time_step/old_time_step) + : + old_velocity_values[q]); + + for (unsigned int i=0; iget_dof_indices (local_dof_indices); @@ -1216,6 +1271,7 @@ void BoussinesqFlowProblem::assemble_temperature_system () template void BoussinesqFlowProblem::solve () { + computing_timer.enter_section (" Solve Stokes system"); pcout << " Solving..." << std::endl; { @@ -1245,31 +1301,36 @@ void BoussinesqFlowProblem::solve () stokes_constraints.distribute (localized_stokes_solution); stokes_solution = localized_stokes_solution; } + computing_timer.exit_section(); + + + computing_timer.enter_section (" Assemble temperature rhs"); old_time_step = time_step; - time_step = 1./(1.6*dim*std::sqrt(1.*dim)) / + time_step = 1./(1.9*dim*std::sqrt(1.*dim)) / temperature_degree * GridTools::minimal_cell_diameter(triangulation) / - std::max (get_maximal_velocity(), 1.e-5); + std::max (get_maximal_velocity(), 0.01); pcout << " " << "Time step: " << time_step - << std::endl; - + << std::endl; + temperature_solution = old_temperature_solution; assemble_temperature_system (); - { + computing_timer.exit_section (); + + computing_timer.enter_section (" Solve temperature system"); + + { SolverControl solver_control (temperature_matrix.m(), 1e-8*temperature_rhs.l2_norm()); SolverCG cg (solver_control); - TrilinosWrappers::PreconditionIC preconditioner; - preconditioner.initialize (temperature_matrix); - cg.solve (temperature_matrix, temperature_solution, - temperature_rhs, preconditioner); + temperature_rhs, *T_preconditioner); TrilinosWrappers::Vector localized_temperature_solution (temperature_solution); temperature_constraints.distribute (localized_temperature_solution); @@ -1277,12 +1338,12 @@ void BoussinesqFlowProblem::solve () pcout << " " << solver_control.last_step() - << " CG iterations for temperature." - << std::endl; + << " CG iterations for temperature" << std::endl; + computing_timer.exit_section(); double min_temperature = localized_temperature_solution(0), max_temperature = localized_temperature_solution(0); - for (unsigned int i=0; i (min_temperature, localized_temperature_solution(i)); @@ -1312,60 +1373,60 @@ void BoussinesqFlowProblem::output_results () const Assert (joint_dof_handler.n_dofs() == stokes_dof_handler.n_dofs() + temperature_dof_handler.n_dofs(), ExcInternalError()); - + Vector joint_solution (joint_dof_handler.n_dofs()); TrilinosWrappers::BlockVector localized_stokes_solution (stokes_solution); TrilinosWrappers::Vector localized_temperature_solution (temperature_solution); - { - std::vector local_joint_dof_indices (joint_fe.dofs_per_cell); - std::vector local_stokes_dof_indices (stokes_fe.dofs_per_cell); - std::vector local_temperature_dof_indices (temperature_fe.dofs_per_cell); - - typename DoFHandler::active_cell_iterator - joint_cell = joint_dof_handler.begin_active(), - joint_endc = joint_dof_handler.end(), - stokes_cell = stokes_dof_handler.begin_active(), - temperature_cell = temperature_dof_handler.begin_active(); - for (; joint_cell!=joint_endc; ++joint_cell, ++stokes_cell, ++temperature_cell) + if (Utilities::Trilinos::get_this_mpi_process(trilinos_communicator) == 0) + { + { - joint_cell->get_dof_indices (local_joint_dof_indices); - stokes_cell->get_dof_indices (local_stokes_dof_indices); - temperature_cell->get_dof_indices (local_temperature_dof_indices); - - for (unsigned int i=0; i local_joint_dof_indices (joint_fe.dofs_per_cell); + std::vector local_stokes_dof_indices (stokes_fe.dofs_per_cell); + std::vector local_temperature_dof_indices (temperature_fe.dofs_per_cell); + + typename DoFHandler::active_cell_iterator + joint_cell = joint_dof_handler.begin_active(), + joint_endc = joint_dof_handler.end(), + stokes_cell = stokes_dof_handler.begin_active(), + temperature_cell = temperature_dof_handler.begin_active(); + for (; joint_cell!=joint_endc; ++joint_cell, ++stokes_cell, ++temperature_cell) + { + joint_cell->get_dof_indices (local_joint_dof_indices); + stokes_cell->get_dof_indices (local_stokes_dof_indices); + temperature_cell->get_dof_indices (local_temperature_dof_indices); + + for (unsigned int i=0; i joint_solution_names (dim, "velocity"); - joint_solution_names.push_back ("p"); - joint_solution_names.push_back ("T"); - DataOut data_out; + std::vector joint_solution_names (dim, "velocity"); + joint_solution_names.push_back ("p"); + joint_solution_names.push_back ("T"); + + DataOut data_out; - if (Utilities::Trilinos::get_this_mpi_process(trilinos_communicator) == 0) - { data_out.attach_dof_handler (joint_dof_handler); std::vector @@ -1394,6 +1455,7 @@ void BoussinesqFlowProblem::output_results () const template void BoussinesqFlowProblem::refine_mesh (const unsigned int max_grid_level) { + computing_timer.enter_section ("Refine mesh structure, part 1"); Vector estimated_error_per_cell (triangulation.n_active_cells()); TrilinosWrappers::Vector localized_temperature_solution (temperature_solution); @@ -1413,29 +1475,45 @@ void BoussinesqFlowProblem::refine_mesh (const unsigned int max_grid_level) cell != triangulation.end(); ++cell) cell->clear_refine_flag (); - std::vector x_solution (2); - x_solution[0] = temperature_solution; - x_solution[1] = old_temperature_solution; + std::vector x_temperature (2); + x_temperature[0] = temperature_solution; + x_temperature[1] = old_temperature_solution; + TrilinosWrappers::BlockVector x_stokes = stokes_solution; - SolutionTransfer soltrans(temperature_dof_handler); + SolutionTransfer + temperature_trans(temperature_dof_handler); + SolutionTransfer + stokes_trans(stokes_dof_handler); triangulation.prepare_coarsening_and_refinement(); - soltrans.prepare_for_coarsening_and_refinement(x_solution); + temperature_trans.prepare_for_coarsening_and_refinement(x_temperature); + stokes_trans.prepare_for_coarsening_and_refinement(x_stokes); triangulation.execute_coarsening_and_refinement (); + computing_timer.exit_section(); + setup_dofs (); - std::vector tmp (2); - tmp[0] = temperature_solution; - tmp[1] = temperature_solution; - soltrans.interpolate(x_solution, tmp); + computing_timer.enter_section ("Refine mesh structure, part 2"); + + std::vector tmp (2); + tmp[0].reinit (temperature_solution); + tmp[1].reinit (temperature_solution); + temperature_trans.interpolate(x_temperature, tmp); temperature_solution = tmp[0]; old_temperature_solution = tmp[1]; - rebuild_stokes_matrix = true; - rebuild_temperature_matrices = true; - rebuild_stokes_preconditioner = true; + TrilinosWrappers::BlockVector x_stokes_new = stokes_solution; + stokes_trans.interpolate (x_stokes, x_stokes_new); + stokes_solution = x_stokes_new; + + rebuild_stokes_matrix = true; + rebuild_stokes_preconditioner = true; + rebuild_temperature_matrices = true; + rebuild_temperature_preconditioner = true; + + computing_timer.exit_section(); } @@ -1447,11 +1525,13 @@ void BoussinesqFlowProblem::run () const unsigned int initial_refinement = (dim == 2 ? 3 : 2); const unsigned int n_pre_refinement_steps = (dim == 2 ? 4 : 3); - GridGenerator::half_hyper_shell (triangulation, - Point(), 0.5, 1.0); + //GridGenerator::half_hyper_shell (triangulation, + // Point(), 0.5, 1.0); - static HalfHyperShellBoundary boundary; - triangulation.set_boundary (0, boundary); + //static HyperShellBoundary boundary; + //triangulation.set_boundary (0, boundary); + GridGenerator::hyper_cube (triangulation); + global_Omega_diameter = GridTools::diameter (triangulation); triangulation.refine_global (initial_refinement); @@ -1461,11 +1541,72 @@ void BoussinesqFlowProblem::run () start_time_iteration: - VectorTools::project (temperature_dof_handler, - temperature_constraints, - QGauss(temperature_degree+2), - EquationData::TemperatureInitialValues(), - old_temperature_solution); + //VectorTools::project (temperature_dof_handler, + // temperature_constraints, + // QGauss(temperature_degree+2), + // EquationData::TemperatureInitialValues(), + // old_temperature_solution); + assemble_temperature_matrix (); + // Create right hand side for projection + { + QGauss quadrature(temperature_degree+2); + UpdateFlags update_flags = UpdateFlags(update_values | + update_quadrature_points | + update_JxW_values); + FEValues fe_values (temperature_fe, quadrature, update_flags); + + const unsigned int dofs_per_cell = fe_values.dofs_per_cell, + n_q_points = fe_values.n_quadrature_points; + + std::vector dofs (dofs_per_cell); + Vector cell_vector (dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell = temperature_dof_handler.begin_active(), + endc = temperature_dof_handler.end(); + + std::vector rhs_values(n_q_points); + + TrilinosWrappers::MPI::Vector rhs (temperature_mass_matrix.matrix->RowMap()), + sol (temperature_mass_matrix.matrix->RowMap()); + + for (; cell!=endc; ++cell) + if (cell->subdomain_id() == + Utilities::Trilinos::get_this_mpi_process(trilinos_communicator)) + { + fe_values.reinit(cell); + + const std::vector &weights = fe_values.get_JxW_values (); + EquationData::TemperatureInitialValues().value_list + (fe_values.get_quadrature_points(), rhs_values); + + cell_vector = 0; + for (unsigned int point=0; pointget_dof_indices (dofs); + + temperature_constraints.distribute_local_to_global (cell_vector, + dofs, + rhs); + } + + ReductionControl control(5*rhs.size(), 0., 1e-12, false, false); + GrowingVectorMemory memory; + SolverCG cg(control,memory); + + TrilinosWrappers::PreconditionIC prec; + prec.initialize(temperature_mass_matrix); + // solve + cg.solve (temperature_mass_matrix, sol, rhs, prec); + + // distribute solution + old_temperature_solution = sol; + temperature_constraints.distribute (old_temperature_solution); + } timestep_number = 0; time_step = old_time_step = 0; @@ -1502,10 +1643,12 @@ void BoussinesqFlowProblem::run () time += time_step; ++timestep_number; + old_stokes_solution = stokes_solution; old_old_temperature_solution = old_temperature_solution; old_temperature_solution = temperature_solution; } - while (time <= 100); + while (time <= 0); + } @@ -1519,7 +1662,7 @@ int main (int argc, char *argv[]) Utilities::System::MPI_InitFinalize mpi_initialization(argc, argv); - BoussinesqFlowProblem<2> flow_problem; + BoussinesqFlowProblem<3> flow_problem; flow_problem.run (); } catch (std::exception &exc)