// @sect3{The <code>HeatEquation</code> implementation}
//
- // The next step then is the implementation of the main class.
+ // It is time now for the implementation of the main class. Let's
+ // start with the constructor which selects a linear element, a time
+ // step constant at 1/500 (remember that one period of the source
+ // on the right hand side was set to 0.2 above, so we resolve each
+ // period with 100 time steps) and chooses the Crank Nicolson method
+ // by setting $\theta=1/2$.
template<int dim>
HeatEquation<dim>::HeatEquation ()
:
dof_handler(triangulation),
time_step(1. / 500),
theta(0.5)
- {
- }
+ {}
+ // @sect4{<code>HeatEquation::setup_system</code>}
+ //
+ // The next function is the one that sets up the DoFHandler object,
+ // computes the constraints, and sets the linear algebra objects
+ // to their correct sizes. We also compute the mass and Laplace
+ // matrix here by simply calling two functions in the library.
+ //
+ // Note that we compute these matrices taking into account already the
+ // constraints due to hanging nodes. These are all homogenous, i.e.,
+ // they only consist of constraints of the form $U_i = \alpha_{ij} U_j
+ // + \alpha_{ik} U_k$ (whereas inhomogenous constraints would also
+ // have a term not proportional to $U$, i.e., $U_i = \alpha_{ij} U_j
+ // + \alpha_{ik} U_k + c_i$). For this kind of constraint, we can
+ // eliminate hanging nodes independently in the matrix and the
+ // right hand side vectors, but this is not the case for inhomogenous
+ // constraints for which we can eliminate constrained degrees of freedom
+ // only by looking at both the system matrix and corresponding right
+ // right hand side at the same time. This may become a problem when
+ // dealing with non-zero Dirichlet boundary conditions, though we
+ // do not do this here in the current program.
template<int dim>
void HeatEquation<dim>::setup_system()
{
}
-
+ // @sect4{<code>HeatEquation::solve_time_step</code>}
+ //
+ // The next function is the one that solves the actual linear system
+ // for a single time step. There is nothing surprising here:
template<int dim>
void HeatEquation<dim>::solve_time_step()
{
PreconditionSSOR<> preconditioner;
preconditioner.initialize(system_matrix, 1.0);
- cg.solve(system_matrix, solution, system_rhs, preconditioner);
+ cg.solve(system_matrix, solution, system_rhs,
+ preconditioner);
constraints.distribute(solution);
+ // @sect4{<code>HeatEquation::output_results</code>}
+ //
+ // Neither is there anything new in generating graphical output:
template<int dim>
void HeatEquation<dim>::output_results() const
{
}
- // @sect4{BoussinesqFlowProblem::refine_mesh}
+ // @sect4{<code>HeatEquation::refine_mesh</code>}
//
- // This function takes care of the adaptive mesh refinement. The three tasks
+ // This function is the interesting part of the program. It takes care of
+ // the adaptive mesh refinement. The three tasks
// this function performs is to first find out which cells to
// refine/coarsen, then to actually do the refinement and eventually
// transfer the solution vectors between the two different grids. The first
// task is simply achieved by using the well-established Kelly error
- // estimator on the temperature (it is the temperature we're mainly
- // interested in for this program, and we need to be accurate in regions of
- // high temperature gradients, also to not have too much numerical
- // diffusion). The second task is to actually do the remeshing. That
- // involves only basic functions as well, such as the
+ // estimator on the solution. The second task is to actually do the
+ // remeshing. That involves only basic functions as well, such as the
// <code>refine_and_coarsen_fixed_fraction</code> that refines those cells
- // with the largest estimated error that together make up 80 per cent of the
+ // with the largest estimated error that together make up 60 per cent of the
// error, and coarsens those cells with the smallest error that make up for
- // a combined 10 per cent of the error.
+ // a combined 40 per cent of the error. Note that for problems such as the
+ // current one where the areas where something is going on are shifting
+ // around, we want to aggressively coarsen so that we can move cells
+ // around to where it is necessary.
//
- // If implemented like this, we would get a program that will not make much
- // progress: Remember that we expect temperature fields that are nearly
- // discontinuous (the diffusivity $\kappa$ is very small after all) and
- // consequently we can expect that a freely adapted mesh will refine further
- // and further into the areas of large gradients. This decrease in mesh size
- // will then be accompanied by a decrease in time step, requiring an
- // exceedingly large number of time steps to solve to a given final time. It
- // will also lead to meshes that are much better at resolving
- // discontinuities after several mesh refinement cycles than in the
- // beginning.
- //
- // In particular to prevent the decrease in time step size and the
- // correspondingly large number of time steps, we limit the maximal
- // refinement depth of the mesh. To this end, after the refinement indicator
- // has been applied to the cells, we simply loop over all cells on the
- // finest level and unselect them from refinement if they would result in
- // too high a mesh level.
+ // As already discussed in the introduction, too small a mesh leads to
+ // too small a time step, whereas too large a mesh leads to too little
+ // resolution. Consequently, after the first two steps, we have two
+ // loops that limit refinement and coarsening to an allowable range of
+ // cells:
template <int dim>
void HeatEquation<dim>::refine_mesh (const unsigned int min_grid_level,
const unsigned int max_grid_level)
GridRefinement::refine_and_coarsen_fixed_fraction (triangulation,
estimated_error_per_cell,
0.6, 0.4);
+
if (triangulation.n_levels() > max_grid_level)
for (typename Triangulation<dim>::active_cell_iterator
cell = triangulation.begin_active(max_grid_level);
// SolutionTransfer class and we have to prepare the solution vectors that
// should be transferred to the new grid (we will lose the old grid once
// we have done the refinement so the transfer has to happen concurrently
- // with refinement). What we definitely need are the current and the old
- // temperature (BDF-2 time stepping requires two old solutions). Since the
- // SolutionTransfer objects only support to transfer one object per dof
- // handler, we need to collect the two temperature solutions in one data
- // structure. Moreover, we choose to transfer the Stokes solution, too,
- // since we need the velocity at two previous time steps, of which only
- // one is calculated on the fly.
+ // with refinement). At the point where we call this function, we will
+ // have just computed the solution, so we no longer need the old_solution
+ // variable (it will be overwritten by the solution just after the mesh
+ // may have been refined, i.e., at the end of the time step; see below).
+ // In other words, we only need the one solution vector, and we copy it
+ // to a temporary object where it is safe from being reset when we further
+ // down below call <code>setup_system()</code>.
//
- // Consequently, we initialize two SolutionTransfer objects for the Stokes
- // and temperature DoFHandler objects, by attaching them to the old dof
- // handlers. With this at place, we can prepare the triangulation and the
- // data vectors for refinement (in this order).
- std::vector<Vector<double> > x_solution (2);
- x_solution[0] = solution;
- x_solution[1] = old_solution;
-
+ // Consequently, we initialize a SolutionTransfer object by attaching
+ // it to the old DoF handler. We then prepare the triangulation and the
+ // data vector for refinement (in this order).
SolutionTransfer<dim> solution_trans(dof_handler);
+ Vector<double> previous_solution;
+ previous_solution = solution;
triangulation.prepare_coarsening_and_refinement();
- solution_trans.prepare_for_coarsening_and_refinement(x_solution);
+ solution_trans.prepare_for_coarsening_and_refinement(previous_solution);
// Now everything is ready, so do the refinement and recreate the dof
// structure on the new grid, and initialize the matrix structures and the
- // new vectors in the <code>setup_dofs</code> function. Next, we actually
- // perform the interpolation of the solutions between the grids. We create
- // another copy of temporary vectors for temperature (now corresponding to
- // the new grid), and let the interpolate function do the job. Then, the
- // resulting array of vectors is written into the respective vector member
- // variables. For the Stokes vector, everything is just the same –
- // except that we do not need another temporary vector since we just
- // interpolate a single vector. In the end, we have to tell the program
- // that the matrices and preconditioners need to be regenerated, since the
- // mesh has changed.
+ // new vectors in the <code>setup_system</code> function. Next, we actually
+ // perform the interpolation of the solution from old to new grid.
triangulation.execute_coarsening_and_refinement ();
setup_system ();
- std::vector<Vector<double> > tmp (2);
- tmp[0].reinit (solution);
- tmp[1].reinit (solution);
- solution_trans.interpolate(x_solution, tmp);
-
- solution = tmp[0];
- old_solution = tmp[1];
+ solution_trans.interpolate(previous_solution, solution);
}
+ // @sect4{<code>HeatEquation::run</code>}
+ //
+ // This is the main driver of the program, where we loop over all
+ // time steps. At the top of the function, we set the number of
+ // initial global mesh refinements and the number of initial cycles of
+ // adaptive mesh refinement by repeating the first time step a few
+ // times. Then we create a mesh, initialize the various objects we will
+ // work with, set a label for where we should start when re-running
+ // the first time step, and interpolate the initial solution onto
+ // out mesh (we choose the zero function here, which of course we could
+ // do in a simpler way by just setting the solution vector to zero). We
+ // also output the initial time step once.
template<int dim>
void HeatEquation<dim>::run()
{
start_time_iteration:
+ tmp.reinit (solution.size());
+ forcing_terms.reinit (solution.size());
+
+
VectorTools::interpolate(dof_handler,
ZeroFunction<dim>(),
old_solution);
output_results();
+ // Then we start the main loop until the computed time exceeds our
+ // end time of 0.5. The first task is to build the right hand
+ // side of the linear system we need to solve in each time step.
+ // Recall that it contains the term $MU^{n-1}-(1-\theta)k_n AU^{n-1}$.
+ // We put these terms into the variable system_rhs, with the
+ // help of a temporary vector:
while (time <= 0.5)
{
time += time_step;
std::cout << "Time step " << timestep_number << " at t=" << time
<< std::endl;
- tmp.reinit (solution.size());
- forcing_terms.reinit (solution.size());
-
mass_matrix.vmult(system_rhs, old_solution);
laplace_matrix.vmult(tmp, old_solution);
system_rhs.add(-(1 - theta) * time_step, tmp);
+ // The second piece is to compute the contributions of the source
+ // terms. This corresponds to the term $k_n
+ // \left[ (1-\theta)F^{n-1} + \theta F^n \right]$. The following
+ // code calls VectorTools::create_right_hand_side to compute the
+ // vectors $F$, where we set the time of the right hand side
+ // (source) function before we evaluate it. The result of this
+ // all ends up in the forcing_terms variable:
RightHandSide<dim> rhs_function;
rhs_function.set_time(time);
VectorTools::create_right_hand_side(dof_handler,
forcing_terms.add(time_step * (1 - theta), tmp);
+ // Next, we add the forcing terms to the ones that
+ // come from the time stepping, and also build the matrix
+ // $M+k_n\theta A$ that we have to invert in each time step.
+ // The final piece of these operations is to eliminate
+ // hanging node constrained degrees of freedom from the
+ // linear system:
system_rhs += forcing_terms;
+ system_matrix.copy_from(mass_matrix);
+ system_matrix.add(theta * time_step, laplace_matrix);
+
+ constraints.condense (system_matrix, system_rhs);
+
+ // There is one more operation we need to do before we
+ // can solve it: boundary values. To this end, we create
+ // a boundary value object, set the proper time to the one
+ // of the current time step, and evaluate it as we have
+ // done many times before. The result is used to also
+ // set the correct boundary values in the linear system:
{
BoundaryValues<dim> boundary_values_function;
boundary_values_function.set_time(time);
boundary_values_function,
boundary_values);
- system_matrix.copy_from(mass_matrix);
- system_matrix.add(theta * time_step, laplace_matrix);
MatrixTools::apply_boundary_values(boundary_values,
system_matrix,
solution,
system_rhs);
}
- constraints.condense (system_matrix, system_rhs);
-
+ // With this out of the way, all we have to do is solve the
+ // system, generate graphical data, and...
solve_time_step();
output_results();
+ // ...take care of mesh refinement. Here, what we want to do is
+ // (i) refine the requested number of times at the very beginning
+ // of the solution procedure, after which we jump to the top to
+ // restart the time iteration, (ii) refine every fifth time
+ // step after that.
+ //
+ // The time loop and, indeed, the main part of the program ends
+ // with starting into the next time step by setting old_solution
+ // to the solution we have just computed.
if ((timestep_number == 1) &&
(pre_refinement_step < n_adaptive_pre_refinement_steps))
{
initial_global_refinement + n_adaptive_pre_refinement_steps);
++pre_refinement_step;
+ tmp.reinit (solution.size());
+ forcing_terms.reinit (solution.size());
+
std::cout << std::endl;
goto start_time_iteration;
}
else if ((timestep_number > 0) && (timestep_number % 5 == 0))
- refine_mesh (initial_global_refinement,
- initial_global_refinement + n_adaptive_pre_refinement_steps);
+ {
+ refine_mesh (initial_global_refinement,
+ initial_global_refinement + n_adaptive_pre_refinement_steps);
+ tmp.reinit (solution.size());
+ forcing_terms.reinit (solution.size());
+ }
old_solution = solution;
}
}
}
+
+// @sect3{The <code>main</code> function}
+//
+// Having made it this far, there is, again, nothing
+// much to discuss for the main function of this
+// program: it looks like all such functions since step-6.
int main()
{
try
projects / dealii-svn.git / commitdiff