Skip to content
Snippets Groups Projects
Commit b3f6d3c1 authored by Emilie Marchandise's avatar Emilie Marchandise
Browse files

No commit message 

No commit message
parent 27761919
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Solver/CMakeLists.txt
+ 1
− 0
View file @ b3f6d3c1
......@@ -19,6 +19,7 @@ set(SRC
dgConservationLawShallowWater2d.cpp
dgConservationLawWaveEquation.cpp
dgConservationLawPerfectGas.cpp
SlopeLimiter.cpp
function.cpp
functionDerivator.cpp
luaFunction.cpp
......
Solver/SlopeLimiter.cc → Solver/SlopeLimiter.cpp
+ 139
− 0
View file @ b3f6d3c1
#include "dgLimiter.h"
#include "dgGroupOfElements.h"
#include "dgSystemOfEquations.h"
#include "dgConservationLaw.h"
// //----------------------------------------------------------------------------------
......@@ -10,43 +9,39 @@
// std::vector<dgGroupOfFaces*> &fGroups)
// {
// const DGIntMatrix &topo = domain->getFaceToRegion();
// const int fSize = solution->FunctionSpaceSize();
// const ClosedFunctionSpace* fsC = DGModelDataManager::instance().getFunctionSpace(problemName,0);
// // //--- Correct unphysical values
// // double Amin = 1.e-5;
// // for (int iE=0;iE<solution->getNbElements();iE++){
// // fullMatrix soljElem = solution->touchElementCoefs(iE);
// // for (int k=0;k< soljElem.size1() ;k++){
// // if (soljElem(k,0) < Amin) soljElem(k,0) = Amin;
// // }
// // }
// //WARNING FOR ONLY 1 GROUP OF FACES
// dgGroupOfFaces* group = fGroups[0];
// fullMatrix<double> SolLeft = solution._dataProxys[0];
// fullMatrix<double> SolRight = solution._dataProxys[0];
// int nbFields = solution->getNbFields();
// int totNbElems = solution->getNbElements();
// // first compute max and min of all fields for all stencils
// //----------------------------------------------------------
// fullMatrix MIN (solution->getNbElements(), solution->getNbFields());
// fullMatrix MAX (solution->getNbElements(), solution->getNbFields());
// fullMatrix<double> MIN (totNbElems,nbFields );
// fullMatrix<double> MAX (totNbElems, nbFields);
// MIN.set_all ( 1.e22);
// MAX.set_all (-1.e22);
// fullMatrix MEANL (solution->getNbElements(), solution->getNbFields());
// fullMatrix MEANR (solution->getNbElements(), solution->getNbFields());
// fullMatrix<double> MEANL (totNbElems, nbFields);
// fullMatrix<double> MEANR (totNbElems, nbFields);
// MEANL.set_all ( 0.01);
// MEANR.set_all ( 0.01);
// dgGroupOfElements* group = fGroups[0];
// for(int iFace=0 ; iFace<group->getNbElements();iFace++)
// {
// const int iElementL = group->getElementLeft(iFace);
// const int iElementR = group->getElementRight(iFace);
// int iElementL = group->getElementLeftId(iFace);
// int iElementR = group->getElementRightId(iFace);
// if (iElementR >= 0)
// {
// const fullMatrix TempR = fullMatrix(solution->seeElementCoefs (iElementR));
// const fullMatrix TempL = fullMatrix(solution->seeElementCoefs (iElementL));
// for (int k=0; k<nbFields; ++k)
// fullMatrix<double> TempL, TempR;
// TempL.setAsProxy(SolLeft, nbFields*iELementL, nbFields );
// TempR.setAsProxy(SolRight, nbFields*iELementR, nbFields );
// printf("TempL size 1 = %d %d ", TempL.size1(), TempL.size2());
// exit(1);
// int fSize = TempL.size2();
// for (int k=0; k< fSize; ++k)
// {
// double AVGL = 0;
// double AVGR = 0;
......@@ -68,7 +63,8 @@
// }
// }
// else{
// const fullMatrix TempL = fullMatrix(solution->seeElementCoefs (iElementL));
// fullMatrix<double> TempL, ;
// TempL.setAsProxy(SolLeft, nbFields*iELementL, nbFields );
// for (int k=0; k<nbFields; ++k){
// for (int i=0; i<fSize; ++i){
// MIN ( iElementL , k ) = std::min ( TempL(i,k) , MIN ( iElementL , k ) );
......@@ -90,53 +86,53 @@
// // then limit the solution
// //----------------------------------------------------------
// for (int iElement=0 ; iElement<solution->getNbElements(); ++iElement)
// {
// fullMatrix Temp = fullMatrix(solution->touchElementCoefs (iElement));
// 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)
// {
// AVG += Temp(i,k);
// locMax = std::max (locMax, Temp (i,k));
// locMin = std::min (locMin, Temp (i,k));
// }
// AVG /= (double) fSize;
// // for (int iElement=0 ; iElement<totNbElems; ++iElement)
// // {
// // fullMatrix Temp = fullMatrix(solution->touchElementCoefs (iElement));
// // 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)
// // {
// // 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;
// // //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;
// // 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 (detectDISC(iElement) == 0.0) slopeLimiterValue = 1.0; // do not limit
// //if (slopeLimiterValue != 1.0) printf("limit (%g) elem =%d \n", slopeLimiterValue, iElement);
// // if (detectDISC(iElement) == 0.0) slopeLimiterValue = 1.0; // do not limit
// // //if (slopeLimiterValue != 1.0) printf("limit (%g) elem =%d \n", slopeLimiterValue, iElement);
// for (int i=0; i<fSize; ++i) Temp(i,k) *= slopeLimiterValue;
// // for (int i=0; i<fSize; ++i) Temp(i,k) *= slopeLimiterValue;
// for (int i=0; i<fSize; ++i) Temp(i,k) += AVG;
// // for (int i=0; i<fSize; ++i) Temp(i,k) += AVG;
// //MINMOD DG (IF CHANGE TO THIS ADD LIMITER LINE IN RUNGEKUTTANNOPERATOR.cc)
// //--------------------------------
// // if (detectDISC(iElement) != 0.0){
// // const size_t nn = Temp.size1();
// // double alpha = 0.5; //0.5 1.0 1.3 (alpha=0.5 => MUSC Schem VAN LEER)
// // Temp(0,k) = AVG-minmod(AVG-Temp(0,k), alpha*(AVG-MEANL(iElement,k)), alpha*(MEANR(iElement,k)-AVG));
// // Temp(nn-1,k) = AVG+minmod(Temp(nn-1,k)-AVG, alpha*(AVG-MEANL(iElement,k)), alpha*(MEANR(iElement,k)-AVG));
// // }
// // //MINMOD DG (IF CHANGE TO THIS ADD LIMITER LINE IN RUNGEKUTTANNOPERATOR.cc)
// // //--------------------------------
// // // if (detectDISC(iElement) != 0.0){
// // // const size_t nn = Temp.size1();
// // // double alpha = 0.5; //0.5 1.0 1.3 (alpha=0.5 => MUSC Schem VAN LEER)
// // // Temp(0,k) = AVG-minmod(AVG-Temp(0,k), alpha*(AVG-MEANL(iElement,k)), alpha*(MEANR(iElement,k)-AVG));
// // // Temp(nn-1,k) = AVG+minmod(Temp(nn-1,k)-AVG, alpha*(AVG-MEANL(iElement,k)), alpha*(MEANR(iElement,k)-AVG));
// // // }
// }
// }
// // }
// // }
// return true;
......
Solver/dgAlgorithm.cpp
+ 2
− 2
View file @ b3f6d3c1
......@@ -313,14 +313,14 @@ void dgAlgorithm::rungeKutta (const dgConservationLaw &claw, // conservation l
}
}
Unp._data->axpy(*(K._data),a[j]);
if (limiter) limiter->apply(Unp);
if (limiter) limiter->apply(Unp, eGroups, fGroups);
}
for (int i=0;i<sol._dataSize;i++){
// printf("tempSol[%d] = %g\n",i,(*tempSol._data)(i));
// memcp
(*sol._data)(i)=(*Unp._data)(i);
//if (limiter) limiter->apply(sol);
//if (limiter) limiter->apply(sol, eGroups, fGroups);
}
}
......
Solver/dgLimiter.h
+ 11
− 12
View file @ b3f6d3c1
......@@ -3,24 +3,23 @@
#include "fullMatrix.h"
#include <vector>
struct dgDofContainer;
class dgGroupOfElements;
class dgGroupOfFaces;
class dgLimiter
{
class dgLimiter{
public:
dgLimiter () {}
virtual ~dgLimiter() {}
virtual bool apply ( dgDofContainer &sol ) =0;
virtual bool apply ( dgDofContainer &sol, std::vector<dgGroupOfElements*> &eGroups,
std::vector<dgGroupOfFaces*> &fGroup ) = 0;
};
class SlopeLimiter : public dgLimiter
{
class SlopeLimiter : public dgLimiter{
public :
SlopeLimiter () : dgLimiter () {}
~SlopeLimiter() {}
bool apply ( dgDofContainer &sol );
virtual bool apply ( dgDofContainer &solution,
std::vector<dgGroupOfElements*> &eGroups,
std::vector<dgGroupOfFaces*> &fGroup);
};
#endif
Solver/dgSystemOfEquations.cpp
+ 3
− 2
View file @ b3f6d3c1
......@@ -5,7 +5,7 @@
#include "MElement.h"
#include "PView.h"
#include "PViewData.h"
#include "dgLimiter.h"
class dgConservationLawL2Projection : public dgConservationLaw {
std::string _functionName;
......@@ -77,7 +77,8 @@ void dgSystemOfEquations::setup(){
double dgSystemOfEquations::RK44(double dt){
_algo->rungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups, dt, *_solution, *_rightHandSide);
//dgLimiter *sl = new SlopeLimiter();
_algo->rungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups, dt, *_solution, *_rightHandSide, NULL);
return _solution->_data->norm();
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment