From f148869d4a823565e07abb159c3e4f1a6344b6be Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Sun, 23 Aug 2015 21:05:11 +0000
Subject: [PATCH] small errors corrected
---
AMDiS/AMDISConfig.cmake.in | 7 +
AMDiS/CMakeLists.txt | 15 +-
AMDiS/src/AMDiS.h | 1 +
AMDiS/src/CouplingProblemStat.h | 3 +-
AMDiS/src/Functors.h | 2 +
AMDiS/src/Global.cc | 24 +--
AMDiS/src/ProblemStat.h | 6 +
AMDiS/src/SystemVector.h | 6 +-
AMDiS/src/config/Config_clang.h | 4 +-
AMDiS/src/config/Config_gcc.h | 4 +-
AMDiS/src/config/Config_intel.h | 4 +-
AMDiS/src/config/Config_msc.h | 4 +-
AMDiS/src/io/GridWriter.hh | 4 +-
AMDiS/src/time/RosenbrockAdaptInstationary.cc | 20 ++-
AMDiS/src/time/RosenbrockAdaptInstationary.h | 10 +-
AMDiS/src/time/RosenbrockStationary.h | 166 ++++++++++--------
extensions/BackgroundMesh.cc | 12 +-
extensions/BoundaryFunctions.h | 4 +-
extensions/Helpers.cc | 20 +--
.../base_problems/CahnHilliardNavierStokes.cc | 2 +-
.../CahnHilliardNavierStokes_RB.cc | 2 +-
.../CahnHilliardNavierStokes_TwoPhase_RB.cc | 4 +-
.../base_problems/LinearElasticityPhase.cc | 2 +-
.../base_problems/NavierStokesCahnHilliard.cc | 2 +-
.../base_problems/PhaseFieldCrystal_Phase.h | 25 ++-
.../base_problems/PhaseFieldCrystal_Phase.hh | 22 +++
extensions/base_problems/QuasiCrystal_RB.cc | 2 +-
extensions/preconditioner/MTLPreconPfc.h | 27 ++-
.../time/ExtendedRosenbrockStationary.cc | 4 +-
.../time/ExtendedRosenbrockStationary.h | 2 +-
30 files changed, 262 insertions(+), 148 deletions(-)
diff --git a/AMDiS/AMDISConfig.cmake.in b/AMDiS/AMDISConfig.cmake.in
index 95f6e667..643aa24c 100644
--- a/AMDiS/AMDISConfig.cmake.in
+++ b/AMDiS/AMDISConfig.cmake.in
@@ -96,6 +96,7 @@ list(APPEND AMDIS_INCLUDE_DIRS ${AMDIS_INCLUDE_DIR}/muparser)
set(AMDIS_NEED_ZOLTAN @ENABLE_ZOLTAN@)
set(AMDIS_OPENMP @ENABLE_OPENMP@)
+set(AMDIS_NEED_CXX11 @ENABLE_CXX11@)
set(AMDIS_NEED_UMFPACK @ENABLE_UMFPACK@)
set(AMDIS_NEED_HYPRE @ENABLE_HYPRE@)
set(AMDIS_NEED_SEQ_PETSC @ENABLE_SEQ_PETSC@)
@@ -114,6 +115,12 @@ if(CMAKE_BUILD_TYPE STREQUAL "")
set(CMAKE_BUILD_TYPE "Release")
endif()
+if (AMDIS_NEED_CXX11)
+ list(APPEND AMDIS_COMPILEFLAGS "-DHAS_CXX11=1" "-std=c++11")
+else ()
+ list(APPEND AMDIS_COMPILEFLAGS "-DHAS_CXX11=0")
+endif (AMDIS_NEED_CXX11)
+
#load mpi-compiler for use with parmetis and parallel_domain
if(AMDIS_HAS_PARALLEL_DOMAIN)
find_package(MPI REQUIRED)
diff --git a/AMDiS/CMakeLists.txt b/AMDiS/CMakeLists.txt
index 333d3e7a..4cecdd49 100644
--- a/AMDiS/CMakeLists.txt
+++ b/AMDiS/CMakeLists.txt
@@ -44,6 +44,7 @@ endif()
SET(ENABLE_PARALLEL_DOMAIN "OFF" CACHE STRING "use parallel domain decomposition. please set to one of: PMTL, PETSC, OFF" )
option(USE_PETSC_DEV false)
option(ENABLE_ZOLTAN false)
+option(ENABLE_CXX11 "Enable C++11 compiler features" false)
option(ENABLE_SEQ_PETSC "Use sequential PETSc solvers" false)
option(ENABLE_UMFPACK "Use of UMFPACK solver" false)
option(ENABLE_HYPRE "Use HYPRE AMG solver" false)
@@ -201,7 +202,7 @@ SET(AMDIS_SRC ${SOURCE_DIR}/AdaptBase.cc
${SOURCE_DIR}/time/RosenbrockMethod.cc
${SOURCE_DIR}/time/RosenbrockStationary.cc
)
-
+
if(ENABLE_PARALLEL_DOMAIN)
### Search for MPI compiler and get compile flags and include paths. ###
@@ -621,6 +622,18 @@ if(WIN32)
endif(WIN32)
+if (ENABLE_CXX11)
+ list(APPEND COMPILEFLAGS "-DHAS_CXX11=1")
+ set_property(TARGET "amdis" PROPERTY CXX_STANDARD 11)
+ set_property(TARGET "compositeFEM" PROPERTY CXX_STANDARD 11)
+ set_property(TARGET "reinit" PROPERTY CXX_STANDARD 11)
+ set_property(TARGET "muparser" PROPERTY CXX_STANDARD 11)
+else ()
+ list(APPEND COMPILEFLAGS "-DHAS_CXX11=0")
+endif (ENABLE_CXX11)
+
+
+
if(CMAKE_SYSTEM_NAME MATCHES Darwin)
if(CMAKE_COMPILER_IS_GNUCXX OR ${CMAKE_CXX_COMPILER_ID} STREQUAL Clang)
SET(CMAKE_SHARED_LIBRARY_CREATE_C_FLAGS "${CMAKE_SHARED_LIBRARY_CREATE_C_FLAGS} -flat_namespace -single_module -undefined dynamic_lookup")
diff --git a/AMDiS/src/AMDiS.h b/AMDiS/src/AMDiS.h
index 5cdfde0f..cdeedbf8 100644
--- a/AMDiS/src/AMDiS.h
+++ b/AMDiS/src/AMDiS.h
@@ -62,6 +62,7 @@
#include "Element.h"
#include "ElementDofIterator.h"
#include "Error.h"
+#include "Expressions.h"
#include "FiniteElemSpace.h"
#include "FirstOrderTerm.h"
#include "FixVec.h"
diff --git a/AMDiS/src/CouplingProblemStat.h b/AMDiS/src/CouplingProblemStat.h
index 3c205d1b..4b214e03 100644
--- a/AMDiS/src/CouplingProblemStat.h
+++ b/AMDiS/src/CouplingProblemStat.h
@@ -30,6 +30,8 @@
#include <list>
#include "AMDiS_fwd.h"
#include "ProblemStat.h"
+#include "CoarseningManager.h"
+#include "RefinementManager.h"
#include "Initfile.h"
#include "utility/to_string.hpp"
@@ -67,7 +69,6 @@ namespace AMDiS {
virtual void addProblem(ProblemStatType* prob)
{
problems.push_back(prob);
- nComponents += prob->getNumComponents();
};
/// Initialisation of the problem.
diff --git a/AMDiS/src/Functors.h b/AMDiS/src/Functors.h
index 8bb9c4cb..0d20973f 100644
--- a/AMDiS/src/Functors.h
+++ b/AMDiS/src/Functors.h
@@ -43,6 +43,8 @@ private:
T1 val;
};
+typedef Const<double, WorldVector<double> > Constant;
+
template<typename T=double>
struct Factor : public AbstractFunction<T,T>
{
diff --git a/AMDiS/src/Global.cc b/AMDiS/src/Global.cc
index 4b433c0d..adf3c6aa 100644
--- a/AMDiS/src/Global.cc
+++ b/AMDiS/src/Global.cc
@@ -74,10 +74,10 @@ namespace AMDiS {
FUNCNAME("Msg::change_out()");
if (fp) {
- if (out && *out != std::cout && *out != std::cerr) {
- dynamic_cast< std::ofstream*>(out)->close();
- delete out;
- }
+// if (out && *out != std::cout && *out != std::cerr) {
+// dynamic_cast< std::ofstream*>(out)->close();
+// delete out;
+// }
out = fp;
} else {
@@ -92,10 +92,10 @@ namespace AMDiS {
FUNCNAME("Msg::change_error_out()");
if (fp) {
- if (error && *error != std::cout && *error != std::cerr) {
- dynamic_cast< std::ofstream*>(error)->close();
- delete error;
- }
+// if (error && *error != std::cout && *error != std::cerr) {
+// dynamic_cast< std::ofstream*>(error)->close();
+// delete error;
+// }
error = fp;
} else {
@@ -111,10 +111,10 @@ namespace AMDiS {
std::ofstream *fp;
if (filename && (fp = new std::ofstream(filename, type))) {
- if (error && *error != std::cout && *error != std::cerr) {
- dynamic_cast< std::ofstream*>(error)->close();
- delete error;
- }
+// if (error && *error != std::cout && *error != std::cerr) {
+// dynamic_cast< std::ofstream*>(error)->close();
+// delete error;
+// }
error = fp;
} else {
diff --git a/AMDiS/src/ProblemStat.h b/AMDiS/src/ProblemStat.h
index 4780a739..77e31619 100644
--- a/AMDiS/src/ProblemStat.h
+++ b/AMDiS/src/ProblemStat.h
@@ -142,6 +142,12 @@ namespace AMDiS {
void buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
bool assembleMatrix = true,
bool assembleVector = true) override;
+
+ /// assemble all operators of matrix and vector side
+ void assemble(AdaptInfo* adaptInfo)
+ {
+ buildAfterCoarsen(adaptInfo, 0, true, true);
+ }
bool dualMeshTraverseRequired();
diff --git a/AMDiS/src/SystemVector.h b/AMDiS/src/SystemVector.h
index 17b4c8b2..f30b4996 100644
--- a/AMDiS/src/SystemVector.h
+++ b/AMDiS/src/SystemVector.h
@@ -343,7 +343,7 @@ namespace AMDiS {
int size = x->getNumVectors();
for (int i = 0; i < size; i++)
result += x->getDOFVector(i)->squareNrm2();
- return sqrt(result);
+ return std::sqrt(result);
}
/// L2 norm of system vector.
@@ -353,7 +353,7 @@ namespace AMDiS {
int size = x->getNumVectors();
for (int i = 0; i < size; i++)
result += x->getDOFVector(i)->L2NormSquare();
- return sqrt(result);
+ return std::sqrt(result);
}
/// H1 norm of system vector.
@@ -363,7 +363,7 @@ namespace AMDiS {
int size = x->getNumVectors();
for (int i = 0; i < size; i++)
result += x->getDOFVector(i)->H1NormSquare();
- return sqrt(result);
+ return std::sqrt(result);
}
inline void mv(const Matrix<DOFMatrix*> &matrix,
diff --git a/AMDiS/src/config/Config_clang.h b/AMDiS/src/config/Config_clang.h
index a76f79fe..22fd55f9 100644
--- a/AMDiS/src/config/Config_clang.h
+++ b/AMDiS/src/config/Config_clang.h
@@ -48,7 +48,9 @@ typedef size_t aligned_size_t __attribute__ ((aligned(CACHE_LINE)));
// C++11 features
// --------------
#if __cplusplus > 199711L
-#define HAS_CXX11 1
+#ifndef HAS_CXX11
+ #define HAS_CXX11 1
+#endif
// __has_feature(cxx_rvalue_references)
#if CLANG_VERSION >= 20900 && !defined(HAS_VARIADIC_TEMPLATES)
diff --git a/AMDiS/src/config/Config_gcc.h b/AMDiS/src/config/Config_gcc.h
index 71550484..4f14aba2 100644
--- a/AMDiS/src/config/Config_gcc.h
+++ b/AMDiS/src/config/Config_gcc.h
@@ -48,7 +48,9 @@ typedef size_t aligned_size_t __attribute__ ((aligned(CACHE_LINE)));
// C++11 features
// --------------
#if __cplusplus > 199711L
-#define HAS_CXX11 1
+#ifndef HAS_CXX11
+ #define HAS_CXX11 1
+#endif
#if GCC_VERSION >= 40300 && !defined(HAS_VARIADIC_TEMPLATES)
#define HAS_VARIADIC_TEMPLATES 1
diff --git a/AMDiS/src/config/Config_intel.h b/AMDiS/src/config/Config_intel.h
index c4623f0f..9569dbb4 100644
--- a/AMDiS/src/config/Config_intel.h
+++ b/AMDiS/src/config/Config_intel.h
@@ -57,7 +57,9 @@ typedef __declspec(align(CACHE_LINE)) size_t aligned_size_t;
#include <functional>
#if defined(_GLIBCXX_TUPLE) || defined(_GLIBCXX_TYPE_TRAITS) || defined(_GLIBCXX_ARRAY) || (__cplusplus > 199711L)
-#define HAS_CXX11 1
+#ifndef HAS_CXX11
+ #define HAS_CXX11 1
+#endif
#if INTEL_VERSION >= 1201 && !defined(HAS_VARIADIC_TEMPLATES)
#define HAS_VARIADIC_TEMPLATES 1
diff --git a/AMDiS/src/config/Config_msc.h b/AMDiS/src/config/Config_msc.h
index 27080c31..dd8428fd 100644
--- a/AMDiS/src/config/Config_msc.h
+++ b/AMDiS/src/config/Config_msc.h
@@ -54,7 +54,9 @@ typedef __declspec(align(CACHE_LINE)) size_t aligned_size_t;
// C++11 features
// --------------
#if __cplusplus > 199711L
-#define HAS_CXX11 1
+#ifndef HAS_CXX11
+ #define HAS_CXX11 1
+#endif
#if MSC_VERSION >= 1800 && !defined(HAS_VARIADIC_TEMPLATES)
#define HAS_VARIADIC_TEMPLATES 1
diff --git a/AMDiS/src/io/GridWriter.hh b/AMDiS/src/io/GridWriter.hh
index 6d50b4e5..f4afb9f6 100644
--- a/AMDiS/src/io/GridWriter.hh
+++ b/AMDiS/src/io/GridWriter.hh
@@ -77,7 +77,7 @@ namespace AMDiS { namespace io {
basis[i][j] = p[i + 1][j] ;
lengthBasis[i] += basis[i][j] * basis[i][j];
}
- lengthBasis[i] = sqrt(lengthBasis[i]);
+ lengthBasis[i] = std::sqrt(lengthBasis[i]);
}
// norm basis, get steps
@@ -191,7 +191,7 @@ namespace AMDiS { namespace io {
basis[i][j] = p[i + 1][j] ;
lengthBasis[i] += basis[i][j] * basis[i][j];
}
- lengthBasis[i] = sqrt(lengthBasis[i]);
+ lengthBasis[i] = std::sqrt(lengthBasis[i]);
}
// norm basis, get steps
diff --git a/AMDiS/src/time/RosenbrockAdaptInstationary.cc b/AMDiS/src/time/RosenbrockAdaptInstationary.cc
index 0101a52b..57bcabd8 100644
--- a/AMDiS/src/time/RosenbrockAdaptInstationary.cc
+++ b/AMDiS/src/time/RosenbrockAdaptInstationary.cc
@@ -26,13 +26,14 @@
namespace AMDiS {
RosenbrockAdaptInstationary::RosenbrockAdaptInstationary(std::string name,
- RosenbrockStationary *problemStat,
+ ProblemIterationInterface *problemStat,
+ RosenbrockInterface *rbprob,
AdaptInfo *info,
ProblemTimeInterface *problemInstat,
AdaptInfo *initialInfo,
time_t initialTimestamp)
: AdaptInstationary(name, problemStat, info, problemInstat, initialInfo, initialTimestamp),
- rosenbrockStat(problemStat),
+ rosenbrockStat(rbprob),
firstTimestep(true),
lastTimestepRejected(false),
succRejection(false),
@@ -45,19 +46,20 @@ namespace AMDiS {
minusInvTauGamma(-1.0),
dbgTimestepStudy(false)
{
- initConstructor(problemStat);
+ initConstructor(rbprob);
rosenbrockStat->setOldTime(adaptInfo->getStartTime());
}
RosenbrockAdaptInstationary::RosenbrockAdaptInstationary(std::string name,
- RosenbrockStationary &problemStat,
+ ProblemIterationInterface &problemStat,
+ RosenbrockInterface &rbprob,
AdaptInfo &info,
ProblemTimeInterface &problemInstat,
AdaptInfo &initialInfo,
time_t initialTimestamp)
: AdaptInstationary(name, problemStat, info, problemInstat, initialInfo, initialTimestamp),
- rosenbrockStat(&problemStat),
+ rosenbrockStat(&rbprob),
firstTimestep(true),
lastTimestepRejected(false),
succRejection(false),
@@ -70,12 +72,12 @@ namespace AMDiS {
minusInvTauGamma(-1.0),
dbgTimestepStudy(false)
{
- initConstructor(&problemStat);
+ initConstructor(&rbprob);
rosenbrockStat->setOldTime(adaptInfo->getStartTime());
}
- void RosenbrockAdaptInstationary::initConstructor(RosenbrockStationary *problemStat)
+ void RosenbrockAdaptInstationary::initConstructor(RosenbrockInterface *problemStat)
{
FUNCNAME_DBG("RosenbrockAdaptInstationary::initConstructor()");
@@ -283,11 +285,11 @@ namespace AMDiS {
{
double errorEst = 0.0;
for (int i = 0; i < adaptInfo->getSize(); i++) {
- double weight = adaptInfo->getTimeTolerance(i) + rosenbrockStat->getSolution(i)->L2Norm() * adaptInfo->getTimeRelativeTolerance(i);
+ double weight = adaptInfo->getTimeTolerance(i) + adaptInfo->getEstSum(i) * adaptInfo->getTimeRelativeTolerance(i);
errorEst += sqr(adaptInfo->getTimeEstCombined(i)) / sqr(weight);
}
- errorEst = sqrt(errorEst / adaptInfo->getSize());
+ errorEst = std::sqrt(errorEst / adaptInfo->getSize());
adaptInfo->setTimeEst(errorEst);
return errorEst;
}
diff --git a/AMDiS/src/time/RosenbrockAdaptInstationary.h b/AMDiS/src/time/RosenbrockAdaptInstationary.h
index 6266f368..a75dc4a9 100644
--- a/AMDiS/src/time/RosenbrockAdaptInstationary.h
+++ b/AMDiS/src/time/RosenbrockAdaptInstationary.h
@@ -42,7 +42,8 @@ namespace AMDiS {
* dependent problem problemInstat. TODO: Make obsolete!
*/
RosenbrockAdaptInstationary(std::string name,
- RosenbrockStationary *problemStat,
+ ProblemIterationInterface *problemStat,
+ RosenbrockInterface *rbprob,
AdaptInfo *info,
ProblemTimeInterface *problemInstat,
AdaptInfo *initialInfo,
@@ -54,7 +55,8 @@ namespace AMDiS {
* dependent problem problemInstat.
*/
RosenbrockAdaptInstationary(std::string name,
- RosenbrockStationary &problemStat,
+ ProblemIterationInterface &problemStat,
+ RosenbrockInterface &rbprob,
AdaptInfo &info,
ProblemTimeInterface &problemInstat,
AdaptInfo &initialInfo,
@@ -71,7 +73,7 @@ namespace AMDiS {
* removed, remove also this function.
* TODO: Remove if obsolete constructor will be removed.
*/
- void initConstructor(RosenbrockStationary *problemStat);
+ void initConstructor(RosenbrockInterface *problemStat);
void reset();
@@ -80,7 +82,7 @@ namespace AMDiS {
RosenbrockMethod *rosenbrockMethod;
/// Pointer to the stationary problem;
- RosenbrockStationary *rosenbrockStat;
+ RosenbrockInterface *rosenbrockStat;
/// Indicates, if this is the very first timestep.
bool firstTimestep;
diff --git a/AMDiS/src/time/RosenbrockStationary.h b/AMDiS/src/time/RosenbrockStationary.h
index 10612a61..95a40689 100644
--- a/AMDiS/src/time/RosenbrockStationary.h
+++ b/AMDiS/src/time/RosenbrockStationary.h
@@ -52,6 +52,88 @@ namespace AMDiS {
int row;
int col;
};
+
+
+ class RosenbrockInterface
+ {
+ public:
+ RosenbrockInterface()
+ : rm(NULL),
+ oldTime(0.0),
+ tauPtr(NULL),
+ tauGamma(NULL),
+ minusTauGamma(NULL),
+ invTauGamma(NULL),
+ minusInvTauGamma(NULL)
+ {}
+
+ virtual ~RosenbrockInterface() {}
+
+ /// update solution vector and oldTime value
+ virtual void acceptTimestep(AdaptInfo* adaptInfo) = 0;
+
+
+ double* getTau()
+ {
+ return tauPtr;
+ }
+
+ double* getTauGamma()
+ {
+ return tauGamma;
+ }
+
+ double* getMinusTauGamma()
+ {
+ return minusTauGamma;
+ }
+
+ double* getInvTauGamma()
+ {
+ return invTauGamma;
+ }
+
+ double* getMinusInvTauGamma()
+ {
+ return minusInvTauGamma;
+ }
+
+ virtual void setRosenbrockMethod(RosenbrockMethod *method)
+ {
+ rm = method;
+ }
+
+ double* getOldTime()
+ {
+ return &oldTime;
+ }
+
+ void setTau(double *ptr)
+ {
+ tauPtr = ptr;
+ }
+
+ void setTauGamma(double *ptr0, double *ptr1, double *ptr2, double *ptr3)
+ {
+ tauGamma = ptr0;
+ minusTauGamma = ptr1;
+ invTauGamma = ptr2;
+ minusInvTauGamma = ptr3;
+ }
+
+ void setOldTime(double t)
+ {
+ oldTime = t;
+ }
+
+ protected:
+ RosenbrockMethod *rm;
+ double oldTime;
+
+ double *tauPtr;
+ double *tauGamma, *minusTauGamma, *invTauGamma, *minusInvTauGamma;
+ };
+
/// Realization of a Rosenbrock time-discretization of M*d_t(X) = F[x]
/**
@@ -65,21 +147,16 @@ namespace AMDiS {
* and new solution
* X^{k+1} = X^k + sum_{j=1}^{s} m_j * Y_j^k
**/
- class RosenbrockStationary : public ProblemStat
+ class RosenbrockStationary : public ProblemStat,
+ public RosenbrockInterface
{
public:
RosenbrockStationary(std::string name, int componentShift_ = 0)
: ProblemStat(name),
- first(true),
- componentShift(componentShift_),
- minusOne(-1.0),
- stageTime(0.0),
- oldTime(0.0),
- tauPtr(NULL),
- tauGamma(NULL),
- minusTauGamma(NULL),
- invTauGamma(NULL),
- minusInvTauGamma(NULL)
+ first(true),
+ componentShift(componentShift_),
+ minusOne(-1.0),
+ stageTime(0.0)
{}
@@ -91,7 +168,7 @@ namespace AMDiS {
virtual void estimateTimeError(AdaptInfo* adaptInfo);
/// update solution vector and oldTime value
- void acceptTimestep(AdaptInfo* adaptInfo);
+ void acceptTimestep(AdaptInfo* adaptInfo) override;
/// Add operators of function F
virtual void addOperator(Operator &op, int row, int col,
@@ -120,42 +197,12 @@ namespace AMDiS {
{
return timeRhsVec->getDOFVector(i);
}
-
- double* getTauGamma()
- {
- return tauGamma;
- }
-
- double* getTau()
- {
- return tauPtr;
- }
-
- double* getMinusTauGamma()
- {
- return minusTauGamma;
- }
-
- double* getInvTauGamma()
- {
- return invTauGamma;
- }
-
- double* getMinusInvTauGamma()
- {
- return minusInvTauGamma;
- }
double* getStageTime()
{
return &stageTime;
}
- double* getOldTime()
- {
- return &oldTime;
- }
-
double* getTauGammaI()
{
return &tauGammaI;
@@ -164,33 +211,15 @@ namespace AMDiS {
// setting methods
// _________________________________________________________________________
- void setRosenbrockMethod(RosenbrockMethod *method)
- {
- rm = method;
- init();
- }
-
- void setTau(double *ptr)
- {
- tauPtr = ptr;
- }
-
- void setTauGamma(double *ptr0, double *ptr1, double *ptr2, double *ptr3)
- {
- tauGamma = ptr0;
- minusTauGamma = ptr1;
- invTauGamma = ptr2;
- minusInvTauGamma = ptr3;
- }
-
- void setOldTime(double t)
+ void setStageTime(double t)
{
- oldTime = t;
+ stageTime = t;
}
- void setStageTime(double t)
+ void setRosenbrockMethod(RosenbrockMethod *method) override
{
- stageTime = t;
+ rm = method;
+ init();
}
@@ -225,11 +254,10 @@ namespace AMDiS {
stageSolution->set(0.0);
unVec->set(0.0);
for (int i = 0; i < rm->getStages(); i++)
- stageSolutions[i]->set(0.0);
+ stageSolutions[i]->set(0.0);
}
protected:
- RosenbrockMethod *rm;
SystemVector *stageSolution, *unVec, *timeRhsVec, *newUn, *tmp, *lowSol;
@@ -241,12 +269,8 @@ namespace AMDiS {
double minusOne;
double stageTime; // t_n + alpha_i*tau
- double oldTime;
double tauGammaI; // tau*gamma_i
- double *tauPtr;
- double *tauGamma, *minusTauGamma, *invTauGamma, *minusInvTauGamma;
-
std::vector<RosenbrockBoundary> boundaries;
};
diff --git a/extensions/BackgroundMesh.cc b/extensions/BackgroundMesh.cc
index 0dc74990..c90036e0 100644
--- a/extensions/BackgroundMesh.cc
+++ b/extensions/BackgroundMesh.cc
@@ -208,7 +208,7 @@ namespace AMDiS { namespace extensions {
std::vector<int> N_;
int N_sum = 0;
for (int i = 0; i < Global::getGeo(WORLD); i++) {
- N_.push_back(boost::math::iround(pow(static_cast<double>(N_const),lengths[i]/sum_lengths)));
+ N_.push_back(boost::math::iround(std::pow(static_cast<double>(N_const),lengths[i]/sum_lengths)));
N_sum += N_[i];
}
@@ -306,8 +306,8 @@ namespace AMDiS { namespace extensions {
break;
case 2:
for (double phi = 0.0; phi < 2.0*m_pi; phi += incr) {
- newIdx[0] = idx[0] + boost::math::iround(radius*cos(phi));
- newIdx[1] = idx[1] + boost::math::iround(radius*sin(phi));
+ newIdx[0] = idx[0] + boost::math::iround(radius*std::cos(phi));
+ newIdx[1] = idx[1] + boost::math::iround(radius*std::sin(phi));
if (newIdx[0] >= 0 && newIdx[0] < N[0]
&& newIdx[1] >= 0 && newIdx[1] < N[1])
newNrs.insert(idx2nr(newIdx));
@@ -316,9 +316,9 @@ namespace AMDiS { namespace extensions {
case 3:
for (double phi = -m_pi; phi < m_pi; phi += incr) {
for (double theta = 0.0; theta < m_pi; theta += incr) {
- newIdx[0] = idx[0] + boost::math::iround(radius*sin(theta)*cos(phi));
- newIdx[1] = idx[1] + boost::math::iround(radius*sin(theta)*sin(phi));
- newIdx[2] = idx[2] + boost::math::iround(radius*cos(theta));
+ newIdx[0] = idx[0] + boost::math::iround(radius*std::sin(theta)*std::cos(phi));
+ newIdx[1] = idx[1] + boost::math::iround(radius*std::sin(theta)*std::sin(phi));
+ newIdx[2] = idx[2] + boost::math::iround(radius*std::cos(theta));
if (newIdx[0] >= 0 && newIdx[0] < N[0]
&& newIdx[1] >= 0 && newIdx[1] < N[1]
&& newIdx[2] >= 0 && newIdx[2] < N[2])
diff --git a/extensions/BoundaryFunctions.h b/extensions/BoundaryFunctions.h
index 8e391251..8432668d 100644
--- a/extensions/BoundaryFunctions.h
+++ b/extensions/BoundaryFunctions.h
@@ -25,8 +25,6 @@
namespace AMDiS { namespace extensions {
-using namespace std;
-
// parabolic inflow
class ParabolicInflow: public AbstractFunction<double, WorldVector<double> >
@@ -254,7 +252,7 @@ class ConstantTimeFct : public AbstractFunction<double, WorldVector<double> >, p
public:
ConstantTimeFct(double c, double endTime_)
: AbstractFunction<double, WorldVector<double> >(1), constant(c), endTime(endTime_) {};
- double operator()(const WorldVector<double>& x) const { return constant*max(0.0, std::min(1.0, (endTime-(*timePtr))/endTime)); }
+ double operator()(const WorldVector<double>& x) const { return constant*std::max(0.0, std::min(1.0, (endTime-(*timePtr))/endTime)); }
protected:
double constant;
diff --git a/extensions/Helpers.cc b/extensions/Helpers.cc
index 9c2e2048..9ff8f67e 100644
--- a/extensions/Helpers.cc
+++ b/extensions/Helpers.cc
@@ -80,19 +80,19 @@ namespace Helpers {
Rz = 1.0;
S = 1.0;
if (dow == 3) {
- Rx[1][1] = cos(rotate[0]); Rx[1][2] = -sin(rotate[0]);
- Rx[2][1] = sin(rotate[0]); Rx[2][2] = cos(rotate[0]);
+ Rx[1][1] = std::cos(rotate[0]); Rx[1][2] = -std::sin(rotate[0]);
+ Rx[2][1] = std::sin(rotate[0]); Rx[2][2] = std::cos(rotate[0]);
- Ry[0][0] = cos(rotate[1]); Ry[0][2] = sin(rotate[1]);
- Ry[2][0] = -sin(rotate[1]); Ry[2][2] = cos(rotate[1]);
+ Ry[0][0] = std::cos(rotate[1]); Ry[0][2] = std::sin(rotate[1]);
+ Ry[2][0] = -std::sin(rotate[1]); Ry[2][2] = std::cos(rotate[1]);
- Rz[0][0] = cos(rotate[2]); Rz[0][1] = -sin(rotate[2]);
- Rz[1][0] = sin(rotate[2]); Rz[1][1] = cos(rotate[2]);
+ Rz[0][0] = std::cos(rotate[2]); Rz[0][1] = -std::sin(rotate[2]);
+ Rz[1][0] = std::sin(rotate[2]); Rz[1][1] = std::cos(rotate[2]);
S[0][0] = scale[0]; S[1][1] = scale[1]; S[2][2] = scale[2];
} else if (dow == 2) {
- Rx[0][0] = cos(rotate[0]); Rx[0][1] = -sin(rotate[0]);
- Rx[1][0] = sin(rotate[0]); Rx[1][1] = cos(rotate[0]);
+ Rx[0][0] = std::cos(rotate[0]); Rx[0][1] = -std::sin(rotate[0]);
+ Rx[1][0] = std::sin(rotate[0]); Rx[1][1] = std::cos(rotate[0]);
Ry = 1.0;
Rz = 1.0;
@@ -425,8 +425,8 @@ namespace Helpers {
{
unsigned num=values.size();
unsigned cols=std::min(numCols,(int)num), rows=std::min(numRows,(int)num);
- unsigned step=(unsigned)floor(num/(double)cols);
- cols=(unsigned)ceil(num/(double)step);
+ unsigned step=(unsigned)std::floor(num/(double)cols);
+ cols=(unsigned)std::ceil(num/(double)step);
double maxVal=*max_element(values.begin(),values.end());
double minVal=*min_element(values.begin(),values.end());
diff --git a/extensions/base_problems/CahnHilliardNavierStokes.cc b/extensions/base_problems/CahnHilliardNavierStokes.cc
index 6681a858..8101570e 100644
--- a/extensions/base_problems/CahnHilliardNavierStokes.cc
+++ b/extensions/base_problems/CahnHilliardNavierStokes.cc
@@ -75,7 +75,7 @@ CahnHilliardNavierStokes::CahnHilliardNavierStokes(const std::string &name_) :
// Parameters for NS-Coupling
// cahn-hiliard-force: sigma*mu*grad(c)
Initfile::get(name + "->sigma", sigma); // surface tension
- surfaceTension = -sigma*3.0/(2.0*sqrt(2.0)) / eps;
+ surfaceTension = -sigma*3.0/(2.0*std::sqrt(2.0)) / eps;
force.set(0.0);
diff --git a/extensions/base_problems/CahnHilliardNavierStokes_RB.cc b/extensions/base_problems/CahnHilliardNavierStokes_RB.cc
index 05ad0b35..669d008a 100644
--- a/extensions/base_problems/CahnHilliardNavierStokes_RB.cc
+++ b/extensions/base_problems/CahnHilliardNavierStokes_RB.cc
@@ -76,7 +76,7 @@ CahnHilliardNavierStokes_RB::CahnHilliardNavierStokes_RB(const std::string &name
// Parameters for NS-Coupling
// cahn-hiliard-force: sigma*mu*grad(c)
Initfile::get(name + "->sigma", sigma); // surface tension
- surfaceTension = -sigma*3.0/(2.0*sqrt(2.0)) / eps;
+ surfaceTension = -sigma*3.0/(2.0*std::sqrt(2.0)) / eps;
force.set(0.0);
diff --git a/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.cc b/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.cc
index 24579b77..2e616b7e 100644
--- a/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.cc
+++ b/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.cc
@@ -74,9 +74,9 @@ CahnHilliardNavierStokes_TwoPhase_RB::CahnHilliardNavierStokes_TwoPhase_RB(const
// Parameters for NS-Coupling
// cahn-hiliard-force: sigma*mu*grad(c)
Initfile::get(name + "->sigma", sigma); // surface tension
- surfaceTension = sigma*3.0/(2.0*sqrt(2.0)) / eps;
+ surfaceTension = sigma*3.0/(2.0*std::sqrt(2.0)) / eps;
Initfile::get(name + "->surfaceTension", surfaceTension); // surface tension
- surfaceTensionFactor = 2.0*sqrt(2.0)*eps/(3.0*sigma);
+ surfaceTensionFactor = 2.0*std::sqrt(2.0)*eps/(3.0*sigma);
force.set(0.0);
diff --git a/extensions/base_problems/LinearElasticityPhase.cc b/extensions/base_problems/LinearElasticityPhase.cc
index 1723149e..39e1d46a 100644
--- a/extensions/base_problems/LinearElasticityPhase.cc
+++ b/extensions/base_problems/LinearElasticityPhase.cc
@@ -69,7 +69,7 @@ void LinearElasticityPhase::initData()
super::initData();
- dbc_factor = beta/pow(epsilon, alpha);
+ dbc_factor = beta/std::pow(epsilon, alpha);
MSG("dbc_factor = %f\n", dbc_factor);
}
diff --git a/extensions/base_problems/NavierStokesCahnHilliard.cc b/extensions/base_problems/NavierStokesCahnHilliard.cc
index e2b127d8..4309bd4b 100644
--- a/extensions/base_problems/NavierStokesCahnHilliard.cc
+++ b/extensions/base_problems/NavierStokesCahnHilliard.cc
@@ -64,7 +64,7 @@ NavierStokesCahnHilliard::NavierStokesCahnHilliard(const std::string &name_, boo
a=1.0/sigma*sqr(eps)*a_fac* (std::max(density1,density2)/tau);
b=1.0;
ab=a*b;
- surfaceTension = -sigma/eps * 3.0/2.0/sqrt(2.0) * a ;
+ surfaceTension = -sigma/eps * 3.0/2.0/std::sqrt(2.0) * a ;
theta1 = 1.0-theta;
minusTheta1 = -theta1;
diff --git a/extensions/base_problems/PhaseFieldCrystal_Phase.h b/extensions/base_problems/PhaseFieldCrystal_Phase.h
index bf17b241..996e6000 100644
--- a/extensions/base_problems/PhaseFieldCrystal_Phase.h
+++ b/extensions/base_problems/PhaseFieldCrystal_Phase.h
@@ -41,25 +41,32 @@ namespace AMDiS { namespace base_problems {
PhaseFieldCrystal_Phase(const std::string &name_);
~PhaseFieldCrystal_Phase();
+
+ void initData() override;
DOFVector<double> *getPhase()
{
- return phase;
+ return phase;
}
-
- void setPhase(DOFVector<double>* phase_)
+
+ void setPhase(DOFVector<double>* phase_, DOFVector<double>* phaseDisturbed_ = NULL)
{
- if (phase == NULL)
- phase = new DOFVector<double>(self::prob->getFeSpace(1), "phase");
- if (phaseDisturbed == NULL)
- phaseDisturbed = new DOFVector<double>(self::prob->getFeSpace(1), "phaseDisturbed");
+ if (phase_)
+ phase = phase_;
+ if (phaseDisturbed_)
+ phaseDisturbed = phaseDisturbed_;
+ }
- phase->interpol(phase_);
- *phaseDisturbed << max(valueOf(phase), 1.e-6);
+ void interpolPhase(DOFVector<double>* phase_)
+ {
+ *phase << valueOf(phase_);
+ *phaseDisturbed << max(valueOf(phase), 1.e-6);
};
virtual void fillOperators() override;
+ void transferInitialSolution(AdaptInfo* adaptInfo) override;
+
double calcEnergy();
private: // variables
diff --git a/extensions/base_problems/PhaseFieldCrystal_Phase.hh b/extensions/base_problems/PhaseFieldCrystal_Phase.hh
index 183fdbf6..db901a51 100644
--- a/extensions/base_problems/PhaseFieldCrystal_Phase.hh
+++ b/extensions/base_problems/PhaseFieldCrystal_Phase.hh
@@ -24,10 +24,22 @@ namespace detail
}
}
+
+ template<typename P>
+ void PhaseFieldCrystal_Phase<P>::initData()
+ {
+ phase = new DOFVector<double>(self::prob->getFeSpace(1), "phase");
+ phaseDisturbed = new DOFVector<double>(self::prob->getFeSpace(1), "phaseDisturbed");
+
+ super::initData();
+ }
+
template<typename P>
void PhaseFieldCrystal_Phase<P>::fillOperators()
{
+ assert( phase && phaseDisturbed );
+
const FiniteElemSpace* feSpace0 = self::prob->getFeSpace(0);
const FiniteElemSpace* feSpace1 = self::prob->getFeSpace(1);
const FiniteElemSpace* feSpace2 = self::prob->getFeSpace(2);
@@ -85,6 +97,16 @@ namespace detail
addZOT(op22, valueOf(phaseDisturbed));
self::prob->addMatrixOperator(op22, 2, 2); // nu
}
+
+
+ template<typename P>
+ void PhaseFieldCrystal_Phase<P>::transferInitialSolution(AdaptInfo* adaptInfo)
+ {
+ super::transferInitialSolution(adaptInfo);
+
+ assert( phase );
+ *phaseDisturbed << max(valueOf(phase), 1.e-5);
+ }
template<typename P>
diff --git a/extensions/base_problems/QuasiCrystal_RB.cc b/extensions/base_problems/QuasiCrystal_RB.cc
index f95df93f..9c26feb6 100644
--- a/extensions/base_problems/QuasiCrystal_RB.cc
+++ b/extensions/base_problems/QuasiCrystal_RB.cc
@@ -6,7 +6,7 @@ using namespace std;
QuasiCrystal_RB::QuasiCrystal_RB(const std::string &name_) :
super(name_),
- q(2.0*cos(m_pi/12.0)),
+ q(2.0*std::cos(m_pi/12.0)),
r(-0.015),
c(50.0),
minusC(-50.0),
diff --git a/extensions/preconditioner/MTLPreconPfc.h b/extensions/preconditioner/MTLPreconPfc.h
index 508876b8..7b56f812 100644
--- a/extensions/preconditioner/MTLPreconPfc.h
+++ b/extensions/preconditioner/MTLPreconPfc.h
@@ -104,9 +104,16 @@ struct MTLPreconPfcBase : BlockPreconditioner<MatrixType>
if (useDirectSolver && solverM)
(*solverM)(x,b);
else {
+ std::string solver = "cg";
+ Parameters::get("precon_pfc_M->solver", solver);
itl::basic_iteration<double> iter(b, nIterM, rtolM, tolM);
x = 0.0;
- itl::cg(getM(), x, b, *PDiagM, *PId, iter);
+ if (solver == "cg")
+ itl::cg(getM(), x, b, *PDiagM, *PId, iter);
+ else if (solver == "bicgstab")
+ itl::bicgstab(getM(), x, b, *PDiagM, iter);
+ else if (solver == "bicgstab_ell")
+ itl::bicgstab_ell(getM(), x, b, *PDiagM, *PId, iter, 3);
}
}
@@ -123,9 +130,16 @@ struct MTLPreconPfcBase : BlockPreconditioner<MatrixType>
if (useDirectSolver && solverMpL)
(*solverMpL)(x,b);
else {
+ std::string solver = "cg";
+ Parameters::get("precon_pfc_MpL->solver", solver);
itl::basic_iteration<double> iter(b, nIterMpL, rtolMpL, tolMpL);
x = 0.0;
- itl::cg(MpL, x, b, *PDiagMpL, *PId, iter);
+ if (solver == "cg")
+ itl::cg(MpL, x, b, *PDiagMpL, *PId, iter);
+ else if (solver == "bicgstab")
+ itl::bicgstab(MpL, x, b, *PDiagMpL, iter);
+ else if (solver == "bicgstab_ell")
+ itl::bicgstab_ell(MpL, x, b, *PDiagMpL, *PId, iter, 3);
}
}
@@ -142,9 +156,16 @@ struct MTLPreconPfcBase : BlockPreconditioner<MatrixType>
if (useDirectSolver && solverMpL2)
(*solverMpL2)(x,b);
else {
+ std::string solver = "cg";
+ Parameters::get("precon_pfc_MpL2->solver", solver);
itl::basic_iteration<double> iter(b, nIterMpL2, rtolMpL2, tolMpL2);
x = 0.0;
- itl::cg(MpL2, x, b, *PDiagMpL2, *PId, iter);
+ if (solver == "cg")
+ itl::cg(MpL2, x, b, *PDiagMpL2, *PId, iter);
+ else if (solver == "bicgstab")
+ itl::bicgstab(MpL2, x, b, *PDiagMpL2, iter);
+ else if (solver == "bicgstab_ell")
+ itl::bicgstab_ell(MpL2, x, b, *PDiagMpL2, *PId, iter, 3);
}
}
diff --git a/extensions/time/ExtendedRosenbrockStationary.cc b/extensions/time/ExtendedRosenbrockStationary.cc
index 8afa4c83..5dc7118c 100644
--- a/extensions/time/ExtendedRosenbrockStationary.cc
+++ b/extensions/time/ExtendedRosenbrockStationary.cc
@@ -223,7 +223,7 @@ namespace AMDiS { namespace extensions {
double ExtendedRosenbrockStationary::getNewTimestep(double tol, double tau, bool restrict)
{
- double fac = pow((tol / errorEst), 1.0 / order);
+ double fac = std::pow((tol / errorEst), 1.0 / order);
if (restrict)
fac = std::min(3.0, fac);
return 0.95 * (fac * tau);
@@ -231,7 +231,7 @@ namespace AMDiS { namespace extensions {
double ExtendedRosenbrockStationary::getNewTimestep(double tol, double oldTau, double tau, bool restrict)
{
- double fac = tau / oldTau * pow((tol * oldErrorEst / sqr(errorEst)), 1.0 / order);
+ double fac = tau / oldTau * std::pow((tol * oldErrorEst / sqr(errorEst)), 1.0 / order);
if (restrict)
fac = std::min(3.0, fac);
return 0.95 * (fac * tau);
diff --git a/extensions/time/ExtendedRosenbrockStationary.h b/extensions/time/ExtendedRosenbrockStationary.h
index d01ccabb..7d2ff005 100644
--- a/extensions/time/ExtendedRosenbrockStationary.h
+++ b/extensions/time/ExtendedRosenbrockStationary.h
@@ -130,7 +130,7 @@ namespace AMDiS { namespace extensions {
void reduceOrder(double tauRej, double tau)
{
- double newOrder = log(errorEst / oldErrorEst) / log(tau / tauRej);
+ double newOrder = std::log(errorEst / oldErrorEst) / std::log(tau / tauRej);
if (newOrder < order && newOrder > 0.0)
order = newOrder;
}
--
GitLab