]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix application of boundary values: Need to reset the matrix in every time step.
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 14 Jan 2011 10:56:56 +0000 (10:56 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 14 Jan 2011 10:56:56 +0000 (10:56 +0000)
git-svn-id: https://svn.dealii.org/trunk@23189 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/examples/step-23/doc/results.dox
deal.II/examples/step-23/step-23.cc

index 99ac7a42f98b6f57233f2888d618b5fd062456dd..ea51e0f20a52405c556b9b29ad52f878d92a03e8 100644 (file)
@@ -31,7 +31,19 @@ inconvenience this causes.
 <h3>General</h3>
 
 <ol>
-<li> Extended: <code>base/tensor.h</code> has an additional collection of contractions between three tensors (<i>ie</i>. <code>contract3</code>). This can be useful for writing matrix/vector assembly in a more compact form than before.<br>
+<li> 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.<br>
+(Martin Stoll, Martin Kronbichler, 2011/01/14)
+</ol>
+
+<ol>
+<li> Extended: <code>base/tensor.h</code> has an additional collection of
+contractions between three tensors (<i>ie</i>. <code>contract3</code>).
+This can be useful for writing matrix/vector assembly in a more compact
+form than before.<br>
 (Toby D. Young, 2011/01/12)
 </ol>
 
index ef5118cab5db149b6ddd79486bb70b2efaf2ed03..02fa8d70587138577025ec441c4cf51227378f19 100644 (file)
@@ -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.
+
+  <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>
index 96804c531507cd2f88ad142dac1925274864d259..23db0571d8b532f2668e0b7df91e0fdd14f19513 100644 (file)
@@ -116,28 +116,31 @@ using namespace dealii;
 
                                 // @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
@@ -167,9 +170,10 @@ class WaveEquation
     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;
@@ -420,36 +424,30 @@ void WaveEquation<dim>::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<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
@@ -498,7 +496,7 @@ void WaveEquation<dim>::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<dim>::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<dim>::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<dim>::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);
       }

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.