From: Guido Kanschat Date: Mon, 5 Jul 2010 18:14:15 +0000 (+0000) Subject: working on new example X-Git-Tag: v8.0.0~5852 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9c8ac43c0dd25daa60804514998ae6c04c8deeab;p=dealii.git working on new example git-svn-id: https://svn.dealii.org/trunk@21454 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/doxygen/Makefile b/deal.II/examples/doxygen/Makefile index 1ed195f86c..5bd22b1a7c 100644 --- a/deal.II/examples/doxygen/Makefile +++ b/deal.II/examples/doxygen/Makefile @@ -2,7 +2,7 @@ # $Id$ # Version: $Name$ # -# Copyright (C) 1998 - 2007 by the deal.II authors +# Copyright (C) 1998 - 2007, 2010 by the deal.II authors # # This file is subject to QPL and may not be distributed # without copyright and license information. Please refer @@ -13,7 +13,8 @@ include ../../common/Make.global_options all: product_matrix$(EXEEXT) block_matrix_array$(EXEEXT) \ - compressed_block_sparsity_pattern$(EXEEXT) + compressed_block_sparsity_pattern$(EXEEXT) \ + theta_timestepping$(EXEEXT) ###################################################################### # Compilation of source code @@ -42,6 +43,10 @@ compressed_block_sparsity_pattern$(EXEEXT): compressed_block_sparsity_pattern.g. @echo ============================ Linking $@ @$(CXX) -o $@ $^ $(LIBS) $(LDFLAGS) +theta_timestepping$(EXEEXT): theta_timestepping.g.$(OBJEXT) $(lib-deal2-2d.g) $(lib-lac.g) $(lib-base.g) + @echo ============================ Linking $@ + @$(CXX) -o $@ $^ $(LIBS) $(LDFLAGS) + ###################################################################### # Pseudo target for cleaning directory ###################################################################### diff --git a/deal.II/examples/doxygen/theta_timestepping.cc b/deal.II/examples/doxygen/theta_timestepping.cc index 8d9bf626aa..6f42cbd37f 100644 --- a/deal.II/examples/doxygen/theta_timestepping.cc +++ b/deal.II/examples/doxygen/theta_timestepping.cc @@ -22,43 +22,106 @@ using namespace dealii; using namespace Algorithms; + class Explicit - : Operator > + : public Operator > { public: + Explicit(const FullMatrix& matrix); void operator() (NamedData*>& out, const NamedData*>& in); - void initialize_timestep_data(const TimeStepData&); + void initialize_timestep_data(const TimestepData&); private: - const TimeStepData* timestep_data; + const TimestepData* timestep_data; + SmartPointer, Explicit> matrix; + FullMatrix m; }; class Implicit - : Operator > + : public Operator > { public: + Implicit(const FullMatrix& matrix); void operator() (NamedData*>& out, const NamedData*>& in); - void initialize_timestep_data(const TimeStepData&); + void initialize_timestep_data(const TimestepData&); private: - const TimeStepData* timestep_data; + const TimestepData* timestep_data; + SmartPointer, Implicit> matrix; + FullMatrix m; }; int main() { - Explicit op_explicit; - Implicit op_implicit; + FullMatrix matrix(2); + matrix(0,0) = -.1; + matrix(1,1) = -.1; + matrix(0,1) = 31.4; + matrix(1,0) = -31.4; + + Explicit op_explicit(matrix); + Implicit op_implicit(matrix); ThetaTimestepping > solver(op_explicit, op_implicit); op_explicit.initialize_timestep_data(solver.explicit_data()); op_implicit.initialize_timestep_data(solver.implicit_data()); - + solver.notify(Events::initial); + Vector value(2); NamedData*> indata; NamedData*> outdata; - outdata.add(&value, "value"); + Vector* p = &value; + outdata.add(p, "value"); solver(outdata, indata); } + + +Explicit::Explicit(const FullMatrix& M) + : + matrix(&M) +{} + + +void +Explicit::initialize_timestep_data(const TimestepData& t) +{ + timestep_data = &t; +} + + +void +Explicit::operator() (NamedData*>& out, const NamedData*>& in) +{ + this->notifications.print(deallog); + deallog << std::endl; + unsigned int i = in.find("Previous time"); + m.vmult(*out(0), *in(i)); +} + + +Implicit::Implicit(const FullMatrix& M) + : + matrix(&M) +{} + + +void +Implicit::initialize_timestep_data(const TimestepData& t) +{ + timestep_data = &t; +} + + +void +Implicit::operator() (NamedData*>& out, const NamedData*>& in) +{ + this->notifications.print(deallog); + deallog << std::endl; + unsigned int i = in.find("Previous time"); + m.vmult(*out(0), *in(i)); +} + + diff --git a/deal.II/examples/step-20/doc/intro.dox b/deal.II/examples/step-20/doc/intro.dox index 91c7b3ca77..28e6b48b88 100644 --- a/deal.II/examples/step-20/doc/intro.dox +++ b/deal.II/examples/step-20/doc/intro.dox @@ -525,43 +525,6 @@ this looks almost as expensive as solving with $S$ right away. However, note that in the inner iteration, we do not have to calculate $M^{-1}$, but only the inverse of its diagonal, which is cheap. -To implement something like this, let us first generalize the -InverseMatrix class so that it can work not only with -SparseMatrix objects, but with any matrix type. This looks like so: - -@code -template -class InverseMatrix -{ - public: - InverseMatrix (const Matrix &m); - - void vmult (Vector &dst, - const Vector &src) const; - - private: - const SmartPointer matrix; - - //... -}; - - -template -void InverseMatrix::vmult (Vector &dst, - const Vector &src) const -{ - SolverControl solver_control (src.size(), 1e-8*src.l2_norm()); - SolverCG<> cg (solver_control, vector_memory); - - dst = 0; - - cg.solve (*matrix, dst, src, PreconditionIdentity()); -} -@endcode - -Essentially, the only change we have made is the introduction of a template -argument that generalizes the use of SparseMatrix. - The next step is to define a class that represents the approximate Schur complement. This should look very much like the Schur complement class itself, except that it doesn't need the object representing $M^{-1}$ any more: @@ -662,9 +625,7 @@ chose to introduce this material here anyway to demonstrate how to work with block matrices and to develop solvers and preconditioners, rather than using black box components from the library. -For those interested in looking up the corresponding library classes: the -InverseMatrix is roughly equivalent to the -PreconditionLACSolver class in the library. Likewise, the Schur +For those interested in looking up the corresponding library classes: the Schur complement class corresponds to the SchurMatrix class. diff --git a/deal.II/examples/step-20/step-20.cc b/deal.II/examples/step-20/step-20.cc index febeb0602f..6cc260a81b 100644 --- a/deal.II/examples/step-20/step-20.cc +++ b/deal.II/examples/step-20/step-20.cc @@ -34,6 +34,7 @@ // inverse of a matrix by calling an // iterative solver. #include +#include #include #include @@ -726,86 +727,6 @@ void MixedLaplaceProblem::assemble_system () // rather only comment on // implementational aspects. - - // @sect4{The SchurComplement class template} - - // The next class is the Schur - // complement class. Its rationale - // has also been discussed in length - // in the introduction. The only - // things we would like to note is - // that the class, too, is derived - // from the Subscriptor class and - // that as mentioned above it stores - // pointers to the entire block - // matrix and the inverse of the mass - // matrix block using - // SmartPointer objects. - // - // The vmult function requires - // two temporary vectors that we do - // not want to re-allocate and free - // every time we call this - // function. Since here, we have full - // control over the use of these - // vectors (unlike above, where a - // class called by the vmult - // function required these vectors, - // not the vmult function - // itself), we allocate them - // directly, rather than going - // through the VectorMemory - // mechanism. However, again, these - // member variables do not carry any - // state between successive calls to - // the member functions of this class - // (i.e., we never care what values - // they were set to the last time a - // member function was called), we - // mark these vectors as mutable. - // - // The rest of the (short) - // implementation of this class is - // straightforward if you know the - // order of matrix-vector - // multiplications performed by the - // vmult function: -class SchurComplement : public Subscriptor -{ - public: - SchurComplement (const BlockSparseMatrix &A, - const IterativeInverse > &Minv); - - void vmult (Vector &dst, - const Vector &src) const; - - private: - const SmartPointer > system_matrix; - const SmartPointer > > m_inverse; - - mutable Vector tmp1, tmp2; -}; - - -SchurComplement::SchurComplement (const BlockSparseMatrix &A, - const IterativeInverse > &Minv) - : - system_matrix (&A), - m_inverse (&Minv), - tmp1 (A.block(0,0).m()), - tmp2 (A.block(0,0).m()) -{} - - -void SchurComplement::vmult (Vector &dst, - const Vector &src) const -{ - system_matrix->block(0,1).vmult (tmp1, src); - m_inverse->vmult (tmp2, tmp1); - system_matrix->block(1,0).vmult (dst, tmp2); -} - - // @sect4{The ApproximateSchurComplement class template} // The third component of our solver @@ -901,7 +822,11 @@ void MixedLaplaceProblem::solve () m_inverse.solver.select("cg"); ReductionControl inner_control(1000, 0., 1.e-13); m_inverse.solver.control = inner_control; - + + SchurComplement schur_complement (m_inverse, + system_matrix.block(1,0), + system_matrix.block(1,0)); + Vector tmp (solution.block(0).size()); // Now on to the first @@ -925,8 +850,6 @@ void MixedLaplaceProblem::solve () schur_rhs -= system_rhs.block(1); - SchurComplement - schur_complement (system_matrix, m_inverse); ApproximateSchurComplement approximate_schur_complement (system_matrix);