Skip to content
Snippets Groups Projects
Commit ac707cca authored by Jean-François Remacle's avatar Jean-François Remacle
Browse files

No commit message

No commit message
parent fda545c6
No related branches found
No related tags found
No related merge requests found
......@@ -18,6 +18,7 @@ set(SRC
dgSystemOfEquations.cpp
dgConservationLawShallowWater2d.cpp
dgConservationLawWaveEquation.cpp
dgConservationLawPerfectGas.cpp
function.cpp
functionDerivator.cpp
luaFunction.cpp
......
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 ***'
......@@ -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
......
......@@ -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};
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};
......@@ -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");
}
......@@ -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
#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 &sol;
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();
}
......@@ -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;
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment