diff --git a/AMDiS/CMakeLists.txt b/AMDiS/CMakeLists.txt
index 16e74b64068286014adc49f78df89b7f93aabd8a..66bf10c4abab75d00d93a6c419ac64440a952276 100644
--- a/AMDiS/CMakeLists.txt
+++ b/AMDiS/CMakeLists.txt
@@ -145,7 +145,6 @@ 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
@@ -323,6 +322,11 @@ INSTALL(FILES ${HEADERS}
DESTINATION include/amdis/parallel/)
list(APPEND deb_add_dirs "include/amdis/parallel")
+FILE(GLOB HEADERS "${SOURCE_DIR}/nonlin/*.h")
+INSTALL(FILES ${HEADERS}
+ DESTINATION include/amdis/nonlin/)
+list(APPEND deb_add_dirs "include/amdis/nonlin")
+
FILE(GLOB HEADERS "${SOURCE_DIR}/time/*.h")
INSTALL(FILES ${HEADERS}
DESTINATION include/amdis/time/)
diff --git a/AMDiS/src/AMDiS.h b/AMDiS/src/AMDiS.h
index b950d31a3f8f09f77f3408152b58c4dd5559e7d7..fae091ce4769c8838dcdac13838cffa4a311cc97 100644
--- a/AMDiS/src/AMDiS.h
+++ b/AMDiS/src/AMDiS.h
@@ -122,9 +122,13 @@
#include "io/ValueWriter.h"
#include "io/VtkWriter.h"
+#include "nonlin/ProblemNonLin.h"
+#include "nonlin/NonLinSolver.h"
+
#include "time/RosenbrockAdaptInstationary.h"
#include "time/RosenbrockStationary.h"
+
#if HAVE_PARALLEL_DOMAIN_AMDIS
#include "parallel/InteriorBoundary.h"
#include "parallel/MpiHelper.h"
diff --git a/AMDiS/src/CreatorMap.cc b/AMDiS/src/CreatorMap.cc
index 70091af41569a29112569ae36ca91eb1b85268ce..d21ba5c13814f819ac4331c512b2f884c602aca9 100644
--- a/AMDiS/src/CreatorMap.cc
+++ b/AMDiS/src/CreatorMap.cc
@@ -26,6 +26,7 @@
#include "DOFMatrix.h"
#include "UmfPackSolver.h"
#include "time/RosenbrockMethod.h"
+#include "nonlin/NonLinSolver.h"
namespace AMDiS {
@@ -142,4 +143,11 @@ namespace AMDiS {
addCreator("ros3p", new Ros3p::Creator);
addCreator("rodasp", new Rodasp::Creator);
}
+
+
+ template<>
+ void CreatorMap<NonLinSolver>::addDefaultCreators()
+ {
+ addCreator("newton", new Newton::Creator);
+ }
}
diff --git a/AMDiS/src/ProblemStat.cc b/AMDiS/src/ProblemStat.cc
index 97e56a52a574d86d6bc91d3998e1966a10ad886e..8e29c3d161c65b4083c3648a2a56e9cf44aa46fc 100644
--- a/AMDiS/src/ProblemStat.cc
+++ b/AMDiS/src/ProblemStat.cc
@@ -634,8 +634,6 @@ namespace AMDiS {
{
FUNCNAME("ProblemStat::buildAfterCoarsen()");
- VtkWriter::writeFile(rhs->getDOFVector(0), "test.vtu");
-
// buildAfterCoarsen_sebastianMode(adaptInfo, flag);
// return;
@@ -780,7 +778,6 @@ namespace AMDiS {
if (asmMatrix) {
solverMatrix.setMatrix(*systemMatrix);
- createPrecon();
INFO(info, 8)("fillin of assembled matrix: %d\n", nnz);
}
@@ -971,8 +968,6 @@ namespace AMDiS {
solverMatrix.setMatrix(*systemMatrix);
- createPrecon();
-
INFO(info, 8)("fillin of assembled matrix: %d\n", nnz);
INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n",
TIME_USED(first, clock()));
@@ -1232,9 +1227,6 @@ namespace AMDiS {
solverMatrix.setMatrix(*systemMatrix);
- createPrecon();
-
-
INFO(info, 8)("fillin of assembled matrix: %d\n", nnz);
INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n",
@@ -1242,24 +1234,6 @@ namespace AMDiS {
}
- void ProblemStatSeq::createPrecon()
- {
- /*string preconType("no");
- Parameters::get(name + "->solver->left precon", preconType);
-
- CreatorInterface<ITL_BasePreconditioner> *preconCreator =
- CreatorMap<ITL_BasePreconditioner>::getCreator(preconType);
-
- solver->setLeftPrecon(preconCreator->create(solverMatrix.getMatrix()));
-
- preconType= "no";
- Parameters::get(name + "->solver->right precon", preconType);
-
- preconCreator = CreatorMap<ITL_BasePreconditioner>::getCreator(preconType);
- solver->setRightPrecon(preconCreator->create(solverMatrix.getMatrix()));*/
- }
-
-
void ProblemStatSeq::writeFiles(AdaptInfo *adaptInfo, bool force)
{
FUNCNAME("ProblemStat::writeFiles()");
@@ -1570,153 +1544,43 @@ namespace AMDiS {
Mesh *mesh = feSpace->getMesh();
const BasisFunction *basisFcts = feSpace->getBasisFcts();
-#ifdef _OPENMP
- TraverseParallelStack stack(0, 1);
-#else
TraverseStack stack;
-#endif
-
- // == Initialize matrix and vector. If we have to assemble in parallel, ===
- // == temporary thread owned matrix and vector must be created. ===
-
-#ifdef _OPENMP
-#pragma omp parallel
-#endif
- {
- BoundaryType *bound =
- useGetBound ? new BoundaryType[basisFcts->getNumber()] : NULL;
- // Create for every thread its private matrix and vector, on that
- // the thread will assemble its part of the mesh.
- DOFMatrix *tmpMatrix = NULL;
- DOFVector<double> *tmpVector = NULL;
-
-#ifdef _OPENMP
- if (matrix) {
- tmpMatrix = new DOFMatrix(matrix->getRowFeSpace(), matrix->getColFeSpace(), "tmp");
-
- // Copy the global matrix to the private matrix, because we need the
- // operators defined on the global matrix in the private one. Only the
- // values have to be set to zero.
- *tmpMatrix = *matrix;
- tmpMatrix->clear();
+ BoundaryType *bound =
+ useGetBound ? new BoundaryType[basisFcts->getNumber()] : NULL;
- tmpMatrix->getBaseMatrix().change_dim(matrix->getRowFeSpace()->getAdmin()->getUsedSize(),
- matrix->getColFeSpace()->getAdmin()->getUsedSize());
- tmpMatrix->startInsertion(10);
- }
+ matrix->startInsertion(matrix->getNnz());
+ vector->set(0.0);
- if (vector) {
- tmpVector = new DOFVector<double>(vector->getFeSpace(), "tmp");
+ // == Traverse mesh and assemble. ==
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag);
+ while (elInfo) {
+ if (useGetBound)
+ basisFcts->getBound(elInfo, bound);
- // Copy the global vector to the private vector, because we need the
- // operatirs defined on the global vector in the private one. But set
- // the values to zero of the private vector after copying.
- *tmpVector = *vector;
- tmpVector->set(0.0);
- }
-#else
if (matrix) {
- tmpMatrix = matrix;
- tmpMatrix->startInsertion(matrix->getNnz());
- }
- if (vector) {
- tmpVector = vector;
- tmpVector->set(0.0);
- }
-#endif
-
- // == Traverse mesh (either sequentially or in parallel) and assemble. ==
-
-
- // Because we are using the parallel traverse stack, each thread will
- // traverse only a part of the mesh.
- ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag);
-
- // After creating privat copies of the DOFMatrix and the DOFVector, all threads
- // have to wait at this barrier. Especially for small problems this is required,
- // because otherwise one thread may be finished with assembling, before another
- // has made his private copy.
-#ifdef _OPENMP
-#pragma omp barrier
-#endif
- while (elInfo) {
- if (useGetBound)
- basisFcts->getBound(elInfo, bound);
-
- if (matrix) {
- tmpMatrix->assemble(1.0, elInfo, bound);
-
- // Take the matrix boundary manager from the public matrix,
- // but assemble the boundary conditions on the thread private matrix.
- if (matrix->getBoundaryManager())
- matrix->getBoundaryManager()->fillBoundaryConditions(elInfo, tmpMatrix);
- }
+ matrix->assemble(1.0, elInfo, bound);
- if (vector)
- tmpVector->assemble(1.0, elInfo, bound, NULL);
-
- elInfo = stack.traverseNext(elInfo);
- }
-
- // == Finally, if we have assembled in parallel, we have to add the thread ==
- // == private matrix and vector to the global one. ==
-
-#ifdef _OPENMP
- if (tmpMatrix)
- tmpMatrix->finishInsertion();
-
- // After mesh traverse, all thread have to added their private matrices and
- // vectors to the global public matrix and public vector. Therefore, this is
- // a critical section, which is allowed to be executed by on thread only at
- // the same time.
-
- if (matrix) {
-#pragma omp critical
- {
- matrix->getBaseMatrix() += tmpMatrix->getBaseMatrix();
- }
- }
-
-#pragma omp barrier
-
-#pragma omp master
- {
- if (matrix)
- matrix->startInsertion();
- }
-
-#pragma omp barrier
-
- if (matrix) {
- // Remove rows corresponding to DOFs on a Dirichlet boundary.
-#pragma omp critical
- matrix->removeRowsWithDBC(tmpMatrix->getApplyDBCs());
-
- delete tmpMatrix;
- }
-
- if (vector) {
-#pragma omp critical
- *vector += *tmpVector;
-
- delete tmpVector;
+ // Take the matrix boundary manager from the public matrix,
+ // but assemble the boundary conditions on the thread private matrix.
+ if (matrix->getBoundaryManager())
+ matrix->getBoundaryManager()->fillBoundaryConditions(elInfo, matrix);
}
+
+ if (vector)
+ vector->assemble(1.0, elInfo, bound, NULL);
-#else
- if (matrix)
- matrix->removeRowsWithDBC(matrix->getApplyDBCs());
-#endif
+ elInfo = stack.traverseNext(elInfo);
+ }
- if (useGetBound)
- delete [] bound;
+ // == Finally, if we have assembled in parallel, we have to add the thread ==
+ // == private matrix and vector to the global one. ==
- } // pragma omp parallel
+ if (matrix)
+ matrix->removeRowsWithDBC(matrix->getApplyDBCs());
-// INFO(info, 8)("Assemble matrix with %d mat ops and %d vec ops needed %.5f seconds\n",
-// matrix ? matrix->getOperators().size() : 0,
-// vector ? vector->getOperators().size() : 0,
-// TIME_USED(first, clock()));
+ if (useGetBound)
+ delete [] bound;
}
@@ -1727,6 +1591,7 @@ namespace AMDiS {
{
FUNCNAME("ProblemStat::assembleBoundaryConditions()");
+
// === Initialization of vectors ===
if (rhs->getBoundaryManager())
@@ -1734,55 +1599,6 @@ namespace AMDiS {
if (solution->getBoundaryManager())
solution->getBoundaryManager()->initVector(solution);
-#ifdef _OPENMP
- // === Parallel boundary assemblage ===
-
- TraverseParallelStack stack;
-
-#pragma omp parallel
- {
- // Each thread assembles on its own dof-vectors.
- DOFVector<double> *tmpRhsVec = new DOFVector<double>(rhs->getFeSpace(), "tmpRhs");
- DOFVector<double> *tmpSolVec = new DOFVector<double>(solution->getFeSpace(), "tmpSol");
- tmpRhsVec->set(0.0);
- tmpSolVec->set(0.0);
-
- // (Parallel) traverse of mesh.
- ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag);
- while (elInfo) {
- if (rhs->getBoundaryManager())
- rhs->getBoundaryManager()-> fillBoundaryConditions(elInfo, tmpRhsVec);
-
- if (solution->getBoundaryManager())
- solution->getBoundaryManager()->fillBoundaryConditions(elInfo, tmpSolVec);
-
- elInfo = stack.traverseNext(elInfo);
- }
-
- // After (parallel) mesh traverse, the result is applied to the final
- // vectors. This section is not allowed to be executed by more than one
- // thread at the same time.
-#pragma omp critical
- {
- DOFVector<double>::Iterator rhsIt(rhs, USED_DOFS);
- DOFVector<double>::Iterator solIt(solution, USED_DOFS);
- DOFVector<double>::Iterator tmpRhsIt(tmpRhsVec, USED_DOFS);
- DOFVector<double>::Iterator tmpSolIt(tmpSolVec, USED_DOFS);
- for (rhsIt.reset(), solIt.reset(), tmpRhsIt.reset(), tmpSolIt.reset();
- !rhsIt.end();
- ++rhsIt, ++solIt, ++tmpRhsIt, ++tmpSolIt) {
- *rhsIt += *tmpRhsIt;
- *solIt += *tmpSolIt;
- }
- } // pragma omp critical
-
- delete tmpRhsVec;
- delete tmpSolVec;
- } // pragma omp parallel
-
-#else
- // === Sequentiell boundary assemblage ===
-
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag);
while (elInfo) {
@@ -1794,7 +1610,6 @@ namespace AMDiS {
elInfo = stack.traverseNext(elInfo);
}
-#endif
// === Finalize vectors ===
diff --git a/AMDiS/src/ProblemStat.h b/AMDiS/src/ProblemStat.h
index 75ee265f8298c93bf17759ff36d0751301b936b9..7a541914015f8a4595d50a503317a15bc141f7a4 100644
--- a/AMDiS/src/ProblemStat.h
+++ b/AMDiS/src/ProblemStat.h
@@ -48,6 +48,7 @@ namespace AMDiS {
Flag operatorType;
};
+
/** \brief
* This class defines the stationary problem definition in sequential
* computations. For parallel computations, see
@@ -195,9 +196,6 @@ namespace AMDiS {
void dualAssemble(AdaptInfo *adaptInfo, Flag flag,
bool asmMatrix = true, bool asmVector = true);
- void createPrecon();
-
-
/** \brief
* Determines the execution order of the single adaption steps. If adapt is
* true, mesh adaption will be performed. This allows to avoid mesh adaption,
diff --git a/AMDiS/src/nonlin/Newton.h b/AMDiS/src/nonlin/Newton.h
index 2be9aaeb5d3a118e6aba5b2fed0fddeb4e11fe00..7cf430aba217222efec1cfed081b989fb93d7214 100644
--- a/AMDiS/src/nonlin/Newton.h
+++ b/AMDiS/src/nonlin/Newton.h
@@ -26,6 +26,7 @@
#include "CreatorInterface.h"
#include "NonLinSolver.h"
#include "OEMSolver.h"
+#include "io/VtkWriter.h"
namespace AMDiS {
@@ -45,114 +46,93 @@ namespace AMDiS {
public:
virtual ~Creator() {}
- /** \brief
- * Returns a new Newton object.
- */
+ /// Returns a new Newton object.
NonLinSolver* create()
{
- return new Newton(this->name, this->linearSolver, this->nonLinUpdater);
+ return new Newton(this->name, this->linearSolver);
}
};
- /** \brief
- * Calls constructor of base class NonLinSolver
- */
- Newton(const std::string& name_,
- OEMSolver *linSolver_,
- NonLinUpdater *updater)
- : NonLinSolver(name_, linSolver_, updater),
+ /// Calls constructor of base class NonLinSolver
+ Newton(const std::string& name, OEMSolver *linSolver)
+ : NonLinSolver(name, linSolver),
b(NULL)
{}
private:
- /** \brief
- * realisation of NonLinSolver::init
- */
- void init()
- {}
+ /// Realisation of NonLinSolver::init
+ void init() {}
- /** \brief
- * realisation of NonLinSolver::nlsolve
- */
+ /// realisation of NonLinSolver::nlsolve
int nlsolve(SolverMatrix<Matrix<DOFMatrix*> >& mat,
- SystemVector& x, SystemVector& rhs)
+ SystemVector& x, SystemVector& rhs,
+ AdaptInfo *adaptInfo,
+ ProblemStat *prob)
{
- //DOFVector<T> *d = problem->getRHS();
- //DOFVector<T> *x = problem->getSolution();;
-
- b = new SystemVector(x);
- *b = rhs;
+ FUNCNAME("Newton::nlsolve()");
- // // copy operators from fh to b
- // std::vector<Operator*>::iterator op;
- // std::vector<double*>::iterator fac;
- // for(op = d->getOperatorsBegin(),
- // fac = d->getOperatorFactorBegin();
- // op != d->getOperatorsEnd();
- // ++op, ++fac)
- // {
- // b->addOperator(*op, *fac);
- // }
+ if (b == NULL)
+ b = new SystemVector(x);
- double err = 0.0, err_old = -1.0;
- int iter, n;
+ double err = 0.0, errOld = -1.0;
+ int iter, n;
- INFO(this->info,2)("iter. | this->residual | red. | n |\n");
+ MSG("iter. | this->residual | red. | n |\n");
for (iter = 1; iter <= this->maxIter; iter++) {
- /*--- Assemble DF(x) and F(x) ----------------------------------------------*/
- this->nonLinUpdater->update(/*x,*/true, b);
- /*--- Initial guess is zero ------------------------------------------------*/
- rhs.set(0.0);
- /*--- solve linear system --------------------------------------------------*/
+ // Assemble DF(x) and F(x)
+ prob->buildAfterCoarsen(adaptInfo, 0, true, true);
+
+ // Initial guess is zero
+ b->set(0.0);
+
+ // Solve linear system
n = solveLinearSystem(mat, *b, rhs);
- /*--- x = x - d ------------------------------------------------------------*/
- x -= rhs;
+
+ // x = x + d
+ x += *b;
if (this->usedNorm == NO_NORM || this->usedNorm == L2_NORM)
- err = L2Norm(&rhs); // sollte hier nicht b genommen werden (s. NewtonS) ?
+ err = L2Norm(b);
else
- err = H1Norm(&rhs); // sollte hier nicht b genommen werden (s. NewtonS) ?
+ err = H1Norm(b);
- if (iter == 1) this->initial_residual = err;
-
- if (err_old <= 0) {
- INFO(this->info,2)("%5d | %12.5e | -------- | %4d |\n", iter, err, n);
- } else {
- INFO(this->info,2)("%5d | %12.5e | %8.2e | %4d |\n",
- iter, err, err/err_old, n);
- }
+ if (iter == 1)
+ this->initialResidual = err;
- if ((this->residual = err) < this->tolerance) {
- INFO(this->info,4)("finished successfully\n");
- return iter;
- }
- err_old = err;
+ if (errOld <= 0)
+ MSG("%5d | %12.5e | -------- | %4d |\n", iter, err, n);
+ else
+ MSG("%5d | %12.5e | %8.2e | %4d |\n", iter, err, err/errOld, n);
+
+ residual = err;
+ if (err < this->tolerance) {
+ MSG("Finished successfully!\n");
+ return iter;
+ }
+ errOld = err;
}
- if (this->info < 2)
- INFO(this->info,1)("iter. %d, residual: %12.5e\n", iter, err);
- INFO(this->info,1)("tolerance %e not reached\n", this->tolerance);
+ MSG("iter. %d, residual: %12.5e\n", iter, err);
+ MSG("tolerance %e not reached\n", this->tolerance);
this->residual = err;
return iter;
}
- /** \brief
- * realisation of NonLinSolver::exit
- */
+ /// Realisation of NonLinSolver::exit
void exit()
{
- if (b != NULL)
+ if (b != NULL) {
delete b;
+ b = NULL;
+ }
}
private:
- /** \brief
- * internal used data
- */
+ /// Internal used data
SystemVector *b;
};
diff --git a/AMDiS/src/nonlin/NonLinSolver.h b/AMDiS/src/nonlin/NonLinSolver.h
index 16ee8a2875684f4c88af0b757f494029516c3732..559f9fbd9c15ee2e4dde700fd7199850e3ab0ff6 100644
--- a/AMDiS/src/nonlin/NonLinSolver.h
+++ b/AMDiS/src/nonlin/NonLinSolver.h
@@ -30,6 +30,7 @@
#include "MatrixVector.h"
#include "DOFMatrix.h"
#include "OEMSolver.h"
+#include "ProblemStat.h"
namespace AMDiS {
@@ -37,27 +38,23 @@ 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.
+ * Pure virtual base class for specific non linear solvers.
*/
class NonLinSolver
{
public:
/** \brief
- * constructor.
- * \param name name of this solver
+ * Constructor
+ *
+ * \param name Name of this solver
*/
- NonLinSolver(const std::string &name_,
- OEMSolver *solver,
- NonLinUpdater *updater)
- : name(name_),
+ NonLinSolver(const std::string &name, OEMSolver *solver)
+ : name(name),
linSolver(solver),
- nonLinUpdater(updater),
- tolerance(1.e-4),
+ tolerance(1.e-8),
maxIter(50),
info(8),
- initial_residual(0.0),
+ initialResidual(0.0),
residual(0.0),
usedNorm(NO_NORM)
{
@@ -67,20 +64,22 @@ namespace AMDiS {
Parameters::get(name + "->norm", usedNorm);
}
- /** \brief
- * destructor
- */
+ /// Destructor
virtual ~NonLinSolver() {}
/** \brief
* solves the non linear system. Uses sub class methods
*/
inline int solve(SolverMatrix<Matrix<DOFMatrix*> >& mat,
- SystemVector &x, SystemVector &rhs)
+ SystemVector &solution, SystemVector &rhs,
+ AdaptInfo *adaptInfo,
+ ProblemStat *prob)
{
init();
- int result = nlsolve(mat, x, rhs);
+ solution.set(0.0);
+ int result = nlsolve(mat, solution, rhs, adaptInfo, prob);
exit();
+
return result;
}
@@ -99,68 +98,63 @@ namespace AMDiS {
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.
- */
+ /// 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.
- */
+ /// Solves the non linear system. Must be overriden in sub classes.
virtual int nlsolve(SolverMatrix<Matrix<DOFMatrix*> >& matVec,
- SystemVector& x, SystemVector& rhs) = 0;
+ SystemVector& x, SystemVector& rhs,
+ AdaptInfo *adaptInfo,
+ ProblemStat *prob) = 0;
- /** \brief
- * Frees needed memory. Must be overriden in sub classes.
- */
+ /// Frees needed memory. Must be overriden in sub classes.
virtual void exit() = 0;
virtual int solveLinearSystem(SolverMatrix<Matrix<DOFMatrix*> >& mat,
- SystemVector &fh, SystemVector &x)
+ SystemVector &x, SystemVector &b)
{
FUNCNAME("NonLinSolver::solveLinearSystem()");
TEST_EXIT(linSolver)("no solver\n");
- return linSolver->solveSystem(mat, x, fh);
+ return linSolver->solveSystem(mat, x, b);
}
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 */
+ /// Name of the solver.
+ std::string name;
+
+ /// Linear solver object.
+ OEMSolver *linSolver;
+
+ /// Solver tolerance.
+ double tolerance;
+
+ /// Maximal number of iterations.
+ int maxIter;
+
+ /// Info level.
+ int info;
+
+ /// Initial residual.
+ double initialResidual;
+
+ /// Current residual.
+ double residual;
+
+ /// Used norm for convergent test.
+ Norm usedNorm;
};
- // ============================================================================
- // ===== class NonLinSolverCreator ============================================
- // ============================================================================
- /** \brief
- * Interface for creators of concrete NonLinSolvers.
- */
+ /// Interface for creators of concrete NonLinSolvers.
class NonLinSolverCreator : public CreatorInterface<NonLinSolver>
{
public:
virtual ~NonLinSolverCreator() {}
- void setName(std::string name_)
+ void setName(std::string n)
{
- name = name_;
+ name = n;
}
void setLinearSolver(OEMSolver *solver)
@@ -168,15 +162,10 @@ namespace AMDiS {
linearSolver = solver;
}
- void setNonLinUpdater(NonLinUpdater *updater)
- {
- nonLinUpdater = updater;
- }
-
protected:
std::string name;
+
OEMSolver *linearSolver;
- NonLinUpdater *nonLinUpdater;
};
}
diff --git a/AMDiS/src/nonlin/NonLinUpdater.cc b/AMDiS/src/nonlin/NonLinUpdater.cc
deleted file mode 100644
index 451c70a4ef2c98927e55df36bfc4193307ee11b7..0000000000000000000000000000000000000000
--- a/AMDiS/src/nonlin/NonLinUpdater.cc
+++ /dev/null
@@ -1,104 +0,0 @@
-//
-// 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
deleted file mode 100644
index 3092845eb255845cbdc2ba410418b8efc081fc42..0000000000000000000000000000000000000000
--- a/AMDiS/src/nonlin/NonLinUpdater.h
+++ /dev/null
@@ -1,52 +0,0 @@
-// ============================================================================
-// == ==
-// == 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
index 1c66a26962d2944d825a6a1ad4a968d1621b06e1..002de8b99c7594dfc823f669629f7732d4103032 100644
--- a/AMDiS/src/nonlin/ProblemNonLin.cc
+++ b/AMDiS/src/nonlin/ProblemNonLin.cc
@@ -13,7 +13,6 @@
#include "ProblemNonLin.h"
#include "NonLinSolver.h"
#include "CreatorMap.h"
-#include "NonLinUpdater.h"
#include "Traverse.h"
#include "AdaptInfo.h"
@@ -23,43 +22,25 @@ namespace AMDiS {
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();
- }
- }
+ FUNCNAME("ProblemNonLin::initialize()");
- if (updater_ == NULL)
- WARNING("no updater created\n");
+ ProblemStat::initialize(initFlag, adoptProblem, adoptFlag);
// === create nonlinear solver ===
- if (nonLinSolver_) {
- WARNING("nonlinear solver already created\n");
+ 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();
+ TEST_EXIT(nonLinSolver == NULL)("Nonlinear solver already created!\n");
+ nonLinSolver = dynamic_cast<ProblemNonLin*>(adoptProblem)->getNonLinSolver();
}
}
- if (nonLinSolver_ == NULL)
- WARNING("no nonlinear solver created\n");
+ if (nonLinSolver == NULL)
+ WARNING("No nonlinear solver created!\n");
}
@@ -74,51 +55,53 @@ namespace AMDiS {
nonLinSolverCreator->setLinearSolver(solver);
nonLinSolverCreator->setName(name + "->nonlin solver");
- nonLinSolverCreator->setNonLinUpdater(updater_);
- nonLinSolver_ = nonLinSolverCreator->create();
+ nonLinSolver = nonLinSolverCreator->create();
}
- void ProblemNonLin::solve(AdaptInfo *adaptInfo)
+ void ProblemNonLin::solve(AdaptInfo *adaptInfo, bool)
{
- TEST_EXIT(nonLinSolver_)("no non-linear solver!\n");
- nonLinSolver_->solve(solverMatrix, *solution, *rhs);
+ TEST_EXIT(nonLinSolver)("no non-linear solver!\n");
+
+ nonLinSolver->solve(solverMatrix, *solution, *rhs, adaptInfo, this);
}
- void ProblemNonLin::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag, bool, bool)
+ void ProblemNonLin::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag f, bool b0, bool b1)
{
FUNCNAME("ProblemNonLin::buildAfterCoarsen()");
+
+ ProblemStat::buildAfterCoarsen(adaptInfo, f, b0, b1);
+
+#if 0
int nMeshes = static_cast<int>(meshes.size());
for (int i = 0; i < nMeshes; i++)
meshes[i]->dofCompress();
- for (int i = 0; i < nComponents; i++) {
+ 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;
+ componentSpaces[i]->getName().c_str());
// 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);
+ TraverseStack stack;
+ ElInfo *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()) {
+ if (solution->getDOFVector(i)->getBoundaryManager())
solution->getDOFVector(i)->
getBoundaryManager()->fillBoundaryConditions(elInfo, solution->getDOFVector(i));
- }
+
elInfo = stack.traverseNext(elInfo);
}
}
+#endif
}
}
diff --git a/AMDiS/src/nonlin/ProblemNonLin.h b/AMDiS/src/nonlin/ProblemNonLin.h
index ded55285fdc05f9c8481e78c14103876c12c48d5..91a59d66bf9d3ca71bea0f34fd8beb6183a887cd 100644
--- a/AMDiS/src/nonlin/ProblemNonLin.h
+++ b/AMDiS/src/nonlin/ProblemNonLin.h
@@ -24,7 +24,6 @@
#define AMDIS_PROBLEMNONLIN_H_
#include "ProblemStat.h"
-#include "NonLinUpdater.h"
#include "NonLinSolver.h"
#include "DOFVector.h"
#include "SystemVector.h"
@@ -44,112 +43,64 @@ namespace AMDiS {
/// Constructs a ProblemNonLin with given name.
ProblemNonLin(const std::string& name_)
: ProblemStat(name_.c_str()),
- nonLinSolver_(NULL),
- updater_(NULL)
+ nonLinSolver(NULL)
{
- u0_.resize(nComponents);
- u0_.set(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,
+ /// 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;
+
+ TEST_EXIT(index < nComponents)("Invalid index!\n");
+ u0[index] = u0Fct;
solution->getDOFVector(index)->interpol(u0Fct);
}
- /** \brief
- * Destructor.
- */
- virtual ~ProblemNonLin() {}
+ /// Destructor
+ virtual ~ProblemNonLin()
+ {}
- /** \brief
- * Initialisation of the problem.
- */
+ /// Initialization 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().
- */
+ /// Used in \ref initialize().
virtual void createNonLinSolver();
/** \brief
- * Overrides ProblemVec::solve(). Uses the non linear solver
- * \ref nonLinSolver_.
+ * Overrides ProblemStat::solve(). Uses the non linear solver
+ * \ref nonLinSolver.
*/
- virtual void solve(AdaptInfo *adaptInfo);
+ virtual void solve(AdaptInfo *adaptInfo, bool fixedMatrix = false);
- /** \brief
- * Overrides Problem::buildAfterCoarsen(). Sets dirichlet boundaries
- * in \ref A_ and \ref rhs_.
- */
+ /// Overrides ProblemStat::buildAfterCoarsen().
virtual void buildAfterCoarsen(AdaptInfo *adaptInfo, Flag,
bool assembleMatrix = true,
bool assembleVector = true);
- /** \brief
- * Returns /ref nonLinSolver_.
- */
+ /// Returns \ref nonLinSolver.
inline NonLinSolver *getNonLinSolver()
{
- return nonLinSolver_;
- }
-
- /** \brief
- * Returns \ref updater.
- */
- inline NonLinUpdater* getUpdater()
- {
- return updater_;
+ return nonLinSolver;
}
- /** \brief
- * Sets \ref nonLinSolver_.
- */
- inline void setNonLinSolver(NonLinSolver *nonLinSolver)
+ /// Sets \ref nonLinSolver.
+ inline void setNonLinSolver(NonLinSolver *s)
{
- nonLinSolver_ = nonLinSolver;
+ nonLinSolver = s;
}
- /** \brief
- * Sets \ref updater.
- */
- inline void setUpdater(NonLinUpdater *updater)
- {
- updater_ = updater;
- }
-
-
protected:
- /** \brief
- * Initial guess for the solution.
- */
- Vector<AbstractFunction<double, WorldVector<double> >*> u0_;
+ /// 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_;
+ /// Non linear solver object. Used in \ref solve().
+ NonLinSolver *nonLinSolver;
};
}
diff --git a/demo/init/nonlin.dat.2d b/demo/init/nonlin.dat.2d
new file mode 100644
index 0000000000000000000000000000000000000000..66bf1a8241ab5f9fe2e83ba0b51b285690c09666
--- /dev/null
+++ b/demo/init/nonlin.dat.2d
@@ -0,0 +1,33 @@
+dimension of world: 2
+
+nonlinMesh->macro file name: ./macro/macro.stand.2d
+nonlinMesh->global refinements: 10
+
+nonlin->mesh: nonlinMesh
+nonlin->dim: 2
+nonlin->polynomial degree[0]: 1
+nonlin->components: 1
+
+nonlin->nonlin solver: newton
+nonlin->solver: umfpack
+nonlin->solver->max iteration: 1000
+nonlin->solver->tolerance: 1.e-8
+nonlin->solver->info: 2
+nonlin->solver->left precon: no
+nonlin->solver->right precon: no
+
+nonlin->estimator[0]: residual
+nonlin->estimator[0]->error norm: 1 % 1: H1_NORM, 2: L2_NORM
+nonlin->estimator[0]->C0: 0.1 % constant of element residual
+nonlin->estimator[0]->C1: 0.1 % constant of jump residual
+
+nonlin->marker[0]->strategy: 0 % 0: no adaption 1: GR 2: MS 3: ES 4:GERS
+nonlin->marker[0]->MSGamma: 0.5
+
+nonlin->adapt[0]->tolerance: 1e-4
+nonlin->adapt[0]->refine bisections: 2
+
+nonlin->adapt->max iteration: 0
+
+nonlin->output->filename: output/nonlin
+nonlin->output->ParaView format: 1
diff --git a/demo/src/nonlin.cc b/demo/src/nonlin.cc
new file mode 100644
index 0000000000000000000000000000000000000000..0c0cdc4d6c78b515b6e128143350c493fdc07626
--- /dev/null
+++ b/demo/src/nonlin.cc
@@ -0,0 +1,131 @@
+#include "AMDiS.h"
+
+using namespace AMDiS;
+using namespace std;
+
+// ===========================================================================
+// ===== function definitions ================================================
+// ===========================================================================
+
+/// Dirichlet boundary function
+class G : public AbstractFunction<double, WorldVector<double> >
+{
+public:
+
+ G()
+ : AbstractFunction<double, WorldVector<double> >()
+ {}
+
+ /// Implementation of AbstractFunction::operator().
+ double operator()(const WorldVector<double>& x) const
+ {
+ return 0.0;
+ }
+};
+
+
+/// RHS function
+class F : public AbstractFunction<double, WorldVector<double> >
+{
+public:
+
+ F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}
+
+ /// Implementation of AbstractFunction::operator().
+ double operator()(const WorldVector<double>& x) const
+ {
+ int dow = Global::getGeo(WORLD);
+ double r2 = (x * x);
+ double ux = exp(-10.0 * r2);
+ return -(400.0 * r2 - 20.0 * dow) * ux + pow(ux, 4.0);
+ }
+};
+
+
+class NonlinFctLeft : public AbstractFunction<double, double>
+{
+public:
+ NonlinFctLeft()
+ : AbstractFunction<double, double>()
+ {}
+
+ double operator()(const double& u) const
+ {
+ return 4.0 * u * u * u;
+ }
+};
+
+
+class NonlinFctRight : public AbstractFunction<double, double>
+{
+public:
+ NonlinFctRight()
+ : AbstractFunction<double, double>()
+ {}
+
+ double operator()(const double& u) const
+ {
+ return -(u * u * u * u);
+ }
+};
+
+// ===========================================================================
+// ===== main program ========================================================
+// ===========================================================================
+
+int main(int argc, char* argv[])
+{
+ FUNCNAME("main");
+
+ // ===== check for init file =====
+ TEST_EXIT(argc >= 2)("usage: nonlin initfile\n");
+
+
+ // ===== init parameters =====
+ Parameters::init(true, argv[1]);
+
+
+ // ===== create and init the scalar problem =====
+ ProblemNonLin nonlin("nonlin");
+ nonlin.initialize(INIT_ALL);
+
+ // === create adapt info ===
+ AdaptInfo adaptInfo("nonlin->adapt");
+
+
+ // === create adapt ===
+ AdaptStationary adapt("nonlin->adapt", nonlin, adaptInfo);
+
+
+ // === create matrix operators ===
+ Operator mat01(nonlin.getFeSpace());
+ mat01.addSecondOrderTerm(new Simple_SOT());
+ nonlin.addMatrixOperator(mat01, 0, 0);
+
+ Operator mat02(nonlin.getFeSpace());
+ mat02.addZeroOrderTerm(new VecAtQP_ZOT(nonlin.getSolution(0), new NonlinFctLeft()));
+ nonlin.addMatrixOperator(mat02, 0, 0);
+
+ // === create rhs operators ===
+ Operator vec01(nonlin.getFeSpace());
+ vec01.setUhOld(nonlin.getSolution(0));
+ vec01.addSecondOrderTerm(new Simple_SOT(-1.0));
+ nonlin.addVectorOperator(vec01, 0);
+
+ Operator vec02(nonlin.getFeSpace());
+ vec02.addZeroOrderTerm(new VecAtQP_ZOT(nonlin.getSolution(0), new NonlinFctRight()));
+ nonlin.addVectorOperator(vec02, 0);
+
+ Operator vec03(nonlin.getFeSpace());
+ vec03.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(nonlin.getFeSpace()->getBasisFcts()->getDegree())));
+ nonlin.addVectorOperator(vec03, 0);
+
+ // ===== add boundary conditions =====
+ nonlin.addDirichletBC(1, 0, 0, new G);
+
+ // ===== start adaption loop =====
+ adapt.adapt();
+
+
+ nonlin.writeFiles(adaptInfo, true);
+}