private:
void setup_dofs ();
- void assemble_preconditioner ();
- void assemble_system ();
- void assemble_rhs_T ();
+ void assemble_stokes_preconditioner ();
+ void build_stokes_preconditioner ();
+ void assemble_stokes_system ();
+ void assemble_temperature_system ();
double get_maximal_velocity () const;
double get_maximal_temperature () const;
void solve ();
void output_results () const;
void refine_mesh (const unsigned int max_grid_level);
- const unsigned int degree;
+ Triangulation<dim> triangulation;
- Triangulation<dim> triangulation;
- FESystem<dim> fe;
- DoFHandler<dim> dof_handler;
+ const unsigned int degree;
+
+ FESystem<dim> stokes_fe;
+ DoFHandler<dim> stokes_dof_handler;
+ ConstraintMatrix stokes_constraints;
+
+ BlockSparsityPattern stokes_sparsity_pattern;
+ BlockSparseMatrix<double> stokes_matrix;
+ BlockSparsityPattern stokes_preconditioner_sparsity_pattern;
+ BlockSparseMatrix<double> stokes_preconditioner_matrix;
+
+ BlockVector<double> stokes_solution;
+ BlockVector<double> stokes_rhs;
- ConstraintMatrix hanging_node_constraints;
+
+ FE_Q<dim> temperature_fe;
+ DoFHandler<dim> temperature_dof_handler;
+ ConstraintMatrix temperature_constraints;
+
+ SparsityPattern temperature_sparsity_pattern;
+ SparseMatrix<double> temperature_matrix;
- BlockSparsityPattern sparsity_pattern;
- BlockSparseMatrix<double> system_matrix;
+ Vector<double> temperature_solution;
+ Vector<double> old_temperature_solution;
+ Vector<double> old_old_temperature_solution;
+ Vector<double> temperature_rhs;
- BlockSparsityPattern preconditioner_sparsity_pattern;
- BlockSparseMatrix<double> preconditioner_matrix;
double time_step;
double old_time_step;
unsigned int timestep_number;
- BlockVector<double> solution;
- BlockVector<double> old_solution;
- BlockVector<double> old_old_solution;
- BlockVector<double> system_rhs;
-
boost::shared_ptr<PreconditionerTrilinosAmg> Amg_preconditioner;
boost::shared_ptr<SparseILU<double> > Mp_preconditioner;
// @sect4{Initial values}
template <int dim>
-class InitialValues : public Function<dim>
+class TemperatureInitialValues : public Function<dim>
{
public:
- InitialValues () : Function<dim>(dim+2) {}
+ TemperatureInitialValues () : Function<dim>(1) {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
template <int dim>
double
-InitialValues<dim>::value (const Point<dim> &,
+TemperatureInitialValues<dim>::value (const Point<dim> &,
const unsigned int) const
{
return 0;
template <int dim>
void
-InitialValues<dim>::vector_value (const Point<dim> &p,
+TemperatureInitialValues<dim>::vector_value (const Point<dim> &p,
Vector<double> &values) const
{
for (unsigned int c=0; c<this->n_components; ++c)
- values(c) = InitialValues<dim>::value (p, c);
+ values(c) = TemperatureInitialValues<dim>::value (p, c);
}
// in the problem description in the
// introduction.
template <int dim>
-class RightHandSide : public Function<dim>
+class TemperatureRightHandSide : public Function<dim>
{
public:
- RightHandSide () : Function<dim>(dim+2) {}
+ TemperatureRightHandSide () : Function<dim>(1) {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
template <int dim>
double
-RightHandSide<dim>::value (const Point<dim> &p,
- const unsigned int component) const
+TemperatureRightHandSide<dim>::value (const Point<dim> &p,
+ const unsigned int /*component*/) const
{
- // the right hand side only acts on the
- // temperature component
- if (component == dim+1)
- {
- static const Point<dim> source_centers[3]
- = { (dim == 2 ? Point<dim>(.3,.1) : Point<dim>(.3,.5,.1)),
- (dim == 2 ? Point<dim>(.45,.1) : Point<dim>(.45,.5,.1)),
- (dim == 2 ? Point<dim>(.75,.1) : Point<dim>(.75,.5,.1)) };
- static const double source_radius
- = (dim == 2 ? 1./32 : 1./8);
+ static const Point<dim> source_centers[3]
+ = { (dim == 2 ? Point<dim>(.3,.1) : Point<dim>(.3,.5,.1)),
+ (dim == 2 ? Point<dim>(.45,.1) : Point<dim>(.45,.5,.1)),
+ (dim == 2 ? Point<dim>(.75,.1) : Point<dim>(.75,.5,.1)) };
+ static const double source_radius
+ = (dim == 2 ? 1./32 : 1./8);
- return ((source_centers[0].distance (p) < source_radius)
- ||
- (source_centers[1].distance (p) < source_radius)
- ||
- (source_centers[2].distance (p) < source_radius)
- ?
- 1
- :
- 0);
- }
- else
- return 0;
+ return ((source_centers[0].distance (p) < source_radius)
+ ||
+ (source_centers[1].distance (p) < source_radius)
+ ||
+ (source_centers[2].distance (p) < source_radius)
+ ?
+ 1
+ :
+ 0);
}
template <int dim>
void
-RightHandSide<dim>::vector_value (const Point<dim> &p,
- Vector<double> &values) const
+TemperatureRightHandSide<dim>::vector_value (const Point<dim> &p,
+ Vector<double> &values) const
{
for (unsigned int c=0; c<this->n_components; ++c)
- values(c) = RightHandSide<dim>::value (p, c);
+ values(c) = TemperatureRightHandSide<dim>::value (p, c);
}
const BlockVector<double> &src) const;
private:
- const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
+ const SmartPointer<const BlockSparseMatrix<double> > stokes_matrix;
const SmartPointer<const InverseMatrix<SparseMatrix<double>,
PreconditionerMp > > m_inverse;
const PreconditionerA &a_preconditioner;
const PreconditionerA &Apreconditioner
)
:
- system_matrix (&S),
+ stokes_matrix (&S),
m_inverse (&Mpinv),
a_preconditioner (Apreconditioner),
tmp (S.block(1,1).m())
const BlockVector<double> &src) const
{
a_preconditioner.vmult (dst.block(0), src.block(0));
- system_matrix->block(1,0).residual(tmp, dst.block(0), src.block(1));
+ stokes_matrix->block(1,0).residual(tmp, dst.block(0), src.block(1));
tmp *= -1;
m_inverse->vmult (dst.block(1), tmp);
}
template <int dim>
BoussinesqFlowProblem<dim>::BoussinesqFlowProblem (const unsigned int degree)
:
- degree (degree),
triangulation (Triangulation<dim>::maximum_smoothing),
- fe (FE_Q<dim>(degree+1), dim,
- FE_Q<dim>(degree), 1,
- FE_Q<dim>(degree), 1),
- dof_handler (triangulation),
+ degree (degree),
+ stokes_fe (FE_Q<dim>(degree+1), dim,
+ FE_Q<dim>(degree), 1),
+ stokes_dof_handler (triangulation),
+ temperature_fe (degree),
+ temperature_dof_handler (triangulation),
time_step (0),
old_time_step (0),
timestep_number (0),
template <int dim>
void BoussinesqFlowProblem<dim>::setup_dofs ()
{
- dof_handler.distribute_dofs (fe);
- DoFRenumbering::Cuthill_McKee (dof_handler);
- std::vector<unsigned int> block_component (dim+2,0);
- block_component[dim] = 1;
- block_component[dim+1] = 2;
- DoFRenumbering::component_wise (dof_handler, block_component);
-
- // The next step is, as usual,
- // to write some information
- // to the screen. The information
- // that is most interesting during
- // the calculations is the
- // number of degrees of freedom
- // in the individual components,
- // so we count them. The function
- // to do this is the same as the
- // one used in step-22, which
- // uses the grouping of all
- // velocity components into
- // one block as introduced
- // above.
- std::vector<unsigned int> dofs_per_block (3);
- DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component);
- const unsigned int n_u = dofs_per_block[0],
- n_p = dofs_per_block[1],
- n_T = dofs_per_block[2];
+ std::vector<unsigned int> stokes_block_component (dim+1,0);
+ stokes_block_component[dim] = 1;
+
+ {
+ stokes_dof_handler.distribute_dofs (stokes_fe);
+ DoFRenumbering::Cuthill_McKee (stokes_dof_handler);
+ DoFRenumbering::component_wise (stokes_dof_handler, stokes_block_component);
+
+ stokes_constraints.clear ();
+ DoFTools::make_hanging_node_constraints (stokes_dof_handler,
+ stokes_constraints);
+ std::set<unsigned char> no_normal_flux_boundaries;
+ no_normal_flux_boundaries.insert (0);
+ VectorTools::compute_no_normal_flux_constraints (stokes_dof_handler, 0,
+ no_normal_flux_boundaries,
+ stokes_constraints);
+ stokes_constraints.close ();
+ }
+ {
+ temperature_dof_handler.distribute_dofs (temperature_fe);
+ DoFRenumbering::Cuthill_McKee (temperature_dof_handler);
+
+ temperature_constraints.clear ();
+ DoFTools::make_hanging_node_constraints (temperature_dof_handler,
+ temperature_constraints);
+ temperature_constraints.close ();
+ }
+
+ std::vector<unsigned int> stokes_dofs_per_block (2);
+ DoFTools::count_dofs_per_block (stokes_dof_handler, stokes_dofs_per_block,
+ stokes_block_component);
+
+ const unsigned int n_u = stokes_dofs_per_block[0],
+ n_p = stokes_dofs_per_block[1],
+ n_T = temperature_dof_handler.n_dofs();
std::cout << "Number of active cells: "
<< triangulation.n_active_cells()
<< " levels)"
<< std::endl
<< "Number of degrees of freedom: "
- << dof_handler.n_dofs()
+ << n_u + n_p + n_T
<< " (" << n_u << '+' << n_p << '+'<< n_T <<')'
<< std::endl
<< std::endl;
- hanging_node_constraints.clear ();
- DoFTools::make_hanging_node_constraints (dof_handler,
- hanging_node_constraints);
- std::set<unsigned char> no_normal_flux_boundaries;
- no_normal_flux_boundaries.insert (0);
- VectorTools::compute_no_normal_flux_constraints (dof_handler, 0,
- no_normal_flux_boundaries,
- hanging_node_constraints);
- hanging_node_constraints.close ();
// The next step is to
// Observe that we use a
// coupling argument for
// telling the function
- // <code>make_sparsity_pattern</code>
+ // <code>make_stokes_sparsity_pattern</code>
// which components actually
// will hold data and which
// we're going to neglect.
// system matrix structure to
// the sparsity pattern.
{
- system_matrix.clear ();
+ stokes_matrix.clear ();
- BlockCompressedSetSparsityPattern csp (3,3);
+ BlockCompressedSetSparsityPattern csp (2,2);
csp.block(0,0).reinit (n_u, n_u);
csp.block(0,1).reinit (n_u, n_p);
- csp.block(0,2).reinit (n_u, n_T);
csp.block(1,0).reinit (n_p, n_u);
csp.block(1,1).reinit (n_p, n_p);
- csp.block(1,2).reinit (n_p, n_T);
- csp.block(2,0).reinit (n_T, n_u);
- csp.block(2,1).reinit (n_T, n_p);
- csp.block(2,2).reinit (n_T, n_T);
csp.collect_sizes ();
- Table<2,DoFTools::Coupling> coupling (dim+2, dim+2);
+ Table<2,DoFTools::Coupling> coupling (dim+1, dim+1);
// build the sparsity pattern. note
// that all dim velocities couple with
{
case 2:
{
- static const bool coupling_matrix[4][4]
- = {{ 1, 1, 1, 0 },
- { 1, 1, 1, 0 },
+ static const bool coupling_matrix[3][3]
+ = {{ 1, 1, 1},
+ { 1, 1, 1},
- { 1, 1, 0, 0 },
-
- { 0, 0, 0, 1 }};
- for (unsigned int c=0; c<dim+2; ++c)
- for (unsigned int d=0; d<dim+2; ++d)
+ { 1, 1, 0}};
+ for (unsigned int c=0; c<dim+1; ++c)
+ for (unsigned int d=0; d<dim+1; ++d)
if (coupling_matrix[c][d] == true)
coupling[c][d] = DoFTools::always;
else
case 3:
{
- static const bool coupling_matrix[5][5]
- = {{ 1, 1, 1, 1, 0 },
- { 1, 1, 1, 1, 0 },
- { 1, 1, 1, 1, 0 },
+ static const bool coupling_matrix[4][4]
+ = {{ 1, 1, 1, 1},
+ { 1, 1, 1, 1},
+ { 1, 1, 1, 1},
- { 1, 1, 1, 0, 0 },
-
- { 0, 0, 0, 0, 1 }};
- for (unsigned int c=0; c<dim+2; ++c)
- for (unsigned int d=0; d<dim+2; ++d)
+ { 1, 1, 1, 0}};
+ for (unsigned int c=0; c<dim+1; ++c)
+ for (unsigned int d=0; d<dim+1; ++d)
if (coupling_matrix[c][d] == true)
coupling[c][d] = DoFTools::always;
else
Assert (false, ExcNotImplemented());
}
- DoFTools::make_sparsity_pattern (dof_handler, coupling, csp);
- hanging_node_constraints.condense (csp);
- sparsity_pattern.copy_from (csp);
+ DoFTools::make_sparsity_pattern (stokes_dof_handler, coupling, csp);
+ stokes_constraints.condense (csp);
+ stokes_sparsity_pattern.copy_from (csp);
- system_matrix.reinit (sparsity_pattern);
+ stokes_matrix.reinit (stokes_sparsity_pattern);
}
-
+
{
Amg_preconditioner.reset ();
Mp_preconditioner.reset ();
- preconditioner_matrix.clear ();
+ stokes_preconditioner_matrix.clear ();
- BlockCompressedSetSparsityPattern csp (3,3);
+ BlockCompressedSetSparsityPattern csp (2,2);
csp.block(0,0).reinit (n_u, n_u);
csp.block(0,1).reinit (n_u, n_p);
- csp.block(0,2).reinit (n_u, n_T);
csp.block(1,0).reinit (n_p, n_u);
csp.block(1,1).reinit (n_p, n_p);
- csp.block(1,2).reinit (n_p, n_T);
- csp.block(2,0).reinit (n_T, n_u);
- csp.block(2,1).reinit (n_T, n_p);
- csp.block(2,2).reinit (n_T, n_T);
csp.collect_sizes ();
- Table<2,DoFTools::Coupling> coupling (dim+2, dim+2);
+ Table<2,DoFTools::Coupling> coupling (dim+1, dim+1);
// build the sparsity pattern. note
// that all dim velocities couple with
{
case 2:
{
- static const bool coupling_matrix[4][4]
- = {{ 1, 0, 0, 0 },
- { 0, 1, 0, 0 },
+ static const bool coupling_matrix[3][3]
+ = {{ 1, 0, 0},
+ { 0, 1, 0},
- { 0, 0, 1, 0 },
-
- { 0, 0, 0, 0 }};
- for (unsigned int c=0; c<dim+2; ++c)
- for (unsigned int d=0; d<dim+2; ++d)
+ { 0, 0, 1}};
+ for (unsigned int c=0; c<dim+1; ++c)
+ for (unsigned int d=0; d<dim+1; ++d)
if (coupling_matrix[c][d] == true)
coupling[c][d] = DoFTools::always;
else
case 3:
{
- static const bool coupling_matrix[5][5]
- = {{ 1, 0, 0, 0, 0 },
- { 0, 1, 0, 0, 0 },
- { 0, 0, 1, 0, 0 },
+ static const bool coupling_matrix[4][4]
+ = {{ 1, 0, 0, 0},
+ { 0, 1, 0, 0},
+ { 0, 0, 1, 0},
- { 0, 0, 0, 1, 0 },
-
- { 0, 0, 0, 0, 0 }};
- for (unsigned int c=0; c<dim+2; ++c)
- for (unsigned int d=0; d<dim+2; ++d)
+ { 0, 0, 0, 1}};
+ for (unsigned int c=0; c<dim+1; ++c)
+ for (unsigned int d=0; d<dim+1; ++d)
if (coupling_matrix[c][d] == true)
coupling[c][d] = DoFTools::always;
else
default:
Assert (false, ExcNotImplemented());
- }
-
- DoFTools::make_sparsity_pattern (dof_handler, coupling, csp);
- hanging_node_constraints.condense (csp);
- preconditioner_sparsity_pattern.copy_from (csp);
+ }
+
- preconditioner_matrix.reinit (preconditioner_sparsity_pattern);
+ DoFTools::make_sparsity_pattern (stokes_dof_handler, coupling, csp);
+ stokes_constraints.condense (csp);
+ stokes_preconditioner_sparsity_pattern.copy_from (csp);
+
+ stokes_preconditioner_matrix.reinit (stokes_preconditioner_sparsity_pattern);
}
- // As last action in this function,
- // we need to set the vectors
- // for the solution, the old
- // solution (required for
- // time stepping) and the system
- // right hand side to the
- // three-block structure given
- // by velocity, pressure and
- // temperature.
- solution.reinit (3);
- solution.block(0).reinit (n_u);
- solution.block(1).reinit (n_p);
- solution.block(2).reinit (n_T);
- solution.collect_sizes ();
-
- old_solution.reinit (3);
- old_solution.block(0).reinit (n_u);
- old_solution.block(1).reinit (n_p);
- old_solution.block(2).reinit (n_T);
- old_solution.collect_sizes ();
-
- old_old_solution.reinit (3);
- old_old_solution.block(0).reinit (n_u);
- old_old_solution.block(1).reinit (n_p);
- old_old_solution.block(2).reinit (n_T);
- old_old_solution.collect_sizes ();
+ {
+ temperature_matrix.clear ();
+
+ CompressedSetSparsityPattern csp (n_T, n_T);
+ DoFTools::make_sparsity_pattern (temperature_dof_handler, csp);
+ temperature_constraints.condense (csp);
+ temperature_sparsity_pattern.copy_from (csp);
+
+ temperature_matrix.reinit (temperature_sparsity_pattern);
+ }
+
+
- system_rhs.reinit (3);
- system_rhs.block(0).reinit (n_u);
- system_rhs.block(1).reinit (n_p);
- system_rhs.block(2).reinit (n_T);
- system_rhs.collect_sizes ();
+ // As last action in this function,
+ // we need to set the vectors
+ // for the solution, the old
+ // solution (required for
+ // time stepping) and the system
+ // right hand side to the
+ // three-block structure given
+ // by velocity, pressure and
+ // temperature.
+ stokes_solution.reinit (2);
+ stokes_solution.block(0).reinit (n_u);
+ stokes_solution.block(1).reinit (n_p);
+ stokes_solution.collect_sizes ();
+
+ stokes_rhs.reinit (2);
+ stokes_rhs.block(0).reinit (n_u);
+ stokes_rhs.block(1).reinit (n_p);
+ stokes_rhs.collect_sizes ();
+
+ temperature_solution.reinit (n_T);
+ old_temperature_solution.reinit (n_T);
+ old_old_temperature_solution.reinit (n_T);
+
+ temperature_rhs.reinit (n_T);
}
template <int dim>
void
-BoussinesqFlowProblem<dim>::assemble_preconditioner ()
+BoussinesqFlowProblem<dim>::assemble_stokes_preconditioner ()
{
- preconditioner_matrix = 0;
+ stokes_preconditioner_matrix = 0;
QGauss<dim> quadrature_formula(degree+2);
- FEValues<dim> fe_values (fe, quadrature_formula,
- update_JxW_values |
- update_values |
- update_gradients);
- const unsigned int dofs_per_cell = fe.dofs_per_cell;
+ FEValues<dim> stokes_fe_values (stokes_fe, quadrature_formula,
+ update_JxW_values |
+ update_values |
+ update_gradients);
+ const unsigned int dofs_per_cell = stokes_fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
const FEValuesExtractors::Scalar pressure (dim);
typename DoFHandler<dim>::active_cell_iterator
- cell = dof_handler.begin_active(),
- endc = dof_handler.end();
+ cell = stokes_dof_handler.begin_active(),
+ endc = stokes_dof_handler.end();
for (; cell!=endc; ++cell)
{
- fe_values.reinit (cell);
+ stokes_fe_values.reinit (cell);
local_matrix = 0;
for (unsigned int q=0; q<n_q_points; ++q)
{
for (unsigned int k=0; k<dofs_per_cell; ++k)
{
- phi_grad_u[k] = fe_values[velocities].gradient(k,q);
- phi_p[k] = fe_values[pressure].value (k, q);
+ phi_grad_u[k] = stokes_fe_values[velocities].gradient(k,q);
+ phi_p[k] = stokes_fe_values[pressure].value (k, q);
}
for (unsigned int i=0; i<dofs_per_cell; ++i)
local_matrix(i,j) += (scalar_product (phi_grad_u[i], phi_grad_u[j])
+
phi_p[i] * phi_p[j])
- * fe_values.JxW(q);
+ * stokes_fe_values.JxW(q);
}
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
- preconditioner_matrix.add (local_dof_indices[i],
- local_dof_indices[j],
- local_matrix(i,j));
+ stokes_preconditioner_matrix.add (local_dof_indices[i],
+ local_dof_indices[j],
+ local_matrix(i,j));
}
- hanging_node_constraints.condense (preconditioner_matrix);
+ stokes_constraints.condense (stokes_preconditioner_matrix);
+}
+
+
+
+template <int dim>
+void
+BoussinesqFlowProblem<dim>::build_stokes_preconditioner ()
+{
+ std::cout << " Rebuilding Stokes preconditioner..." << std::flush;
+
+
+ // This last step of the assembly
+ // function sets up the preconditioners
+ // used for the solution of the
+ // system. We are going to use an
+ // ILU preconditioner for the
+ // velocity block (to be used
+ // by BlockSchurPreconditioner class)
+ // as well as an ILU preconditioner
+ // for the inversion of the
+ // pressure mass matrix. Recall that
+ // the velocity-velocity block sits
+ // at position (0,0) in the
+ // global system matrix, and
+ // the pressure mass matrix in
+ // (1,1). The
+ // storage of these objects is
+ // as in step-22, that is, we
+ // include them using a
+ // shared pointer structure from the
+ // boost library.
+ //
+ // When all work is done, we
+ // change the flags
+ // <code>rebuild_preconditioner</code>
+ // and
+ // <code>rebuild_matrices</code>
+ // to false.
+ assemble_stokes_preconditioner ();
+
+ Amg_preconditioner = boost::shared_ptr<PreconditionerTrilinosAmg>
+ (new PreconditionerTrilinosAmg());
+
+ const unsigned int n_u = stokes_preconditioner_matrix.block(0,0).m();
+ std::vector<double> null_space (dim * n_u, 0.);
+
+ std::vector<bool> precondition_dof_list (stokes_dof_handler.n_dofs(), false);
+
+ for (unsigned int component=0; component < dim; ++component)
+ {
+ std::vector<bool> precondition_mask (dim + 1, false);
+ precondition_mask[component] = true;
+ DoFTools::extract_dofs (stokes_dof_handler, precondition_mask,
+ precondition_dof_list);
+
+ // TODO: The current implementation
+ // assumes that we are working on
+ // the first components of a system when
+ // writing into the null vector.
+ // Change this to the general case,
+ // probably use something similar as
+ // for block vectors.
+ unsigned int counter = 0;
+ for (unsigned int i=0; i<stokes_dof_handler.n_dofs(); ++i)
+ {
+ if (precondition_dof_list[i])
+ {
+ Assert(i < n_u,
+ ExcMessage("Could not correctly locate "
+ "preconditioner dofs in system!"));
+ null_space [component * n_u + i] = 1.;
+ ++counter;
+ }
+ }
+ Assert (counter == n_u / dim,
+ ExcDimensionMismatch(counter, n_u / dim));
+ }
+
+ Amg_preconditioner->initialize(stokes_preconditioner_matrix.block(0,0),
+ null_space, dim, true, true, false);
+
+ // TODO: we could throw away the (0,0)
+ // block here since things have been
+ // copied over to Trilinos. we need to
+ // keep the (1,1) block, though
+
+ Mp_preconditioner
+ = boost::shared_ptr<SparseILU<double> >
+ (new SparseILU<double>);
+ Mp_preconditioner->initialize (stokes_preconditioner_matrix.block(1,1),
+ SparseILU<double>::AdditionalData());
+
+ std::cout << std::endl;
+
+ rebuild_preconditioner = false;
}
- // @sect4{BoussinesqFlowProblem::assemble_system}
+ // @sect4{BoussinesqFlowProblem::assemble_stokes_system}
//
// The assembly of the Boussinesq
// system is acutally a two-step
// the local dofs compared to
// the global system.
template <int dim>
-void BoussinesqFlowProblem<dim>::assemble_system ()
+void BoussinesqFlowProblem<dim>::assemble_stokes_system ()
{
+ std::cout << " Assembling..." << std::flush;
+
if (rebuild_matrices == true)
- system_matrix=0;
+ stokes_matrix=0;
- system_rhs=0;
+ stokes_rhs=0;
QGauss<dim> quadrature_formula(degree+2);
QGauss<dim-1> face_quadrature_formula(degree+2);
- FEValues<dim> fe_values (fe, quadrature_formula,
- update_values |
- update_quadrature_points |
- update_JxW_values |
- (rebuild_matrices == true
- ?
- update_gradients
- :
- UpdateFlags(0)));
- FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula,
- update_values |
- update_normal_vectors |
- update_quadrature_points |
- update_JxW_values);
-
- const unsigned int dofs_per_cell = fe.dofs_per_cell;
+ FEValues<dim> stokes_fe_values (stokes_fe, quadrature_formula,
+ update_values |
+ update_quadrature_points |
+ update_JxW_values |
+ (rebuild_matrices == true
+ ?
+ update_gradients
+ :
+ UpdateFlags(0)));
+
+ FEValues<dim> temperature_fe_values (temperature_fe, quadrature_formula,
+ update_values);
+
+ FEFaceValues<dim> stokes_fe_face_values (stokes_fe, face_quadrature_formula,
+ update_values |
+ update_normal_vectors |
+ update_quadrature_points |
+ update_JxW_values);
+
+ const unsigned int dofs_per_cell = stokes_fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
const unsigned int n_face_q_points = face_quadrature_formula.size();
std::vector<unsigned int> local_dof_indices (dofs_per_cell);
- // These few declarations provide
- // the structures for the evaluation
- // of inhomogeneous Neumann boundary
- // conditions from the function
- // declaration made above.
- // The vector <code>old_solution_values</code>
- // evaluates the solution
- // at the old time level, since
- // the temperature from the
- // old time level enters the
- // Stokes system as a source
- // term in the momentum equation.
- //
- // Then, we create a variable
- // to hold the Rayleigh number,
- // the measure of buoyancy.
- //
- // The set of vectors we create
- // next hold the evaluations of
- // the basis functions that will
- // be used for creating the
- // matrices. This gives faster
- // access to that data, which
- // increases the performance
- // of the assembly. See step-22
- // for details.
- //
- // The last few declarations
- // are used to extract the
- // individual blocks (velocity,
- // pressure, temperature) from
- // the total FE system.
+ // These few declarations provide
+ // the structures for the evaluation
+ // of inhomogeneous Neumann boundary
+ // conditions from the function
+ // declaration made above.
+ // The vector <code>old_solution_values</code>
+ // evaluates the solution
+ // at the old time level, since
+ // the temperature from the
+ // old time level enters the
+ // Stokes system as a source
+ // term in the momentum equation.
+ //
+ // Then, we create a variable
+ // to hold the Rayleigh number,
+ // the measure of buoyancy.
+ //
+ // The set of vectors we create
+ // next hold the evaluations of
+ // the basis functions that will
+ // be used for creating the
+ // matrices. This gives faster
+ // access to that data, which
+ // increases the performance
+ // of the assembly. See step-22
+ // for details.
+ //
+ // The last few declarations
+ // are used to extract the
+ // individual blocks (velocity,
+ // pressure, temperature) from
+ // the total FE system.
const PressureBoundaryValues<dim> pressure_boundary_values;
std::vector<double> boundary_values (n_face_q_points);
- std::vector<Vector<double> > old_solution_values(n_q_points,
- Vector<double>(dim+2));
+ std::vector<double> old_temperature_values(n_q_points);
const double Rayleigh_number = 10;
const FEValuesExtractors::Vector velocities (0);
const FEValuesExtractors::Scalar pressure (dim);
- const FEValuesExtractors::Scalar temperature (dim+1);
-
- // Now starts the loop over
- // all cells in the problem.
- // The first commands are all
- // very familiar, doing the
- // evaluations of the element
- // basis functions, resetting
- // the local arrays and
- // getting the values of the
- // old solution at the
- // quadrature point. Then we
- // are ready to loop over
- // the quadrature points
- // on the cell.
+
+ // Now start the loop over
+ // all cells in the problem.
+ // The first commands are all
+ // very familiar, doing the
+ // evaluations of the element
+ // basis functions, resetting
+ // the local arrays and
+ // getting the values of the
+ // old solution at the
+ // quadrature point. Then we
+ // are ready to loop over
+ // the quadrature points
+ // on the cell.
typename DoFHandler<dim>::active_cell_iterator
- cell = dof_handler.begin_active(),
- endc = dof_handler.end();
- for (; cell!=endc; ++cell)
+ cell = stokes_dof_handler.begin_active(),
+ endc = stokes_dof_handler.end();
+ typename DoFHandler<dim>::active_cell_iterator
+ temperature_cell = temperature_dof_handler.begin_active();
+
+ for (; cell!=endc; ++cell, ++temperature_cell)
{
- fe_values.reinit (cell);
+ stokes_fe_values.reinit (cell);
+ temperature_fe_values.reinit (temperature_cell);
+
local_matrix = 0;
local_rhs = 0;
- fe_values.get_function_values (old_solution, old_solution_values);
+ temperature_fe_values.get_function_values (old_temperature_solution, old_temperature_values);
for (unsigned int q=0; q<n_q_points; ++q)
{
- const double old_temperature = old_solution_values[q](dim+1);
-
- // Extract the basis relevant
- // terms in the inner products
- // once in advance as shown
- // in step-22 in order to
- // accelerate assembly.
- //
- // Once this is done, we
- // start the loop over the
- // rows and columns of the
- // local matrix and feed
- // the matrix with the relevant
- // products. The right hand
- // side is filled with the
- // forcing term driven by
- // temperature in direction
- // of gravity (which is
- // vertical in our example).
- // Note that the right hand
- // side term is always generated,
- // whereas the matrix
- // contributions are only
- // updated when it is
- // requested by the
- // <code>rebuild_matrices</code>
- // flag.
+ const double old_temperature = old_temperature_values[q];
+
+ // Extract the basis relevant
+ // terms in the inner products
+ // once in advance as shown
+ // in step-22 in order to
+ // accelerate assembly.
+ //
+ // Once this is done, we
+ // start the loop over the
+ // rows and columns of the
+ // local matrix and feed
+ // the matrix with the relevant
+ // products. The right hand
+ // side is filled with the
+ // forcing term driven by
+ // temperature in direction
+ // of gravity (which is
+ // vertical in our example).
+ // Note that the right hand
+ // side term is always generated,
+ // whereas the matrix
+ // contributions are only
+ // updated when it is
+ // requested by the
+ // <code>rebuild_matrices</code>
+ // flag.
for (unsigned int k=0; k<dofs_per_cell; ++k)
{
- phi_u[k] = fe_values[velocities].value (k,q);
+ phi_u[k] = stokes_fe_values[velocities].value (k,q);
if (rebuild_matrices)
{
- grads_phi_u[k] = fe_values[velocities].symmetric_gradient(k,q);
- div_phi_u[k] = fe_values[velocities].divergence (k, q);
- phi_p[k] = fe_values[pressure].value (k, q);
+ grads_phi_u[k] = stokes_fe_values[velocities].symmetric_gradient(k,q);
+ div_phi_u[k] = stokes_fe_values[velocities].divergence (k, q);
+ phi_p[k] = stokes_fe_values[pressure].value (k, q);
}
}
+ // define viscosity
+ const double eta = 1;
+ if (rebuild_matrices)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ local_matrix(i,j) += (eta * grads_phi_u[i] * grads_phi_u[j]
+ - div_phi_u[i] * phi_p[j]
+ - phi_p[i] * div_phi_u[j])
+ * stokes_fe_values.JxW(q);
+
+ const Point<dim> gravity = ( (dim == 2) ? (Point<dim> (0,1)) :
+ (Point<dim> (0,0,1)) );
for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- const Tensor<1,dim> phi_i_u = fe_values[velocities].value (i, q);
-
- // define viscosity
- const double eta = 1;
-
- // build Stokes part of
- // the matrix. we have
- // to rebuild the
- // temperature part of
- // it in each time step
- // due to the fact that
- // it depends on the
- // time step size
- if (rebuild_matrices)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- local_matrix(i,j) += (eta * grads_phi_u[i] * grads_phi_u[j]
- - div_phi_u[i] * phi_p[j]
- - phi_p[i] * div_phi_u[j])
- * fe_values.JxW(q);
-
- const Point<dim> gravity = ( (dim == 2) ? (Point<dim> (0,1)) :
- (Point<dim> (0,0,1)) );
-
- local_rhs(i) += (Rayleigh_number *
- gravity * phi_u[i] * old_temperature)*
- fe_values.JxW(q);
- }
+ local_rhs(i) += (Rayleigh_number *
+ gravity * phi_u[i] * old_temperature)*
+ stokes_fe_values.JxW(q);
}
- // Next follows the assembly
- // of the face terms, result
- // from Neumann boundary
- // conditions. Since these
- // terms only enter the right
- // hand side vector and not
- // the matrix, there is no
- // substantial benefit from
- // extracting the data
- // before using it, so
- // we remain in the lines
- // of step-20 at this point.
+
+ // Next follows the assembly
+ // of the face terms, result
+ // from Neumann boundary
+ // conditions. Since these
+ // terms only enter the right
+ // hand side vector and not
+ // the matrix, there is no
+ // substantial benefit from
+ // extracting the data
+ // before using it, so
+ // we remain in the lines
+ // of step-20 at this point.
for (unsigned int face_no=0;
- face_no<GeometryInfo<dim>::faces_per_cell;
- ++face_no)
- if (cell->at_boundary(face_no))
- {
- fe_face_values.reinit (cell, face_no);
-
- pressure_boundary_values
- .value_list (fe_face_values.get_quadrature_points(),
- boundary_values);
-
- for (unsigned int q=0; q<n_face_q_points; ++q)
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- const Tensor<1,dim>
- phi_i_u = fe_face_values[velocities].value (i, q);
-
- local_rhs(i) += -(phi_i_u *
- fe_face_values.normal_vector(q) *
- boundary_values[q] *
- fe_face_values.JxW(q));
- }
- }
-
- // The last step in the loop
- // over all cells is to
- // enter the local contributions
- // into the global matrix and
- // vector structures to the
- // positions specified in
- // <code>local_dof_indices</code>.
- // Again, we only add the
- // matrix data when it is
- // requested.
+ face_no<GeometryInfo<dim>::faces_per_cell;
+ ++face_no)
+ if (cell->at_boundary(face_no))
+ {
+ stokes_fe_face_values.reinit (cell, face_no);
+
+ pressure_boundary_values
+ .value_list (stokes_fe_face_values.get_quadrature_points(),
+ boundary_values);
+
+ for (unsigned int q=0; q<n_face_q_points; ++q)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ const Tensor<1,dim>
+ phi_i_u = stokes_fe_face_values[velocities].value (i, q);
+
+ local_rhs(i) += -(phi_i_u *
+ stokes_fe_face_values.normal_vector(q) *
+ boundary_values[q] *
+ stokes_fe_face_values.JxW(q));
+ }
+ }
+
+ // The last step in the loop
+ // over all cells is to
+ // enter the local contributions
+ // into the global matrix and
+ // vector structures to the
+ // positions specified in
+ // <code>local_dof_indices</code>.
+ // Again, we only add the
+ // matrix data when it is
+ // requested.
cell->get_dof_indices (local_dof_indices);
if (rebuild_matrices == true)
- {
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- system_matrix.add (local_dof_indices[i],
- local_dof_indices[j],
- local_matrix(i,j));
- }
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ stokes_matrix.add (local_dof_indices[i],
+ local_dof_indices[j],
+ local_matrix(i,j));
for (unsigned int i=0; i<dofs_per_cell; ++i)
- system_rhs(local_dof_indices[i]) += local_rhs(i);
+ stokes_rhs(local_dof_indices[i]) += local_rhs(i);
}
- // Back at the outermost
- // level of this function,
- // we continue the work
- // by condensing hanging
- // node constraints to the
- // right hand side and,
- // possibly, to the matrix.
- if (rebuild_matrices == true)
- hanging_node_constraints.condense (system_matrix);
-
- hanging_node_constraints.condense (system_rhs);
+ // Back at the outermost
+ // level of this function,
+ // we continue the work
+ // by condensing hanging
+ // node constraints to the
+ // right hand side and,
+ // possibly, to the matrix.
+ stokes_constraints.condense (stokes_rhs);
if (rebuild_matrices == true)
{
+ stokes_constraints.condense (stokes_matrix);
+
// std::map<unsigned int,double> boundary_values;
// typename DoFHandler<dim>::active_cell_iterator
// component_mask);
// MatrixTools::apply_boundary_values (boundary_values,
-// system_matrix,
+// stokes_matrix,
// solution,
// system_rhs);
- }
-
- // This last step of the assembly
- // function sets up the preconditioners
- // used for the solution of the
- // system. We are going to use an
- // ILU preconditioner for the
- // velocity block (to be used
- // by BlockSchurPreconditioner class)
- // as well as an ILU preconditioner
- // for the inversion of the
- // pressure mass matrix. Recall that
- // the velocity-velocity block sits
- // at position (0,0) in the
- // global system matrix, and
- // the pressure mass matrix in
- // (1,1). The
- // storage of these objects is
- // as in step-22, that is, we
- // include them using a
- // shared pointer structure from the
- // boost library.
- //
- // When all work is done, we
- // change the flags
- // <code>rebuild_preconditioner</code>
- // and
- // <code>rebuild_matrices</code>
- // to false.
- if (rebuild_preconditioner == true)
- {
- Assert (rebuild_matrices == true,
- ExcMessage ("There is no point in rebuilding the preconditioner "
- "without a rebuilt matrix!"));
-
- std::cout << " Rebuilding preconditioner..." << std::flush;
-
- assemble_preconditioner ();
-
- Amg_preconditioner = boost::shared_ptr<PreconditionerTrilinosAmg>
- (new PreconditionerTrilinosAmg());
-
- const unsigned int n_u = preconditioner_matrix.block(0,0).m();
- std::vector<double> null_space (dim * n_u, 0.);
-
- std::vector<bool> precondition_dof_list (dof_handler.n_dofs(), false);
- std::vector<bool> precondition_mask (dim + 2, false);
-
- for (unsigned int component=0; component < dim; component++)
- {
- precondition_mask[component] = true;
- DoFTools::extract_dofs (dof_handler, precondition_mask,
- precondition_dof_list);
- precondition_mask[component] = false;
-
- // TODO: The current implementation
- // assumes that we are working on
- // the first components of a system when
- // writing into the null vector.
- // Change this to the general case,
- // probably use something similar as
- // for block vectors.
- unsigned int counter = 0;
- for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)
- {
- if (precondition_dof_list[i])
- {
- Assert(i < n_u,
- ExcMessage("Could not correctly locate "
- "preconditioner dofs in system!"));
- null_space [component * n_u + i] = 1.;
- ++counter;
- }
- }
- Assert (counter == n_u / dim,
- ExcDimensionMismatch(counter, n_u / dim));
- }
-
- Amg_preconditioner->initialize(preconditioner_matrix.block(0,0),
- null_space, dim, true, true, false);
-
- // TODO: we could throw away the (0,0)
- // block here since things have been
- // copied over to Trilinos. we need to
- // keep the (1,1) block, though
-
- Mp_preconditioner
- = boost::shared_ptr<SparseILU<double> >
- (new SparseILU<double>);
- Mp_preconditioner->initialize (preconditioner_matrix.block(1,1),
- SparseILU<double>::AdditionalData());
-
- std::cout << std::endl;
-
- rebuild_preconditioner = false;
}
rebuild_matrices = false;
+
+ std::cout << std::endl;
}
template <int dim>
double compute_viscosity(
- const std::vector<Vector<double> > present_solution,
- const std::vector<Vector<double> > old_solution,
- const std::vector<Vector<double> > old_old_solution,
- const std::vector<std::vector<Tensor<1,dim> > > old_solution_grads,
- const std::vector<std::vector<Tensor<1,dim> > > old_old_solution_grads,
- const std::vector<std::vector<Tensor<2,dim> > > old_solution_hessians,
- const std::vector<std::vector<Tensor<2,dim> > > old_old_solution_hessians,
- const std::vector<double> gamma_values,
- const double kappa,
- const double global_u_infty,
- const double global_T_infty,
- const double global_Omega_diameter,
- const double cell_diameter,
- const double old_time_step
- )
+ const std::vector<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<Tensor<2,dim> > &old_temperature_hessians,
+ const std::vector<Tensor<2,dim> > &old_old_temperature_hessians,
+ const std::vector<Vector<double> > &present_stokes_values,
+ const std::vector<double> &gamma_values,
+ const double kappa,
+ const double global_u_infty,
+ const double global_T_infty,
+ const double global_Omega_diameter,
+ const double cell_diameter,
+ const double old_time_step
+)
{
const double beta = 0.1;
const double alpha = 1;
if (global_u_infty == 0)
return 5e-3 * cell_diameter;
- const unsigned int n_q_points = old_solution.size();
+ const unsigned int n_q_points = old_temperature.size();
- // Stage 1: calculate residual
+ // Stage 1: calculate residual
double max_residual = 0;
double max_velocity = 0;
for (unsigned int q=0; q < n_q_points; ++q)
{
- const double dT_dt = (old_solution[q](dim+1) - old_old_solution[q](dim+1))
- / old_time_step;
- double u_grad_T = 0.;
+ Tensor<1,dim> u;
for (unsigned int d=0; d<dim; ++d)
- u_grad_T += present_solution[q](d)*((old_solution_grads[q][dim+1][d] +
- old_old_solution_grads[q][dim+1][d]) / 2);
+ u[d] = present_stokes_values[q](d);
+
+ 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 = kappa
- * (trace(old_solution_hessians[q][dim+1]) +
- trace(old_old_solution_hessians[q][dim+1])) / 2;
+ * (trace(old_temperature_hessians[q]) +
+ trace(old_old_temperature_hessians[q])) / 2;
const double residual
= std::abs((dT_dt + u_grad_T - kappa_Delta_T - gamma_values[q]) *
- std::pow((old_solution[q](dim+1)+old_old_solution[q](dim+1)) / 2,
+ std::pow((old_temperature[q]+old_old_temperature[q]) / 2,
alpha-1.));
- max_residual = std::max (residual, max_residual);
-
- double velocity_squared = 0;
- for (unsigned int d=0; d<dim; ++d)
- velocity_squared += present_solution[q](d) * present_solution[q](d);
-
- max_velocity = std::max (std::sqrt (velocity_squared),
- max_velocity);
+ max_residual = std::max (residual, max_residual);
+ max_velocity = std::max (std::sqrt (u*u), max_velocity);
}
const double global_scaling = global_u_infty * global_T_infty /
- std::pow(global_Omega_diameter, alpha - 2.);
-
+ std::pow(global_Omega_diameter, alpha - 2.);
+
return (beta *
max_velocity *
std::min (cell_diameter,
- // @sect4{BoussinesqFlowProblem::assemble_rhs_T}
+ // @sect4{BoussinesqFlowProblem::assemble_temperature_system}
//
// This function does the second
// part of the assembly work, the
// flags at face level are the
// same as in step-12.
template <int dim>
-void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
+void BoussinesqFlowProblem<dim>::assemble_temperature_system ()
{
const bool use_bdf2_scheme = (timestep_number != 0);
- system_matrix.block(2,2) = 0;
+ temperature_matrix = 0;
+ temperature_rhs = 0;
QGauss<dim> quadrature_formula(degree+2);
- FEValues<dim> fe_values (fe, quadrature_formula,
- update_values | update_gradients |
- update_hessians |
- update_quadrature_points | update_JxW_values);
-
- const unsigned int dofs_per_cell = fe.dofs_per_cell;
+ FEValues<dim> temperature_fe_values (temperature_fe, quadrature_formula,
+ update_values | update_gradients |
+ update_hessians |
+ update_quadrature_points | update_JxW_values);
+ FEValues<dim> stokes_fe_values (stokes_fe, quadrature_formula,
+ update_values);
+
+ const unsigned int dofs_per_cell = temperature_fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
Vector<double> local_rhs (dofs_per_cell);
std::vector<unsigned int> local_dof_indices (dofs_per_cell);
- // Here comes the declaration
- // of vectors to hold the old
- // and present solution values
- // and gradients
- // for both the cell as well as faces
- // to the cell. Next comes the
- // declaration of an object
- // to hold the temperature
- // boundary values and a
- // well-known extractor for
- // accessing the temperature
- // part of the FE system.
- std::vector<Vector<double> > present_solution_values (n_q_points,
- Vector<double>(dim+2));
- std::vector<Vector<double> > old_solution_values(n_q_points,
- Vector<double>(dim+2));
- std::vector<Vector<double> > old_old_solution_values(n_q_points,
- Vector<double>(dim+2));
- std::vector<std::vector<Tensor<1,dim> > >
- old_solution_grads(n_q_points,
- std::vector<Tensor<1,dim> >(dim+2));
- std::vector<std::vector<Tensor<1,dim> > >
- old_old_solution_grads(n_q_points,
- std::vector<Tensor<1,dim> >(dim+2));
- std::vector<std::vector<Tensor<2,dim> > > old_solution_hessians(
- n_q_points,
- std::vector<Tensor<2,dim> >(dim+2));
- std::vector<std::vector<Tensor<2,dim> > > old_old_solution_hessians(
- n_q_points,
- std::vector<Tensor<2,dim> >(dim+2));
+ // Here comes the declaration
+ // of vectors to hold the old
+ // and present solution values
+ // and gradients
+ // for both the cell as well as faces
+ // to the cell. Next comes the
+ // declaration of an object
+ // to hold the temperature
+ // boundary values and a
+ // well-known extractor for
+ // accessing the temperature
+ // part of the FE system.
+ std::vector<Vector<double> > present_stokes_values (n_q_points,
+ Vector<double>(dim+1));
+
+
+ std::vector<double> old_temperature_values (n_q_points);
+ std::vector<double> old_old_temperature_values(n_q_points);
+ std::vector<Tensor<1,dim> > old_temperature_grads(n_q_points);
+ std::vector<Tensor<1,dim> > old_old_temperature_grads(n_q_points);
+ std::vector<Tensor<2,dim> > old_temperature_hessians(n_q_points);
+ std::vector<Tensor<2,dim> > old_old_temperature_hessians(n_q_points);
+
TemperatureBoundaryValues<dim> temperature_boundary_values;
- RightHandSide<dim> right_hand_side;
+ TemperatureRightHandSide<dim> temperature_right_hand_side;
std::vector<double> gamma_values (n_q_points);
- const FEValuesExtractors::Scalar temperature (dim+1);
-
std::vector<double> phi_T (dofs_per_cell);
std::vector<Tensor<1,dim> > grad_phi_T (dofs_per_cell);
const double global_T_infty = get_maximal_temperature();
const double global_Omega_diameter = GridTools::diameter (triangulation);
- // Now, let's start the loop
- // over all cells in the
- // triangulation. The first
- // actions within the loop
- // are, 0as usual, the evaluation
- // of the FE basis functions
- // and the old and present
- // solution at the quadrature
- // points.
+ // Now, let's start the loop
+ // over all cells in the
+ // triangulation. The first
+ // actions within the loop
+ // are, 0as usual, the evaluation
+ // of the FE basis functions
+ // and the old and present
+ // solution at the quadrature
+ // points.
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = temperature_dof_handler.begin_active(),
+ endc = temperature_dof_handler.end();
typename DoFHandler<dim>::active_cell_iterator
- cell = dof_handler.begin_active(),
- endc = dof_handler.end();
+ stokes_cell = stokes_dof_handler.begin_active();
- for (; cell!=endc; ++cell)
+ for (; cell!=endc; ++cell, ++stokes_cell)
{
local_rhs = 0;
local_matrix = 0;
- fe_values.reinit (cell);
+ temperature_fe_values.reinit (cell);
+ stokes_fe_values.reinit (stokes_cell);
- fe_values.get_function_values (solution, present_solution_values);
- fe_values.get_function_values (old_solution, old_solution_values);
- fe_values.get_function_values (old_old_solution, old_old_solution_values);
+ temperature_fe_values.get_function_values (old_temperature_solution, old_temperature_values);
+ temperature_fe_values.get_function_values (old_old_temperature_solution, old_old_temperature_values);
- fe_values.get_function_gradients (old_solution, old_solution_grads);
- fe_values.get_function_gradients (old_old_solution, old_old_solution_grads);
+ temperature_fe_values.get_function_gradients (old_temperature_solution, old_temperature_grads);
+ temperature_fe_values.get_function_gradients (old_old_temperature_solution, old_old_temperature_grads);
- fe_values.get_function_hessians (old_solution, old_solution_hessians);
- fe_values.get_function_hessians (old_old_solution, old_old_solution_hessians);
+ temperature_fe_values.get_function_hessians (old_temperature_solution, old_temperature_hessians);
+ temperature_fe_values.get_function_hessians (old_old_temperature_solution, old_old_temperature_hessians);
- right_hand_side.value_list (fe_values.get_quadrature_points(),
- gamma_values, dim+1);
+ temperature_right_hand_side.value_list (temperature_fe_values.get_quadrature_points(),
+ gamma_values);
+ stokes_fe_values.get_function_values (stokes_solution, present_stokes_values);
+
// build matrix contributions
// define diffusion. take the
// impirically) to keep the
// scheme stable
const double kappa = 1e-6;
- const double nu =
- compute_viscosity (present_solution_values, old_solution_values,
- old_old_solution_values, old_solution_grads, old_old_solution_grads,
- old_solution_hessians, old_old_solution_hessians, gamma_values,
- kappa, global_u_infty, global_T_infty,
- global_Omega_diameter, cell->diameter(),
- old_time_step);
+ const double nu
+ = compute_viscosity (old_temperature_values,
+ old_old_temperature_values,
+ old_temperature_grads,
+ old_old_temperature_grads,
+ old_temperature_hessians,
+ old_old_temperature_hessians,
+ present_stokes_values,
+ gamma_values,
+ kappa, global_u_infty, global_T_infty,
+ global_Omega_diameter, cell->diameter(),
+ old_time_step);
for (unsigned int q=0; q<n_q_points; ++q)
{
for (unsigned int k=0; k<dofs_per_cell; ++k)
{
- grad_phi_T[k] = fe_values[temperature].gradient(k,q);
- phi_T[k] = fe_values[temperature].value (k, q);
+ grad_phi_T[k] = temperature_fe_values.shape_grad (k,q);
+ phi_T[k] = temperature_fe_values.shape_value (k, q);
}
- const double old_T = old_solution_values[q](dim+1);
- const double old_old_T = old_old_solution_values[q](dim+1);
+ const double old_T = old_temperature_values[q];
+ const double old_old_T = old_old_temperature_values[q];
- const Tensor<1,dim> old_grad_T = old_solution_grads[q][dim+1];
- const Tensor<1,dim> old_old_grad_T = old_old_solution_grads[q][dim+1];
+ const Tensor<1,dim> old_grad_T = old_temperature_grads[q];
+ const Tensor<1,dim> old_old_grad_T = old_old_temperature_grads[q];
Tensor<1,dim> present_u;
for (unsigned int d=0; d<dim; ++d)
- present_u[d] = present_solution_values[q](d);
+ present_u[d] = present_stokes_values[q](d);
if (use_bdf2_scheme == true)
{
+
time_step *
kappa * grad_phi_T[i] * grad_phi_T[j])
- * fe_values.JxW(q);
+ * temperature_fe_values.JxW(q);
for (unsigned int i=0; i<dofs_per_cell; ++i)
local_rhs(i) += ((time_step + old_time_step) / old_time_step *
time_step *
gamma_values[q] * phi_T[i])
*
- fe_values.JxW(q);
+ temperature_fe_values.JxW(q);
}
else
{
+
time_step *
kappa * grad_phi_T[i] * grad_phi_T[j])
- * fe_values.JxW(q);
+ * temperature_fe_values.JxW(q);
for (unsigned int i=0; i<dofs_per_cell; ++i)
local_rhs(i) += (old_T * phi_T[i]
time_step *
gamma_values[q] * phi_T[i])
*
- fe_values.JxW(q);
+ temperature_fe_values.JxW(q);
}
}
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
- system_matrix.add (local_dof_indices[i],
- local_dof_indices[j],
- local_matrix(i,j));
+ temperature_matrix.add (local_dof_indices[i],
+ local_dof_indices[j],
+ local_matrix(i,j));
for (unsigned int i=0; i<dofs_per_cell; ++i)
- system_rhs(local_dof_indices[i]) += local_rhs(i);
+ temperature_rhs(local_dof_indices[i]) += local_rhs(i);
}
}
template <int dim>
void BoussinesqFlowProblem<dim>::solve ()
{
- solution = old_solution;
-
- // Use the BlockMatrixArray structure
- // for extracting only the upper left
- // 2x2 blocks from the matrix that will
- // be used for the solution of the
- // blocked system.
+ std::cout << " Solving..." << std::endl;
+
+ // Use the BlockMatrixArray structure
+ // for extracting only the upper left
+ // 2x2 blocks from the matrix that will
+ // be used for the solution of the
+ // blocked system.
{
- GrowingVectorMemory<Vector<double> > simple_mem;
- BlockMatrixArray<double> stokes_submatrix(2, 2, simple_mem);
-
- stokes_submatrix.enter(system_matrix.block(0,0),0,0);
- stokes_submatrix.enter(system_matrix.block(0,1),0,1);
- stokes_submatrix.enter(system_matrix.block(1,0),1,0);
-
- // Define some temporary vectors
- // for the solution process.
- // TODO: Can we somehow avoid copying
- // these vectors back and forth? I.e.
- // accessing the block vectors in a
- // similar way as the matrix with the
- // BlockMatrixArray class?
- std::vector<unsigned int> block_sizes(2);
- block_sizes[0] = solution.block(0).size();
- block_sizes[1] = solution.block(1).size();
-
- BlockVector<double> up_rhs(block_sizes);
- BlockVector<double> up(block_sizes);
-
- up_rhs.block(0) = system_rhs.block(0);
- up_rhs.block(1) = system_rhs.block(1);
-
- // Set up inverse matrix for
- // pressure mass matrix
+ // Set up inverse matrix for
+ // pressure mass matrix
InverseMatrix<SparseMatrix<double>,SparseILU<double> >
- mp_inverse (preconditioner_matrix.block(1,1), *Mp_preconditioner);
+ mp_inverse (stokes_preconditioner_matrix.block(1,1), *Mp_preconditioner);
- // Set up block Schur preconditioner
- /*BlockSchurPreconditioner<typename InnerPreconditioner<dim>::type,
- SparseILU<double> >
- preconditioner (system_matrix, mp_inverse, *A_preconditioner);*/
+ // Set up block Schur preconditioner
+ /*BlockSchurPreconditioner<typename InnerPreconditioner<dim>::type,
+ SparseILU<double> >
+ preconditioner (stokes_matrix, mp_inverse, *A_preconditioner);*/
BlockSchurPreconditioner<PreconditionerTrilinosAmg, SparseILU<double> >
- preconditioner (system_matrix, mp_inverse, *Amg_preconditioner);
+ preconditioner (stokes_matrix, mp_inverse, *Amg_preconditioner);
- // Set up GMRES solver and
- // solve.
- SolverControl solver_control (system_matrix.m(),
- 1e-6*system_rhs.l2_norm());
+ // Set up GMRES solver and
+ // solve.
+ SolverControl solver_control (stokes_matrix.m(),
+ 1e-6*stokes_rhs.l2_norm());
SolverGMRES<BlockVector<double> > gmres(solver_control,
- SolverGMRES<BlockVector<double> >::AdditionalData(100));
+ SolverGMRES<BlockVector<double> >::AdditionalData(100));
- gmres.solve(stokes_submatrix, up, up_rhs, preconditioner);
+ gmres.solve(stokes_matrix, stokes_solution, stokes_rhs, preconditioner);
std::cout << " "
<< solver_control.last_step()
<< " GMRES iterations for Stokes subsystem."
<< std::endl;
- solution.block(0) = up.block(0);
- solution.block(1) = up.block(1);
-
- // Produce a constistent solution
- // field (we can't do this on the 'up'
- // vector since it does not have the
- // temperature component, but
- // hanging_node_constraints has
- // constraints also for the
- // temperature vector)
- hanging_node_constraints.distribute (solution);
+ // Produce a constistent solution
+ // field (we can't do this on the 'up'
+ // vector since it does not have the
+ // temperature component, but
+ // hanging_node_constraints has
+ // constraints also for the
+ // temperature vector)
+ stokes_constraints.distribute (stokes_solution);
}
// TODO: determine limit of stability for
old_time_step = time_step;
time_step = GridTools::minimal_cell_diameter(triangulation) /
std::max (get_maximal_velocity(), .05) / 4;
+
+ temperature_solution = old_temperature_solution;
- assemble_rhs_T ();
+
+ assemble_temperature_system ();
{
- SolverControl solver_control (system_matrix.block(2,2).m(),
- 1e-8*system_rhs.block(2).l2_norm());
+ SolverControl solver_control (temperature_matrix.m(),
+ 1e-8*temperature_rhs.l2_norm());
SolverCG<> cg (solver_control);
PreconditionSSOR<> preconditioner;
- preconditioner.initialize (system_matrix.block(2,2), 1.2);
+ preconditioner.initialize (temperature_matrix, 1.2);
- cg.solve (system_matrix.block(2,2), solution.block(2),
- system_rhs.block(2), preconditioner);
+ cg.solve (temperature_matrix, temperature_solution,
+ temperature_rhs, preconditioner);
// produce a consistent temperature field
- hanging_node_constraints.distribute (solution);
+ temperature_constraints.distribute (temperature_solution);
std::cout << " "
<< solver_control.last_step()
<< " CG iterations for temperature."
<< std::endl;
std::cout << " Max temperature: "
- << *std::max_element (solution.block(2).begin(),
- solution.block(2).end())
+ << *std::max_element (temperature_solution.begin(),
+ temperature_solution.end())
<< std::endl;
}
}
if (timestep_number % 10 != 0)
return;
- std::vector<std::string> solution_names (dim, "velocity");
- solution_names.push_back ("p");
- solution_names.push_back ("T");
+//TODO!!
+// Assert (false, ExcNotImplemented());
+
+// std::vector<std::string> solution_names (dim, "velocity");
+// solution_names.push_back ("p");
+// solution_names.push_back ("T");
- DataOut<dim> data_out;
+// DataOut<dim> data_out;
- data_out.attach_dof_handler (dof_handler);
+// data_out.attach_dof_handler (dof_handler);
- std::vector<DataComponentInterpretation::DataComponentInterpretation>
- data_component_interpretation
- (dim+2, DataComponentInterpretation::component_is_scalar);
- for (unsigned int i=0; i<dim; ++i)
- data_component_interpretation[i]
- = DataComponentInterpretation::component_is_part_of_vector;
+// std::vector<DataComponentInterpretation::DataComponentInterpretation>
+// data_component_interpretation
+// (dim+2, DataComponentInterpretation::component_is_scalar);
+// for (unsigned int i=0; i<dim; ++i)
+// data_component_interpretation[i]
+// = DataComponentInterpretation::component_is_part_of_vector;
- data_out.add_data_vector (solution, solution_names,
- DataOut<dim>::type_dof_data,
- data_component_interpretation);
- data_out.build_patches (degree);
+// data_out.add_data_vector (solution, solution_names,
+// DataOut<dim>::type_dof_data,
+// data_component_interpretation);
+// data_out.build_patches (degree);
- std::ostringstream filename;
- filename << "solution-" << Utilities::int_to_string(timestep_number, 4) << ".vtk";
+// std::ostringstream filename;
+// filename << "solution-" << Utilities::int_to_string(timestep_number, 4) << ".vtk";
- std::ofstream output (filename.str().c_str());
- data_out.write_vtk (output);
+// std::ofstream output (filename.str().c_str());
+// data_out.write_vtk (output);
}
{
Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
- std::vector<bool> component_mask (dim+2, false);
- component_mask[dim+1] = true;
- KellyErrorEstimator<dim>::estimate (dof_handler,
- QGauss<dim-1>(3),
+ KellyErrorEstimator<dim>::estimate (temperature_dof_handler,
+ QGauss<dim-1>(degree+1),
typename FunctionMap<dim>::type(),
- solution,
- estimated_error_per_cell,
- component_mask);
+ temperature_solution,
+ estimated_error_per_cell);
GridRefinement::refine_and_coarsen_fixed_fraction (triangulation,
estimated_error_per_cell,
if (cell->has_children() == false)
cell->clear_refine_flag ();
- SolutionTransfer<dim, double> soltrans(dof_handler);
+ SolutionTransfer<dim, double> soltrans(temperature_dof_handler);
triangulation.prepare_coarsening_and_refinement();
std::vector<Vector<double> > x_solution (2);
- x_solution[0].reinit (dof_handler.n_dofs());
- x_solution[0] = solution;
- x_solution[1].reinit (dof_handler.n_dofs());
- x_solution[1] = old_solution;
+ x_solution[0].reinit (temperature_dof_handler.n_dofs());
+ x_solution[0] = temperature_solution;
+ x_solution[1].reinit (temperature_dof_handler.n_dofs());
+ x_solution[1] = old_temperature_solution;
soltrans.prepare_for_coarsening_and_refinement(x_solution);
setup_dofs ();
std::vector<Vector<double> > tmp (2);
- tmp[0].reinit (dof_handler.n_dofs());
- tmp[1].reinit (dof_handler.n_dofs());
+ tmp[0].reinit (temperature_dof_handler.n_dofs());
+ tmp[1].reinit (temperature_dof_handler.n_dofs());
soltrans.interpolate(x_solution, tmp);
+ temperature_solution = tmp[0];
+ //TODO: this should be [1]!!!!!!!
+ old_temperature_solution = tmp[0];
+
rebuild_matrices = true;
rebuild_preconditioner = true;
-
- solution = tmp[0];
- old_solution = tmp[0];
}
const unsigned int n_q_points
= quadrature_formula.size();
- FEValues<dim> fe_values (fe, quadrature_formula,
+ FEValues<dim> fe_values (stokes_fe, quadrature_formula,
update_values);
- std::vector<Vector<double> > solution_values(n_q_points,
- Vector<double>(dim+2));
+ std::vector<Vector<double> > stokes_values(n_q_points,
+ Vector<double>(dim+1));
double max_velocity = 0;
typename DoFHandler<dim>::active_cell_iterator
- cell = dof_handler.begin_active(),
- endc = dof_handler.end();
+ cell = stokes_dof_handler.begin_active(),
+ endc = stokes_dof_handler.end();
for (; cell!=endc; ++cell)
{
fe_values.reinit (cell);
- fe_values.get_function_values (solution, solution_values);
+ fe_values.get_function_values (stokes_solution, stokes_values);
for (unsigned int q=0; q<n_q_points; ++q)
{
Tensor<1,dim> velocity;
for (unsigned int i=0; i<dim; ++i)
- velocity[i] = solution_values[q](i);
+ velocity[i] = stokes_values[q](i);
max_velocity = std::max (max_velocity,
velocity.norm());
double BoussinesqFlowProblem<dim>::get_maximal_temperature () const
{
QGauss<dim> quadrature_formula(degree+2);
- const unsigned int n_q_points
- = quadrature_formula.size();
+ const unsigned int n_q_points = quadrature_formula.size();
- FEValues<dim> fe_values (fe, quadrature_formula,
+ FEValues<dim> fe_values (temperature_fe, quadrature_formula,
update_values);
- std::vector<Vector<double> > old_solution_values(n_q_points,
- Vector<double>(dim+2));
- std::vector<Vector<double> > old_old_solution_values(n_q_points,
- Vector<double>(dim+2));
+ std::vector<double> old_temperature_values(n_q_points);
+ std::vector<double> old_old_temperature_values(n_q_points);
+
double max_temperature = 0;
typename DoFHandler<dim>::active_cell_iterator
- cell = dof_handler.begin_active(),
- endc = dof_handler.end();
+ cell = temperature_dof_handler.begin_active(),
+ endc = temperature_dof_handler.end();
for (; cell!=endc; ++cell)
{
fe_values.reinit (cell);
- fe_values.get_function_values (old_solution, old_solution_values);
- fe_values.get_function_values (old_old_solution, old_old_solution_values);
+ fe_values.get_function_values (old_temperature_solution, old_temperature_values);
+ fe_values.get_function_values (old_old_temperature_solution, old_old_temperature_values);
for (unsigned int q=0; q<n_q_points; ++q)
{
double temperature =
- (1. + time_step/old_time_step) * old_solution_values[q](dim+1)-
- time_step/old_time_step * old_old_solution_values[q](dim+1);
+ (1. + time_step/old_time_step) * old_temperature_values[q]-
+ time_step/old_time_step * old_old_temperature_values[q];
max_temperature = std::max (max_temperature,
temperature);
start_time_iteration:
- VectorTools::project (dof_handler,
- hanging_node_constraints,
+ VectorTools::project (temperature_dof_handler,
+ temperature_constraints,
QGauss<dim>(degree+2),
- InitialValues<dim>(),
- old_solution);
+ TemperatureInitialValues<dim>(),
+ old_temperature_solution);
timestep_number = 0;
double time = 0;
<< ", dt=" << time_step
<< std::endl;
- std::cout << " Assembling..." << std::endl;
- assemble_system ();
+ assemble_stokes_system ();
+
+ if (rebuild_preconditioner == true)
+ build_stokes_preconditioner ();
- std::cout << " Solving..." << std::endl;
solve ();
output_results ();
time += time_step;
++timestep_number;
- old_old_solution = old_solution;
- old_solution = solution;
+ old_old_temperature_solution = old_temperature_solution;
+ old_temperature_solution = temperature_solution;
}
while (time <= 10);
}