From 127822edec427cd5433ba498f977e789a6ab37f9 Mon Sep 17 00:00:00 2001 From: bangerth Date: Fri, 23 Sep 2011 00:04:16 +0000 Subject: [PATCH] Wrap everything into a namespace Step32 as was done for all other tutorial programs. git-svn-id: https://svn.dealii.org/trunk@24379 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-32/step-32.cc | 5104 ++++++++++++++------------- 1 file changed, 2555 insertions(+), 2549 deletions(-) diff --git a/deal.II/examples/step-32/step-32.cc b/deal.II/examples/step-32/step-32.cc index 5981134d3c..8f7c094490 100644 --- a/deal.II/examples/step-32/step-32.cc +++ b/deal.II/examples/step-32/step-32.cc @@ -94,7 +94,9 @@ -using namespace dealii; +namespace Step32 +{ + using namespace dealii; // @sect3{Equation data} @@ -103,150 +105,150 @@ using namespace dealii; // these are exhaustively discussed in the // description of the testcase in the // introduction: -namespace EquationData -{ - const double eta = 1e21; /* Pa s */ - const double kappa = 1e-6; - const double reference_density = 3300; /* kg / m^3 */ - const double reference_temperature = 293; /* K */ - const double expansion_coefficient = 2e-5; /* 1/K */ - const double specific_heat = 1250; /* J / K / kg */ //?? - const double radiogenic_heating = 7.4e-12; /* W / kg */ //?? + namespace EquationData + { + const double eta = 1e21; /* Pa s */ + const double kappa = 1e-6; + const double reference_density = 3300; /* kg / m^3 */ + const double reference_temperature = 293; /* K */ + const double expansion_coefficient = 2e-5; /* 1/K */ + const double specific_heat = 1250; /* J / K / kg */ //?? + const double radiogenic_heating = 7.4e-12; /* W / kg */ //?? - const double R0 = 6371000.-2890000.; /* m */ - const double R1 = 6371000.- 35000.; /* m */ + const double R0 = 6371000.-2890000.; /* m */ + const double R1 = 6371000.- 35000.; /* m */ - const double T0 = 4000+273; /* K */ - const double T1 = 700+273; /* K */ + const double T0 = 4000+273; /* K */ + const double T1 = 700+273; /* K */ - const double year_in_seconds = 60*60*24*365.2425; + const double year_in_seconds = 60*60*24*365.2425; //TODO: document in intro.dox - // scale not by R1-R0, but by a - // typical length scale, say 10km, - // of variation ("plume - // diameter"). this choice also - // roughly equilibrates the sizes - // of the velocity and pressure - // components of the solution - // vectors - const double pressure_scaling = eta / 10000; - - - double density (const double temperature) - { - return (reference_density * - (1 - expansion_coefficient * (temperature - - reference_temperature))); - } + // scale not by R1-R0, but by a + // typical length scale, say 10km, + // of variation ("plume + // diameter"). this choice also + // roughly equilibrates the sizes + // of the velocity and pressure + // components of the solution + // vectors + const double pressure_scaling = eta / 10000; + + + double density (const double temperature) + { + return (reference_density * + (1 - expansion_coefficient * (temperature - + reference_temperature))); + } - template - Tensor<1,dim> gravity_vector (const Point &p) - { + template + Tensor<1,dim> gravity_vector (const Point &p) + { // interpolate the following values with a physically realistic model: // const double g0 = 10.7; /* m / s^2 */ // const double g1 = 9.81; /* m / s^2 */ - const double r = p.norm(); - return -(1.245e-6 * r + 7.714e13/r/r) * p / p.norm(); - } + const double r = p.norm(); + return -(1.245e-6 * r + 7.714e13/r/r) * p / p.norm(); + } - template - double adiabatic_pressure (const Point &p) - { - // the static, adiabatic pressure - // satisfies - // dP/dr = -g rho - - // assuming a constant density, - // we can integrate the pressure - // equation in depth to get that - // the adiabatic pressure equals - // $P(r) = rho_0 \int_r^{R_1} g(r) dr$ - // - // using the model for the - // gravity vector above, this - // yields the following formula: - const double r = p.norm(); - return reference_density * (1./2 * 1.245e-6 * (R1*R1 - r*r) - 7.714e13 * (1./R1 - 1./r)); - } + template + double adiabatic_pressure (const Point &p) + { + // the static, adiabatic pressure + // satisfies + // dP/dr = -g rho + + // assuming a constant density, + // we can integrate the pressure + // equation in depth to get that + // the adiabatic pressure equals + // $P(r) = rho_0 \int_r^{R_1} g(r) dr$ + // + // using the model for the + // gravity vector above, this + // yields the following formula: + const double r = p.norm(); + return reference_density * (1./2 * 1.245e-6 * (R1*R1 - r*r) - 7.714e13 * (1./R1 - 1./r)); + } - template - double adiabatic_temperature (const Point &p) - { - // the static, adiabatic - // temperature satisfies - // $dT/dr = -T alpha/c_P g$ + template + double adiabatic_temperature (const Point &p) + { + // the static, adiabatic + // temperature satisfies + // $dT/dr = -T alpha/c_P g$ - // let's assume constant gravity, - // then we get by integration - const double r = p.norm(); + // let's assume constant gravity, + // then we get by integration + const double r = p.norm(); - return T1 * std::exp(-expansion_coefficient * 9.81 / specific_heat * (r-R1)); - } + return T1 * std::exp(-expansion_coefficient * 9.81 / specific_heat * (r-R1)); + } - template - class TemperatureInitialValues : public Function - { - public: - TemperatureInitialValues () : Function(1) {} + template + class TemperatureInitialValues : public Function + { + public: + TemperatureInitialValues () : Function(1) {} - virtual double value (const Point &p, - const unsigned int component = 0) const; + virtual double value (const Point &p, + const unsigned int component = 0) const; - virtual void vector_value (const Point &p, - Vector &value) const; - }; + virtual void vector_value (const Point &p, + Vector &value) const; + }; - template - double - TemperatureInitialValues::value (const Point &p, - const unsigned int) const - { - const double r = p.norm(); - const double h = R1-R0; - - // s = fraction of the way from - // the inner to the outer - // boundary; 0<=s<=1 - const double s = (r-R0)/h; - - /* now compute an angular variation of the linear temperature field by - stretching the variable s appropriately. note that the following - formula leaves the end points s=0 and s=1 fixed, but stretches the - region in between depending on the angle phi=atan2(x,y). - - For a plot, see - http://www.wolframalpha.com/input/?i=plot+%28%282*sqrt%28x^2%2By^2%29-1%29%2B0.2*%282*sqrt%28x^2%2By^2%29-1%29*%281-%282*sqrt%28x^2%2By^2%29-1%29%29*sin%286*atan2%28x%2Cy%29%29%29%2C+x%3D-1+to+1%2C+y%3D-1+to+1 - */ - const double scale = (dim==3)?std::max(0.0,cos(3.14159*abs(p(2)/R1))):1.0; - const double phi = std::atan2(p(0),p(1)); - const double s_mod = s - + - 0.2 * s * (1-s) * std::sin(6*phi) * scale; - - return T0*(1.0-s_mod) + T1*s_mod; - } + template + double + TemperatureInitialValues::value (const Point &p, + const unsigned int) const + { + const double r = p.norm(); + const double h = R1-R0; + + // s = fraction of the way from + // the inner to the outer + // boundary; 0<=s<=1 + const double s = (r-R0)/h; + + /* now compute an angular variation of the linear temperature field by + stretching the variable s appropriately. note that the following + formula leaves the end points s=0 and s=1 fixed, but stretches the + region in between depending on the angle phi=atan2(x,y). + + For a plot, see + http://www.wolframalpha.com/input/?i=plot+%28%282*sqrt%28x^2%2By^2%29-1%29%2B0.2*%282*sqrt%28x^2%2By^2%29-1%29*%281-%282*sqrt%28x^2%2By^2%29-1%29%29*sin%286*atan2%28x%2Cy%29%29%29%2C+x%3D-1+to+1%2C+y%3D-1+to+1 + */ + const double scale = (dim==3)?std::max(0.0,cos(3.14159*abs(p(2)/R1))):1.0; + const double phi = std::atan2(p(0),p(1)); + const double s_mod = s + + + 0.2 * s * (1-s) * std::sin(6*phi) * scale; + + return T0*(1.0-s_mod) + T1*s_mod; + } - template - void - TemperatureInitialValues::vector_value (const Point &p, - Vector &values) const - { - for (unsigned int c=0; cn_components; ++c) - values(c) = TemperatureInitialValues::value (p, c); + template + void + TemperatureInitialValues::vector_value (const Point &p, + Vector &values) const + { + for (unsigned int c=0; cn_components; ++c) + values(c) = TemperatureInitialValues::value (p, c); + } } -} @@ -275,73 +277,73 @@ namespace EquationData // respective variable in the // BlockSchurPreconditioner class another // name. -namespace LinearSolvers -{ - template - class RightPrecond : public Subscriptor + namespace LinearSolvers { - public: - RightPrecond ( - const TrilinosWrappers::BlockSparseMatrix &S, - const TrilinosWrappers::BlockSparseMatrix &Spre, - const PreconditionerMp &Mppreconditioner, - const PreconditionerA &Apreconditioner, - const bool do_solve_A_in = true) - : - stokes_matrix (&S), - stokes_preconditioner_matrix (&Spre), - mp_preconditioner (Mppreconditioner), - a_preconditioner (Apreconditioner), - do_solve_A (do_solve_A_in) - {} - - void solve_S(TrilinosWrappers::MPI::Vector &dst, - const TrilinosWrappers::MPI::Vector &src) const - { - SolverControl cn(5000, 1e-5); - - TrilinosWrappers::SolverCG solver(cn); - - solver.solve(stokes_preconditioner_matrix->block(1,1), - dst, src, - mp_preconditioner); - - dst*=-1.0; - } - - void solve_A(TrilinosWrappers::MPI::Vector &dst, - const TrilinosWrappers::MPI::Vector &src) const - { - SolverControl cn(5000, src.l2_norm()*1e-2); - TrilinosWrappers::SolverCG solver(cn); - solver.solve(stokes_matrix->block(0,0), dst, src, a_preconditioner); - } - - void vmult (TrilinosWrappers::MPI::BlockVector &dst, - const TrilinosWrappers::MPI::BlockVector &src) const - { - TrilinosWrappers::MPI::Vector utmp(src.block(0)); - - solve_S(dst.block(1), src.block(1)); - - stokes_matrix->block(0,1).vmult(utmp, dst.block(1)); //B^T - utmp*=-1.0; - utmp.add(src.block(0)); - - if (do_solve_A == true) - solve_A(dst.block(0), utmp); - else - a_preconditioner.vmult (dst.block(0), utmp); - } - - private: - const SmartPointer stokes_matrix; - const SmartPointer stokes_preconditioner_matrix; - const PreconditionerMp &mp_preconditioner; - const PreconditionerA &a_preconditioner; - const bool do_solve_A; - }; -} + template + class RightPrecond : public Subscriptor + { + public: + RightPrecond ( + const TrilinosWrappers::BlockSparseMatrix &S, + const TrilinosWrappers::BlockSparseMatrix &Spre, + const PreconditionerMp &Mppreconditioner, + const PreconditionerA &Apreconditioner, + const bool do_solve_A_in = true) + : + stokes_matrix (&S), + stokes_preconditioner_matrix (&Spre), + mp_preconditioner (Mppreconditioner), + a_preconditioner (Apreconditioner), + do_solve_A (do_solve_A_in) + {} + + void solve_S(TrilinosWrappers::MPI::Vector &dst, + const TrilinosWrappers::MPI::Vector &src) const + { + SolverControl cn(5000, 1e-5); + + TrilinosWrappers::SolverCG solver(cn); + + solver.solve(stokes_preconditioner_matrix->block(1,1), + dst, src, + mp_preconditioner); + + dst*=-1.0; + } + + void solve_A(TrilinosWrappers::MPI::Vector &dst, + const TrilinosWrappers::MPI::Vector &src) const + { + SolverControl cn(5000, src.l2_norm()*1e-2); + TrilinosWrappers::SolverCG solver(cn); + solver.solve(stokes_matrix->block(0,0), dst, src, a_preconditioner); + } + + void vmult (TrilinosWrappers::MPI::BlockVector &dst, + const TrilinosWrappers::MPI::BlockVector &src) const + { + TrilinosWrappers::MPI::Vector utmp(src.block(0)); + + solve_S(dst.block(1), src.block(1)); + + stokes_matrix->block(0,1).vmult(utmp, dst.block(1)); //B^T + utmp*=-1.0; + utmp.add(src.block(0)); + + if (do_solve_A == true) + solve_A(dst.block(0), utmp); + else + a_preconditioner.vmult (dst.block(0), utmp); + } + + private: + const SmartPointer stokes_matrix; + const SmartPointer stokes_preconditioner_matrix; + const PreconditionerMp &mp_preconditioner; + const PreconditionerA &a_preconditioner; + const bool do_solve_A; + }; + } @@ -407,398 +409,398 @@ namespace LinearSolvers // itself), and provide some additional // vector fields that are used to improve // performance of assembly. -namespace Assembly -{ - namespace Scratch + namespace Assembly { - template - struct StokesPreconditioner + namespace Scratch { - StokesPreconditioner (const FiniteElement &stokes_fe, - const Quadrature &stokes_quadrature, - const Mapping &mapping, - const UpdateFlags update_flags); - StokesPreconditioner (const StokesPreconditioner &data); - - FEValues stokes_fe_values; - - std::vector > grads_phi_u; - std::vector phi_p; - }; - - template - StokesPreconditioner:: - StokesPreconditioner (const FiniteElement &stokes_fe, - const Quadrature &stokes_quadrature, - const Mapping &mapping, - const UpdateFlags update_flags) - : - stokes_fe_values (mapping, stokes_fe, stokes_quadrature, - update_flags), - grads_phi_u (stokes_fe.dofs_per_cell), - phi_p (stokes_fe.dofs_per_cell) - {} + template + struct StokesPreconditioner + { + StokesPreconditioner (const FiniteElement &stokes_fe, + const Quadrature &stokes_quadrature, + const Mapping &mapping, + const UpdateFlags update_flags); + StokesPreconditioner (const StokesPreconditioner &data); + FEValues stokes_fe_values; + std::vector > grads_phi_u; + std::vector phi_p; + }; - template - StokesPreconditioner:: - StokesPreconditioner (const StokesPreconditioner &scratch) - : - stokes_fe_values (scratch.stokes_fe_values.get_mapping(), - scratch.stokes_fe_values.get_fe(), - scratch.stokes_fe_values.get_quadrature(), - scratch.stokes_fe_values.get_update_flags()), - grads_phi_u (scratch.grads_phi_u), - phi_p (scratch.phi_p) - {} - - - - // Observe that we derive the - // StokesSystem scratch array from the - // StokesPreconditioner array. We do this - // because all the objects that are - // necessary for the assembly of the - // preconditioner are also needed for the - // actual matrix system and right hand - // side, plus some extra data. This makes - // the program more compact. Note also - // that the assembly of the Stokes system - // and the temperature right hand side - // further down requires data from - // temperature and velocity, - // respectively, so we actually need two - // FEValues objects for those two cases. - template - struct StokesSystem : public StokesPreconditioner - { - StokesSystem (const FiniteElement &stokes_fe, - const Mapping &mapping, - const Quadrature &stokes_quadrature, - const UpdateFlags stokes_update_flags, - const FiniteElement &temperature_fe, - const UpdateFlags temperature_update_flags); + template + StokesPreconditioner:: + StokesPreconditioner (const FiniteElement &stokes_fe, + const Quadrature &stokes_quadrature, + const Mapping &mapping, + const UpdateFlags update_flags) + : + stokes_fe_values (mapping, stokes_fe, stokes_quadrature, + update_flags), + grads_phi_u (stokes_fe.dofs_per_cell), + phi_p (stokes_fe.dofs_per_cell) + {} - StokesSystem (const StokesSystem &data); - FEValues temperature_fe_values; - std::vector > phi_u; - std::vector > grads_phi_u; - std::vector div_phi_u; + template + StokesPreconditioner:: + StokesPreconditioner (const StokesPreconditioner &scratch) + : + stokes_fe_values (scratch.stokes_fe_values.get_mapping(), + scratch.stokes_fe_values.get_fe(), + scratch.stokes_fe_values.get_quadrature(), + scratch.stokes_fe_values.get_update_flags()), + grads_phi_u (scratch.grads_phi_u), + phi_p (scratch.phi_p) + {} - std::vector old_temperature_values; - }; - template - StokesSystem:: - StokesSystem (const FiniteElement &stokes_fe, - const Mapping &mapping, - const Quadrature &stokes_quadrature, - const UpdateFlags stokes_update_flags, - const FiniteElement &temperature_fe, - const UpdateFlags temperature_update_flags) - : - StokesPreconditioner (stokes_fe, stokes_quadrature, - mapping, - stokes_update_flags), - temperature_fe_values (mapping, temperature_fe, stokes_quadrature, - temperature_update_flags), - phi_u (stokes_fe.dofs_per_cell), - grads_phi_u (stokes_fe.dofs_per_cell), - div_phi_u (stokes_fe.dofs_per_cell), - old_temperature_values (stokes_quadrature.size()) - {} + // Observe that we derive the + // StokesSystem scratch array from the + // StokesPreconditioner array. We do this + // because all the objects that are + // necessary for the assembly of the + // preconditioner are also needed for the + // actual matrix system and right hand + // side, plus some extra data. This makes + // the program more compact. Note also + // that the assembly of the Stokes system + // and the temperature right hand side + // further down requires data from + // temperature and velocity, + // respectively, so we actually need two + // FEValues objects for those two cases. + template + struct StokesSystem : public StokesPreconditioner + { + StokesSystem (const FiniteElement &stokes_fe, + const Mapping &mapping, + const Quadrature &stokes_quadrature, + const UpdateFlags stokes_update_flags, + const FiniteElement &temperature_fe, + const UpdateFlags temperature_update_flags); + StokesSystem (const StokesSystem &data); - template - StokesSystem:: - StokesSystem (const StokesSystem &scratch) - : - StokesPreconditioner (scratch), - temperature_fe_values (scratch.temperature_fe_values.get_mapping(), - scratch.temperature_fe_values.get_fe(), - scratch.temperature_fe_values.get_quadrature(), - scratch.temperature_fe_values.get_update_flags()), - phi_u (scratch.phi_u), - grads_phi_u (scratch.grads_phi_u), - div_phi_u (scratch.div_phi_u), - old_temperature_values (scratch.old_temperature_values) - {} + FEValues temperature_fe_values; + std::vector > phi_u; + std::vector > grads_phi_u; + std::vector div_phi_u; + std::vector old_temperature_values; + }; - template - struct TemperatureMatrix - { - TemperatureMatrix (const FiniteElement &temperature_fe, - const Mapping &mapping, - const Quadrature &temperature_quadrature); - TemperatureMatrix (const TemperatureMatrix &data); - FEValues temperature_fe_values; + template + StokesSystem:: + StokesSystem (const FiniteElement &stokes_fe, + const Mapping &mapping, + const Quadrature &stokes_quadrature, + const UpdateFlags stokes_update_flags, + const FiniteElement &temperature_fe, + const UpdateFlags temperature_update_flags) + : + StokesPreconditioner (stokes_fe, stokes_quadrature, + mapping, + stokes_update_flags), + temperature_fe_values (mapping, temperature_fe, stokes_quadrature, + temperature_update_flags), + phi_u (stokes_fe.dofs_per_cell), + grads_phi_u (stokes_fe.dofs_per_cell), + div_phi_u (stokes_fe.dofs_per_cell), + old_temperature_values (stokes_quadrature.size()) + {} - std::vector phi_T; - std::vector > grad_phi_T; - }; - template - TemperatureMatrix:: - TemperatureMatrix (const FiniteElement &temperature_fe, - const Mapping &mapping, - const Quadrature &temperature_quadrature) - : - temperature_fe_values (mapping, - temperature_fe, temperature_quadrature, - update_values | update_gradients | - update_JxW_values), - phi_T (temperature_fe.dofs_per_cell), - grad_phi_T (temperature_fe.dofs_per_cell) - {} + template + StokesSystem:: + StokesSystem (const StokesSystem &scratch) + : + StokesPreconditioner (scratch), + temperature_fe_values (scratch.temperature_fe_values.get_mapping(), + scratch.temperature_fe_values.get_fe(), + scratch.temperature_fe_values.get_quadrature(), + scratch.temperature_fe_values.get_update_flags()), + phi_u (scratch.phi_u), + grads_phi_u (scratch.grads_phi_u), + div_phi_u (scratch.div_phi_u), + old_temperature_values (scratch.old_temperature_values) + {} - template - TemperatureMatrix:: - TemperatureMatrix (const TemperatureMatrix &scratch) - : - temperature_fe_values (scratch.temperature_fe_values.get_mapping(), - scratch.temperature_fe_values.get_fe(), - scratch.temperature_fe_values.get_quadrature(), - scratch.temperature_fe_values.get_update_flags()), - phi_T (scratch.phi_T), - grad_phi_T (scratch.grad_phi_T) - {} + template + struct TemperatureMatrix + { + TemperatureMatrix (const FiniteElement &temperature_fe, + const Mapping &mapping, + const Quadrature &temperature_quadrature); + TemperatureMatrix (const TemperatureMatrix &data); - template - struct TemperatureRHS - { - TemperatureRHS (const FiniteElement &temperature_fe, - const FiniteElement &stokes_fe, - const Mapping &mapping, - const Quadrature &quadrature); - TemperatureRHS (const TemperatureRHS &data); + FEValues temperature_fe_values; - FEValues temperature_fe_values; - FEValues stokes_fe_values; + std::vector phi_T; + std::vector > grad_phi_T; + }; - std::vector phi_T; - std::vector > grad_phi_T; + template + TemperatureMatrix:: + TemperatureMatrix (const FiniteElement &temperature_fe, + const Mapping &mapping, + const Quadrature &temperature_quadrature) + : + temperature_fe_values (mapping, + temperature_fe, temperature_quadrature, + update_values | update_gradients | + update_JxW_values), + phi_T (temperature_fe.dofs_per_cell), + grad_phi_T (temperature_fe.dofs_per_cell) + {} - std::vector > old_velocity_values; - std::vector > old_old_velocity_values; - std::vector > old_strain_rates; - std::vector > old_old_strain_rates; + template + TemperatureMatrix:: + TemperatureMatrix (const TemperatureMatrix &scratch) + : + temperature_fe_values (scratch.temperature_fe_values.get_mapping(), + scratch.temperature_fe_values.get_fe(), + scratch.temperature_fe_values.get_quadrature(), + scratch.temperature_fe_values.get_update_flags()), + phi_T (scratch.phi_T), + grad_phi_T (scratch.grad_phi_T) + {} - std::vector old_temperature_values; - std::vector old_old_temperature_values; - std::vector > old_temperature_grads; - std::vector > old_old_temperature_grads; - std::vector old_temperature_laplacians; - std::vector old_old_temperature_laplacians; - }; - template - TemperatureRHS:: - TemperatureRHS (const FiniteElement &temperature_fe, - const FiniteElement &stokes_fe, - const Mapping &mapping, - const Quadrature &quadrature) - : - temperature_fe_values (mapping, - temperature_fe, quadrature, - update_values | - update_gradients | - update_hessians | - update_quadrature_points | - update_JxW_values), - stokes_fe_values (mapping, - stokes_fe, quadrature, - update_values | update_gradients), - phi_T (temperature_fe.dofs_per_cell), - grad_phi_T (temperature_fe.dofs_per_cell), - - old_velocity_values (quadrature.size()), - old_old_velocity_values (quadrature.size()), - old_strain_rates (quadrature.size()), - old_old_strain_rates (quadrature.size()), - - old_temperature_values (quadrature.size()), - old_old_temperature_values(quadrature.size()), - old_temperature_grads(quadrature.size()), - old_old_temperature_grads(quadrature.size()), - old_temperature_laplacians(quadrature.size()), - old_old_temperature_laplacians(quadrature.size()) - {} + template + struct TemperatureRHS + { + TemperatureRHS (const FiniteElement &temperature_fe, + const FiniteElement &stokes_fe, + const Mapping &mapping, + const Quadrature &quadrature); + TemperatureRHS (const TemperatureRHS &data); + + FEValues temperature_fe_values; + FEValues stokes_fe_values; + + std::vector phi_T; + std::vector > grad_phi_T; + + std::vector > old_velocity_values; + std::vector > old_old_velocity_values; + + std::vector > old_strain_rates; + std::vector > old_old_strain_rates; + + std::vector old_temperature_values; + std::vector old_old_temperature_values; + std::vector > old_temperature_grads; + std::vector > old_old_temperature_grads; + std::vector old_temperature_laplacians; + std::vector old_old_temperature_laplacians; + }; + + template + TemperatureRHS:: + TemperatureRHS (const FiniteElement &temperature_fe, + const FiniteElement &stokes_fe, + const Mapping &mapping, + const Quadrature &quadrature) + : + temperature_fe_values (mapping, + temperature_fe, quadrature, + update_values | + update_gradients | + update_hessians | + update_quadrature_points | + update_JxW_values), + stokes_fe_values (mapping, + stokes_fe, quadrature, + update_values | update_gradients), + phi_T (temperature_fe.dofs_per_cell), + grad_phi_T (temperature_fe.dofs_per_cell), + + old_velocity_values (quadrature.size()), + old_old_velocity_values (quadrature.size()), + old_strain_rates (quadrature.size()), + old_old_strain_rates (quadrature.size()), + + old_temperature_values (quadrature.size()), + old_old_temperature_values(quadrature.size()), + old_temperature_grads(quadrature.size()), + old_old_temperature_grads(quadrature.size()), + old_temperature_laplacians(quadrature.size()), + old_old_temperature_laplacians(quadrature.size()) + {} - template - TemperatureRHS:: - TemperatureRHS (const TemperatureRHS &scratch) - : - temperature_fe_values (scratch.temperature_fe_values.get_mapping(), - scratch.temperature_fe_values.get_fe(), - scratch.temperature_fe_values.get_quadrature(), - scratch.temperature_fe_values.get_update_flags()), - stokes_fe_values (scratch.stokes_fe_values.get_mapping(), - scratch.stokes_fe_values.get_fe(), - scratch.stokes_fe_values.get_quadrature(), - scratch.stokes_fe_values.get_update_flags()), - phi_T (scratch.phi_T), - grad_phi_T (scratch.grad_phi_T), - - old_velocity_values (scratch.old_velocity_values), - old_old_velocity_values (scratch.old_old_velocity_values), - old_strain_rates (scratch.old_strain_rates), - old_old_strain_rates (scratch.old_old_strain_rates), - - old_temperature_values (scratch.old_temperature_values), - old_old_temperature_values (scratch.old_old_temperature_values), - old_temperature_grads (scratch.old_temperature_grads), - old_old_temperature_grads (scratch.old_old_temperature_grads), - old_temperature_laplacians (scratch.old_temperature_laplacians), - old_old_temperature_laplacians (scratch.old_old_temperature_laplacians) - {} - } + template + TemperatureRHS:: + TemperatureRHS (const TemperatureRHS &scratch) + : + temperature_fe_values (scratch.temperature_fe_values.get_mapping(), + scratch.temperature_fe_values.get_fe(), + scratch.temperature_fe_values.get_quadrature(), + scratch.temperature_fe_values.get_update_flags()), + stokes_fe_values (scratch.stokes_fe_values.get_mapping(), + scratch.stokes_fe_values.get_fe(), + scratch.stokes_fe_values.get_quadrature(), + scratch.stokes_fe_values.get_update_flags()), + phi_T (scratch.phi_T), + grad_phi_T (scratch.grad_phi_T), + + old_velocity_values (scratch.old_velocity_values), + old_old_velocity_values (scratch.old_old_velocity_values), + old_strain_rates (scratch.old_strain_rates), + old_old_strain_rates (scratch.old_old_strain_rates), + + old_temperature_values (scratch.old_temperature_values), + old_old_temperature_values (scratch.old_old_temperature_values), + old_temperature_grads (scratch.old_temperature_grads), + old_old_temperature_grads (scratch.old_old_temperature_grads), + old_temperature_laplacians (scratch.old_temperature_laplacians), + old_old_temperature_laplacians (scratch.old_old_temperature_laplacians) + {} + } - // The CopyData arrays are similar to the - // Scratch arrays. They provide a - // constructor, a copy operation, and - // some arrays for local matrix, local - // vectors and the relation between local - // and global degrees of freedom (a.k.a. - // local_dof_indices). - namespace CopyData - { - template - struct StokesPreconditioner + // The CopyData arrays are similar to the + // Scratch arrays. They provide a + // constructor, a copy operation, and + // some arrays for local matrix, local + // vectors and the relation between local + // and global degrees of freedom (a.k.a. + // local_dof_indices). + namespace CopyData { - StokesPreconditioner (const FiniteElement &stokes_fe); - StokesPreconditioner (const StokesPreconditioner &data); - - FullMatrix local_matrix; - std::vector local_dof_indices; - }; - - template - StokesPreconditioner:: - StokesPreconditioner (const FiniteElement &stokes_fe) - : - local_matrix (stokes_fe.dofs_per_cell, - stokes_fe.dofs_per_cell), - local_dof_indices (stokes_fe.dofs_per_cell) - {} - - - - template - StokesPreconditioner:: - StokesPreconditioner (const StokesPreconditioner &data) - : - local_matrix (data.local_matrix), - local_dof_indices (data.local_dof_indices) - {} - + template + struct StokesPreconditioner + { + StokesPreconditioner (const FiniteElement &stokes_fe); + StokesPreconditioner (const StokesPreconditioner &data); + + FullMatrix local_matrix; + std::vector local_dof_indices; + }; + + template + StokesPreconditioner:: + StokesPreconditioner (const FiniteElement &stokes_fe) + : + local_matrix (stokes_fe.dofs_per_cell, + stokes_fe.dofs_per_cell), + local_dof_indices (stokes_fe.dofs_per_cell) + {} - template - struct StokesSystem : public StokesPreconditioner - { - StokesSystem (const FiniteElement &stokes_fe); - StokesSystem (const StokesSystem &data); - Vector local_rhs; - }; + template + StokesPreconditioner:: + StokesPreconditioner (const StokesPreconditioner &data) + : + local_matrix (data.local_matrix), + local_dof_indices (data.local_dof_indices) + {} - template - StokesSystem:: - StokesSystem (const FiniteElement &stokes_fe) - : - StokesPreconditioner (stokes_fe), - local_rhs (stokes_fe.dofs_per_cell) - {} + template + struct StokesSystem : public StokesPreconditioner + { + StokesSystem (const FiniteElement &stokes_fe); + StokesSystem (const StokesSystem &data); - template - StokesSystem:: - StokesSystem (const StokesSystem &data) - : - StokesPreconditioner (data), - local_rhs (data.local_rhs) - {} + Vector local_rhs; + }; + template + StokesSystem:: + StokesSystem (const FiniteElement &stokes_fe) + : + StokesPreconditioner (stokes_fe), + local_rhs (stokes_fe.dofs_per_cell) + {} - template - struct TemperatureMatrix - { - TemperatureMatrix (const FiniteElement &temperature_fe); - TemperatureMatrix (const TemperatureMatrix &data); - FullMatrix local_mass_matrix; - FullMatrix local_stiffness_matrix; - std::vector local_dof_indices; - }; + template + StokesSystem:: + StokesSystem (const StokesSystem &data) + : + StokesPreconditioner (data), + local_rhs (data.local_rhs) + {} - template - TemperatureMatrix:: - TemperatureMatrix (const FiniteElement &temperature_fe) - : - local_mass_matrix (temperature_fe.dofs_per_cell, - temperature_fe.dofs_per_cell), - local_stiffness_matrix (temperature_fe.dofs_per_cell, - temperature_fe.dofs_per_cell), - local_dof_indices (temperature_fe.dofs_per_cell) - {} - template - TemperatureMatrix:: - TemperatureMatrix (const TemperatureMatrix &data) - : - local_mass_matrix (data.local_mass_matrix), - local_stiffness_matrix (data.local_stiffness_matrix), - local_dof_indices (data.local_dof_indices) - {} + template + struct TemperatureMatrix + { + TemperatureMatrix (const FiniteElement &temperature_fe); + TemperatureMatrix (const TemperatureMatrix &data); + + FullMatrix local_mass_matrix; + FullMatrix local_stiffness_matrix; + std::vector local_dof_indices; + }; + + template + TemperatureMatrix:: + TemperatureMatrix (const FiniteElement &temperature_fe) + : + local_mass_matrix (temperature_fe.dofs_per_cell, + temperature_fe.dofs_per_cell), + local_stiffness_matrix (temperature_fe.dofs_per_cell, + temperature_fe.dofs_per_cell), + local_dof_indices (temperature_fe.dofs_per_cell) + {} - template - struct TemperatureRHS - { - TemperatureRHS (const FiniteElement &temperature_fe); - TemperatureRHS (const TemperatureRHS &data); + template + TemperatureMatrix:: + TemperatureMatrix (const TemperatureMatrix &data) + : + local_mass_matrix (data.local_mass_matrix), + local_stiffness_matrix (data.local_stiffness_matrix), + local_dof_indices (data.local_dof_indices) + {} - Vector local_rhs; - std::vector local_dof_indices; - FullMatrix matrix_for_bc; - }; - template - TemperatureRHS:: - TemperatureRHS (const FiniteElement &temperature_fe) - : - local_rhs (temperature_fe.dofs_per_cell), - local_dof_indices (temperature_fe.dofs_per_cell), - matrix_for_bc (temperature_fe.dofs_per_cell, - temperature_fe.dofs_per_cell) - {} + template + struct TemperatureRHS + { + TemperatureRHS (const FiniteElement &temperature_fe); + TemperatureRHS (const TemperatureRHS &data); + + Vector local_rhs; + std::vector local_dof_indices; + FullMatrix matrix_for_bc; + }; + + template + TemperatureRHS:: + TemperatureRHS (const FiniteElement &temperature_fe) + : + local_rhs (temperature_fe.dofs_per_cell), + local_dof_indices (temperature_fe.dofs_per_cell), + matrix_for_bc (temperature_fe.dofs_per_cell, + temperature_fe.dofs_per_cell) + {} - template - TemperatureRHS:: - TemperatureRHS (const TemperatureRHS &data) - : - local_rhs (data.local_rhs), - local_dof_indices (data.local_dof_indices), - matrix_for_bc (data.matrix_for_bc) - {} + template + TemperatureRHS:: + TemperatureRHS (const TemperatureRHS &data) + : + local_rhs (data.local_rhs), + local_dof_indices (data.local_dof_indices), + matrix_for_bc (data.matrix_for_bc) + {} + } } -} @@ -886,321 +888,321 @@ namespace Assembly // program (this choice is made in the // from this variable in the results // section of this tutorial program. -template -class BoussinesqFlowProblem -{ - public: - struct Parameters; - BoussinesqFlowProblem (Parameters ¶meters); - void run (); - - private: - void setup_dofs (); - void assemble_stokes_preconditioner (); - void build_stokes_preconditioner (); - void assemble_stokes_system (); - void assemble_temperature_matrix (); - void assemble_temperature_system (const double maximal_velocity); - void project_temperature_field (); - double get_maximal_velocity () const; - double get_cfl_number () const; - double get_entropy_variation (const double average_temperature) const; - std::pair get_extrapolated_temperature_range () const; - void solve (); - void output_results (); - void refine_mesh (const unsigned int max_grid_level); + template + class BoussinesqFlowProblem + { + public: + struct Parameters; + BoussinesqFlowProblem (Parameters ¶meters); + void run (); - double - 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 > &old_strain_rates, - const std::vector > &old_old_strain_rates, - const double global_u_infty, - const double global_T_variation, - const double average_temperature, - const double global_entropy_variation, - const double cell_diameter) const; - - public: - struct Parameters - { - Parameters (const std::string ¶meter_filename); + private: + void setup_dofs (); + void assemble_stokes_preconditioner (); + void build_stokes_preconditioner (); + void assemble_stokes_system (); + void assemble_temperature_matrix (); + void assemble_temperature_system (const double maximal_velocity); + void project_temperature_field (); + double get_maximal_velocity () const; + double get_cfl_number () const; + double get_entropy_variation (const double average_temperature) const; + std::pair get_extrapolated_temperature_range () const; + void solve (); + void output_results (); + void refine_mesh (const unsigned int max_grid_level); + + double + 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 > &old_strain_rates, + const std::vector > &old_old_strain_rates, + const double global_u_infty, + const double global_T_variation, + const double average_temperature, + const double global_entropy_variation, + const double cell_diameter) const; - static void declare_parameters (ParameterHandler &prm); - void parse_parameters (ParameterHandler &prm); + public: + struct Parameters + { + Parameters (const std::string ¶meter_filename); - double end_time; + static void declare_parameters (ParameterHandler &prm); + void parse_parameters (ParameterHandler &prm); - unsigned int initial_global_refinement; - unsigned int initial_adaptive_refinement; + double end_time; - bool generate_graphical_output; - unsigned int graphical_output_interval; + unsigned int initial_global_refinement; + unsigned int initial_adaptive_refinement; - unsigned int adaptive_refinement_interval; + bool generate_graphical_output; + unsigned int graphical_output_interval; - double stabilization_alpha; - double stabilization_c_R; - double stabilization_beta; + unsigned int adaptive_refinement_interval; - unsigned int stokes_velocity_degree; - bool use_locally_conservative_discretization; + double stabilization_alpha; + double stabilization_c_R; + double stabilization_beta; - unsigned int temperature_degree; - }; + unsigned int stokes_velocity_degree; + bool use_locally_conservative_discretization; - private: - Parameters ¶meters; - ConditionalOStream pcout; + unsigned int temperature_degree; + }; - parallel::distributed::Triangulation triangulation; - double global_Omega_diameter; + private: + Parameters ¶meters; + ConditionalOStream pcout; - const MappingQ mapping; + parallel::distributed::Triangulation triangulation; + double global_Omega_diameter; - const FESystem stokes_fe; + const MappingQ mapping; - DoFHandler stokes_dof_handler; - ConstraintMatrix stokes_constraints; + const FESystem stokes_fe; - TrilinosWrappers::BlockSparseMatrix stokes_matrix; - TrilinosWrappers::BlockSparseMatrix stokes_preconditioner_matrix; + DoFHandler stokes_dof_handler; + ConstraintMatrix stokes_constraints; - TrilinosWrappers::MPI::BlockVector stokes_solution; - TrilinosWrappers::MPI::BlockVector old_stokes_solution; - TrilinosWrappers::MPI::BlockVector stokes_rhs; + TrilinosWrappers::BlockSparseMatrix stokes_matrix; + TrilinosWrappers::BlockSparseMatrix stokes_preconditioner_matrix; + TrilinosWrappers::MPI::BlockVector stokes_solution; + TrilinosWrappers::MPI::BlockVector old_stokes_solution; + TrilinosWrappers::MPI::BlockVector stokes_rhs; - FE_Q temperature_fe; - DoFHandler temperature_dof_handler; - ConstraintMatrix temperature_constraints; - TrilinosWrappers::SparseMatrix temperature_mass_matrix; - TrilinosWrappers::SparseMatrix temperature_stiffness_matrix; - TrilinosWrappers::SparseMatrix temperature_matrix; + FE_Q temperature_fe; + DoFHandler temperature_dof_handler; + ConstraintMatrix temperature_constraints; - TrilinosWrappers::MPI::Vector temperature_solution; - TrilinosWrappers::MPI::Vector old_temperature_solution; - TrilinosWrappers::MPI::Vector old_old_temperature_solution; - TrilinosWrappers::MPI::Vector temperature_rhs; + TrilinosWrappers::SparseMatrix temperature_mass_matrix; + TrilinosWrappers::SparseMatrix temperature_stiffness_matrix; + TrilinosWrappers::SparseMatrix temperature_matrix; + TrilinosWrappers::MPI::Vector temperature_solution; + TrilinosWrappers::MPI::Vector old_temperature_solution; + TrilinosWrappers::MPI::Vector old_old_temperature_solution; + TrilinosWrappers::MPI::Vector temperature_rhs; - double time_step; - double old_time_step; - unsigned int timestep_number; - std_cxx1x::shared_ptr Amg_preconditioner; - std_cxx1x::shared_ptr Mp_preconditioner; - std_cxx1x::shared_ptr T_preconditioner; + double time_step; + double old_time_step; + unsigned int timestep_number; - bool rebuild_stokes_matrix; - bool rebuild_stokes_preconditioner; - bool rebuild_temperature_matrices; - bool rebuild_temperature_preconditioner; + std_cxx1x::shared_ptr Amg_preconditioner; + std_cxx1x::shared_ptr Mp_preconditioner; + std_cxx1x::shared_ptr T_preconditioner; - TimerOutput computing_timer; + bool rebuild_stokes_matrix; + bool rebuild_stokes_preconditioner; + bool rebuild_temperature_matrices; + bool rebuild_temperature_preconditioner; - void setup_stokes_matrix (const std::vector &stokes_partitioning); - void setup_stokes_preconditioner (const std::vector &stokes_partitioning); - void setup_temperature_matrices (const IndexSet &temperature_partitioning); + TimerOutput computing_timer; - void - local_assemble_stokes_preconditioner (const typename DoFHandler::active_cell_iterator &cell, - Assembly::Scratch::StokesPreconditioner &scratch, - Assembly::CopyData::StokesPreconditioner &data); + void setup_stokes_matrix (const std::vector &stokes_partitioning); + void setup_stokes_preconditioner (const std::vector &stokes_partitioning); + void setup_temperature_matrices (const IndexSet &temperature_partitioning); - void - copy_local_to_global_stokes_preconditioner (const Assembly::CopyData::StokesPreconditioner &data); + void + local_assemble_stokes_preconditioner (const typename DoFHandler::active_cell_iterator &cell, + Assembly::Scratch::StokesPreconditioner &scratch, + Assembly::CopyData::StokesPreconditioner &data); + void + copy_local_to_global_stokes_preconditioner (const Assembly::CopyData::StokesPreconditioner &data); - void - local_assemble_stokes_system (const typename DoFHandler::active_cell_iterator &cell, - Assembly::Scratch::StokesSystem &scratch, - Assembly::CopyData::StokesSystem &data); - void - copy_local_to_global_stokes_system (const Assembly::CopyData::StokesSystem &data); + void + local_assemble_stokes_system (const typename DoFHandler::active_cell_iterator &cell, + Assembly::Scratch::StokesSystem &scratch, + Assembly::CopyData::StokesSystem &data); + void + copy_local_to_global_stokes_system (const Assembly::CopyData::StokesSystem &data); - void - local_assemble_temperature_matrix (const typename DoFHandler::active_cell_iterator &cell, - Assembly::Scratch::TemperatureMatrix &scratch, - Assembly::CopyData::TemperatureMatrix &data); - void - copy_local_to_global_temperature_matrix (const Assembly::CopyData::TemperatureMatrix &data); + void + local_assemble_temperature_matrix (const typename DoFHandler::active_cell_iterator &cell, + Assembly::Scratch::TemperatureMatrix &scratch, + Assembly::CopyData::TemperatureMatrix &data); + void + copy_local_to_global_temperature_matrix (const Assembly::CopyData::TemperatureMatrix &data); - void - local_assemble_temperature_rhs (const std::pair global_T_range, - const double global_max_velocity, - const double global_entropy_variation, - const typename DoFHandler::active_cell_iterator &cell, - Assembly::Scratch::TemperatureRHS &scratch, - Assembly::CopyData::TemperatureRHS &data); - void - copy_local_to_global_temperature_rhs (const Assembly::CopyData::TemperatureRHS &data); + void + local_assemble_temperature_rhs (const std::pair global_T_range, + const double global_max_velocity, + const double global_entropy_variation, + const typename DoFHandler::active_cell_iterator &cell, + Assembly::Scratch::TemperatureRHS &scratch, + Assembly::CopyData::TemperatureRHS &data); - class Postprocessor; -}; + void + copy_local_to_global_temperature_rhs (const Assembly::CopyData::TemperatureRHS &data); + + class Postprocessor; + }; // @sect3{BoussinesqFlowProblem class implementation} // @sect4{BoussinesqFlowProblem::Parameters} -template -BoussinesqFlowProblem::Parameters::Parameters (const std::string ¶meter_filename) - : - end_time (1e8), - initial_global_refinement (2), - initial_adaptive_refinement (2), - adaptive_refinement_interval (10), - stabilization_alpha (2), - stabilization_c_R (0.11), - stabilization_beta (0.078), - stokes_velocity_degree (2), - use_locally_conservative_discretization (true), - temperature_degree (2) -{ - ParameterHandler prm; - BoussinesqFlowProblem::Parameters::declare_parameters (prm); + template + BoussinesqFlowProblem::Parameters::Parameters (const std::string ¶meter_filename) + : + end_time (1e8), + initial_global_refinement (2), + initial_adaptive_refinement (2), + adaptive_refinement_interval (10), + stabilization_alpha (2), + stabilization_c_R (0.11), + stabilization_beta (0.078), + stokes_velocity_degree (2), + use_locally_conservative_discretization (true), + temperature_degree (2) + { + ParameterHandler prm; + BoussinesqFlowProblem::Parameters::declare_parameters (prm); - std::ifstream parameter_file (parameter_filename.c_str()); + std::ifstream parameter_file (parameter_filename.c_str()); - if (!parameter_file) - { - parameter_file.close (); + if (!parameter_file) + { + parameter_file.close (); - std::ostringstream message; - message << "Input parameter file <" - << parameter_filename << "> not found. Creating a" - << std::endl - << "template file of the same name." - << std::endl; + std::ostringstream message; + message << "Input parameter file <" + << parameter_filename << "> not found. Creating a" + << std::endl + << "template file of the same name." + << std::endl; - std::ofstream parameter_out (parameter_filename.c_str()); - prm.print_parameters (parameter_out, - ParameterHandler::Text); + std::ofstream parameter_out (parameter_filename.c_str()); + prm.print_parameters (parameter_out, + ParameterHandler::Text); - AssertThrow (false, ExcMessage (message.str().c_str())); - } + AssertThrow (false, ExcMessage (message.str().c_str())); + } - const bool success = prm.read_input (parameter_file); - AssertThrow (success, ExcMessage ("Invalid input parameter file.")); + const bool success = prm.read_input (parameter_file); + AssertThrow (success, ExcMessage ("Invalid input parameter file.")); - parse_parameters (prm); -} + parse_parameters (prm); + } -template -void -BoussinesqFlowProblem::Parameters:: -declare_parameters (ParameterHandler &prm) -{ - prm.declare_entry ("End time", "1e8", - Patterns::Double (0), - "The end time of the simulation in years."); - prm.declare_entry ("Initial global refinement", "2", - Patterns::Integer (0), - "The number of global refinement steps performed on " - "the initial coarse mesh, before the problem is first " - "solved there."); - prm.declare_entry ("Initial adaptive refinement", "2", - Patterns::Integer (0), - "The number of adaptive refinement steps performed after " - "initial global refinement."); - prm.declare_entry ("Time steps between mesh refinement", "10", - Patterns::Integer (1), - "The number of time steps after which the mesh is to be " - "adapted based on computed error indicators."); - prm.declare_entry ("Generate graphical output", "false", - Patterns::Bool (), - "Whether graphical output is to be generated or not. " - "You may not want to get graphical output if the number " - "of processors is large."); - prm.declare_entry ("Time steps between graphical output", "50", - Patterns::Integer (1), - "The number of time steps between each generation of " - "graphical output files."); - - prm.enter_subsection ("Stabilization parameters"); + template + void + BoussinesqFlowProblem::Parameters:: + declare_parameters (ParameterHandler &prm) { - prm.declare_entry ("alpha", "2", - Patterns::Double (1, 2), - "The exponent in the entropy viscosity stabilization."); - prm.declare_entry ("c_R", "0.11", - Patterns::Double (0), - "The c_R factor in the entropy viscosity " - "stabilization."); - prm.declare_entry ("beta", "0.078", - Patterns::Double (0), - "The beta factor in the artificial viscosity " - "stabilization. An appropriate value for 2d is 0.052 " - "and 0.078 for 3d."); - } - prm.leave_subsection (); + prm.declare_entry ("End time", "1e8", + Patterns::Double (0), + "The end time of the simulation in years."); + prm.declare_entry ("Initial global refinement", "2", + Patterns::Integer (0), + "The number of global refinement steps performed on " + "the initial coarse mesh, before the problem is first " + "solved there."); + prm.declare_entry ("Initial adaptive refinement", "2", + Patterns::Integer (0), + "The number of adaptive refinement steps performed after " + "initial global refinement."); + prm.declare_entry ("Time steps between mesh refinement", "10", + Patterns::Integer (1), + "The number of time steps after which the mesh is to be " + "adapted based on computed error indicators."); + prm.declare_entry ("Generate graphical output", "false", + Patterns::Bool (), + "Whether graphical output is to be generated or not. " + "You may not want to get graphical output if the number " + "of processors is large."); + prm.declare_entry ("Time steps between graphical output", "50", + Patterns::Integer (1), + "The number of time steps between each generation of " + "graphical output files."); + + prm.enter_subsection ("Stabilization parameters"); + { + prm.declare_entry ("alpha", "2", + Patterns::Double (1, 2), + "The exponent in the entropy viscosity stabilization."); + prm.declare_entry ("c_R", "0.11", + Patterns::Double (0), + "The c_R factor in the entropy viscosity " + "stabilization."); + prm.declare_entry ("beta", "0.078", + Patterns::Double (0), + "The beta factor in the artificial viscosity " + "stabilization. An appropriate value for 2d is 0.052 " + "and 0.078 for 3d."); + } + prm.leave_subsection (); - prm.enter_subsection ("Discretization"); - { - prm.declare_entry ("Stokes velocity polynomial degree", "2", - Patterns::Integer (1), - "The polynomial degree to use for the velocity variables " - "in the Stokes system."); - prm.declare_entry ("Temperature polynomial degree", "2", - Patterns::Integer (1), - "The polynomial degree to use for the temperature variable."); - prm.declare_entry ("Use locally conservative discretization", "true", - Patterns::Bool (), - "Whether to use a Stokes discretization that is locally " - "conservative at the expense of a larger number of degrees " - "of freedom, or to go with a cheaper discretization " - "that does not locally conserve mass (although it is " - "globally conservative."); + prm.enter_subsection ("Discretization"); + { + prm.declare_entry ("Stokes velocity polynomial degree", "2", + Patterns::Integer (1), + "The polynomial degree to use for the velocity variables " + "in the Stokes system."); + prm.declare_entry ("Temperature polynomial degree", "2", + Patterns::Integer (1), + "The polynomial degree to use for the temperature variable."); + prm.declare_entry ("Use locally conservative discretization", "true", + Patterns::Bool (), + "Whether to use a Stokes discretization that is locally " + "conservative at the expense of a larger number of degrees " + "of freedom, or to go with a cheaper discretization " + "that does not locally conserve mass (although it is " + "globally conservative."); + } + prm.leave_subsection (); } - prm.leave_subsection (); -} -template -void -BoussinesqFlowProblem::Parameters:: -parse_parameters (ParameterHandler &prm) -{ - end_time = prm.get_double ("End time"); - initial_global_refinement = prm.get_integer ("Initial global refinement"); - initial_adaptive_refinement = prm.get_integer ("Initial adaptive refinement"); + template + void + BoussinesqFlowProblem::Parameters:: + parse_parameters (ParameterHandler &prm) + { + end_time = prm.get_double ("End time"); + initial_global_refinement = prm.get_integer ("Initial global refinement"); + initial_adaptive_refinement = prm.get_integer ("Initial adaptive refinement"); - adaptive_refinement_interval= prm.get_integer ("Time steps between mesh refinement"); + adaptive_refinement_interval= prm.get_integer ("Time steps between mesh refinement"); - generate_graphical_output = prm.get_bool ("Generate graphical output"); - graphical_output_interval = prm.get_integer ("Time steps between graphical output"); + generate_graphical_output = prm.get_bool ("Generate graphical output"); + graphical_output_interval = prm.get_integer ("Time steps between graphical output"); - prm.enter_subsection ("Stabilization parameters"); - { - stabilization_alpha = prm.get_double ("alpha"); - stabilization_c_R = prm.get_double ("c_R"); - stabilization_beta = prm.get_double ("beta"); - } - prm.leave_subsection (); + prm.enter_subsection ("Stabilization parameters"); + { + stabilization_alpha = prm.get_double ("alpha"); + stabilization_c_R = prm.get_double ("c_R"); + stabilization_beta = prm.get_double ("beta"); + } + prm.leave_subsection (); - prm.enter_subsection ("Discretization"); - { - stokes_velocity_degree = prm.get_integer ("Stokes velocity polynomial degree"); - temperature_degree = prm.get_integer ("Temperature polynomial degree"); - use_locally_conservative_discretization - = prm.get_bool ("Use locally conservative discretization"); + prm.enter_subsection ("Discretization"); + { + stokes_velocity_degree = prm.get_integer ("Stokes velocity polynomial degree"); + temperature_degree = prm.get_integer ("Temperature polynomial degree"); + use_locally_conservative_discretization + = prm.get_bool ("Use locally conservative discretization"); + } + prm.leave_subsection (); } - prm.leave_subsection (); -} @@ -1248,49 +1250,49 @@ parse_parameters (ParameterHandler &prm) // get a summary table in the end of the // program which shows us wallclock times // (as opposed to CPU times). -template -BoussinesqFlowProblem::BoussinesqFlowProblem (Parameters ¶meters_) - : - parameters (parameters_), - pcout (std::cout, - (Utilities::System:: - get_this_mpi_process(MPI_COMM_WORLD) - == 0)), - - triangulation (MPI_COMM_WORLD, - typename Triangulation::MeshSmoothing - (Triangulation::smoothing_on_refinement | - Triangulation::smoothing_on_coarsening)), - - mapping (4), - - stokes_fe (FE_Q(parameters.stokes_velocity_degree), - dim, - (parameters.use_locally_conservative_discretization - ? - static_cast &> - (FE_DGP(parameters.stokes_velocity_degree-1)) - : - static_cast &> - (FE_Q(parameters.stokes_velocity_degree-1))), - 1), - - stokes_dof_handler (triangulation), - - temperature_fe (parameters.temperature_degree), - temperature_dof_handler (triangulation), - - time_step (0), - old_time_step (0), - timestep_number (0), - rebuild_stokes_matrix (true), - rebuild_stokes_preconditioner (true), - rebuild_temperature_matrices (true), - rebuild_temperature_preconditioner (true), - - computing_timer (pcout, TimerOutput::summary, - TimerOutput::wall_times) -{} + template + BoussinesqFlowProblem::BoussinesqFlowProblem (Parameters ¶meters_) + : + parameters (parameters_), + pcout (std::cout, + (Utilities::System:: + get_this_mpi_process(MPI_COMM_WORLD) + == 0)), + + triangulation (MPI_COMM_WORLD, + typename Triangulation::MeshSmoothing + (Triangulation::smoothing_on_refinement | + Triangulation::smoothing_on_coarsening)), + + mapping (4), + + stokes_fe (FE_Q(parameters.stokes_velocity_degree), + dim, + (parameters.use_locally_conservative_discretization + ? + static_cast &> + (FE_DGP(parameters.stokes_velocity_degree-1)) + : + static_cast &> + (FE_Q(parameters.stokes_velocity_degree-1))), + 1), + + stokes_dof_handler (triangulation), + + temperature_fe (parameters.temperature_degree), + temperature_dof_handler (triangulation), + + time_step (0), + old_time_step (0), + timestep_number (0), + rebuild_stokes_matrix (true), + rebuild_stokes_preconditioner (true), + rebuild_temperature_matrices (true), + rebuild_temperature_preconditioner (true), + + computing_timer (pcout, TimerOutput::summary, + TimerOutput::wall_times) + {} @@ -1342,46 +1344,46 @@ BoussinesqFlowProblem::BoussinesqFlowProblem (Parameters ¶meters_) // global maximum (output), and an // integer value one that says that // we only work on one double. -template -double BoussinesqFlowProblem::get_maximal_velocity () const -{ - const QIterated quadrature_formula (QTrapez<1>(), - parameters.stokes_velocity_degree); - const unsigned int n_q_points = quadrature_formula.size(); - - FEValues fe_values (mapping, stokes_fe, quadrature_formula, update_values); - std::vector > velocity_values(n_q_points); + template + double BoussinesqFlowProblem::get_maximal_velocity () const + { + const QIterated quadrature_formula (QTrapez<1>(), + parameters.stokes_velocity_degree); + const unsigned int n_q_points = quadrature_formula.size(); - const FEValuesExtractors::Vector velocities (0); + FEValues fe_values (mapping, stokes_fe, quadrature_formula, update_values); + std::vector > velocity_values(n_q_points); - double max_local_velocity = 0; + const FEValuesExtractors::Vector velocities (0); - typename DoFHandler::active_cell_iterator - cell = stokes_dof_handler.begin_active(), - endc = stokes_dof_handler.end(); - for (; cell!=endc; ++cell) - if (cell->subdomain_id() == - Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)) - { - fe_values.reinit (cell); - fe_values[velocities].get_function_values (stokes_solution, - velocity_values); - - for (unsigned int q=0; q::active_cell_iterator + cell = stokes_dof_handler.begin_active(), + endc = stokes_dof_handler.end(); + for (; cell!=endc; ++cell) + if (cell->subdomain_id() == + Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)) + { + fe_values.reinit (cell); + fe_values[velocities].get_function_values (stokes_solution, + velocity_values); + + for (unsigned int q=0; q::get_maximal_velocity () const // compute the cfl number, i.e., maximal // velocity on a cell divided by the cell // diameter -template -double BoussinesqFlowProblem::get_cfl_number () const -{ - const QIterated quadrature_formula (QTrapez<1>(), - parameters.stokes_velocity_degree); - const unsigned int n_q_points = quadrature_formula.size(); - - FEValues fe_values (mapping, stokes_fe, quadrature_formula, update_values); - std::vector > velocity_values(n_q_points); + template + double BoussinesqFlowProblem::get_cfl_number () const + { + const QIterated quadrature_formula (QTrapez<1>(), + parameters.stokes_velocity_degree); + const unsigned int n_q_points = quadrature_formula.size(); - const FEValuesExtractors::Vector velocities (0); + FEValues fe_values (mapping, stokes_fe, quadrature_formula, update_values); + std::vector > velocity_values(n_q_points); - double max_local_cfl = 0; + const FEValuesExtractors::Vector velocities (0); - typename DoFHandler::active_cell_iterator - cell = stokes_dof_handler.begin_active(), - endc = stokes_dof_handler.end(); - for (; cell!=endc; ++cell) - if (cell->subdomain_id() == - Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)) - { - fe_values.reinit (cell); - fe_values[velocities].get_function_values (stokes_solution, - velocity_values); - - double max_local_velocity = 1e-10; - for (unsigned int q=0; qdiameter()); - } + double max_local_cfl = 0; - double max_cfl_number = 0.; + typename DoFHandler::active_cell_iterator + cell = stokes_dof_handler.begin_active(), + endc = stokes_dof_handler.end(); + for (; cell!=endc; ++cell) + if (cell->subdomain_id() == + Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)) + { + fe_values.reinit (cell); + fe_values[velocities].get_function_values (stokes_solution, + velocity_values); + + double max_local_velocity = 1e-10; + for (unsigned int q=0; qdiameter()); + } + + double max_cfl_number = 0.; #ifdef DEAL_II_COMPILER_SUPPORTS_MPI - MPI_Allreduce (&max_local_cfl, &max_cfl_number, 1, MPI_DOUBLE, - MPI_MAX, MPI_COMM_WORLD); + MPI_Allreduce (&max_local_cfl, &max_cfl_number, 1, MPI_DOUBLE, + MPI_MAX, MPI_COMM_WORLD); #else - max_cfl_number = max_local_cfl; + max_cfl_number = max_local_cfl; #endif - return max_cfl_number; -} + return max_cfl_number; + } -template -double -BoussinesqFlowProblem::get_entropy_variation (const double average_temperature) const -{ - // only do this if we really need entropy - // variation - if (parameters.stabilization_alpha != 2) - return 1.; - - // record maximal entropy on Gauss quadrature - // points - const QGauss quadrature_formula (parameters.temperature_degree+1); - const unsigned int n_q_points = quadrature_formula.size(); - - FEValues fe_values (temperature_fe, quadrature_formula, - update_values | update_JxW_values); - std::vector old_temperature_values(n_q_points); - std::vector old_old_temperature_values(n_q_points); - - double min_entropy = std::numeric_limits::max(), - max_entropy = -std::numeric_limits::max(), - area = 0, - entropy_integrated = 0; - - typename DoFHandler::active_cell_iterator - cell = temperature_dof_handler.begin_active(), - endc = temperature_dof_handler.end(); - for (; cell!=endc; ++cell) - if (cell->subdomain_id() == - Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)) - { - fe_values.reinit (cell); - 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 + double + BoussinesqFlowProblem::get_entropy_variation (const double average_temperature) const + { + // only do this if we really need entropy + // variation + if (parameters.stabilization_alpha != 2) + return 1.; + + // record maximal entropy on Gauss quadrature + // points + const QGauss quadrature_formula (parameters.temperature_degree+1); + const unsigned int n_q_points = quadrature_formula.size(); + + FEValues fe_values (temperature_fe, quadrature_formula, + update_values | update_JxW_values); + std::vector old_temperature_values(n_q_points); + std::vector old_old_temperature_values(n_q_points); + + double min_entropy = std::numeric_limits::max(), + max_entropy = -std::numeric_limits::max(), + area = 0, + entropy_integrated = 0; - // do MPI data exchange: we need to sum over - // the two integrals (area, - // entropy_integrated), and get the extrema - // for maximum and minimum. combine - // MPI_Allreduce for two values since that is - // an expensive operation - double local_for_sum[2], global_for_sum[2]; - double local_for_max[2], global_for_max[2]; - local_for_sum[0] = entropy_integrated; - local_for_sum[1] = area; - local_for_max[0] = -min_entropy; - local_for_max[1] = max_entropy; + typename DoFHandler::active_cell_iterator + cell = temperature_dof_handler.begin_active(), + endc = temperature_dof_handler.end(); + for (; cell!=endc; ++cell) + if (cell->subdomain_id() == + Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)) + { + fe_values.reinit (cell); + 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::get_entropy_variation (const double average_temperat // partition. Two global communication // steps synchronize the data among the // processors. -template -std::pair -BoussinesqFlowProblem::get_extrapolated_temperature_range () const -{ - const QIterated quadrature_formula (QTrapez<1>(), - parameters.temperature_degree); - const unsigned int n_q_points = quadrature_formula.size(); - - FEValues fe_values (mapping, temperature_fe, quadrature_formula, - update_values); - std::vector old_temperature_values(n_q_points); - std::vector old_old_temperature_values(n_q_points); - - // This presets the minimum with a bigger - // and the maximum with a smaller number - // than one that is going to appear. Will - // be overwritten in the cell loop or in - // the communication step at the - // latest. - double min_local_temperature = std::numeric_limits::max(), - max_local_temperature = -std::numeric_limits::max(); - - if (timestep_number != 0) - { - typename DoFHandler::active_cell_iterator - cell = temperature_dof_handler.begin_active(), - endc = temperature_dof_handler.end(); - for (; cell!=endc; ++cell) - if (cell->subdomain_id() == - Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)) - { - fe_values.reinit (cell); - 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::active_cell_iterator - cell = temperature_dof_handler.begin_active(), - endc = temperature_dof_handler.end(); - for (; cell!=endc; ++cell) - if (cell->subdomain_id() == - Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)) - { - fe_values.reinit (cell); - fe_values.get_function_values (old_temperature_solution, - old_temperature_values); - - for (unsigned int q=0; q + std::pair + BoussinesqFlowProblem::get_extrapolated_temperature_range () const + { + const QIterated quadrature_formula (QTrapez<1>(), + parameters.temperature_degree); + const unsigned int n_q_points = quadrature_formula.size(); + + FEValues fe_values (mapping, temperature_fe, quadrature_formula, + update_values); + std::vector old_temperature_values(n_q_points); + std::vector old_old_temperature_values(n_q_points); + + // This presets the minimum with a bigger + // and the maximum with a smaller number + // than one that is going to appear. Will + // be overwritten in the cell loop or in + // the communication step at the + // latest. + double min_local_temperature = std::numeric_limits::max(), + max_local_temperature = -std::numeric_limits::max(); + + if (timestep_number != 0) + { + typename DoFHandler::active_cell_iterator + cell = temperature_dof_handler.begin_active(), + endc = temperature_dof_handler.end(); + for (; cell!=endc; ++cell) + if (cell->subdomain_id() == + Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)) + { + fe_values.reinit (cell); + 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::active_cell_iterator + cell = temperature_dof_handler.begin_active(), + endc = temperature_dof_handler.end(); + for (; cell!=endc; ++cell) + if (cell->subdomain_id() == + Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)) + { + fe_values.reinit (cell); + fe_values.get_function_values (old_temperature_solution, + old_temperature_values); + + for (unsigned int q=0; q -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 > &old_velocity_values, - const std::vector > &old_old_velocity_values, - const std::vector > &old_strain_rates, - const std::vector > &old_old_strain_rates, - const double global_u_infty, - const double global_T_variation, - const double average_temperature, - const double global_entropy_variation, - const double cell_diameter) const -{ - if (global_u_infty == 0) - return 5e-3 * cell_diameter; + 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 > &old_velocity_values, + const std::vector > &old_old_velocity_values, + const std::vector > &old_strain_rates, + const std::vector > &old_old_strain_rates, + const double global_u_infty, + const double global_T_variation, + const double average_temperature, + const double global_entropy_variation, + const double cell_diameter) const + { + if (global_u_infty == 0) + return 5e-3 * cell_diameter; - const unsigned int n_q_points = old_temperature.size(); + const unsigned int n_q_points = old_temperature.size(); - double max_residual = 0; - double max_velocity = 0; + double max_residual = 0; + double max_velocity = 0; - for (unsigned int q=0; q < n_q_points; ++q) - { - const Tensor<1,dim> u = (old_velocity_values[q] + - old_old_velocity_values[q]) / 2; - - const SymmetricTensor<2,dim> strain_rate = (old_strain_rates[q] + - old_old_strain_rates[q]) / 2; - - const double T = (old_temperature[q] + old_old_temperature[q]) / 2; - const double dT_dt = (old_temperature[q] - old_old_temperature[q]) - / old_time_step; - const double u_grad_T = u * (old_temperature_grads[q] + - old_old_temperature_grads[q]) / 2; - - const double kappa_Delta_T = EquationData::kappa - * (old_temperature_laplacians[q] + - old_old_temperature_laplacians[q]) / 2; - const double gamma - = ((EquationData::radiogenic_heating * EquationData::density(T) - + - 2 * EquationData::eta * strain_rate * strain_rate) / - (EquationData::density(T) * EquationData::specific_heat)); - - double residual - = std::abs(dT_dt + u_grad_T - kappa_Delta_T - gamma); - if (parameters.stabilization_alpha == 2) - residual *= std::abs(T - average_temperature); - - max_residual = std::max (residual, max_residual); - max_velocity = std::max (std::sqrt (u*u), max_velocity); - } - - const double max_viscosity = (parameters.stabilization_beta * - max_velocity * cell_diameter); - if (timestep_number == 0) - return max_viscosity; - else - { - Assert (old_time_step > 0, ExcInternalError()); - - double entropy_viscosity; - if (parameters.stabilization_alpha == 2) - entropy_viscosity = (parameters.stabilization_c_R * - cell_diameter * cell_diameter * - max_residual / - global_entropy_variation); - else - entropy_viscosity = (parameters.stabilization_c_R * - cell_diameter * global_Omega_diameter * - max_velocity * max_residual / - (global_u_infty * global_T_variation)); + for (unsigned int q=0; q < n_q_points; ++q) + { + const Tensor<1,dim> u = (old_velocity_values[q] + + old_old_velocity_values[q]) / 2; + + const SymmetricTensor<2,dim> strain_rate = (old_strain_rates[q] + + old_old_strain_rates[q]) / 2; + + const double T = (old_temperature[q] + old_old_temperature[q]) / 2; + const double dT_dt = (old_temperature[q] - old_old_temperature[q]) + / old_time_step; + const double u_grad_T = u * (old_temperature_grads[q] + + old_old_temperature_grads[q]) / 2; + + const double kappa_Delta_T = EquationData::kappa + * (old_temperature_laplacians[q] + + old_old_temperature_laplacians[q]) / 2; + const double gamma + = ((EquationData::radiogenic_heating * EquationData::density(T) + + + 2 * EquationData::eta * strain_rate * strain_rate) / + (EquationData::density(T) * EquationData::specific_heat)); + + double residual + = std::abs(dT_dt + u_grad_T - kappa_Delta_T - gamma); + if (parameters.stabilization_alpha == 2) + residual *= std::abs(T - average_temperature); + + max_residual = std::max (residual, max_residual); + max_velocity = std::max (std::sqrt (u*u), max_velocity); + } - return std::min (max_viscosity, entropy_viscosity); - } -} + const double max_viscosity = (parameters.stabilization_beta * + max_velocity * cell_diameter); + if (timestep_number == 0) + return max_viscosity; + else + { + Assert (old_time_step > 0, ExcInternalError()); + + double entropy_viscosity; + if (parameters.stabilization_alpha == 2) + entropy_viscosity = (parameters.stabilization_c_R * + cell_diameter * cell_diameter * + max_residual / + global_entropy_variation); + else + entropy_viscosity = (parameters.stabilization_c_R * + cell_diameter * global_Omega_diameter * + max_velocity * max_residual / + (global_u_infty * global_T_variation)); + + return std::min (max_viscosity, entropy_viscosity); + } + } @@ -1788,116 +1790,116 @@ compute_viscosity (const std::vector &old_temperature, // exactly the entry that we would fill the // matrix with — otherwise, this will // not work. -template -void BoussinesqFlowProblem::project_temperature_field () -{ - assemble_temperature_matrix (); + template + void BoussinesqFlowProblem::project_temperature_field () + { + assemble_temperature_matrix (); - QGauss quadrature(parameters.temperature_degree+2); - UpdateFlags update_flags = UpdateFlags(update_values | - update_quadrature_points | - update_JxW_values); - FEValues fe_values (mapping, temperature_fe, quadrature, update_flags); + QGauss quadrature(parameters.temperature_degree+2); + UpdateFlags update_flags = UpdateFlags(update_values | + update_quadrature_points | + update_JxW_values); + FEValues fe_values (mapping, temperature_fe, quadrature, update_flags); - const unsigned int dofs_per_cell = fe_values.dofs_per_cell, - n_q_points = fe_values.n_quadrature_points; + const unsigned int dofs_per_cell = fe_values.dofs_per_cell, + n_q_points = fe_values.n_quadrature_points; - std::vector local_dof_indices (dofs_per_cell); - Vector cell_vector (dofs_per_cell); - FullMatrix matrix_for_bc (dofs_per_cell, dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); + Vector cell_vector (dofs_per_cell); + FullMatrix matrix_for_bc (dofs_per_cell, dofs_per_cell); - std::vector rhs_values(n_q_points); + std::vector rhs_values(n_q_points); - TrilinosWrappers::MPI::Vector - rhs (temperature_mass_matrix.row_partitioner()), + TrilinosWrappers::MPI::Vector + rhs (temperature_mass_matrix.row_partitioner()), solution (temperature_mass_matrix.row_partitioner()); - typename DoFHandler::active_cell_iterator - cell = temperature_dof_handler.begin_active(), - endc = temperature_dof_handler.end(); - - for (; cell!=endc; ++cell) - if (cell->subdomain_id() == - Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)) - { - cell->get_dof_indices (local_dof_indices); - fe_values.reinit(cell); - - EquationData::TemperatureInitialValues().value_list - (fe_values.get_quadrature_points(), rhs_values); - - cell_vector = 0; - matrix_for_bc = 0; - for (unsigned int point=0; point::active_cell_iterator + cell = temperature_dof_handler.begin_active(), + endc = temperature_dof_handler.end(); - rhs.compress (Add); - - SolverControl solver_control(5*rhs.size(), 1e-12*rhs.l2_norm()); - SolverCG cg(solver_control); - - TrilinosWrappers::PreconditionJacobi preconditioner_mass; - preconditioner_mass.initialize(temperature_mass_matrix, 1.3); - - cg.solve (temperature_mass_matrix, solution, rhs, preconditioner_mass); - - temperature_constraints.distribute (solution); - // Having so computed the current - // temperature field, let us set - // the member variable that holds - // the temperature nodes. Strictly - // speaking, we really only need to - // set - // old_temperature_solution - // since the first thing we will do - // is to compute the Stokes - // solution that only requires the - // previous time step's temperature - // field. That said, nothing good - // can come from not initializing - // the other vectors as well - // (especially since it's a - // relatively cheap operation and - // we only have to do it once at - // the beginning of the program) if - // we ever want to extend our - // numerical method or physical - // model, and so we initialize - // temperature_solution - // and - // old_old_temperature_solution - // as well. The second and third - // argument to the reinit function - // indicates that the elements of - // the first argument shall be - // copied and that while the - // solution vector is - // strictly distributed, the - // vectors here initialized can - // (and indeed do) have ghost - // elements: - temperature_solution.reinit(solution, false, true); - old_temperature_solution.reinit(solution, false, true); - old_old_temperature_solution.reinit(solution, false, true); -} + for (; cell!=endc; ++cell) + if (cell->subdomain_id() == + Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)) + { + cell->get_dof_indices (local_dof_indices); + fe_values.reinit(cell); + + EquationData::TemperatureInitialValues().value_list + (fe_values.get_quadrature_points(), rhs_values); + + cell_vector = 0; + matrix_for_bc = 0; + for (unsigned int point=0; point cg(solver_control); + + TrilinosWrappers::PreconditionJacobi preconditioner_mass; + preconditioner_mass.initialize(temperature_mass_matrix, 1.3); + + cg.solve (temperature_mass_matrix, solution, rhs, preconditioner_mass); + + temperature_constraints.distribute (solution); + // Having so computed the current + // temperature field, let us set + // the member variable that holds + // the temperature nodes. Strictly + // speaking, we really only need to + // set + // old_temperature_solution + // since the first thing we will do + // is to compute the Stokes + // solution that only requires the + // previous time step's temperature + // field. That said, nothing good + // can come from not initializing + // the other vectors as well + // (especially since it's a + // relatively cheap operation and + // we only have to do it once at + // the beginning of the program) if + // we ever want to extend our + // numerical method or physical + // model, and so we initialize + // temperature_solution + // and + // old_old_temperature_solution + // as well. The second and third + // argument to the reinit function + // indicates that the elements of + // the first argument shall be + // copied and that while the + // solution vector is + // strictly distributed, the + // vectors here initialized can + // (and indeed do) have ghost + // elements: + temperature_solution.reinit(solution, false, true); + old_temperature_solution.reinit(solution, false, true); + old_old_temperature_solution.reinit(solution, false, true); + } @@ -1979,88 +1981,88 @@ void BoussinesqFlowProblem::project_temperature_field () // patterns, so we can safely release the // sp variable once the matrix // has been given the sparsity structure. -template -void BoussinesqFlowProblem:: -setup_stokes_matrix (const std::vector &stokes_partitioning) -{ - stokes_matrix.clear (); - - TrilinosWrappers::BlockSparsityPattern sp (stokes_partitioning, - MPI_COMM_WORLD); - - Table<2,DoFTools::Coupling> coupling (dim+1, dim+1); - - for (unsigned int c=0; c + void BoussinesqFlowProblem:: + setup_stokes_matrix (const std::vector &stokes_partitioning) + { + stokes_matrix.clear (); - DoFTools::make_sparsity_pattern (stokes_dof_handler, - coupling, sp, - stokes_constraints, false, - Utilities::System:: - get_this_mpi_process(MPI_COMM_WORLD)); - sp.compress(); + TrilinosWrappers::BlockSparsityPattern sp (stokes_partitioning, + MPI_COMM_WORLD); - stokes_matrix.reinit (sp); -} + Table<2,DoFTools::Coupling> coupling (dim+1, dim+1); + for (unsigned int c=0; c -void BoussinesqFlowProblem:: -setup_stokes_preconditioner (const std::vector &stokes_partitioning) -{ - Amg_preconditioner.reset (); - Mp_preconditioner.reset (); - - stokes_preconditioner_matrix.clear (); - - TrilinosWrappers::BlockSparsityPattern sp (stokes_partitioning, - MPI_COMM_WORLD); + stokes_matrix.reinit (sp); + } - Table<2,DoFTools::Coupling> coupling (dim+1, dim+1); - for (unsigned int c=0; c + void BoussinesqFlowProblem:: + setup_stokes_preconditioner (const std::vector &stokes_partitioning) + { + Amg_preconditioner.reset (); + Mp_preconditioner.reset (); + + stokes_preconditioner_matrix.clear (); + + TrilinosWrappers::BlockSparsityPattern sp (stokes_partitioning, + MPI_COMM_WORLD); + + Table<2,DoFTools::Coupling> coupling (dim+1, dim+1); + for (unsigned int c=0; c -void BoussinesqFlowProblem:: -setup_temperature_matrices (const IndexSet &temperature_partitioner) -{ - T_preconditioner.reset (); - temperature_mass_matrix.clear (); - temperature_stiffness_matrix.clear (); - temperature_matrix.clear (); - - TrilinosWrappers::SparsityPattern sp (temperature_partitioner, - MPI_COMM_WORLD); - DoFTools::make_sparsity_pattern (temperature_dof_handler, sp, - temperature_constraints, false, - Utilities::System:: - get_this_mpi_process(MPI_COMM_WORLD)); - sp.compress(); - - temperature_matrix.reinit (sp); - temperature_mass_matrix.reinit (sp); - temperature_stiffness_matrix.reinit (sp); -} + template + void BoussinesqFlowProblem:: + setup_temperature_matrices (const IndexSet &temperature_partitioner) + { + T_preconditioner.reset (); + temperature_mass_matrix.clear (); + temperature_stiffness_matrix.clear (); + temperature_matrix.clear (); + + TrilinosWrappers::SparsityPattern sp (temperature_partitioner, + MPI_COMM_WORLD); + DoFTools::make_sparsity_pattern (temperature_dof_handler, sp, + temperature_constraints, false, + Utilities::System:: + get_this_mpi_process(MPI_COMM_WORLD)); + sp.compress(); + + temperature_matrix.reinit (sp); + temperature_mass_matrix.reinit (sp); + temperature_stiffness_matrix.reinit (sp); + } @@ -2142,130 +2144,130 @@ setup_temperature_matrices (const IndexSet &temperature_partitioner) // get_this_mpi_process // (MPI_COMM_WORLD) == 0), // hardly a pretty solution. -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; - stokes_dof_handler.distribute_dofs (stokes_fe); - DoFRenumbering::component_wise (stokes_dof_handler, stokes_sub_blocks); - - temperature_dof_handler.distribute_dofs (temperature_fe); - - std::vector stokes_dofs_per_block (2); - DoFTools::count_dofs_per_block (stokes_dof_handler, stokes_dofs_per_block, - stokes_sub_blocks); - - const unsigned int n_u = stokes_dofs_per_block[0], - n_p = stokes_dofs_per_block[1], - n_T = temperature_dof_handler.n_dofs(); - - // print dof numbers with 1000s - // separator since they are frequently - // large - std::locale s = pcout.get_stream().getloc(); - pcout.get_stream().imbue(std::locale("")); - pcout << "Number of active cells: " - << triangulation.n_global_active_cells() - << " (on " - << triangulation.n_levels() - << " levels)" - << std::endl - << "Number of degrees of freedom: " - << n_u + n_p + n_T - << " (" << n_u << '+' << n_p << '+'<< n_T <<')' - << std::endl - << std::endl; - pcout.get_stream().imbue(s); - - - - std::vector stokes_partitioning, stokes_relevant_partitioning; - IndexSet temperature_partitioning (n_T), temperature_relevant_partitioning (n_T); - IndexSet stokes_relevant_set; + template + void BoussinesqFlowProblem::setup_dofs () { - const unsigned int my_id = - Utilities::System::get_this_mpi_process(MPI_COMM_WORLD); - IndexSet stokes_index_set = stokes_dof_handler.locally_owned_dofs(); - stokes_partitioning.push_back(stokes_index_set.get_view(0,n_u)); - stokes_partitioning.push_back(stokes_index_set.get_view(n_u,n_u+n_p)); - - DoFTools::extract_locally_relevant_dofs (stokes_dof_handler, - stokes_relevant_set); - stokes_relevant_partitioning.push_back(stokes_relevant_set.get_view(0,n_u)); - stokes_relevant_partitioning.push_back(stokes_relevant_set.get_view(n_u,n_u+n_p)); - - temperature_partitioning = temperature_dof_handler.locally_owned_dofs(); - DoFTools::extract_locally_relevant_dofs (temperature_dof_handler, - temperature_relevant_partitioning); - } + computing_timer.enter_section("Setup dof systems"); + + std::vector stokes_sub_blocks (dim+1,0); + stokes_sub_blocks[dim] = 1; + stokes_dof_handler.distribute_dofs (stokes_fe); + DoFRenumbering::component_wise (stokes_dof_handler, stokes_sub_blocks); + + temperature_dof_handler.distribute_dofs (temperature_fe); + + std::vector stokes_dofs_per_block (2); + DoFTools::count_dofs_per_block (stokes_dof_handler, stokes_dofs_per_block, + stokes_sub_blocks); + + const unsigned int n_u = stokes_dofs_per_block[0], + n_p = stokes_dofs_per_block[1], + n_T = temperature_dof_handler.n_dofs(); + + // print dof numbers with 1000s + // separator since they are frequently + // large + std::locale s = pcout.get_stream().getloc(); + pcout.get_stream().imbue(std::locale("")); + pcout << "Number of active cells: " + << triangulation.n_global_active_cells() + << " (on " + << triangulation.n_levels() + << " levels)" + << std::endl + << "Number of degrees of freedom: " + << n_u + n_p + n_T + << " (" << n_u << '+' << n_p << '+'<< n_T <<')' + << std::endl + << std::endl; + pcout.get_stream().imbue(s); + + + + std::vector stokes_partitioning, stokes_relevant_partitioning; + IndexSet temperature_partitioning (n_T), temperature_relevant_partitioning (n_T); + IndexSet stokes_relevant_set; + { + const unsigned int my_id = + Utilities::System::get_this_mpi_process(MPI_COMM_WORLD); + IndexSet stokes_index_set = stokes_dof_handler.locally_owned_dofs(); + stokes_partitioning.push_back(stokes_index_set.get_view(0,n_u)); + stokes_partitioning.push_back(stokes_index_set.get_view(n_u,n_u+n_p)); + + DoFTools::extract_locally_relevant_dofs (stokes_dof_handler, + stokes_relevant_set); + stokes_relevant_partitioning.push_back(stokes_relevant_set.get_view(0,n_u)); + stokes_relevant_partitioning.push_back(stokes_relevant_set.get_view(n_u,n_u+n_p)); + + temperature_partitioning = temperature_dof_handler.locally_owned_dofs(); + DoFTools::extract_locally_relevant_dofs (temperature_dof_handler, + temperature_relevant_partitioning); + } - { + { - stokes_constraints.clear (); + stokes_constraints.clear (); // IndexSet stokes_la; // DoFTools::extract_locally_active_dofs (stokes_dof_handler, // stokes_la); - stokes_constraints.reinit(stokes_relevant_set); - - DoFTools::make_hanging_node_constraints (stokes_dof_handler, - stokes_constraints); - - std::vector velocity_mask (dim+1, true); - velocity_mask[dim] = false; - VectorTools::interpolate_boundary_values (stokes_dof_handler, - 0, - ZeroFunction(dim+1), - stokes_constraints, - velocity_mask); - - std::set no_normal_flux_boundaries; - no_normal_flux_boundaries.insert (1); - VectorTools::compute_no_normal_flux_constraints (stokes_dof_handler, 0, - no_normal_flux_boundaries, - stokes_constraints, - mapping); - stokes_constraints.close (); - } - { - temperature_constraints.clear (); - temperature_constraints.reinit(temperature_relevant_partitioning);//temp_locally_active); - - DoFTools::make_hanging_node_constraints (temperature_dof_handler, - temperature_constraints); - VectorTools::interpolate_boundary_values (temperature_dof_handler, - 0, - EquationData::TemperatureInitialValues(), - temperature_constraints); - VectorTools::interpolate_boundary_values (temperature_dof_handler, - 1, - EquationData::TemperatureInitialValues(), - temperature_constraints); - temperature_constraints.close (); - } + stokes_constraints.reinit(stokes_relevant_set); + + DoFTools::make_hanging_node_constraints (stokes_dof_handler, + stokes_constraints); + + std::vector velocity_mask (dim+1, true); + velocity_mask[dim] = false; + VectorTools::interpolate_boundary_values (stokes_dof_handler, + 0, + ZeroFunction(dim+1), + stokes_constraints, + velocity_mask); + + std::set no_normal_flux_boundaries; + no_normal_flux_boundaries.insert (1); + VectorTools::compute_no_normal_flux_constraints (stokes_dof_handler, 0, + no_normal_flux_boundaries, + stokes_constraints, + mapping); + stokes_constraints.close (); + } + { + temperature_constraints.clear (); + temperature_constraints.reinit(temperature_relevant_partitioning);//temp_locally_active); + + DoFTools::make_hanging_node_constraints (temperature_dof_handler, + temperature_constraints); + VectorTools::interpolate_boundary_values (temperature_dof_handler, + 0, + EquationData::TemperatureInitialValues(), + temperature_constraints); + VectorTools::interpolate_boundary_values (temperature_dof_handler, + 1, + EquationData::TemperatureInitialValues(), + temperature_constraints); + temperature_constraints.close (); + } - setup_stokes_matrix (stokes_partitioning); - setup_stokes_preconditioner (stokes_partitioning); - setup_temperature_matrices (temperature_partitioning); + setup_stokes_matrix (stokes_partitioning); + setup_stokes_preconditioner (stokes_partitioning); + setup_temperature_matrices (temperature_partitioning); - stokes_rhs.reinit (stokes_partitioning, MPI_COMM_WORLD); - stokes_solution.reinit (stokes_relevant_partitioning, MPI_COMM_WORLD); - old_stokes_solution.reinit (stokes_solution); + stokes_rhs.reinit (stokes_partitioning, MPI_COMM_WORLD); + stokes_solution.reinit (stokes_relevant_partitioning, MPI_COMM_WORLD); + old_stokes_solution.reinit (stokes_solution); - temperature_rhs.reinit (temperature_partitioning, MPI_COMM_WORLD); - temperature_solution.reinit (temperature_relevant_partitioning, MPI_COMM_WORLD); - old_temperature_solution.reinit (temperature_solution); - old_old_temperature_solution.reinit (temperature_solution); + temperature_rhs.reinit (temperature_partitioning, MPI_COMM_WORLD); + temperature_solution.reinit (temperature_relevant_partitioning, MPI_COMM_WORLD); + old_temperature_solution.reinit (temperature_solution); + old_old_temperature_solution.reinit (temperature_solution); - rebuild_stokes_matrix = true; - rebuild_stokes_preconditioner = true; - rebuild_temperature_matrices = true; - rebuild_temperature_preconditioner = true; + rebuild_stokes_matrix = true; + rebuild_stokes_preconditioner = true; + rebuild_temperature_matrices = true; + rebuild_temperature_preconditioner = true; - computing_timer.exit_section(); -} + computing_timer.exit_section(); + } @@ -2320,60 +2322,60 @@ void BoussinesqFlowProblem::setup_dofs () // scratch data object, which is allocated // only once before we start the loop over // all cells: -template -void -BoussinesqFlowProblem:: -local_assemble_stokes_preconditioner (const typename DoFHandler::active_cell_iterator &cell, - Assembly::Scratch::StokesPreconditioner &scratch, - Assembly::CopyData::StokesPreconditioner &data) -{ - const unsigned int dofs_per_cell = stokes_fe.dofs_per_cell; - const unsigned int n_q_points = scratch.stokes_fe_values.n_quadrature_points; + template + void + BoussinesqFlowProblem:: + local_assemble_stokes_preconditioner (const typename DoFHandler::active_cell_iterator &cell, + Assembly::Scratch::StokesPreconditioner &scratch, + Assembly::CopyData::StokesPreconditioner &data) + { + const unsigned int dofs_per_cell = stokes_fe.dofs_per_cell; + const unsigned int n_q_points = scratch.stokes_fe_values.n_quadrature_points; - const FEValuesExtractors::Vector velocities (0); - const FEValuesExtractors::Scalar pressure (dim); + const FEValuesExtractors::Vector velocities (0); + const FEValuesExtractors::Scalar pressure (dim); - scratch.stokes_fe_values.reinit (cell); - cell->get_dof_indices (data.local_dof_indices); + scratch.stokes_fe_values.reinit (cell); + cell->get_dof_indices (data.local_dof_indices); - data.local_matrix = 0; + data.local_matrix = 0; - for (unsigned int q=0; q -void -BoussinesqFlowProblem:: -copy_local_to_global_stokes_preconditioner (const Assembly::CopyData::StokesPreconditioner &data) -{ - stokes_constraints.distribute_local_to_global (data.local_matrix, - data.local_dof_indices, - stokes_preconditioner_matrix); -} + template + void + BoussinesqFlowProblem:: + copy_local_to_global_stokes_preconditioner (const Assembly::CopyData::StokesPreconditioner &data) + { + stokes_constraints.distribute_local_to_global (data.local_matrix, + data.local_dof_indices, + stokes_preconditioner_matrix); + } @@ -2429,46 +2431,46 @@ copy_local_to_global_stokes_preconditioner (const Assembly::CopyData::StokesPrec // send it to the owner at the end of // assembly, by calling the // compress() command. -template -void -BoussinesqFlowProblem::assemble_stokes_preconditioner () -{ - stokes_preconditioner_matrix = 0; - - const QGauss quadrature_formula(parameters.stokes_velocity_degree+1); - - typedef - FilteredIterator::active_cell_iterator> - SubdomainFilter; - - WorkStream:: - run (SubdomainFilter (IteratorFilters::SubdomainEqualTo - (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)), - stokes_dof_handler.begin_active()), - SubdomainFilter (IteratorFilters::SubdomainEqualTo - (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)), - stokes_dof_handler.end()), - std_cxx1x::bind (&BoussinesqFlowProblem:: - local_assemble_stokes_preconditioner, - this, - std_cxx1x::_1, - std_cxx1x::_2, - std_cxx1x::_3), - std_cxx1x::bind (&BoussinesqFlowProblem:: - copy_local_to_global_stokes_preconditioner, - this, - std_cxx1x::_1), - Assembly::Scratch:: - StokesPreconditioner (stokes_fe, quadrature_formula, - mapping, - update_JxW_values | - update_values | - update_gradients), - Assembly::CopyData:: - StokesPreconditioner (stokes_fe)); - - stokes_preconditioner_matrix.compress(); -} + template + void + BoussinesqFlowProblem::assemble_stokes_preconditioner () + { + stokes_preconditioner_matrix = 0; + + const QGauss quadrature_formula(parameters.stokes_velocity_degree+1); + + typedef + FilteredIterator::active_cell_iterator> + SubdomainFilter; + + WorkStream:: + run (SubdomainFilter (IteratorFilters::SubdomainEqualTo + (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)), + stokes_dof_handler.begin_active()), + SubdomainFilter (IteratorFilters::SubdomainEqualTo + (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)), + stokes_dof_handler.end()), + std_cxx1x::bind (&BoussinesqFlowProblem:: + local_assemble_stokes_preconditioner, + this, + std_cxx1x::_1, + std_cxx1x::_2, + std_cxx1x::_3), + std_cxx1x::bind (&BoussinesqFlowProblem:: + copy_local_to_global_stokes_preconditioner, + this, + std_cxx1x::_1), + Assembly::Scratch:: + StokesPreconditioner (stokes_fe, quadrature_formula, + mapping, + update_JxW_values | + update_values | + update_gradients), + Assembly::CopyData:: + StokesPreconditioner (stokes_fe)); + + stokes_preconditioner_matrix.compress(); + } @@ -2481,45 +2483,45 @@ BoussinesqFlowProblem::assemble_stokes_preconditioner () // preconditioner for the pressure mass // matrix instead of IC, as discussed in the // introduction. -template -void -BoussinesqFlowProblem::build_stokes_preconditioner () -{ - if (rebuild_stokes_preconditioner == false) - return; + template + void + BoussinesqFlowProblem::build_stokes_preconditioner () + { + if (rebuild_stokes_preconditioner == false) + return; - computing_timer.enter_section (" Build Stokes preconditioner"); - pcout << " Rebuilding Stokes preconditioner..." << std::flush; + computing_timer.enter_section (" Build Stokes preconditioner"); + pcout << " Rebuilding Stokes preconditioner..." << std::flush; - assemble_stokes_preconditioner (); + assemble_stokes_preconditioner (); - std::vector > constant_modes; - std::vector velocity_components (dim+1,true); - velocity_components[dim] = false; - DoFTools::extract_constant_modes (stokes_dof_handler, velocity_components, - constant_modes); + std::vector > constant_modes; + std::vector velocity_components (dim+1,true); + velocity_components[dim] = false; + DoFTools::extract_constant_modes (stokes_dof_handler, velocity_components, + constant_modes); - Mp_preconditioner = std_cxx1x::shared_ptr - (new TrilinosWrappers::PreconditionILU()); - Amg_preconditioner = std_cxx1x::shared_ptr - (new TrilinosWrappers::PreconditionAMG()); + Mp_preconditioner = std_cxx1x::shared_ptr + (new TrilinosWrappers::PreconditionILU()); + Amg_preconditioner = std_cxx1x::shared_ptr + (new TrilinosWrappers::PreconditionAMG()); - TrilinosWrappers::PreconditionAMG::AdditionalData Amg_data; - Amg_data.constant_modes = constant_modes; - Amg_data.elliptic = true; - Amg_data.higher_order_elements = true; - Amg_data.smoother_sweeps = 2; - Amg_data.aggregation_threshold = 0.02; + TrilinosWrappers::PreconditionAMG::AdditionalData Amg_data; + Amg_data.constant_modes = constant_modes; + Amg_data.elliptic = true; + Amg_data.higher_order_elements = true; + Amg_data.smoother_sweeps = 2; + Amg_data.aggregation_threshold = 0.02; - Mp_preconditioner->initialize (stokes_preconditioner_matrix.block(1,1)); - Amg_preconditioner->initialize (stokes_preconditioner_matrix.block(0,0), - Amg_data); + Mp_preconditioner->initialize (stokes_preconditioner_matrix.block(1,1)); + Amg_preconditioner->initialize (stokes_preconditioner_matrix.block(0,0), + Amg_data); - rebuild_stokes_preconditioner = false; + rebuild_stokes_preconditioner = false; - pcout << std::endl; - computing_timer.exit_section(); -} + pcout << std::endl; + computing_timer.exit_section(); + } // @sect5{Stokes system assembly} @@ -2542,152 +2544,152 @@ BoussinesqFlowProblem::build_stokes_preconditioner () // functions in the end of the assembly in // order to send non-local data to the // owner process. -template -void -BoussinesqFlowProblem:: -local_assemble_stokes_system (const typename DoFHandler::active_cell_iterator &cell, - Assembly::Scratch::StokesSystem &scratch, - Assembly::CopyData::StokesSystem &data) -{ - const unsigned int dofs_per_cell = scratch.stokes_fe_values.get_fe().dofs_per_cell; - const unsigned int n_q_points = scratch.stokes_fe_values.n_quadrature_points; + template + void + BoussinesqFlowProblem:: + local_assemble_stokes_system (const typename DoFHandler::active_cell_iterator &cell, + Assembly::Scratch::StokesSystem &scratch, + Assembly::CopyData::StokesSystem &data) + { + const unsigned int dofs_per_cell = scratch.stokes_fe_values.get_fe().dofs_per_cell; + const unsigned int n_q_points = scratch.stokes_fe_values.n_quadrature_points; - const FEValuesExtractors::Vector velocities (0); - const FEValuesExtractors::Scalar pressure (dim); + const FEValuesExtractors::Vector velocities (0); + const FEValuesExtractors::Scalar pressure (dim); - scratch.stokes_fe_values.reinit (cell); + scratch.stokes_fe_values.reinit (cell); - typename DoFHandler::active_cell_iterator - temperature_cell (&triangulation, - cell->level(), - cell->index(), - &temperature_dof_handler); - scratch.temperature_fe_values.reinit (temperature_cell); + typename DoFHandler::active_cell_iterator + temperature_cell (&triangulation, + cell->level(), + cell->index(), + &temperature_dof_handler); + scratch.temperature_fe_values.reinit (temperature_cell); - if (rebuild_stokes_matrix) - data.local_matrix = 0; - data.local_rhs = 0; + if (rebuild_stokes_matrix) + data.local_matrix = 0; + data.local_rhs = 0; - scratch.temperature_fe_values.get_function_values (old_temperature_solution, - scratch.old_temperature_values); + scratch.temperature_fe_values.get_function_values (old_temperature_solution, + scratch.old_temperature_values); - for (unsigned int q=0; q - gravity = EquationData::gravity_vector (scratch.stokes_fe_values - .quadrature_point(q)); - - for (unsigned int i=0; i + gravity = EquationData::gravity_vector (scratch.stokes_fe_values + .quadrature_point(q)); + + for (unsigned int i=0; iget_dof_indices (data.local_dof_indices); -} + cell->get_dof_indices (data.local_dof_indices); + } -template -void -BoussinesqFlowProblem:: -copy_local_to_global_stokes_system (const Assembly::CopyData::StokesSystem &data) -{ - if (rebuild_stokes_matrix == true) - stokes_constraints.distribute_local_to_global (data.local_matrix, - data.local_rhs, - data.local_dof_indices, - stokes_matrix, - stokes_rhs); - else - stokes_constraints.distribute_local_to_global (data.local_rhs, - data.local_dof_indices, - stokes_rhs); -} + template + void + BoussinesqFlowProblem:: + copy_local_to_global_stokes_system (const Assembly::CopyData::StokesSystem &data) + { + if (rebuild_stokes_matrix == true) + stokes_constraints.distribute_local_to_global (data.local_matrix, + data.local_rhs, + data.local_dof_indices, + stokes_matrix, + stokes_rhs); + else + stokes_constraints.distribute_local_to_global (data.local_rhs, + data.local_dof_indices, + stokes_rhs); + } -template -void BoussinesqFlowProblem::assemble_stokes_system () -{ - computing_timer.enter_section (" Assemble Stokes system"); - - if (rebuild_stokes_matrix == true) - stokes_matrix=0; - - stokes_rhs=0; - - const QGauss quadrature_formula(parameters.stokes_velocity_degree+1); - - typedef - FilteredIterator::active_cell_iterator> - SubdomainFilter; - - WorkStream:: - run (SubdomainFilter (IteratorFilters::SubdomainEqualTo - (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)), - stokes_dof_handler.begin_active()), - SubdomainFilter (IteratorFilters::SubdomainEqualTo - (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)), - stokes_dof_handler.end()), - std_cxx1x::bind (&BoussinesqFlowProblem:: - local_assemble_stokes_system, - this, - std_cxx1x::_1, - std_cxx1x::_2, - std_cxx1x::_3), - std_cxx1x::bind (&BoussinesqFlowProblem:: - copy_local_to_global_stokes_system, - this, - std_cxx1x::_1), - Assembly::Scratch:: - StokesSystem (stokes_fe, mapping, quadrature_formula, - (update_values | - update_quadrature_points | - update_JxW_values | - (rebuild_stokes_matrix == true - ? - update_gradients - : - UpdateFlags(0))), - temperature_fe, - update_values), - Assembly::CopyData:: - StokesSystem (stokes_fe)); - - stokes_matrix.compress(); - stokes_rhs.compress(Add); - - rebuild_stokes_matrix = false; - - pcout << std::endl; - computing_timer.exit_section(); -} + template + void BoussinesqFlowProblem::assemble_stokes_system () + { + computing_timer.enter_section (" Assemble Stokes system"); + + if (rebuild_stokes_matrix == true) + stokes_matrix=0; + + stokes_rhs=0; + + const QGauss quadrature_formula(parameters.stokes_velocity_degree+1); + + typedef + FilteredIterator::active_cell_iterator> + SubdomainFilter; + + WorkStream:: + run (SubdomainFilter (IteratorFilters::SubdomainEqualTo + (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)), + stokes_dof_handler.begin_active()), + SubdomainFilter (IteratorFilters::SubdomainEqualTo + (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)), + stokes_dof_handler.end()), + std_cxx1x::bind (&BoussinesqFlowProblem:: + local_assemble_stokes_system, + this, + std_cxx1x::_1, + std_cxx1x::_2, + std_cxx1x::_3), + std_cxx1x::bind (&BoussinesqFlowProblem:: + copy_local_to_global_stokes_system, + this, + std_cxx1x::_1), + Assembly::Scratch:: + StokesSystem (stokes_fe, mapping, quadrature_formula, + (update_values | + update_quadrature_points | + update_JxW_values | + (rebuild_stokes_matrix == true + ? + update_gradients + : + UpdateFlags(0))), + temperature_fe, + update_values), + Assembly::CopyData:: + StokesSystem (stokes_fe)); + + stokes_matrix.compress(); + stokes_rhs.compress(Add); + + rebuild_stokes_matrix = false; + + pcout << std::endl; + computing_timer.exit_section(); + } // @sect5{Temperature matrix assembly} @@ -2705,106 +2707,106 @@ void BoussinesqFlowProblem::assemble_stokes_system () // // The two following functions perform // similar services as the ones above. -template -void BoussinesqFlowProblem:: -local_assemble_temperature_matrix (const typename DoFHandler::active_cell_iterator &cell, - Assembly::Scratch::TemperatureMatrix &scratch, - Assembly::CopyData::TemperatureMatrix &data) -{ - const unsigned int dofs_per_cell = scratch.temperature_fe_values.get_fe().dofs_per_cell; - const unsigned int n_q_points = scratch.temperature_fe_values.n_quadrature_points; + template + void BoussinesqFlowProblem:: + local_assemble_temperature_matrix (const typename DoFHandler::active_cell_iterator &cell, + Assembly::Scratch::TemperatureMatrix &scratch, + Assembly::CopyData::TemperatureMatrix &data) + { + const unsigned int dofs_per_cell = scratch.temperature_fe_values.get_fe().dofs_per_cell; + const unsigned int n_q_points = scratch.temperature_fe_values.n_quadrature_points; - scratch.temperature_fe_values.reinit (cell); - cell->get_dof_indices (data.local_dof_indices); + scratch.temperature_fe_values.reinit (cell); + cell->get_dof_indices (data.local_dof_indices); - data.local_mass_matrix = 0; - data.local_stiffness_matrix = 0; + data.local_mass_matrix = 0; + data.local_stiffness_matrix = 0; - for (unsigned int q=0; q -void -BoussinesqFlowProblem:: -copy_local_to_global_temperature_matrix (const Assembly::CopyData::TemperatureMatrix &data) -{ - temperature_constraints.distribute_local_to_global (data.local_mass_matrix, - data.local_dof_indices, - temperature_mass_matrix); - temperature_constraints.distribute_local_to_global (data.local_stiffness_matrix, - data.local_dof_indices, - temperature_stiffness_matrix); -} + template + void + BoussinesqFlowProblem:: + copy_local_to_global_temperature_matrix (const Assembly::CopyData::TemperatureMatrix &data) + { + temperature_constraints.distribute_local_to_global (data.local_mass_matrix, + data.local_dof_indices, + temperature_mass_matrix); + temperature_constraints.distribute_local_to_global (data.local_stiffness_matrix, + data.local_dof_indices, + temperature_stiffness_matrix); + } -template -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; - - const QGauss quadrature_formula(parameters.temperature_degree+2); - - typedef - FilteredIterator::active_cell_iterator> - SubdomainFilter; - - WorkStream:: - run (SubdomainFilter (IteratorFilters::SubdomainEqualTo - (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)), - temperature_dof_handler.begin_active()), - SubdomainFilter (IteratorFilters::SubdomainEqualTo - (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)), - temperature_dof_handler.end()), - std_cxx1x::bind (&BoussinesqFlowProblem:: - local_assemble_temperature_matrix, - this, - std_cxx1x::_1, - std_cxx1x::_2, - std_cxx1x::_3), - std_cxx1x::bind (&BoussinesqFlowProblem:: - copy_local_to_global_temperature_matrix, - this, - std_cxx1x::_1), - Assembly::Scratch:: - TemperatureMatrix (temperature_fe, mapping, quadrature_formula), - Assembly::CopyData:: - TemperatureMatrix (temperature_fe)); - - temperature_mass_matrix.compress(); - temperature_stiffness_matrix.compress(); - - rebuild_temperature_matrices = false; - rebuild_temperature_preconditioner = true; - - computing_timer.exit_section(); -} + template + 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; + + const QGauss quadrature_formula(parameters.temperature_degree+2); + + typedef + FilteredIterator::active_cell_iterator> + SubdomainFilter; + + WorkStream:: + run (SubdomainFilter (IteratorFilters::SubdomainEqualTo + (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)), + temperature_dof_handler.begin_active()), + SubdomainFilter (IteratorFilters::SubdomainEqualTo + (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)), + temperature_dof_handler.end()), + std_cxx1x::bind (&BoussinesqFlowProblem:: + local_assemble_temperature_matrix, + this, + std_cxx1x::_1, + std_cxx1x::_2, + std_cxx1x::_3), + std_cxx1x::bind (&BoussinesqFlowProblem:: + copy_local_to_global_temperature_matrix, + this, + std_cxx1x::_1), + Assembly::Scratch:: + TemperatureMatrix (temperature_fe, mapping, quadrature_formula), + Assembly::CopyData:: + TemperatureMatrix (temperature_fe)); + + temperature_mass_matrix.compress(); + temperature_stiffness_matrix.compress(); + + rebuild_temperature_matrices = false; + rebuild_temperature_preconditioner = true; + + computing_timer.exit_section(); + } // @sect5{Temperature right hand side assembly} @@ -2831,170 +2833,170 @@ void BoussinesqFlowProblem::assemble_temperature_matrix () // dofs. That will account for the correct // balance of the right hand side vector // with the matrix system of temperature. -template -void BoussinesqFlowProblem:: -local_assemble_temperature_rhs (const std::pair global_T_range, - const double global_max_velocity, - const double global_entropy_variation, - const typename DoFHandler::active_cell_iterator &cell, - Assembly::Scratch::TemperatureRHS &scratch, - Assembly::CopyData::TemperatureRHS &data) -{ - const bool use_bdf2_scheme = (timestep_number != 0); - - const unsigned int dofs_per_cell = scratch.temperature_fe_values.get_fe().dofs_per_cell; - const unsigned int n_q_points = scratch.temperature_fe_values.n_quadrature_points; - - const FEValuesExtractors::Vector velocities (0); - - data.local_rhs = 0; - data.matrix_for_bc = 0; - cell->get_dof_indices (data.local_dof_indices); - - scratch.temperature_fe_values.reinit (cell); - - typename DoFHandler::active_cell_iterator - stokes_cell (&triangulation, - cell->level(), - cell->index(), - &stokes_dof_handler); - scratch.stokes_fe_values.reinit (stokes_cell); - - scratch.temperature_fe_values.get_function_values (old_temperature_solution, - scratch.old_temperature_values); - scratch.temperature_fe_values.get_function_values (old_old_temperature_solution, - scratch.old_old_temperature_values); - - scratch.temperature_fe_values.get_function_gradients (old_temperature_solution, - scratch.old_temperature_grads); - scratch.temperature_fe_values.get_function_gradients (old_old_temperature_solution, - scratch.old_old_temperature_grads); - - scratch.temperature_fe_values.get_function_laplacians (old_temperature_solution, - scratch.old_temperature_laplacians); - scratch.temperature_fe_values.get_function_laplacians (old_old_temperature_solution, - scratch.old_old_temperature_laplacians); - - scratch.stokes_fe_values[velocities].get_function_values (stokes_solution, - scratch.old_velocity_values); - scratch.stokes_fe_values[velocities].get_function_values (old_stokes_solution, - scratch.old_old_velocity_values); - scratch.stokes_fe_values[velocities].get_function_symmetric_gradients (stokes_solution, - scratch.old_strain_rates); - scratch.stokes_fe_values[velocities].get_function_symmetric_gradients (old_stokes_solution, - scratch.old_old_strain_rates); - - const double nu - = compute_viscosity (scratch.old_temperature_values, - scratch.old_old_temperature_values, - scratch.old_temperature_grads, - scratch.old_old_temperature_grads, - scratch.old_temperature_laplacians, - scratch.old_old_temperature_laplacians, - scratch.old_velocity_values, - scratch.old_old_velocity_values, - scratch.old_strain_rates, - scratch.old_old_strain_rates, - global_max_velocity, - global_T_range.second - global_T_range.first, - 0.5 * (global_T_range.second + global_T_range.first), - global_entropy_variation, - cell->diameter()); - - for (unsigned int q=0; q ext_grad_T - = (use_bdf2_scheme ? - (scratch.old_temperature_grads[q] * - (1+time_step/old_time_step) - - - scratch.old_old_temperature_grads[q] * - time_step / old_time_step) - : - scratch.old_temperature_grads[q]); - - const Tensor<1,dim> extrapolated_u - = (use_bdf2_scheme ? - (scratch.old_velocity_values[q] * (1+time_step/old_time_step) - - scratch.old_old_velocity_values[q] * time_step/old_time_step) - : - scratch.old_velocity_values[q]); - const SymmetricTensor<2,dim> extrapolated_strain_rate - = (use_bdf2_scheme ? - (scratch.old_strain_rates[q] * (1+time_step/old_time_step) - - scratch.old_old_strain_rates[q] * time_step/old_time_step) - : - scratch.old_strain_rates[q]); - - const double gamma - = ((EquationData::radiogenic_heating * EquationData::density(old_Ts) //?????? why old_Ts? - + - 2 * EquationData::eta * extrapolated_strain_rate * extrapolated_strain_rate) / - (EquationData::density(old_Ts) * EquationData::specific_heat)); - - for (unsigned int i=0; i + void BoussinesqFlowProblem:: + local_assemble_temperature_rhs (const std::pair global_T_range, + const double global_max_velocity, + const double global_entropy_variation, + const typename DoFHandler::active_cell_iterator &cell, + Assembly::Scratch::TemperatureRHS &scratch, + Assembly::CopyData::TemperatureRHS &data) + { + const bool use_bdf2_scheme = (timestep_number != 0); + const unsigned int dofs_per_cell = scratch.temperature_fe_values.get_fe().dofs_per_cell; + const unsigned int n_q_points = scratch.temperature_fe_values.n_quadrature_points; -template -void -BoussinesqFlowProblem:: -copy_local_to_global_temperature_rhs (const Assembly::CopyData::TemperatureRHS &data) -{ - temperature_constraints.distribute_local_to_global (data.local_rhs, - data.local_dof_indices, - temperature_rhs, - data.matrix_for_bc); -} + const FEValuesExtractors::Vector velocities (0); + + data.local_rhs = 0; + data.matrix_for_bc = 0; + cell->get_dof_indices (data.local_dof_indices); + + scratch.temperature_fe_values.reinit (cell); + + typename DoFHandler::active_cell_iterator + stokes_cell (&triangulation, + cell->level(), + cell->index(), + &stokes_dof_handler); + scratch.stokes_fe_values.reinit (stokes_cell); + + scratch.temperature_fe_values.get_function_values (old_temperature_solution, + scratch.old_temperature_values); + scratch.temperature_fe_values.get_function_values (old_old_temperature_solution, + scratch.old_old_temperature_values); + + scratch.temperature_fe_values.get_function_gradients (old_temperature_solution, + scratch.old_temperature_grads); + scratch.temperature_fe_values.get_function_gradients (old_old_temperature_solution, + scratch.old_old_temperature_grads); + + scratch.temperature_fe_values.get_function_laplacians (old_temperature_solution, + scratch.old_temperature_laplacians); + scratch.temperature_fe_values.get_function_laplacians (old_old_temperature_solution, + scratch.old_old_temperature_laplacians); + + scratch.stokes_fe_values[velocities].get_function_values (stokes_solution, + scratch.old_velocity_values); + scratch.stokes_fe_values[velocities].get_function_values (old_stokes_solution, + scratch.old_old_velocity_values); + scratch.stokes_fe_values[velocities].get_function_symmetric_gradients (stokes_solution, + scratch.old_strain_rates); + scratch.stokes_fe_values[velocities].get_function_symmetric_gradients (old_stokes_solution, + scratch.old_old_strain_rates); + + const double nu + = compute_viscosity (scratch.old_temperature_values, + scratch.old_old_temperature_values, + scratch.old_temperature_grads, + scratch.old_old_temperature_grads, + scratch.old_temperature_laplacians, + scratch.old_old_temperature_laplacians, + scratch.old_velocity_values, + scratch.old_old_velocity_values, + scratch.old_strain_rates, + scratch.old_old_strain_rates, + global_max_velocity, + global_T_range.second - global_T_range.first, + 0.5 * (global_T_range.second + global_T_range.first), + global_entropy_variation, + cell->diameter()); + + for (unsigned int q=0; q ext_grad_T + = (use_bdf2_scheme ? + (scratch.old_temperature_grads[q] * + (1+time_step/old_time_step) + - + scratch.old_old_temperature_grads[q] * + time_step / old_time_step) + : + scratch.old_temperature_grads[q]); + + const Tensor<1,dim> extrapolated_u + = (use_bdf2_scheme ? + (scratch.old_velocity_values[q] * (1+time_step/old_time_step) - + scratch.old_old_velocity_values[q] * time_step/old_time_step) + : + scratch.old_velocity_values[q]); + const SymmetricTensor<2,dim> extrapolated_strain_rate + = (use_bdf2_scheme ? + (scratch.old_strain_rates[q] * (1+time_step/old_time_step) - + scratch.old_old_strain_rates[q] * time_step/old_time_step) + : + scratch.old_strain_rates[q]); + + const double gamma + = ((EquationData::radiogenic_heating * EquationData::density(old_Ts) //?????? why old_Ts? + + + 2 * EquationData::eta * extrapolated_strain_rate * extrapolated_strain_rate) / + (EquationData::density(old_Ts) * EquationData::specific_heat)); + + for (unsigned int i=0; i + void + BoussinesqFlowProblem:: + copy_local_to_global_temperature_rhs (const Assembly::CopyData::TemperatureRHS &data) + { + temperature_constraints.distribute_local_to_global (data.local_rhs, + data.local_dof_indices, + temperature_rhs, + data.matrix_for_bc); + } @@ -3024,82 +3026,82 @@ copy_local_to_global_temperature_rhs (const Assembly::CopyData::TemperatureRHS -void BoussinesqFlowProblem::assemble_temperature_system (const double maximal_velocity) -{ - const bool use_bdf2_scheme = (timestep_number != 0); + template + void BoussinesqFlowProblem::assemble_temperature_system (const double maximal_velocity) + { + const bool use_bdf2_scheme = (timestep_number != 0); - if (use_bdf2_scheme == true) - { - temperature_matrix.copy_from (temperature_mass_matrix); - temperature_matrix *= (2*time_step + old_time_step) / - (time_step + old_time_step); - temperature_matrix.add (time_step, temperature_stiffness_matrix); - } - else - { - temperature_matrix.copy_from (temperature_mass_matrix); - temperature_matrix.add (time_step, temperature_stiffness_matrix); - } - temperature_matrix.compress(); + if (use_bdf2_scheme == true) + { + temperature_matrix.copy_from (temperature_mass_matrix); + temperature_matrix *= (2*time_step + old_time_step) / + (time_step + old_time_step); + temperature_matrix.add (time_step, temperature_stiffness_matrix); + } + else + { + temperature_matrix.copy_from (temperature_mass_matrix); + temperature_matrix.add (time_step, temperature_stiffness_matrix); + } + temperature_matrix.compress(); - if (rebuild_temperature_preconditioner == true) - { - T_preconditioner = std_cxx1x::shared_ptr - (new TrilinosWrappers::PreconditionIC()); - T_preconditioner->initialize (temperature_matrix); - rebuild_temperature_preconditioner = false; - } + if (rebuild_temperature_preconditioner == true) + { + T_preconditioner = std_cxx1x::shared_ptr + (new TrilinosWrappers::PreconditionIC()); + T_preconditioner->initialize (temperature_matrix); + rebuild_temperature_preconditioner = false; + } - temperature_rhs = 0; - - const QGauss quadrature_formula(parameters.temperature_degree+2); - const std::pair - global_T_range = get_extrapolated_temperature_range(); - - // use midpoint between maximum and minimum - // temperature for definition of average - // temperature in entropy viscosity. Could - // also use the integral average, but the - // results are not very sensitive to this - // choice. - const double average_temperature = 0.5 * (global_T_range.first + - global_T_range.second); - const double global_entropy_variation = - get_entropy_variation (average_temperature); - - typedef - FilteredIterator::active_cell_iterator> - SubdomainFilter; - - WorkStream:: - run (SubdomainFilter (IteratorFilters::SubdomainEqualTo - (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)), - temperature_dof_handler.begin_active()), - SubdomainFilter (IteratorFilters::SubdomainEqualTo - (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)), - temperature_dof_handler.end()), - std_cxx1x::bind (&BoussinesqFlowProblem:: - local_assemble_temperature_rhs, - this, - global_T_range, - maximal_velocity, - global_entropy_variation, - std_cxx1x::_1, - std_cxx1x::_2, - std_cxx1x::_3), - std_cxx1x::bind (&BoussinesqFlowProblem:: - copy_local_to_global_temperature_rhs, - this, - std_cxx1x::_1), - Assembly::Scratch:: - TemperatureRHS (temperature_fe, stokes_fe, mapping, - quadrature_formula), - Assembly::CopyData:: - TemperatureRHS (temperature_fe)); - - temperature_rhs.compress(Add); -} + temperature_rhs = 0; + + const QGauss quadrature_formula(parameters.temperature_degree+2); + const std::pair + global_T_range = get_extrapolated_temperature_range(); + + // use midpoint between maximum and minimum + // temperature for definition of average + // temperature in entropy viscosity. Could + // also use the integral average, but the + // results are not very sensitive to this + // choice. + const double average_temperature = 0.5 * (global_T_range.first + + global_T_range.second); + const double global_entropy_variation = + get_entropy_variation (average_temperature); + + typedef + FilteredIterator::active_cell_iterator> + SubdomainFilter; + + WorkStream:: + run (SubdomainFilter (IteratorFilters::SubdomainEqualTo + (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)), + temperature_dof_handler.begin_active()), + SubdomainFilter (IteratorFilters::SubdomainEqualTo + (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)), + temperature_dof_handler.end()), + std_cxx1x::bind (&BoussinesqFlowProblem:: + local_assemble_temperature_rhs, + this, + global_T_range, + maximal_velocity, + global_entropy_variation, + std_cxx1x::_1, + std_cxx1x::_2, + std_cxx1x::_3), + std_cxx1x::bind (&BoussinesqFlowProblem:: + copy_local_to_global_temperature_rhs, + this, + std_cxx1x::_1), + Assembly::Scratch:: + TemperatureRHS (temperature_fe, stokes_fe, mapping, + quadrature_formula), + Assembly::CopyData:: + TemperatureRHS (temperature_fe)); + + temperature_rhs.compress(Add); + } @@ -3150,333 +3152,333 @@ void BoussinesqFlowProblem::assemble_temperature_system (const double maxim // Apart from these two changes, everything // is the same as in step-31, so we don't // need to further comment on it. -template -void BoussinesqFlowProblem::solve () -{ - computing_timer.enter_section (" Solve Stokes system"); - + template + void BoussinesqFlowProblem::solve () { - pcout << " Solving Stokes system... " << std::flush; - - TrilinosWrappers::MPI::BlockVector - distributed_stokes_solution (stokes_rhs); -// distributed_stokes_solution = stokes_solution; - distributed_stokes_solution.block(0).reinit(stokes_solution.block(0),false,true); - distributed_stokes_solution.block(1).reinit(stokes_solution.block(1),false,true); - - - const unsigned int - start = (distributed_stokes_solution.block(0).size() + - distributed_stokes_solution.block(1).local_range().first), - end = (distributed_stokes_solution.block(0).size() + - distributed_stokes_solution.block(1).local_range().second); - for (unsigned int i=start; i mem; - - // step 1: try if the simple and fast solver - // succeeds in 30 steps or less. - unsigned int n_iterations = 0; - double reduction = 0; - const double solver_tolerance = 1e-7 * stokes_rhs.l2_norm(); - SolverControl solver_control (30, solver_tolerance); - - try - { - const LinearSolvers::RightPrecond - preconditioner (stokes_matrix, stokes_preconditioner_matrix, - *Mp_preconditioner, *Amg_preconditioner, - false); - - SolverFGMRES - solver(solver_control, mem, - SolverFGMRES:: - AdditionalData(30, true)); - solver.solve(stokes_matrix, distributed_stokes_solution, stokes_rhs, - preconditioner); - - n_iterations = solver_control.last_step(); - reduction = solver_control.last_value()/solver_control.initial_value(); - } - - // step 2: take the stronger solver in case - // the simple solver failed - catch (SolverControl::NoConvergence) - { - const LinearSolvers::RightPrecond - preconditioner (stokes_matrix, stokes_preconditioner_matrix, - *Mp_preconditioner, *Amg_preconditioner, - true); - - SolverControl solver_control_refined (stokes_matrix.m(), solver_tolerance); - SolverFGMRES - solver(solver_control_refined, mem, - SolverFGMRES:: - AdditionalData(50, true)); - solver.solve(stokes_matrix, distributed_stokes_solution, stokes_rhs, - preconditioner); - - n_iterations = (solver_control.last_step() + - solver_control_refined.last_step()); - reduction = (solver_control_refined.last_value()/ - std::max(solver_control.initial_value(), - solver_control_refined.initial_value())); - } + computing_timer.enter_section (" Solve Stokes system"); + { + pcout << " Solving Stokes system... " << std::flush; - stokes_constraints.distribute (distributed_stokes_solution); - stokes_solution.block(0).reinit(distributed_stokes_solution.block(0), - false, true); - stokes_solution.block(1).reinit(distributed_stokes_solution.block(1), - false, true); - - pcout << n_iterations << " iterations." - << " Reduced residual by " << reduction - << std::endl; - - TrilinosWrappers::MPI::Vector tmp; - tmp.reinit (stokes_rhs.block(1)); - pcout << " Relative divergence residual: " - << stokes_matrix.block(1,0).residual (tmp, - distributed_stokes_solution.block(0), - stokes_rhs.block(1)) - / - distributed_stokes_solution.block(0).l2_norm() / EquationData::pressure_scaling - << std::endl; - - pcout << " Relative vector sizes: " - << distributed_stokes_solution.block(0).linfty_norm() << ' ' - << distributed_stokes_solution.block(1).linfty_norm() << std::endl; - } - computing_timer.exit_section(); + TrilinosWrappers::MPI::BlockVector + distributed_stokes_solution (stokes_rhs); +// distributed_stokes_solution = stokes_solution; + distributed_stokes_solution.block(0).reinit(stokes_solution.block(0),false,true); + distributed_stokes_solution.block(1).reinit(stokes_solution.block(1),false,true); + + + const unsigned int + start = (distributed_stokes_solution.block(0).size() + + distributed_stokes_solution.block(1).local_range().first), + end = (distributed_stokes_solution.block(0).size() + + distributed_stokes_solution.block(1).local_range().second); + for (unsigned int i=start; i mem; + + // step 1: try if the simple and fast solver + // succeeds in 30 steps or less. + unsigned int n_iterations = 0; + double reduction = 0; + const double solver_tolerance = 1e-7 * stokes_rhs.l2_norm(); + SolverControl solver_control (30, solver_tolerance); + + try + { + const LinearSolvers::RightPrecond + preconditioner (stokes_matrix, stokes_preconditioner_matrix, + *Mp_preconditioner, *Amg_preconditioner, + false); + + SolverFGMRES + solver(solver_control, mem, + SolverFGMRES:: + AdditionalData(30, true)); + solver.solve(stokes_matrix, distributed_stokes_solution, stokes_rhs, + preconditioner); + + n_iterations = solver_control.last_step(); + reduction = solver_control.last_value()/solver_control.initial_value(); + } + + // step 2: take the stronger solver in case + // the simple solver failed + catch (SolverControl::NoConvergence) + { + const LinearSolvers::RightPrecond + preconditioner (stokes_matrix, stokes_preconditioner_matrix, + *Mp_preconditioner, *Amg_preconditioner, + true); + + SolverControl solver_control_refined (stokes_matrix.m(), solver_tolerance); + SolverFGMRES + solver(solver_control_refined, mem, + SolverFGMRES:: + AdditionalData(50, true)); + solver.solve(stokes_matrix, distributed_stokes_solution, stokes_rhs, + preconditioner); + + n_iterations = (solver_control.last_step() + + solver_control_refined.last_step()); + reduction = (solver_control_refined.last_value()/ + std::max(solver_control.initial_value(), + solver_control_refined.initial_value())); + } + + + stokes_constraints.distribute (distributed_stokes_solution); + stokes_solution.block(0).reinit(distributed_stokes_solution.block(0), + false, true); + stokes_solution.block(1).reinit(distributed_stokes_solution.block(1), + false, true); + + pcout << n_iterations << " iterations." + << " Reduced residual by " << reduction + << std::endl; + + TrilinosWrappers::MPI::Vector tmp; + tmp.reinit (stokes_rhs.block(1)); + pcout << " Relative divergence residual: " + << stokes_matrix.block(1,0).residual (tmp, + distributed_stokes_solution.block(0), + stokes_rhs.block(1)) + / + distributed_stokes_solution.block(0).l2_norm() / EquationData::pressure_scaling + << std::endl; + + pcout << " Relative vector sizes: " + << distributed_stokes_solution.block(0).linfty_norm() << ' ' + << distributed_stokes_solution.block(1).linfty_norm() << std::endl; + } + computing_timer.exit_section(); - computing_timer.enter_section (" Assemble temperature rhs"); - { - old_time_step = time_step; - const double cfl_number = get_cfl_number(); - - // we found out that we need - // approximately a quarter the time step - // size in 3d - double scaling = (dim==3)?0.25:1.0; - time_step = (scaling/(2.1*dim*std::sqrt(1.*dim)) / - (parameters.temperature_degree * - cfl_number)); - - const double maximal_velocity = get_maximal_velocity(); - pcout << " Maximal velocity: " - << maximal_velocity *EquationData::year_in_seconds * 100 - << " cm/year" - << std::endl; - pcout << " " << "Time step: " - << time_step/EquationData::year_in_seconds - << " years" - << std::endl; - - temperature_solution = old_temperature_solution; - assemble_temperature_system (maximal_velocity); - } - computing_timer.exit_section (); + computing_timer.enter_section (" Assemble temperature rhs"); + { + old_time_step = time_step; + const double cfl_number = get_cfl_number(); + + // we found out that we need + // approximately a quarter the time step + // size in 3d + double scaling = (dim==3)?0.25:1.0; + time_step = (scaling/(2.1*dim*std::sqrt(1.*dim)) / + (parameters.temperature_degree * + cfl_number)); + + const double maximal_velocity = get_maximal_velocity(); + pcout << " Maximal velocity: " + << maximal_velocity *EquationData::year_in_seconds * 100 + << " cm/year" + << std::endl; + pcout << " " << "Time step: " + << time_step/EquationData::year_in_seconds + << " years" + << std::endl; + + temperature_solution = old_temperature_solution; + assemble_temperature_system (maximal_velocity); + } + computing_timer.exit_section (); - computing_timer.enter_section (" Solve temperature system"); - { - SolverControl solver_control (temperature_matrix.m(), - 1e-12*temperature_rhs.l2_norm()); - SolverCG cg (solver_control); + computing_timer.enter_section (" Solve temperature system"); + { + SolverControl solver_control (temperature_matrix.m(), + 1e-12*temperature_rhs.l2_norm()); + SolverCG cg (solver_control); - TrilinosWrappers::MPI::Vector - distributed_temperature_solution (temperature_rhs); + TrilinosWrappers::MPI::Vector + distributed_temperature_solution (temperature_rhs); // distributed_temperature_solution = temperature_solution; - distributed_temperature_solution.reinit(temperature_solution, false, true); + distributed_temperature_solution.reinit(temperature_solution, false, true); - cg.solve (temperature_matrix, distributed_temperature_solution, - temperature_rhs, *T_preconditioner); + cg.solve (temperature_matrix, distributed_temperature_solution, + temperature_rhs, *T_preconditioner); - temperature_constraints.distribute (distributed_temperature_solution); + temperature_constraints.distribute (distributed_temperature_solution); // temperature_solution = distributed_temperature_solution; - temperature_solution.reinit(distributed_temperature_solution, false, true); - - pcout << " " - << solver_control.last_step() - << " CG iterations for temperature" << std::endl; - computing_timer.exit_section(); - - // extract temperature range - std::vector temperature (2), global_temperature (2); - temperature[0] = std::numeric_limits::max(), - temperature[1] = -std::numeric_limits::max(); - for (unsigned int i=0; i (temperature[0], - distributed_temperature_solution.trilinos_vector()[0][i]); - temperature[1] = std::max (temperature[1], - distributed_temperature_solution.trilinos_vector()[0][i]); - } + temperature_solution.reinit(distributed_temperature_solution, false, true); + + pcout << " " + << solver_control.last_step() + << " CG iterations for temperature" << std::endl; + computing_timer.exit_section(); + + // extract temperature range + std::vector temperature (2), global_temperature (2); + temperature[0] = std::numeric_limits::max(), + temperature[1] = -std::numeric_limits::max(); + for (unsigned int i=0; i (temperature[0], + distributed_temperature_solution.trilinos_vector()[0][i]); + temperature[1] = std::max (temperature[1], + distributed_temperature_solution.trilinos_vector()[0][i]); + } #ifdef DEAL_II_COMPILER_SUPPORTS_MPI - temperature[0] *= -1.0; - MPI_Allreduce (&temperature[0], &global_temperature[0], - 2, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); - global_temperature[0] *= -1.0; + temperature[0] *= -1.0; + MPI_Allreduce (&temperature[0], &global_temperature[0], + 2, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); + global_temperature[0] *= -1.0; #else - global_temperature = local_temperature; + global_temperature = local_temperature; #endif - pcout << " Temperature range: " - << global_temperature[0] << ' ' << global_temperature[1] - << std::endl; + pcout << " Temperature range: " + << global_temperature[0] << ' ' << global_temperature[1] + << std::endl; + } } -} // @sect4{BoussinesqFlowProblem::output_results} -template -class BoussinesqFlowProblem::Postprocessor : public DataPostprocessor -{ - public: - Postprocessor (const unsigned int partition, - const double minimal_pressure); + template + class BoussinesqFlowProblem::Postprocessor : public DataPostprocessor + { + public: + Postprocessor (const unsigned int partition, + const double minimal_pressure); - virtual - void - compute_derived_quantities_vector (const std::vector > &uh, - const std::vector > > &duh, - const std::vector > > &dduh, - const std::vector > &normals, - const std::vector > &evaluation_points, - std::vector > &computed_quantities) const; + virtual + void + compute_derived_quantities_vector (const std::vector > &uh, + const std::vector > > &duh, + const std::vector > > &dduh, + const std::vector > &normals, + const std::vector > &evaluation_points, + std::vector > &computed_quantities) const; - virtual std::vector get_names () const; + virtual std::vector get_names () const; - virtual unsigned int n_output_variables() const; + virtual unsigned int n_output_variables() const; - virtual - std::vector - get_data_component_interpretation () const; + virtual + std::vector + get_data_component_interpretation () const; - virtual UpdateFlags get_needed_update_flags () const; + virtual UpdateFlags get_needed_update_flags () const; - private: - const unsigned int partition; - const double minimal_pressure; -}; + private: + const unsigned int partition; + const double minimal_pressure; + }; -template -BoussinesqFlowProblem::Postprocessor:: -Postprocessor (const unsigned int partition, - const double minimal_pressure) - : - partition (partition), - minimal_pressure (minimal_pressure) -{} + template + BoussinesqFlowProblem::Postprocessor:: + Postprocessor (const unsigned int partition, + const double minimal_pressure) + : + partition (partition), + minimal_pressure (minimal_pressure) + {} -template -std::vector -BoussinesqFlowProblem::Postprocessor::get_names() const -{ - std::vector solution_names (dim, "velocity"); - solution_names.push_back ("p"); - solution_names.push_back ("T"); - solution_names.push_back ("friction_heating"); - solution_names.push_back ("partition"); - solution_names.push_back ("non_adiabatic_pressure"); - solution_names.push_back ("non_adiabatic_temperature"); - return solution_names; -} + template + std::vector + BoussinesqFlowProblem::Postprocessor::get_names() const + { + std::vector solution_names (dim, "velocity"); + solution_names.push_back ("p"); + solution_names.push_back ("T"); + solution_names.push_back ("friction_heating"); + solution_names.push_back ("partition"); + solution_names.push_back ("non_adiabatic_pressure"); + solution_names.push_back ("non_adiabatic_temperature"); + return solution_names; + } -template -unsigned int -BoussinesqFlowProblem::Postprocessor::n_output_variables() const -{ - return get_names().size(); -} + template + unsigned int + BoussinesqFlowProblem::Postprocessor::n_output_variables() const + { + return get_names().size(); + } -template -std::vector -BoussinesqFlowProblem::Postprocessor:: -get_data_component_interpretation () const -{ + template std::vector - interpretation (dim, - DataComponentInterpretation::component_is_part_of_vector); + BoussinesqFlowProblem::Postprocessor:: + get_data_component_interpretation () const + { + std::vector + interpretation (dim, + DataComponentInterpretation::component_is_part_of_vector); - interpretation.push_back (DataComponentInterpretation::component_is_scalar); - interpretation.push_back (DataComponentInterpretation::component_is_scalar); - interpretation.push_back (DataComponentInterpretation::component_is_scalar); - interpretation.push_back (DataComponentInterpretation::component_is_scalar); - interpretation.push_back (DataComponentInterpretation::component_is_scalar); - interpretation.push_back (DataComponentInterpretation::component_is_scalar); + interpretation.push_back (DataComponentInterpretation::component_is_scalar); + interpretation.push_back (DataComponentInterpretation::component_is_scalar); + interpretation.push_back (DataComponentInterpretation::component_is_scalar); + interpretation.push_back (DataComponentInterpretation::component_is_scalar); + interpretation.push_back (DataComponentInterpretation::component_is_scalar); + interpretation.push_back (DataComponentInterpretation::component_is_scalar); - return interpretation; -} + return interpretation; + } -template -UpdateFlags -BoussinesqFlowProblem::Postprocessor::get_needed_update_flags() const -{ - return update_values | update_gradients | update_q_points; -} + template + UpdateFlags + BoussinesqFlowProblem::Postprocessor::get_needed_update_flags() const + { + return update_values | update_gradients | update_q_points; + } -template -void -BoussinesqFlowProblem::Postprocessor:: -compute_derived_quantities_vector (const std::vector > &uh, - const std::vector > > &duh, - const std::vector > > &/*dduh*/, - const std::vector > &/*normals*/, - const std::vector > &evaluation_points, - std::vector > &computed_quantities) const -{ - const unsigned int n_quadrature_points = uh.size(); - Assert (duh.size() == n_quadrature_points, ExcInternalError()); - Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError()); - Assert (uh[0].size() == dim+2, ExcInternalError()); - Assert (computed_quantities[0].size()==n_output_variables(),ExcInternalError()); + template + void + BoussinesqFlowProblem::Postprocessor:: + compute_derived_quantities_vector (const std::vector > &uh, + const std::vector > > &duh, + const std::vector > > &/*dduh*/, + const std::vector > &/*normals*/, + const std::vector > &evaluation_points, + std::vector > &computed_quantities) const + { + const unsigned int n_quadrature_points = uh.size(); + Assert (duh.size() == n_quadrature_points, ExcInternalError()); + Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError()); + Assert (uh[0].size() == dim+2, ExcInternalError()); + Assert (computed_quantities[0].size()==n_output_variables(),ExcInternalError()); - for (unsigned int q=0; q grad_u; - for (unsigned int d=0; d strain_rate = symmetrize (grad_u); - computed_quantities[q](dim+2) = 2 * EquationData::eta * - strain_rate * strain_rate; - - computed_quantities[q](dim+3) = partition; - - computed_quantities[q](dim+4) = pressure - - EquationData::adiabatic_pressure (evaluation_points[q]); - - computed_quantities[q](dim+5) = temperature - - EquationData::adiabatic_temperature (evaluation_points[q]); - } -} + for (unsigned int q=0; q grad_u; + for (unsigned int d=0; d strain_rate = symmetrize (grad_u); + computed_quantities[q](dim+2) = 2 * EquationData::eta * + strain_rate * strain_rate; + + computed_quantities[q](dim+3) = partition; + + computed_quantities[q](dim+4) = pressure - + EquationData::adiabatic_pressure (evaluation_points[q]); + + computed_quantities[q](dim+5) = temperature - + EquationData::adiabatic_temperature (evaluation_points[q]); + } + } // This function does mostly what the @@ -3545,121 +3547,121 @@ compute_derived_quantities_vector (const std::vector > // therefore rescale things into // centimeters per year, the unit // commonly used in geophysics. -template -void BoussinesqFlowProblem::output_results () -{ - computing_timer.enter_section ("Postprocessing"); + template + void BoussinesqFlowProblem::output_results () + { + computing_timer.enter_section ("Postprocessing"); - const FESystem joint_fe (stokes_fe, 1, - temperature_fe, 1); + const FESystem joint_fe (stokes_fe, 1, + temperature_fe, 1); - DoFHandler joint_dof_handler (triangulation); - joint_dof_handler.distribute_dofs (joint_fe); - Assert (joint_dof_handler.n_dofs() == - stokes_dof_handler.n_dofs() + temperature_dof_handler.n_dofs(), - ExcInternalError()); + DoFHandler joint_dof_handler (triangulation); + joint_dof_handler.distribute_dofs (joint_fe); + Assert (joint_dof_handler.n_dofs() == + stokes_dof_handler.n_dofs() + temperature_dof_handler.n_dofs(), + ExcInternalError()); - TrilinosWrappers::MPI::Vector joint_solution; - joint_solution.reinit (joint_dof_handler.locally_owned_dofs(), MPI_COMM_WORLD); + TrilinosWrappers::MPI::Vector joint_solution; + joint_solution.reinit (joint_dof_handler.locally_owned_dofs(), MPI_COMM_WORLD); - { - 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); + { + 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 (joint_cell->is_locally_owned()) - { - 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::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 (joint_cell->is_locally_owned()) + { + 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 data_out; - data_out.attach_dof_handler (joint_dof_handler); - data_out.add_data_vector (locally_relevant_joint_solution, postprocessor); - data_out.build_patches (); - - static int out_index=0; - const std::string filename = ("solution-" + - Utilities::int_to_string (out_index, 5) + - "." + - Utilities::int_to_string - (triangulation.locally_owned_subdomain(), 4) + - ".vtu"); - std::ofstream output (filename.c_str()); - data_out.write_vtu (output); - - if (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD) == 0) - { - std::vector filenames; - for (unsigned int i=0; i data_out; + data_out.attach_dof_handler (joint_dof_handler); + data_out.add_data_vector (locally_relevant_joint_solution, postprocessor); + data_out.build_patches (); + + static int out_index=0; + const std::string filename = ("solution-" + + Utilities::int_to_string (out_index, 5) + + "." + + Utilities::int_to_string + (triangulation.locally_owned_subdomain(), 4) + + ".vtu"); + std::ofstream output (filename.c_str()); + data_out.write_vtu (output); + + if (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD) == 0) + { + std::vector filenames; + for (unsigned int i=0; i::output_results () // complete set of refinement // indicators, and the rest of the // function proceeds as before. -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()); - - KellyErrorEstimator::estimate (temperature_dof_handler, - QGauss(parameters.temperature_degree+1), - typename FunctionMap::type(), - temperature_solution, - estimated_error_per_cell, - std::vector(), - 0, - 0, - triangulation.locally_owned_subdomain()); - - parallel::distributed::GridRefinement:: - refine_and_coarsen_fixed_fraction (triangulation, - estimated_error_per_cell, - 0.3, 0.1); - - // limit maximum refinement level - if (triangulation.n_levels() > max_grid_level) - for (typename Triangulation::active_cell_iterator - cell = triangulation.begin_active(max_grid_level); - cell != triangulation.end(); ++cell) - cell->clear_refine_flag (); - - std::vector x_temperature (2); - x_temperature[0] = &temperature_solution; - x_temperature[1] = &old_temperature_solution; - std::vector x_stokes (2); - x_stokes[0] = &stokes_solution; - x_stokes[1] = &old_stokes_solution; - - parallel::distributed::SolutionTransfer - temperature_trans(temperature_dof_handler); - parallel::distributed::SolutionTransfer - stokes_trans(stokes_dof_handler); - - triangulation.prepare_coarsening_and_refinement(); - 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 (); - - computing_timer.enter_section ("Refine mesh structure, part 2"); - + template + void BoussinesqFlowProblem::refine_mesh (const unsigned int max_grid_level) { - TrilinosWrappers::MPI::Vector - distributed_temp1 (temperature_rhs); - TrilinosWrappers::MPI::Vector - distributed_temp2 (temperature_rhs); + computing_timer.enter_section ("Refine mesh structure, part 1"); + Vector estimated_error_per_cell (triangulation.n_active_cells()); + + KellyErrorEstimator::estimate (temperature_dof_handler, + QGauss(parameters.temperature_degree+1), + typename FunctionMap::type(), + temperature_solution, + estimated_error_per_cell, + std::vector(), + 0, + 0, + triangulation.locally_owned_subdomain()); + + parallel::distributed::GridRefinement:: + refine_and_coarsen_fixed_fraction (triangulation, + estimated_error_per_cell, + 0.3, 0.1); + + // limit maximum refinement level + if (triangulation.n_levels() > max_grid_level) + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(max_grid_level); + cell != triangulation.end(); ++cell) + cell->clear_refine_flag (); + + std::vector x_temperature (2); + x_temperature[0] = &temperature_solution; + x_temperature[1] = &old_temperature_solution; + std::vector x_stokes (2); + x_stokes[0] = &stokes_solution; + x_stokes[1] = &old_stokes_solution; + + parallel::distributed::SolutionTransfer + temperature_trans(temperature_dof_handler); + parallel::distributed::SolutionTransfer + stokes_trans(stokes_dof_handler); + + triangulation.prepare_coarsening_and_refinement(); + 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(); - std::vector tmp (2); - tmp[0] = &(distributed_temp1); - tmp[1] = &(distributed_temp2); - temperature_trans.interpolate(tmp); + setup_dofs (); - // temperature_solution = distributed_temp1; - temperature_solution.reinit(distributed_temp1, false, true); - // old_temperature_solution = distributed_temp2; - old_temperature_solution.reinit(distributed_temp2, false, true); - } + computing_timer.enter_section ("Refine mesh structure, part 2"); - { - TrilinosWrappers::MPI::BlockVector - distributed_stokes (stokes_rhs); - TrilinosWrappers::MPI::BlockVector - old_distributed_stokes (stokes_rhs); - std::vector stokes_tmp (2); - stokes_tmp[0] = &(distributed_stokes); - stokes_tmp[1] = &(old_distributed_stokes); - - stokes_trans.interpolate (stokes_tmp); - // stokes_solution = distributed_stokes; - stokes_solution.block(0).reinit(distributed_stokes.block(0), false, true); - stokes_solution.block(1).reinit(distributed_stokes.block(1), false, true); - old_stokes_solution.block(0).reinit(old_distributed_stokes.block(0), false, true); - old_stokes_solution.block(1).reinit(old_distributed_stokes.block(1), false, true); - } + { + TrilinosWrappers::MPI::Vector + distributed_temp1 (temperature_rhs); + TrilinosWrappers::MPI::Vector + distributed_temp2 (temperature_rhs); + + std::vector tmp (2); + tmp[0] = &(distributed_temp1); + tmp[1] = &(distributed_temp2); + temperature_trans.interpolate(tmp); + + // temperature_solution = distributed_temp1; + temperature_solution.reinit(distributed_temp1, false, true); + // old_temperature_solution = distributed_temp2; + old_temperature_solution.reinit(distributed_temp2, false, true); + } - computing_timer.exit_section(); -} + { + TrilinosWrappers::MPI::BlockVector + distributed_stokes (stokes_rhs); + TrilinosWrappers::MPI::BlockVector + old_distributed_stokes (stokes_rhs); + std::vector stokes_tmp (2); + stokes_tmp[0] = &(distributed_stokes); + stokes_tmp[1] = &(old_distributed_stokes); + + stokes_trans.interpolate (stokes_tmp); + // stokes_solution = distributed_stokes; + stokes_solution.block(0).reinit(distributed_stokes.block(0), false, true); + stokes_solution.block(1).reinit(distributed_stokes.block(1), false, true); + old_stokes_solution.block(0).reinit(old_distributed_stokes.block(0), false, true); + old_stokes_solution.block(1).reinit(old_distributed_stokes.block(1), false, true); + } + + computing_timer.exit_section(); + } @@ -3854,111 +3856,112 @@ void BoussinesqFlowProblem::refine_mesh (const unsigned int max_grid_level) // function instead of the library function // VectorTools::project, the // rest is as before. -template -void BoussinesqFlowProblem::run () -{ - GridGenerator::hyper_shell (triangulation, - Point(), - EquationData::R0, - EquationData::R1, - (dim==3) ? 96 : 12, - true); - static HyperShellBoundary boundary; - triangulation.set_boundary (0, boundary); - triangulation.set_boundary (1, boundary); + template + void BoussinesqFlowProblem::run () + { + GridGenerator::hyper_shell (triangulation, + Point(), + EquationData::R0, + EquationData::R1, + (dim==3) ? 96 : 12, + true); + static HyperShellBoundary boundary; + triangulation.set_boundary (0, boundary); + triangulation.set_boundary (1, boundary); - global_Omega_diameter = GridTools::diameter (triangulation); + global_Omega_diameter = GridTools::diameter (triangulation); - triangulation.refine_global (parameters.initial_global_refinement); + triangulation.refine_global (parameters.initial_global_refinement); - setup_dofs(); + setup_dofs(); - unsigned int pre_refinement_step = 0; + unsigned int pre_refinement_step = 0; -start_time_iteration: + start_time_iteration: - project_temperature_field (); + project_temperature_field (); - timestep_number = 0; - time_step = old_time_step = 0; + timestep_number = 0; + time_step = old_time_step = 0; - double time = 0; + double time = 0; - do - { - pcout << "Timestep " << timestep_number - << ": t=" << time/EquationData::year_in_seconds - << " years" - << std::endl; - - assemble_stokes_system (); - build_stokes_preconditioner (); - assemble_temperature_matrix (); - - solve (); - - pcout << std::endl; - - if ((timestep_number == 0) && - (pre_refinement_step < parameters.initial_adaptive_refinement)) - { - refine_mesh (parameters.initial_global_refinement + - parameters.initial_adaptive_refinement); - ++pre_refinement_step; - goto start_time_iteration; - } - else if ((timestep_number > 0) - && - (timestep_number % parameters.adaptive_refinement_interval == 0)) - refine_mesh (parameters.initial_global_refinement + - parameters.initial_adaptive_refinement); - - if ((parameters.generate_graphical_output == true) - && - (timestep_number % parameters.graphical_output_interval == 0)) - output_results (); - - time += time_step; - ++timestep_number; - - // if we are at the end of - // time, stop now - if (time > parameters.end_time * EquationData::year_in_seconds) - break; - - // otherwise prepare for the - // next time step - TrilinosWrappers::MPI::BlockVector old_old_stokes_solution; - old_old_stokes_solution = old_stokes_solution; - old_stokes_solution = stokes_solution; - old_old_temperature_solution = old_temperature_solution; - old_temperature_solution = temperature_solution; - if (old_time_step > 0) - { - stokes_solution.sadd (1.+time_step/old_time_step, -time_step/old_time_step, - old_old_stokes_solution); - temperature_solution.sadd (1.+time_step/old_time_step, - -time_step/old_time_step, - old_old_temperature_solution); - } - - // every 100 time steps output - // a summary of the current - // timing information - if (timestep_number % 100 == 0) - computing_timer.print_summary (); - } - while (true); - - // if we are generating graphical - // output, do so also for the last - // time step unless we had just - // done so before we left the - // do-while loop - if ((parameters.generate_graphical_output == true) - && - !((timestep_number-1) % parameters.graphical_output_interval == 0)) - output_results (); + do + { + pcout << "Timestep " << timestep_number + << ": t=" << time/EquationData::year_in_seconds + << " years" + << std::endl; + + assemble_stokes_system (); + build_stokes_preconditioner (); + assemble_temperature_matrix (); + + solve (); + + pcout << std::endl; + + if ((timestep_number == 0) && + (pre_refinement_step < parameters.initial_adaptive_refinement)) + { + refine_mesh (parameters.initial_global_refinement + + parameters.initial_adaptive_refinement); + ++pre_refinement_step; + goto start_time_iteration; + } + else if ((timestep_number > 0) + && + (timestep_number % parameters.adaptive_refinement_interval == 0)) + refine_mesh (parameters.initial_global_refinement + + parameters.initial_adaptive_refinement); + + if ((parameters.generate_graphical_output == true) + && + (timestep_number % parameters.graphical_output_interval == 0)) + output_results (); + + time += time_step; + ++timestep_number; + + // if we are at the end of + // time, stop now + if (time > parameters.end_time * EquationData::year_in_seconds) + break; + + // otherwise prepare for the + // next time step + TrilinosWrappers::MPI::BlockVector old_old_stokes_solution; + old_old_stokes_solution = old_stokes_solution; + old_stokes_solution = stokes_solution; + old_old_temperature_solution = old_temperature_solution; + old_temperature_solution = temperature_solution; + if (old_time_step > 0) + { + stokes_solution.sadd (1.+time_step/old_time_step, -time_step/old_time_step, + old_old_stokes_solution); + temperature_solution.sadd (1.+time_step/old_time_step, + -time_step/old_time_step, + old_old_temperature_solution); + } + + // every 100 time steps output + // a summary of the current + // timing information + if (timestep_number % 100 == 0) + computing_timer.print_summary (); + } + while (true); + + // if we are generating graphical + // output, do so also for the last + // time step unless we had just + // done so before we left the + // do-while loop + if ((parameters.generate_graphical_output == true) + && + !((timestep_number-1) % parameters.graphical_output_interval == 0)) + output_results (); + } } @@ -3968,6 +3971,9 @@ start_time_iteration: // This is copied verbatim from step-31: int main (int argc, char *argv[]) { + using namespace Step32; + using namespace dealii; + Utilities::System::MPI_InitFinalize mpi_initialization(argc, argv); try -- 2.39.5