From 92d109ed776c9f7ae3f2d8a9f6d965ae6240012f Mon Sep 17 00:00:00 2001 From: David Wells Date: Sat, 25 Jun 2016 19:27:45 -0400 Subject: [PATCH] Get rid of a temporary matrix in step-25. One can just use the action of the Laplace and mass matrices to make this easier to understand. --- examples/step-25/step-25.cc | 39 ++++++++++++------------------------- 1 file changed, 12 insertions(+), 27 deletions(-) diff --git a/examples/step-25/step-25.cc b/examples/step-25/step-25.cc index c7f18d23c2..b63edf7b1a 100644 --- a/examples/step-25/step-25.cc +++ b/examples/step-25/step-25.cc @@ -355,43 +355,28 @@ namespace Step25 compute_nl_matrix (old_solution, solution, tmp_matrix); system_matrix.add (std::pow(time_step*theta,2), tmp_matrix); - // Then, we compute the right-hand side vector $-F_h(U^{n,l})$. - // - // We have to first build up the matrix - // $M+k^2\theta^2 A$, which we put into tmp_matrix - // use it to compute a contribution to the right hand side vector, and - // then build the matrix $M-k^2\theta(1-\theta) A$. We could - // build it in the same way as before, i.e., using code like - // @code - // tmp_matrix.copy_from (mass_matrix); - // tmp_matrix.add (-std::pow(time_step,2)*theta*(1-theta), laplace_matrix); - // @endcode - // but we can save the expense of the copy_from operation - // by starting from what is already in the tmp_matrix - // variable (i.e., $M+k^2\theta^2 A$) and subtracting from this - // $k^2\theta^2 A+k^2\theta(1-\theta) A=k^2\theta A$ when computing the - // second matrix: - system_rhs = 0; - - tmp_matrix.copy_from (mass_matrix); - tmp_matrix.add (std::pow(time_step*theta,2), laplace_matrix); + // Next we compute the right-hand side vector. This is just the + // combination of matrix-vector products implied by the description of + // $-F_h(U^{n,l})$ in the introduction. + system_rhs = 0.; Vector tmp_vector (solution.size()); - tmp_matrix.vmult (tmp_vector, solution); - system_rhs += tmp_vector; - - tmp_matrix.add(-std::pow(time_step, 2) * theta, laplace_matrix); + mass_matrix.vmult (system_rhs, solution); + laplace_matrix.vmult (tmp_vector, solution); + system_rhs.add (std::pow(time_step*theta, 2), tmp_vector); - tmp_matrix.vmult (tmp_vector, old_solution); - system_rhs -= tmp_vector; + mass_matrix.vmult (tmp_vector, old_solution); + system_rhs.add (-1.0, tmp_vector); + laplace_matrix.vmult (tmp_vector, old_solution); + system_rhs.add (std::pow(time_step,2)*theta*(1 - theta), tmp_vector); system_rhs.add (-time_step, M_x_velocity); compute_nl_term (old_solution, solution, tmp_vector); system_rhs.add (std::pow(time_step,2)*theta, tmp_vector); - system_rhs *= -1; + system_rhs *= -1.; } // @sect4{SineGordonProblem::compute_nl_term} -- 2.39.5