From 549802b450b200debce765abffd3db3c7f2ee06b Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 14 Jan 2011 10:56:56 +0000 Subject: [PATCH] Fix application of boundary values: Need to reset the matrix in every time step. git-svn-id: https://svn.dealii.org/trunk@23189 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 14 ++- deal.II/examples/step-23/doc/results.dox | 57 ++++++---- deal.II/examples/step-23/step-23.cc | 137 +++++++++++++---------- 3 files changed, 124 insertions(+), 84 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 99ac7a42f9..ea51e0f20a 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -31,7 +31,19 @@ inconvenience this causes.

General

    -
  1. Extended: base/tensor.h has an additional collection of contractions between three tensors (ie. contract3). This can be useful for writing matrix/vector assembly in a more compact form than before.
    +
  2. Fixed: Boundary conditions in the step-23 tutorial program are now +applied correctly. Matrix columns get eliminated with the used method +and introduce some contribution to the right hand side coming from +inhomogeneous boundary values. The old implementation did not reset the +matrix columns before applying new boundary values.
    +(Martin Stoll, Martin Kronbichler, 2011/01/14) +
+ +
    +
  1. Extended: base/tensor.h has an additional collection of +contractions between three tensors (ie. contract3). +This can be useful for writing matrix/vector assembly in a more compact +form than before.
    (Toby D. Young, 2011/01/12)
diff --git a/deal.II/examples/step-23/doc/results.dox b/deal.II/examples/step-23/doc/results.dox index ef5118cab5..02fa8d7058 100644 --- a/deal.II/examples/step-23/doc/results.dox +++ b/deal.II/examples/step-23/doc/results.dox @@ -8,53 +8,53 @@ Number of degrees of freedom: 16641 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 @@ -108,4 +108,17 @@ If you want to explore a bit, try out some of the following things: 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. + +
  • 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 diff --git a/deal.II/examples/step-23/step-23.cc b/deal.II/examples/step-23/step-23.cc index 96804c5315..23db0571d8 100644 --- a/deal.II/examples/step-23/step-23.cc +++ b/deal.II/examples/step-23/step-23.cc @@ -116,28 +116,31 @@ using namespace dealii; // @sect3{The WaveEquation 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 - // system_rhs 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 - // solve_u and + // system_rhs 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 solve_u and // solve_v. // // Finally, the variable @@ -167,9 +170,10 @@ class WaveEquation ConstraintMatrix constraints; SparsityPattern sparsity_pattern; - SparseMatrix system_matrix; SparseMatrix mass_matrix; SparseMatrix laplace_matrix; + SparseMatrix matrix_u; + SparseMatrix matrix_v; Vector solution_u, solution_v; Vector old_solution_u, old_solution_v; @@ -420,36 +424,30 @@ void WaveEquation::setup_system () // 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(3), mass_matrix); MatrixCreator::create_laplace_matrix (dof_handler, QGauss(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 @@ -498,7 +496,7 @@ void WaveEquation::solve_u () 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() @@ -513,7 +511,7 @@ void WaveEquation::solve_v () 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() @@ -676,26 +674,42 @@ void WaveEquation::run () 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); @@ -718,8 +732,9 @@ void WaveEquation::run () 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); } -- 2.39.5