From f8e50c2d6e60d1d4079e9bda509de5b1bed9ea7e Mon Sep 17 00:00:00 2001 From: akva Date: Wed, 18 Jan 2012 13:43:33 +0000 Subject: [PATCH] changed: reorder operations in GlbNorm to make them work with multithreaded assembly git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1405 e10b68d5-8a6e-419e-a041-bce267b0401d --- src/ASM/ElmMats.h | 3 --- src/ASM/ElmNorm.h | 41 +++++++++++++++++++++++++++++++++------- src/ASM/GlbNorm.C | 7 +++++++ src/ASM/GlbNorm.h | 2 +- src/ASM/GlobalIntegral.h | 2 +- src/ASM/IntegrandBase.h | 3 +++ src/SIM/SIMbase.C | 36 ++++++++++++++++++++++------------- 7 files changed, 69 insertions(+), 25 deletions(-) diff --git a/src/ASM/ElmMats.h b/src/ASM/ElmMats.h index f5d86f62..84d88fbd 100644 --- a/src/ASM/ElmMats.h +++ b/src/ASM/ElmMats.h @@ -70,9 +70,6 @@ public: return b.front(); } - //! \brief Virtual destruction method to clean up after numerical integration. - virtual void destruct() { delete this; } - std::vector A; //!< The element coefficient matrices std::vector b; //!< The element right-hand-side vectors diff --git a/src/ASM/ElmNorm.h b/src/ASM/ElmNorm.h index 4b9d09fa..38d7e99e 100644 --- a/src/ASM/ElmNorm.h +++ b/src/ASM/ElmNorm.h @@ -14,8 +14,8 @@ #ifndef _ELM_NORM_H #define _ELM_NORM_H -#include #include "LocalIntegral.h" +#include /*! @@ -29,7 +29,17 @@ class ElmNorm : public LocalIntegral { public: //! \brief The constructor assigns the internal pointer. + //! \param[in] p Pointer to element norm values + //! \param[in] n Number of norm values ElmNorm(double* p, size_t n) : ptr(p), nnv(n) {} + //! \brief Alternative constructor using the internal buffer \a buf. + //! \param[in] n Number of norm values + //! + //! \details This constructor is used when the element norms are not requested + //! by the application, but are only used to assembly the global norms. + //! To avoid the need for a global array of element norms in that case, + //! an internal array is then used instead. + ElmNorm(size_t n) : buf(n,0.0), nnv(n) { ptr = &buf.front(); } //! \brief Empty destructor. virtual ~ElmNorm() {} @@ -41,15 +51,32 @@ public: //! \brief Returns the number of norm values. size_t size() const { return nnv; } + //! \brief Returns whether the element norms are stored externally or not. + bool externalStorage() const { return buf.empty(); } + //! \brief Virtual destruction method to clean up after numerical integration. - //! \details For this class we only clear the two solution vector containers - //! (to save memory). We do NOT delete the object itself, since it might be - //! needed in a second integration loop over element boundaries, for instance. - virtual void destruct() { vec.clear(); psol.clear(); } + //! \details Unless the internal buffer \a buf is used, these method only + //! clears the two solution vector containers (to save memory). + //! The object itself is then NOT deleted, since it might be used in a second + //! integration loop over element boundaries, for instance. + virtual void destruct() + { + if (buf.empty()) + { + vec.clear(); + psol.clear(); + } + else + delete this; // The internal buffer has been used, delete the whole thing + } private: - double* ptr; //!< Pointer to the actual norm values - size_t nnv; //!< Number of norm values + RealArray buf; //!< Internal buffer used when element norms are not requested + double* ptr; //!< Pointer to the actual norm values + size_t nnv; //!< Number of norm values + +public: + Vectors psol; //!< Element-level projected solution vectors }; #endif diff --git a/src/ASM/GlbNorm.C b/src/ASM/GlbNorm.C index 7ebf83f9..2fce3352 100644 --- a/src/ASM/GlbNorm.C +++ b/src/ASM/GlbNorm.C @@ -34,6 +34,13 @@ bool GlbNorm::assemble (const LocalIntegral* elmObj, int elmId) const ElmNorm* ptr = dynamic_cast(elmObj); if (!ptr) return false; + // If the element norms are requested (i.e. the internal buffer + // is not used) the actuall summation of element norms into the + // global norm is postponed to the very end (when invoked with + // elmId=0). This must be done like this to allow more than one + // loop over the elements during the norm integration. + if (elmId > 0 && ptr->externalStorage()) return true; + ElmNorm& elVals = *const_cast(ptr); for (size_t i = 0; i < elVals.size(); i++) { diff --git a/src/ASM/GlbNorm.h b/src/ASM/GlbNorm.h index 64a94ba4..476b3167 100644 --- a/src/ASM/GlbNorm.h +++ b/src/ASM/GlbNorm.h @@ -40,7 +40,7 @@ public: //! \brief Adds element norm quantities into the global norm object. //! \param[in] elmObj Pointer to the element norms to add into \a *this //! \param[in] elmId Global number of the element associated with \a *elmObj - virtual bool assemble(const LocalIntegral* elmObj, int elmId); + virtual bool assemble(const LocalIntegral* elmObj, int elmId = 0); private: //! \brief Applies the operation \a myOp on the given \a value. diff --git a/src/ASM/GlobalIntegral.h b/src/ASM/GlobalIntegral.h index 44304345..9e137197 100644 --- a/src/ASM/GlobalIntegral.h +++ b/src/ASM/GlobalIntegral.h @@ -34,7 +34,7 @@ public: //! \brief Adds a LocalIntegral object into a corresponding global object. //! \param[in] elmObj The local integral object to add into \a *this. //! \param[in] elmId Global number of the element associated with elmObj - virtual bool assemble(const LocalIntegral* elmObj, int elmId) = 0; + virtual bool assemble(const LocalIntegral* elmObj, int elmId) { return true; } }; #endif diff --git a/src/ASM/IntegrandBase.h b/src/ASM/IntegrandBase.h index 1d09c67e..92ee9746 100644 --- a/src/ASM/IntegrandBase.h +++ b/src/ASM/IntegrandBase.h @@ -22,6 +22,8 @@ class NormBase; class AnaSol; class VTF; +typedef std::vector LintegralVec; + /*! \brief Base class representing a system level integrated quantity. */ @@ -259,6 +261,7 @@ protected: Vectors prjsol; //!< Projected secondary solution vectors for current patch unsigned short int nrcmp; //!< Number of projected solution components + LintegralVec* lints; //!< Vector of local integrals used during norm integration }; #endif diff --git a/src/SIM/SIMbase.C b/src/SIM/SIMbase.C index 50106411..0587452d 100644 --- a/src/SIM/SIMbase.C +++ b/src/SIM/SIMbase.C @@ -1139,6 +1139,8 @@ bool SIMbase::solutionNorms (const TimeDomain& time, const Vectors& psol, const Vectors& ssol, Vector& gNorm, Matrix* eNorm) { + PROFILE1("Norm integration"); + NormBase* norm = myProblem->getNormIntegrand(mySol); if (!norm) { @@ -1148,8 +1150,6 @@ bool SIMbase::solutionNorms (const TimeDomain& time, return false; } - PROFILE1("Norm integration"); - myProblem->initIntegration(time); const Vector& primsol = psol.front(); @@ -1162,6 +1162,14 @@ bool SIMbase::solutionNorms (const TimeDomain& time, size_t nNorms = norm->getNoFields(); gNorm.resize(nNorms,true); +#ifdef USE_OPENMP + // When assembling in parallel, we must always do the norm summation + // at the end in a serial loop, to avoid that the threads try to update the + // same memory address simultaneously. + Matrix dummy; + if (!eNorm) eNorm = &dummy; +#endif + // Initialize norm integral classes GlbNorm globalNorm(gNorm,GlbNorm::SQRT); LintegralVec elementNorms; @@ -1171,7 +1179,9 @@ bool SIMbase::solutionNorms (const TimeDomain& time, elementNorms.reserve(eNorm->cols()); for (i = 0; i < eNorm->cols(); i++) elementNorms.push_back(new ElmNorm(eNorm->ptr(i),nNorms)); + norm->setLocalIntegrals(&elementNorms); } + norm->setLocalIntegrals(&elementNorms); // Loop over the different material regions, integrating solution norm terms // for the patch domain associated with each material @@ -1206,12 +1216,6 @@ bool SIMbase::solutionNorms (const TimeDomain& time, lp = i+1; } - // Integrate norm contributions due to Neumann boundary conditions, if any. - // Note: Currently, only the global norms are considered here. - // The corresponding element-level norms are not stored. This is mainly - // because the current design only allows one loop over the elements - // in the element-norm calculations. Consider rivising this later. - if (norm->hasBoundaryTerms()) for (i = 0; i < myProps.size() && ok; i++) if (myProps[i].pcode == Property::NEUMANN) @@ -1242,8 +1246,18 @@ bool SIMbase::solutionNorms (const TimeDomain& time, else ok = false; + // Clean up the dynamically allocated norm objects. This will also perform + // the actual global norm assembly, in case the element norms are stored, + // and always when multi-threading is used. + for (i = 0; i < elementNorms.size(); i++) + { + globalNorm.assemble(elementNorms[i]); + delete elementNorms[i]; + } + delete norm; + // Add problem-dependent external norm contributions - if (gNorm.size() >= 2) + if (gNorm.size() > 1) gNorm(2) += this->externalEnergy(psol); #ifdef PARALLEL_PETSC @@ -1256,10 +1270,6 @@ bool SIMbase::solutionNorms (const TimeDomain& time, } #endif - delete norm; - for (i = 0; i < elementNorms.size(); i++) - delete elementNorms[i]; - return ok; }