# $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
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
@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
######################################################################
using namespace dealii;
using namespace Algorithms;
+
class Explicit
- : Operator<Vector<double> >
+ : public Operator<Vector<double> >
{
public:
+ Explicit(const FullMatrix<double>& matrix);
void operator() (NamedData<Vector<double>*>& out,
const NamedData<Vector<double>*>& in);
- void initialize_timestep_data(const TimeStepData&);
+ void initialize_timestep_data(const TimestepData&);
private:
- const TimeStepData* timestep_data;
+ const TimestepData* timestep_data;
+ SmartPointer<const FullMatrix<double>, Explicit> matrix;
+ FullMatrix<double> m;
};
class Implicit
- : Operator<Vector<double> >
+ : public Operator<Vector<double> >
{
public:
+ Implicit(const FullMatrix<double>& matrix);
void operator() (NamedData<Vector<double>*>& out,
const NamedData<Vector<double>*>& in);
- void initialize_timestep_data(const TimeStepData&);
+ void initialize_timestep_data(const TimestepData&);
private:
- const TimeStepData* timestep_data;
+ const TimestepData* timestep_data;
+ SmartPointer<const FullMatrix<double>, Implicit> matrix;
+ FullMatrix<double> m;
};
int main()
{
- Explicit op_explicit;
- Implicit op_implicit;
+ FullMatrix<double> 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<Vector<double> > 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<double> value(2);
NamedData<Vector<double>*> indata;
NamedData<Vector<double>*> outdata;
- outdata.add(&value, "value");
+ Vector<double>* p = &value;
+ outdata.add(p, "value");
solver(outdata, indata);
}
+
+
+Explicit::Explicit(const FullMatrix<double>& M)
+ :
+ matrix(&M)
+{}
+
+
+void
+Explicit::initialize_timestep_data(const TimestepData& t)
+{
+ timestep_data = &t;
+}
+
+
+void
+Explicit::operator() (NamedData<Vector<double>*>& out, const NamedData<Vector<double>*>& 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<double>& M)
+ :
+ matrix(&M)
+{}
+
+
+void
+Implicit::initialize_timestep_data(const TimestepData& t)
+{
+ timestep_data = &t;
+}
+
+
+void
+Implicit::operator() (NamedData<Vector<double>*>& out, const NamedData<Vector<double>*>& in)
+{
+ this->notifications.print(deallog);
+ deallog << std::endl;
+ unsigned int i = in.find("Previous time");
+ m.vmult(*out(0), *in(i));
+}
+
+
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
-<code>InverseMatrix</code> class so that it can work not only with
-<code>SparseMatrix</code> objects, but with any matrix type. This looks like so:
-
-@code
-template <class Matrix>
-class InverseMatrix
-{
- public:
- InverseMatrix (const Matrix &m);
-
- void vmult (Vector<double> &dst,
- const Vector<double> &src) const;
-
- private:
- const SmartPointer<const Matrix> matrix;
-
- //...
-};
-
-
-template <class Matrix>
-void InverseMatrix<Matrix>::vmult (Vector<double> &dst,
- const Vector<double> &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 <code>SparseMatrix</code>.
-
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:
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
-<code>InverseMatrix</code> is roughly equivalent to the
-<code>PreconditionLACSolver</code> class in the library. Likewise, the Schur
+For those interested in looking up the corresponding library classes: the Schur
complement class corresponds to the <code>SchurMatrix</code> class.
// inverse of a matrix by calling an
// iterative solver.
#include <lac/iterative_inverse.h>
+#include <lac/schur_matrix.h>
#include <grid/tria.h>
#include <grid/grid_generator.h>
// rather only comment on
// implementational aspects.
-
- // @sect4{The <code>SchurComplement</code> 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 <code>Subscriptor</code> class and
- // that as mentioned above it stores
- // pointers to the entire block
- // matrix and the inverse of the mass
- // matrix block using
- // <code>SmartPointer</code> objects.
- //
- // The <code>vmult</code> 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 <code>vmult</code>
- // function required these vectors,
- // not the <code>vmult</code> function
- // itself), we allocate them
- // directly, rather than going
- // through the <code>VectorMemory</code>
- // 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 <code>mutable</code>.
- //
- // The rest of the (short)
- // implementation of this class is
- // straightforward if you know the
- // order of matrix-vector
- // multiplications performed by the
- // <code>vmult</code> function:
-class SchurComplement : public Subscriptor
-{
- public:
- SchurComplement (const BlockSparseMatrix<double> &A,
- const IterativeInverse<Vector<double> > &Minv);
-
- void vmult (Vector<double> &dst,
- const Vector<double> &src) const;
-
- private:
- const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
- const SmartPointer<const IterativeInverse<Vector<double> > > m_inverse;
-
- mutable Vector<double> tmp1, tmp2;
-};
-
-
-SchurComplement::SchurComplement (const BlockSparseMatrix<double> &A,
- const IterativeInverse<Vector<double> > &Minv)
- :
- system_matrix (&A),
- m_inverse (&Minv),
- tmp1 (A.block(0,0).m()),
- tmp2 (A.block(0,0).m())
-{}
-
-
-void SchurComplement::vmult (Vector<double> &dst,
- const Vector<double> &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 <code>ApproximateSchurComplement</code> class template}
// The third component of our solver
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<double> tmp (solution.block(0).size());
// Now on to the first
schur_rhs -= system_rhs.block(1);
- SchurComplement
- schur_complement (system_matrix, m_inverse);
ApproximateSchurComplement
approximate_schur_complement (system_matrix);