-using namespace dealii;
+namespace Step32
+{
+ using namespace dealii;
// @sect3{Equation data}
// 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 <int dim>
- Tensor<1,dim> gravity_vector (const Point<dim> &p)
- {
+ template <int dim>
+ Tensor<1,dim> gravity_vector (const Point<dim> &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 <int dim>
- double adiabatic_pressure (const Point<dim> &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 <int dim>
+ double adiabatic_pressure (const Point<dim> &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 <int dim>
- double adiabatic_temperature (const Point<dim> &p)
- {
- // the static, adiabatic
- // temperature satisfies
- // $dT/dr = -T alpha/c_P g$
+ template <int dim>
+ double adiabatic_temperature (const Point<dim> &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 <int dim>
- class TemperatureInitialValues : public Function<dim>
- {
- public:
- TemperatureInitialValues () : Function<dim>(1) {}
+ template <int dim>
+ class TemperatureInitialValues : public Function<dim>
+ {
+ public:
+ TemperatureInitialValues () : Function<dim>(1) {}
- virtual double value (const Point<dim> &p,
- const unsigned int component = 0) const;
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
- virtual void vector_value (const Point<dim> &p,
- Vector<double> &value) const;
- };
+ virtual void vector_value (const Point<dim> &p,
+ Vector<double> &value) const;
+ };
- template <int dim>
- double
- TemperatureInitialValues<dim>::value (const Point<dim> &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 <int dim>
+ double
+ TemperatureInitialValues<dim>::value (const Point<dim> &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 <int dim>
- void
- TemperatureInitialValues<dim>::vector_value (const Point<dim> &p,
- Vector<double> &values) const
- {
- for (unsigned int c=0; c<this->n_components; ++c)
- values(c) = TemperatureInitialValues<dim>::value (p, c);
+ template <int dim>
+ void
+ TemperatureInitialValues<dim>::vector_value (const Point<dim> &p,
+ Vector<double> &values) const
+ {
+ for (unsigned int c=0; c<this->n_components; ++c)
+ values(c) = TemperatureInitialValues<dim>::value (p, c);
+ }
}
-}
// respective variable in the
// BlockSchurPreconditioner class another
// name.
-namespace LinearSolvers
-{
- template <class PreconditionerA, class PreconditionerMp>
- 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<const TrilinosWrappers::BlockSparseMatrix> stokes_matrix;
- const SmartPointer<const TrilinosWrappers::BlockSparseMatrix> stokes_preconditioner_matrix;
- const PreconditionerMp &mp_preconditioner;
- const PreconditionerA &a_preconditioner;
- const bool do_solve_A;
- };
-}
+ template <class PreconditionerA, class PreconditionerMp>
+ 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<const TrilinosWrappers::BlockSparseMatrix> stokes_matrix;
+ const SmartPointer<const TrilinosWrappers::BlockSparseMatrix> stokes_preconditioner_matrix;
+ const PreconditionerMp &mp_preconditioner;
+ const PreconditionerA &a_preconditioner;
+ const bool do_solve_A;
+ };
+ }
// itself), and provide some additional
// vector fields that are used to improve
// performance of assembly.
-namespace Assembly
-{
- namespace Scratch
+ namespace Assembly
{
- template <int dim>
- struct StokesPreconditioner
+ namespace Scratch
{
- StokesPreconditioner (const FiniteElement<dim> &stokes_fe,
- const Quadrature<dim> &stokes_quadrature,
- const Mapping<dim> &mapping,
- const UpdateFlags update_flags);
- StokesPreconditioner (const StokesPreconditioner &data);
-
- FEValues<dim> stokes_fe_values;
-
- std::vector<SymmetricTensor<2,dim> > grads_phi_u;
- std::vector<double> phi_p;
- };
-
- template <int dim>
- StokesPreconditioner<dim>::
- StokesPreconditioner (const FiniteElement<dim> &stokes_fe,
- const Quadrature<dim> &stokes_quadrature,
- const Mapping<dim> &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 <int dim>
+ struct StokesPreconditioner
+ {
+ StokesPreconditioner (const FiniteElement<dim> &stokes_fe,
+ const Quadrature<dim> &stokes_quadrature,
+ const Mapping<dim> &mapping,
+ const UpdateFlags update_flags);
+ StokesPreconditioner (const StokesPreconditioner &data);
+ FEValues<dim> stokes_fe_values;
+ std::vector<SymmetricTensor<2,dim> > grads_phi_u;
+ std::vector<double> phi_p;
+ };
- template <int dim>
- StokesPreconditioner<dim>::
- 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 <int dim>
- struct StokesSystem : public StokesPreconditioner<dim>
- {
- StokesSystem (const FiniteElement<dim> &stokes_fe,
- const Mapping<dim> &mapping,
- const Quadrature<dim> &stokes_quadrature,
- const UpdateFlags stokes_update_flags,
- const FiniteElement<dim> &temperature_fe,
- const UpdateFlags temperature_update_flags);
+ template <int dim>
+ StokesPreconditioner<dim>::
+ StokesPreconditioner (const FiniteElement<dim> &stokes_fe,
+ const Quadrature<dim> &stokes_quadrature,
+ const Mapping<dim> &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<dim> &data);
- FEValues<dim> temperature_fe_values;
- std::vector<Tensor<1,dim> > phi_u;
- std::vector<SymmetricTensor<2,dim> > grads_phi_u;
- std::vector<double> div_phi_u;
+ template <int dim>
+ StokesPreconditioner<dim>::
+ 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<double> old_temperature_values;
- };
- template <int dim>
- StokesSystem<dim>::
- StokesSystem (const FiniteElement<dim> &stokes_fe,
- const Mapping<dim> &mapping,
- const Quadrature<dim> &stokes_quadrature,
- const UpdateFlags stokes_update_flags,
- const FiniteElement<dim> &temperature_fe,
- const UpdateFlags temperature_update_flags)
- :
- StokesPreconditioner<dim> (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 <int dim>
+ struct StokesSystem : public StokesPreconditioner<dim>
+ {
+ StokesSystem (const FiniteElement<dim> &stokes_fe,
+ const Mapping<dim> &mapping,
+ const Quadrature<dim> &stokes_quadrature,
+ const UpdateFlags stokes_update_flags,
+ const FiniteElement<dim> &temperature_fe,
+ const UpdateFlags temperature_update_flags);
+ StokesSystem (const StokesSystem<dim> &data);
- template <int dim>
- StokesSystem<dim>::
- StokesSystem (const StokesSystem<dim> &scratch)
- :
- StokesPreconditioner<dim> (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<dim> temperature_fe_values;
+ std::vector<Tensor<1,dim> > phi_u;
+ std::vector<SymmetricTensor<2,dim> > grads_phi_u;
+ std::vector<double> div_phi_u;
+ std::vector<double> old_temperature_values;
+ };
- template <int dim>
- struct TemperatureMatrix
- {
- TemperatureMatrix (const FiniteElement<dim> &temperature_fe,
- const Mapping<dim> &mapping,
- const Quadrature<dim> &temperature_quadrature);
- TemperatureMatrix (const TemperatureMatrix &data);
- FEValues<dim> temperature_fe_values;
+ template <int dim>
+ StokesSystem<dim>::
+ StokesSystem (const FiniteElement<dim> &stokes_fe,
+ const Mapping<dim> &mapping,
+ const Quadrature<dim> &stokes_quadrature,
+ const UpdateFlags stokes_update_flags,
+ const FiniteElement<dim> &temperature_fe,
+ const UpdateFlags temperature_update_flags)
+ :
+ StokesPreconditioner<dim> (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<double> phi_T;
- std::vector<Tensor<1,dim> > grad_phi_T;
- };
- template <int dim>
- TemperatureMatrix<dim>::
- TemperatureMatrix (const FiniteElement<dim> &temperature_fe,
- const Mapping<dim> &mapping,
- const Quadrature<dim> &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 <int dim>
+ StokesSystem<dim>::
+ StokesSystem (const StokesSystem<dim> &scratch)
+ :
+ StokesPreconditioner<dim> (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 <int dim>
- TemperatureMatrix<dim>::
- 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 <int dim>
+ struct TemperatureMatrix
+ {
+ TemperatureMatrix (const FiniteElement<dim> &temperature_fe,
+ const Mapping<dim> &mapping,
+ const Quadrature<dim> &temperature_quadrature);
+ TemperatureMatrix (const TemperatureMatrix &data);
- template <int dim>
- struct TemperatureRHS
- {
- TemperatureRHS (const FiniteElement<dim> &temperature_fe,
- const FiniteElement<dim> &stokes_fe,
- const Mapping<dim> &mapping,
- const Quadrature<dim> &quadrature);
- TemperatureRHS (const TemperatureRHS &data);
+ FEValues<dim> temperature_fe_values;
- FEValues<dim> temperature_fe_values;
- FEValues<dim> stokes_fe_values;
+ std::vector<double> phi_T;
+ std::vector<Tensor<1,dim> > grad_phi_T;
+ };
- std::vector<double> phi_T;
- std::vector<Tensor<1,dim> > grad_phi_T;
+ template <int dim>
+ TemperatureMatrix<dim>::
+ TemperatureMatrix (const FiniteElement<dim> &temperature_fe,
+ const Mapping<dim> &mapping,
+ const Quadrature<dim> &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<Tensor<1,dim> > old_velocity_values;
- std::vector<Tensor<1,dim> > old_old_velocity_values;
- std::vector<SymmetricTensor<2,dim> > old_strain_rates;
- std::vector<SymmetricTensor<2,dim> > old_old_strain_rates;
+ template <int dim>
+ TemperatureMatrix<dim>::
+ 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<double> old_temperature_values;
- std::vector<double> old_old_temperature_values;
- std::vector<Tensor<1,dim> > old_temperature_grads;
- std::vector<Tensor<1,dim> > old_old_temperature_grads;
- std::vector<double> old_temperature_laplacians;
- std::vector<double> old_old_temperature_laplacians;
- };
- template <int dim>
- TemperatureRHS<dim>::
- TemperatureRHS (const FiniteElement<dim> &temperature_fe,
- const FiniteElement<dim> &stokes_fe,
- const Mapping<dim> &mapping,
- const Quadrature<dim> &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 <int dim>
+ struct TemperatureRHS
+ {
+ TemperatureRHS (const FiniteElement<dim> &temperature_fe,
+ const FiniteElement<dim> &stokes_fe,
+ const Mapping<dim> &mapping,
+ const Quadrature<dim> &quadrature);
+ TemperatureRHS (const TemperatureRHS &data);
+
+ FEValues<dim> temperature_fe_values;
+ FEValues<dim> stokes_fe_values;
+
+ std::vector<double> phi_T;
+ std::vector<Tensor<1,dim> > grad_phi_T;
+
+ std::vector<Tensor<1,dim> > old_velocity_values;
+ std::vector<Tensor<1,dim> > old_old_velocity_values;
+
+ std::vector<SymmetricTensor<2,dim> > old_strain_rates;
+ std::vector<SymmetricTensor<2,dim> > old_old_strain_rates;
+
+ std::vector<double> old_temperature_values;
+ std::vector<double> old_old_temperature_values;
+ std::vector<Tensor<1,dim> > old_temperature_grads;
+ std::vector<Tensor<1,dim> > old_old_temperature_grads;
+ std::vector<double> old_temperature_laplacians;
+ std::vector<double> old_old_temperature_laplacians;
+ };
+
+ template <int dim>
+ TemperatureRHS<dim>::
+ TemperatureRHS (const FiniteElement<dim> &temperature_fe,
+ const FiniteElement<dim> &stokes_fe,
+ const Mapping<dim> &mapping,
+ const Quadrature<dim> &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 <int dim>
- TemperatureRHS<dim>::
- 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 <int dim>
+ TemperatureRHS<dim>::
+ 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.
- // <code>local_dof_indices</code>).
- namespace CopyData
- {
- template <int dim>
- 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.
+ // <code>local_dof_indices</code>).
+ namespace CopyData
{
- StokesPreconditioner (const FiniteElement<dim> &stokes_fe);
- StokesPreconditioner (const StokesPreconditioner &data);
-
- FullMatrix<double> local_matrix;
- std::vector<unsigned int> local_dof_indices;
- };
-
- template <int dim>
- StokesPreconditioner<dim>::
- StokesPreconditioner (const FiniteElement<dim> &stokes_fe)
- :
- local_matrix (stokes_fe.dofs_per_cell,
- stokes_fe.dofs_per_cell),
- local_dof_indices (stokes_fe.dofs_per_cell)
- {}
-
-
-
- template <int dim>
- StokesPreconditioner<dim>::
- StokesPreconditioner (const StokesPreconditioner &data)
- :
- local_matrix (data.local_matrix),
- local_dof_indices (data.local_dof_indices)
- {}
-
+ template <int dim>
+ struct StokesPreconditioner
+ {
+ StokesPreconditioner (const FiniteElement<dim> &stokes_fe);
+ StokesPreconditioner (const StokesPreconditioner &data);
+
+ FullMatrix<double> local_matrix;
+ std::vector<unsigned int> local_dof_indices;
+ };
+
+ template <int dim>
+ StokesPreconditioner<dim>::
+ StokesPreconditioner (const FiniteElement<dim> &stokes_fe)
+ :
+ local_matrix (stokes_fe.dofs_per_cell,
+ stokes_fe.dofs_per_cell),
+ local_dof_indices (stokes_fe.dofs_per_cell)
+ {}
- template <int dim>
- struct StokesSystem : public StokesPreconditioner<dim>
- {
- StokesSystem (const FiniteElement<dim> &stokes_fe);
- StokesSystem (const StokesSystem<dim> &data);
- Vector<double> local_rhs;
- };
+ template <int dim>
+ StokesPreconditioner<dim>::
+ StokesPreconditioner (const StokesPreconditioner &data)
+ :
+ local_matrix (data.local_matrix),
+ local_dof_indices (data.local_dof_indices)
+ {}
- template <int dim>
- StokesSystem<dim>::
- StokesSystem (const FiniteElement<dim> &stokes_fe)
- :
- StokesPreconditioner<dim> (stokes_fe),
- local_rhs (stokes_fe.dofs_per_cell)
- {}
+ template <int dim>
+ struct StokesSystem : public StokesPreconditioner<dim>
+ {
+ StokesSystem (const FiniteElement<dim> &stokes_fe);
+ StokesSystem (const StokesSystem<dim> &data);
- template <int dim>
- StokesSystem<dim>::
- StokesSystem (const StokesSystem<dim> &data)
- :
- StokesPreconditioner<dim> (data),
- local_rhs (data.local_rhs)
- {}
+ Vector<double> local_rhs;
+ };
+ template <int dim>
+ StokesSystem<dim>::
+ StokesSystem (const FiniteElement<dim> &stokes_fe)
+ :
+ StokesPreconditioner<dim> (stokes_fe),
+ local_rhs (stokes_fe.dofs_per_cell)
+ {}
- template <int dim>
- struct TemperatureMatrix
- {
- TemperatureMatrix (const FiniteElement<dim> &temperature_fe);
- TemperatureMatrix (const TemperatureMatrix &data);
- FullMatrix<double> local_mass_matrix;
- FullMatrix<double> local_stiffness_matrix;
- std::vector<unsigned int> local_dof_indices;
- };
+ template <int dim>
+ StokesSystem<dim>::
+ StokesSystem (const StokesSystem<dim> &data)
+ :
+ StokesPreconditioner<dim> (data),
+ local_rhs (data.local_rhs)
+ {}
- template <int dim>
- TemperatureMatrix<dim>::
- TemperatureMatrix (const FiniteElement<dim> &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 <int dim>
- TemperatureMatrix<dim>::
- 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 <int dim>
+ struct TemperatureMatrix
+ {
+ TemperatureMatrix (const FiniteElement<dim> &temperature_fe);
+ TemperatureMatrix (const TemperatureMatrix &data);
+
+ FullMatrix<double> local_mass_matrix;
+ FullMatrix<double> local_stiffness_matrix;
+ std::vector<unsigned int> local_dof_indices;
+ };
+
+ template <int dim>
+ TemperatureMatrix<dim>::
+ TemperatureMatrix (const FiniteElement<dim> &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 <int dim>
- struct TemperatureRHS
- {
- TemperatureRHS (const FiniteElement<dim> &temperature_fe);
- TemperatureRHS (const TemperatureRHS &data);
+ template <int dim>
+ TemperatureMatrix<dim>::
+ 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<double> local_rhs;
- std::vector<unsigned int> local_dof_indices;
- FullMatrix<double> matrix_for_bc;
- };
- template <int dim>
- TemperatureRHS<dim>::
- TemperatureRHS (const FiniteElement<dim> &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 <int dim>
+ struct TemperatureRHS
+ {
+ TemperatureRHS (const FiniteElement<dim> &temperature_fe);
+ TemperatureRHS (const TemperatureRHS &data);
+
+ Vector<double> local_rhs;
+ std::vector<unsigned int> local_dof_indices;
+ FullMatrix<double> matrix_for_bc;
+ };
+
+ template <int dim>
+ TemperatureRHS<dim>::
+ TemperatureRHS (const FiniteElement<dim> &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 <int dim>
- TemperatureRHS<dim>::
- TemperatureRHS (const TemperatureRHS &data)
- :
- local_rhs (data.local_rhs),
- local_dof_indices (data.local_dof_indices),
- matrix_for_bc (data.matrix_for_bc)
- {}
+ template <int dim>
+ TemperatureRHS<dim>::
+ TemperatureRHS (const TemperatureRHS &data)
+ :
+ local_rhs (data.local_rhs),
+ local_dof_indices (data.local_dof_indices),
+ matrix_for_bc (data.matrix_for_bc)
+ {}
+ }
}
-}
// program (this choice is made in the
// from this variable in the results
// section of this tutorial program.
-template <int dim>
-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<double,double> get_extrapolated_temperature_range () const;
- void solve ();
- void output_results ();
- void refine_mesh (const unsigned int max_grid_level);
+ template <int dim>
+ class BoussinesqFlowProblem
+ {
+ public:
+ struct Parameters;
+ BoussinesqFlowProblem (Parameters ¶meters);
+ void run ();
- double
- compute_viscosity(const std::vector<double> &old_temperature,
- const std::vector<double> &old_old_temperature,
- const std::vector<Tensor<1,dim> > &old_temperature_grads,
- const std::vector<Tensor<1,dim> > &old_old_temperature_grads,
- const std::vector<double> &old_temperature_laplacians,
- const std::vector<double> &old_old_temperature_laplacians,
- const std::vector<Tensor<1,dim> > &old_velocity_values,
- const std::vector<Tensor<1,dim> > &old_old_velocity_values,
- const std::vector<SymmetricTensor<2,dim> > &old_strain_rates,
- const std::vector<SymmetricTensor<2,dim> > &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<double,double> 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<double> &old_temperature,
+ const std::vector<double> &old_old_temperature,
+ const std::vector<Tensor<1,dim> > &old_temperature_grads,
+ const std::vector<Tensor<1,dim> > &old_old_temperature_grads,
+ const std::vector<double> &old_temperature_laplacians,
+ const std::vector<double> &old_old_temperature_laplacians,
+ const std::vector<Tensor<1,dim> > &old_velocity_values,
+ const std::vector<Tensor<1,dim> > &old_old_velocity_values,
+ const std::vector<SymmetricTensor<2,dim> > &old_strain_rates,
+ const std::vector<SymmetricTensor<2,dim> > &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<dim> triangulation;
- double global_Omega_diameter;
+ private:
+ Parameters ¶meters;
+ ConditionalOStream pcout;
- const MappingQ<dim> mapping;
+ parallel::distributed::Triangulation<dim> triangulation;
+ double global_Omega_diameter;
- const FESystem<dim> stokes_fe;
+ const MappingQ<dim> mapping;
- DoFHandler<dim> stokes_dof_handler;
- ConstraintMatrix stokes_constraints;
+ const FESystem<dim> stokes_fe;
- TrilinosWrappers::BlockSparseMatrix stokes_matrix;
- TrilinosWrappers::BlockSparseMatrix stokes_preconditioner_matrix;
+ DoFHandler<dim> 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<dim> temperature_fe;
- DoFHandler<dim> temperature_dof_handler;
- ConstraintMatrix temperature_constraints;
- TrilinosWrappers::SparseMatrix temperature_mass_matrix;
- TrilinosWrappers::SparseMatrix temperature_stiffness_matrix;
- TrilinosWrappers::SparseMatrix temperature_matrix;
+ FE_Q<dim> temperature_fe;
+ DoFHandler<dim> 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<TrilinosWrappers::PreconditionAMG> Amg_preconditioner;
- std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionILU> Mp_preconditioner;
- std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionIC> 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<TrilinosWrappers::PreconditionAMG> Amg_preconditioner;
+ std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionILU> Mp_preconditioner;
+ std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionIC> 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<IndexSet> &stokes_partitioning);
- void setup_stokes_preconditioner (const std::vector<IndexSet> &stokes_partitioning);
- void setup_temperature_matrices (const IndexSet &temperature_partitioning);
+ TimerOutput computing_timer;
- void
- local_assemble_stokes_preconditioner (const typename DoFHandler<dim>::active_cell_iterator &cell,
- Assembly::Scratch::StokesPreconditioner<dim> &scratch,
- Assembly::CopyData::StokesPreconditioner<dim> &data);
+ void setup_stokes_matrix (const std::vector<IndexSet> &stokes_partitioning);
+ void setup_stokes_preconditioner (const std::vector<IndexSet> &stokes_partitioning);
+ void setup_temperature_matrices (const IndexSet &temperature_partitioning);
- void
- copy_local_to_global_stokes_preconditioner (const Assembly::CopyData::StokesPreconditioner<dim> &data);
+ void
+ local_assemble_stokes_preconditioner (const typename DoFHandler<dim>::active_cell_iterator &cell,
+ Assembly::Scratch::StokesPreconditioner<dim> &scratch,
+ Assembly::CopyData::StokesPreconditioner<dim> &data);
+ void
+ copy_local_to_global_stokes_preconditioner (const Assembly::CopyData::StokesPreconditioner<dim> &data);
- void
- local_assemble_stokes_system (const typename DoFHandler<dim>::active_cell_iterator &cell,
- Assembly::Scratch::StokesSystem<dim> &scratch,
- Assembly::CopyData::StokesSystem<dim> &data);
- void
- copy_local_to_global_stokes_system (const Assembly::CopyData::StokesSystem<dim> &data);
+ void
+ local_assemble_stokes_system (const typename DoFHandler<dim>::active_cell_iterator &cell,
+ Assembly::Scratch::StokesSystem<dim> &scratch,
+ Assembly::CopyData::StokesSystem<dim> &data);
+ void
+ copy_local_to_global_stokes_system (const Assembly::CopyData::StokesSystem<dim> &data);
- void
- local_assemble_temperature_matrix (const typename DoFHandler<dim>::active_cell_iterator &cell,
- Assembly::Scratch::TemperatureMatrix<dim> &scratch,
- Assembly::CopyData::TemperatureMatrix<dim> &data);
- void
- copy_local_to_global_temperature_matrix (const Assembly::CopyData::TemperatureMatrix<dim> &data);
+ void
+ local_assemble_temperature_matrix (const typename DoFHandler<dim>::active_cell_iterator &cell,
+ Assembly::Scratch::TemperatureMatrix<dim> &scratch,
+ Assembly::CopyData::TemperatureMatrix<dim> &data);
+ void
+ copy_local_to_global_temperature_matrix (const Assembly::CopyData::TemperatureMatrix<dim> &data);
- void
- local_assemble_temperature_rhs (const std::pair<double,double> global_T_range,
- const double global_max_velocity,
- const double global_entropy_variation,
- const typename DoFHandler<dim>::active_cell_iterator &cell,
- Assembly::Scratch::TemperatureRHS<dim> &scratch,
- Assembly::CopyData::TemperatureRHS<dim> &data);
- void
- copy_local_to_global_temperature_rhs (const Assembly::CopyData::TemperatureRHS<dim> &data);
+ void
+ local_assemble_temperature_rhs (const std::pair<double,double> global_T_range,
+ const double global_max_velocity,
+ const double global_entropy_variation,
+ const typename DoFHandler<dim>::active_cell_iterator &cell,
+ Assembly::Scratch::TemperatureRHS<dim> &scratch,
+ Assembly::CopyData::TemperatureRHS<dim> &data);
- class Postprocessor;
-};
+ void
+ copy_local_to_global_temperature_rhs (const Assembly::CopyData::TemperatureRHS<dim> &data);
+
+ class Postprocessor;
+ };
// @sect3{BoussinesqFlowProblem class implementation}
// @sect4{BoussinesqFlowProblem::Parameters}
-template <int dim>
-BoussinesqFlowProblem<dim>::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<dim>::Parameters::declare_parameters (prm);
+ template <int dim>
+ BoussinesqFlowProblem<dim>::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<dim>::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 <int dim>
-void
-BoussinesqFlowProblem<dim>::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 <int dim>
+ void
+ BoussinesqFlowProblem<dim>::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 <int dim>
-void
-BoussinesqFlowProblem<dim>::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 <int dim>
+ void
+ BoussinesqFlowProblem<dim>::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 ();
-}
// get a summary table in the end of the
// program which shows us wallclock times
// (as opposed to CPU times).
-template <int dim>
-BoussinesqFlowProblem<dim>::BoussinesqFlowProblem (Parameters ¶meters_)
- :
- parameters (parameters_),
- pcout (std::cout,
- (Utilities::System::
- get_this_mpi_process(MPI_COMM_WORLD)
- == 0)),
-
- triangulation (MPI_COMM_WORLD,
- typename Triangulation<dim>::MeshSmoothing
- (Triangulation<dim>::smoothing_on_refinement |
- Triangulation<dim>::smoothing_on_coarsening)),
-
- mapping (4),
-
- stokes_fe (FE_Q<dim>(parameters.stokes_velocity_degree),
- dim,
- (parameters.use_locally_conservative_discretization
- ?
- static_cast<const FiniteElement<dim> &>
- (FE_DGP<dim>(parameters.stokes_velocity_degree-1))
- :
- static_cast<const FiniteElement<dim> &>
- (FE_Q<dim>(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 <int dim>
+ BoussinesqFlowProblem<dim>::BoussinesqFlowProblem (Parameters ¶meters_)
+ :
+ parameters (parameters_),
+ pcout (std::cout,
+ (Utilities::System::
+ get_this_mpi_process(MPI_COMM_WORLD)
+ == 0)),
+
+ triangulation (MPI_COMM_WORLD,
+ typename Triangulation<dim>::MeshSmoothing
+ (Triangulation<dim>::smoothing_on_refinement |
+ Triangulation<dim>::smoothing_on_coarsening)),
+
+ mapping (4),
+
+ stokes_fe (FE_Q<dim>(parameters.stokes_velocity_degree),
+ dim,
+ (parameters.use_locally_conservative_discretization
+ ?
+ static_cast<const FiniteElement<dim> &>
+ (FE_DGP<dim>(parameters.stokes_velocity_degree-1))
+ :
+ static_cast<const FiniteElement<dim> &>
+ (FE_Q<dim>(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)
+ {}
// global maximum (output), and an
// integer value one that says that
// we only work on one double.
-template <int dim>
-double BoussinesqFlowProblem<dim>::get_maximal_velocity () const
-{
- const QIterated<dim> quadrature_formula (QTrapez<1>(),
- parameters.stokes_velocity_degree);
- const unsigned int n_q_points = quadrature_formula.size();
-
- FEValues<dim> fe_values (mapping, stokes_fe, quadrature_formula, update_values);
- std::vector<Tensor<1,dim> > velocity_values(n_q_points);
+ template <int dim>
+ double BoussinesqFlowProblem<dim>::get_maximal_velocity () const
+ {
+ const QIterated<dim> quadrature_formula (QTrapez<1>(),
+ parameters.stokes_velocity_degree);
+ const unsigned int n_q_points = quadrature_formula.size();
- const FEValuesExtractors::Vector velocities (0);
+ FEValues<dim> fe_values (mapping, stokes_fe, quadrature_formula, update_values);
+ std::vector<Tensor<1,dim> > velocity_values(n_q_points);
- double max_local_velocity = 0;
+ const FEValuesExtractors::Vector velocities (0);
- typename DoFHandler<dim>::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<n_q_points; ++q)
- max_local_velocity = std::max (max_local_velocity,
- velocity_values[q].norm());
- }
+ double max_local_velocity = 0;
- double max_velocity = 0.;
+ typename DoFHandler<dim>::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<n_q_points; ++q)
+ max_local_velocity = std::max (max_local_velocity,
+ velocity_values[q].norm());
+ }
+
+ double max_velocity = 0.;
#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
- MPI_Allreduce (&max_local_velocity, &max_velocity, 1, MPI_DOUBLE,
- MPI_MAX, MPI_COMM_WORLD);
+ MPI_Allreduce (&max_local_velocity, &max_velocity, 1, MPI_DOUBLE,
+ MPI_MAX, MPI_COMM_WORLD);
#else
- max_velocity = max_local_velocity;
+ max_velocity = max_local_velocity;
#endif
- return max_velocity;
-}
+ return max_velocity;
+ }
// compute the cfl number, i.e., maximal
// velocity on a cell divided by the cell
// diameter
-template <int dim>
-double BoussinesqFlowProblem<dim>::get_cfl_number () const
-{
- const QIterated<dim> quadrature_formula (QTrapez<1>(),
- parameters.stokes_velocity_degree);
- const unsigned int n_q_points = quadrature_formula.size();
-
- FEValues<dim> fe_values (mapping, stokes_fe, quadrature_formula, update_values);
- std::vector<Tensor<1,dim> > velocity_values(n_q_points);
+ template <int dim>
+ double BoussinesqFlowProblem<dim>::get_cfl_number () const
+ {
+ const QIterated<dim> quadrature_formula (QTrapez<1>(),
+ parameters.stokes_velocity_degree);
+ const unsigned int n_q_points = quadrature_formula.size();
- const FEValuesExtractors::Vector velocities (0);
+ FEValues<dim> fe_values (mapping, stokes_fe, quadrature_formula, update_values);
+ std::vector<Tensor<1,dim> > velocity_values(n_q_points);
- double max_local_cfl = 0;
+ const FEValuesExtractors::Vector velocities (0);
- typename DoFHandler<dim>::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; q<n_q_points; ++q)
- max_local_velocity = std::max (max_local_velocity,
- velocity_values[q].norm());
- max_local_cfl = std::max(max_local_cfl,
- max_local_velocity / cell->diameter());
- }
+ double max_local_cfl = 0;
- double max_cfl_number = 0.;
+ typename DoFHandler<dim>::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; q<n_q_points; ++q)
+ max_local_velocity = std::max (max_local_velocity,
+ velocity_values[q].norm());
+ max_local_cfl = std::max(max_local_cfl,
+ max_local_velocity / cell->diameter());
+ }
+
+ 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 <int dim>
-double
-BoussinesqFlowProblem<dim>::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<dim> quadrature_formula (parameters.temperature_degree+1);
- const unsigned int n_q_points = quadrature_formula.size();
-
- FEValues<dim> fe_values (temperature_fe, quadrature_formula,
- update_values | update_JxW_values);
- std::vector<double> old_temperature_values(n_q_points);
- std::vector<double> old_old_temperature_values(n_q_points);
-
- double min_entropy = std::numeric_limits<double>::max(),
- max_entropy = -std::numeric_limits<double>::max(),
- area = 0,
- entropy_integrated = 0;
-
- typename DoFHandler<dim>::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<n_q_points; ++q)
- {
- const double T = (old_temperature_values[q] +
- old_old_temperature_values[q]) / 2;
- const double entropy = ((T-average_temperature) *
- (T-average_temperature));
-
- min_entropy = std::min (min_entropy, entropy);
- max_entropy = std::max (max_entropy, entropy);
- area += fe_values.JxW(q);
- entropy_integrated += fe_values.JxW(q) * entropy;
- }
- }
+ template <int dim>
+ double
+ BoussinesqFlowProblem<dim>::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<dim> quadrature_formula (parameters.temperature_degree+1);
+ const unsigned int n_q_points = quadrature_formula.size();
+
+ FEValues<dim> fe_values (temperature_fe, quadrature_formula,
+ update_values | update_JxW_values);
+ std::vector<double> old_temperature_values(n_q_points);
+ std::vector<double> old_old_temperature_values(n_q_points);
+
+ double min_entropy = std::numeric_limits<double>::max(),
+ max_entropy = -std::numeric_limits<double>::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<dim>::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<n_q_points; ++q)
+ {
+ const double T = (old_temperature_values[q] +
+ old_old_temperature_values[q]) / 2;
+ const double entropy = ((T-average_temperature) *
+ (T-average_temperature));
+
+ min_entropy = std::min (min_entropy, entropy);
+ max_entropy = std::max (max_entropy, entropy);
+ area += fe_values.JxW(q);
+ entropy_integrated += fe_values.JxW(q) * entropy;
+ }
+ }
+
+ // 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;
#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
- MPI_Allreduce (&local_for_sum[0], &global_for_sum[0], 2, MPI_DOUBLE,
- MPI_SUM, MPI_COMM_WORLD);
- MPI_Allreduce (&local_for_max[0], &global_for_max[0], 2, MPI_DOUBLE,
- MPI_MAX, MPI_COMM_WORLD);
+ MPI_Allreduce (&local_for_sum[0], &global_for_sum[0], 2, MPI_DOUBLE,
+ MPI_SUM, MPI_COMM_WORLD);
+ MPI_Allreduce (&local_for_max[0], &global_for_max[0], 2, MPI_DOUBLE,
+ MPI_MAX, MPI_COMM_WORLD);
#else
- global_for_sum[0] = local_for_sum[0];
- global_for_sum[1] = local_for_sum[1];
- global_for_max[0] = local_for_max[0];
- global_for_max[1] = local_for_max[1];
+ global_for_sum[0] = local_for_sum[0];
+ global_for_sum[1] = local_for_sum[1];
+ global_for_max[0] = local_for_max[0];
+ global_for_max[1] = local_for_max[1];
#endif
- const double average_entropy = global_for_sum[0] / global_for_sum[1];
- const double entropy_diff = std::max(global_for_max[1] - average_entropy,
- average_entropy - (-global_for_max[0]));
- return entropy_diff;
-}
+ const double average_entropy = global_for_sum[0] / global_for_sum[1];
+ const double entropy_diff = std::max(global_for_max[1] - average_entropy,
+ average_entropy - (-global_for_max[0]));
+ return entropy_diff;
+ }
// partition. Two global communication
// steps synchronize the data among the
// processors.
-template <int dim>
-std::pair<double,double>
-BoussinesqFlowProblem<dim>::get_extrapolated_temperature_range () const
-{
- const QIterated<dim> quadrature_formula (QTrapez<1>(),
- parameters.temperature_degree);
- const unsigned int n_q_points = quadrature_formula.size();
-
- FEValues<dim> fe_values (mapping, temperature_fe, quadrature_formula,
- update_values);
- std::vector<double> old_temperature_values(n_q_points);
- std::vector<double> 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<double>::max(),
- max_local_temperature = -std::numeric_limits<double>::max();
-
- if (timestep_number != 0)
- {
- typename DoFHandler<dim>::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<n_q_points; ++q)
- {
- const double temperature =
- (1. + time_step/old_time_step) * old_temperature_values[q]-
- time_step/old_time_step * old_old_temperature_values[q];
-
- min_local_temperature = std::min (min_local_temperature,
- temperature);
- max_local_temperature = std::max (max_local_temperature,
- temperature);
- }
- }
- }
- else
- {
- typename DoFHandler<dim>::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<n_q_points; ++q)
- {
- const double temperature = old_temperature_values[q];
-
- min_local_temperature = std::min (min_local_temperature,
- temperature);
- max_local_temperature = std::max (max_local_temperature,
- temperature);
- }
- }
- }
+ template <int dim>
+ std::pair<double,double>
+ BoussinesqFlowProblem<dim>::get_extrapolated_temperature_range () const
+ {
+ const QIterated<dim> quadrature_formula (QTrapez<1>(),
+ parameters.temperature_degree);
+ const unsigned int n_q_points = quadrature_formula.size();
+
+ FEValues<dim> fe_values (mapping, temperature_fe, quadrature_formula,
+ update_values);
+ std::vector<double> old_temperature_values(n_q_points);
+ std::vector<double> 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<double>::max(),
+ max_local_temperature = -std::numeric_limits<double>::max();
+
+ if (timestep_number != 0)
+ {
+ typename DoFHandler<dim>::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<n_q_points; ++q)
+ {
+ const double temperature =
+ (1. + time_step/old_time_step) * old_temperature_values[q]-
+ time_step/old_time_step * old_old_temperature_values[q];
+
+ min_local_temperature = std::min (min_local_temperature,
+ temperature);
+ max_local_temperature = std::max (max_local_temperature,
+ temperature);
+ }
+ }
+ }
+ else
+ {
+ typename DoFHandler<dim>::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<n_q_points; ++q)
+ {
+ const double temperature = old_temperature_values[q];
+
+ min_local_temperature = std::min (min_local_temperature,
+ temperature);
+ max_local_temperature = std::max (max_local_temperature,
+ temperature);
+ }
+ }
+ }
- double min_temperature, max_temperature;
+ double min_temperature, max_temperature;
#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
- MPI_Allreduce (&max_local_temperature, &max_temperature, 1, MPI_DOUBLE,
- MPI_MAX, MPI_COMM_WORLD);
- MPI_Allreduce (&min_local_temperature, &min_temperature, 1, MPI_DOUBLE,
- MPI_MIN, MPI_COMM_WORLD);
+ MPI_Allreduce (&max_local_temperature, &max_temperature, 1, MPI_DOUBLE,
+ MPI_MAX, MPI_COMM_WORLD);
+ MPI_Allreduce (&min_local_temperature, &min_temperature, 1, MPI_DOUBLE,
+ MPI_MIN, MPI_COMM_WORLD);
#else
- min_temperature = min_local_temperature;
- max_temperature = max_local_temperature;
+ min_temperature = min_local_temperature;
+ max_temperature = max_local_temperature;
#endif
- return std::make_pair(min_temperature, max_temperature);
-}
+ return std::make_pair(min_temperature, max_temperature);
+ }
// The function that calculates the
// viscosity is purely local, so this is
// the same code as in step-31.
-template <int dim>
-double
-BoussinesqFlowProblem<dim>::
-compute_viscosity (const std::vector<double> &old_temperature,
- const std::vector<double> &old_old_temperature,
- const std::vector<Tensor<1,dim> > &old_temperature_grads,
- const std::vector<Tensor<1,dim> > &old_old_temperature_grads,
- const std::vector<double> &old_temperature_laplacians,
- const std::vector<double> &old_old_temperature_laplacians,
- const std::vector<Tensor<1,dim> > &old_velocity_values,
- const std::vector<Tensor<1,dim> > &old_old_velocity_values,
- const std::vector<SymmetricTensor<2,dim> > &old_strain_rates,
- const std::vector<SymmetricTensor<2,dim> > &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 <int dim>
+ double
+ BoussinesqFlowProblem<dim>::
+ compute_viscosity (const std::vector<double> &old_temperature,
+ const std::vector<double> &old_old_temperature,
+ const std::vector<Tensor<1,dim> > &old_temperature_grads,
+ const std::vector<Tensor<1,dim> > &old_old_temperature_grads,
+ const std::vector<double> &old_temperature_laplacians,
+ const std::vector<double> &old_old_temperature_laplacians,
+ const std::vector<Tensor<1,dim> > &old_velocity_values,
+ const std::vector<Tensor<1,dim> > &old_old_velocity_values,
+ const std::vector<SymmetricTensor<2,dim> > &old_strain_rates,
+ const std::vector<SymmetricTensor<2,dim> > &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);
+ }
+ }
// exactly the entry that we would fill the
// matrix with — otherwise, this will
// not work.
-template <int dim>
-void BoussinesqFlowProblem<dim>::project_temperature_field ()
-{
- assemble_temperature_matrix ();
+ template <int dim>
+ void BoussinesqFlowProblem<dim>::project_temperature_field ()
+ {
+ assemble_temperature_matrix ();
- QGauss<dim> quadrature(parameters.temperature_degree+2);
- UpdateFlags update_flags = UpdateFlags(update_values |
- update_quadrature_points |
- update_JxW_values);
- FEValues<dim> fe_values (mapping, temperature_fe, quadrature, update_flags);
+ QGauss<dim> quadrature(parameters.temperature_degree+2);
+ UpdateFlags update_flags = UpdateFlags(update_values |
+ update_quadrature_points |
+ update_JxW_values);
+ FEValues<dim> 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<unsigned int> local_dof_indices (dofs_per_cell);
- Vector<double> cell_vector (dofs_per_cell);
- FullMatrix<double> matrix_for_bc (dofs_per_cell, dofs_per_cell);
+ std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+ Vector<double> cell_vector (dofs_per_cell);
+ FullMatrix<double> matrix_for_bc (dofs_per_cell, dofs_per_cell);
- std::vector<double> rhs_values(n_q_points);
+ std::vector<double> 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<dim>::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<dim>().value_list
- (fe_values.get_quadrature_points(), rhs_values);
-
- cell_vector = 0;
- matrix_for_bc = 0;
- for (unsigned int point=0; point<n_q_points; ++point)
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- cell_vector(i) += rhs_values[point] *
- fe_values.shape_value(i,point) *
- fe_values.JxW(point);
- if (temperature_constraints.is_inhomogeneously_constrained(local_dof_indices[i]))
- {
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- matrix_for_bc(j,i) += fe_values.shape_value(i,point) *
- fe_values.shape_value(j,point) *
- fe_values.JxW(point);
- }
- }
-
- temperature_constraints.distribute_local_to_global (cell_vector,
- local_dof_indices,
- rhs,
- matrix_for_bc);
- }
+ typename DoFHandler<dim>::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<TrilinosWrappers::MPI::Vector> 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
- // <code>old_temperature_solution</code>
- // 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
- // <code>temperature_solution</code>
- // and
- // <code>old_old_temperature_solution</code>
- // 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
- // <code>solution</code> 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<dim>().value_list
+ (fe_values.get_quadrature_points(), rhs_values);
+
+ cell_vector = 0;
+ matrix_for_bc = 0;
+ for (unsigned int point=0; point<n_q_points; ++point)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ cell_vector(i) += rhs_values[point] *
+ fe_values.shape_value(i,point) *
+ fe_values.JxW(point);
+ if (temperature_constraints.is_inhomogeneously_constrained(local_dof_indices[i]))
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ matrix_for_bc(j,i) += fe_values.shape_value(i,point) *
+ fe_values.shape_value(j,point) *
+ fe_values.JxW(point);
+ }
+ }
+
+ temperature_constraints.distribute_local_to_global (cell_vector,
+ local_dof_indices,
+ rhs,
+ matrix_for_bc);
+ }
+
+ rhs.compress (Add);
+
+ SolverControl solver_control(5*rhs.size(), 1e-12*rhs.l2_norm());
+ SolverCG<TrilinosWrappers::MPI::Vector> 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
+ // <code>old_temperature_solution</code>
+ // 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
+ // <code>temperature_solution</code>
+ // and
+ // <code>old_old_temperature_solution</code>
+ // 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
+ // <code>solution</code> 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);
+ }
// patterns, so we can safely release the
// <code>sp</code> variable once the matrix
// has been given the sparsity structure.
-template <int dim>
-void BoussinesqFlowProblem<dim>::
-setup_stokes_matrix (const std::vector<IndexSet> &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<dim+1; ++c)
- for (unsigned int d=0; d<dim+1; ++d)
- if (! ((c==dim) && (d==dim)))
- coupling[c][d] = DoFTools::always;
- else
- coupling[c][d] = DoFTools::none;
+ template <int dim>
+ void BoussinesqFlowProblem<dim>::
+ setup_stokes_matrix (const std::vector<IndexSet> &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<dim+1; ++c)
+ for (unsigned int d=0; d<dim+1; ++d)
+ if (! ((c==dim) && (d==dim)))
+ coupling[c][d] = DoFTools::always;
+ else
+ coupling[c][d] = DoFTools::none;
+ DoFTools::make_sparsity_pattern (stokes_dof_handler,
+ coupling, sp,
+ stokes_constraints, false,
+ Utilities::System::
+ get_this_mpi_process(MPI_COMM_WORLD));
+ sp.compress();
-template <int dim>
-void BoussinesqFlowProblem<dim>::
-setup_stokes_preconditioner (const std::vector<IndexSet> &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<dim+1; ++c)
- for (unsigned int d=0; d<dim+1; ++d)
- if (c == d)
- coupling[c][d] = DoFTools::always;
- else
- coupling[c][d] = DoFTools::none;
- DoFTools::make_sparsity_pattern (stokes_dof_handler,
- coupling, sp,
- stokes_constraints, false,
- Utilities::System::
- get_this_mpi_process(MPI_COMM_WORLD));
- sp.compress();
- stokes_preconditioner_matrix.reinit (sp);
-}
+ template <int dim>
+ void BoussinesqFlowProblem<dim>::
+ setup_stokes_preconditioner (const std::vector<IndexSet> &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<dim+1; ++c)
+ for (unsigned int d=0; d<dim+1; ++d)
+ if (c == d)
+ coupling[c][d] = DoFTools::always;
+ else
+ coupling[c][d] = DoFTools::none;
+
+ DoFTools::make_sparsity_pattern (stokes_dof_handler,
+ coupling, sp,
+ stokes_constraints, false,
+ Utilities::System::
+ get_this_mpi_process(MPI_COMM_WORLD));
+ sp.compress();
+
+ stokes_preconditioner_matrix.reinit (sp);
+ }
-template <int dim>
-void BoussinesqFlowProblem<dim>::
-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 <int dim>
+ void BoussinesqFlowProblem<dim>::
+ 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);
+ }
// get_this_mpi_process
// (MPI_COMM_WORLD) == 0)</code>,
// hardly a pretty solution.
-template <int dim>
-void BoussinesqFlowProblem<dim>::setup_dofs ()
-{
- computing_timer.enter_section("Setup dof systems");
-
- std::vector<unsigned int> 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<unsigned int> 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<IndexSet> stokes_partitioning, stokes_relevant_partitioning;
- IndexSet temperature_partitioning (n_T), temperature_relevant_partitioning (n_T);
- IndexSet stokes_relevant_set;
+ template <int dim>
+ void BoussinesqFlowProblem<dim>::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<unsigned int> 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<unsigned int> 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<IndexSet> 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<bool> velocity_mask (dim+1, true);
- velocity_mask[dim] = false;
- VectorTools::interpolate_boundary_values (stokes_dof_handler,
- 0,
- ZeroFunction<dim>(dim+1),
- stokes_constraints,
- velocity_mask);
-
- std::set<unsigned char> 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<dim>(),
- temperature_constraints);
- VectorTools::interpolate_boundary_values (temperature_dof_handler,
- 1,
- EquationData::TemperatureInitialValues<dim>(),
- temperature_constraints);
- temperature_constraints.close ();
- }
+ stokes_constraints.reinit(stokes_relevant_set);
+
+ DoFTools::make_hanging_node_constraints (stokes_dof_handler,
+ stokes_constraints);
+
+ std::vector<bool> velocity_mask (dim+1, true);
+ velocity_mask[dim] = false;
+ VectorTools::interpolate_boundary_values (stokes_dof_handler,
+ 0,
+ ZeroFunction<dim>(dim+1),
+ stokes_constraints,
+ velocity_mask);
+
+ std::set<unsigned char> 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<dim>(),
+ temperature_constraints);
+ VectorTools::interpolate_boundary_values (temperature_dof_handler,
+ 1,
+ EquationData::TemperatureInitialValues<dim>(),
+ 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();
+ }
// scratch data object, which is allocated
// only once before we start the loop over
// all cells:
-template <int dim>
-void
-BoussinesqFlowProblem<dim>::
-local_assemble_stokes_preconditioner (const typename DoFHandler<dim>::active_cell_iterator &cell,
- Assembly::Scratch::StokesPreconditioner<dim> &scratch,
- Assembly::CopyData::StokesPreconditioner<dim> &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 <int dim>
+ void
+ BoussinesqFlowProblem<dim>::
+ local_assemble_stokes_preconditioner (const typename DoFHandler<dim>::active_cell_iterator &cell,
+ Assembly::Scratch::StokesPreconditioner<dim> &scratch,
+ Assembly::CopyData::StokesPreconditioner<dim> &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<n_q_points; ++q)
- {
- for (unsigned int k=0; k<dofs_per_cell; ++k)
- {
- scratch.grads_phi_u[k] = scratch.stokes_fe_values[velocities].symmetric_gradient(k,q);
- scratch.phi_p[k] = scratch.stokes_fe_values[pressure].value (k, q);
- }
-
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- if (stokes_fe.system_to_component_index(i).first
- ==
- stokes_fe.system_to_component_index(j).first)
- data.local_matrix(i,j) += (EquationData::eta *
- (scratch.grads_phi_u[i] *
- scratch.grads_phi_u[j])
- +
- (1./EquationData::eta) *
- EquationData::pressure_scaling *
- EquationData::pressure_scaling *
- (scratch.phi_p[i] * scratch.phi_p[j]))
- * scratch.stokes_fe_values.JxW(q);
- }
-}
+ for (unsigned int q=0; q<n_q_points; ++q)
+ {
+ for (unsigned int k=0; k<dofs_per_cell; ++k)
+ {
+ scratch.grads_phi_u[k] = scratch.stokes_fe_values[velocities].symmetric_gradient(k,q);
+ scratch.phi_p[k] = scratch.stokes_fe_values[pressure].value (k, q);
+ }
+
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ if (stokes_fe.system_to_component_index(i).first
+ ==
+ stokes_fe.system_to_component_index(j).first)
+ data.local_matrix(i,j) += (EquationData::eta *
+ (scratch.grads_phi_u[i] *
+ scratch.grads_phi_u[j])
+ +
+ (1./EquationData::eta) *
+ EquationData::pressure_scaling *
+ EquationData::pressure_scaling *
+ (scratch.phi_p[i] * scratch.phi_p[j]))
+ * scratch.stokes_fe_values.JxW(q);
+ }
+ }
-template <int dim>
-void
-BoussinesqFlowProblem<dim>::
-copy_local_to_global_stokes_preconditioner (const Assembly::CopyData::StokesPreconditioner<dim> &data)
-{
- stokes_constraints.distribute_local_to_global (data.local_matrix,
- data.local_dof_indices,
- stokes_preconditioner_matrix);
-}
+ template <int dim>
+ void
+ BoussinesqFlowProblem<dim>::
+ copy_local_to_global_stokes_preconditioner (const Assembly::CopyData::StokesPreconditioner<dim> &data)
+ {
+ stokes_constraints.distribute_local_to_global (data.local_matrix,
+ data.local_dof_indices,
+ stokes_preconditioner_matrix);
+ }
// send it to the owner at the end of
// assembly, by calling the
// <code>compress()</code> command.
-template <int dim>
-void
-BoussinesqFlowProblem<dim>::assemble_stokes_preconditioner ()
-{
- stokes_preconditioner_matrix = 0;
-
- const QGauss<dim> quadrature_formula(parameters.stokes_velocity_degree+1);
-
- typedef
- FilteredIterator<typename DoFHandler<dim>::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<dim>::
- local_assemble_stokes_preconditioner,
- this,
- std_cxx1x::_1,
- std_cxx1x::_2,
- std_cxx1x::_3),
- std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
- copy_local_to_global_stokes_preconditioner,
- this,
- std_cxx1x::_1),
- Assembly::Scratch::
- StokesPreconditioner<dim> (stokes_fe, quadrature_formula,
- mapping,
- update_JxW_values |
- update_values |
- update_gradients),
- Assembly::CopyData::
- StokesPreconditioner<dim> (stokes_fe));
-
- stokes_preconditioner_matrix.compress();
-}
+ template <int dim>
+ void
+ BoussinesqFlowProblem<dim>::assemble_stokes_preconditioner ()
+ {
+ stokes_preconditioner_matrix = 0;
+
+ const QGauss<dim> quadrature_formula(parameters.stokes_velocity_degree+1);
+
+ typedef
+ FilteredIterator<typename DoFHandler<dim>::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<dim>::
+ local_assemble_stokes_preconditioner,
+ this,
+ std_cxx1x::_1,
+ std_cxx1x::_2,
+ std_cxx1x::_3),
+ std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
+ copy_local_to_global_stokes_preconditioner,
+ this,
+ std_cxx1x::_1),
+ Assembly::Scratch::
+ StokesPreconditioner<dim> (stokes_fe, quadrature_formula,
+ mapping,
+ update_JxW_values |
+ update_values |
+ update_gradients),
+ Assembly::CopyData::
+ StokesPreconditioner<dim> (stokes_fe));
+
+ stokes_preconditioner_matrix.compress();
+ }
// preconditioner for the pressure mass
// matrix instead of IC, as discussed in the
// introduction.
-template <int dim>
-void
-BoussinesqFlowProblem<dim>::build_stokes_preconditioner ()
-{
- if (rebuild_stokes_preconditioner == false)
- return;
+ template <int dim>
+ void
+ BoussinesqFlowProblem<dim>::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<std::vector<bool> > constant_modes;
- std::vector<bool> velocity_components (dim+1,true);
- velocity_components[dim] = false;
- DoFTools::extract_constant_modes (stokes_dof_handler, velocity_components,
- constant_modes);
+ std::vector<std::vector<bool> > constant_modes;
+ std::vector<bool> 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<TrilinosWrappers::PreconditionILU>
- (new TrilinosWrappers::PreconditionILU());
- Amg_preconditioner = std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionAMG>
- (new TrilinosWrappers::PreconditionAMG());
+ Mp_preconditioner = std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionILU>
+ (new TrilinosWrappers::PreconditionILU());
+ Amg_preconditioner = std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionAMG>
+ (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}
// functions in the end of the assembly in
// order to send non-local data to the
// owner process.
-template <int dim>
-void
-BoussinesqFlowProblem<dim>::
-local_assemble_stokes_system (const typename DoFHandler<dim>::active_cell_iterator &cell,
- Assembly::Scratch::StokesSystem<dim> &scratch,
- Assembly::CopyData::StokesSystem<dim> &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 <int dim>
+ void
+ BoussinesqFlowProblem<dim>::
+ local_assemble_stokes_system (const typename DoFHandler<dim>::active_cell_iterator &cell,
+ Assembly::Scratch::StokesSystem<dim> &scratch,
+ Assembly::CopyData::StokesSystem<dim> &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<dim>::active_cell_iterator
- temperature_cell (&triangulation,
- cell->level(),
- cell->index(),
- &temperature_dof_handler);
- scratch.temperature_fe_values.reinit (temperature_cell);
+ typename DoFHandler<dim>::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<n_q_points; ++q)
- {
- const double old_temperature = scratch.old_temperature_values[q];
-
- for (unsigned int k=0; k<dofs_per_cell; ++k)
- {
- scratch.phi_u[k] = scratch.stokes_fe_values[velocities].value (k,q);
- if (rebuild_stokes_matrix)
- {
- scratch.grads_phi_u[k] = scratch.stokes_fe_values[velocities].symmetric_gradient(k,q);
- scratch.div_phi_u[k] = scratch.stokes_fe_values[velocities].divergence (k, q);
- scratch.phi_p[k] = scratch.stokes_fe_values[pressure].value (k, q);
- }
- }
-
- if (rebuild_stokes_matrix)
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- data.local_matrix(i,j) += (EquationData::eta * 2 *
- (scratch.grads_phi_u[i] * scratch.grads_phi_u[j])
- - (EquationData::pressure_scaling *
- scratch.div_phi_u[i] * scratch.phi_p[j])
- - (EquationData::pressure_scaling *
- scratch.phi_p[i] * scratch.div_phi_u[j]))
- * scratch.stokes_fe_values.JxW(q);
-
- const Tensor<1,dim>
- gravity = EquationData::gravity_vector (scratch.stokes_fe_values
- .quadrature_point(q));
-
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- data.local_rhs(i) += (EquationData::density(old_temperature) *
- gravity *
- scratch.phi_u[i]) *
- scratch.stokes_fe_values.JxW(q);
- }
+ for (unsigned int q=0; q<n_q_points; ++q)
+ {
+ const double old_temperature = scratch.old_temperature_values[q];
+
+ for (unsigned int k=0; k<dofs_per_cell; ++k)
+ {
+ scratch.phi_u[k] = scratch.stokes_fe_values[velocities].value (k,q);
+ if (rebuild_stokes_matrix)
+ {
+ scratch.grads_phi_u[k] = scratch.stokes_fe_values[velocities].symmetric_gradient(k,q);
+ scratch.div_phi_u[k] = scratch.stokes_fe_values[velocities].divergence (k, q);
+ scratch.phi_p[k] = scratch.stokes_fe_values[pressure].value (k, q);
+ }
+ }
+
+ if (rebuild_stokes_matrix)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ data.local_matrix(i,j) += (EquationData::eta * 2 *
+ (scratch.grads_phi_u[i] * scratch.grads_phi_u[j])
+ - (EquationData::pressure_scaling *
+ scratch.div_phi_u[i] * scratch.phi_p[j])
+ - (EquationData::pressure_scaling *
+ scratch.phi_p[i] * scratch.div_phi_u[j]))
+ * scratch.stokes_fe_values.JxW(q);
+
+ const Tensor<1,dim>
+ gravity = EquationData::gravity_vector (scratch.stokes_fe_values
+ .quadrature_point(q));
+
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ data.local_rhs(i) += (EquationData::density(old_temperature) *
+ gravity *
+ scratch.phi_u[i]) *
+ scratch.stokes_fe_values.JxW(q);
+ }
- cell->get_dof_indices (data.local_dof_indices);
-}
+ cell->get_dof_indices (data.local_dof_indices);
+ }
-template <int dim>
-void
-BoussinesqFlowProblem<dim>::
-copy_local_to_global_stokes_system (const Assembly::CopyData::StokesSystem<dim> &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 <int dim>
+ void
+ BoussinesqFlowProblem<dim>::
+ copy_local_to_global_stokes_system (const Assembly::CopyData::StokesSystem<dim> &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 <int dim>
-void BoussinesqFlowProblem<dim>::assemble_stokes_system ()
-{
- computing_timer.enter_section (" Assemble Stokes system");
-
- if (rebuild_stokes_matrix == true)
- stokes_matrix=0;
-
- stokes_rhs=0;
-
- const QGauss<dim> quadrature_formula(parameters.stokes_velocity_degree+1);
-
- typedef
- FilteredIterator<typename DoFHandler<dim>::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<dim>::
- local_assemble_stokes_system,
- this,
- std_cxx1x::_1,
- std_cxx1x::_2,
- std_cxx1x::_3),
- std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
- copy_local_to_global_stokes_system,
- this,
- std_cxx1x::_1),
- Assembly::Scratch::
- StokesSystem<dim> (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<dim> (stokes_fe));
-
- stokes_matrix.compress();
- stokes_rhs.compress(Add);
-
- rebuild_stokes_matrix = false;
-
- pcout << std::endl;
- computing_timer.exit_section();
-}
+ template <int dim>
+ void BoussinesqFlowProblem<dim>::assemble_stokes_system ()
+ {
+ computing_timer.enter_section (" Assemble Stokes system");
+
+ if (rebuild_stokes_matrix == true)
+ stokes_matrix=0;
+
+ stokes_rhs=0;
+
+ const QGauss<dim> quadrature_formula(parameters.stokes_velocity_degree+1);
+
+ typedef
+ FilteredIterator<typename DoFHandler<dim>::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<dim>::
+ local_assemble_stokes_system,
+ this,
+ std_cxx1x::_1,
+ std_cxx1x::_2,
+ std_cxx1x::_3),
+ std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
+ copy_local_to_global_stokes_system,
+ this,
+ std_cxx1x::_1),
+ Assembly::Scratch::
+ StokesSystem<dim> (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<dim> (stokes_fe));
+
+ stokes_matrix.compress();
+ stokes_rhs.compress(Add);
+
+ rebuild_stokes_matrix = false;
+
+ pcout << std::endl;
+ computing_timer.exit_section();
+ }
// @sect5{Temperature matrix assembly}
//
// The two following functions perform
// similar services as the ones above.
-template <int dim>
-void BoussinesqFlowProblem<dim>::
-local_assemble_temperature_matrix (const typename DoFHandler<dim>::active_cell_iterator &cell,
- Assembly::Scratch::TemperatureMatrix<dim> &scratch,
- Assembly::CopyData::TemperatureMatrix<dim> &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 <int dim>
+ void BoussinesqFlowProblem<dim>::
+ local_assemble_temperature_matrix (const typename DoFHandler<dim>::active_cell_iterator &cell,
+ Assembly::Scratch::TemperatureMatrix<dim> &scratch,
+ Assembly::CopyData::TemperatureMatrix<dim> &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<n_q_points; ++q)
- {
- for (unsigned int k=0; k<dofs_per_cell; ++k)
- {
- scratch.grad_phi_T[k] = scratch.temperature_fe_values.shape_grad (k,q);
- scratch.phi_T[k] = scratch.temperature_fe_values.shape_value (k, q);
- }
-
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- {
- data.local_mass_matrix(i,j)
- += (scratch.phi_T[i] * scratch.phi_T[j]
- *
- scratch.temperature_fe_values.JxW(q));
- data.local_stiffness_matrix(i,j)
- += (EquationData::kappa * scratch.grad_phi_T[i] * scratch.grad_phi_T[j]
- *
- scratch.temperature_fe_values.JxW(q));
- }
- }
-}
+ for (unsigned int q=0; q<n_q_points; ++q)
+ {
+ for (unsigned int k=0; k<dofs_per_cell; ++k)
+ {
+ scratch.grad_phi_T[k] = scratch.temperature_fe_values.shape_grad (k,q);
+ scratch.phi_T[k] = scratch.temperature_fe_values.shape_value (k, q);
+ }
+
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ data.local_mass_matrix(i,j)
+ += (scratch.phi_T[i] * scratch.phi_T[j]
+ *
+ scratch.temperature_fe_values.JxW(q));
+ data.local_stiffness_matrix(i,j)
+ += (EquationData::kappa * scratch.grad_phi_T[i] * scratch.grad_phi_T[j]
+ *
+ scratch.temperature_fe_values.JxW(q));
+ }
+ }
+ }
-template <int dim>
-void
-BoussinesqFlowProblem<dim>::
-copy_local_to_global_temperature_matrix (const Assembly::CopyData::TemperatureMatrix<dim> &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 <int dim>
+ void
+ BoussinesqFlowProblem<dim>::
+ copy_local_to_global_temperature_matrix (const Assembly::CopyData::TemperatureMatrix<dim> &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 <int dim>
-void BoussinesqFlowProblem<dim>::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<dim> quadrature_formula(parameters.temperature_degree+2);
-
- typedef
- FilteredIterator<typename DoFHandler<dim>::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<dim>::
- local_assemble_temperature_matrix,
- this,
- std_cxx1x::_1,
- std_cxx1x::_2,
- std_cxx1x::_3),
- std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
- copy_local_to_global_temperature_matrix,
- this,
- std_cxx1x::_1),
- Assembly::Scratch::
- TemperatureMatrix<dim> (temperature_fe, mapping, quadrature_formula),
- Assembly::CopyData::
- TemperatureMatrix<dim> (temperature_fe));
-
- temperature_mass_matrix.compress();
- temperature_stiffness_matrix.compress();
-
- rebuild_temperature_matrices = false;
- rebuild_temperature_preconditioner = true;
-
- computing_timer.exit_section();
-}
+ template <int dim>
+ void BoussinesqFlowProblem<dim>::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<dim> quadrature_formula(parameters.temperature_degree+2);
+
+ typedef
+ FilteredIterator<typename DoFHandler<dim>::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<dim>::
+ local_assemble_temperature_matrix,
+ this,
+ std_cxx1x::_1,
+ std_cxx1x::_2,
+ std_cxx1x::_3),
+ std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
+ copy_local_to_global_temperature_matrix,
+ this,
+ std_cxx1x::_1),
+ Assembly::Scratch::
+ TemperatureMatrix<dim> (temperature_fe, mapping, quadrature_formula),
+ Assembly::CopyData::
+ TemperatureMatrix<dim> (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}
// dofs. That will account for the correct
// balance of the right hand side vector
// with the matrix system of temperature.
-template <int dim>
-void BoussinesqFlowProblem<dim>::
-local_assemble_temperature_rhs (const std::pair<double,double> global_T_range,
- const double global_max_velocity,
- const double global_entropy_variation,
- const typename DoFHandler<dim>::active_cell_iterator &cell,
- Assembly::Scratch::TemperatureRHS<dim> &scratch,
- Assembly::CopyData::TemperatureRHS<dim> &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<dim>::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<n_q_points; ++q)
- {
- for (unsigned int k=0; k<dofs_per_cell; ++k)
- {
- scratch.phi_T[k] = scratch.temperature_fe_values.shape_value (k, q);
- scratch.grad_phi_T[k] = scratch.temperature_fe_values.shape_grad (k,q);
- }
-
-
- const double old_Ts
- = (use_bdf2_scheme ?
- (scratch.old_temperature_values[q] *
- (time_step + old_time_step) / old_time_step
- -
- scratch.old_old_temperature_values[q] *
- (time_step * time_step) /
- (old_time_step * (time_step + old_time_step)))
- :
- scratch.old_temperature_values[q]);
-
- const Tensor<1,dim> 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<dofs_per_cell; ++i)
- {
- data.local_rhs(i) += (old_Ts * scratch.phi_T[i]
- -
- time_step *
- extrapolated_u * ext_grad_T * scratch.phi_T[i]
- -
- time_step *
- nu * ext_grad_T * scratch.grad_phi_T[i]
- +
- time_step *
- gamma * scratch.phi_T[i])
- *
- scratch.temperature_fe_values.JxW(q);
-
- if (temperature_constraints.is_inhomogeneously_constrained(data.local_dof_indices[i]))
- {
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- data.matrix_for_bc(j,i) += (scratch.phi_T[i] * scratch.phi_T[j] *
- (use_bdf2_scheme ?
- ((2*time_step + old_time_step) /
- (time_step + old_time_step)) : 1.)
- +
- scratch.grad_phi_T[i] *
- scratch.grad_phi_T[j] *
- EquationData::kappa *
- time_step)
- *
- scratch.temperature_fe_values.JxW(q);
- }
- }
- }
-}
+ template <int dim>
+ void BoussinesqFlowProblem<dim>::
+ local_assemble_temperature_rhs (const std::pair<double,double> global_T_range,
+ const double global_max_velocity,
+ const double global_entropy_variation,
+ const typename DoFHandler<dim>::active_cell_iterator &cell,
+ Assembly::Scratch::TemperatureRHS<dim> &scratch,
+ Assembly::CopyData::TemperatureRHS<dim> &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 <int dim>
-void
-BoussinesqFlowProblem<dim>::
-copy_local_to_global_temperature_rhs (const Assembly::CopyData::TemperatureRHS<dim> &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<dim>::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<n_q_points; ++q)
+ {
+ for (unsigned int k=0; k<dofs_per_cell; ++k)
+ {
+ scratch.phi_T[k] = scratch.temperature_fe_values.shape_value (k, q);
+ scratch.grad_phi_T[k] = scratch.temperature_fe_values.shape_grad (k,q);
+ }
+
+
+ const double old_Ts
+ = (use_bdf2_scheme ?
+ (scratch.old_temperature_values[q] *
+ (time_step + old_time_step) / old_time_step
+ -
+ scratch.old_old_temperature_values[q] *
+ (time_step * time_step) /
+ (old_time_step * (time_step + old_time_step)))
+ :
+ scratch.old_temperature_values[q]);
+
+ const Tensor<1,dim> 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<dofs_per_cell; ++i)
+ {
+ data.local_rhs(i) += (old_Ts * scratch.phi_T[i]
+ -
+ time_step *
+ extrapolated_u * ext_grad_T * scratch.phi_T[i]
+ -
+ time_step *
+ nu * ext_grad_T * scratch.grad_phi_T[i]
+ +
+ time_step *
+ gamma * scratch.phi_T[i])
+ *
+ scratch.temperature_fe_values.JxW(q);
+
+ if (temperature_constraints.is_inhomogeneously_constrained(data.local_dof_indices[i]))
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ data.matrix_for_bc(j,i) += (scratch.phi_T[i] * scratch.phi_T[j] *
+ (use_bdf2_scheme ?
+ ((2*time_step + old_time_step) /
+ (time_step + old_time_step)) : 1.)
+ +
+ scratch.grad_phi_T[i] *
+ scratch.grad_phi_T[j] *
+ EquationData::kappa *
+ time_step)
+ *
+ scratch.temperature_fe_values.JxW(q);
+ }
+ }
+ }
+ }
+
+
+ template <int dim>
+ void
+ BoussinesqFlowProblem<dim>::
+ copy_local_to_global_temperature_rhs (const Assembly::CopyData::TemperatureRHS<dim> &data)
+ {
+ temperature_constraints.distribute_local_to_global (data.local_rhs,
+ data.local_dof_indices,
+ temperature_rhs,
+ data.matrix_for_bc);
+ }
// too big a problem because we remesh every
// fifth time step (and regenerate the
// preconditioner then).
-template <int dim>
-void BoussinesqFlowProblem<dim>::assemble_temperature_system (const double maximal_velocity)
-{
- const bool use_bdf2_scheme = (timestep_number != 0);
+ template <int dim>
+ void BoussinesqFlowProblem<dim>::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<TrilinosWrappers::PreconditionIC>
- (new TrilinosWrappers::PreconditionIC());
- T_preconditioner->initialize (temperature_matrix);
- rebuild_temperature_preconditioner = false;
- }
+ if (rebuild_temperature_preconditioner == true)
+ {
+ T_preconditioner = std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionIC>
+ (new TrilinosWrappers::PreconditionIC());
+ T_preconditioner->initialize (temperature_matrix);
+ rebuild_temperature_preconditioner = false;
+ }
- temperature_rhs = 0;
-
- const QGauss<dim> quadrature_formula(parameters.temperature_degree+2);
- const std::pair<double,double>
- 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<typename DoFHandler<dim>::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<dim>::
- 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<dim>::
- copy_local_to_global_temperature_rhs,
- this,
- std_cxx1x::_1),
- Assembly::Scratch::
- TemperatureRHS<dim> (temperature_fe, stokes_fe, mapping,
- quadrature_formula),
- Assembly::CopyData::
- TemperatureRHS<dim> (temperature_fe));
-
- temperature_rhs.compress(Add);
-}
+ temperature_rhs = 0;
+
+ const QGauss<dim> quadrature_formula(parameters.temperature_degree+2);
+ const std::pair<double,double>
+ 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<typename DoFHandler<dim>::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<dim>::
+ 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<dim>::
+ copy_local_to_global_temperature_rhs,
+ this,
+ std_cxx1x::_1),
+ Assembly::Scratch::
+ TemperatureRHS<dim> (temperature_fe, stokes_fe, mapping,
+ quadrature_formula),
+ Assembly::CopyData::
+ TemperatureRHS<dim> (temperature_fe));
+
+ temperature_rhs.compress(Add);
+ }
// Apart from these two changes, everything
// is the same as in step-31, so we don't
// need to further comment on it.
-template <int dim>
-void BoussinesqFlowProblem<dim>::solve ()
-{
- computing_timer.enter_section (" Solve Stokes system");
-
+ template <int dim>
+ void BoussinesqFlowProblem<dim>::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<end; ++i)
- if (stokes_constraints.is_constrained (i))
- distributed_stokes_solution(i) = 0;
-
-
- PrimitiveVectorMemory< TrilinosWrappers::MPI::BlockVector > 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<TrilinosWrappers::PreconditionAMG,
- TrilinosWrappers::PreconditionILU>
- preconditioner (stokes_matrix, stokes_preconditioner_matrix,
- *Mp_preconditioner, *Amg_preconditioner,
- false);
-
- SolverFGMRES<TrilinosWrappers::MPI::BlockVector>
- solver(solver_control, mem,
- SolverFGMRES<TrilinosWrappers::MPI::BlockVector>::
- 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<TrilinosWrappers::PreconditionAMG,
- TrilinosWrappers::PreconditionILU>
- preconditioner (stokes_matrix, stokes_preconditioner_matrix,
- *Mp_preconditioner, *Amg_preconditioner,
- true);
-
- SolverControl solver_control_refined (stokes_matrix.m(), solver_tolerance);
- SolverFGMRES<TrilinosWrappers::MPI::BlockVector>
- solver(solver_control_refined, mem,
- SolverFGMRES<TrilinosWrappers::MPI::BlockVector>::
- 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<end; ++i)
+ if (stokes_constraints.is_constrained (i))
+ distributed_stokes_solution(i) = 0;
+
+
+ PrimitiveVectorMemory< TrilinosWrappers::MPI::BlockVector > 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<TrilinosWrappers::PreconditionAMG,
+ TrilinosWrappers::PreconditionILU>
+ preconditioner (stokes_matrix, stokes_preconditioner_matrix,
+ *Mp_preconditioner, *Amg_preconditioner,
+ false);
+
+ SolverFGMRES<TrilinosWrappers::MPI::BlockVector>
+ solver(solver_control, mem,
+ SolverFGMRES<TrilinosWrappers::MPI::BlockVector>::
+ 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<TrilinosWrappers::PreconditionAMG,
+ TrilinosWrappers::PreconditionILU>
+ preconditioner (stokes_matrix, stokes_preconditioner_matrix,
+ *Mp_preconditioner, *Amg_preconditioner,
+ true);
+
+ SolverControl solver_control_refined (stokes_matrix.m(), solver_tolerance);
+ SolverFGMRES<TrilinosWrappers::MPI::BlockVector>
+ solver(solver_control_refined, mem,
+ SolverFGMRES<TrilinosWrappers::MPI::BlockVector>::
+ 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<TrilinosWrappers::MPI::Vector> cg (solver_control);
+ computing_timer.enter_section (" Solve temperature system");
+ {
+ SolverControl solver_control (temperature_matrix.m(),
+ 1e-12*temperature_rhs.l2_norm());
+ SolverCG<TrilinosWrappers::MPI::Vector> 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<double> temperature (2), global_temperature (2);
- temperature[0] = std::numeric_limits<double>::max(),
- temperature[1] = -std::numeric_limits<double>::max();
- for (unsigned int i=0; i<distributed_temperature_solution.local_size(); ++i)
- {
- temperature[0] = std::min<double> (temperature[0],
- distributed_temperature_solution.trilinos_vector()[0][i]);
- temperature[1] = std::max<double> (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<double> temperature (2), global_temperature (2);
+ temperature[0] = std::numeric_limits<double>::max(),
+ temperature[1] = -std::numeric_limits<double>::max();
+ for (unsigned int i=0; i<distributed_temperature_solution.local_size(); ++i)
+ {
+ temperature[0] = std::min<double> (temperature[0],
+ distributed_temperature_solution.trilinos_vector()[0][i]);
+ temperature[1] = std::max<double> (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 <int dim>
-class BoussinesqFlowProblem<dim>::Postprocessor : public DataPostprocessor<dim>
-{
- public:
- Postprocessor (const unsigned int partition,
- const double minimal_pressure);
+ template <int dim>
+ class BoussinesqFlowProblem<dim>::Postprocessor : public DataPostprocessor<dim>
+ {
+ public:
+ Postprocessor (const unsigned int partition,
+ const double minimal_pressure);
- virtual
- void
- compute_derived_quantities_vector (const std::vector<Vector<double> > &uh,
- const std::vector<std::vector<Tensor<1,dim> > > &duh,
- const std::vector<std::vector<Tensor<2,dim> > > &dduh,
- const std::vector<Point<dim> > &normals,
- const std::vector<Point<dim> > &evaluation_points,
- std::vector<Vector<double> > &computed_quantities) const;
+ virtual
+ void
+ compute_derived_quantities_vector (const std::vector<Vector<double> > &uh,
+ const std::vector<std::vector<Tensor<1,dim> > > &duh,
+ const std::vector<std::vector<Tensor<2,dim> > > &dduh,
+
const std::vector<Point<dim> > &normals,
+
const std::vector<Point<dim> > &evaluation_points,
+ std::vector<Vector<double> > &computed_quantities) const;
- virtual std::vector<std::string> get_names () const;
+ virtual std::vector<std::string> get_names () const;
- virtual unsigned int n_output_variables() const;
+ virtual unsigned int n_output_variables() const;
- virtual
- std::vector<DataComponentInterpretation::DataComponentInterpretation>
- get_data_component_interpretation () const;
+ virtual
+ std::vector<DataComponentInterpretation::DataComponentInterpretation>
+ 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 <int dim>
-BoussinesqFlowProblem<dim>::Postprocessor::
-Postprocessor (const unsigned int partition,
- const double minimal_pressure)
- :
- partition (partition),
- minimal_pressure (minimal_pressure)
-{}
+ template <int dim>
+ BoussinesqFlowProblem<dim>::Postprocessor::
+ Postprocessor (const unsigned int partition,
+ const double minimal_pressure)
+ :
+ partition (partition),
+ minimal_pressure (minimal_pressure)
+ {}
-template <int dim>
-std::vector<std::string>
-BoussinesqFlowProblem<dim>::Postprocessor::get_names() const
-{
- std::vector<std::string> 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 <int dim>
+ std::vector<std::string>
+ BoussinesqFlowProblem<dim>::Postprocessor::get_names() const
+ {
+ std::vector<std::string> 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 <int dim>
-unsigned int
-BoussinesqFlowProblem<dim>::Postprocessor::n_output_variables() const
-{
- return get_names().size();
-}
+ template <int dim>
+ unsigned int
+ BoussinesqFlowProblem<dim>::Postprocessor::n_output_variables() const
+ {
+ return get_names().size();
+ }
-template <int dim>
-std::vector<DataComponentInterpretation::DataComponentInterpretation>
-BoussinesqFlowProblem<dim>::Postprocessor::
-get_data_component_interpretation () const
-{
+ template <int dim>
std::vector<DataComponentInterpretation::DataComponentInterpretation>
- interpretation (dim,
- DataComponentInterpretation::component_is_part_of_vector);
+ BoussinesqFlowProblem<dim>::Postprocessor::
+ get_data_component_interpretation () const
+ {
+ std::vector<DataComponentInterpretation::DataComponentInterpretation>
+ 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 <int dim>
-UpdateFlags
-BoussinesqFlowProblem<dim>::Postprocessor::get_needed_update_flags() const
-{
- return update_values | update_gradients | update_q_points;
-}
+ template <int dim>
+ UpdateFlags
+ BoussinesqFlowProblem<dim>::Postprocessor::get_needed_update_flags() const
+ {
+ return update_values | update_gradients | update_q_points;
+ }
-template <int dim>
-void
-BoussinesqFlowProblem<dim>::Postprocessor::
-compute_derived_quantities_vector (const std::vector<Vector<double> > &uh,
- const std::vector<std::vector<Tensor<1,dim> > > &duh,
- const std::vector<std::vector<Tensor<2,dim> > > &/*dduh*/,
- const std::vector<Point<dim> > &/*normals*/,
- const std::vector<Point<dim> > &evaluation_points,
- std::vector<Vector<double> > &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 <int dim>
+ void
+ BoussinesqFlowProblem<dim>::Postprocessor::
+ compute_derived_quantities_vector (const std::vector<Vector<double> > &uh,
+ const std::vector<std::vector<Tensor<1,dim> > > &duh,
+ const std::vector<std::vector<Tensor<2,dim> > > &/*dduh*/,
+
const std::vector<Point<dim> > &/*normals*/,
+
const std::vector<Point<dim> > &evaluation_points,
+ std::vector<Vector<double> > &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<n_quadrature_points; ++q)
- {
- // velocity; rescale in cm/year
- for (unsigned int d=0; d<dim; ++d)
- computed_quantities[q](d)
- = (uh[q](d) * EquationData::year_in_seconds * 100);
-
- // pressure
- const double pressure = (uh[q](dim)-minimal_pressure) * EquationData::pressure_scaling;
- computed_quantities[q](dim) = pressure;
-
- // temperature
- const double temperature = uh[q](dim+1);
- computed_quantities[q](dim+1) = temperature;
-
- // friction heating
- Tensor<2,dim> grad_u;
- for (unsigned int d=0; d<dim; ++d)
- grad_u[d] = duh[q][d];
- const SymmetricTensor<2,dim> 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<n_quadrature_points; ++q)
+ {
+ // velocity; rescale in cm/year
+ for (unsigned int d=0; d<dim; ++d)
+ computed_quantities[q](d)
+ = (uh[q](d) * EquationData::year_in_seconds * 100);
+
+ // pressure
+ const double pressure = (uh[q](dim)-minimal_pressure) * EquationData::pressure_scaling;
+ computed_quantities[q](dim) = pressure;
+
+ // temperature
+ const double temperature = uh[q](dim+1);
+ computed_quantities[q](dim+1) = temperature;
+
+ // friction heating
+ Tensor<2,dim> grad_u;
+ for (unsigned int d=0; d<dim; ++d)
+ grad_u[d] = duh[q][d];
+ const SymmetricTensor<2,dim> 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
// therefore rescale things into
// centimeters per year, the unit
// commonly used in geophysics.
-template <int dim>
-void BoussinesqFlowProblem<dim>::output_results ()
-{
- computing_timer.enter_section ("Postprocessing");
+ template <int dim>
+ void BoussinesqFlowProblem<dim>::output_results ()
+ {
+ computing_timer.enter_section ("Postprocessing");
- const FESystem<dim> joint_fe (stokes_fe, 1,
- temperature_fe, 1);
+ const FESystem<dim> joint_fe (stokes_fe, 1,
+ temperature_fe, 1);
- DoFHandler<dim> 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<dim> 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<unsigned int> local_joint_dof_indices (joint_fe.dofs_per_cell);
- std::vector<unsigned int> local_stokes_dof_indices (stokes_fe.dofs_per_cell);
- std::vector<unsigned int> local_temperature_dof_indices (temperature_fe.dofs_per_cell);
+ {
+ std::vector<unsigned int> local_joint_dof_indices (joint_fe.dofs_per_cell);
+ std::vector<unsigned int> local_stokes_dof_indices (stokes_fe.dofs_per_cell);
+ std::vector<unsigned int> local_temperature_dof_indices (temperature_fe.dofs_per_cell);
- typename DoFHandler<dim>::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<joint_fe.dofs_per_cell; ++i)
- if (joint_fe.system_to_base_index(i).first.first == 0)
- {
- Assert (joint_fe.system_to_base_index(i).second
- <
- local_stokes_dof_indices.size(),
- ExcInternalError());
-
- joint_solution(local_joint_dof_indices[i])
- = stokes_solution(local_stokes_dof_indices
- [joint_fe.system_to_base_index(i).second]);
- }
- else
- {
- Assert (joint_fe.system_to_base_index(i).first.first == 1,
- ExcInternalError());
- Assert (joint_fe.system_to_base_index(i).second
- <
- local_temperature_dof_indices.size(),
- ExcInternalError());
- joint_solution(local_joint_dof_indices[i])
- = temperature_solution(local_temperature_dof_indices
- [joint_fe.system_to_base_index(i).second]);
- }
- }
- }
+ typename DoFHandler<dim>::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<joint_fe.dofs_per_cell; ++i)
+ if (joint_fe.system_to_base_index(i).first.first == 0)
+ {
+ Assert (joint_fe.system_to_base_index(i).second
+ <
+ local_stokes_dof_indices.size(),
+ ExcInternalError());
+
+ joint_solution(local_joint_dof_indices[i])
+ = stokes_solution(local_stokes_dof_indices
+ [joint_fe.system_to_base_index(i).second]);
+ }
+ else
+ {
+ Assert (joint_fe.system_to_base_index(i).first.first == 1,
+ ExcInternalError());
+ Assert (joint_fe.system_to_base_index(i).second
+ <
+ local_temperature_dof_indices.size(),
+ ExcInternalError());
+ joint_solution(local_joint_dof_indices[i])
+ = temperature_solution(local_temperature_dof_indices
+ [joint_fe.system_to_base_index(i).second]);
+ }
+ }
+ }
- IndexSet locally_relevant_joint_dofs(joint_dof_handler.n_dofs());
- DoFTools::extract_locally_relevant_dofs (joint_dof_handler, locally_relevant_joint_dofs);
- TrilinosWrappers::MPI::Vector locally_relevant_joint_solution;
- locally_relevant_joint_solution.reinit (locally_relevant_joint_dofs, MPI_COMM_WORLD);
- locally_relevant_joint_solution = joint_solution;
-
- Postprocessor postprocessor (Utilities::System::
- get_this_mpi_process(MPI_COMM_WORLD),
- stokes_solution.block(1).minimal_value());
-
- DataOut<dim> 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<std::string> filenames;
- for (unsigned int i=0; i<Utilities::System::get_n_mpi_processes(MPI_COMM_WORLD); ++i)
- filenames.push_back (std::string("solution-") +
- Utilities::int_to_string (out_index, 5) +
- "." +
- Utilities::int_to_string(i, 4) +
- ".vtu");
- const std::string
- pvtu_master_filename = ("solution-" +
- Utilities::int_to_string (out_index, 5) +
- ".pvtu");
- std::ofstream pvtu_master (pvtu_master_filename.c_str());
- data_out.write_pvtu_record (pvtu_master, filenames);
-
- const std::string
- visit_master_filename = ("solution-" +
- Utilities::int_to_string (out_index, 5) +
- ".visit");
- std::ofstream visit_master (visit_master_filename.c_str());
- data_out.write_visit_record (visit_master, filenames);
- }
+ IndexSet locally_relevant_joint_dofs(joint_dof_handler.n_dofs());
+ DoFTools::extract_locally_relevant_dofs (joint_dof_handler, locally_relevant_joint_dofs);
+ TrilinosWrappers::MPI::Vector locally_relevant_joint_solution;
+ locally_relevant_joint_solution.reinit (locally_relevant_joint_dofs, MPI_COMM_WORLD);
+ locally_relevant_joint_solution = joint_solution;
+
+ Postprocessor postprocessor (Utilities::System::
+ get_this_mpi_process(MPI_COMM_WORLD),
+ stokes_solution.block(1).minimal_value());
+
+ DataOut<dim> 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<std::string> filenames;
+ for (unsigned int i=0; i<Utilities::System::get_n_mpi_processes(MPI_COMM_WORLD); ++i)
+ filenames.push_back (std::string("solution-") +
+ Utilities::int_to_string (out_index, 5) +
+ "." +
+ Utilities::int_to_string(i, 4) +
+ ".vtu");
+ const std::string
+ pvtu_master_filename = ("solution-" +
+ Utilities::int_to_string (out_index, 5) +
+ ".pvtu");
+ std::ofstream pvtu_master (pvtu_master_filename.c_str());
+ data_out.write_pvtu_record (pvtu_master, filenames);
+
+ const std::string
+ visit_master_filename = ("solution-" +
+ Utilities::int_to_string (out_index, 5) +
+ ".visit");
+ std::ofstream visit_master (visit_master_filename.c_str());
+ data_out.write_visit_record (visit_master, filenames);
+ }
- computing_timer.exit_section ();
- out_index++;
-}
+ computing_timer.exit_section ();
+ out_index++;
+ }
// complete set of refinement
// indicators, and the rest of the
// function proceeds as before.
-template <int dim>
-void BoussinesqFlowProblem<dim>::refine_mesh (const unsigned int max_grid_level)
-{
- computing_timer.enter_section ("Refine mesh structure, part 1");
- Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
-
- KellyErrorEstimator<dim>::estimate (temperature_dof_handler,
- QGauss<dim-1>(parameters.temperature_degree+1),
- typename FunctionMap<dim>::type(),
- temperature_solution,
- estimated_error_per_cell,
- std::vector<bool>(),
- 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<dim>::active_cell_iterator
- cell = triangulation.begin_active(max_grid_level);
- cell != triangulation.end(); ++cell)
- cell->clear_refine_flag ();
-
- std::vector<const TrilinosWrappers::MPI::Vector *> x_temperature (2);
- x_temperature[0] = &temperature_solution;
- x_temperature[1] = &old_temperature_solution;
- std::vector<const TrilinosWrappers::MPI::BlockVector *> x_stokes (2);
- x_stokes[0] = &stokes_solution;
- x_stokes[1] = &old_stokes_solution;
-
- parallel::distributed::SolutionTransfer<dim,TrilinosWrappers::MPI::Vector>
- temperature_trans(temperature_dof_handler);
- parallel::distributed::SolutionTransfer<dim,TrilinosWrappers::MPI::BlockVector>
- 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 <int dim>
+ void BoussinesqFlowProblem<dim>::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<float> estimated_error_per_cell (triangulation.n_active_cells());
+
+ KellyErrorEstimator<dim>::estimate (temperature_dof_handler,
+ QGauss<dim-1>(parameters.temperature_degree+1),
+ typename FunctionMap<dim>::type(),
+ temperature_solution,
+ estimated_error_per_cell,
+ std::vector<bool>(),
+ 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<dim>::active_cell_iterator
+ cell = triangulation.begin_active(max_grid_level);
+ cell != triangulation.end(); ++cell)
+ cell->clear_refine_flag ();
+
+ std::vector<const TrilinosWrappers::MPI::Vector *> x_temperature (2);
+ x_temperature[0] = &temperature_solution;
+ x_temperature[1] = &old_temperature_solution;
+ std::vector<const TrilinosWrappers::MPI::BlockVector *> x_stokes (2);
+ x_stokes[0] = &stokes_solution;
+ x_stokes[1] = &old_stokes_solution;
+
+ parallel::distributed::SolutionTransfer<dim,TrilinosWrappers::MPI::Vector>
+ temperature_trans(temperature_dof_handler);
+ parallel::distributed::SolutionTransfer<dim,TrilinosWrappers::MPI::BlockVector>
+ 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<TrilinosWrappers::MPI::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<TrilinosWrappers::MPI::BlockVector *> 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<TrilinosWrappers::MPI::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<TrilinosWrappers::MPI::BlockVector *> 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();
+ }
// function instead of the library function
// <code>VectorTools::project</code>, the
// rest is as before.
-template <int dim>
-void BoussinesqFlowProblem<dim>::run ()
-{
- GridGenerator::hyper_shell (triangulation,
- Point<dim>(),
- EquationData::R0,
- EquationData::R1,
- (dim==3) ? 96 : 12,
- true);
- static HyperShellBoundary<dim> boundary;
- triangulation.set_boundary (0, boundary);
- triangulation.set_boundary (1, boundary);
+ template <int dim>
+ void BoussinesqFlowProblem<dim>::run ()
+ {
+ GridGenerator::hyper_shell (triangulation,
+ Point<dim>(),
+ EquationData::R0,
+ EquationData::R1,
+ (dim==3) ? 96 : 12,
+ true);
+ static HyperShellBoundary<dim> 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 ();
+ }
}
// 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
projects / dealii-svn.git / commitdiff