diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt
index add58f8e9a3e9bd53b8808d1f652049259752aa1..fcb4bc40ae8e3759861cda8ef515b955bee16db0 100644
--- a/Solver/CMakeLists.txt
+++ b/Solver/CMakeLists.txt
@@ -18,6 +18,7 @@ set(SRC
dgSystemOfEquations.cpp
dgConservationLawShallowWater2d.cpp
dgConservationLawWaveEquation.cpp
+ dgConservationLawPerfectGas.cpp
function.cpp
functionDerivator.cpp
luaFunction.cpp
diff --git a/Solver/TESTCASES/BackwardFacingStep.lua b/Solver/TESTCASES/BackwardFacingStep.lua
new file mode 100644
index 0000000000000000000000000000000000000000..00bc95c90faa79457141cc7fdd07894e68e56332
--- /dev/null
+++ b/Solver/TESTCASES/BackwardFacingStep.lua
@@ -0,0 +1,66 @@
+MACH = 0.1;
+RHO = 1.0;
+PRES = 1./(MACH*RHO*RHO*1.4*1.4)
+V = 1.0
+SOUND = V/MACH
+
+--[[
+ Function for initial conditions
+--]]
+function free_stream( x, f )
+ FCT = fullMatrix(f)
+ XYZ = fullMatrix(x)
+ for i=0,XYZ:size1()-1 do
+ FCT:set(i,0,RHO)
+ FCT:set(i,1,RHO*V)
+ FCT:set(i,2,0.0)
+ FCT:set(i,3, 0.5*RHO*V*V+PRES/0.4)
+ end
+end
+
+--[[
+ Example of a lua program driving the DG code
+--]]
+
+order = 1
+print'*** Loading the mesh and the model ***'
+myModel = GModel ()
+myModel:load ('step.geo')
+myModel:load ('step.msh')
+print'*** Create a dg solver ***'
+DG = dgSystemOfEquations (myModel)
+DG:setOrder(order)
+DG:setConservationLaw('PerfectGas2d')
+DG:addBoundaryCondition('Walls','Wall')
+
+FS = createFunction.lua(4, 'free_stream', 'XYZ')
+
+DG:addBoundaryCondition('LeftRight','FreeStream',FS)
+DG:setup()
+
+
+print'*** setting the initial solution ***'
+
+DG:L2Projection(FS)
+
+print'*** export ***'
+
+DG:exportSolution('solution_0')
+
+print'*** solve ***'
+
+LC = 0.1
+dt = .3*LC/(SOUND+V);
+print('DT=',dt)
+
+for i=1,1000 do
+ norm = DG:RK44(dt)
+ print('*** ITER ***',i,norm)
+ if (i % 20 == 0) then
+ DG:exportSolution(string.format("solution-%03d", i))
+ end
+end
+
+print'*** done ***'
+
+
diff --git a/Solver/TESTCASES/WavePulse.lua b/Solver/TESTCASES/WavePulse.lua
index 1e0e93c4a9e3d21bf9fe20cffe062782daef20c8..3d33870d1f04b3cde3d4ca3ef53945793572daee 100644
--- a/Solver/TESTCASES/WavePulse.lua
+++ b/Solver/TESTCASES/WavePulse.lua
@@ -48,8 +48,7 @@ for i=1,1000 do
norm = DG:RK44(dt)
print('*** ITER ***',i,norm)
if (i % 10 == 0) then
- AA = string.format("output/solution-%03d", i)
- DG:exportSolution(AA)
+ DG:exportSolution(string.format("solution-%03d", i))
end
end
diff --git a/Solver/TESTCASES/square.geo b/Solver/TESTCASES/square.geo
index 4eb2bc249f14600a30f0e76a9360771bd9a76109..a214e082b0b163d953a2f4e77243284a87f3a9de 100644
--- a/Solver/TESTCASES/square.geo
+++ b/Solver/TESTCASES/square.geo
@@ -10,3 +10,5 @@ Line Loop(5) = {2, 3, 4, 1};
Plane Surface(6) = {5};
Physical Line("Border") = {1, 2, 3, 4};
Physical Surface("Inside") = {6};
+Transfinite Line {1, 2, 4, 3} = 2 Using Progression 1;
+Transfinite Surface {6};
diff --git a/Solver/TESTCASES/step.geo b/Solver/TESTCASES/step.geo
new file mode 100644
index 0000000000000000000000000000000000000000..62afa73c6a7e652daf7fd6af7765e2162eb72279
--- /dev/null
+++ b/Solver/TESTCASES/step.geo
@@ -0,0 +1,18 @@
+Point(1) = {0, 0, 0, .1};
+Point(2) = {3, 0, 0, .1};
+Point(3) = {3, 1, 0, .1};
+Point(4) = {5, 1, 0, .1};
+Point(5) = {5, 2, 0, .1};
+Point(6) = {0, 2, 0, .1};
+Line(1) = {6, 5};
+Line(2) = {5, 4};
+Line(3) = {4, 3};
+Line(4) = {3, 2};
+Line(5) = {2, 1};
+Line(6) = {1, 6};
+Line Loop(7) = {2, 3, 4, 5, 6, 1};
+Plane Surface(8) = {7};
+Physical Line("Walls") = {1, 3, 4, 5};
+Physical Line("LeftRight") = {2, 6};
+//Physical Line("LeftRight") = {1,2,3,4,5,6};
+Physical Surface("Body") = {8};
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index da3d0ebfa66d2dfd1b5e66fb84f30aec8cc7aa2c..0a56bb059136b495dbdc6ab18d189f2607a5d16f 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -418,6 +418,7 @@ void dgAlgorithm::residual( const dgConservationLaw &claw,
residu[i]->scale(0);
residualVolume(claw,*eGroups[i],*solution[i],*residu[i]);
}
+ // residu[0]->print("Volume");
//interface term
for(size_t i=0;i<fGroups.size() ; i++) {
dgGroupOfFaces &faces = *fGroups[i];
@@ -433,6 +434,7 @@ void dgAlgorithm::residual( const dgConservationLaw &claw,
residualInterface(claw,faces,solInterface,*solution[iGroupLeft], *solution[iGroupRight],residuInterface);
faces.mapFromInterface(nbFields, residuInterface, *residu[iGroupLeft], *residu[iGroupLeft]);
}
+ // residu[0]->print("Interfaces");
//boundaries
for(size_t i=0;i<bGroups.size() ; i++) {
dgGroupOfFaces &faces = *bGroups[i];
@@ -447,4 +449,5 @@ void dgAlgorithm::residual( const dgConservationLaw &claw,
residualBoundary(claw,faces,solInterface,*solution[iGroupLeft],residuInterface);
faces.mapFromInterface(nbFields, residuInterface, *residu[iGroupLeft], *residu[iGroupRight]);
}
+ // residu[0]->print("Boundaries");
}
diff --git a/Solver/dgConservationLaw.h b/Solver/dgConservationLaw.h
index 93df46b5757f789455fa5fbde78a2f24b4f6f10c..3e3c75a64244289adfc888266616b6ca2a80177f 100644
--- a/Solver/dgConservationLaw.h
+++ b/Solver/dgConservationLaw.h
@@ -54,7 +54,12 @@ class dgConservationLaw {
dgConservationLaw *dgNewConservationLawAdvection(const std::string vname);
dgConservationLaw *dgNewConservationLawShallowWater2d();
dgConservationLaw *dgNewConservationLawWaveEquation();
+dgConservationLaw *dgNewPerfectGasLaw2d();
+
dgBoundaryCondition *dgNewBoundaryConditionWaveEquationWall();
dgBoundaryCondition *dgNewBoundaryConditionShallowWater2dWall();
+dgBoundaryCondition *dgNewBoundaryConditionPerfectGasLaw2dWall();
+dgBoundaryCondition *dgNewBoundaryConditionPerfectGasLaw2dFreeStream(std::string&);
+
#endif
diff --git a/Solver/dgConservationLawPerfectGas.cpp b/Solver/dgConservationLawPerfectGas.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..eb0dae9ed034c3ebb5fde2e985b083d44027ac98
--- /dev/null
+++ b/Solver/dgConservationLawPerfectGas.cpp
@@ -0,0 +1,328 @@
+#include "dgConservationLaw.h"
+#include "function.h"
+
+const double GAMMA = 1.4;
+static inline void _ROE2D (const double &_GAMMA,
+ const double &nx,
+ const double &ny,
+ const double *solL,
+ const double *solR,
+ double *FLUX){
+
+ const double gamma = _GAMMA;
+ const double GM1 = gamma - 1.;
+
+ double uL[4], uR[4];
+ const double overRhoL = 1./solL[0];
+ uL[0] = solL[0];
+ uL[1] = solL[1]*overRhoL;
+ uL[2] = solL[2]*overRhoL;
+ uL[3] = GM1*(solL[3] - 0.5* (solL[1]*solL[1]+solL[2]*solL[2])*overRhoL);
+
+ const double overRhoR = 1./solR[0];
+ uR[0] = solR[0];
+ 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 ---*/
+ 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]);
+
+ 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 ---*/
+
+ 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]);
+
+ 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 ---*/
+
+ double sqr_rhoL = sqrt(uL[0]);
+ double sqr_rhoR = sqrt(uR[0]);
+
+ double invz1 = 1./ (sqr_rhoL + sqr_rhoR);
+
+ // double rho = sqr_rhoL * sqr_rhoR;
+ 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;
+
+ double dU[4];
+ for (int k=0;k<4;k++) dU[k] = solR[k] - solL[k];
+
+ double un = u*nx + v*ny ;
+ double tet = ny*u - nx*v;
+
+ double u2 = u*u + v*v;
+ double c2 = GM1 * ( H - 0.5 * u2);
+ double c = sqrt(c2);
+
+ double g1 = gamma - 1.;
+ double oC = 1./c;
+ double oC2 = oC*oC;
+
+ double g1OnC2 = g1*oC2;
+ double TtOnT = (1.0 - u2*g1OnC2);
+
+ // 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++)
+ for (int j=0;j<4;j++)
+ 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,
+ u*nx ,ny ,0.5*(nx*oC + u*oC2),0.5*(-nx*oC + u*oC2),
+ 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
+
+
+ const double A = 0.5;
+ double eps = 1.e-6;
+ double absUn = (fabs(un) + eps*c)*A;
+ double lA[4] = {absUn,
+ absUn,
+ fabs(un + c)*A,
+ (fabs(un - c)+eps*c)*A};
+
+ // add wave contributions to flux
+ int index = 0;
+
+ for (int k=0;k<4;k++) {
+ double lflux = FLUX[k];
+ for (int j=0;j<4;j++)
+ lflux -= lA[j]*dW[j]*R[index++];
+ FLUX[k] = -lflux;
+ }
+}
+
+class dgPerfectGasLaw2d : public dgConservationLaw {
+
+ // perfect gas law, GAMMA is the only parameter
+
+ class advection : public dataCacheDouble {
+ dataCacheDouble /
+ public:
+ advection(dataCacheMap &cacheMap):
+ sol(cacheMap.get("Solution",this))
+ {};
+ void _eval () {
+ const int nQP = sol().size1();
+ if(_value.size1() != nQP)
+ _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;
+ const double q11 = sol(k,1)*sol(k,1)*invrho;
+ const double q22 = sol(k,2)*sol(k,2)*invrho;
+
+ const double p = GM1*sol(k,3) - 0.5*GM1*(q11+q22);
+ const double qq = invrho*(sol(k,3)+p);
+
+ _value(k,0) = sol(k,1);
+ _value(k,1) = q11+p;
+ _value(k,2) = q12;
+ _value(k,3) = sol(k,1)*qq;
+
+ _value(k,0+4) = sol(k,2);
+ _value(k,1+4) = q12;
+ _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;*/
+ }
+ }
+ };
+
+ class riemann : public dataCacheDouble {
+ dataCacheDouble &normals, &solL, &solR;
+ public:
+ riemann(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight):
+ normals(cacheMapLeft.get("Normals", this)),
+ solL(cacheMapLeft.get("Solution", this)),
+ solR(cacheMapRight.get("Solution", this))
+ {};
+ void _eval () {
+ int nQP = solL().size1();
+ if(_value.size1() != nQP)
+ _value = fullMatrix<double>(nQP,8);
+
+ for(int i=0; i< nQP; i++) {
+ const double solLeft [4] = {solL(i,0),solL(i,1),solL(i,2),solL(i,3)};
+ const double solRight[4] = {solR(i,0),solR(i,1),solR(i,2),solR(i,3)};
+ double FLUX[4] ;
+ const double nx = normals(0,i);
+ const double ny = normals(1,i);
+ _ROE2D (GAMMA,nx,ny,solLeft,solRight,FLUX);
+
+ /*
+ printf("RSOLL %g %g %g %g\n",solLeft[0],solLeft[1],solLeft[2], solLeft[3]);
+ printf("RSOLR %g %g %g %g\n",solRight[0],solRight[1],solRight[2], solRight[3]);
+ printf("RN %g %g\n",nx,ny);
+ printf("RFLUX %g %g %g %g\n",FLUX[0],FLUX[1],FLUX[2],FLUX[3]);
+ */
+ _value(i,0) = FLUX[0];
+ _value(i,1) = FLUX[1];
+ _value(i,2) = FLUX[2];
+ _value(i,3) = FLUX[3];
+ _value(i,4) = -_value(i,0);
+ _value(i,5) = -_value(i,1);
+ _value(i,6) = -_value(i,2);
+ _value(i,7) = -_value(i,3);
+ }
+ }
+ };
+ public:
+ dataCacheDouble *newConvectiveFlux( dataCacheMap &cacheMap) const {
+ return new advection(cacheMap);
+ }
+ dataCacheDouble *newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const {
+ return new riemann(cacheMapLeft, cacheMapRight);
+ }
+ dataCacheDouble *newDiffusiveFlux( dataCacheMap &cacheMap) const {
+ return 0;
+ }
+ dataCacheDouble *newSourceTerm (dataCacheMap &cacheMap) const {
+ return 0;
+ }
+ dgPerfectGasLaw2d()
+ {
+ _nbf = 4; // H U(=Hu) V(=Hv)
+ }
+};
+
+class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition {
+ class term : public dataCacheDouble {
+ dataCacheDouble &sol,&normals;
+ public:
+ term(dataCacheMap &cacheMap):
+ sol(cacheMap.get("Solution",this)),
+ normals(cacheMap.get("Normals",this)){}
+ void _eval () {
+ int nQP = sol().size1();
+ if(_value.size1() != nQP)
+ _value = fullMatrix<double>(nQP,4);
+
+ for(int i=0; i< nQP; i++) {
+ const double nx = normals(0,i);
+ const double ny = normals(1,i);
+ const double solLeft [4] = {sol(i,0),sol(i,1),sol(i,2),sol(i,3)};
+ const double vn = (solLeft [1] * nx + solLeft [2] * ny);
+ const double solRight[4] = {sol(i,0),
+ sol(i,1) - 2 * vn * nx,
+ sol(i,2) - 2 * vn * ny,
+ sol(i,3)};
+ // double FLUX[4] ;
+ // _ROE2D (GAMMA,nx,ny,solLeft,solRight,FLUX);
+ const double q11 = sol(i,1)*sol(i,1)/sol(i,0);
+ const double q22 = sol(i,2)*sol(i,2)/sol(i,0);
+ const double p = (GAMMA-1)*sol(i,3) - 0.5*(GAMMA-1)*(q11+q22);
+ _value(i,0) = 0;//FLUX[0];
+ _value(i,1) = -p*nx;//FLUX[1];
+ _value(i,2) = -p*ny;//FLUX[2];
+ _value(i,3) = 0.0;//FLUX[3];
+ }
+ }
+ };
+ public:
+ dgBoundaryConditionPerfectGasLaw2dWall(){}
+ dataCacheDouble *newBoundaryTerm(dataCacheMap &cacheMapLeft) const {
+ return new term(cacheMapLeft);
+ }
+};
+
+class dgBoundaryConditionPerfectGasLaw2dFreeStream : public dgBoundaryCondition {
+ class term : public dataCacheDouble {
+ dataCacheDouble &sol,&normals,&freeStream;
+ public:
+ term(dataCacheMap &cacheMap, const std::string & freeStreamName):
+ sol(cacheMap.get("Solution",this)),
+ normals(cacheMap.get("Normals",this)),
+ freeStream(cacheMap.get(freeStreamName,this)){}
+ void _eval () {
+ int nQP = sol().size1();
+ if(_value.size1() != nQP)
+ _value = fullMatrix<double>(nQP,4);
+
+ for(int i=0; i< nQP; i++) {
+ const double nx = normals(0,i);
+ const double ny = normals(1,i);
+ const double solLeft [4] = {sol(i,0),sol(i,1),sol(i,2),sol(i,3)};
+ const double solRight[4] = {freeStream(i,0),
+ freeStream(i,1),
+ freeStream(i,2),
+ freeStream(i,3)};
+ double FLUX[4] ;
+ _ROE2D (GAMMA,nx,ny,solLeft,solRight,FLUX);
+ /*
+ printf("SOLL %g %g %g %g\n",solLeft[0],solLeft[1],solLeft[2], solLeft[3]);
+ printf("SOLR %g %g %g %g\n",solRight[0],solRight[1],solRight[2], solRight[3]);
+ printf("N %g %g\n",nx,ny);
+ printf("FLUX %g %g %g %g\n",FLUX[0],FLUX[1],FLUX[2],FLUX[3]);
+ */
+ _value(i,0) = FLUX[0];
+ _value(i,1) = FLUX[1];
+ _value(i,2) = FLUX[2];
+ _value(i,3) = FLUX[3];
+ }
+ }
+ };
+ public:
+ std::string _freeStreamName;
+ dgBoundaryConditionPerfectGasLaw2dFreeStream(std::string & freeStreamName)
+ : _freeStreamName(freeStreamName){}
+ dataCacheDouble *newBoundaryTerm(dataCacheMap &cacheMapLeft) const {
+ return new term(cacheMapLeft,_freeStreamName);
+ }
+};
+
+
+dgBoundaryCondition *dgNewBoundaryConditionPerfectGasLaw2dWall() {
+ return new dgBoundaryConditionPerfectGasLaw2dWall();
+}
+
+dgBoundaryCondition *dgNewBoundaryConditionPerfectGasLaw2dFreeStream(std::string &freeStreamName) {
+ return new dgBoundaryConditionPerfectGasLaw2dFreeStream(freeStreamName);
+}
+
+dgConservationLaw *dgNewPerfectGasLaw2d() {
+ return new dgPerfectGasLaw2d();
+}
diff --git a/Solver/dgSystemOfEquations.cpp b/Solver/dgSystemOfEquations.cpp
index 7f038fa95ebc87241b542ccbaaae16fb450f9e6f..5609f35610f5fb867b70e1497c8644fbfa6a4e97 100644
--- a/Solver/dgSystemOfEquations.cpp
+++ b/Solver/dgSystemOfEquations.cpp
@@ -66,6 +66,8 @@ int dgSystemOfEquations::setConservationLaw(lua_State *L){
_claw = dgNewConservationLawWaveEquation();
else if (_cLawName == "ShallowWater2d")
_claw = dgNewConservationLawShallowWater2d();
+ else if (_cLawName == "PerfectGas2d")
+ _claw = dgNewPerfectGasLaw2d();
else if (_cLawName == "AdvectionDiffusion"){
std::string advFunction = luaL_checkstring(L,2);
_claw = dgNewConservationLawAdvection(advFunction);
@@ -103,6 +105,17 @@ int dgSystemOfEquations::addBoundaryCondition (lua_State *L){
}
else throw;
}
+ else if (_cLawName == "PerfectGas2d"){
+ if (bcName == "Wall"){
+ _claw->addBoundaryCondition(physicalName,dgNewBoundaryConditionPerfectGasLaw2dWall());
+ }
+ else if (bcName == "FreeStream"){
+ std::string freeStreamName(luaL_checkstring(L, 3));
+ _claw->addBoundaryCondition(physicalName,
+ dgNewBoundaryConditionPerfectGasLaw2dFreeStream(freeStreamName));
+ }
+ else throw;
+ }
else throw;
}