Commit 615a8b2c authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

First general support for Rosenbrock methods for time discretization.

parent bb398e1b
......@@ -62,6 +62,7 @@ SET(AMDIS_SRC ${SOURCE_DIR}/DOFIndexed.cc ${SOURCE_DIR}/GNUPlotWriter.cc ${SOURC
${SOURCE_DIR}/PovrayWriter.cc ${SOURCE_DIR}/DataCollector.cc ${SOURCE_DIR}/ScalableQuadrature.cc
${SOURCE_DIR}/SubElInfo.cc ${SOURCE_DIR}/parallel/InteriorBoundary.cc ${SOURCE_DIR}/ElementDofIterator.cc
${SOURCE_DIR}/Debug.cc
${SOURCE_DIR}/time/RosenbrockAdaptInstationary.cc ${SOURCE_DIR}/time/RosenbrockStationary.cc
)
if(ENABLE_PARMETIS)
......
......@@ -223,7 +223,9 @@ $(SOURCE_DIR)/SubElInfo.h $(SOURCE_DIR)/SubElInfo.cc \
$(SOURCE_DIR)/SolutionDataStorage.h $(SOURCE_DIR)/SolutionDataStorage.hh \
$(SOURCE_DIR)/ElementDofIteartor.h $(SOURCE_DIR)/ElementDofIterator.cc \
$(SOURCE_DIR)/parallel/InteriorBoundary.h $(SOURCE_DIR)/parallel/InteriorBoundary.cc \
$(SOURCE_DIR)/Debug.h $(SOURCE_DIR)/Debug.cc
$(SOURCE_DIR)/Debug.h $(SOURCE_DIR)/Debug.cc \
$(SOURCE_DIR)/time/RosenbrockAdaptInstationary.h $(SOURCE_DIR)/time/RosenbrockAdaptInstationary.cc \
$(SOURCE_DIR)/time/RosenbrockStationary.h $(SOURCE_DIR)/time/RosenbrockStationary.cc
COMPOSITE_SOURCE_DIR = ../compositeFEM/src
......
......@@ -239,7 +239,11 @@ am__libamdis_la_SOURCES_DIST = $(SOURCE_DIR)/parallel/StdMpi.h \
$(SOURCE_DIR)/ElementDofIterator.cc \
$(SOURCE_DIR)/parallel/InteriorBoundary.h \
$(SOURCE_DIR)/parallel/InteriorBoundary.cc \
$(SOURCE_DIR)/Debug.h $(SOURCE_DIR)/Debug.cc
$(SOURCE_DIR)/Debug.h $(SOURCE_DIR)/Debug.cc \
$(SOURCE_DIR)/time/RosenbrockAdaptInstationary.h \
$(SOURCE_DIR)/time/RosenbrockAdaptInstationary.cc \
$(SOURCE_DIR)/time/RosenbrockStationary.h \
$(SOURCE_DIR)/time/RosenbrockStationary.cc
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@am__objects_1 = libamdis_la-StdMpi.lo \
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ libamdis_la-ParallelDomainBase.lo \
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ libamdis_la-ParallelDomainDbg.lo \
......@@ -306,7 +310,9 @@ am_libamdis_la_OBJECTS = $(am__objects_2) libamdis_la-DOFIndexed.lo \
libamdis_la-PngWriter.lo libamdis_la-PovrayWriter.lo \
libamdis_la-DataCollector.lo libamdis_la-ScalableQuadrature.lo \
libamdis_la-SubElInfo.lo libamdis_la-ElementDofIterator.lo \
libamdis_la-InteriorBoundary.lo libamdis_la-Debug.lo
libamdis_la-InteriorBoundary.lo libamdis_la-Debug.lo \
libamdis_la-RosenbrockAdaptInstationary.lo \
libamdis_la-RosenbrockStationary.lo
libamdis_la_OBJECTS = $(am_libamdis_la_OBJECTS)
libcompositeFEM_la_LIBADD =
am_libcompositeFEM_la_OBJECTS = libcompositeFEM_la-CFE_Integration.lo \
......@@ -650,7 +656,9 @@ $(SOURCE_DIR)/SubElInfo.h $(SOURCE_DIR)/SubElInfo.cc \
$(SOURCE_DIR)/SolutionDataStorage.h $(SOURCE_DIR)/SolutionDataStorage.hh \
$(SOURCE_DIR)/ElementDofIteartor.h $(SOURCE_DIR)/ElementDofIterator.cc \
$(SOURCE_DIR)/parallel/InteriorBoundary.h $(SOURCE_DIR)/parallel/InteriorBoundary.cc \
$(SOURCE_DIR)/Debug.h $(SOURCE_DIR)/Debug.cc
$(SOURCE_DIR)/Debug.h $(SOURCE_DIR)/Debug.cc \
$(SOURCE_DIR)/time/RosenbrockAdaptInstationary.h $(SOURCE_DIR)/time/RosenbrockAdaptInstationary.cc \
$(SOURCE_DIR)/time/RosenbrockStationary.h $(SOURCE_DIR)/time/RosenbrockStationary.cc
COMPOSITE_SOURCE_DIR = ../compositeFEM/src
libcompositeFEM_la_CXXFLAGS = $(libamdis_la_CXXFLAGS)
......@@ -830,6 +838,8 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-RefinementManager3d.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-ResidualEstimator.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-RobinBC.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-RosenbrockAdaptInstationary.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-RosenbrockStationary.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-ScalableQuadrature.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-SecondOrderAssembler.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-SecondOrderTerm.Plo@am__quote@
......@@ -1615,6 +1625,20 @@ libamdis_la-Debug.lo: $(SOURCE_DIR)/Debug.cc
@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCXX_FALSE@ $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -c -o libamdis_la-Debug.lo `test -f '$(SOURCE_DIR)/Debug.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/Debug.cc
libamdis_la-RosenbrockAdaptInstationary.lo: $(SOURCE_DIR)/time/RosenbrockAdaptInstationary.cc
@am__fastdepCXX_TRUE@ if $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -MT libamdis_la-RosenbrockAdaptInstationary.lo -MD -MP -MF "$(DEPDIR)/libamdis_la-RosenbrockAdaptInstationary.Tpo" -c -o libamdis_la-RosenbrockAdaptInstationary.lo `test -f '$(SOURCE_DIR)/time/RosenbrockAdaptInstationary.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/time/RosenbrockAdaptInstationary.cc; \
@am__fastdepCXX_TRUE@ then mv -f "$(DEPDIR)/libamdis_la-RosenbrockAdaptInstationary.Tpo" "$(DEPDIR)/libamdis_la-RosenbrockAdaptInstationary.Plo"; else rm -f "$(DEPDIR)/libamdis_la-RosenbrockAdaptInstationary.Tpo"; exit 1; fi
@AMDEP_TRUE@@am__fastdepCXX_FALSE@ source='$(SOURCE_DIR)/time/RosenbrockAdaptInstationary.cc' object='libamdis_la-RosenbrockAdaptInstationary.lo' libtool=yes @AMDEPBACKSLASH@
@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCXX_FALSE@ $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -c -o libamdis_la-RosenbrockAdaptInstationary.lo `test -f '$(SOURCE_DIR)/time/RosenbrockAdaptInstationary.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/time/RosenbrockAdaptInstationary.cc
libamdis_la-RosenbrockStationary.lo: $(SOURCE_DIR)/time/RosenbrockStationary.cc
@am__fastdepCXX_TRUE@ if $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -MT libamdis_la-RosenbrockStationary.lo -MD -MP -MF "$(DEPDIR)/libamdis_la-RosenbrockStationary.Tpo" -c -o libamdis_la-RosenbrockStationary.lo `test -f '$(SOURCE_DIR)/time/RosenbrockStationary.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/time/RosenbrockStationary.cc; \
@am__fastdepCXX_TRUE@ then mv -f "$(DEPDIR)/libamdis_la-RosenbrockStationary.Tpo" "$(DEPDIR)/libamdis_la-RosenbrockStationary.Plo"; else rm -f "$(DEPDIR)/libamdis_la-RosenbrockStationary.Tpo"; exit 1; fi
@AMDEP_TRUE@@am__fastdepCXX_FALSE@ source='$(SOURCE_DIR)/time/RosenbrockStationary.cc' object='libamdis_la-RosenbrockStationary.lo' libtool=yes @AMDEPBACKSLASH@
@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCXX_FALSE@ $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -c -o libamdis_la-RosenbrockStationary.lo `test -f '$(SOURCE_DIR)/time/RosenbrockStationary.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/time/RosenbrockStationary.cc
libcompositeFEM_la-CFE_Integration.lo: $(COMPOSITE_SOURCE_DIR)/CFE_Integration.cc
@am__fastdepCXX_TRUE@ if $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libcompositeFEM_la_CXXFLAGS) $(CXXFLAGS) -MT libcompositeFEM_la-CFE_Integration.lo -MD -MP -MF "$(DEPDIR)/libcompositeFEM_la-CFE_Integration.Tpo" -c -o libcompositeFEM_la-CFE_Integration.lo `test -f '$(COMPOSITE_SOURCE_DIR)/CFE_Integration.cc' || echo '$(srcdir)/'`$(COMPOSITE_SOURCE_DIR)/CFE_Integration.cc; \
@am__fastdepCXX_TRUE@ then mv -f "$(DEPDIR)/libcompositeFEM_la-CFE_Integration.Tpo" "$(DEPDIR)/libcompositeFEM_la-CFE_Integration.Plo"; else rm -f "$(DEPDIR)/libcompositeFEM_la-CFE_Integration.Tpo"; exit 1; fi
......
......@@ -97,6 +97,9 @@
#include "VtkWriter.h"
#include "ZeroOrderTerm.h"
#include "time/RosenbrockAdaptInstationary.h"
#include "time/RosenbrockStationary.h"
#if HAVE_PARALLEL_DOMAIN_AMDIS
#include "parallel/InteriorBoundary.h"
#include "parallel/MpiHelper.h"
......
......@@ -84,6 +84,8 @@ namespace AMDiS {
class Q1Psi;
class RCNeighbourList;
class RefinementManager;
class RosenbrockAdaptInstationary;
class RosenbrockStationary;
class RobinBC;
class SubElInfo;
class SurfaceOperator;
......
#include "time/RosenbrockAdaptInstationary.h"
#include "AdaptInfo.h"
#include "ProblemTimeInterface.h"
namespace AMDiS {
void RosenbrockAdaptInstationary::oneTimestep()
{
FUNCNAME("RosenbrockAdaptInstationary::oneTimestep()");
// estimate before first adaption
if (adaptInfo->getTime() <= adaptInfo->getStartTime())
problemIteration->oneIteration(adaptInfo, ESTIMATE);
bool rejected = false;
do {
// increment time
adaptInfo->setTime(adaptInfo->getTime() + adaptInfo->getTimestep());
problemTime->setTime(adaptInfo);
INFO(info, 6)("time = %e timestep = %e\n",
adaptInfo->getTime(), adaptInfo->getTimestep());
adaptInfo->setSpaceIteration(0);
// do the iteration
problemIteration->beginIteration(adaptInfo);
problemIteration->oneIteration(adaptInfo, FULL_ITERATION);
problemIteration->endIteration(adaptInfo);
double errorEst = adaptInfo->getTimeEstSum(0);
double newTimestep = 0.0;
if (errorEst < timeTol) {
double fac = 1.0;
if (firstTimestep || succRejection) {
fac = pow((timeTol / errorEst), 1.0 / 3.0);
} else {
fac = adaptInfo->getTimestep() / tauAcc *
pow((timeTol * estAcc / (errorEst * errorEst)), 1.0 / 3.0);
}
fac = std::min(fac, 3.0);
newTimestep = fac * adaptInfo->getTimestep();
tauAcc = adaptInfo->getTimestep();
estAcc = errorEst;
lastTimestepRejected = false;
succRejection = false;
} else {
double p = 3.0;
if (lastTimestepRejected) {
succRejection = true;
double reducedP = log(errorEst / estRej) / log(adaptInfo->getTimestep() / tauRej);
if (reducedP < p && reducedP > 0.0)
p = reducedP;
}
newTimestep = pow((timeTol / errorEst), 1.0 / p) * adaptInfo->getTimestep();
tauRej = adaptInfo->getTimestep();
estRej = errorEst;
lastTimestepRejected = true;
}
newTimestep *= 0.95;
if (errorEst < timeTol || firstTimestep) {
INFO(info, 6)("Accepted timestep at time = %f with timestep = %f\n",
adaptInfo->getTime() - adaptInfo->getTimestep(),
adaptInfo->getTimestep());
adaptInfo->setTimestep(newTimestep);
if (firstTimestep)
firstTimestep = false;
rejected = false;
for (int i = 0; i < 3; i++)
INFO(info, 6)("time estimator for component %d = %f\n",
i, adaptInfo->getTimeEstSum(i));
} else {
INFO(info, 6)("Rejected timestep at time = %f with timestep = %f\n",
adaptInfo->getTime() - adaptInfo->getTimestep(),
adaptInfo->getTimestep());
adaptInfo->setTime(adaptInfo->getTime() - adaptInfo->getTimestep());
adaptInfo->setTimestep(newTimestep);
rejected = true;
}
} while (rejected);
rosenbrockStat->acceptTimestep();
adaptInfo->setLastProcessedTimestep(adaptInfo->getTimestep());
adaptInfo->incTimestepNumber();
}
}
// ============================================================================
// == ==
// == AMDiS - Adaptive multidimensional simulations ==
// == ==
// ============================================================================
// == ==
// == TU Dresden ==
// == ==
// == Institut fr Wissenschaftliches Rechnen ==
// == Zellescher Weg 12-14 ==
// == 01069 Dresden ==
// == germany ==
// == ==
// ============================================================================
// == ==
// == https://gforge.zih.tu-dresden.de/projects/amdis/ ==
// == ==
// ============================================================================
/** \file RosenbrockAdaptInstationary.h */
#ifndef AMDIS_ROSENBROCKADAPTINSTATIONARY_H
#define AMDIS_ROSENBROCKADAPTINSTATIONARY_H
#include "AMDiS_fwd.h"
#include "time/RosenbrockStationary.h"
#include "AdaptInstationary.h"
namespace AMDiS {
class RosenbrockAdaptInstationary : public AdaptInstationary
{
public:
RosenbrockAdaptInstationary(std::string name,
RosenbrockStationary *problemStat,
AdaptInfo *info,
ProblemTimeInterface *problemInstat,
AdaptInfo *initialInfo,
time_t initialTimestamp = 0)
: AdaptInstationary(name, problemStat, info, problemInstat, initialInfo, initialTimestamp),
rosenbrockStat(problemStat),
firstTimestep(true),
lastTimestepRejected(false),
succRejection(false)
{
timeTol = 0.0;
GET_PARAMETER(1, "pfc->time tol", "%f", &timeTol);
}
void oneTimestep();
protected:
RosenbrockStationary *rosenbrockStat;
double timeTol;
bool firstTimestep, lastTimestepRejected, succRejection;
double tauAcc, estAcc, tauRej, estRej;
};
}
#endif // AMDIS_ROSENBROCKADAPTINSTATIONARY_H
#include "time/RosenbrockStationary.h"
#include "ProblemVec.h"
#include "SystemVector.h"
#include "OEMSolver.h"
namespace AMDiS {
void RosenbrockStationary::acceptTimestep()
{
*solution = *newUn;
*unVec = *newUn;
}
void RosenbrockStationary::initRB()
{
stages = 3;
stageSolution = new SystemVector(*solution);
unVec = new SystemVector(*solution);
newUn = new SystemVector(*solution);
tmp = new SystemVector(*solution);
lowSol = new SystemVector(*solution);
stageSolutions.resize(stages);
for (int i = 0; i < stages; i++)
stageSolutions[i] = new SystemVector(*solution);
phiSum = new DOFVector<double>(feSpaces[0], "phiSum");
tmpDof = new DOFVector<double>(feSpaces[0], "phiSum");
a.resize(stages);
for (int i = 0; i < stages; i++) {
a[i].resize(stages);
for (int j = 0; j < stages; j++)
a[i][j] = 0.0;
}
a[0][0] = 0.0;
a[1][0] = 1.605996252195329e+00;
a[1][1] = 0.0;
a[2][0] = 1.605996252195329e+00;
a[2][1] = 0.0;
a[2][2] = 0.0;
c.resize(stages);
for (int i = 0; i < stages; i++) {
c[i].resize(stages);
for (int j = 0; j < stages; j++)
c[i][j] = 0.0;
}
c[0][0] = 2.294280360279042e+00;
c[1][0] = -8.874044410657833e-01;
c[1][1] = 2.294280360279042e+00;
c[2][0] = -2.398747971635036e+01;
c[2][1] = -5.263722371562129e+00;
c[2][2] = 2.294280360279042e+00;
m.resize(stages);
m[0] = 2.236727045296590e+00;
m[1] = 2.250067730969644e+00;
m[2] = -2.092514044390320e-01;
m2.resize(stages);
m2[0] = 2.059356167645940e+00;
m2[1] = 1.694014319346528e-01;
m2[2] = 0.0;
}
void RosenbrockStationary::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
bool asmMatrix, bool asmVector)
{
if (first) {
first = false;
*unVec = *solution;
}
*newUn = *unVec;
*lowSol = *unVec;
for (int i = 0; i < stages; i++) {
*stageSolution = *unVec;
for (int j = 0; j < i; j++) {
*tmp = *(stageSolutions[j]);
*tmp *= a[i][j];
*stageSolution += *tmp;
}
phiSum->set(0.0);
for (int j = 0; j < i; j++) {
*tmpDof = *(stageSolutions[j]->getDOFVector(0));
*tmpDof *= (c[i][j] / *tauPtr);
*phiSum += *tmpDof;
}
ProblemVec::buildAfterCoarsen(adaptInfo, flag, (i == 0), asmVector);
solver->setMultipleRhs(i != 0);
ProblemVec::solve(adaptInfo);
*(stageSolutions[i]) = *solution;
*tmp = *solution;
*tmp *= m[i];
*newUn += *tmp;
*tmp = *solution;
*tmp *= m2[i];
*lowSol += *tmp;
}
for (int i = 0; i < nComponents; i++) {
(*(lowSol->getDOFVector(i))) -= (*(newUn->getDOFVector(i)));
adaptInfo->setTimeEstSum(lowSol->getDOFVector(i)->l2norm(), i);
}
}
void RosenbrockStationary::solve(AdaptInfo *adaptInfo, bool fixedMatrix)
{}
}
// ============================================================================
// == ==
// == AMDiS - Adaptive multidimensional simulations ==
// == ==
// ============================================================================
// == ==
// == TU Dresden ==
// == ==
// == Institut fr Wissenschaftliches Rechnen ==
// == Zellescher Weg 12-14 ==
// == 01069 Dresden ==
// == germany ==
// == ==
// ============================================================================
// == ==
// == https://gforge.zih.tu-dresden.de/projects/amdis/ ==
// == ==
// ============================================================================
/** \file RosenbrockStationary.h */
#ifndef AMDIS_ROSENBROCKSTATIONARY_H
#define AMDIS_ROSENBROCKSTATIONARY_H
#include "AMDiS_fwd.h"
#include "ProblemVec.h"
#include "SystemVector.h"
namespace AMDiS {
class RosenbrockStationary : public ProblemVec
{
public:
RosenbrockStationary(std::string name)
: ProblemVec(name),
first(true)
{}
void initRB();
DOFVector<double>* getStageSolution(int i)
{
return stageSolution->getDOFVector(i);
}
DOFVector<double>* getUnVec(int i)
{
return unVec->getDOFVector(i);
}
DOFVector<double>* getPhiSum()
{
return phiSum;
}
void setTauPtr(double *ptr)
{
tauPtr = ptr;
}
void acceptTimestep();
void buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
bool asmMatrix, bool asmVector);
void solve(AdaptInfo *adaptInfo, bool fixedMatrix);
protected:
int stages;
SystemVector *stageSolution, *unVec, *newUn, *tmp, *lowSol;
std::vector<SystemVector*> stageSolutions;
DOFVector<double> *phiSum, *tmpDof;
std::vector<std::vector<double> > a, c;
std::vector<double> m, m2;
double *tauPtr;
bool first;
};
}
#endif // AMDIS_ROSENBROCKSTATIONARY_H
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment