Add ObstacleTNNMGStep class.

This implements the TNNMG method for quadratic obstacle problems.
A much more general implementation can be found in the dune-tnnmg
module (please contact C. Gräser or O. Sander).

[[Imported from SVN: r5826]]

dune/solvers/iterationsteps/Makefile.am

@@ -12,6 +12,7 @@ iterationsteps_HEADERS = amgstep.hh \
                          mmgstep.hh mmgstep.cc \
                          multigridstep.hh \
                          multigridstep.cc \
+                         obstacletnnmgstep.hh \
                          projectedblockgsstep.hh \
                          projectedblockgsstep.cc \
                          projectedlinegsstep.cc projectedlinegsstep.hh \

dune/solvers/iterationsteps/obstacletnnmgstep.hh

+#ifndef OBSTACLE_TNNMG_STEP_HH
+#define OBSTACLE_TNNMG_STEP_HH
+
+// std includes
+#include <vector>
+#include <algorithm>
+
+// dune-common includes
+#include <dune/common/bitsetvector.hh>
+
+// dune-solver includes
+#include <dune/solvers/iterationsteps/iterationstep.hh>
+#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
+#include <dune/solvers/iterationsteps/projectedblockgsstep.hh>
+#include <dune/solvers/iterationsteps/multigridstep.hh>
+
+#include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/common/boxconstraint.hh>
+
+
+
+
+
+/**
+ * \brief TNNMG step for quadratic obstacle problems
+ *
+ * This is a dune-solvers iteration step implementing
+ * the TNNMG (truncated nonsmooth Newton multigrid) method
+ * for convex quadratic obstacle problems for use in the
+ * LoopSolver class.
+ *
+ * The method basically consists of a projected Gauss-Seidel
+ * method as nonlinear pre- and post-smoother and a linear
+ * multigrid method for a truncated linear system. For
+ * details please refer to
+ *
+ *   'Multigrid Methods for Obstacle Problems', C. Gräser and R. Kornhuber
+ *   Journal of Computational Mathematics 27 (1), 2009, pp. 1-44 
+ *
+ * A much more general and flexible implementation that can also be used for
+ * energy functionals with nonquadratic (anisotropic) smooth part
+ * and other nonsmooth part, e.g. incorporating local simplex
+ * constraints, can be found in the dune-tnnmg module.
+ *
+ * \tparam MatrixType Type for the matrix describing the quadratic part
+ * \tparam VectorType Type for solution and rhs vectors
+ */
+template<class MatrixType, class VectorType>
+class ObstacleTNNMGStep
+    : public IterationStep<VectorType>
+{
+        typedef IterationStep<VectorType> Base;
+        typedef typename VectorType::field_type field_type;
+
+        enum {blockSize = VectorType::block_type::dimension};
+    public:
+
+        typedef VectorType Vector;
+        typedef MatrixType Matrix;
+        typedef BoxConstraint<field_type, blockSize> BlockBoxConstraint;
+        typedef typename std::vector<BlockBoxConstraint> BoxConstraintVector;
+        typedef typename Dune::BitSetVector<VectorType::block_type::dimension> BitVector;
+
+        typedef MultigridTransfer<Vector, BitVector, Matrix> Transfer;
+        typedef typename std::vector<Transfer*> TransferPointerVector;
+
+    protected:
+
+        typedef ProjectedBlockGSStep<Matrix, Vector> NonlinearSmoother;
+        typedef TruncatedBlockGSStep<Matrix, Vector> LinearSmoother;
+        typedef MultigridStep<Matrix, Vector, BitVector> LinearMultigridStep;
+        typedef EnergyNorm<Matrix, Vector> LinearBaseSolverEnergyNorm;
+        typedef LoopSolver<Vector> LinearBaseSolver;
+
+    public:
+
+        /**
+         * \brief Construct iteration step
+         *
+         * \param mat Matrix describing the quadratic part of the energy
+         * \param x Solution vector containing the initial iterate
+         * \param rhs Vector containing the linear part of the energy
+         * \param ignore BitVector marking components to ignore (e.g. Dirichlet of hanging nodes)
+         * \param obstacles std::vector containing local obstacles for each global index
+         * \param transfer A vector of pointers to transfer operators describing the subspace hierarchy
+         * \param truncationTol Components will be truncated for the coarse grid correction if they are
+         *                      this close to the obstacle (the default will work well in general).
+         *
+         */
+        ObstacleTNNMGStep(
+                const Matrix& mat,
+                Vector& x,
+                const Vector& rhs,
+                const BitVector& ignore,
+                const BoxConstraintVector& obstacles,
+                const TransferPointerVector& transfer,
+                double truncationTol=1e-10) :
+            Base(x),
+            mat_(mat),
+            rhs_(rhs),
+            obstacles_(obstacles),
+            transfer_(transfer),
+            baseSolverNorm_(baseSolverStep_),
+            baseSolver_(&baseSolverStep_, 100, 1e-8, &baseSolverNorm_, LinearBaseSolver::QUIET),
+            nonlinearSmoother_(mat_, *x_, rhs_),
+            preSmoothingSteps_(3),
+            postSmoothingSteps_(3),
+            linearIterationSteps_(1),
+            linearPreSmoothingSteps_(3),
+            linearPostSmoothingSteps_(3),
+            truncationTol_(truncationTol)
+        {
+            ignoreNodes_ = &ignore;
+        }
+
+
+        /**
+         * \brief Set type of multigrid cycle
+         *
+         * If you don't call this method manually, 3 nonlinear pre- and post-smoothing
+         * steps and a linear V(3,3)-cycle will be used.
+         *
+         * \param preSmoothingSteps Number of nonlinear pre-smoothing steps.
+         * \param linearIterationSteps Number of linear multigrid steps for truncated problem
+         * \param postSmoothingSteps Number of nonlinear post-smoothing steps.
+         * \param linearPreSmoothingSteps Number of linear pre-smoothing steps during the linear multigrid V-cycle
+         * \param linearPostSmoothingSteps Number of linear post-smoothing steps during the linear multigrid V-cycle
+         */
+        void setSmoothingSteps(int preSmoothingSteps, int linearIterationSteps, int postSmoothingSteps, int linearPreSmoothingSteps, int linearPostSmoothingSteps)
+        {
+            preSmoothingSteps_ = preSmoothingSteps;
+            postSmoothingSteps_ = postSmoothingSteps;
+            linearIterationSteps_ = linearIterationSteps;
+            linearPreSmoothingSteps_ = linearPreSmoothingSteps;
+            linearPostSmoothingSteps_ = linearPostSmoothingSteps;
+        }
+
+
+        /**
+         * \brief Setup the iteration step
+         *
+         * This method must be called before the iterate method is called.
+         */
+        virtual void preprocess()
+        {
+            truncatedResidual_.resize(x_->size());
+            residual_.resize(x_->size());
+            coarseCorrection_.resize(x_->size());
+            temp_.resize(x_->size());
+
+            linearMGStep_.setNumberOfLevels(transfer_.size()+1);
+            linearMGStep_.setTransferOperators(transfer_);
+            linearMGStep_.setSmoother(&linearSmoother_);
+            linearMGStep_.basesolver_ = &baseSolver_;
+            linearMGStep_.setMGType(1, linearPreSmoothingSteps_, linearPostSmoothingSteps_);
+
+            hasObstacle_.resize(x_->size());
+            hasObstacle_.setAll();
+
+            nonlinearSmoother_.obstacles_ = &obstacles_;
+            nonlinearSmoother_.hasObstacle_ = &hasObstacle_;
+            nonlinearSmoother_.ignoreNodes_ = ignoreNodes_;
+            outStream_.str("");
+            outStream_ << "      step size";
+            for(int j=0; j<blockSize; ++j)
+                outStream_ << "  trunc" << std::setw(2) << j;
+        }
+
+
+        /**
+         * \brief Apply one iteration step
+         *
+         * You can query the current iterate using the getSol() method
+         * afterwards. The preprocess() method must be called before this
+         * one.
+         */
+        virtual void iterate()
+        {
+            // clear previous output data
+            outStream_.str("");
+
+            // use reference to simplify the following code
+            Vector& x = *x_;
+
+            // nonlinear pre-smoothing
+            for(int i=0; i<preSmoothingSteps_; ++i)
+                nonlinearSmoother_.iterate();
+            x = nonlinearSmoother_.getSol();
+
+            // compute residual
+            residual_ = rhs_;
+            mat_.mmv(x, residual_);
+
+            // initialize matrix, residual, and correction vector
+            // for linear coarse grid correction
+            Matrix truncatedMat = mat_;
+            truncatedResidual_ = residual_;
+            coarseCorrection_ = 0;
+
+            // mark indices for truncation and count truncated indices per component
+            typename Dune::array<int, blockSize> truncationCount;
+            truncationCount.fill(0);
+            BitVector truncate(x.size(), false);
+            for(int i=0; i<x.size(); ++i)
+            {
+                for (int j = 0; j < blockSize; ++j)
+                {
+                    if ((x[i][j]<=obstacles_[i].lower(j)+truncationTol_) or (x[i][j]>=obstacles_[i].upper(j)-truncationTol_))
+                    {
+                        truncate[i][j] = true;
+                        ++truncationCount[j];
+                    }
+                }
+            }
+
+            // truncate matrix and residual
+            typename Matrix::row_type::Iterator it;
+            typename Matrix::row_type::Iterator end;
+            for(int i=0; i<x.size(); ++i)
+            {
+                it = truncatedMat[i].begin();
+                end = truncatedMat[i].end();
+                for(; it!=end; ++it)
+                {
+                    int j = it.index();
+                    typename Matrix::block_type& Aij = *it;
+                    for (int ii = 0; ii < blockSize; ++ii)
+                        for (int jj = 0; jj < blockSize; ++jj)
+                            if (truncate[i][ii] or truncate[j][jj])
+                                Aij[ii][jj] = 0;
+                }
+                for (int ii = 0; ii < blockSize; ++ii)
+                    if (truncate[i][ii])
+                        truncatedResidual_[i][ii] = 0;
+            }
+
+            // apply linear multigrid to approximatively solve the truncated linear system
+            linearMGStep_.setProblem(truncatedMat, coarseCorrection_, truncatedResidual_);
+            linearMGStep_.ignoreNodes_ = ignoreNodes_;
+            linearMGStep_.preprocess();
+            for(int i=0; i<linearIterationSteps_; ++i)
+                linearMGStep_.iterate();
+            coarseCorrection_ = linearMGStep_.getSol();
+
+            // post process multigrid correction
+            // There will be spurious entries in non-truncated
+            // components since we did not truncate the transfer
+            // operators. These can safely be removed.
+            for(int i=0; i<x.size(); ++i)
+                for (int ii = 0; ii < blockSize; ++ii)
+                    if (truncate[i][ii])
+                        coarseCorrection_[i][ii] = 0;
+
+            // project correction into the defect domain
+            // and compute directional constraints for line search
+            BoxConstraint<field_type, 1> directionalConstraints;
+            for(int i=0; i<x.size(); ++i)
+            {
+                // compute local defect constraints
+                BlockBoxConstraint defectConstraint = obstacles_[i];
+                defectConstraint -= x[i];
+
+                // project correction into local defect constraints
+                for (int j = 0; j < blockSize; ++j)
+                {
+                    if (coarseCorrection_[i][j] < defectConstraint.lower(j))
+                        coarseCorrection_[i][j] = defectConstraint.lower(j);
+                    if (coarseCorrection_[i][j] > defectConstraint.upper(j))
+                        coarseCorrection_[i][j] = defectConstraint.upper(j);
+                }
+
+                for (int j = 0; j < blockSize; ++j)
+                {
+                    // compute relative defect constraints
+                    defectConstraint.lower(j) /= coarseCorrection_[i][j];
+                    defectConstraint.upper(j) /= coarseCorrection_[i][j];
+
+                    // orient relative defect constraints
+                    if (defectConstraint.lower(j) > defectConstraint.upper(j))
+                        std::swap(defectConstraint.lower(j), defectConstraint.upper(j));
+
+                    // update directional constraints
+                    directionalConstraints.lower(0) = std::max(directionalConstraints.lower(0), defectConstraint.lower(j));
+                    directionalConstraints.upper(0) = std::min(directionalConstraints.upper(0), defectConstraint.upper(j));
+                }
+            }
+
+            // compute damping parameter
+            mat_.mv(coarseCorrection_, temp_);
+            field_type damping = 0;
+            field_type coarseCorrectionNorm2 = coarseCorrection_*temp_;
+            if (coarseCorrectionNorm2 > 0)
+                damping = (residual_*coarseCorrection_) / coarseCorrectionNorm2;
+            else
+                damping = 0;
+            damping = std::max(damping, directionalConstraints.lower(0));
+            damping = std::min(damping, directionalConstraints.upper(0));
+
+            // apply correction
+            x.axpy(damping, coarseCorrection_);
+
+            // nonlinear post-smoothing
+            for(int i=0; i<postSmoothingSteps_; ++i)
+                nonlinearSmoother_.iterate();
+            x = nonlinearSmoother_.getSol();
+
+            // fill stream with output data
+            outStream_.setf(std::ios::scientific);
+            outStream_ << std::setw(15) << damping;
+            for(int j=0; j<blockSize; ++j)
+                outStream_ << std::setw(9) << truncationCount[j];
+        }
+
+
+        /**
+         * \brief Get the current iterate
+         */
+        VectorType getSol()
+        {
+            return *x_;
+        }
+
+
+        /**
+         * \brief Get output accumulated during last iterate() call.
+         *
+         * This method is used to report TNNMG-specific data to the
+         * LoopSolver class during the iterative solution.
+         */
+        std::string getOutput() const
+        {
+            std::string s = outStream_.str();
+            outStream_.str("");
+            return s;
+        }
+
+    protected:
+
+        // problem description given by the user
+        const Matrix& mat_;
+        using Base::x_;
+        using Base::ignoreNodes_;
+        const Vector& rhs_;
+        const BoxConstraintVector& obstacles_;
+
+        // transfer operators given by the user
+        const TransferPointerVector& transfer_;
+
+        // vectors needed during iteration
+        Vector residual_;
+        Vector temp_;
+        Vector truncatedResidual_;
+        Vector coarseCorrection_;
+
+        // solver components
+        NonlinearSmoother nonlinearSmoother_;
+        typename Dune::BitSetVector<1> hasObstacle_;
+        double truncationTol_;
+
+        LinearMultigridStep linearMGStep_;
+        LinearSmoother linearSmoother_;
+        LinearSmoother baseSolverStep_;
+        LinearBaseSolverEnergyNorm baseSolverNorm_;
+        LinearBaseSolver baseSolver_;
+
+        int preSmoothingSteps_;
+        int postSmoothingSteps_;
+        int linearIterationSteps_;
+        int linearPreSmoothingSteps_;
+        int linearPostSmoothingSteps_;
+
+        mutable std::ostringstream outStream_;
+};
+
+
+
+
+
+
+
+
+#endif