From: bangerth Date: Tue, 19 Aug 2008 03:48:14 +0000 (+0000) Subject: Use two DoFHandlers, matrices, sets of vectors, etc, so that we don't always have... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=68d9977653c511327592e54f338deb7a6b398be2;p=dealii-svn.git Use two DoFHandlers, matrices, sets of vectors, etc, so that we don't always have to rebuild all FEValues data for all components, independent of what we're currently doing. git-svn-id: https://svn.dealii.org/trunk@16593 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc index 410b887042..b1ccf2d3dd 100644 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@ -319,38 +319,50 @@ class BoussinesqFlowProblem private: void setup_dofs (); - void assemble_preconditioner (); - void assemble_system (); - void assemble_rhs_T (); + void assemble_stokes_preconditioner (); + void build_stokes_preconditioner (); + void assemble_stokes_system (); + void assemble_temperature_system (); double get_maximal_velocity () const; double get_maximal_temperature () const; void solve (); void output_results () const; void refine_mesh (const unsigned int max_grid_level); - const unsigned int degree; + Triangulation triangulation; - Triangulation triangulation; - FESystem fe; - DoFHandler dof_handler; + const unsigned int degree; + + FESystem stokes_fe; + DoFHandler stokes_dof_handler; + ConstraintMatrix stokes_constraints; + + BlockSparsityPattern stokes_sparsity_pattern; + BlockSparseMatrix stokes_matrix; + BlockSparsityPattern stokes_preconditioner_sparsity_pattern; + BlockSparseMatrix stokes_preconditioner_matrix; + + BlockVector stokes_solution; + BlockVector stokes_rhs; - ConstraintMatrix hanging_node_constraints; + + FE_Q temperature_fe; + DoFHandler temperature_dof_handler; + ConstraintMatrix temperature_constraints; + + SparsityPattern temperature_sparsity_pattern; + SparseMatrix temperature_matrix; - BlockSparsityPattern sparsity_pattern; - BlockSparseMatrix system_matrix; + Vector temperature_solution; + Vector old_temperature_solution; + Vector old_old_temperature_solution; + Vector temperature_rhs; - BlockSparsityPattern preconditioner_sparsity_pattern; - BlockSparseMatrix preconditioner_matrix; double time_step; double old_time_step; unsigned int timestep_number; - BlockVector solution; - BlockVector old_solution; - BlockVector old_old_solution; - BlockVector system_rhs; - boost::shared_ptr Amg_preconditioner; boost::shared_ptr > Mp_preconditioner; @@ -441,10 +453,10 @@ TemperatureBoundaryValues::value (const Point &p, // @sect4{Initial values} template -class InitialValues : public Function +class TemperatureInitialValues : public Function { public: - InitialValues () : Function(dim+2) {} + TemperatureInitialValues () : Function(1) {} virtual double value (const Point &p, const unsigned int component = 0) const; @@ -456,7 +468,7 @@ class InitialValues : public Function template double -InitialValues::value (const Point &, +TemperatureInitialValues::value (const Point &, const unsigned int) const { return 0; @@ -465,11 +477,11 @@ InitialValues::value (const Point &, template void -InitialValues::vector_value (const Point &p, +TemperatureInitialValues::vector_value (const Point &p, Vector &values) const { for (unsigned int c=0; cn_components; ++c) - values(c) = InitialValues::value (p, c); + values(c) = TemperatureInitialValues::value (p, c); } @@ -488,10 +500,10 @@ InitialValues::vector_value (const Point &p, // in the problem description in the // introduction. template -class RightHandSide : public Function +class TemperatureRightHandSide : public Function { public: - RightHandSide () : Function(dim+2) {} + TemperatureRightHandSide () : Function(1) {} virtual double value (const Point &p, const unsigned int component = 0) const; @@ -503,42 +515,35 @@ class RightHandSide : public Function template double -RightHandSide::value (const Point &p, - const unsigned int component) const +TemperatureRightHandSide::value (const Point &p, + const unsigned int /*component*/) const { - // the right hand side only acts on the - // temperature component - if (component == dim+1) - { - 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); + 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); - } - else - return 0; + return ((source_centers[0].distance (p) < source_radius) + || + (source_centers[1].distance (p) < source_radius) + || + (source_centers[2].distance (p) < source_radius) + ? + 1 + : + 0); } template void -RightHandSide::vector_value (const Point &p, - Vector &values) const +TemperatureRightHandSide::vector_value (const Point &p, + Vector &values) const { for (unsigned int c=0; cn_components; ++c) - values(c) = RightHandSide::value (p, c); + values(c) = TemperatureRightHandSide::value (p, c); } @@ -726,7 +731,7 @@ class BlockSchurPreconditioner : public Subscriptor const BlockVector &src) const; private: - const SmartPointer > system_matrix; + const SmartPointer > stokes_matrix; const SmartPointer, PreconditionerMp > > m_inverse; const PreconditionerA &a_preconditioner; @@ -742,7 +747,7 @@ BlockSchurPreconditioner::BlockSchurPrecondit const PreconditionerA &Apreconditioner ) : - system_matrix (&S), + stokes_matrix (&S), m_inverse (&Mpinv), a_preconditioner (Apreconditioner), tmp (S.block(1,1).m()) @@ -775,7 +780,7 @@ void BlockSchurPreconditioner::vmult ( const BlockVector &src) const { a_preconditioner.vmult (dst.block(0), src.block(0)); - system_matrix->block(1,0).residual(tmp, dst.block(0), src.block(1)); + stokes_matrix->block(1,0).residual(tmp, dst.block(0), src.block(1)); tmp *= -1; m_inverse->vmult (dst.block(1), tmp); } @@ -801,12 +806,13 @@ void BlockSchurPreconditioner::vmult ( template BoussinesqFlowProblem::BoussinesqFlowProblem (const unsigned int degree) : - degree (degree), triangulation (Triangulation::maximum_smoothing), - fe (FE_Q(degree+1), dim, - FE_Q(degree), 1, - FE_Q(degree), 1), - dof_handler (triangulation), + degree (degree), + stokes_fe (FE_Q(degree+1), dim, + FE_Q(degree), 1), + stokes_dof_handler (triangulation), + temperature_fe (degree), + temperature_dof_handler (triangulation), time_step (0), old_time_step (0), timestep_number (0), @@ -860,32 +866,41 @@ BoussinesqFlowProblem::BoussinesqFlowProblem (const unsigned int degree) template void BoussinesqFlowProblem::setup_dofs () { - dof_handler.distribute_dofs (fe); - DoFRenumbering::Cuthill_McKee (dof_handler); - std::vector block_component (dim+2,0); - block_component[dim] = 1; - block_component[dim+1] = 2; - DoFRenumbering::component_wise (dof_handler, block_component); - - // The next step is, as usual, - // to write some information - // to the screen. The information - // that is most interesting during - // the calculations is the - // number of degrees of freedom - // in the individual components, - // so we count them. The function - // to do this is the same as the - // one used in step-22, which - // uses the grouping of all - // velocity components into - // one block as introduced - // above. - std::vector dofs_per_block (3); - DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component); - const unsigned int n_u = dofs_per_block[0], - n_p = dofs_per_block[1], - n_T = dofs_per_block[2]; + std::vector stokes_block_component (dim+1,0); + stokes_block_component[dim] = 1; + + { + stokes_dof_handler.distribute_dofs (stokes_fe); + DoFRenumbering::Cuthill_McKee (stokes_dof_handler); + DoFRenumbering::component_wise (stokes_dof_handler, stokes_block_component); + + stokes_constraints.clear (); + DoFTools::make_hanging_node_constraints (stokes_dof_handler, + stokes_constraints); + std::set no_normal_flux_boundaries; + no_normal_flux_boundaries.insert (0); + VectorTools::compute_no_normal_flux_constraints (stokes_dof_handler, 0, + no_normal_flux_boundaries, + stokes_constraints); + stokes_constraints.close (); + } + { + temperature_dof_handler.distribute_dofs (temperature_fe); + DoFRenumbering::Cuthill_McKee (temperature_dof_handler); + + temperature_constraints.clear (); + DoFTools::make_hanging_node_constraints (temperature_dof_handler, + temperature_constraints); + temperature_constraints.close (); + } + + std::vector stokes_dofs_per_block (2); + DoFTools::count_dofs_per_block (stokes_dof_handler, stokes_dofs_per_block, + stokes_block_component); + + const unsigned int n_u = stokes_dofs_per_block[0], + n_p = stokes_dofs_per_block[1], + n_T = temperature_dof_handler.n_dofs(); std::cout << "Number of active cells: " << triangulation.n_active_cells() @@ -894,20 +909,11 @@ void BoussinesqFlowProblem::setup_dofs () << " levels)" << std::endl << "Number of degrees of freedom: " - << dof_handler.n_dofs() + << n_u + n_p + n_T << " (" << n_u << '+' << n_p << '+'<< n_T <<')' << std::endl << std::endl; - hanging_node_constraints.clear (); - DoFTools::make_hanging_node_constraints (dof_handler, - hanging_node_constraints); - std::set no_normal_flux_boundaries; - no_normal_flux_boundaries.insert (0); - VectorTools::compute_no_normal_flux_constraints (dof_handler, 0, - no_normal_flux_boundaries, - hanging_node_constraints); - hanging_node_constraints.close (); // The next step is to @@ -960,7 +966,7 @@ void BoussinesqFlowProblem::setup_dofs () // Observe that we use a // coupling argument for // telling the function - // make_sparsity_pattern + // make_stokes_sparsity_pattern // which components actually // will hold data and which // we're going to neglect. @@ -970,23 +976,18 @@ void BoussinesqFlowProblem::setup_dofs () // system matrix structure to // the sparsity pattern. { - system_matrix.clear (); + stokes_matrix.clear (); - BlockCompressedSetSparsityPattern csp (3,3); + BlockCompressedSetSparsityPattern csp (2,2); csp.block(0,0).reinit (n_u, n_u); csp.block(0,1).reinit (n_u, n_p); - csp.block(0,2).reinit (n_u, n_T); csp.block(1,0).reinit (n_p, n_u); csp.block(1,1).reinit (n_p, n_p); - csp.block(1,2).reinit (n_p, n_T); - csp.block(2,0).reinit (n_T, n_u); - csp.block(2,1).reinit (n_T, n_p); - csp.block(2,2).reinit (n_T, n_T); csp.collect_sizes (); - Table<2,DoFTools::Coupling> coupling (dim+2, dim+2); + Table<2,DoFTools::Coupling> coupling (dim+1, dim+1); // build the sparsity pattern. note // that all dim velocities couple with @@ -997,15 +998,13 @@ void BoussinesqFlowProblem::setup_dofs () { case 2: { - static const bool coupling_matrix[4][4] - = {{ 1, 1, 1, 0 }, - { 1, 1, 1, 0 }, + static const bool coupling_matrix[3][3] + = {{ 1, 1, 1}, + { 1, 1, 1}, - { 1, 1, 0, 0 }, - - { 0, 0, 0, 1 }}; - for (unsigned int c=0; c::setup_dofs () case 3: { - static const bool coupling_matrix[5][5] - = {{ 1, 1, 1, 1, 0 }, - { 1, 1, 1, 1, 0 }, - { 1, 1, 1, 1, 0 }, + static const bool coupling_matrix[4][4] + = {{ 1, 1, 1, 1}, + { 1, 1, 1, 1}, + { 1, 1, 1, 1}, - { 1, 1, 1, 0, 0 }, - - { 0, 0, 0, 0, 1 }}; - for (unsigned int c=0; c::setup_dofs () Assert (false, ExcNotImplemented()); } - DoFTools::make_sparsity_pattern (dof_handler, coupling, csp); - hanging_node_constraints.condense (csp); - sparsity_pattern.copy_from (csp); + DoFTools::make_sparsity_pattern (stokes_dof_handler, coupling, csp); + stokes_constraints.condense (csp); + stokes_sparsity_pattern.copy_from (csp); - system_matrix.reinit (sparsity_pattern); + stokes_matrix.reinit (stokes_sparsity_pattern); } - + { Amg_preconditioner.reset (); Mp_preconditioner.reset (); - preconditioner_matrix.clear (); + stokes_preconditioner_matrix.clear (); - BlockCompressedSetSparsityPattern csp (3,3); + BlockCompressedSetSparsityPattern csp (2,2); csp.block(0,0).reinit (n_u, n_u); csp.block(0,1).reinit (n_u, n_p); - csp.block(0,2).reinit (n_u, n_T); csp.block(1,0).reinit (n_p, n_u); csp.block(1,1).reinit (n_p, n_p); - csp.block(1,2).reinit (n_p, n_T); - csp.block(2,0).reinit (n_T, n_u); - csp.block(2,1).reinit (n_T, n_p); - csp.block(2,2).reinit (n_T, n_T); csp.collect_sizes (); - Table<2,DoFTools::Coupling> coupling (dim+2, dim+2); + Table<2,DoFTools::Coupling> coupling (dim+1, dim+1); // build the sparsity pattern. note // that all dim velocities couple with @@ -1075,15 +1067,13 @@ void BoussinesqFlowProblem::setup_dofs () { case 2: { - static const bool coupling_matrix[4][4] - = {{ 1, 0, 0, 0 }, - { 0, 1, 0, 0 }, + static const bool coupling_matrix[3][3] + = {{ 1, 0, 0}, + { 0, 1, 0}, - { 0, 0, 1, 0 }, - - { 0, 0, 0, 0 }}; - for (unsigned int c=0; c::setup_dofs () case 3: { - static const bool coupling_matrix[5][5] - = {{ 1, 0, 0, 0, 0 }, - { 0, 1, 0, 0, 0 }, - { 0, 0, 1, 0, 0 }, + static const bool coupling_matrix[4][4] + = {{ 1, 0, 0, 0}, + { 0, 1, 0, 0}, + { 0, 0, 1, 0}, - { 0, 0, 0, 1, 0 }, - - { 0, 0, 0, 0, 0 }}; - for (unsigned int c=0; c::setup_dofs () default: Assert (false, ExcNotImplemented()); - } - - DoFTools::make_sparsity_pattern (dof_handler, coupling, csp); - hanging_node_constraints.condense (csp); - preconditioner_sparsity_pattern.copy_from (csp); + } + - preconditioner_matrix.reinit (preconditioner_sparsity_pattern); + DoFTools::make_sparsity_pattern (stokes_dof_handler, coupling, csp); + stokes_constraints.condense (csp); + stokes_preconditioner_sparsity_pattern.copy_from (csp); + + stokes_preconditioner_matrix.reinit (stokes_preconditioner_sparsity_pattern); } - // As last action in this function, - // we need to set the vectors - // for the solution, the old - // solution (required for - // time stepping) and the system - // right hand side to the - // three-block structure given - // by velocity, pressure and - // temperature. - solution.reinit (3); - solution.block(0).reinit (n_u); - solution.block(1).reinit (n_p); - solution.block(2).reinit (n_T); - solution.collect_sizes (); - - old_solution.reinit (3); - old_solution.block(0).reinit (n_u); - old_solution.block(1).reinit (n_p); - old_solution.block(2).reinit (n_T); - old_solution.collect_sizes (); - - old_old_solution.reinit (3); - old_old_solution.block(0).reinit (n_u); - old_old_solution.block(1).reinit (n_p); - old_old_solution.block(2).reinit (n_T); - old_old_solution.collect_sizes (); + { + temperature_matrix.clear (); + + CompressedSetSparsityPattern csp (n_T, n_T); + DoFTools::make_sparsity_pattern (temperature_dof_handler, csp); + temperature_constraints.condense (csp); + temperature_sparsity_pattern.copy_from (csp); + + temperature_matrix.reinit (temperature_sparsity_pattern); + } + + - system_rhs.reinit (3); - system_rhs.block(0).reinit (n_u); - system_rhs.block(1).reinit (n_p); - system_rhs.block(2).reinit (n_T); - system_rhs.collect_sizes (); + // As last action in this function, + // we need to set the vectors + // for the solution, the old + // solution (required for + // time stepping) and the system + // right hand side to the + // three-block structure given + // by velocity, pressure and + // temperature. + stokes_solution.reinit (2); + stokes_solution.block(0).reinit (n_u); + stokes_solution.block(1).reinit (n_p); + stokes_solution.collect_sizes (); + + stokes_rhs.reinit (2); + stokes_rhs.block(0).reinit (n_u); + stokes_rhs.block(1).reinit (n_p); + stokes_rhs.collect_sizes (); + + temperature_solution.reinit (n_T); + old_temperature_solution.reinit (n_T); + old_old_temperature_solution.reinit (n_T); + + temperature_rhs.reinit (n_T); } @@ -1173,16 +1167,16 @@ double scalar_product (const Tensor<2,dim> &t1, template void -BoussinesqFlowProblem::assemble_preconditioner () +BoussinesqFlowProblem::assemble_stokes_preconditioner () { - preconditioner_matrix = 0; + stokes_preconditioner_matrix = 0; QGauss quadrature_formula(degree+2); - FEValues fe_values (fe, quadrature_formula, - update_JxW_values | - update_values | - update_gradients); - const unsigned int dofs_per_cell = fe.dofs_per_cell; + FEValues stokes_fe_values (stokes_fe, quadrature_formula, + update_JxW_values | + update_values | + update_gradients); + const unsigned int dofs_per_cell = stokes_fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); @@ -1196,19 +1190,19 @@ BoussinesqFlowProblem::assemble_preconditioner () const FEValuesExtractors::Scalar pressure (dim); typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); + cell = stokes_dof_handler.begin_active(), + endc = stokes_dof_handler.end(); for (; cell!=endc; ++cell) { - fe_values.reinit (cell); + stokes_fe_values.reinit (cell); local_matrix = 0; for (unsigned int q=0; q::assemble_preconditioner () local_matrix(i,j) += (scalar_product (phi_grad_u[i], phi_grad_u[j]) + phi_p[i] * phi_p[j]) - * fe_values.JxW(q); + * stokes_fe_values.JxW(q); } cell->get_dof_indices (local_dof_indices); for (unsigned int i=0; i +void +BoussinesqFlowProblem::build_stokes_preconditioner () +{ + std::cout << " Rebuilding Stokes preconditioner..." << std::flush; + + + // This last step of the assembly + // function sets up the preconditioners + // used for the solution of the + // system. We are going to use an + // ILU preconditioner for the + // velocity block (to be used + // by BlockSchurPreconditioner class) + // as well as an ILU preconditioner + // for the inversion of the + // pressure mass matrix. Recall that + // the velocity-velocity block sits + // at position (0,0) in the + // global system matrix, and + // the pressure mass matrix in + // (1,1). The + // storage of these objects is + // as in step-22, that is, we + // include them using a + // shared pointer structure from the + // boost library. + // + // When all work is done, we + // change the flags + // rebuild_preconditioner + // and + // rebuild_matrices + // to false. + assemble_stokes_preconditioner (); + + Amg_preconditioner = boost::shared_ptr + (new PreconditionerTrilinosAmg()); + + const unsigned int n_u = stokes_preconditioner_matrix.block(0,0).m(); + std::vector null_space (dim * n_u, 0.); + + std::vector precondition_dof_list (stokes_dof_handler.n_dofs(), false); + + for (unsigned int component=0; component < dim; ++component) + { + std::vector precondition_mask (dim + 1, false); + precondition_mask[component] = true; + DoFTools::extract_dofs (stokes_dof_handler, precondition_mask, + precondition_dof_list); + + // TODO: The current implementation + // assumes that we are working on + // the first components of a system when + // writing into the null vector. + // Change this to the general case, + // probably use something similar as + // for block vectors. + unsigned int counter = 0; + for (unsigned int i=0; iinitialize(stokes_preconditioner_matrix.block(0,0), + null_space, dim, true, true, false); + + // TODO: we could throw away the (0,0) + // block here since things have been + // copied over to Trilinos. we need to + // keep the (1,1) block, though + + Mp_preconditioner + = boost::shared_ptr > + (new SparseILU); + Mp_preconditioner->initialize (stokes_preconditioner_matrix.block(1,1), + SparseILU::AdditionalData()); + + std::cout << std::endl; + + rebuild_preconditioner = false; } - // @sect4{BoussinesqFlowProblem::assemble_system} + // @sect4{BoussinesqFlowProblem::assemble_stokes_system} // // The assembly of the Boussinesq // system is acutally a two-step @@ -1302,32 +1391,38 @@ BoussinesqFlowProblem::assemble_preconditioner () // the local dofs compared to // the global system. template -void BoussinesqFlowProblem::assemble_system () +void BoussinesqFlowProblem::assemble_stokes_system () { + std::cout << " Assembling..." << std::flush; + if (rebuild_matrices == true) - system_matrix=0; + stokes_matrix=0; - system_rhs=0; + stokes_rhs=0; QGauss quadrature_formula(degree+2); QGauss face_quadrature_formula(degree+2); - FEValues fe_values (fe, quadrature_formula, - update_values | - update_quadrature_points | - update_JxW_values | - (rebuild_matrices == true - ? - update_gradients - : - UpdateFlags(0))); - FEFaceValues fe_face_values (fe, face_quadrature_formula, - update_values | - update_normal_vectors | - update_quadrature_points | - update_JxW_values); - - const unsigned int dofs_per_cell = fe.dofs_per_cell; + FEValues stokes_fe_values (stokes_fe, quadrature_formula, + update_values | + update_quadrature_points | + update_JxW_values | + (rebuild_matrices == true + ? + update_gradients + : + UpdateFlags(0))); + + FEValues temperature_fe_values (temperature_fe, quadrature_formula, + update_values); + + FEFaceValues stokes_fe_face_values (stokes_fe, face_quadrature_formula, + update_values | + update_normal_vectors | + update_quadrature_points | + update_JxW_values); + + const unsigned int dofs_per_cell = stokes_fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); const unsigned int n_face_q_points = face_quadrature_formula.size(); @@ -1337,43 +1432,42 @@ void BoussinesqFlowProblem::assemble_system () std::vector local_dof_indices (dofs_per_cell); - // These few declarations provide - // the structures for the evaluation - // of inhomogeneous Neumann boundary - // conditions from the function - // declaration made above. - // The vector old_solution_values - // evaluates the solution - // at the old time level, since - // the temperature from the - // old time level enters the - // Stokes system as a source - // term in the momentum equation. - // - // Then, we create a variable - // to hold the Rayleigh number, - // the measure of buoyancy. - // - // The set of vectors we create - // next hold the evaluations of - // the basis functions that will - // be used for creating the - // matrices. This gives faster - // access to that data, which - // increases the performance - // of the assembly. See step-22 - // for details. - // - // The last few declarations - // are used to extract the - // individual blocks (velocity, - // pressure, temperature) from - // the total FE system. + // These few declarations provide + // the structures for the evaluation + // of inhomogeneous Neumann boundary + // conditions from the function + // declaration made above. + // The vector old_solution_values + // evaluates the solution + // at the old time level, since + // the temperature from the + // old time level enters the + // Stokes system as a source + // term in the momentum equation. + // + // Then, we create a variable + // to hold the Rayleigh number, + // the measure of buoyancy. + // + // The set of vectors we create + // next hold the evaluations of + // the basis functions that will + // be used for creating the + // matrices. This gives faster + // access to that data, which + // increases the performance + // of the assembly. See step-22 + // for details. + // + // The last few declarations + // are used to extract the + // individual blocks (velocity, + // pressure, temperature) from + // the total FE system. const PressureBoundaryValues pressure_boundary_values; std::vector boundary_values (n_face_q_points); - std::vector > old_solution_values(n_q_points, - Vector(dim+2)); + std::vector old_temperature_values(n_q_points); const double Rayleigh_number = 10; @@ -1384,178 +1478,167 @@ void BoussinesqFlowProblem::assemble_system () const FEValuesExtractors::Vector velocities (0); const FEValuesExtractors::Scalar pressure (dim); - const FEValuesExtractors::Scalar temperature (dim+1); - - // Now starts the loop over - // all cells in the problem. - // The first commands are all - // very familiar, doing the - // evaluations of the element - // basis functions, resetting - // the local arrays and - // getting the values of the - // old solution at the - // quadrature point. Then we - // are ready to loop over - // the quadrature points - // on the cell. + + // Now start the loop over + // all cells in the problem. + // The first commands are all + // very familiar, doing the + // evaluations of the element + // basis functions, resetting + // the local arrays and + // getting the values of the + // old solution at the + // quadrature point. Then we + // are ready to loop over + // the quadrature points + // on the cell. typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) + cell = stokes_dof_handler.begin_active(), + endc = stokes_dof_handler.end(); + typename DoFHandler::active_cell_iterator + temperature_cell = temperature_dof_handler.begin_active(); + + for (; cell!=endc; ++cell, ++temperature_cell) { - fe_values.reinit (cell); + stokes_fe_values.reinit (cell); + temperature_fe_values.reinit (temperature_cell); + local_matrix = 0; local_rhs = 0; - fe_values.get_function_values (old_solution, old_solution_values); + temperature_fe_values.get_function_values (old_temperature_solution, old_temperature_values); for (unsigned int q=0; qrebuild_matrices - // flag. + const double old_temperature = old_temperature_values[q]; + + // Extract the basis relevant + // terms in the inner products + // once in advance as shown + // in step-22 in order to + // accelerate assembly. + // + // Once this is done, we + // start the loop over the + // rows and columns of the + // local matrix and feed + // the matrix with the relevant + // products. The right hand + // side is filled with the + // forcing term driven by + // temperature in direction + // of gravity (which is + // vertical in our example). + // Note that the right hand + // side term is always generated, + // whereas the matrix + // contributions are only + // updated when it is + // requested by the + // rebuild_matrices + // flag. for (unsigned int k=0; k gravity = ( (dim == 2) ? (Point (0,1)) : + (Point (0,0,1)) ); for (unsigned int i=0; i phi_i_u = fe_values[velocities].value (i, q); - - // define viscosity - const double eta = 1; - - // build Stokes part of - // the matrix. we have - // to rebuild the - // temperature part of - // it in each time step - // due to the fact that - // it depends on the - // time step size - if (rebuild_matrices) - for (unsigned int j=0; j gravity = ( (dim == 2) ? (Point (0,1)) : - (Point (0,0,1)) ); - - local_rhs(i) += (Rayleigh_number * - gravity * phi_u[i] * old_temperature)* - fe_values.JxW(q); - } + local_rhs(i) += (Rayleigh_number * + gravity * phi_u[i] * old_temperature)* + stokes_fe_values.JxW(q); } - // Next follows the assembly - // of the face terms, result - // from Neumann boundary - // conditions. Since these - // terms only enter the right - // hand side vector and not - // the matrix, there is no - // substantial benefit from - // extracting the data - // before using it, so - // we remain in the lines - // of step-20 at this point. + + // Next follows the assembly + // of the face terms, result + // from Neumann boundary + // conditions. Since these + // terms only enter the right + // hand side vector and not + // the matrix, there is no + // substantial benefit from + // extracting the data + // before using it, so + // we remain in the lines + // of step-20 at this point. for (unsigned int face_no=0; - face_no::faces_per_cell; - ++face_no) - if (cell->at_boundary(face_no)) - { - fe_face_values.reinit (cell, face_no); - - pressure_boundary_values - .value_list (fe_face_values.get_quadrature_points(), - boundary_values); - - for (unsigned int q=0; q - phi_i_u = fe_face_values[velocities].value (i, q); - - local_rhs(i) += -(phi_i_u * - fe_face_values.normal_vector(q) * - boundary_values[q] * - fe_face_values.JxW(q)); - } - } - - // The last step in the loop - // over all cells is to - // enter the local contributions - // into the global matrix and - // vector structures to the - // positions specified in - // local_dof_indices. - // Again, we only add the - // matrix data when it is - // requested. + face_no::faces_per_cell; + ++face_no) + if (cell->at_boundary(face_no)) + { + stokes_fe_face_values.reinit (cell, face_no); + + pressure_boundary_values + .value_list (stokes_fe_face_values.get_quadrature_points(), + boundary_values); + + for (unsigned int q=0; q + phi_i_u = stokes_fe_face_values[velocities].value (i, q); + + local_rhs(i) += -(phi_i_u * + stokes_fe_face_values.normal_vector(q) * + boundary_values[q] * + stokes_fe_face_values.JxW(q)); + } + } + + // The last step in the loop + // over all cells is to + // enter the local contributions + // into the global matrix and + // vector structures to the + // positions specified in + // local_dof_indices. + // Again, we only add the + // matrix data when it is + // requested. cell->get_dof_indices (local_dof_indices); if (rebuild_matrices == true) - { - for (unsigned int i=0; i boundary_values; // typename DoFHandler::active_cell_iterator @@ -1585,107 +1668,14 @@ void BoussinesqFlowProblem::assemble_system () // component_mask); // MatrixTools::apply_boundary_values (boundary_values, -// system_matrix, +// stokes_matrix, // solution, // system_rhs); - } - - // This last step of the assembly - // function sets up the preconditioners - // used for the solution of the - // system. We are going to use an - // ILU preconditioner for the - // velocity block (to be used - // by BlockSchurPreconditioner class) - // as well as an ILU preconditioner - // for the inversion of the - // pressure mass matrix. Recall that - // the velocity-velocity block sits - // at position (0,0) in the - // global system matrix, and - // the pressure mass matrix in - // (1,1). The - // storage of these objects is - // as in step-22, that is, we - // include them using a - // shared pointer structure from the - // boost library. - // - // When all work is done, we - // change the flags - // rebuild_preconditioner - // and - // rebuild_matrices - // to false. - if (rebuild_preconditioner == true) - { - Assert (rebuild_matrices == true, - ExcMessage ("There is no point in rebuilding the preconditioner " - "without a rebuilt matrix!")); - - std::cout << " Rebuilding preconditioner..." << std::flush; - - assemble_preconditioner (); - - Amg_preconditioner = boost::shared_ptr - (new PreconditionerTrilinosAmg()); - - const unsigned int n_u = preconditioner_matrix.block(0,0).m(); - std::vector null_space (dim * n_u, 0.); - - std::vector precondition_dof_list (dof_handler.n_dofs(), false); - std::vector precondition_mask (dim + 2, false); - - for (unsigned int component=0; component < dim; component++) - { - precondition_mask[component] = true; - DoFTools::extract_dofs (dof_handler, precondition_mask, - precondition_dof_list); - precondition_mask[component] = false; - - // TODO: The current implementation - // assumes that we are working on - // the first components of a system when - // writing into the null vector. - // Change this to the general case, - // probably use something similar as - // for block vectors. - unsigned int counter = 0; - for (unsigned int i=0; iinitialize(preconditioner_matrix.block(0,0), - null_space, dim, true, true, false); - - // TODO: we could throw away the (0,0) - // block here since things have been - // copied over to Trilinos. we need to - // keep the (1,1) block, though - - Mp_preconditioner - = boost::shared_ptr > - (new SparseILU); - Mp_preconditioner->initialize (preconditioner_matrix.block(1,1), - SparseILU::AdditionalData()); - - std::cout << std::endl; - - rebuild_preconditioner = false; } rebuild_matrices = false; + + std::cout << std::endl; } @@ -1693,21 +1683,21 @@ void BoussinesqFlowProblem::assemble_system () template double compute_viscosity( - const std::vector > present_solution, - const std::vector > old_solution, - const std::vector > old_old_solution, - const std::vector > > old_solution_grads, - const std::vector > > old_old_solution_grads, - const std::vector > > old_solution_hessians, - const std::vector > > old_old_solution_hessians, - const std::vector gamma_values, - const double kappa, - const double global_u_infty, - const double global_T_infty, - const double global_Omega_diameter, - const double cell_diameter, - const double old_time_step - ) + 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_hessians, + const std::vector > &old_old_temperature_hessians, + const std::vector > &present_stokes_values, + const std::vector &gamma_values, + const double kappa, + const double global_u_infty, + const double global_T_infty, + const double global_Omega_diameter, + const double cell_diameter, + const double old_time_step +) { const double beta = 0.1; const double alpha = 1; @@ -1715,43 +1705,39 @@ double compute_viscosity( if (global_u_infty == 0) return 5e-3 * cell_diameter; - const unsigned int n_q_points = old_solution.size(); + const unsigned int n_q_points = old_temperature.size(); - // Stage 1: calculate residual + // Stage 1: calculate residual double max_residual = 0; double max_velocity = 0; for (unsigned int q=0; q < n_q_points; ++q) { - const double dT_dt = (old_solution[q](dim+1) - old_old_solution[q](dim+1)) - / old_time_step; - double u_grad_T = 0.; + Tensor<1,dim> u; for (unsigned int d=0; d -void BoussinesqFlowProblem::assemble_rhs_T () +void BoussinesqFlowProblem::assemble_temperature_system () { const bool use_bdf2_scheme = (timestep_number != 0); - system_matrix.block(2,2) = 0; + temperature_matrix = 0; + temperature_rhs = 0; QGauss quadrature_formula(degree+2); - FEValues fe_values (fe, quadrature_formula, - update_values | update_gradients | - update_hessians | - update_quadrature_points | update_JxW_values); - - const unsigned int dofs_per_cell = fe.dofs_per_cell; + FEValues temperature_fe_values (temperature_fe, quadrature_formula, + update_values | update_gradients | + update_hessians | + update_quadrature_points | update_JxW_values); + FEValues stokes_fe_values (stokes_fe, quadrature_formula, + update_values); + + const unsigned int dofs_per_cell = temperature_fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); Vector local_rhs (dofs_per_cell); @@ -1804,43 +1793,34 @@ void BoussinesqFlowProblem::assemble_rhs_T () std::vector local_dof_indices (dofs_per_cell); - // Here comes the declaration - // of vectors to hold the old - // and present solution values - // and gradients - // for both the cell as well as faces - // to the cell. Next comes the - // declaration of an object - // to hold the temperature - // boundary values and a - // well-known extractor for - // accessing the temperature - // part of the FE system. - std::vector > present_solution_values (n_q_points, - Vector(dim+2)); - std::vector > old_solution_values(n_q_points, - Vector(dim+2)); - std::vector > old_old_solution_values(n_q_points, - Vector(dim+2)); - std::vector > > - old_solution_grads(n_q_points, - std::vector >(dim+2)); - std::vector > > - old_old_solution_grads(n_q_points, - std::vector >(dim+2)); - std::vector > > old_solution_hessians( - n_q_points, - std::vector >(dim+2)); - std::vector > > old_old_solution_hessians( - n_q_points, - std::vector >(dim+2)); + // Here comes the declaration + // of vectors to hold the old + // and present solution values + // and gradients + // for both the cell as well as faces + // to the cell. Next comes the + // declaration of an object + // to hold the temperature + // boundary values and a + // well-known extractor for + // accessing the temperature + // part of the FE system. + std::vector > present_stokes_values (n_q_points, + Vector(dim+1)); + + + std::vector old_temperature_values (n_q_points); + std::vector old_old_temperature_values(n_q_points); + std::vector > old_temperature_grads(n_q_points); + std::vector > old_old_temperature_grads(n_q_points); + std::vector > old_temperature_hessians(n_q_points); + std::vector > old_old_temperature_hessians(n_q_points); + TemperatureBoundaryValues temperature_boundary_values; - RightHandSide right_hand_side; + TemperatureRightHandSide temperature_right_hand_side; std::vector gamma_values (n_q_points); - const FEValuesExtractors::Scalar temperature (dim+1); - std::vector phi_T (dofs_per_cell); std::vector > grad_phi_T (dofs_per_cell); @@ -1848,39 +1828,43 @@ void BoussinesqFlowProblem::assemble_rhs_T () const double global_T_infty = get_maximal_temperature(); const double global_Omega_diameter = GridTools::diameter (triangulation); - // Now, let's start the loop - // over all cells in the - // triangulation. The first - // actions within the loop - // are, 0as usual, the evaluation - // of the FE basis functions - // and the old and present - // solution at the quadrature - // points. + // Now, let's start the loop + // over all cells in the + // triangulation. The first + // actions within the loop + // are, 0as usual, the evaluation + // of the FE basis functions + // and the old and present + // solution at the quadrature + // points. + typename DoFHandler::active_cell_iterator + cell = temperature_dof_handler.begin_active(), + endc = temperature_dof_handler.end(); typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); + stokes_cell = stokes_dof_handler.begin_active(); - for (; cell!=endc; ++cell) + for (; cell!=endc; ++cell, ++stokes_cell) { local_rhs = 0; local_matrix = 0; - fe_values.reinit (cell); + temperature_fe_values.reinit (cell); + stokes_fe_values.reinit (stokes_cell); - fe_values.get_function_values (solution, present_solution_values); - fe_values.get_function_values (old_solution, old_solution_values); - fe_values.get_function_values (old_old_solution, old_old_solution_values); + temperature_fe_values.get_function_values (old_temperature_solution, old_temperature_values); + temperature_fe_values.get_function_values (old_old_temperature_solution, old_old_temperature_values); - fe_values.get_function_gradients (old_solution, old_solution_grads); - fe_values.get_function_gradients (old_old_solution, old_old_solution_grads); + temperature_fe_values.get_function_gradients (old_temperature_solution, old_temperature_grads); + temperature_fe_values.get_function_gradients (old_old_temperature_solution, old_old_temperature_grads); - fe_values.get_function_hessians (old_solution, old_solution_hessians); - fe_values.get_function_hessians (old_old_solution, old_old_solution_hessians); + temperature_fe_values.get_function_hessians (old_temperature_solution, old_temperature_hessians); + temperature_fe_values.get_function_hessians (old_old_temperature_solution, old_old_temperature_hessians); - right_hand_side.value_list (fe_values.get_quadrature_points(), - gamma_values, dim+1); + temperature_right_hand_side.value_list (temperature_fe_values.get_quadrature_points(), + gamma_values); + stokes_fe_values.get_function_values (stokes_solution, present_stokes_values); + // build matrix contributions // define diffusion. take the @@ -1890,32 +1874,37 @@ void BoussinesqFlowProblem::assemble_rhs_T () // impirically) to keep the // scheme stable const double kappa = 1e-6; - const double nu = - compute_viscosity (present_solution_values, old_solution_values, - old_old_solution_values, old_solution_grads, old_old_solution_grads, - old_solution_hessians, old_old_solution_hessians, gamma_values, - kappa, global_u_infty, global_T_infty, - global_Omega_diameter, cell->diameter(), - old_time_step); + const double nu + = compute_viscosity (old_temperature_values, + old_old_temperature_values, + old_temperature_grads, + old_old_temperature_grads, + old_temperature_hessians, + old_old_temperature_hessians, + present_stokes_values, + gamma_values, + kappa, global_u_infty, global_T_infty, + global_Omega_diameter, cell->diameter(), + old_time_step); for (unsigned int q=0; q old_grad_T = old_solution_grads[q][dim+1]; - const Tensor<1,dim> old_old_grad_T = old_old_solution_grads[q][dim+1]; + const Tensor<1,dim> old_grad_T = old_temperature_grads[q]; + const Tensor<1,dim> old_old_grad_T = old_old_temperature_grads[q]; Tensor<1,dim> present_u; for (unsigned int d=0; d::assemble_rhs_T () + time_step * kappa * grad_phi_T[i] * grad_phi_T[j]) - * fe_values.JxW(q); + * temperature_fe_values.JxW(q); for (unsigned int i=0; i::assemble_rhs_T () time_step * gamma_values[q] * phi_T[i]) * - fe_values.JxW(q); + temperature_fe_values.JxW(q); } else { @@ -1964,7 +1953,7 @@ void BoussinesqFlowProblem::assemble_rhs_T () + time_step * kappa * grad_phi_T[i] * grad_phi_T[j]) - * fe_values.JxW(q); + * temperature_fe_values.JxW(q); for (unsigned int i=0; i::assemble_rhs_T () time_step * gamma_values[q] * phi_T[i]) * - fe_values.JxW(q); + temperature_fe_values.JxW(q); } } @@ -1987,12 +1976,12 @@ void BoussinesqFlowProblem::assemble_rhs_T () for (unsigned int i=0; i::assemble_rhs_T () template void BoussinesqFlowProblem::solve () { - solution = old_solution; - - // Use the BlockMatrixArray structure - // for extracting only the upper left - // 2x2 blocks from the matrix that will - // be used for the solution of the - // blocked system. + std::cout << " Solving..." << std::endl; + + // Use the BlockMatrixArray structure + // for extracting only the upper left + // 2x2 blocks from the matrix that will + // be used for the solution of the + // blocked system. { - GrowingVectorMemory > simple_mem; - BlockMatrixArray stokes_submatrix(2, 2, simple_mem); - - stokes_submatrix.enter(system_matrix.block(0,0),0,0); - stokes_submatrix.enter(system_matrix.block(0,1),0,1); - stokes_submatrix.enter(system_matrix.block(1,0),1,0); - - // Define some temporary vectors - // for the solution process. - // TODO: Can we somehow avoid copying - // these vectors back and forth? I.e. - // accessing the block vectors in a - // similar way as the matrix with the - // BlockMatrixArray class? - std::vector block_sizes(2); - block_sizes[0] = solution.block(0).size(); - block_sizes[1] = solution.block(1).size(); - - BlockVector up_rhs(block_sizes); - BlockVector up(block_sizes); - - up_rhs.block(0) = system_rhs.block(0); - up_rhs.block(1) = system_rhs.block(1); - - // Set up inverse matrix for - // pressure mass matrix + // Set up inverse matrix for + // pressure mass matrix InverseMatrix,SparseILU > - mp_inverse (preconditioner_matrix.block(1,1), *Mp_preconditioner); + mp_inverse (stokes_preconditioner_matrix.block(1,1), *Mp_preconditioner); - // Set up block Schur preconditioner - /*BlockSchurPreconditioner::type, - SparseILU > - preconditioner (system_matrix, mp_inverse, *A_preconditioner);*/ + // Set up block Schur preconditioner + /*BlockSchurPreconditioner::type, + SparseILU > + preconditioner (stokes_matrix, mp_inverse, *A_preconditioner);*/ BlockSchurPreconditioner > - preconditioner (system_matrix, mp_inverse, *Amg_preconditioner); + preconditioner (stokes_matrix, mp_inverse, *Amg_preconditioner); - // Set up GMRES solver and - // solve. - SolverControl solver_control (system_matrix.m(), - 1e-6*system_rhs.l2_norm()); + // Set up GMRES solver and + // solve. + SolverControl solver_control (stokes_matrix.m(), + 1e-6*stokes_rhs.l2_norm()); SolverGMRES > gmres(solver_control, - SolverGMRES >::AdditionalData(100)); + SolverGMRES >::AdditionalData(100)); - gmres.solve(stokes_submatrix, up, up_rhs, preconditioner); + gmres.solve(stokes_matrix, stokes_solution, stokes_rhs, preconditioner); std::cout << " " << solver_control.last_step() << " GMRES iterations for Stokes subsystem." << std::endl; - solution.block(0) = up.block(0); - solution.block(1) = up.block(1); - - // Produce a constistent solution - // field (we can't do this on the 'up' - // vector since it does not have the - // temperature component, but - // hanging_node_constraints has - // constraints also for the - // temperature vector) - hanging_node_constraints.distribute (solution); + // Produce a constistent solution + // field (we can't do this on the 'up' + // vector since it does not have the + // temperature component, but + // hanging_node_constraints has + // constraints also for the + // temperature vector) + stokes_constraints.distribute (stokes_solution); } // TODO: determine limit of stability for @@ -2082,29 +2044,32 @@ void BoussinesqFlowProblem::solve () old_time_step = time_step; time_step = GridTools::minimal_cell_diameter(triangulation) / std::max (get_maximal_velocity(), .05) / 4; + + temperature_solution = old_temperature_solution; - assemble_rhs_T (); + + assemble_temperature_system (); { - SolverControl solver_control (system_matrix.block(2,2).m(), - 1e-8*system_rhs.block(2).l2_norm()); + SolverControl solver_control (temperature_matrix.m(), + 1e-8*temperature_rhs.l2_norm()); SolverCG<> cg (solver_control); PreconditionSSOR<> preconditioner; - preconditioner.initialize (system_matrix.block(2,2), 1.2); + preconditioner.initialize (temperature_matrix, 1.2); - cg.solve (system_matrix.block(2,2), solution.block(2), - system_rhs.block(2), preconditioner); + cg.solve (temperature_matrix, temperature_solution, + temperature_rhs, preconditioner); // produce a consistent temperature field - hanging_node_constraints.distribute (solution); + temperature_constraints.distribute (temperature_solution); std::cout << " " << solver_control.last_step() << " CG iterations for temperature." << std::endl; std::cout << " Max temperature: " - << *std::max_element (solution.block(2).begin(), - solution.block(2).end()) + << *std::max_element (temperature_solution.begin(), + temperature_solution.end()) << std::endl; } } @@ -2118,31 +2083,34 @@ void BoussinesqFlowProblem::output_results () const if (timestep_number % 10 != 0) return; - std::vector solution_names (dim, "velocity"); - solution_names.push_back ("p"); - solution_names.push_back ("T"); +//TODO!! +// Assert (false, ExcNotImplemented()); + +// std::vector solution_names (dim, "velocity"); +// solution_names.push_back ("p"); +// solution_names.push_back ("T"); - DataOut data_out; +// DataOut data_out; - data_out.attach_dof_handler (dof_handler); +// data_out.attach_dof_handler (dof_handler); - std::vector - data_component_interpretation - (dim+2, DataComponentInterpretation::component_is_scalar); - for (unsigned int i=0; i +// data_component_interpretation +// (dim+2, DataComponentInterpretation::component_is_scalar); +// for (unsigned int i=0; i::type_dof_data, - data_component_interpretation); - data_out.build_patches (degree); +// data_out.add_data_vector (solution, solution_names, +// DataOut::type_dof_data, +// data_component_interpretation); +// data_out.build_patches (degree); - std::ostringstream filename; - filename << "solution-" << Utilities::int_to_string(timestep_number, 4) << ".vtk"; +// std::ostringstream filename; +// filename << "solution-" << Utilities::int_to_string(timestep_number, 4) << ".vtk"; - std::ofstream output (filename.str().c_str()); - data_out.write_vtk (output); +// std::ofstream output (filename.str().c_str()); +// data_out.write_vtk (output); } @@ -2153,14 +2121,11 @@ void BoussinesqFlowProblem::refine_mesh (const unsigned int max_grid_level) { Vector estimated_error_per_cell (triangulation.n_active_cells()); - std::vector component_mask (dim+2, false); - component_mask[dim+1] = true; - KellyErrorEstimator::estimate (dof_handler, - QGauss(3), + KellyErrorEstimator::estimate (temperature_dof_handler, + QGauss(degree+1), typename FunctionMap::type(), - solution, - estimated_error_per_cell, - component_mask); + temperature_solution, + estimated_error_per_cell); GridRefinement::refine_and_coarsen_fixed_fraction (triangulation, estimated_error_per_cell, @@ -2172,15 +2137,15 @@ void BoussinesqFlowProblem::refine_mesh (const unsigned int max_grid_level) if (cell->has_children() == false) cell->clear_refine_flag (); - SolutionTransfer soltrans(dof_handler); + SolutionTransfer soltrans(temperature_dof_handler); triangulation.prepare_coarsening_and_refinement(); std::vector > x_solution (2); - x_solution[0].reinit (dof_handler.n_dofs()); - x_solution[0] = solution; - x_solution[1].reinit (dof_handler.n_dofs()); - x_solution[1] = old_solution; + x_solution[0].reinit (temperature_dof_handler.n_dofs()); + x_solution[0] = temperature_solution; + x_solution[1].reinit (temperature_dof_handler.n_dofs()); + x_solution[1] = old_temperature_solution; soltrans.prepare_for_coarsening_and_refinement(x_solution); @@ -2188,15 +2153,16 @@ void BoussinesqFlowProblem::refine_mesh (const unsigned int max_grid_level) setup_dofs (); std::vector > tmp (2); - tmp[0].reinit (dof_handler.n_dofs()); - tmp[1].reinit (dof_handler.n_dofs()); + tmp[0].reinit (temperature_dof_handler.n_dofs()); + tmp[1].reinit (temperature_dof_handler.n_dofs()); soltrans.interpolate(x_solution, tmp); + temperature_solution = tmp[0]; + //TODO: this should be [1]!!!!!!! + old_temperature_solution = tmp[0]; + rebuild_matrices = true; rebuild_preconditioner = true; - - solution = tmp[0]; - old_solution = tmp[0]; } @@ -2209,25 +2175,25 @@ double BoussinesqFlowProblem::get_maximal_velocity () const const unsigned int n_q_points = quadrature_formula.size(); - FEValues fe_values (fe, quadrature_formula, + FEValues fe_values (stokes_fe, quadrature_formula, update_values); - std::vector > solution_values(n_q_points, - Vector(dim+2)); + std::vector > stokes_values(n_q_points, + Vector(dim+1)); double max_velocity = 0; typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); + cell = stokes_dof_handler.begin_active(), + endc = stokes_dof_handler.end(); for (; cell!=endc; ++cell) { fe_values.reinit (cell); - fe_values.get_function_values (solution, solution_values); + fe_values.get_function_values (stokes_solution, stokes_values); for (unsigned int q=0; q velocity; for (unsigned int i=0; i double BoussinesqFlowProblem::get_maximal_temperature () const { QGauss quadrature_formula(degree+2); - const unsigned int n_q_points - = quadrature_formula.size(); + const unsigned int n_q_points = quadrature_formula.size(); - FEValues fe_values (fe, quadrature_formula, + FEValues fe_values (temperature_fe, quadrature_formula, update_values); - std::vector > old_solution_values(n_q_points, - Vector(dim+2)); - std::vector > old_old_solution_values(n_q_points, - Vector(dim+2)); + std::vector old_temperature_values(n_q_points); + std::vector old_old_temperature_values(n_q_points); + double max_temperature = 0; typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); + cell = temperature_dof_handler.begin_active(), + endc = temperature_dof_handler.end(); for (; cell!=endc; ++cell) { fe_values.reinit (cell); - fe_values.get_function_values (old_solution, old_solution_values); - fe_values.get_function_values (old_old_solution, old_old_solution_values); + fe_values.get_function_values (old_temperature_solution, old_temperature_values); + fe_values.get_function_values (old_old_temperature_solution, old_old_temperature_values); for (unsigned int q=0; q::run () start_time_iteration: - VectorTools::project (dof_handler, - hanging_node_constraints, + VectorTools::project (temperature_dof_handler, + temperature_constraints, QGauss(degree+2), - InitialValues(), - old_solution); + TemperatureInitialValues(), + old_temperature_solution); timestep_number = 0; double time = 0; @@ -2314,10 +2278,11 @@ void BoussinesqFlowProblem::run () << ", dt=" << time_step << std::endl; - std::cout << " Assembling..." << std::endl; - assemble_system (); + assemble_stokes_system (); + + if (rebuild_preconditioner == true) + build_stokes_preconditioner (); - std::cout << " Solving..." << std::endl; solve (); output_results (); @@ -2338,8 +2303,8 @@ void BoussinesqFlowProblem::run () time += time_step; ++timestep_number; - old_old_solution = old_solution; - old_solution = solution; + old_old_temperature_solution = old_temperature_solution; + old_temperature_solution = temperature_solution; } while (time <= 10); }