From 757ba5c5fe8c515f63f8b8a198f562739a857cba Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Sun, 26 Jun 2011 09:10:12 +0000
Subject: [PATCH] Added first try of nonlinear solvers.
---
AMDiS/CMakeLists.txt | 2 +
AMDiS/src/OEMSolver.h | 6 +-
AMDiS/src/nonlin/NonLinSolver.h | 185 ++++++++++++++++++++++++++++++
AMDiS/src/nonlin/NonLinUpdater.cc | 104 +++++++++++++++++
AMDiS/src/nonlin/NonLinUpdater.h | 52 +++++++++
AMDiS/src/nonlin/ProblemNonLin.cc | 124 ++++++++++++++++++++
AMDiS/src/nonlin/ProblemNonLin.h | 158 +++++++++++++++++++++++++
7 files changed, 628 insertions(+), 3 deletions(-)
create mode 100644 AMDiS/src/nonlin/NonLinSolver.h
create mode 100644 AMDiS/src/nonlin/NonLinUpdater.cc
create mode 100644 AMDiS/src/nonlin/NonLinUpdater.h
create mode 100644 AMDiS/src/nonlin/ProblemNonLin.cc
create mode 100644 AMDiS/src/nonlin/ProblemNonLin.h
diff --git a/AMDiS/CMakeLists.txt b/AMDiS/CMakeLists.txt
index a238f467..16e74b64 100644
--- a/AMDiS/CMakeLists.txt
+++ b/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
diff --git a/AMDiS/src/OEMSolver.h b/AMDiS/src/OEMSolver.h
index 329bad38..37a7df8e 100644
--- a/AMDiS/src/OEMSolver.h
+++ b/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);
diff --git a/AMDiS/src/nonlin/NonLinSolver.h b/AMDiS/src/nonlin/NonLinSolver.h
new file mode 100644
index 00000000..671a1d1c
--- /dev/null
+++ b/AMDiS/src/nonlin/NonLinSolver.h
@@ -0,0 +1,185 @@
+// ============================================================================
+// == ==
+// == 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
+
diff --git a/AMDiS/src/nonlin/NonLinUpdater.cc b/AMDiS/src/nonlin/NonLinUpdater.cc
new file mode 100644
index 00000000..451c70a4
--- /dev/null
+++ b/AMDiS/src/nonlin/NonLinUpdater.cc
@@ -0,0 +1,104 @@
+//
+// 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;
+ }
+
+}
diff --git a/AMDiS/src/nonlin/NonLinUpdater.h b/AMDiS/src/nonlin/NonLinUpdater.h
new file mode 100644
index 00000000..3092845e
--- /dev/null
+++ b/AMDiS/src/nonlin/NonLinUpdater.h
@@ -0,0 +1,52 @@
+// ============================================================================
+// == ==
+// == 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
diff --git a/AMDiS/src/nonlin/ProblemNonLin.cc b/AMDiS/src/nonlin/ProblemNonLin.cc
new file mode 100644
index 00000000..1c66a269
--- /dev/null
+++ b/AMDiS/src/nonlin/ProblemNonLin.cc
@@ -0,0 +1,124 @@
+//
+// 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);
+ }
+ }
+ }
+
+}
diff --git a/AMDiS/src/nonlin/ProblemNonLin.h b/AMDiS/src/nonlin/ProblemNonLin.h
new file mode 100644
index 00000000..ded55285
--- /dev/null
+++ b/AMDiS/src/nonlin/ProblemNonLin.h
@@ -0,0 +1,158 @@
+// ============================================================================
+// == ==
+// == 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;
+ }
+
+ /** \brief
+ * Sets \ref updater.
+ */
+ inline void setUpdater(NonLinUpdater *updater)
+ {
+ updater_ = updater;
+ }
+
+
+ protected:
+ /** \brief
+ * Initial guess for the solution.
+ */
+ Vector<AbstractFunction<double, WorldVector<double> >*> u0_;
+
+ /** \brief
+ * Non linear solver. Used in \ref solve().
+ */
+ NonLinSolver *nonLinSolver_;
+
+ /** \brief
+ * Updator
+ */
+ NonLinUpdater *updater_;
+ };
+
+}
+
+#endif
+
--
GitLab