Commit e59d6928 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

First working nonlinear solver (but very restricted).

parent deb831ee
...@@ -145,7 +145,6 @@ SET(AMDIS_SRC ${SOURCE_DIR}/AdaptBase.cc ...@@ -145,7 +145,6 @@ SET(AMDIS_SRC ${SOURCE_DIR}/AdaptBase.cc
${SOURCE_DIR}/io/ValueReader.cc ${SOURCE_DIR}/io/ValueReader.cc
${SOURCE_DIR}/io/ValueWriter.cc ${SOURCE_DIR}/io/ValueWriter.cc
${SOURCE_DIR}/io/VtkWriter.cc ${SOURCE_DIR}/io/VtkWriter.cc
${SOURCE_DIR}/nonlin/NonLinUpdater.cc
${SOURCE_DIR}/nonlin/ProblemNonLin.cc ${SOURCE_DIR}/nonlin/ProblemNonLin.cc
${SOURCE_DIR}/parallel/InteriorBoundary.cc ${SOURCE_DIR}/parallel/InteriorBoundary.cc
${SOURCE_DIR}/time/RosenbrockAdaptInstationary.cc ${SOURCE_DIR}/time/RosenbrockAdaptInstationary.cc
...@@ -323,6 +322,11 @@ INSTALL(FILES ${HEADERS} ...@@ -323,6 +322,11 @@ INSTALL(FILES ${HEADERS}
DESTINATION include/amdis/parallel/) DESTINATION include/amdis/parallel/)
list(APPEND deb_add_dirs "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") FILE(GLOB HEADERS "${SOURCE_DIR}/time/*.h")
INSTALL(FILES ${HEADERS} INSTALL(FILES ${HEADERS}
DESTINATION include/amdis/time/) DESTINATION include/amdis/time/)
......
...@@ -122,9 +122,13 @@ ...@@ -122,9 +122,13 @@
#include "io/ValueWriter.h" #include "io/ValueWriter.h"
#include "io/VtkWriter.h" #include "io/VtkWriter.h"
#include "nonlin/ProblemNonLin.h"
#include "nonlin/NonLinSolver.h"
#include "time/RosenbrockAdaptInstationary.h" #include "time/RosenbrockAdaptInstationary.h"
#include "time/RosenbrockStationary.h" #include "time/RosenbrockStationary.h"
#if HAVE_PARALLEL_DOMAIN_AMDIS #if HAVE_PARALLEL_DOMAIN_AMDIS
#include "parallel/InteriorBoundary.h" #include "parallel/InteriorBoundary.h"
#include "parallel/MpiHelper.h" #include "parallel/MpiHelper.h"
......
...@@ -26,6 +26,7 @@ ...@@ -26,6 +26,7 @@
#include "DOFMatrix.h" #include "DOFMatrix.h"
#include "UmfPackSolver.h" #include "UmfPackSolver.h"
#include "time/RosenbrockMethod.h" #include "time/RosenbrockMethod.h"
#include "nonlin/NonLinSolver.h"
namespace AMDiS { namespace AMDiS {
...@@ -142,4 +143,11 @@ namespace AMDiS { ...@@ -142,4 +143,11 @@ namespace AMDiS {
addCreator("ros3p", new Ros3p::Creator); addCreator("ros3p", new Ros3p::Creator);
addCreator("rodasp", new Rodasp::Creator); addCreator("rodasp", new Rodasp::Creator);
} }
template<>
void CreatorMap<NonLinSolver>::addDefaultCreators()
{
addCreator("newton", new Newton::Creator);
}
} }
...@@ -634,8 +634,6 @@ namespace AMDiS { ...@@ -634,8 +634,6 @@ namespace AMDiS {
{ {
FUNCNAME("ProblemStat::buildAfterCoarsen()"); FUNCNAME("ProblemStat::buildAfterCoarsen()");
VtkWriter::writeFile(rhs->getDOFVector(0), "test.vtu");
// buildAfterCoarsen_sebastianMode(adaptInfo, flag); // buildAfterCoarsen_sebastianMode(adaptInfo, flag);
// return; // return;
...@@ -780,7 +778,6 @@ namespace AMDiS { ...@@ -780,7 +778,6 @@ namespace AMDiS {
if (asmMatrix) { if (asmMatrix) {
solverMatrix.setMatrix(*systemMatrix); solverMatrix.setMatrix(*systemMatrix);
createPrecon();
INFO(info, 8)("fillin of assembled matrix: %d\n", nnz); INFO(info, 8)("fillin of assembled matrix: %d\n", nnz);
} }
...@@ -971,8 +968,6 @@ namespace AMDiS { ...@@ -971,8 +968,6 @@ namespace AMDiS {
solverMatrix.setMatrix(*systemMatrix); solverMatrix.setMatrix(*systemMatrix);
createPrecon();
INFO(info, 8)("fillin of assembled matrix: %d\n", nnz); INFO(info, 8)("fillin of assembled matrix: %d\n", nnz);
INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n", INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n",
TIME_USED(first, clock())); TIME_USED(first, clock()));
...@@ -1232,9 +1227,6 @@ namespace AMDiS { ...@@ -1232,9 +1227,6 @@ namespace AMDiS {
solverMatrix.setMatrix(*systemMatrix); solverMatrix.setMatrix(*systemMatrix);
createPrecon();
INFO(info, 8)("fillin of assembled matrix: %d\n", nnz); INFO(info, 8)("fillin of assembled matrix: %d\n", nnz);
INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n", INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n",
...@@ -1242,24 +1234,6 @@ namespace AMDiS { ...@@ -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) void ProblemStatSeq::writeFiles(AdaptInfo *adaptInfo, bool force)
{ {
FUNCNAME("ProblemStat::writeFiles()"); FUNCNAME("ProblemStat::writeFiles()");
...@@ -1570,91 +1544,31 @@ namespace AMDiS { ...@@ -1570,91 +1544,31 @@ namespace AMDiS {
Mesh *mesh = feSpace->getMesh(); Mesh *mesh = feSpace->getMesh();
const BasisFunction *basisFcts = feSpace->getBasisFcts(); const BasisFunction *basisFcts = feSpace->getBasisFcts();
#ifdef _OPENMP
TraverseParallelStack stack(0, 1);
#else
TraverseStack stack; 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 = BoundaryType *bound =
useGetBound ? new BoundaryType[basisFcts->getNumber()] : NULL; useGetBound ? new BoundaryType[basisFcts->getNumber()] : NULL;
// Create for every thread its private matrix and vector, on that matrix->startInsertion(matrix->getNnz());
// the thread will assemble its part of the mesh. vector->set(0.0);
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();
tmpMatrix->getBaseMatrix().change_dim(matrix->getRowFeSpace()->getAdmin()->getUsedSize(),
matrix->getColFeSpace()->getAdmin()->getUsedSize());
tmpMatrix->startInsertion(10);
}
if (vector) {
tmpVector = new DOFVector<double>(vector->getFeSpace(), "tmp");
// 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 mesh and assemble. ==
// traverse only a part of the mesh.
ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag); 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) { while (elInfo) {
if (useGetBound) if (useGetBound)
basisFcts->getBound(elInfo, bound); basisFcts->getBound(elInfo, bound);
if (matrix) { if (matrix) {
tmpMatrix->assemble(1.0, elInfo, bound); matrix->assemble(1.0, elInfo, bound);
// Take the matrix boundary manager from the public matrix, // Take the matrix boundary manager from the public matrix,
// but assemble the boundary conditions on the thread private matrix. // but assemble the boundary conditions on the thread private matrix.
if (matrix->getBoundaryManager()) if (matrix->getBoundaryManager())
matrix->getBoundaryManager()->fillBoundaryConditions(elInfo, tmpMatrix); matrix->getBoundaryManager()->fillBoundaryConditions(elInfo, matrix);
} }
if (vector) if (vector)
tmpVector->assemble(1.0, elInfo, bound, NULL); vector->assemble(1.0, elInfo, bound, NULL);
elInfo = stack.traverseNext(elInfo); elInfo = stack.traverseNext(elInfo);
} }
...@@ -1662,61 +1576,11 @@ namespace AMDiS { ...@@ -1662,61 +1576,11 @@ namespace AMDiS {
// == Finally, if we have assembled in parallel, we have to add the thread == // == Finally, if we have assembled in parallel, we have to add the thread ==
// == private matrix and vector to the global one. == // == 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;
}
#else
if (matrix) if (matrix)
matrix->removeRowsWithDBC(matrix->getApplyDBCs()); matrix->removeRowsWithDBC(matrix->getApplyDBCs());
#endif
if (useGetBound) if (useGetBound)
delete [] bound; delete [] bound;
} // pragma omp parallel
// 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()));
} }
...@@ -1727,6 +1591,7 @@ namespace AMDiS { ...@@ -1727,6 +1591,7 @@ namespace AMDiS {
{ {
FUNCNAME("ProblemStat::assembleBoundaryConditions()"); FUNCNAME("ProblemStat::assembleBoundaryConditions()");
// === Initialization of vectors === // === Initialization of vectors ===
if (rhs->getBoundaryManager()) if (rhs->getBoundaryManager())
...@@ -1734,55 +1599,6 @@ namespace AMDiS { ...@@ -1734,55 +1599,6 @@ namespace AMDiS {
if (solution->getBoundaryManager()) if (solution->getBoundaryManager())
solution->getBoundaryManager()->initVector(solution); 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; TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag); ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag);
while (elInfo) { while (elInfo) {
...@@ -1794,7 +1610,6 @@ namespace AMDiS { ...@@ -1794,7 +1610,6 @@ namespace AMDiS {
elInfo = stack.traverseNext(elInfo); elInfo = stack.traverseNext(elInfo);
} }
#endif
// === Finalize vectors === // === Finalize vectors ===
......
...@@ -48,6 +48,7 @@ namespace AMDiS { ...@@ -48,6 +48,7 @@ namespace AMDiS {
Flag operatorType; Flag operatorType;
}; };
/** \brief /** \brief
* This class defines the stationary problem definition in sequential * This class defines the stationary problem definition in sequential
* computations. For parallel computations, see * computations. For parallel computations, see
...@@ -195,9 +196,6 @@ namespace AMDiS { ...@@ -195,9 +196,6 @@ namespace AMDiS {
void dualAssemble(AdaptInfo *adaptInfo, Flag flag, void dualAssemble(AdaptInfo *adaptInfo, Flag flag,
bool asmMatrix = true, bool asmVector = true); bool asmMatrix = true, bool asmVector = true);
void createPrecon();
/** \brief /** \brief
* Determines the execution order of the single adaption steps. If adapt is * Determines the execution order of the single adaption steps. If adapt is
* true, mesh adaption will be performed. This allows to avoid mesh adaption, * true, mesh adaption will be performed. This allows to avoid mesh adaption,
......
...@@ -26,6 +26,7 @@ ...@@ -26,6 +26,7 @@
#include "CreatorInterface.h" #include "CreatorInterface.h"
#include "NonLinSolver.h" #include "NonLinSolver.h"
#include "OEMSolver.h" #include "OEMSolver.h"
#include "io/VtkWriter.h"
namespace AMDiS { namespace AMDiS {
...@@ -45,114 +46,93 @@ namespace AMDiS { ...@@ -45,114 +46,93 @@ namespace AMDiS {
public: public:
virtual ~Creator() {} virtual ~Creator() {}
/** \brief /// Returns a new Newton object.
* Returns a new Newton object.
*/
NonLinSolver* create() 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
* Calls constructor of base class NonLinSolver Newton(const std::string& name, OEMSolver *linSolver)
*/ : NonLinSolver(name, linSolver),
Newton(const std::string& name_,
OEMSolver *linSolver_,
NonLinUpdater *updater)
: NonLinSolver(name_, linSolver_, updater),
b(NULL) b(NULL)
{} {}
private: private:
/** \brief /// Realisation of NonLinSolver::init
* realisation of NonLinSolver::init void init() {}
*/
void init()
{}
/** \brief /// realisation of NonLinSolver::nlsolve
* realisation of NonLinSolver::nlsolve
*/
int nlsolve(SolverMatrix<Matrix<DOFMatrix*> >& mat, int nlsolve(SolverMatrix<Matrix<DOFMatrix*> >& mat,
SystemVector& x, SystemVector& rhs) SystemVector& x, SystemVector& rhs,
AdaptInfo *adaptInfo,
ProblemStat *prob)
{ {
//DOFVector<T> *d = problem->getRHS(); FUNCNAME("Newton::nlsolve()");
//DOFVector<T> *x = problem->getSolution();;
if (b == NULL)
b = new SystemVector(x); b = new SystemVector(x);
*b = rhs;
double err = 0.0, errOld = -1.0;
// // 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);
// }
double err = 0.0, err_old = -1.0;
int iter, n; 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++) { for (iter = 1; iter <= this->maxIter; iter++) {
/*--- Assemble DF(x) and F(x) ----------------------------------------------*/ // Assemble DF(x) and F(x)
this->nonLinUpdater->update(/*x,*/true, b); prob->buildAfterCoarsen(adaptInfo, 0, true, true);
/*--- Initial guess is zero ------------------------------------------------*/
rhs.set(0.0); // Initial guess is zero
/*--- solve linear system --------------------------------------------------*/ b->set(0.0);
// Solve linear system
n = solveLinearSystem(mat, *b, rhs); n = solveLinearSystem(mat, *b, rhs);
/*--- x = x - d ------------------------------------------------------------*/
x -= rhs; // x = x + d
x += *b;
if (this->usedNorm == NO_NORM || this->usedNorm == L2_NORM) if (this->usedNorm == NO_NORM || this->usedNorm == L2_NORM)
err = L2Norm(&rhs); // sollte hier nicht b genommen werden (s. NewtonS) ? err = L2Norm(b);
else else
err = H1Norm(&rhs); // sollte hier nicht b genommen werden (s. NewtonS) ? err = H1Norm(b);
if (iter == 1) this->initial_residual = err; if (iter == 1)
this->initialResidual = err;
if (err_old <= 0) { if (errOld <= 0)
INFO(this->info,2)("%5d | %12.5e | -------- | %4d |\n", iter, err, n); MSG("%5d | %12.5e | -------- | %4d |\n", iter, err, n);
} else { else
INFO(this->info,2)("%5d | %12.5e | %8.2e | %4d |\n", MSG("%5d | %12.5e | %8.2e | %4d |\n", iter, err, err/errOld, n);
iter, err, err/err_old, n);
}
if ((this->residual = err) < this->tolerance) { residual = err;
INFO(this->info,4)("finished successfully\n"); if (err < this->tolerance) {
MSG("Finished successfully!\n");
return iter; return iter;
} }
err_old = err; errOld = err;
} }
if (this->info < 2)