From 57578b87b60a738634d32c80fefccb3987943df2 Mon Sep 17 00:00:00 2001 From: bangerth Date: Thu, 7 Sep 2006 23:02:51 +0000 Subject: [PATCH] Add somethings git-svn-id: https://svn.dealii.org/trunk@13854 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-23/doc/intro.dox | 8 + deal.II/examples/step-23/step-23.cc | 314 ++++++++++++++++++------- 2 files changed, 238 insertions(+), 84 deletions(-) diff --git a/deal.II/examples/step-23/doc/intro.dox b/deal.II/examples/step-23/doc/intro.dox index bfb31d8df4..78172918f5 100644 --- a/deal.II/examples/step-23/doc/intro.dox +++ b/deal.II/examples/step-23/doc/intro.dox @@ -193,6 +193,10 @@ time stepping schemes do not conserve it after time discretization). Since $v^n$ also appears in the equation for $u^n$, the Crank-Nicolson scheme is also implicit. +In the program, we will leave $\theta$ as a parameter, so that it will +be easy to play with it. The results section will show some numerical +evidence comparing the different schemes. + The equations above (called the semidiscretized equations because we have only discretized the time, but not space), can be simplified a bit by eliminating $v^n$ from the first equation and @@ -306,6 +310,10 @@ boundary values for the first equation above (i.e. the one for $U^n$), whereas we do not have to do that for the second one. +

Who are Courant, Friedrichs, and Levy?

+ +CFL condition +

How the program works

Given the above formulation, .... diff --git a/deal.II/examples/step-23/step-23.cc b/deal.II/examples/step-23/step-23.cc index 4afa1129b8..c7a9739b69 100644 --- a/deal.II/examples/step-23/step-23.cc +++ b/deal.II/examples/step-23/step-23.cc @@ -1,10 +1,10 @@ /* $Id: step-4.cc,v 1.34 2006/02/06 21:33:10 wolf Exp $ */ -/* Author: Wolfgang Bangerth, University of Heidelberg, 1999 */ +/* Author: Wolfgang Bangerth, Texas A&M University, 2006 */ /* $Id: step-4.cc,v 1.34 2006/02/06 21:33:10 wolf Exp $ */ /* Version: $Name: $ */ /* */ -/* Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors */ +/* Copyright (C) 2006 by the deal.II authors */ /* */ /* This file is subject to QPL and may not be distributed */ /* without copyright and license information. Please refer */ @@ -12,76 +12,195 @@ /* further information on this license. */ -#include -#include -#include -#include -#include -#include -#include -#include -#include + // @sect3{Include files} + + // We start with the usual assortment + // of include files that we've seen + // in so many of the previous tests: #include #include -#include -#include +#include + #include #include #include #include #include + +#include +#include +#include +#include + +#include +#include +#include #include -#include -#include +#include +#include #include + #include #include #include -#include - + // Here are the only three include + // files of some new interest: The + // first one is already used, for + // example, for the + // VectorTools::interpolate_boundary_values + // and + // VectorTools::apply_boundary_values + // functions. However, we here use + // another function in that class, + // VectorTools::project to compute + // our initial values as the $L^2$ + // projection of the continuous + // initial values. Furthermore, we + // use + // VectorTools::create_right_hand_side + // to generate the integrals + // $(f^n,\phi^n_i)$. These were + // previously always generated by + // hand in + // assemble_system or + // similar functions in application + // code. However, we're too lazy to + // do that here, so simply use a + // library function: +#include + // In a very similar vein, we are + // also too lazy to write the code to + // assemble mass and Laplace + // matrices, although it would have + // only taken copying the relevant + // code from any number of previous + // tutorial programs. Rather, we want + // to focus on the things that are + // truly new to this program and + // therefore use the + // MatrixTools::create_mass_matrix + // and + // MatrixTools::create_laplace_matrix + // functions. They are declared here: +#include + // Finally, here is an include file + // that contains all sorts of tool + // functions that one sometimes + // needs. In particular, we need the + // Utilities::int_to_string class + // that, given an integer argument, + // returns a string representation of + // it. It is particularly useful + // since it allows for a second + // parameter indicating the number of + // digits to which we want the result + // padded with leading zeros. We will + // use this to write output files + // that have the form + // solution-XXX.gnuplot + // where XXX denotes the + // number of the time step and always + // consists of three digits even if + // we are still in the single or + // double digit time steps. +#include + + + // @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, + // $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 we have to + // solve in each time step. These + // will be solved in the two + // functions solve_u and + // solve_v. + // + // Finally, the variable + // theta is used to + // indicate the parameter $\theta$ + // that is used to define which time + // stepping scheme to use. The rest + // is self-explanatory. template -class WaveEquationProblem +class WaveEquation { public: - WaveEquationProblem (); + WaveEquation (); void run (); private: - void make_grid_and_dofs (); - void assemble_system (); + void setup_system (); void solve_u (); void solve_v (); - void output_results (const unsigned int timestep_number) const; + void output_results () const; Triangulation triangulation; FE_Q fe; DoFHandler dof_handler; + ConstraintMatrix constraints; + SparsityPattern sparsity_pattern; SparseMatrix system_matrix; SparseMatrix mass_matrix; SparseMatrix laplace_matrix; - double time_step; - double theta; - Vector solution_u, solution_v; Vector old_solution_u, old_solution_v; Vector system_rhs; + + double time, time_step; + unsigned int timestep_number; + const double theta; }; + // @sect3{Equation data} + + // Before we go on filling in the + // details of the main class, let us + // define the equation data + // corresponding to the problem, + // i.e. initial and boundary values + // as well as a right hand side + // class. We do so using classes + // derived from the Function class + // template that has been used many + // times before, so the following + // should not be a surprise. + // + // Let's start with initial values + // and choose zero for both the value + // $u$ as well as its time + // derivative, the velocity $v$: template -class RightHandSide : public Function +class InitialValuesU : public Function { public: - RightHandSide () : Function() {}; + InitialValuesU () : Function() {}; virtual double value (const Point &p, const unsigned int component = 0) const; @@ -89,10 +208,10 @@ class RightHandSide : public Function template -class InitialValuesU : public Function +class InitialValuesV : public Function { public: - InitialValuesU () : Function() {}; + InitialValuesV () : Function() {}; virtual double value (const Point &p, const unsigned int component = 0) const; @@ -101,10 +220,33 @@ class InitialValuesU : public Function template -class BoundaryValues : public Function +double InitialValuesU::value (const Point &/*p*/, + const unsigned int component) const +{ + Assert (component == 0, ExcInternalError()); + return 0; +} + + + +template +double InitialValuesV::value (const Point &/*p*/, + const unsigned int component) const +{ + Assert (component == 0, ExcInternalError()); + return 0; +} + + + + // Secondly, we have the right hand + // side forcing term. Boring as we + // are, we choose zero here as well: +template +class RightHandSide : public Function { public: - BoundaryValues () : Function() {}; + RightHandSide () : Function() {}; virtual double value (const Point &p, const unsigned int component = 0) const; @@ -112,45 +254,55 @@ class BoundaryValues : public Function - template -double RightHandSide::value (const Point &/*p*/, - const unsigned int /*component*/) const +double RightHandSide::value (const Point &/*p*/, + const unsigned int component) const { -// if (get_time() <= 0.25) -// if ((p[0] <=0) && (p[1] <= 0)) -// return 1; - + Assert (component == 0, ExcInternalError()); return 0; } + + // Finally, we have boundary + // values. They are as described in + // the introduction: template -double InitialValuesU::value (const Point &p, - const unsigned int /*component*/) const +class BoundaryValues : public Function { - // return std::sqrt(p.square()) * std::exp (-p.square()) / 3; - if ((p[0] <=0) && (p[1] <= 0)) - return 1; + public: + BoundaryValues () : Function() {}; + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + - return 0;} template -double BoundaryValues::value (const Point &/*p*/, - const unsigned int /*component*/) const +double BoundaryValues::value (const Point &p, + const unsigned int component) const { - return 0; + Assert (component == 0, ExcInternalError()); + + if ((this->get_time() <= 1) && + (p[0] < 1) && + (p[1] < 1./3) && + (p[1] > -1./3)) + return std::sin (this->get_time() * 2 * deal_II_numbers::PI); + else + return 0; } - + // @sect3{Implementation of the WaveEquation class} template -WaveEquationProblem::WaveEquationProblem () : +WaveEquation::WaveEquation () : fe (1), dof_handler (triangulation), time_step (1./64), @@ -160,17 +312,18 @@ WaveEquationProblem::WaveEquationProblem () : template -void WaveEquationProblem::make_grid_and_dofs () +void WaveEquation::setup_system () { GridGenerator::hyper_cube (triangulation, -1, 1); triangulation.refine_global (7); - std::cout << " Number of active cells: " + std::cout << "Number of active cells: " << triangulation.n_active_cells() << std::endl - << " Total number of cells: " + << "Total number of cells: " << triangulation.n_cells() - << std::endl; + << std::endl + << std::endl; dof_handler.distribute_dofs (fe); @@ -188,18 +341,6 @@ void WaveEquationProblem::make_grid_and_dofs () mass_matrix.reinit (sparsity_pattern); laplace_matrix.reinit (sparsity_pattern); - solution_u.reinit (dof_handler.n_dofs()); - solution_v.reinit (dof_handler.n_dofs()); - old_solution_u.reinit (dof_handler.n_dofs()); - old_solution_v.reinit (dof_handler.n_dofs()); - system_rhs.reinit (dof_handler.n_dofs()); -} - - - -template -void WaveEquationProblem::assemble_system () -{ MatrixCreator::create_mass_matrix (dof_handler, QGauss(3), mass_matrix); MatrixCreator::create_laplace_matrix (dof_handler, QGauss(3), @@ -207,12 +348,20 @@ void WaveEquationProblem::assemble_system () system_matrix.copy_from (mass_matrix); system_matrix.add (theta * theta * time_step * time_step, laplace_matrix); + + solution_u.reinit (dof_handler.n_dofs()); + solution_v.reinit (dof_handler.n_dofs()); + old_solution_u.reinit (dof_handler.n_dofs()); + old_solution_v.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); + + constraints.close (); } template -void WaveEquationProblem::solve_u () +void WaveEquation::solve_u () { SolverControl solver_control (1000, 1e-8*system_rhs.l2_norm()); SolverCG<> cg (solver_control); @@ -226,7 +375,7 @@ void WaveEquationProblem::solve_u () template -void WaveEquationProblem::solve_v () +void WaveEquation::solve_v () { SolverControl solver_control (1000, 1e-8*system_rhs.l2_norm()); SolverCG<> cg (solver_control); @@ -241,7 +390,7 @@ void WaveEquationProblem::solve_v () template -void WaveEquationProblem::output_results (const unsigned int timestep_number) const +void WaveEquation::output_results () const { DataOut data_out; @@ -253,7 +402,7 @@ void WaveEquationProblem::output_results (const unsigned int timestep_numbe std::ostringstream filename; filename << "solution-" - << timestep_number + << Utilities::int_to_string (timestep_number, 3) << ".gnuplot"; std::ofstream output (filename.str().c_str()); data_out.write_gnuplot (output); @@ -263,24 +412,18 @@ void WaveEquationProblem::output_results (const unsigned int timestep_numbe template -void WaveEquationProblem::run () +void WaveEquation::run () { - std::cout << "Solving problem in " << dim << " space dimensions." << std::endl; - - make_grid_and_dofs(); - assemble_system (); + setup_system(); - ConstraintMatrix constraints; - constraints.close(); VectorTools::project (dof_handler, constraints, QGauss(3), InitialValuesU(), old_solution_u); VectorTools::project (dof_handler, constraints, QGauss(3), - ZeroFunction(), + InitialValuesV(), old_solution_v); - - unsigned int timestep_number = 1; - for (double time = time_step; time<=5; time+=time_step, ++timestep_number) + + for (timestep_number=1, time=time_step; time<=5; time+=time_step, ++timestep_number) { std::cout << "Time step " << timestep_number << " at t=" << time @@ -308,10 +451,13 @@ void WaveEquationProblem::run () system_rhs.add (theta * (1-theta) * time_step * time_step, tmp); + BoundaryValues boundary_values_function; + boundary_values_function.set_time (time); + std::map boundary_values; VectorTools::interpolate_boundary_values (dof_handler, 0, - BoundaryValues(), + boundary_values_function, boundary_values); MatrixTools::apply_boundary_values (boundary_values, system_matrix, @@ -341,7 +487,7 @@ void WaveEquationProblem::run () solve_v (); - output_results (timestep_number); + output_results (); old_solution_u = solution_u; old_solution_v = solution_v; @@ -354,8 +500,8 @@ int main () { deallog.depth_console (0); { - WaveEquationProblem<2> wave_equation_problem_2d; - wave_equation_problem_2d.run (); + WaveEquation<2> wave_equation_solver; + wave_equation_solver.run (); } return 0; -- 2.39.5