From da7e203f5263dc74114cb4993cf3df100a5e53ca Mon Sep 17 00:00:00 2001
From: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date: Tue, 2 Feb 2010 18:24:14 +0000
Subject: [PATCH] dg: limiter can use more than one group and works in parallel
---
Solver/TESTCASES/supra.lua | 4 +-
Solver/dgAlgorithm.cpp | 4 +-
Solver/dgDofContainer.cpp | 8 ++
Solver/dgDofContainer.h | 4 +
Solver/dgSlopeLimiter.cpp | 149 ++++++++++++++++++++-----------------
5 files changed, 95 insertions(+), 74 deletions(-)
diff --git a/Solver/TESTCASES/supra.lua b/Solver/TESTCASES/supra.lua
index 5b6a7e476c..7ad6ccc051 100644
--- a/Solver/TESTCASES/supra.lua
+++ b/Solver/TESTCASES/supra.lua
@@ -18,13 +18,13 @@ end
function valueRight(f)
for i=0,f:size1()-1 do
- f:set(i,0,math.sin(2*math.pi*t/T))
+ f:set(i,0,math.sin(2*math.pi*t/T+math.pi/2))
end
end
function valueLeft(f)
for i=0,f:size1()-1 do
- f:set(i,0,-math.sin(2*math.pi*t/T))
+ f:set(i,0,-math.sin(2*math.pi*t/T+math.pi/2))
end
end
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index 4e4c426205..db458c85e6 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -293,9 +293,7 @@ void dgAlgorithm::rungeKutta (const dgConservationLaw &claw, // conservation l
K.axpy(sol);
}
- if (limiter){
- limiter->apply(K, groups);
- }
+ if (limiter) limiter->apply(K, groups);
this->residual(claw,groups,K,resd);
K.scale(0.);
for(int k=0; k < groups.getNbElementGroups(); k++) {
diff --git a/Solver/dgDofContainer.cpp b/Solver/dgDofContainer.cpp
index 34aa06db1c..fe13a95aad 100644
--- a/Solver/dgDofContainer.cpp
+++ b/Solver/dgDofContainer.cpp
@@ -28,6 +28,7 @@ dgDofContainer::dgDofContainer (dgGroupCollection &groups, int nbFields):
_dataProxys.resize(_groups.getNbElementGroups()+_groups.getNbGhostGroups());
for (int i=0;i<_groups.getNbElementGroups();i++){
dgGroupOfElements *group = _groups.getElementGroup(i);
+ _groupId[group] = i;
int nbNodes = group->getNbNodes();
int nbElements = group->getNbElements();
_dataProxys[i] = new fullMatrix<double> (&(*_data)(offset),nbNodes, _nbFields*nbElements);
@@ -51,6 +52,7 @@ dgDofContainer::dgDofContainer (dgGroupCollection &groups, int nbFields):
dgGroupOfElements *group = _groups.getGhostGroup(i);
int nbNodes = group->getNbNodes();
int nbElements = group->getNbElements();
+ _groupId[group] = i+_groups.getNbElementGroups();
_dataProxys[i+_groups.getNbElementGroups()] = new fullMatrix<double> (&(*_ghostData)(offset),nbNodes, _nbFields*nbElements);
offset += nbNodes*_nbFields*nbElements;
}
@@ -142,6 +144,12 @@ void dgDofContainer::scatter() {
sol.setAll(recvProxy);
}
}
+void dgDofContainer::setAll(double v) {
+ for(int i=0;i<_data->size();i++)
+ (*_data)(i)=v;
+ for(int i=0;i<_ghostData->size();i++)
+ (*_ghostData)(i)=v;
+}
void dgDofContainer::scale(double f)
{
_data->scale(f);
diff --git a/Solver/dgDofContainer.h b/Solver/dgDofContainer.h
index 8eaf22b0e3..6c21641777 100644
--- a/Solver/dgDofContainer.h
+++ b/Solver/dgDofContainer.h
@@ -4,6 +4,7 @@
#include <vector>
#include "fullMatrix.h"
class dgGroupCollection;
+class dgGroupOfElements;
struct dgDofContainer {
private:
@@ -21,16 +22,19 @@ private:
int *countSend,*countRecv,*shiftSend,*shiftRecv,*groupShiftRecv;
double *sendBuf, *recvBuf;
std::vector<fullMatrix<double> *> _dataProxys; // proxys
+ std::map<const dgGroupOfElements*,int> _groupId;
public:
void scale(double f);
double norm();
void axpy(dgDofContainer &x, double a=1.);
inline fullMatrix<double> &getGroupProxy(int gId){ return *(_dataProxys[gId]); }
+ inline fullMatrix<double> &getGroupProxy(const dgGroupOfElements* g){ return *(_dataProxys[_groupId[g]]); }
dgDofContainer (dgGroupCollection &groups, int nbFields);
~dgDofContainer ();
int getNbElements() {return _totalNbElements;}
void scatter();
void save(const std::string name);
void load(const std::string name);
+ void setAll(double v);
};
#endif
diff --git a/Solver/dgSlopeLimiter.cpp b/Solver/dgSlopeLimiter.cpp
index adcdd9355d..6f1009f33b 100644
--- a/Solver/dgSlopeLimiter.cpp
+++ b/Solver/dgSlopeLimiter.cpp
@@ -6,89 +6,101 @@
//----------------------------------------------------------------------------------
bool dgSlopeLimiter::apply ( dgDofContainer &solution, dgGroupCollection &groups)
{
- //WARNING: ONLY FOR 1 GROUP OF FACES
- //TODO: make this more general
-
- dgGroupOfFaces* group = groups.getFaceGroup(0);
- fullMatrix<double> &solleft = solution.getGroupProxy(0);
- fullMatrix<double> &solright = solution.getGroupProxy(0);
+ solution.scatter();
int nbFields =_claw->nbFields();
- int totNbElems = solution.getNbElements();
// first compute max and min of all fields for all stencils
//----------------------------------------------------------
+ dgDofContainer MIN(groups, nbFields);
+ dgDofContainer MAX(groups, nbFields);
- fullMatrix<double> MIN (totNbElems,nbFields );
- fullMatrix<double> MAX (totNbElems, nbFields);
MIN.setAll ( 1.e22);
MAX.setAll (-1.e22);
int iElementL, iElementR, fSize;
- for(int iFace=0 ; iFace<group->getNbElements();iFace++) {
-
- iElementL = group->getElementLeftId(iFace);
- iElementR = group->getElementRightId(iFace);
-
- fullMatrix<double> TempL, TempR;
- TempL.setAsProxy(solleft, nbFields*iElementL, nbFields );
- TempR.setAsProxy(solright, nbFields*iElementR, nbFields );
-
- fSize = TempL.size1();
- for (int k=0; k< nbFields; ++k){
- double AVGL = 0;
- double AVGR = 0;
- for (int i=0; i<fSize; ++i) {
- AVGL += TempL(i,k);
- AVGR += TempR(i,k);
- }
- AVGL /= (double) fSize;
- AVGR /= (double) fSize;
- MIN ( iElementL , k ) = std::min ( AVGR , MIN ( iElementL , k ) );
- MAX ( iElementL , k ) = std::max ( AVGR , MAX ( iElementL , k ) );
- MIN ( iElementR , k ) = std::min ( AVGL , MIN ( iElementR , k ) );
- MAX ( iElementR , k ) = std::max ( AVGL , MAX ( iElementR , k ) );
- }
+ fullMatrix<double> TempL, TempR;
+ for( int iGFace=0; iGFace<groups.getNbFaceGroups(); iGFace++) {
+ dgGroupOfFaces* group = groups.getFaceGroup(iGFace);
+ const dgGroupOfElements *groupLeft = &group->getGroupLeft();
+ const dgGroupOfElements *groupRight = &group->getGroupRight();
+
+ fullMatrix<double> &solleft = solution.getGroupProxy(groupLeft);
+ fullMatrix<double> &solright = solution.getGroupProxy(groupRight);
+ fullMatrix<double> &MINLeft = MIN.getGroupProxy(groupLeft);
+ fullMatrix<double> &MAXLeft = MAX.getGroupProxy(groupLeft);
+ fullMatrix<double> &MINRight = MIN.getGroupProxy(groupRight);
+ fullMatrix<double> &MAXRight = MAX.getGroupProxy(groupRight);
+
+ for(int iFace=0 ; iFace<group->getNbElements();iFace++) {
+
+ iElementL = group->getElementLeftId(iFace);
+ iElementR = group->getElementRightId(iFace);
+
+ TempL.setAsProxy(solleft, nbFields*iElementL, nbFields );
+ TempR.setAsProxy(solright, nbFields*iElementR, nbFields );
+
+ fSize = TempL.size1();
+ for (int k=0; k< nbFields; ++k){
+ double AVGL = 0;
+ double AVGR = 0;
+ for (int i=0; i<fSize; ++i) {
+ AVGL += TempL(i,k);
+ AVGR += TempR(i,k);
+ }
+ AVGL /= (double) fSize;
+ AVGR /= (double) fSize;
+ MINLeft ( iElementL , k ) = std::min ( AVGR , MINLeft ( iElementL , k ) );
+ MAXLeft ( iElementL , k ) = std::max ( AVGR , MAXLeft ( iElementL , k ) );
+ MINRight ( iElementR , k ) = std::min ( AVGL , MINRight ( iElementR , k ) );
+ MAXRight ( iElementR , k ) = std::max ( AVGL , MAXRight ( iElementR , k ) );
+ }
+ }
}
//----------------------------------------------------------
// then limit the solution
//----------------------------------------------------------
- for (int iElement=0 ; iElement<totNbElems; ++iElement) {
+ for (int iGroup=0 ; iGroup<groups.getNbElementGroups() ; iGroup++) {
+ dgGroupOfElements &group = *groups.getElementGroup(iGroup);
+ fullMatrix<double> &sol = solution.getGroupProxy(iGroup);
+ fullMatrix<double> &MAXG = MAX.getGroupProxy(iGroup);
+ fullMatrix<double> &MING = MIN.getGroupProxy(iGroup);
fullMatrix<double> Temp;
- Temp.setAsProxy(solleft, nbFields*iElement, nbFields );
- for (int k=0; k<nbFields; ++k)
- {
- double AVG = 0.;
- double locMax = -1.e22;
- double locMin = 1.e22;
- double neighMax = MAX (iElement,k);
- double neighMin = MIN (iElement,k);
- for (int i=0; i<fSize; ++i)
+ for (int iElement=0 ; iElement<group.getNbElements() ; ++iElement) {
+ Temp.setAsProxy(sol, nbFields*iElement, nbFields );
+ for (int k=0; k<nbFields; ++k)
{
- AVG += Temp(i,k);
- locMax = std::max (locMax, Temp (i,k));
- locMin = std::min (locMin, Temp (i,k));
- }
- AVG /= (double) fSize;
-
- //SLOPE LIMITING DG
- //-------------------
- for (int i=0; i<fSize; ++i)Temp(i,k) -= AVG;
-
- double slopeLimiterValue = 1.0;
- if (locMax != AVG && locMax > neighMax) slopeLimiterValue = (neighMax-AVG) / (locMax-AVG);
- if (locMin != AVG && locMin < neighMin) slopeLimiterValue = std::min ( slopeLimiterValue , (AVG-neighMin) / (AVG-locMin) );
- if (AVG < neighMin) slopeLimiterValue = 0;
- if (AVG > neighMax) slopeLimiterValue = 0;
-
- // if (slopeLimiterValue != 1.0) printf("LIMTING %g\n",slopeLimiterValue);
- // slopeLimiterValue = 0.0;
-
- for (int i=0; i<fSize; ++i) Temp(i,k) = AVG + Temp(i,k)*slopeLimiterValue;
+ double AVG = 0.;
+ double locMax = -1.e22;
+ double locMin = 1.e22;
+ double neighMax = MAXG (iElement,k);
+ double neighMin = MING (iElement,k);
+ for (int i=0; i<fSize; ++i)
+ {
+ AVG += Temp(i,k);
+ locMax = std::max (locMax, Temp (i,k));
+ locMin = std::min (locMin, Temp (i,k));
+ }
+ AVG /= (double) fSize;
- }
- }
+ //SLOPE LIMITING DG
+ //-------------------
+ for (int i=0; i<fSize; ++i)Temp(i,k) -= AVG;
+
+ double slopeLimiterValue = 1.0;
+ if (locMax != AVG && locMax > neighMax) slopeLimiterValue = (neighMax-AVG) / (locMax-AVG);
+ if (locMin != AVG && locMin < neighMin) slopeLimiterValue = std::min ( slopeLimiterValue , (AVG-neighMin) / (AVG-locMin) );
+ if (AVG < neighMin) slopeLimiterValue = 0;
+ if (AVG > neighMax) slopeLimiterValue = 0;
+
+ // if (slopeLimiterValue != 1.0) printf("LIMTING %g\n",slopeLimiterValue);
+ // slopeLimiterValue = 0.0;
+
+ for (int i=0; i<fSize; ++i) Temp(i,k) = AVG + Temp(i,k)*slopeLimiterValue;
+ }
+ }
+ }
// --- CLIPPING: check unphysical values
for (int iG = 0; iG < groups.getNbElementGroups(); iG++){
dgGroupOfElements* egroup = groups.getElementGroup(iG);
@@ -100,12 +112,11 @@ bool dgSlopeLimiter::apply ( dgDofContainer &solution, dgGroupCollection &groups
dataCacheDouble *solutionEClipped = _claw->newClipToPhysics(cacheMap);
if (solutionEClipped){
for (int iElement=0 ; iElement<egroup->getNbElements() ;++iElement) {
- solutionE.setAsProxy(solGroup, iElement*nbFields, nbFields );
- solutionE.set((*solutionEClipped)());
- }
+ solutionE.setAsProxy(solGroup, iElement*nbFields, nbFields );
+ solutionE.set((*solutionEClipped)());
+ }
}
}
return true;
-
}
--
GitLab