No commit message

No commit message

Solver/TESTCASES/Stommel.lua

 
+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

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)

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&);

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
   }
 };
 

Solver/dgConservationLawWaveEquation.cpp

@@ -5,53 +5,51 @@
 // 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 () {
       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);
+	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);
 
@@ -62,21 +60,19 @@ class dgConservationLawWaveEquation : public dgConservationLaw {
         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 {
+    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);
+	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);
 }

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) {

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;
   }