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

No commit message 

No commit message
parent c3dd1649
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 132 additions and 107 deletions
Common/CreateFile.cpp
+ 4
− 4
View file @ e7526c1e
......@@ -171,16 +171,16 @@ void CreateOutputFile(std::string fileName, int format)
break;
case FORMAT_MSH:
GModel::current()->writeMSH
(fileName, CTX::instance()->mesh.mshFileVersion,
CTX::instance()->mesh.binary, CTX::instance()->mesh.saveAll,
CTX::instance()->mesh.saveParametric, CTX::instance()->mesh.scalingFactor);
if (CTX::instance()->mesh.saveDistance){
GModel::current()->writeDistanceMSH
(fileName, CTX::instance()->mesh.mshFileVersion,
CTX::instance()->mesh.binary, CTX::instance()->mesh.saveAll,
CTX::instance()->mesh.saveParametric, CTX::instance()->mesh.scalingFactor);
}
GModel::current()->writeMSH
(fileName, CTX::instance()->mesh.mshFileVersion,
CTX::instance()->mesh.binary, CTX::instance()->mesh.saveAll,
CTX::instance()->mesh.saveParametric, CTX::instance()->mesh.scalingFactor);
break;
case FORMAT_STL:
......
Numeric/polynomialBasis.cpp
+ 0
− 9
View file @ e7526c1e
......@@ -677,15 +677,6 @@ static void getFaceClosure(int iFace, int iSign, int iRotate, std::vector<int> &
default:
int face[4][3] = {{-3, -2, -1}, {1, -6, 4}, {-4, 5, 3}, {6, 2, -5}};
int order1node[4][3] = {{0, 2, 1}, {0, 1, 3}, {0, 3, 2}, {3, 1, 2}};
// int facedofstmp[200];
// face 0 | 0 1 2
// face 1 | 0 4 -3
// face 2 | 1 5 -4
// face 3 | -2 5 -3
// edge 0 | 4 -> 4+order-2
// edge 1 | 4+order-1 -> 4 + 2*(order-1)-1
// edge 2 | 4+2*(order-1)-> 4+ 3*(order-1)-1
// edge k | 4+k*(order-1) -> 4+(k+1)*(order-1)-1
for (int i = 0; i < 3; ++i){
int k = (3 + (iSign * i) + iRotate) % 3; //- iSign * iRotate
closure[i] = order1node[iFace][k];
......
Solver/TESTCASES/ForwardFacingStep.lua
+ 2
− 2
View file @ e7526c1e
MACH = 3;
MACH = 3.0;
GAMMA = 1.4;
U = 3.0
V = 0.0
......@@ -26,6 +26,7 @@ end
--[[
Example of a lua program driving the DG code
--]]
order = 1
print'*** Loading the mesh and the model ***'
myModel = GModel ()
......@@ -56,7 +57,6 @@ print'*** export ***'
DG:exportSolution('output/solution_0')
print'*** solve ***'
CFL = 2
for i=1,5000 do
......
Solver/dgAlgorithm.cpp
+ 1
− 6
View file @ e7526c1e
......@@ -268,7 +268,7 @@ void dgAlgorithm::rungeKutta (const dgConservationLaw &claw, // conservation l
std::vector<dgGroupOfElements*> &eGroups, // group of elements
std::vector<dgGroupOfFaces*> &fGroups, // group of interfacs
std::vector<dgGroupOfFaces*> &bGroups, // group of boundaries
double h, // time-step
double h, // time-step
dgDofContainer &sol,
dgDofContainer &resd,
dgLimiter *limiter,
......@@ -305,9 +305,6 @@ void dgAlgorithm::rungeKutta (const dgConservationLaw &claw, // conservation l
if(j){
K._data->scale(b[j]);
K._data->axpy(*(sol._data));
if (limiter){
// limiter->apply(K, eGroups, fGroups);
}
}
if (limiter){
......@@ -325,7 +322,6 @@ void dgAlgorithm::rungeKutta (const dgConservationLaw &claw, // conservation l
}
}
Unp._data->axpy(*(K._data),a[j]);
//if (limiter) limiter->apply(Unp, eGroups, fGroups);
}
if (limiter) limiter->apply(Unp, eGroups, fGroups);
for (int i=0;i<sol._dataSize;i++){
......@@ -546,7 +542,6 @@ void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (mayb
}
}
}
// ----- 3 ---- do the redistribution at face nodes using BLAS3
residual.gemm(group.getRedistributionMatrix(),NormalFluxQP);
......
Solver/dgConservationLawAdvection.h
+ 2
− 0
View file @ e7526c1e
#ifndef _DG_CONSERVATION_LAW_ADVECTION_H
#define _DG_CONSERVATION_LAW_ADVECTION_H
#include "dgConservationLaw.h"
//Advection diffusion equation
//dc/dt + v div(c) - nu lapl(c) = 0
class dgConservationLawAdvection : public dgConservationLaw {
std::string _vFunctionName,_nuFunctionName;
class advection;
......
Solver/dgConservationLawPerfectGas.cpp
+ 48
− 33
View file @ e7526c1e
......@@ -213,11 +213,14 @@ static inline void _ROE2D (const double &_GAMMA,
const double *solR,
double *FLUX){
//sol cons var: rho, rhou, rhov, rhoE
//u prim var : rho, u, v, p
const double gamma = _GAMMA;
const double GM1 = gamma - 1.;
double uL[4], uR[4];
const double overRhoL = 1./solL[0];
const double overRhoL = 1./solL[0];
uL[0] = solL[0];
uL[1] = solL[1]*overRhoL;
uL[2] = solL[2]*overRhoL;
......@@ -228,38 +231,33 @@ static inline void _ROE2D (const double &_GAMMA,
uR[1] = solR[1]*overRhoR;
uR[2] = solR[2]*overRhoR;
uR[3] = GM1*(solR[3] - 0.5* (solR[1]*solR[1]+solR[2]*solR[2])*overRhoR);
// central contributions
double halfrhoun;
/* --- left contributions ---*/
//* --- central contributions ---*/
double halfrhoun;
// --- left contributions
halfrhoun = 0.5*(uL[0]*(uL[1]*nx + uL[2]*ny));
double HL = gamma/GM1* uL[3]/uL[0]+0.5*(uL[1]*uL[1]+
uL[2]*uL[2]);
double HL = gamma/GM1* uL[3]/uL[0]+0.5*(uL[1]*uL[1]+uL[2]*uL[2]);
FLUX[0] = halfrhoun;
FLUX[1] = halfrhoun*uL[1] + .5*uL[3]*nx;
FLUX[2] = halfrhoun*uL[2] + .5*uL[3]*ny;
FLUX[3] = halfrhoun*HL;
/* --- right contributions ---*/
// --- right contributions
halfrhoun = 0.5*(uR[0]*(uR[1]*nx+uR[2]*ny));
double HR = gamma/GM1* uR[3]/uR[0]+0.5*(uR[1]*uR[1]+uR[2]*uR[2]);
double HR = gamma/GM1* uR[3]/uR[0]+0.5*(uR[1]*uR[1]+uR[2]*uR[2]);
FLUX[0] += halfrhoun;
FLUX[1] += halfrhoun*uR[1] + .5*uR[3]*nx;
FLUX[2] += halfrhoun*uR[2] + .5*uR[3]*ny;
FLUX[3] += halfrhoun*HR;
/* --- add dissipation ---*/
/* --- add rhoe dissipation ---*/
double sqr_rhoL = sqrt(uL[0]);
double sqr_rhoR = sqrt(uR[0]);
double sqr_rhoR = sqrt(uR[0]);
double invz1 = 1./ (sqr_rhoL + sqr_rhoR);
// double rho = sqr_rhoL * sqr_rhoR;
//rhoe average state
double u = ( sqr_rhoL* uL[1] + sqr_rhoR * uR[1] ) * invz1;
double v = ( sqr_rhoL* uL[2] + sqr_rhoR * uR[2] ) * invz1;
double H = ( sqr_rhoL* HL + sqr_rhoR * HR ) * invz1;
......@@ -281,17 +279,14 @@ static inline void _ROE2D (const double &_GAMMA,
double g1OnC2 = g1*oC2;
double TtOnT = (1.0 - u2*g1OnC2);
// matrix of left eigenvectors
// matrix of left eigenvectors
double L[16] = {
nx*TtOnT ,nx*u*g1OnC2 ,nx*v*g1OnC2 , -nx*g1OnC2, // L1
- tet , ny ,-nx , 0, // L3
g1*u2 - c*un , c*nx - g1*u , c*ny - g1*v, g1, // L3
g1*u2 + c*un ,-c*nx - g1*u ,-c*ny - g1*v, g1}; // L4
// characteristic decomposition of differences
double dW[4] = {0,0,0,0};
int idx = 0;
for (int i=0;i<4;i++)
......@@ -299,7 +294,6 @@ static inline void _ROE2D (const double &_GAMMA,
dW[i] += L[idx++]*dU[j];
// matrix of right eigenvectors
double R[16] = {
//R1 //R2 //R3 //R4
nx ,0 ,0.5*oC2 ,0.5*oC2,
......@@ -307,11 +301,9 @@ static inline void _ROE2D (const double &_GAMMA,
v*nx ,- nx ,0.5*(ny*oC + v*oC2),0.5*(-ny*oC + v*oC2),
u2*nx,tet ,0.5*(un*oC + H*oC2),0.5*(-un*oC + H*oC2)};
// eigenvalues
// KH : shouldn't we take into account an entropy correction ?
// absorb half the surface : scaling wrt central term
// TODO : shouldn't we take into account an entropy correction ?
// absorb half the surface : scaling wrt central term
const double A = 0.5;
double eps = 1.e-6;
......@@ -330,8 +322,7 @@ static inline void _ROE2D (const double &_GAMMA,
lflux -= lA[j]*dW[j]*R[index++];
FLUX[k] = -lflux;
}
}
// perfect gas law, GAMMA is the only parameter
};
class dgPerfectGasLaw2d::advection : public dataCacheDouble {
dataCacheDouble &sol;
......@@ -345,7 +336,6 @@ class dgPerfectGasLaw2d::advection : public dataCacheDouble {
_value=fullMatrix<double>(nQP,8);
const double GM1 = GAMMA - 1.0;
for (size_t k = 0 ; k < nQP; k++ ){
// printf("%d %g %g %g %g\n",k,sol(k,0),sol(k,1),sol(k,2),sol(k,3));
const double invrho = 1./sol(k,0);
const double q12 = sol(k,1)*sol(k,2)*invrho;
......@@ -365,10 +355,6 @@ class dgPerfectGasLaw2d::advection : public dataCacheDouble {
_value(k,2+4) = q22+p;
_value(k,3+4) = sol(k,2)*qq;
/* _value(k,8) = 0;
_value(k,9) = 0;
_value(k,10) = 0;
_value(k,11) = 0;*/
}
}
};
......@@ -455,6 +441,32 @@ public:
}
};
class dgPerfectGasLaw2d::clipToPhysics : public dataCacheDouble {
dataCacheDouble &sol;
public:
clipToPhysics(dataCacheMap &cacheMap):
sol(cacheMap.get("Solution",this))
{};
void _eval () {
double rhomin = 1.e-3;
double presmin= 1.e-3;
const int nQP = sol().size1();
for (size_t k = 0 ; k < nQP; k++ ){
_value(k,0) = sol(k,0);
_value(k,1) = sol(k,1);
_value(k,2) = sol(k,1);
_value(k,3) = sol(k,3);
if (sol(k,0) < rhomin)
_value(k,0) = rhomin;
double rhoV2 = sol(k,1)*sol(k,1)+sol(k,2)*sol(k,2);
rhoV2 /= sol(k,0);
const double p = (GAMMA-1)*(sol(k,3) - 0.5*rhoV2);
if (p < presmin)
_value(k,3) = presmin / (GAMMA-1) + 0.5 *rhoV2 ;
}
}
};
class dgPerfectGasLaw2d::riemann : public dataCacheDouble {
dataCacheDouble &normals, &solL, &solR;
......@@ -621,6 +633,9 @@ dataCacheDouble *dgPerfectGasLaw2d::newSourceTerm (dataCacheMap &cacheMap) const
else
return new source(cacheMap,_sourceFunctionName);
}
dataCacheDouble *dgPerfectGasLaw2d::newClipToPhysics( dataCacheMap &cacheMap) const {
return new clipToPhysics(cacheMap);
}
dgPerfectGasLaw2d::dgPerfectGasLaw2d()
{
......
Solver/dgConservationLawPerfectGas.h
+ 7
− 0
View file @ e7526c1e
#ifndef _DG_CONSERVATION_LAW_PERFECT_GAS_H_
#define _DG_CONSERVATION_LAW_PERFECT_GAS_H_
#include "dgConservationLaw.h"
// Compressible Navier-Stokes equations
// d(rho)/dt + d(rhou)/dx + d(rhov)/dy = 0
// d(rhou)/dt + d(rhou^2+p)/dx + d(rhouv)/dy + d(t_xx)/dx + d(t_xy)/dy= 0
// d(rhov)/dt + d(rhouv)/dx + d(rhov^2+p)/dy + d(t_xy)/dx + d(t_yy)/dy= 0
// d(rhoE)/dt + d(rhouH)/dx + d(rhovH)/dy + d(ut_xx+vt_xy-q_x)/dx + d(ut_xy+vt_yy-q_y)/dy= 0
class dgPerfectGasLaw2d : public dgConservationLaw {
class advection;
class diffusion;
......@@ -9,6 +14,7 @@ class dgPerfectGasLaw2d : public dgConservationLaw {
class source;
class maxConvectiveSpeed;
class maxDiffusivity;
class clipToPhysics;
// the name of the functions for
// viscosity (_muFunctionName)
// thermal conductivity (_kappaFunctionName)
......@@ -22,6 +28,7 @@ class dgPerfectGasLaw2d : public dgConservationLaw {
dataCacheDouble *newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const;
dataCacheDouble *newDiffusiveFlux( dataCacheMap &cacheMap) const;
dataCacheDouble *newSourceTerm (dataCacheMap &cacheMap) const;
dataCacheDouble *newClipToPhysics (dataCacheMap &cacheMap) const;
dgPerfectGasLaw2d();
dgBoundaryCondition *newWallBoundary() ;
dgBoundaryCondition *newSlipWallBoundary() ;
......
Solver/dgConservationLawWaveEquation.cpp
+ 1
− 3
View file @ e7526c1e
#include "dgConservationLaw.h"
#include "dgConservationLawWaveEquation.h"
#include "function.h"
// dp/dt - rho*c^2 div(u,v) = 0
// du/dt + 1/rho dp/dx = 0
// dv/dt + 1/rho dp/dy = 0
static double c=1;
static double rho=1;
......
Solver/dgConservationLawWaveEquation.h
+ 4
− 0
View file @ e7526c1e
#ifndef _DG_CONSERVATION_LAW_WAVE_EQUATION_H_
#define _DG_CONSERVATION_LAW_WAVE_EQUATION_H_
#include "dgConservationLaw.h"
// Linear acoustics
// dp/dt - rho*c^2 div(u,v) = 0
// du/dt + 1/rho dp/dx = 0
// dv/dt + 1/rho dp/dy = 0
class methodBinding;
class constructorBinding;
class dgConservationLawWaveEquation : public dgConservationLaw {
......
Solver/dgLimiter.h
+ 5
− 2
View file @ e7526c1e
......@@ -6,17 +6,20 @@
struct dgDofContainer;
class dgGroupOfElements;
class dgGroupOfFaces;
class dgConservationLaw;
class dgLimiter{
protected:
dgConservationLaw *_claw;
public:
dgLimiter () {}
dgLimiter (dgConservationLaw *claw) : _claw(claw) {}
virtual bool apply ( dgDofContainer &sol, std::vector<dgGroupOfElements*> &eGroups,
std::vector<dgGroupOfFaces*> &fGroup ) = 0;
};
class dgSlopeLimiter : public dgLimiter{
public :
dgSlopeLimiter () : dgLimiter () {}
dgSlopeLimiter (dgConservationLaw *claw) : dgLimiter (claw) {}
virtual bool apply ( dgDofContainer &solution,
std::vector<dgGroupOfElements*> &eGroups,
std::vector<dgGroupOfFaces*> &fGroup);
......
Solver/dgSlopeLimiter.cpp
+ 55
− 43
View file @ e7526c1e
#include "dgLimiter.h"
#include "dgGroupOfElements.h"
#include "dgSystemOfEquations.h"
#include "function.h"
//----------------------------------------------------------------------------------
bool dgSlopeLimiter::apply ( dgDofContainer &solution,
......@@ -13,13 +13,11 @@ bool dgSlopeLimiter::apply ( dgDofContainer &solution,
//TODO: make this more general
dgGroupOfFaces* group = fGroups[0];
fullMatrix<double> &SolLeft = *(solution._dataProxys[0]);
fullMatrix<double> &SolRight = *(solution._dataProxys[0]);
int nbFields = solution.getNbFields();
fullMatrix<double> &solleft = solution.getGroupProxy(0); //*(solution._dataProxys[0]);
fullMatrix<double> &solright = solution.getGroupProxy(0); //*(solution._dataProxys[0]);
int nbFields =_claw->nbFields();
int totNbElems = solution.getNbElements();
// --- CLIPPING: check unphysical values
// first compute max and min of all fields for all stencils
//----------------------------------------------------------
......@@ -29,14 +27,14 @@ bool dgSlopeLimiter::apply ( dgDofContainer &solution,
MAX.setAll (-1.e22);
int iElementL, iElementR, fSize;
for(int iFace=0 ; iFace<group->getNbElements();iFace++)
{
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 );
TempL.setAsProxy(solleft, nbFields*iElementL, nbFields );
TempR.setAsProxy(solright, nbFields*iElementR, nbFields );
fSize = TempL.size1();
for (int k=0; k< nbFields; ++k){
......@@ -55,16 +53,13 @@ bool dgSlopeLimiter::apply ( dgDofContainer &solution,
}
}
// printf("fSize = %d\n",fSize);
//----------------------------------------------------------
//----------------------------------------------------------
// then limit the solution
//----------------------------------------------------------
for (int iElement=0 ; iElement<totNbElems; ++iElement)
{
for (int iElement=0 ; iElement<totNbElems; ++iElement) {
fullMatrix<double> Temp;
Temp.setAsProxy(SolRight, nbFields*iElement, nbFields );
Temp.setAsProxy(solleft, nbFields*iElement, nbFields );
for (int k=0; k<nbFields; ++k)
{
double AVG = 0.;
......@@ -90,9 +85,7 @@ bool dgSlopeLimiter::apply ( dgDofContainer &solution,
if (AVG < neighMin) slopeLimiterValue = 0;
if (AVG > neighMax) slopeLimiterValue = 0;
// if (slopeLimiterValue != 1.0){
// printf("LIMTING %g\n",slopeLimiterValue);
// }
// 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;
......@@ -100,31 +93,50 @@ bool dgSlopeLimiter::apply ( dgDofContainer &solution,
}
}
#if 0
double rhomin = 1.e-3;
double presmin= 1.e-3;
for (int iElement=0;iElement<totNbElems;iElement++){
fullMatrix<double> solElem;
solElem.setAsProxy(SolRight, nbFields*iElement, nbFields );
for (int k=0;k< solElem.size1() ;k++){
if (solElem(k,0) < rhomin) {
solElem(k,0) = rhomin;
}
double rhoV2 = 0;
for (int j=0;j<2;j++) {
double rhov = solElem(k,j+1);
rhoV2 += rhov*rhov;
}
rhoV2 /= solElem(k,0);
const double p = (1.4-1)*solElem(k,3) - 0.5*(1.4-1)* rhoV2;
if (p < presmin) {
solElem(k,3) = 0.5 *rhoV2 + presmin / (1.4-1);
// printf("negative pressure %g cliiped !!\n",p);
}
}
}
#endif
// --- CLIPPING: check unphysical values
// fullMatrix<double> solutionQP (group->getNbIntegrationPoints(),group->getNbElements() * nbFields);
// group->getCollocationMatrix().mult(solleft, solutionQP);
// //dataCacheDouble &solutionQPe = cacheMap.provideData("Solution");
// //solutionQPe.set(fullMatrix<double>(group.getNbIntegrationPoints(),nbFields));
// //dataCacheElement &cacheElement = cacheMap.getElement();
// dataCacheMap cacheMap(group->getNbIntegrationPoints());
// dataCacheDouble *clipping = _claw->newClipToPhysics(cacheMap);
// fullMatrix<double> clippedSol;
// for (int iElement=0 ; iElement<group->getNbElements() ;++iElement) {
// clippedSol.setAsProxy(solutionQP, iElement*nbFields, nbFields );
// for (int iPt =0; iPt< group->getNbIntegrationPoints(); iPt++) {
// for (int k=0;k<nbFields;k++)
// clippedSol(iPt,k) = (*clipping)(iPt,k);
// }
// //cacheElement.set(group.getElement(iElement));
// }
// solleft.gemm(group->getRedistributionMatrix(),solutionQP);
//#if 0
// double rhomin = 1.e-3;
// double presmin= 1.e-3;
// for (int iElement=0;iElement<totNbElems;iElement++){
// fullMatrix<double> solElem;
// solElem.setAsProxy(solleft, nbFields*iElement, nbFields );
// for (int k=0;k< solElem.size1() ;k++){
// if (solElem(k,0) < rhomin) {
// solElem(k,0) = rhomin;
// }
// double rhoV2 = 0;
// for (int j=0;j<2;j++) {
// double rhov = solElem(k,j+1);
// rhoV2 += rhov*rhov;
// }
// rhoV2 /= solElem(k,0);
// const double p = (1.4-1)*solElem(k,3) - 0.5*(1.4-1)* rhoV2;
// if (p < presmin) {
// solElem(k,3) = 0.5 *rhoV2 + presmin / (1.4-1);
// }
// }
// }
//#endif
return true;
}
Solver/dgSystemOfEquations.cpp
+ 3
− 4
View file @ e7526c1e
......@@ -103,8 +103,7 @@ double dgSystemOfEquations::RK44(double dt){
return _solution->_data->norm();
}
double dgSystemOfEquations::computeInvSpectralRadius(){
double dgSystemOfEquations::computeInvSpectralRadius(){
double sr = 1.e22;
for (int i=0;i<_elementGroups.size();i++){
std::vector<double> DTS;
......@@ -115,7 +114,7 @@ double dgSystemOfEquations::computeInvSpectralRadius(){
}
double dgSystemOfEquations::RK44_limiter(double dt){
dgLimiter *sl = new dgSlopeLimiter();
dgLimiter *sl = new dgSlopeLimiter(_claw);
_algo->rungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups, dt, *_solution, *_rightHandSide, sl);
delete sl;
return _solution->_data->norm();
......@@ -135,7 +134,7 @@ void dgSystemOfEquations::exportSolution(std::string outputFile){
}
void dgSystemOfEquations::limitSolution(){
dgLimiter *sl = new dgSlopeLimiter();
dgLimiter *sl = new dgSlopeLimiter(_claw);
sl->apply(*_solution,_elementGroups,_faceGroups);
delete sl;
......
Solver/dgSystemOfEquations.h
+ 0
− 1
View file @ e7526c1e
......@@ -23,7 +23,6 @@ public:
dgDofContainer (std::vector<dgGroupOfElements*> &groups, const dgConservationLaw &claw);
~dgDofContainer ();
int getNbElements() {return totalNbElements;}
int getNbFields() {return nbFields;}
};
class binding;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment