Time step 1 at t=0.015625
u-equation: 8 CG iterations.
v-equation: 22 CG iterations.
- Total energy: 0.982265
+ Total energy: 1.17887
Time step 2 at t=0.03125
u-equation: 8 CG iterations.
- v-equation: 23 CG iterations.
- Total energy: 4.10195
+ v-equation: 20 CG iterations.
+ Total energy: 2.9655
Time step 3 at t=0.046875
u-equation: 8 CG iterations.
- v-equation: 20 CG iterations.
- Total energy: 6.95484
+ v-equation: 21 CG iterations.
+ Total energy: 4.33761
Time step 4 at t=0.0625
- u-equation: 8 CG iterations.
+ u-equation: 7 CG iterations.
v-equation: 21 CG iterations.
- Total energy: 7.92609
+ Total energy: 5.35499
Time step 5 at t=0.078125
u-equation: 7 CG iterations.
- v-equation: 22 CG iterations.
- Total energy: 8.9877
-Time step 6 at t=0.09375
- u-equation: 8 CG iterations.
v-equation: 21 CG iterations.
- Total energy: 10.1318
+ Total energy: 6.18652
+Time step 6 at t=0.09375
+ u-equation: 7 CG iterations.
+ v-equation: 20 CG iterations.
+ Total energy: 6.6799
...
Time step 31 at t=0.484375
u-equation: 7 CG iterations.
- v-equation: 21 CG iterations.
- Total energy: 21.6306
+ v-equation: 20 CG iterations.
+ Total energy: 21.9068
Time step 32 at t=0.5
u-equation: 7 CG iterations.
- v-equation: 21 CG iterations.
- Total energy: 23.417
+ v-equation: 20 CG iterations.
+ Total energy: 23.3394
Time step 33 at t=0.515625
u-equation: 7 CG iterations.
- v-equation: 21 CG iterations.
- Total energy: 23.2328
+ v-equation: 20 CG iterations.
+ Total energy: 23.1019
...
Time step 319 at t=4.98438
u-equation: 7 CG iterations.
- v-equation: 21 CG iterations.
- Total energy: 23.2328
+ v-equation: 20 CG iterations.
+ Total energy: 23.1019
Time step 320 at t=5
u-equation: 7 CG iterations.
- v-equation: 21 CG iterations.
- Total energy: 23.2328
+ v-equation: 20 CG iterations.
+ Total energy: 23.1019
@endcode
What we see immediately is that the energy is a constant at least after
be used to pass non-constant coefficient functions to them. The required
changes are therefore relatively small. On the other hand, care must be
taken again to make sure the time step is within the allowed range.
+
+ <li>In the in-code comments, we discussed the fact that the matrices for
+ solving for $U^n$ and $V^n$ need to be reset in every time because of
+ boundary conditions, even though the actual content does not change. It is
+ possible to avoid copying by not eliminating columns in the linear systems,
+ which is implemented by appending a @p false argument to the call:
+ @code
+ MatrixTools::apply_boundary_values (boundary_values,
+ matrix_u,
+ solution_u,
+ system_rhs,
+ false);
+ @endcode
</ul>
// @sect3{The <code>WaveEquation</code> class}
- // Next comes the declaration of the
- // main class. It's public interface
- // of functions is like in most of
- // the other tutorial programs. Worth
- // mentioning is that we now have to
- // store three matrices instead of
- // one: the mass matrix $M$, the
- // Laplace matrix $A$, and the system
- // matrix $M+k^2\theta^2A$ used when
- // solving for $U^n$. Likewise, we
- // need solution vectors for
- // $U^n,V^n$ as well as for the
- // corresponding vectors at the
- // previous time step,
+ // Next comes the declaration of the main
+ // class. It's public interface of functions
+ // is like in most of the other tutorial
+ // programs. Worth mentioning is that we now
+ // have to store four matrices instead of
+ // one: the mass matrix $M$, the Laplace
+ // matrix $A$, the matrix $M+k^2\theta^2A$
+ // used for solving for $U^n$, and a copy of
+ // the mass matrix with boundary conditions
+ // applied used for solving for $V^n$. Note
+ // that it is a bit wasteful to have an
+ // additional copy of the mass matrix
+ // around. We will discuss strategies for how
+ // to avoid this in the section on possible
+ // improvements.
+ //
+ // Likewise, we need solution vectors for
+ // $U^n,V^n$ as well as for the corresponding
+ // vectors at the previous time step,
// $U^{n-1},V^{n-1}$. The
- // <code>system_rhs</code> will be
- // used for whatever right hand side
- // vector we have when solving one of
- // the two linear systems in each
- // time step. These will be solved in
- // the two functions
- // <code>solve_u</code> and
+ // <code>system_rhs</code> will be used for
+ // whatever right hand side vector we have
+ // when solving one of the two linear systems
+ // in each time step. These will be solved in
+ // the two functions <code>solve_u</code> and
// <code>solve_v</code>.
//
// Finally, the variable
ConstraintMatrix constraints;
SparsityPattern sparsity_pattern;
- SparseMatrix<double> system_matrix;
SparseMatrix<double> mass_matrix;
SparseMatrix<double> laplace_matrix;
+ SparseMatrix<double> matrix_u;
+ SparseMatrix<double> matrix_v;
Vector<double> solution_u, solution_v;
Vector<double> old_solution_u, old_solution_v;
// memory on it several times.
//
// After initializing all of these
- // matrices, we call library
- // functions that build the Laplace
- // and mass matrices. All they need
- // is a DoFHandler object and a
- // quadrature formula object that
- // is to be used for numerical
- // integration. Note that in many
- // respects these functions are
- // better than what we would
- // usually do in application
- // programs, for example because
- // they automatically parallelize
- // building the matrices if
- // multiple processors are
- // available in a machine. When we
- // have both of these matrices, we
- // form the third one by copying
- // and adding the first two in
- // appropriate multiples:
- system_matrix.reinit (sparsity_pattern);
+ // matrices, we call library functions that
+ // build the Laplace and mass matrices. All
+ // they need is a DoFHandler object and a
+ // quadrature formula object that is to be
+ // used for numerical integration. Note
+ // that in many respects these functions
+ // are better than what we would usually do
+ // in application programs, for example
+ // because they automatically parallelize
+ // building the matrices if multiple
+ // processors are available in a
+ // machine. The matrices for solving linear
+ // systems will be filled in the run()
+ // method because we need to re-apply
+ // boundary conditions every time step.
mass_matrix.reinit (sparsity_pattern);
laplace_matrix.reinit (sparsity_pattern);
+ matrix_u.reinit (sparsity_pattern);
+ matrix_v.reinit (sparsity_pattern);
MatrixCreator::create_mass_matrix (dof_handler, QGauss<dim>(3),
mass_matrix);
MatrixCreator::create_laplace_matrix (dof_handler, QGauss<dim>(3),
laplace_matrix);
-
- system_matrix.copy_from (mass_matrix);
- system_matrix.add (theta * theta * time_step * time_step, laplace_matrix);
// The rest of the function is spent on
// setting vector sizes to the correct
SolverControl solver_control (1000, 1e-8*system_rhs.l2_norm());
SolverCG<> cg (solver_control);
- cg.solve (system_matrix, solution_u, system_rhs,
+ cg.solve (matrix_u, solution_u, system_rhs,
PreconditionIdentity());
std::cout << " u-equation: " << solver_control.last_step()
SolverControl solver_control (1000, 1e-8*system_rhs.l2_norm());
SolverCG<> cg (solver_control);
- cg.solve (mass_matrix, solution_v, system_rhs,
+ cg.solve (matrix_v, solution_v, system_rhs,
PreconditionIdentity());
std::cout << " v-equation: " << solver_control.last_step()
0,
boundary_values_u_function,
boundary_values);
+
+ // The matrix for solve_u() is the same in
+ // every time steps, so one could think
+ // that it is enough to do this only once
+ // at the beginning of the
+ // simulation. However, since we need to
+ // apply boundary values to the linear
+ // system (which eliminate some matrix rows
+ // and columns and give contributions to
+ // the right hand side), we have to refill
+ // the matrix in every time steps before we
+ // actually apply boundary data. The actual
+ // content is very simple: it is the sum of
+ // the mass matrix and a weighted Laplace
+ // matrix:
+ matrix_u.copy_from (mass_matrix);
+ matrix_u.add (theta * theta * time_step * time_step, laplace_matrix);
MatrixTools::apply_boundary_values (boundary_values,
- system_matrix,
+ matrix_u,
solution_u,
system_rhs);
}
solve_u ();
- // The second step,
- // i.e. solving for $V^n$,
- // works similarly, except that
- // this time the matrix on the
- // left is the mass matrix, and
- // the right hand side is
- // $MV^{n-1} - k\left[ \theta A
- // U^n + (1-\theta)
- // AU^{n-1}\right]$ plus
- // forcing terms. %Boundary
- // values are applied in the
- // same way as before, except
+ // The second step, i.e. solving for
+ // $V^n$, works similarly, except that
+ // this time the matrix on the left is
+ // the mass matrix (which we copy again
+ // in order to be able to apply
+ // boundary conditions, and the right
+ // hand side is $MV^{n-1} - k\left[
+ // \theta A U^n + (1-\theta)
+ // AU^{n-1}\right]$ plus forcing
+ // terms. %Boundary values are applied
+ // in the same way as before, except
// that now we have to use the
// BoundaryValuesV class:
laplace_matrix.vmult (system_rhs, solution_u);
0,
boundary_values_v_function,
boundary_values);
+ matrix_v.copy_from (mass_matrix);
MatrixTools::apply_boundary_values (boundary_values,
- mass_matrix,
+ matrix_v,
solution_v,
system_rhs);
}