diff --git a/Solver/TESTCASES/Stommel.lua b/Solver/TESTCASES/Stommel.lua
index e420a01be00b27003d001b9b569c2fd2bd6a281d..bc3c89380dabf3c86b86db405f0d251dae8d9d13 100644
--- a/Solver/TESTCASES/Stommel.lua
+++ b/Solver/TESTCASES/Stommel.lua
@@ -1,8 +1,9 @@
+order = 3
model = GModel()
model:load ('stommel_square.msh')
dg = dgSystemOfEquations (model)
-dg:setOrder (1)
+dg:setOrder (order)
dg:setConservationLaw('ShallowWater2d')
dg:addBoundaryCondition('Wall','Wall')
dg:setup()
@@ -10,9 +11,11 @@ dg:setup()
dg:exportSolution('output/solution_00000')
for i=1,1000 do
- norm = dg:RK44(50)
- if ( i%100 ==0 ) then
+ norm = dg:RK44(80*(3/(2.*order+1)))
+ if ( i%10 ==0 ) then
print ('iter ', i, norm)
+ end
+ if ( i%100 ==0 ) then
dg:exportSolution(string.format('output/solution-%05d',i))
end
end
diff --git a/Solver/TESTCASES/WavePulse.lua b/Solver/TESTCASES/WavePulse.lua
index 3d33870d1f04b3cde3d4ca3ef53945793572daee..0d5c1e189a8978cfddfb3439c4273b7f55021c6f 100644
--- a/Solver/TESTCASES/WavePulse.lua
+++ b/Solver/TESTCASES/WavePulse.lua
@@ -19,7 +19,7 @@ end
Example of a lua program driving the DG code
--]]
-order = 1
+order = 5
print'*** Loading the mesh and the model ***'
myModel = GModel ()
myModel:load ('square.geo')
@@ -43,7 +43,7 @@ DG:exportSolution('output/solution_000')
print'*** solve ***'
-dt = 0.005;
+dt = 0.00125;
for i=1,1000 do
norm = DG:RK44(dt)
print('*** ITER ***',i,norm)
diff --git a/Solver/dgConservationLaw.h b/Solver/dgConservationLaw.h
index dbd7cdf68b683a7f6d1ff2b49d99384908cadcc1..f8c9f840ecfc6842f2fa95211b20dab1dfe46696 100644
--- a/Solver/dgConservationLaw.h
+++ b/Solver/dgConservationLaw.h
@@ -55,10 +55,10 @@ class dgConservationLaw {
dgConservationLaw *dgNewConservationLawAdvection(const std::string vname,const std::string nuname);
dgConservationLaw *dgNewConservationLawShallowWater2d();
-dgConservationLaw *dgNewConservationLawWaveEquation();
+dgConservationLaw *dgNewConservationLawWaveEquation(int);
dgConservationLaw *dgNewPerfectGasLaw2d();
-dgBoundaryCondition *dgNewBoundaryConditionWaveEquationWall();
+dgBoundaryCondition *dgNewBoundaryConditionWaveEquationWall(int);
dgBoundaryCondition *dgNewBoundaryConditionShallowWater2dWall();
dgBoundaryCondition *dgNewBoundaryConditionPerfectGasLaw2dWall();
dgBoundaryCondition *dgNewBoundaryConditionPerfectGasLaw2dFreeStream(std::string&);
diff --git a/Solver/dgConservationLawPerfectGas.cpp b/Solver/dgConservationLawPerfectGas.cpp
index eb0dae9ed034c3ebb5fde2e985b083d44027ac98..e90235a0d4da4bdd6e2b250e58f6ad6571756d1c 100644
--- a/Solver/dgConservationLawPerfectGas.cpp
+++ b/Solver/dgConservationLawPerfectGas.cpp
@@ -225,7 +225,7 @@ class dgPerfectGasLaw2d : public dgConservationLaw {
}
dgPerfectGasLaw2d()
{
- _nbf = 4; // H U(=Hu) V(=Hv)
+ _nbf = 4; // \rho \rho u \rho v \rho e
}
};
diff --git a/Solver/dgConservationLawWaveEquation.cpp b/Solver/dgConservationLawWaveEquation.cpp
index b5d2e457fb56c5a4d3729c536d27678392de23a7..ab6447ed281d3b2b34a67fc283610e1bfb76dddd 100644
--- a/Solver/dgConservationLawWaveEquation.cpp
+++ b/Solver/dgConservationLawWaveEquation.cpp
@@ -5,78 +5,74 @@
// dv/dt + 1/rho dp/dy = 0
static double c=1;
static double rho=1;
+
class dgConservationLawWaveEquation : public dgConservationLaw {
+ int _DIM;
class hyperbolicFlux : public dataCacheDouble {
dataCacheDouble /
+ const int _DIM,_nbf;
public:
- hyperbolicFlux(dataCacheMap &cacheMap):
- sol(cacheMap.get("Solution",this))
+ hyperbolicFlux(dataCacheMap &cacheMap,int DIM):
+ sol(cacheMap.get("Solution",this)),_DIM(DIM),_nbf(DIM+1)
{};
- void _eval () {
+ void _eval () {
int nQP = sol().size1();
if(_value.size1() != nQP)
- _value=fullMatrix<double>(nQP,9);
+ _value=fullMatrix<double>(nQP,3*_nbf);
+ _value.scale(0.);
for(int i=0; i< nQP; i++) {
- double p = sol(i,0);
- double u = sol(i,1);
- double v = sol(i,2);
- // flux_x
- _value(i,0) = c*c*rho*u;
- _value(i,1) = p/rho;
- _value(i,2) = 0;
- // flux_y
- _value(i,3) = c*c*rho*v;
- _value(i,4) = 0;
- _value(i,5) = p/rho;
- // flux_z
- _value(i,6) = 0;
- _value(i,7) = 0;
- _value(i,8) = 0;
+ const double p = sol(i,0);
+ for (int j=0;j<_DIM;j++){
+ _value(i,0 +j*_nbf) = c*c*rho*sol(i,j+1);
+ _value(i,j+1+j*_nbf) = p/rho;
+ }
}
}
};
class riemann : public dataCacheDouble {
dataCacheDouble &normals, &solL, &solR;
+ const int _DIM,_nbf;
public:
- riemann(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight):
+ riemann(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight, int DIM):
normals(cacheMapLeft.get("Normals", this)),
solL(cacheMapLeft.get("Solution", this)),
- solR(cacheMapRight.get("Solution", this))
+ solR(cacheMapRight.get("Solution", this)),
+ _DIM(DIM),_nbf(DIM+1)
{};
void _eval () {
int nQP = solL().size1();
if(_value.size1() != nQP)
- _value = fullMatrix<double>(nQP,6);
+ _value = fullMatrix<double>(nQP,2*_nbf);
for(int i=0; i< nQP; i++) {
- double nx = normals(0,i);
- double ny = normals(1,i);
- double unL = nx*solL(i,1) + ny*solL(i,2);
- double unR = nx*solR(i,1) + ny*solR(i,2);
- double pR = solR(i,0);
+ const double n[3] = {normals(0,i),normals(1,i),normals(2,i)};
+ double unL=0,unR=0;
+ for (int j=0;j<_DIM;j++){
+ unL += n[j]*solL(i,j+1);
+ unR += n[j]*solR(i,j+1);
+ }
+ double pR = solR(i,0);
double pL = solL(i,0);
- double pRiemann = 0.5 * (pL+pR + (unL-unR)*(rho*c));
+ double pRiemann = 0.5 * (pL+pR + (unL-unR)*(rho*c));
double unRiemann = 0.5 * (unL+unR + (pL-pR)/(rho*c));
double Fp = -rho*c*c*unRiemann;
double Fun = -pRiemann/rho;
_value(i,0) = Fp;
- _value(i,1) = Fun*nx;
- _value(i,2) = Fun*ny;
-
- _value(i,3) = -_value(i,0);
- _value(i,4) = -_value(i,1);
- _value(i,5) = -_value(i,2);
+ for (int j=0;j<_DIM;j++)
+ _value(i,j+1) = Fun*n[j];
+ for (int j=0;j<_nbf;j++)
+ _value(i,j+_nbf) = -_value(i,j);
}
}
};
public:
dataCacheDouble *newConvectiveFlux( dataCacheMap &cacheMap) const {
- return new hyperbolicFlux(cacheMap);
+ return new hyperbolicFlux(cacheMap,_DIM);
}
dataCacheDouble *newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const {
- return new riemann(cacheMapLeft, cacheMapRight);
+ return new riemann(cacheMapLeft, cacheMapRight,_DIM);
}
dataCacheDouble *newDiffusiveFlux( dataCacheMap &cacheMap) const {
return 0;
@@ -84,44 +80,45 @@ class dgConservationLawWaveEquation : public dgConservationLaw {
dataCacheDouble *newSourceTerm (dataCacheMap &cacheMap) const {
return 0;
}
- dgConservationLawWaveEquation()
+ dgConservationLawWaveEquation(int DIM) : _DIM(DIM)
{
- _nbf = 3; // H U(=Hu) V(=Hv)
+ _nbf = 1 + _DIM; // H U(=Hu) V(=Hv)
}
};
class dgBoundaryConditionWaveEquationWall : public dgBoundaryCondition {
+ int _DIM;
class term : public dataCacheDouble {
- dataCacheDouble &sol,&normals;
+ int _DIM;
+ dataCacheDouble &sol,&normals;
public:
- term(dataCacheMap &cacheMap):
+ term(dataCacheMap &cacheMap, int DIM):
sol(cacheMap.get("Solution",this)),
- normals(cacheMap.get("Normals",this)){}
+ normals(cacheMap.get("Normals",this)),
+ _DIM(DIM){}
void _eval () {
int nQP = sol().size1();
if(_value.size1() != nQP)
- _value = fullMatrix<double>(nQP,3);
+ _value = fullMatrix<double>(nQP,_DIM+1);
for(int i=0; i< nQP; i++) {
- double nx = normals(0,i);
- double ny = normals(1,i);
- double p = sol(i,0);
-
+ const double n[3] = {normals(0,i),normals(1,i),normals(2,i)};
+ double p = sol(i,0);
_value(i,0) = 0;
- _value(i,1) = -p/rho*nx;
- _value(i,2) = -p/rho*ny;
+ for (int j=0;j<_DIM;j++)
+ _value(i,j+1) = -p/rho*n[j];
}
}
};
public:
- dgBoundaryConditionWaveEquationWall(){}
+ dgBoundaryConditionWaveEquationWall(int DIM):_DIM(DIM){}
dataCacheDouble *newBoundaryTerm(dataCacheMap &cacheMapLeft) const {
- return new term(cacheMapLeft);
+ return new term(cacheMapLeft,_DIM);
}
};
-dgConservationLaw *dgNewConservationLawWaveEquation() {
- return new dgConservationLawWaveEquation();
+dgConservationLaw *dgNewConservationLawWaveEquation(int DIM) {
+ return new dgConservationLawWaveEquation(DIM);
}
-dgBoundaryCondition *dgNewBoundaryConditionWaveEquationWall() {
- return new dgBoundaryConditionWaveEquationWall();
+dgBoundaryCondition *dgNewBoundaryConditionWaveEquationWall(int DIM) {
+ return new dgBoundaryConditionWaveEquationWall(DIM);
}
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index 62778f75d2fdab635841205bf3232d05e9290836..b2c8a303d62dfa9b7eed5dae9a2a650a4f50486a 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -150,6 +150,7 @@ void dgGroupOfFaces::addFace(const MFace &topoFace, int iElLeft, int iElRight){
for (int j=0; j<geomClosure.size() ; j++)
vertices.push_back( elLeft.getVertex(geomClosure[j]) );
// triangular face
+ printf("the topological face has %d vertices and the geometrical %d\n",topoFace.getNumVertices(),vertices.size());
if (topoFace.getNumVertices() == 3){
switch(vertices.size()){
case 3 : _faces.push_back(new MTriangle (vertices) ); break;
@@ -161,7 +162,9 @@ void dgGroupOfFaces::addFace(const MFace &topoFace, int iElLeft, int iElRight){
}
}
// quad face 2 do
- else throw;
+ else {
+ throw;
+ }
}
void dgGroupOfFaces::addEdge(const MEdge &topoEdge, int iElLeft, int iElRight) {
diff --git a/Solver/dgSystemOfEquations.cpp b/Solver/dgSystemOfEquations.cpp
index 12add17a2b79cd69dcbfdd1945570a14b547fcce..6372d32deae57b24eb77b89ebab6d36459d0a4c8 100644
--- a/Solver/dgSystemOfEquations.cpp
+++ b/Solver/dgSystemOfEquations.cpp
@@ -63,7 +63,7 @@ int dgSystemOfEquations::setConservationLaw(lua_State *L){
//int argc = (int)luaL_checknumber(L,0);
_cLawName = std::string (luaL_checkstring(L, 1));
if (_cLawName == "WaveEquation")
- _claw = dgNewConservationLawWaveEquation();
+ _claw = dgNewConservationLawWaveEquation(_dimension);
else if (_cLawName == "ShallowWater2d")
_claw = dgNewConservationLawShallowWater2d();
else if (_cLawName == "PerfectGas2d")
@@ -94,7 +94,7 @@ int dgSystemOfEquations::addBoundaryCondition (lua_State *L){
//specific boundary conditions
else if (_cLawName == "WaveEquation"){
if (bcName == "Wall"){
- _claw->addBoundaryCondition(physicalName,dgNewBoundaryConditionWaveEquationWall());
+ _claw->addBoundaryCondition(physicalName,dgNewBoundaryConditionWaveEquationWall(_dimension));
}
else throw;
}