Skip to content

GitLab

Commit 3dc9cf95 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

* OpenMP bugfixes for vecphase program

parent 244db4c9
Hide whitespace changes
Inline Side-by-side
AMDiS/src/AbstractFunction.h
View file @ 3dc9cf95
......@@ -49,19 +49,23 @@ namespace AMDiS {
/** \brief
* Constructor.
*/
AbstractFunction(int degree = 0) : degree_(degree) {};
AbstractFunction(int degree = 0) :
degree_(degree)
{};
virtual ~AbstractFunction() {};
/** \brief
* Returns \ref degree_.
*/
inline int getDegree() const { return degree_; };
inline int getDegree() const {
return degree_;
};
/** \brief
* Deligates the evaluation to overriden method f.
*/
virtual const ReturnType& operator()(const ArgumentType& x) const = 0;
virtual ReturnType operator()(const ArgumentType& x) const = 0;
protected:
int degree_;
......@@ -86,7 +90,9 @@ namespace AMDiS {
/** \brief
* Constructor.
*/
BinaryAbstractFunction(int degree = 0) : degree_(degree) {};
BinaryAbstractFunction(int degree = 0) :
degree_(degree)
{};
virtual ~BinaryAbstractFunction() {};
......@@ -100,8 +106,8 @@ namespace AMDiS {
/** \brief
* Deligates the evaluation to overriden method f.
*/
virtual const ReturnType& operator()(const ArgumentType1& x,
const ArgumentType2& y) const = 0;
virtual ReturnType operator()(const ArgumentType1& x,
const ArgumentType2& y) const = 0;
protected:
int degree_;
......@@ -127,7 +133,9 @@ namespace AMDiS {
/** \brief
* Constructor.
*/
TertiaryAbstractFunction(int degree = 0) : degree_(degree) {};
TertiaryAbstractFunction(int degree = 0) :
degree_(degree)
{};
virtual ~TertiaryAbstractFunction() {};
......@@ -141,9 +149,9 @@ namespace AMDiS {
/** \brief
* function evaluation.
*/
virtual const ReturnType& operator()(const ArgumentType1& x,
const ArgumentType2& y,
const ArgumentType3& z) const = 0;
virtual ReturnType operator()(const ArgumentType1& x,
const ArgumentType2& y,
const ArgumentType3& z) const = 0;
protected:
int degree_;
......@@ -171,7 +179,9 @@ namespace AMDiS {
/** \brief
* Constructor.
*/
QuartAbstractFunction(int degree = 0) : degree_(degree) {};
QuartAbstractFunction(int degree = 0) :
degree_(degree)
{};
virtual ~QuartAbstractFunction() {};
......@@ -185,10 +195,10 @@ namespace AMDiS {
/** \brief
* function evaluation.
*/
virtual const ReturnType& operator()(const ArgumentType1& x,
const ArgumentType2& y,
const ArgumentType3& z,
const ArgumentType4& u) const = 0;
virtual ReturnType operator()(const ArgumentType1& x,
const ArgumentType2& y,
const ArgumentType3& z,
const ArgumentType4& u) const = 0;
protected:
int degree_;
......
AMDiS/src/Assembler.h
View file @ 3dc9cf95
......@@ -37,6 +37,7 @@
#include "FirstOrderAssembler.h"
#include "SecondOrderAssembler.h"
#include "ElInfo.h"
#include "OpenMP.h"
namespace AMDiS {
......@@ -218,9 +219,9 @@ namespace AMDiS {
* Checks whether this is a new travese.
*/
inline void checkForNewTraverse() {
if (lastTraverseId != ElInfo::traverseId) {
if (lastTraverseId != ElInfo::traverseId[omp_get_thread_num()]) {
lastVecEl = lastMatEl = NULL;
lastTraverseId = ElInfo::traverseId;
lastTraverseId = ElInfo::traverseId[omp_get_thread_num()];
}
};
......
AMDiS/src/DOFAdmin.cc
View file @ 3dc9cf95
......@@ -15,13 +15,6 @@ namespace AMDiS {
const int DOFAdmin::sizeIncrement = 10;
void DOFAdmin::init()
{
firstHole = size = usedCount = holeCount = sizeUsed = 0;
dofFree.clear();
}
DOFAdmin::DOFAdmin(Mesh* m)
: mesh(m),
nrDOF(mesh->getDim(), NO_INIT),
......@@ -39,6 +32,12 @@ namespace AMDiS {
init();
}
void DOFAdmin::init()
{
firstHole = size = usedCount = holeCount = sizeUsed = 0;
dofFree.clear();
}
DOFAdmin& DOFAdmin::operator=(const DOFAdmin& src)
{
if (this != &src) {
......@@ -75,7 +74,8 @@ namespace AMDiS {
DOFAdmin::DOFAdmin(const DOFAdmin&)
{}
{
}
void DOFAdmin::freeDOFIndex(int dof) {
FUNCNAME("DOFAdmin::freeDOFIndex()");
......@@ -184,29 +184,42 @@ namespace AMDiS {
void DOFAdmin::addDOFIndexed(DOFIndexedBase* dofIndexed) {
FUNCNAME("DOFAdmin::addDOFIndexed()");
TEST_EXIT_DBG(dofIndexed)("no dofIndexed\n");
if (dofIndexed->getSize() < size) {
dofIndexed->resize(size);
}
TEST_EXIT(dofIndexed)("no dofIndexed\n");
dofIndexedList.push_back(dofIndexed);
#ifdef _OPENMP
#pragma omp critical (dofIndexAccess)
#endif
{
if (dofIndexed->getSize() < size) {
dofIndexed->resize(size);
}
dofIndexedList.push_back(dofIndexed);
}
}
void DOFAdmin::removeDOFIndexed(DOFIndexedBase* dofIndexed)
{
FUNCNAME("DOFAdmin::removeDOFIndexed()");
std::list<DOFIndexedBase*>::iterator it;
std::list<DOFIndexedBase*>::iterator end = dofIndexedList.end();
for (it = dofIndexedList.begin(); it != end; ++it) {
if (*it == dofIndexed) {
dofIndexedList.erase(it);
return;
bool removed = false;
#ifdef _OPENMP
#pragma omp critical (dofIndexAccess)
#endif
{
std::list<DOFIndexedBase*>::iterator it;
std::list<DOFIndexedBase*>::iterator end = dofIndexedList.end();
for (it = dofIndexedList.begin(); it != end; ++it) {
if (*it == dofIndexed) {
dofIndexedList.erase(it);
removed = true;
break;
}
}
}
ERROR("DOFIndexed not in list\n");
TEST_EXIT(removed)("DOFIndexed not in list\n");
}
void DOFAdmin::addDOFContainer(DOFContainer* cont)
......@@ -315,7 +328,8 @@ namespace AMDiS {
}
DOFAdmin::~DOFAdmin()
{}
{
}
void DOFAdmin::serialize(std::ostream &out)
{
......
AMDiS/src/DOFAdmin.h
View file @ 3dc9cf95
......@@ -36,6 +36,7 @@
#include "FixVec.h"
#include "MemoryManager.h"
#include "Serializable.h"
#include "OpenMP.h"
#include <vector>
#include <memory>
#include <list>
......
AMDiS/src/DOFMatrix.cc
View file @ 3dc9cf95
......@@ -441,10 +441,10 @@ namespace AMDiS {
if (!op && operators.size() == 0) {
return;
}
Operator *operat = op ? op : operators[0];
operat->getAssembler(omp_get_thread_num())->initElementMatrix(elementMatrix, elInfo);
if (op) {
op->getElementMatrix(elInfo, elementMatrix);
} else {
......
AMDiS/src/ElInfo.cc
View file @ 3dc9cf95
......@@ -14,7 +14,8 @@
namespace AMDiS {
int ElInfo::traverseId = -1;
int ElInfo::traverseId[16] = {-1, -1, -1, -1, -1, -1, -1, -1,
-1, -1, -1, -1, -1, -1, -1, -1};
int ElInfo::traverseIdCounter = 0;
ElInfo::ElInfo(Mesh *aMesh)
......
AMDiS/src/ElInfo.h
View file @ 3dc9cf95
......@@ -586,9 +586,10 @@ namespace AMDiS {
static const int childEdge[3][2][6];
/** \brief
* Used to determine to which traverse an ElInfo belongs
* Used to determine to which traverse an ElInfo belongs.
* Stores the value for each thread, maximum value 16!.
*/
static int traverseId;
static int traverseId[16];
/** \brief
* Used to determine to which traverse an ElInfo belongs
......
AMDiS/src/Error.h
View file @ 3dc9cf95
......@@ -83,9 +83,9 @@ namespace AMDiS {
static int relFct(ElInfo* elinfo);
public:
static const T& errUFct(const DimVec<double>& lambda);
static T errUFct(const DimVec<double>& lambda);
static const WorldVector<T>& grdErrUFct(const DimVec<double>& lambda);
static WorldVector<T> grdErrUFct(const DimVec<double>& lambda);
/** \brief
*
......@@ -95,12 +95,12 @@ namespace AMDiS {
public:
AbstrFctErrU() : AbstractFunction<T, DimVec<double> >(0) {};
inline const T& operator()(const DimVec<double> & x) const {
return Error<T>::errUFct(x);
inline T operator()(const DimVec<double> & x) const {
return Error<T>::errUFct(x);
};
};
static AbstrFctErrU errU;
static AbstrFctErrU errU;
/** \brief
*
......@@ -110,8 +110,8 @@ namespace AMDiS {
public:
AbstrFctGrdErrU() : AbstractFunction<WorldVector<T>, DimVec<double> >(0) {};
inline const WorldVector<T>& operator()(const DimVec<double> & x) const {
return Error<T> :: grdErrUFct (x);
inline WorldVector<T> operator()(const DimVec<double> & x) const {
return Error<T>::grdErrUFct(x);
};
};
......
AMDiS/src/Error.hh
View file @ 3dc9cf95
......@@ -5,26 +5,18 @@
namespace AMDiS {
template<typename T>
const T& Error<T>::errUFct(const DimVec<double>& lambda)
T Error<T>::errUFct(const DimVec<double>& lambda)
{
WorldVector<double> x;
// if (el_parametric) {
// ERROR_EXIT("not yet\n");
// } else {
elinfo->coordToWorld(lambda, &x);
// }
return((*pU)(x));
}
template<typename T>
const WorldVector<T>& Error<T>::grdErrUFct(const DimVec<double>& lambda)
WorldVector<T> Error<T>::grdErrUFct(const DimVec<double>& lambda)
{
WorldVector<double> x;
// if (el_parametric) {
// ERROR_EXIT("not yet\n");
// } else {
elinfo->coordToWorld(lambda, &x);
// }
return((*pGrdU)(x));
}
......@@ -32,23 +24,11 @@ namespace AMDiS {
template<typename T>
int Error<T>::maxErrAtQpFct(ElInfo* el_info)
{
int i;
double err;
const double *u_vec, *uh_vec;
// Parametric *parametric = el_info->getMesh()->getParametric();
elinfo = el_info;
// if (parametric && parametric->initElement(el_info)) {
// el_parametric = parametric;
// } else {
// el_parametric = NULL;
// }
u_vec = quadFast->getQuadrature()->fAtQp(errU, NULL);
//uh_el = errUh->getLocalVector(el_info->getElement(), NULL);
//uh_vec = quadFast->uhAtQp(uh_el, NULL);
uh_vec = errUh->getVecAtQPs(el_info,
NULL,
quadFast,
......@@ -56,11 +36,10 @@ namespace AMDiS {
int numPoints = quadFast->getNumPoints();
for (i = 0; i < numPoints; i++)
{
err = u_vec[i] > uh_vec[i] ? u_vec[i] - uh_vec[i] : uh_vec[i] - u_vec[i];
maxErr = max(maxErr, err);
}
for (int i = 0; i < numPoints; i++) {
err = u_vec[i] > uh_vec[i] ? u_vec[i] - uh_vec[i] : uh_vec[i] - u_vec[i];
maxErr = max(maxErr, err);
}
return;
}
......@@ -70,14 +49,14 @@ namespace AMDiS {
const DOFVector<T>& uh,
const Quadrature* q)
{
FUNCNAME("Error<T>::maxErrAtQp");
FUNCNAME("Error<T>::maxErrAtQp()");
const FiniteElemSpace *fe_space;
if (!(pU = &u))
{
ERROR("no function u specified; doing nothing\n");
return(-1.0);
}
if (!(pU = &u)) {
ERROR("no function u specified; doing nothing\n");
return(-1.0);
}
if (!(errUh = &uh) || !(fe_space = uh->getFESpace()))
{
ERROR("no discrete function or no fe_space for it; doing nothing\n");
......@@ -113,32 +92,11 @@ namespace AMDiS {
{
int i, j;
double err, err_2, h1_err_el, norm_el, norm2, det, exact;
//int dim = el_info->getMesh()->getDim();
//const DimVec<WorldVector<double> > &Lambda = el_info->getGrdLambda();
//const T *uh_el;
const WorldVector<double> *grdu_vec, *grduh_vec;
// const Parametric *parametric = el_info->getMesh()->getParametric();
elinfo = el_info;
// if (parametric && parametric->initElement(el_info)) {
// el_parametric = parametric;
// }
// else {
// el_parametric = NULL;
// }
grdu_vec = quadFast->getQuadrature()->grdFAtQp(grdErrU, NULL);
//uh_el = errUh->getLocalVector(el_info->getElement(), NULL);
// if (el_parametric) {
// el_parametric->grd_lambda(el_info, NULL, 1, NULL, &Lambda, &det);
// }
// else {
det = el_info->getDet();
// }
//grduh_vec = quadFast->grdUhAtQp(Lambda, uh_el, NULL);
grduh_vec = errUh->getGrdAtQPs(elinfo, NULL, quadFast, NULL);
int numPoints = quadFast->getNumPoints();
......@@ -177,7 +135,6 @@ namespace AMDiS {
double Error<T>::H1Err(
const AbstractFunction<WorldVector<T>, WorldVector<double> >& grdU,
const DOFVector<T>& uh,
//const Quadrature* quad,
int relErr,
double* max,
bool writeLeafData,
......@@ -212,7 +169,6 @@ namespace AMDiS {
int deg = grdU.getDegree();
int degree = deg ? deg : 2 * fe_space->getBasisFcts()->getDegree() - 2;
//if (!quad)
q = Quadrature::provideQuadrature(dim, degree);
quadFast = FastQuadrature::provideFastQuadrature(basFct,
*q,
......@@ -250,34 +206,13 @@ namespace AMDiS {
int i;
double err, det, l2_err_el, norm_el, exact;
const double *u_vec, *uh_vec;
// Parametric *parametric = const_cast<Parametric*>(el_info->getMesh()->getParametric());
elinfo = el_info;
// if (parametric && parametric->initElement(el_info)) {
// el_parametric = parametric;
// }
// else {
// el_parametric = NULL;
// }
u_vec = quadFast->getQuadrature()->fAtQp(errU, NULL);
//uh_el = errUh->getLocalVector(el_info->getElement(), NULL);
//uh_vec = quadFast->uhAtQp(uh_el, NULL);
uh_vec = errUh->getVecAtQPs(el_info,
NULL,
quadFast,
NULL);
// if (el_parametric) {
// el_parametric->det(el_info, NULL, 1, NULL, &det);
// }
// else {
// det = el_info->calcDet();
// }
det = el_info->getDet();
int numPoints = quadFast->getNumPoints();
......
AMDiS/src/FileWriter.cc
View file @ 3dc9cf95
......@@ -28,8 +28,6 @@ namespace AMDiS {
solutionVecs_.resize(1);
solutionVecs_[0] = vec;
dataCollectors_.resize(1);
}
......@@ -50,8 +48,6 @@ namespace AMDiS {
feSpace = vecs[0]->getFESpace();
solutionVecs_ = vecs;
dataCollectors_.resize(vecs.size());
}
......@@ -83,8 +79,6 @@ namespace AMDiS {
}
feSpace = vec->getFESpace();
dataCollectors_.resize(nTmpSolutions_);
}
......@@ -97,12 +91,6 @@ namespace AMDiS {
DELETE solutionVecs_[i];
}
}
for (int i = 0; i < static_cast<int>(dataCollectors_.size()); i++) {
if (dataCollectors_[i]) {
DELETE dataCollectors_[i];
}
}
}
......@@ -180,15 +168,18 @@ namespace AMDiS {
ERROR_EXIT("This should not happen!\n");
}
// Containers, which store the data to be written;
std::vector< DataCollector* > dataCollectors(solutionVecs_.size());
if (writeElem) {
for (int i = 0; i < static_cast<int>(dataCollectors_.size()); i++) {
dataCollectors_[i] = NEW DataCollector(feSpace, solutionVecs_[i],
level, flag, writeElem);
for (int i = 0; i < static_cast<int>(dataCollectors.size()); i++) {
dataCollectors[i] = NEW DataCollector(feSpace, solutionVecs_[i],
level, flag, writeElem);
}
} else {
for (int i = 0; i < static_cast<int>(dataCollectors_.size()); i++) {
dataCollectors_[i] = NEW DataCollector(feSpace, solutionVecs_[i],
traverseLevel, flag | traverseFlag, writeElement);
for (int i = 0; i < static_cast<int>(dataCollectors.size()); i++) {
dataCollectors[i] = NEW DataCollector(feSpace, solutionVecs_[i],
traverseLevel, flag | traverseFlag, writeElement);
}
}
......@@ -214,8 +205,8 @@ namespace AMDiS {
ERROR_EXIT("Delay writing only supported for ParaView file format!\n");
}
for (int i = 0; i < static_cast<int>(dataCollectors_.size()); i++) {
dataCollectors_[i]->fillAllData();
for (int i = 0; i < static_cast<int>(dataCollectors.size()); i++) {
dataCollectors[i]->fillAllData();
}
writingIsDelayed_ = true;
......@@ -236,26 +227,26 @@ namespace AMDiS {
if (writeAMDiSFormat) {
TEST_EXIT(mesh)("no mesh\n");
MacroWriter::writeMacro(dataCollectors_[0],
MacroWriter::writeMacro(dataCollectors[0],
const_cast<char*>((fn + amdisMeshExt).c_str()),
adaptInfo ? adaptInfo->getTime() : 0.0);
MSG("macro file written to %s\n", (fn + amdisMeshExt).c_str());
ValueWriter::writeValues(dataCollectors_[0],
ValueWriter::writeValues(dataCollectors[0],
(fn + amdisDataExt).c_str(),
adaptInfo ? adaptInfo->getTime() : 0.0);
MSG("value file written to %s\n", (fn + amdisDataExt).c_str());
}
if (writePeriodicFormat) {
MacroWriter::writePeriodicFile(dataCollectors_[0],
MacroWriter::writePeriodicFile(dataCollectors[0],
(fn + periodicFileExt).c_str());
MSG("periodic file written to %s\n", (fn + periodicFileExt).c_str());
}
if (writeParaViewFormat) {
VtkWriter vtkWriter(&dataCollectors_);
VtkWriter vtkWriter(&dataCollectors);
vtkWriter.setCompression(compression);
vtkWriter.writeFile(const_cast<char*>((fn + paraViewFileExt).c_str()));
......@@ -263,41 +254,20 @@ namespace AMDiS {
}
if (writeParaViewAnimation) {
VtkWriter vtkWriter(&dataCollectors_);
VtkWriter vtkWriter(&dataCollectors);
vtkWriter.updateAnimationFile(fn + paraViewFileExt,
&paraViewAnimationFrames_,
const_cast<char*>((filename + ".pvd").c_str()));
}
for (int i = 0; i < static_cast<int>(dataCollectors_.size()); i++) {
DELETE dataCollectors_[i];
for (int i = 0; i < static_cast<int>(dataCollectors.size()); i++) {
DELETE dataCollectors[i];
}
}
void FileWriter::writeDelayedFiles()
{
if (!writingIsDelayed_) {
return;
}
if (writeParaViewFormat) {
VtkWriter vtkWriter(&dataCollectors_);
vtkWriter.writeFile(const_cast<char*>((delayedFilename_ + paraViewFileExt).c_str()));
}
if (writeParaViewAnimation) {
VtkWriter vtkWriter(&dataCollectors_);
vtkWriter.updateAnimationFile(delayedFilename_ + paraViewFileExt,
&paraViewAnimationFrames_,
const_cast<char*>((filename + ".pvd").c_str()));
}
for (int i = 0; i < static_cast<int>(dataCollectors_.size()); i++) {
DELETE dataCollectors_[i];
}
writingIsDelayed_ = false;
delayedFilename_ = "";
ERROR_EXIT("no more!\n");