Skip to content

GitLab

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

Fixed a bug when using openmp parallelization for assemlbing matrices.

parent da73b714
Hide whitespace changes
Inline Side-by-side
AMDiS/src/BasisFunction.cc
View file @ a6ae4ebc
......@@ -140,6 +140,7 @@ namespace AMDiS {
return ((*val));
}
int BasisFunction::getNumberOfDOFs(int i) const
{
return (*nDOF)[i];
......
AMDiS/src/DOFMatrix.cc
View file @ a6ae4ebc
......@@ -18,6 +18,7 @@
namespace AMDiS {
using namespace mtl;
using boost::lexical_cast;
DOFMatrix *DOFMatrix::traversePtr = NULL;
......@@ -186,7 +187,6 @@ namespace AMDiS {
}
}
for (int i = 0; i < nRow; i++) {
DegreeOfFreedom row = rowIndices[i];
......
AMDiS/src/DOFVector.cc
View file @ a6ae4ebc
......@@ -71,11 +71,8 @@ namespace AMDiS {
}
Mesh *mesh = feSpace->getMesh();
int dim = mesh->getDim();
const BasisFunction *basFcts = feSpace->getBasisFcts();
DOFAdmin *admin = feSpace->getAdmin();
// count number of nodes and dofs per node
......@@ -99,30 +96,24 @@ namespace AMDiS {
numNodes += numPositionNodes;
}
// TEST_EXIT_DBG(numNodes == mesh->getNumberOfNodes())
// ("invalid number of nodes\n");
TEST_EXIT_DBG(numDOFs == basFcts->getNumber())
("number of dofs != number of basis functions\n");
for (int i = 0; i < numDOFs; i++) {
bary.push_back(basFcts->getCoords(i));
}
for (int i = 0; i < numDOFs; i++)
bary.push_back(basFcts->getCoords(i));
double *localUh = new double[basFcts->getNumber()];
// traverse mesh
std::vector<bool> visited(getUsedSize(), false);
TraverseStack stack;
Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
while (elInfo) {
const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
const double *localUh = getLocalVector(elInfo->getElement(), NULL);
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
getLocalVector(elInfo->getElement(), localUh);
int localDOFNr = 0;
for (int i = 0; i < numNodes; i++) { // for all nodes
......@@ -141,6 +132,8 @@ namespace AMDiS {
elInfo = stack.traverseNext(elInfo);
}
delete [] localUh;
return result;
}
......@@ -150,6 +143,10 @@ namespace AMDiS {
{
FUNCNAME("DOFVector<double>::getRecoveryGradient()");
#ifdef _OPENMP
ERROR_EXIT("Using static variable while using OpenMP parallelization!\n");
#endif
// define result vector
static DOFVector<WorldVector<double> > *vec = NULL;
......@@ -160,9 +157,9 @@ namespace AMDiS {
delete vec;
vec = NULL;
}
if (!vec) {
if (!vec)
vec = new DOFVector<WorldVector<double> >(feSpace, "gradient");
}
result = vec;
}
......@@ -183,11 +180,13 @@ namespace AMDiS {
Mesh::CALL_LEAF_EL | Mesh::FILL_DET |
Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS);
double *localUh = new double[basFcts->getNumber()];
while (elInfo) {
double det = elInfo->getDet();
const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
const double *localUh = getLocalVector(elInfo->getElement(), NULL);
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
getLocalVector(elInfo->getElement(), localUh);
basFcts->evalGrdUh(bary, grdLambda, localUh, &grd);
for (int i = 0; i < dim + 1; i++) {
......@@ -199,14 +198,14 @@ namespace AMDiS {
elInfo = stack.traverseNext(elInfo);
}
delete [] localUh;
DOFVector<double>::Iterator volIt(&volume, USED_DOFS);
DOFVector<WorldVector<double> >::Iterator grdIt(result, USED_DOFS);
for (volIt.reset(), grdIt.reset(); !volIt.end(); ++volIt, ++grdIt) {
if (*volIt != 0.0) {
for (volIt.reset(), grdIt.reset(); !volIt.end(); ++volIt, ++grdIt)
if (*volIt != 0.0)
*grdIt *= 1.0 / (*volIt);
}
}
return result;
}
......@@ -242,9 +241,9 @@ namespace AMDiS {
delete [] grd;
grd = new WorldVector<double>[nPoints];
for (int i = 0; i < nPoints; i++) {
for (int i = 0; i < nPoints; i++)
grd[i].set(0.0);
}
result = grd;
}
......@@ -612,6 +611,10 @@ namespace AMDiS {
{
FUNCNAME("DOFVector<double>::getGradient()");
#ifdef _OPENMP
ERROR_EXIT("Using static variable while using OpenMP parallelization!\n");
#endif
Mesh *mesh = feSpace->getMesh();
int dim = mesh->getDim();
int dow = Global::getGeo(WORLD);
......@@ -660,15 +663,11 @@ namespace AMDiS {
numNodes += numPositionNodes;
}
// TEST_EXIT_DBG(numNodes == mesh->getNumberOfNodes())
// ("invalid number of nodes\n");
TEST_EXIT_DBG(numDOFs == basFcts->getNumber())
("number of dofs != number of basis functions\n");
for (int i = 0; i < numDOFs; i++) {
bary.push_back(basFcts->getCoords(i));
}
for (int i = 0; i < numDOFs; i++)
bary.push_back(basFcts->getCoords(i));
// traverse mesh
std::vector<bool> visited(getUsedSize(), false);
......@@ -691,9 +690,8 @@ namespace AMDiS {
if (!visited[dofIndex]) {
basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda, localUh, &grd);
for (int k = 0; k < dow; k++) {
(*(*result)[k])[dofIndex] = grd[k];
}
for (int k = 0; k < dow; k++)
(*(*result)[k])[dofIndex] = grd[k];
visited[dofIndex] = true;
}
......@@ -779,9 +777,8 @@ namespace AMDiS {
if (localvec)
delete [] localvec;
localvec = new double[nPoints];
for (int i = 0; i < nPoints; i++) {
localvec[i] = 0.0;
}
for (int i = 0; i < nPoints; i++)
localvec[i] = 0.0;
result = localvec;
}
......@@ -795,9 +792,9 @@ namespace AMDiS {
result[i] = 0.0;
for (int j = 0; j < nBasFcts; j++) {
double val = 0.0;
for (int k = 0; k < nBasFcts; k++) {
for (int k = 0; k < nBasFcts; k++)
val += m[j][k] * (*(basFcts->getPhi(k)))(quad->getLambda(i));
}
result[i] += localVec[j] * val;
}
}
......
AMDiS/src/DOFVector.hh
View file @ a6ae4ebc
......@@ -962,6 +962,9 @@ namespace AMDiS {
if (d) {
result = d;
} else {
#ifdef _OPENMP
ERROR_EXIT("Using static variable while using OpenMP parallelization!\n");
#endif
if (localVec && nBasFcts > localVecSize) {
delete [] localVec;
localVec = new T[nBasFcts];
......@@ -1002,35 +1005,40 @@ namespace AMDiS {
Element *el = elInfo->getElement();
const Quadrature *quadrature = quadFast ? quadFast->getQuadrature() : quad;
const BasisFunction *basFcts = feSpace->getBasisFcts();
int numPoints = quadrature->getNumPoints();
int nPoints = quadrature->getNumPoints();
int nBasFcts = basFcts->getNumber();
int j;
static T *localvec = NULL;
T *result;
if (vecAtQPs) {
result = vecAtQPs;
} else {
#ifdef _OPENMP
ERROR_EXIT("Using static variable while using OpenMP parallelization!\n");
#endif
if (localvec)
delete [] localvec;
localvec = new T[numPoints];
for (int i = 0; i < numPoints; i++)
localvec = new T[nPoints];
for (int i = 0; i < nPoints; i++)
localvec[i] = 0.0;
result = localvec;
}
const T *localVec = getLocalVector(el, NULL);
T *localVec = new T[nBasFcts];
getLocalVector(el, localVec);
for (int i = 0; i < numPoints; i++) {
for (result[i] = j = 0; j < nBasFcts; j++) {
for (int i = 0; i < nPoints; i++) {
result[i] = 0.0;
for (int j = 0; j < nBasFcts; j++)
result[i] += localVec[j] *
(quadFast ?
(quadFast->getPhi(i, j)) :
((*(basFcts->getPhi(j)))(quad->getLambda(i))));
}
((*(basFcts->getPhi(j)))(quad->getLambda(i))));
}
delete [] localVec;
return const_cast<const T*>(result);
}
......
AMDiS/src/ElementFunction.h
View file @ a6ae4ebc
......@@ -83,24 +83,25 @@ namespace AMDiS {
{
public:
/// constructor.
ElementFunctionDOFVec(const DOFVector<T> *dofVector)
ElementFunctionDOFVec(const DOFVector<T> *vec)
: ElementFunction<T>(),
dofVector_(dofVector)
dofVector(vec)
{}
/// evaluation at given coordinates.
T operator()(const DimVec<double>& bary) const
{
static T t;
const T* localVec =
dofVector_->getLocalVector(this->elInfo_->getElement(), NULL);
t = dofVector_->getFESpace()->getBasisFcts()->evalUh(bary, localVec);
T* localVec = new T[dofVector->getFESpace()->getBasisFcts()->getNumber()];
dofVector->getLocalVector(this->elInfo_->getElement(), localVec);
T t = dofVector->getFESpace()->getBasisFcts()->evalUh(bary, localVec);
delete [] localVec;
return t;
}
protected:
/// DOFVector to be interpolated.
const DOFVector<T> *dofVector_;
const DOFVector<T> *dofVector;
};
}
......
AMDiS/src/Lagrange.cc
View file @ a6ae4ebc
......@@ -2774,11 +2774,11 @@ namespace AMDiS {
DegreeOfFreedom pdof[4], dof9;
const DOFAdmin *admin;
if (n < 1) return;
el = list->getElement(0);
if (n < 1)
return;
el = list->getElement(0);
admin = drv->getFESpace()->getAdmin();
basFct->getLocalIndices(el, admin, pdof);
/****************************************************************************/
......@@ -2824,7 +2824,8 @@ namespace AMDiS {
(*drv)[cdof[5]] + 0.9375*(*drv)[cdof[6]] + 0.1875*(*drv)[cdof[9]];
(*drv)[pdof[9]] += 0.75*(*drv)[cdof[9]];
if (n <= 1) return;
if (n <= 1)
return;
/****************************************************************************/
/* adjust the value on the neihgbour */
......
AMDiS/src/Operator.cc
View file @ a6ae4ebc
......@@ -6,8 +6,13 @@
#include "Quadrature.h"
#include "OpenMP.h"
#include "boost/lexical_cast.hpp"
#include <fstream>
namespace AMDiS {
using boost::lexical_cast;
const Flag OperatorTerm::PW_CONST = 1;
const Flag OperatorTerm::SYMMETRIC = 2;
......@@ -71,7 +76,7 @@ namespace AMDiS {
{
FUNCNAME("OperatorTerm::getVectorAtQPs()");
TEST_EXIT(elInfo->getMesh() == vec->getFESpace()->getMesh())
TEST_EXIT_DBG(elInfo->getMesh() == vec->getFESpace()->getMesh())
("There is something wrong!\n");
return subAssembler->getVectorAtQPs(vec, elInfo, quad);
......@@ -94,23 +99,20 @@ namespace AMDiS {
// Both elements have the same size, so we can use the simple procedure
// to determine the vecAtQPs.
if (vec->getFESpace()->getMesh() == smallElInfo->getMesh()) {
if (vec->getFESpace()->getMesh() == smallElInfo->getMesh())
return subAssembler->getVectorAtQPs(vec, smallElInfo, quad);
} else {
return subAssembler->getVectorAtQPs(vec, largeElInfo, quad);
}
else
return subAssembler->getVectorAtQPs(vec, largeElInfo, quad);
} else {
// The two elements are different. If the vector is defined on the mesh of the
// small element, we can still use the simple procedure to determine the vecAtQPs.
if (vec->getFESpace()->getMesh() == largeElInfo->getMesh()) {
if (vec->getFESpace()->getMesh() == largeElInfo->getMesh())
return subAssembler->getVectorAtQPs(vec, smallElInfo, largeElInfo, quad);
} else {
else
return subAssembler->getVectorAtQPs(vec, smallElInfo, quad);
}
}
}
......@@ -165,7 +167,7 @@ namespace AMDiS {
const WorldMatrix<double>& matrix,
DimMat<double>& LALt,
bool symm,
double factor)
double factor) const
{
int j, k, l;
const int dimOfWorld = Global::getGeo(WORLD);
......@@ -219,7 +221,7 @@ namespace AMDiS {
void OperatorTerm::l1lt(const DimVec<WorldVector<double> >& Lambda,
DimMat<double>& LALt,
double factor)
double factor) const
{
const int dimOfWorld = Global::getGeo(WORLD);
int dim = LALt.getNumRows() - 1;
......@@ -342,15 +344,14 @@ namespace AMDiS {
#pragma omp critical (initAssembler)
#endif
{
if (optimized) {
assembler[rank] = new OptimizedAssembler(this,
quad2, quad1GrdPsi, quad1GrdPhi, quad0,
rowFESpace, colFESpace);
} else {
assembler[rank] = new StandardAssembler(this,
quad2, quad1GrdPsi, quad1GrdPhi, quad0,
rowFESpace, colFESpace);
}
if (optimized)
assembler[rank] =
new OptimizedAssembler(this, quad2, quad1GrdPsi, quad1GrdPhi, quad0,
rowFESpace, colFESpace);
else
assembler[rank] =
new StandardAssembler(this, quad2, quad1GrdPsi, quad1GrdPhi, quad0,
rowFESpace, colFESpace);
}
}
......
AMDiS/src/Operator.h
View file @ a6ae4ebc
......@@ -144,11 +144,11 @@ namespace AMDiS {
protected:
/// Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$.
static void lalt(const DimVec<WorldVector<double> >& Lambda,
const WorldMatrix<double>& matrix,
DimMat<double>& LALt,
bool symm,
double factor);
void lalt(const DimVec<WorldVector<double> >& Lambda,
const WorldMatrix<double>& matrix,
DimMat<double>& LALt,
bool symm,
double factor) const;
/** \brief
* Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$ for \f$ A \f$
......@@ -161,9 +161,9 @@ namespace AMDiS {
double factor);
/// Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$ for A equal to the identity.
static void l1lt(const DimVec<WorldVector<double> >& Lambda,
DimMat<double>& LALt,
double factor);
void l1lt(const DimVec<WorldVector<double> >& Lambda,
DimMat<double>& LALt,
double factor) const;
/// Evaluation of \f$ \Lambda \cdot b\f$.
static inline void lb(const DimVec<WorldVector<double> >& Lambda,
......@@ -628,20 +628,14 @@ namespace AMDiS {
Vec2AtQP_SOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
BinaryAbstractFunction<double, double, double> *af);
/** \brief
* Implementation of \ref OperatorTerm::initElement().
*/
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
/** \brief
* Implements SecondOrderTerm::getLALt().
*/
/// Implements SecondOrderTerm::getLALt().
void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
/** \brief
* Implenetation of SecondOrderTerm::eval().
*/
/// Implenetation of SecondOrderTerm::eval().
void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
......@@ -649,30 +643,21 @@ namespace AMDiS {
double *result,
double factor) const;
/** \brief
* Implenetation of SecondOrderTerm::weakEval().
*/
/// Implenetation of SecondOrderTerm::weakEval().
void weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const;
protected:
/** \brief
* DOFVector to be evaluated at quadrature points.
*/
/// DOFVector to be evaluated at quadrature points.
DOFVectorBase<double>* vec1;
DOFVectorBase<double>* vec2;
/** \brief
* Pointer to an array containing the DOFVector evaluated at quadrature
* points.
*/
/// Pointer to an array containing the DOFVector evaluated at quadrature points.
double* vecAtQPs1;
double* vecAtQPs2;
/** \brief
* Function evaluated at \ref vecAtQPs.
*/
/// Function evaluated at \ref vecAtQPs.
BinaryAbstractFunction<double, double, double> *f;
};
......
AMDiS/src/ProblemVec.cc
View file @ a6ae4ebc
......@@ -770,7 +770,7 @@ namespace AMDiS {
matrix->getBoundaryManager()->exitMatrix(matrix);
if (matrix)
nnz += matrix->getBaseMatrix().nnz();
nnz += matrix->getBaseMatrix().nnz();
}
// And now assemble boundary conditions on the vectors
......
AMDiS/src/Quadrature.cc
View file @ a6ae4ebc
......@@ -39,7 +39,7 @@ namespace AMDiS {
{
FUNCNAME("Quadrature::grdFAtQp()");
static WorldVector<double> *quad_vec_d = NULL;
static WorldVector<double> *quad_vec_d = NULL;
static size_t size = 0;
WorldVector<double> *val;
WorldVector<double> grd;
......@@ -47,7 +47,11 @@ namespace AMDiS {
if (vec) {
val = vec;
} else {
if (static_cast<int>( size) < n_points) {
#ifdef _OPENMP
ERROR_EXIT("Using static variable while using OpenMP parallelization!\n");
#endif
if (static_cast<int>(size) < n_points) {
size_t new_size = std::max(maxNQuadPoints[dim], n_points);
delete [] quad_vec_d;
quad_vec_d = new WorldVector<double>[new_size];
......@@ -67,9 +71,8 @@ namespace AMDiS {
return val;
}
const double
*Quadrature::fAtQp(const AbstractFunction<double, DimVec<double> >& f,
double *vec) const
const double *Quadrature::fAtQp(const AbstractFunction<double, DimVec<double> >& f,
double *vec) const
{
FUNCNAME("Quadrature::fAtQp()");
......@@ -80,6 +83,10 @@ namespace AMDiS {
if (vec) {
val = vec;
} else {
#ifdef _OPENMP
ERROR_EXIT("Using static variable while using OpenMP parallelization!\n");
#endif
if (static_cast<int>( size) < n_points) {
size_t new_size = std::max(maxNQuadPoints[dim], n_points);
delete [] quad_vec;
......@@ -1532,9 +1539,8 @@ namespace AMDiS {
// fill memory
for (int i = 0; i< nPoints; i++) {
lambda = quadrature->getLambda(i);
for (int j = 0; j < nBasFcts; j++) {
for (int j = 0; j < nBasFcts; j++)
(*(basisFunctions->getGrdPhi(j)))(lambda, (*(grdPhi))[i][j]);
}
}
// update flag
......@@ -1551,9 +1557,8 @@ namespace AMDiS {
// fill memory
for (int i = 0; i < nPoints; i++) {
lambda = quadrature->getLambda(i);
for (int j = 0; j < nBasFcts; j++) {
for (int j = 0; j < nBasFcts; j++)
(*(basisFunctions->getD2Phi(j)))(lambda, (*(D2Phi))[i][j]);
}
}
// update flag
......
AMDiS/src/SecondOrderAssembler.cc
View file @ a6ae4ebc
......@@ -186,11 +186,16 @@ namespace AMDiS {
tmpLALt[myRank] = new DimMat<double>*[nPoints];
for (int j = 0; j < nPoints; j++)
tmpLALt[myRank][j] = new DimMat<double>(dim, NO_INIT);
const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
basFcts = owner->getColFESpace()->getBasisFcts();
phiFast = updateFastQuadrature(phiFast, basFcts, INIT_GRD_PHI);
#ifdef _OPENMP
#pragma omp critical (dofIndexAccess)
#endif
{
psiFast = updateFastQuadrature(psiFast, owner->getRowFESpace()->getBasisFcts(),
INIT_GRD_PHI);
phiFast = updateFastQuadrature(phiFast, owner->getRowFESpace()->getBasisFcts(),
INIT_GRD_PHI);
}
firstCall = false;
}
......@@ -211,7 +216,7 @@ namespace AMDiS {
for (int iq = 0; iq < nPoints; iq++) {
(*LALt[iq]) *= elInfo->getDet();
grdPsi = psiFast->getGradient(iq);
grdPhi = phiFast->getGradient(iq);
......@@ -267,9 +272,8 @@ namespace AMDiS {
int myRank = omp_get_thread_num();