Added first try of nonlinear solvers.

AMDiS/CMakeLists.txt

@@ -145,6 +145,8 @@ SET(AMDIS_SRC ${SOURCE_DIR}/AdaptBase.cc
 	      ${SOURCE_DIR}/io/ValueReader.cc
 	      ${SOURCE_DIR}/io/ValueWriter.cc
 	      ${SOURCE_DIR}/io/VtkWriter.cc
+	      ${SOURCE_DIR}/nonlin/NonLinUpdater.cc
+	      ${SOURCE_DIR}/nonlin/ProblemNonLin.cc
 	      ${SOURCE_DIR}/parallel/InteriorBoundary.cc
 	      ${SOURCE_DIR}/time/RosenbrockAdaptInstationary.cc
 	      ${SOURCE_DIR}/time/RosenbrockMethod.cc

AMDiS/src/OEMSolver.h

@@ -98,9 +98,9 @@ namespace AMDiS {
 			    SystemVector& b,
 			    VectorialMapper&) = 0;
 
-    inline int solveSystem(const SolverMatrix< Matrix< DOFMatrix* > >& A,
-		    SystemVector& x,
-		    SystemVector& b) 
+    inline int solveSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
+			   SystemVector& x,
+			   SystemVector& b) 
     {
       VectorialMapper mapper(A);
       return solveSystem(A, x, b, mapper);

AMDiS/src/nonlin/NonLinSolver.h

+// ============================================================================
+// ==                                                                        ==
+// == AMDiS - Adaptive multidimensional simulations                          ==
+// ==                                                                        ==
+// ==  http://www.amdis-fem.org                                              ==
+// ==                                                                        ==
+// ============================================================================
+//
+// Software License for AMDiS
+//
+// Copyright (c) 2010 Dresden University of Technology 
+// All rights reserved.
+// Authors: Simon Vey, Thomas Witkowski et al.
+//
+// This file is part of AMDiS
+//
+// See also license.opensource.txt in the distribution.
+
+
+
+/** \file NonLinSolver.h */
+
+#ifndef AMDIS_NONLINSOLVER_H
+#define AMDIS_NONLINSOLVER_H
+
+#include <string>
+
+#include "Global.h"
+#include "CreatorInterface.h"
+#include "MatrixVector.h"
+#include "DOFMatrix.h"
+#include "OEMSolver.h"
+
+namespace AMDiS {
+
+  /** 
+   * \ingroup Solver
+   *
+   * \brief
+   * Pure virtual base class for Newton, NewtonS and Banach which all
+   * solves non linear equation systems. Sub classes must implement the methods
+   * \ref init, \ref exit and \ref nlsolve.
+   */
+  class NonLinSolver
+  {
+  public:
+    /** \brief
+     * constructor.
+     * \param name name of this solver
+     */
+    NonLinSolver(const std::string &name_, 
+		 OEMSolver *solver,
+		 NonLinUpdater *updater)
+      : name(name_),
+	linSolver(solver),
+	nonLinUpdater(updater),
+	tolerance(1.e-4),
+	maxIter(50),
+	info(8),
+	initial_residual(0.0),
+	residual(0.0),
+	usedNorm(NO_NORM)
+    {
+      Parameters::get(name + "->tolerance", tolerance);
+      Parameters::get(name + "->max iteration", maxIter);
+      Parameters::get(name + "->info", info);
+      Parameters::get(name + "->norm", usedNorm);  
+    }
+
+    /** \brief
+     * destructor
+     */
+    virtual ~NonLinSolver() {}
+
+    /** \brief
+     * solves the non linear system. Uses sub class methods
+     */
+    inline int solve(SolverMatrix<Matrix<DOFMatrix*> >& mat,
+		     SystemVector &x, SystemVector &rhs) 
+    {
+      init();
+      int result = nlsolve(mat, x, rhs);
+      exit();
+      return result;
+    }
+
+    inline void setTolerance(double tol)
+    {
+      tolerance = tol;
+    }
+
+    inline double getTolerance()
+    {
+      return tolerance;
+    }
+
+    inline OEMSolver* getLinearSolver()
+    { 
+      return linSolver; 
+    }
+
+    inline void setNonLinUpdater(NonLinUpdater *up) 
+    {
+      nonLinUpdater = up;
+    }
+
+    inline NonLinUpdater *getNonLinUpdater()
+    {
+      return nonLinUpdater;
+    }
+
+  protected:
+    /** \brief
+     * Allocates needed memory. Must be overriden in sub classes.
+     */  
+    virtual void init() = 0;
+
+    /** \brief
+     * Solves the non linear system. Must be overriden in sub classes.
+     */  
+    virtual int nlsolve(SolverMatrix<Matrix<DOFMatrix*> >& matVec,
+			SystemVector& x, SystemVector& rhs) = 0;
+
+    /** \brief
+     * Frees needed memory. Must be overriden in sub classes.
+     */  
+    virtual void exit() = 0;
+  
+    virtual int solveLinearSystem(SolverMatrix<Matrix<DOFMatrix*> >& mat,
+				  SystemVector &fh, SystemVector &x)
+    {
+      FUNCNAME("NonLinSolver::solveLinearSystem()");
+      TEST_EXIT(linSolver)("no solver\n");
+      return linSolver->solveSystem(mat, x, fh);
+    }
+
+  protected:
+    std::string        name;             /**< \brief name of the solver */
+    OEMSolver *linSolver;        /**< \brief linear solver*/
+    NonLinUpdater *nonLinUpdater;    /**< \brief non linear updater */
+    double             tolerance;        /**< \brief solver tolerance */
+    int                maxIter;          /**< \brief maximal # of iterations */
+    int                info;             /**< \brief info level */
+    double             initial_residual; /**< \brief initial residual */
+    double             residual;         /**< \brief current residual */
+    Norm               usedNorm;         /**< \brief used norm */
+  };
+
+  // ============================================================================
+  // ===== class NonLinSolverCreator ============================================
+  // ============================================================================
+
+  /** \brief
+   * Interface for creators of concrete NonLinSolvers. 
+   */
+  class NonLinSolverCreator : public CreatorInterface<NonLinSolver>
+  { 
+  public:
+    virtual ~NonLinSolverCreator() {}
+
+    void setName(std::string name_) 
+    { 
+      name = name_; 
+    }
+
+    void setLinearSolver(OEMSolver *solver) 
+    { 
+      linearSolver = solver; 
+    }
+
+    void setNonLinUpdater(NonLinUpdater *updater) 
+    {
+      nonLinUpdater = updater;
+    }
+
+  protected:
+    std::string name;
+    OEMSolver *linearSolver;
+    NonLinUpdater *nonLinUpdater;
+  };
+
+}
+
+#endif // AMDIS_NONLINSOLVER_H
+

AMDiS/src/nonlin/NonLinUpdater.cc

+//
+// Software License for AMDiS
+//
+// Copyright (c) 2010 Dresden University of Technology 
+// All rights reserved.
+// Authors: Simon Vey, Thomas Witkowski et al.
+//
+// This file is part of AMDiS
+//
+// See also license.opensource.txt in the distribution.
+
+
+#include "NonLinUpdater.h"
+#include "Global.h"
+#include "DOFMatrix.h"
+#include "DOFVector.h"
+#include "Traverse.h"
+#include "ElInfo.h"
+#include "Flag.h"
+#include "Operator.h"
+#include "FiniteElemSpace.h"
+#include "BasisFunction.h"
+#include "BoundaryManager.h"
+
+namespace AMDiS {
+
+  void NonLinUpdater::update(bool updateDF, SystemVector *F)
+  {
+    FUNCNAME("NonLinUpdater::update()");
+
+    TraverseStack stack;
+    ElInfo *el_info;
+    Flag fill_flag;
+
+    TEST_EXIT(F || df)("neither F nor df are set\n");
+
+    if (!updateDF && (F == NULL)) {
+      WARNING("No RHS vector and no update for system matrix! Updater does nothing!\n");
+      return;
+    }
+
+    int size = df->getNumRows();
+
+    const FiniteElemSpace *feSpace = 
+      F ? 
+      F->getDOFVector(0)->getFeSpace() : 
+      (*df)[0][0]->getFeSpace();
+
+    if (updateDF) {
+      TEST_EXIT(df)("df not set but update tried!\n");
+
+      for (int i = 0; i < size; i++) {
+	for (int j = 0; j < size; j++) {
+	  if ((*df)[i][j]) {
+	    (*df)[i][j]->clear();
+	  }
+	}
+      }
+    }
+
+    BasisFunction *basFcts = const_cast<BasisFunction*>(feSpace->getBasisFcts());
+
+    if (F) {
+      for (int i = 0; i < size; i++) {
+	F->getDOFVector(i)->set(0.0);
+      }
+    }
+
+    fill_flag = 
+      Mesh::CALL_LEAF_EL|
+      Mesh::FILL_COORDS|
+      Mesh::FILL_BOUND|
+      Mesh::FILL_DET|
+      Mesh::FILL_GRD_LAMBDA;
+
+    BoundaryType *bound = new BoundaryType[basFcts->getNumber()];
+
+    el_info = stack.traverseFirst(feSpace->getMesh(), -1, fill_flag);
+    while (el_info) {
+      basFcts->getBound(el_info, bound);
+
+      if (updateDF) {
+	for (int i = 0; i < size; i++) {
+	  for (int j = 0; j < size; j++) {
+	    if ((*df)[i][j]) {
+	      (*df)[i][j]->assemble(1.0, el_info, bound);
+	    }
+	  }
+	}
+      }
+
+      if (F) {
+	for(int i = 0; i < size; i++) {
+	  F->getDOFVector(i)->assemble(1.0, el_info, bound);
+	}
+      }
+    
+      el_info = stack.traverseNext(el_info);
+    }
+
+    delete [] bound;
+  }
+
+}

AMDiS/src/nonlin/NonLinUpdater.h

+// ============================================================================
+// ==                                                                        ==
+// == AMDiS - Adaptive multidimensional simulations                          ==
+// ==                                                                        ==
+// ==  http://www.amdis-fem.org                                              ==
+// ==                                                                        ==
+// ============================================================================
+//
+// Software License for AMDiS
+//
+// Copyright (c) 2010 Dresden University of Technology 
+// All rights reserved.
+// Authors: Simon Vey, Thomas Witkowski et al.
+//
+// This file is part of AMDiS
+//
+// See also license.opensource.txt in the distribution.
+
+
+
+/** \file NonLinUpdater.h */
+
+#ifndef AMDIS_NONLINAPDATOR_H
+#define AMDIS_NONLINAPDATOR_H
+
+#include "SystemVector.h"
+
+namespace AMDiS {
+
+  /** \brief
+   * Used in non linear solvers for matrix and vector assemblage.
+   */
+  class NonLinUpdater
+  {
+  public:
+    NonLinUpdater(Matrix<DOFMatrix*> *m) 
+      : df(m) 
+    {}
+
+    ~NonLinUpdater() 
+    {}
+
+    void update(bool updateDF, SystemVector *F);
+
+  protected:
+    /// Matrix to be updated.
+    Matrix<DOFMatrix*> *df;
+  };
+
+}
+
+#endif

AMDiS/src/nonlin/ProblemNonLin.cc

+//
+// Software License for AMDiS
+//
+// Copyright (c) 2010 Dresden University of Technology 
+// All rights reserved.
+// Authors: Simon Vey, Thomas Witkowski et al.
+//
+// This file is part of AMDiS
+//
+// See also license.opensource.txt in the distribution.
+
+
+#include "ProblemNonLin.h"
+#include "NonLinSolver.h"
+#include "CreatorMap.h"
+#include "NonLinUpdater.h"
+#include "Traverse.h"
+#include "AdaptInfo.h"
+
+namespace AMDiS {
+
+  void ProblemNonLin::initialize(Flag initFlag,
+				 ProblemStat *adoptProblem /*= NULL*/,
+				 Flag adoptFlag /*= INIT_NOTHING*/) 
+  {
+    ProblemStat::initialize(initFlag, adoptProblem, adoptFlag);
+
+    // === create Updater ===
+    if (updater_) { 
+      WARNING("updater already created\n");
+    } else {
+      if (initFlag.isSet(INIT_UPDATER) || 
+	  ((!adoptFlag.isSet(INIT_UPDATER)) && initFlag.isSet(INIT_NONLIN_SOLVER))) {
+	createUpdater();
+      } 
+
+      if (adoptProblem && 
+	  (adoptFlag.isSet(INIT_UPDATER) ||  
+	   ((!initFlag.isSet(INIT_UPDATER))&& adoptFlag.isSet(INIT_NONLIN_SOLVER)))) {
+	TEST_EXIT(updater_==NULL)("updater already created\n");
+	updater_ = dynamic_cast<ProblemNonLin*>(adoptProblem)->getUpdater();
+      }
+    }
+
+    if (updater_ == NULL) 
+      WARNING("no updater created\n");
+
+    // === create nonlinear solver ===
+    if (nonLinSolver_) { 
+      WARNING("nonlinear solver already created\n");
+    } else {
+      if (initFlag.isSet(INIT_NONLIN_SOLVER))
+	createNonLinSolver();     
+
+      if (adoptProblem && adoptFlag.isSet(INIT_NONLIN_SOLVER)) {
+	TEST_EXIT(nonLinSolver_ == NULL)("nonlinear solver already created\n");
+	nonLinSolver_ = dynamic_cast<ProblemNonLin*>(adoptProblem)->getNonLinSolver();
+      }
+    }
+
+    if (nonLinSolver_ == NULL)
+      WARNING("no nonlinear solver created\n");
+  }
+
+
+  void ProblemNonLin::createNonLinSolver() 
+  {
+    // create non-linear solver
+    std::string nonLinSolverType("no");
+    Parameters::get(name + "->nonlin solver", nonLinSolverType);
+
+    NonLinSolverCreator *nonLinSolverCreator = 
+      dynamic_cast<NonLinSolverCreator*>(CreatorMap<NonLinSolver>::getCreator(nonLinSolverType));
+
+    nonLinSolverCreator->setLinearSolver(solver);
+    nonLinSolverCreator->setName(name + "->nonlin solver");
+    nonLinSolverCreator->setNonLinUpdater(updater_);
+    nonLinSolver_ = nonLinSolverCreator->create();
+  }
+
+
+  void ProblemNonLin::solve(AdaptInfo *adaptInfo) 
+  {
+    TEST_EXIT(nonLinSolver_)("no non-linear solver!\n");
+    nonLinSolver_->solve(solverMatrix, *solution, *rhs);
+  }
+
+
+  void ProblemNonLin::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag, bool, bool)
+  {
+    FUNCNAME("ProblemNonLin::buildAfterCoarsen()");
+    int nMeshes = static_cast<int>(meshes.size());
+
+    for (int i = 0; i < nMeshes; i++)
+      meshes[i]->dofCompress();    
+
+    for (int i = 0; i < nComponents; i++) {
+      MSG("%d DOFs for %s\n", 
+	  componentSpaces[i]->getAdmin()->getUsedSize(), 
+	  componentSpaces[i]->getName().c_str());
+    }
+    
+    TraverseStack stack;
+    ElInfo *elInfo;
+    
+    // for all elements ...
+    for (int i = 0; i < nComponents; i++) {
+      elInfo = stack.traverseFirst(componentMeshes[i], -1, 
+				   Mesh::CALL_LEAF_EL | 
+				   Mesh::FILL_BOUND |
+				   Mesh::FILL_COORDS |
+				   Mesh::FILL_DET |
+				   Mesh::FILL_GRD_LAMBDA);
+      while (elInfo) {
+	if (solution->getDOFVector(i)->getBoundaryManager()) {
+	  solution->getDOFVector(i)->
+	    getBoundaryManager()->fillBoundaryConditions(elInfo, solution->getDOFVector(i));
+	}
+	elInfo = stack.traverseNext(elInfo);
+      }
+    }
+  }
+
+}

AMDiS/src/nonlin/ProblemNonLin.h

+// ============================================================================
+// ==                                                                        ==
+// == AMDiS - Adaptive multidimensional simulations                          ==
+// ==                                                                        ==
+// ==  http://www.amdis-fem.org                                              ==
+// ==                                                                        ==
+// ============================================================================
+//
+// Software License for AMDiS
+//
+// Copyright (c) 2010 Dresden University of Technology 
+// All rights reserved.
+// Authors: Simon Vey, Thomas Witkowski et al.
+//
+// This file is part of AMDiS
+//
+// See also license.opensource.txt in the distribution.
+
+
+
+/** \file ProblemNonLin.h */
+
+#ifndef AMDIS_PROBLEMNONLIN_H_
+#define AMDIS_PROBLEMNONLIN_H_
+
+#include "ProblemStat.h"
+#include "NonLinUpdater.h"
+#include "NonLinSolver.h"
+#include "DOFVector.h"
+#include "SystemVector.h"
+#include "MatrixVector.h"
+
+namespace AMDiS {
+
+  /**
+   * \ingroup Problem
+   * 
+   * \brief
+   * Standard implementation for a vector valued non linear problem.
+   */
+  class ProblemNonLin : public ProblemStat
+  {
+  public:
+    /// Constructs a ProblemNonLin with given name.
+    ProblemNonLin(const std::string& name_)  
+      : ProblemStat(name_.c_str()),
+	nonLinSolver_(NULL),
+	updater_(NULL)
+    {
+      u0_.resize(nComponents);
+      u0_.set(NULL);
+    }
+
+    /** \brief
+     * Sets \ref u0_ and interpolates it to \ref solution_.
+     */
+    inline void setU0(AbstractFunction<double, WorldVector<double> > *u0Fct,
+		      int index) 
+    {
+      FUNCNAME("ProblemNonLinVec::setU0()");
+      TEST_EXIT(index < nComponents)("invalid index\n");
+      u0_[index] = u0Fct;
+      solution->getDOFVector(index)->interpol(u0Fct);
+    }
+
+    /** \brief
+     * Destructor.
+     */
+    virtual ~ProblemNonLin() {}
+
+    /** \brief
+     * Initialisation of the problem.
+     */
+    virtual void initialize(Flag initFlag,
+			    ProblemStat *adoptProblem = NULL,
+			    Flag adoptFlag = INIT_NOTHING);
+
+    /** \brief
+     * Used in \ref initialize().
+     */
+    virtual void createUpdater() 
+    {
+      updater_ = new NonLinUpdater(systemMatrix);
+    }
+
+    /** \brief
+     * Used in \ref initialize().
+     */
+    virtual void createNonLinSolver();
+
+    /** \brief
+     * Overrides ProblemVec::solve(). Uses the non linear solver 
+     * \ref nonLinSolver_.
+     */
+    virtual void solve(AdaptInfo *adaptInfo);
+
+    /** \brief
+     * Overrides Problem::buildAfterCoarsen(). Sets dirichlet boundaries
+     * in \ref A_ and \ref rhs_.
+     */
+    virtual void buildAfterCoarsen(AdaptInfo *adaptInfo, Flag,
+				   bool assembleMatrix = true,
+				   bool assembleVector = true);
+
+    /** \brief
+     * Returns /ref nonLinSolver_.
+     */
+    inline NonLinSolver *getNonLinSolver() 
+    {
+      return nonLinSolver_;
+    }
+
+    /** \brief
+     * Returns \ref updater.
+     */
+    inline NonLinUpdater* getUpdater() 
+    {
+      return updater_;
+    }
+
+    /** \brief
+     * Sets \ref nonLinSolver_.
+     */
+    inline void setNonLinSolver(NonLinSolver *nonLinSolver) 
+    {
+      nonLinSolver_ = nonLinSolver;
+    }