No commit message

No commit message

Common/LuaBindings.cpp

@@ -15,6 +15,7 @@
 #include "dgGroupOfElements.h"
 #include "dgDofContainer.h"
 #include "dgConservationLawShallowWater2d.h"
+#include "dgConservationLawShallowWater1d.h"
 #include "dgConservationLawAdvection.h"
 #include "dgConservationLawPerfectGas.h"
 #include "dgConservationLawWaveEquation.h"
@@ -344,6 +345,7 @@ binding::binding(){
   dgConservationLaw::registerBindings(this);
   dgConservationLawAdvectionDiffusionRegisterBindings(this);
    dgConservationLawMaxwellRegisterBindings(this);
+   dgConservationLawShallowWater1dRegisterBindings(this);
   dgConservationLawShallowWater2dRegisterBindings(this);
   dgConservationLawWaveEquationRegisterBindings(this);
   dgDofContainer::registerBindings(this);

Geo/GFaceCompound.cpp

@@ -503,15 +503,14 @@ bool GFaceCompound::parametrize() const
   else if (_mapping == CONFORMAL){
     Msg::Debug("Parametrizing surface %d with 'conformal map'", tag());
     fillNeumannBCS();
-    bool withoutFolding = parametrize_conformal_spectral() ;
-    printStuff();
-    //exit(1);
+    //bool withoutFolding = parametrize_conformal_spectral() ;
+    bool withoutFolding = parametrize_conformal();
     if ( withoutFolding == false ){
+      //printStuff(); exit(1);
       Msg::Warning("$$$ Parametrization switched to harmonic map");
       parametrize(ITERU,HARMONIC); 
       parametrize(ITERV,HARMONIC);
     }
-    //exit(1);
   }
   //Distance function
   //-----------------
@@ -1138,7 +1137,7 @@ bool GFaceCompound::parametrize_conformal_spectral() const
       cross21.addToMatrix(myAssembler, &se);
     }
   }
-  double epsilon = 1.e-6;
+  double epsilon = 1.e-7;
   for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
       MVertex *v = *itv;
       if (std::find(ordered.begin(), ordered.end(), v) == ordered.end() ){
@@ -1167,11 +1166,16 @@ bool GFaceCompound::parametrize_conformal_spectral() const
        myAssembler.assemble(v, 0, 2, v, 0, 2,  1.0);
      }
    }
+
 //    int NB = ordered.size();
-//    for(std::vector<MVertex *>::iterator itv1 = ordered.begin(); itv1 !=ordered.end() ; ++itv1){
-//      for(std::vector<MVertex *>::iterator itv2 = ordered.begin(); itv2 !=ordered.end() ; ++itv2){
-//        myAssembler.assemble(*itv1, 0, 1, *itv2, 0, 1,  -1/NB);
-//        myAssembler.assemble(*itv1, 0, 2, *itv2, 0, 2,  -1/NB);
+//    for (int i = 0; i< NB; i++){
+//      MVertex *v1 = ordered[i];
+//      for (int j = i; j< NB; j++){
+//        MVertex *v2 = ordered[j];
+//        myAssembler.assemble(v1, 0, 1, v2, 0, 1,  -1./NB);
+//        myAssembler.assemble(v1, 0, 2, v2, 0, 2,  -1./NB);
+//        myAssembler.assemble(v2, 0, 1, v1, 0, 1,  -1./NB);
+//        myAssembler.assemble(v2, 0, 2, v1, 0, 2,  -1./NB);
 //      }
 //    }
  
@@ -1188,7 +1192,7 @@ bool GFaceCompound::parametrize_conformal_spectral() const
 
    //-------------------------------
    eigenSolver eig(&myAssembler, "B" , "A", true);
-   bool converged = eig.solve(2, "largest");
+   bool converged = eig.solve(1, "largest");
      
    if(converged) {
      int k = 0;
@@ -1202,19 +1206,19 @@ bool GFaceCompound::parametrize_conformal_spectral() const
      }
      
      //if folding take second sallest eigenvalue
-     bool noFolding = checkFolding(ordered);
-     if (!noFolding ){
-       coordinates.clear();
-       int k = 0;
-       std::vector<std::complex<double> > &ev = eig.getEigenVector(1); 
-       for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
-	 MVertex *v = *itv;
-	 double paramu = ev[k].real();
-	 double paramv = ev[k+1].real();
-	 coordinates[v] = SPoint3(paramu,paramv,0.0);
-	 k = k+2;
-       }
-     }
+//      bool noFolding = checkFolding(ordered);
+//      if (!noFolding ){
+//        coordinates.clear();
+//        int k = 0;
+//        std::vector<std::complex<double> > &ev = eig.getEigenVector(1); 
+//        for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
+// 	 MVertex *v = *itv;
+// 	 double paramu = ev[k].real();
+// 	 double paramv = ev[k+1].real();
+// 	 coordinates[v] = SPoint3(paramu,paramv,0.0);
+// 	 k = k+2;
+//        }
+//      }
 
      lsysA->clear();
      lsysB->clear();

Geo/MLine.cpp

@@ -39,3 +39,13 @@ void MLine::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
   *npts = nbP;
   *pts = GQL;
 }
+
+double MLine::getInnerRadius()
+{
+#if defined(HAVE_MESH)
+  double dist = _v[0]->distance(_v[1]);
+  return dist;
+#else
+  return 0.;
+#endif
+}

Geo/MLine.h

@@ -37,6 +37,7 @@ class MLine : public MElement {
   virtual int getDim() const { return 1; }
   virtual int getNumVertices() const { return 2; }
   virtual MVertex *getVertex(int num){ return _v[num]; }
+  virtual double getInnerRadius(); // length of segment line
   virtual void getVertexInfo(const MVertex * vertex, int &ithVertex) const
   { 
     ithVertex = _v[0] == vertex ? 0 : 1;

Mesh/multiscalePartition.cpp

@@ -167,51 +167,52 @@ static int getAspectRatio(std::vector<MElement *> &elements,
     AR = (int) ceil(2*3.14*area3D/(tot_length*tot_length));
   }
 
-//   std::set<MVertex*> vs;
-//   for(unsigned int i = 0; i < elements.size(); i++){
-//     MElement *e = elements[i];
-//     for(int j = 0; j < e->getNumVertices(); j++){
-//       vs.insert(e->getVertex(j));
-//     }
-//   }
-//   SBoundingBox3d bb;
-//   std::vector<SPoint3> vertices;
-//   for (std::set<MVertex* >::iterator it = vs.begin(); it != vs.end(); it++){
-//     SPoint3 pt((*it)->x(),(*it)->y(), (*it)->z());
-//     vertices.push_back(pt);
-//     bb += pt;
-//   }
-//   double H = norm(SVector3(bb.max(), bb.min()));
-
-//   //SOrientedBoundingBox obbox =  SOrientedBoundingBox::buildOBB(vertices);
-//   //double H = obbox.getMaxSize(); 
-  
-//   double D = H;
-//   if (boundaries.size()  > 0 ) D = 0.0;
-//   for (unsigned int i = 0; i < boundaries.size(); i++){
-//     std::set<MVertex*> vb;
-//     std::vector<MEdge> iBound = boundaries[i];
-//     for (unsigned int j = 0; j < iBound.size(); j++){
-//       MEdge e = iBound[j];
-//       vb.insert(e.getVertex(0));
-//       vb.insert(e.getVertex(1));
-//     }
-//     std::vector<SPoint3> vBounds;
-//     SBoundingBox3d bb;
-//     for (std::set<MVertex* >::iterator it = vb.begin(); it != vb.end(); it++){
-//       SPoint3 pt((*it)->x(),(*it)->y(), (*it)->z());
-//       vBounds.push_back(pt);
-//       bb +=pt;
-//     }
-//     double iD = norm(SVector3(bb.max(), bb.min()));
-//     D = std::max(D, iD);
-    
-//     //SOrientedBoundingBox obboxD = SOrientedBoundingBox::buildOBB(vBounds); 
-//     //D = std::max(D, obboxD.getMaxSize());
-//   }
-//   int AR = (int)ceil(H/D);
-  
-  return AR;
+  //compute AR also with Bounding box
+  std::set<MVertex*> vs;
+  for(unsigned int i = 0; i < elements.size(); i++){
+    MElement *e = elements[i];
+    for(int j = 0; j < e->getNumVertices(); j++){
+      vs.insert(e->getVertex(j));
+    }
+  }
+  SBoundingBox3d bb;
+  std::vector<SPoint3> vertices;
+  for (std::set<MVertex* >::iterator it = vs.begin(); it != vs.end(); it++){
+    SPoint3 pt((*it)->x(),(*it)->y(), (*it)->z());
+    vertices.push_back(pt);
+    bb += pt;
+  }
+  double H = norm(SVector3(bb.max(), bb.min()));
+
+  //SOrientedBoundingBox obbox =  SOrientedBoundingBox::buildOBB(vertices);
+  //double H = obbox.getMaxSize(); 
+  
+  double D = H;
+  if (boundaries.size()  > 0 ) D = 0.0;
+  for (unsigned int i = 0; i < boundaries.size(); i++){
+    std::set<MVertex*> vb;
+    std::vector<MEdge> iBound = boundaries[i];
+    for (unsigned int j = 0; j < iBound.size(); j++){
+      MEdge e = iBound[j];
+      vb.insert(e.getVertex(0));
+      vb.insert(e.getVertex(1));
+    }
+    std::vector<SPoint3> vBounds;
+    SBoundingBox3d bb;
+    for (std::set<MVertex* >::iterator it = vb.begin(); it != vb.end(); it++){
+      SPoint3 pt((*it)->x(),(*it)->y(), (*it)->z());
+      vBounds.push_back(pt);
+      bb +=pt;
+    }
+    double iD = norm(SVector3(bb.max(), bb.min()));
+    D = std::max(D, iD);
+    
+    //SOrientedBoundingBox obboxD = SOrientedBoundingBox::buildOBB(vBounds); 
+    //D = std::max(D, obboxD.getMaxSize());
+  }
+  int AR2 = (int)ceil(H/D);
+  
+  return std::min(AR, AR2);
 }
 
 static void getGenusAndRatio(std::vector<MElement *> &elements, int & genus, int &AR, int &NB)

Solver/CMakeLists.txt

@@ -18,6 +18,7 @@ set(SRC
   dgConservationLaw.cpp
   dgConservationLawAdvection.cpp
   dgSystemOfEquations.cpp
+  dgConservationLawShallowWater1d.cpp	
   dgConservationLawShallowWater2d.cpp
   dgConservationLawWaveEquation.cpp
   dgConservationLawMaxwell.cpp

Solver/TESTCASES/Advection1D.lua

@@ -63,12 +63,13 @@ solution:exportMsh('output/init_limit')
 
 print'*** solve ***'
 local x = os.clock()
-n = 5
 dt = 0.03
-for i=1,100*n do
+for i=1,100 do
   norm = rk:iterate44(law,dt,solution)
-  if (i % n == 0) then 
+   if (i % 20 == 0) then 
        print('|ITER|',i,'|NORM|',norm,'|DT|',dt,'|CPU|',os.clock() - x)
+    end
+    if (i % 200 == 0) then 
      solution:exportMsh(string.format('output/solution-%06d',i))
     end
 end

Solver/TESTCASES/ForwardFacingStep.lua

@@ -41,7 +41,7 @@ FS = functionLua(4, 'free_stream', {'XYZ'}):getName()
 law=dgPerfectGasLaw2d()
 DG:setConservationLaw(law)
 
-law:addBoundaryCondition('Walls',law:newBoundaryWall())
+law:addBoundaryCondition('Walls',law:newNonSlipWallBoundary())
 law:addBoundaryCondition('LeftRight',law:newOutsideValueBoundary(FS))
 --law:addBoundaryCondition('Walls',law:newOutsideValueBoundary(FS))
 

Solver/TESTCASES/ForwardFacingStepNEW.lua

+MACH = 3.0;
+GAMMA = 1.4;
+U = 3.0
+V = 0.0 
+RHO  = 1.4;
+
+PRES = RHO*U*U/(GAMMA*MACH*MACH)
+
+--PRES = 1;
+--PRES = ./(MACH*RHO*RHO*GAMMA*GAMMA) 
+
+SOUND = math.sqrt(U*U+V*V)/MACH
+
+--[[ 
+     Function for initial conditions
+--]]
+function free_stream( XYZ, FCT )
+  for i=0,XYZ:size1()-1 do
+    FCT:set(i,0,RHO) 
+    FCT:set(i,1,RHO*U) 
+    FCT:set(i,2,RHO*V) 
+    FCT:set(i,3, 0.5*RHO*(U*U+V*V)+PRES/(GAMMA-1)) 
+  end
+end
+
+--[[ 
+     Example of a lua program driving the DG code
+--]]
+
+
+print'*** Loading the mesh and the model ***'
+model   = GModel  ()
+model:load ('step.geo')
+model:load ('step.msh')
+order = 1
+dimension = 2
+
+FS = functionLua(4, 'free_stream', {'XYZ'}):getName()
+
+-- boundary condition
+law=dgPerfectGasLaw2d()
+law:addBoundaryCondition('Walls',law:newNonSlipWallBoundary())
+law:addBoundaryCondition('LeftRight',law:newOutsideValueBoundary(FS))
+--law:addBoundaryCondition('Walls',law:newOutsideValueBoundary(FS))
+
+groups = dgGroupCollection(model, dimension, order)
+groups:buildGroupsOfInterfaces()
+
+-- build Runge Kutta and limiter
+rk=dgRungeKutta()
+limiter = dgSlopeLimiter(law)
+rk:setLimiter(limiter) 
+
+-- build solution vector
+solution = dgDofContainer(groups, law:getNbFields())
+solution:L2Projection(FS)
+solution:exportGroupIdMsh()
+
+print'*** setting the initial solution ***'
+solution:exportMsh('output/init')
+limiter:apply(solution)
+solution:exportMsh('output/init_limit')
+
+print'*** solve ***'
+dg = dgSystemOfEquations (model)
+CFL = 2
+local x = os.clock()
+for i=1,5000 do
+    dt = CFL * rk:computeInvSpectralRadius(law,solution);    
+--  norm = rk:iterate44(law,dt,solution)
+    norm = rk:iterateEuler(law,dt,solution)
+    if (i % 10 == 0) then 
+       print('|ITER|',i,'|NORM|',norm,'|DT|',dt,'|CPU|',os.clock() - x)
+    end
+    if (i % 100 == 0) then 
+       solution:exportMsh(string.format('output/solution-%06d', i)) 
+    end
+end
+
+print'*** done ***'
+
+

Solver/TESTCASES/MultirateAdvection.lua

@@ -60,7 +60,7 @@ FS = functionLua(1, 'initial_condition', {'XYZ'}):getName()
 GC=dgGroupCollection(myModel,2,order)
 solTmp=dgDofContainer(GC,1)
 solTmp:L2Projection(FS)
-dt=GC:splitGroupsForMultirate(100,law,solTmp)
+dt=GC:splitGroupsForMultirate(4,law,solTmp)
 GC:buildGroupsOfInterfaces(myModel,2,order)
 solution=dgDofContainer(GC,1)
 solution:L2Projection(FS)

Solver/TESTCASES/edge.geo

+
 Point(1) = {-1, 0, 0, 0.01};
 Point(2) = {1, 0, 0, 0.01};
 Line(1) = {1, 2};
+
 Physical Point("Left") = {1};
 Physical Point("Right") = {2};
 Physical Line("Line") = {1};

Solver/TESTCASES/rect.geo

@@ -15,9 +15,11 @@ Line Loop(8) = {7, -2, 5, 6};
 Plane Surface(9) = {8};
 Line Loop(10) = {3, 4, 1, 2};
 Plane Surface(11) = {10};
+
 Physical Surface("sprut") = {11, 9};
 Physical Line("Walls") = {5, 7, 3, 4, 1};
 Physical Line("Top") = {6};
+
 Recombine Surface {9, 11};
 
 Field[1] = MathEval;

Solver/TESTCASES/validation/Advection1D.lua

+
+--[[ 
+     Function for initial conditions
+--]]
+
+function initial_condition( xyz , f )
+  for i=0,xyz:size1()-1 do
+    x = xyz:get(i,0)
+    y = xyz:get(i,1)
+    z = xyz:get(i,2)
+    if (x>-0.3 and x<0.3) then
+     f:set (i, 0, 1)
+   else
+     f:set (i, 0, 0)
+    end	
+  end
+end
+
+--[[ 
+     Example of a lua program driving the DG code
+--]]
+
+model = GModel  ()
+model:load ('edge.msh')
+order=1
+dimension=1
+
+-- conservation law
+-- advection speed
+v=fullMatrix(3,1);
+v:set(0,0,0.15)
+v:set(1,0,0)
+v:set(2,0,0)
+-- diffusivity
+nu=fullMatrix(1,1);
+nu:set(0,0,0)
+
+law = dgConservationLawAdvectionDiffusion(functionConstant(v):getName(), '') 
+
+-- boundary condition
+outside=fullMatrix(1,1)
+outside:set(0,0,0.)
+bndcondition=law:newOutsideValueBoundary(functionConstant(outside):getName())
+law:addBoundaryCondition('Left',bndcondition)
+law:addBoundaryCondition('Right',bndcondition)
+
+groups = dgGroupCollection(model, dimension, order)
+groups:buildGroupsOfInterfaces()
+
+-- build Runge Kutta and limiter
+rk=dgRungeKutta()
+limiter = dgSlopeLimiter(law)
+rk:setLimiter(limiter) 
+
+-- build solution vector
+FS = functionLua(1, 'initial_condition', {'XYZ'}):getName()
+solution = dgDofContainer(groups, law:getNbFields())
+solution:L2Projection(FS)
+
+solution:exportMsh('output/init')
+limiter:apply(solution)
+solution:exportMsh('output/init_limit')
+
+print'*** solve ***'
+local x = os.clock()
+dt = 0.03
+for i=1,100 do
+  norm = rk:iterate44(law,dt,solution)
+   if (i % 20 == 0) then 
+       print('|ITER|',i,'|NORM|',norm,'|DT|',dt,'|CPU|',os.clock() - x)
+    end
+    if (i % 200 == 0) then 
+     solution:exportMsh(string.format('output/solution-%06d',i))
+    end
+end

Solver/TESTCASES/validation/edge.msh

+$MeshFormat
+2.2 0 8
+$EndMeshFormat
+$PhysicalNames
+3
+0 1 "Left"
+0 2 "Right"
+1 3 "Line"
+$EndPhysicalNames
+$Nodes
+200
+1 -1 0 0
+2 1 0 0
+3 -0.9899497487437465 0 0
+4 -0.9798994974874929 0 0
+5 -0.9698492462312395 0 0
+6 -0.959798994974986 0 0
+7 -0.9497487437187324 0 0
+8 -0.9396984924624789 0 0
+9 -0.9296482412062255 0 0
+10 -0.9195979899499719 0 0
+11 -0.9095477386937184 0 0
+12 -0.8994974874374648 0 0
+13 -0.8894472361812114 0 0
+14 -0.8793969849249579 0 0
+15 -0.8693467336687044 0 0
+16 -0.8592964824124508 0 0
+17 -0.8492462311561974 0 0
+18 -0.8391959798999438 0 0
+19 -0.8291457286436903 0 0
+20 -0.8190954773874368 0 0
+21 -0.8090452261311832 0 0
+22 -0.7989949748749298 0 0
+23 -0.7889447236186762 0 0
+24 -0.7788944723624227 0 0
+25 -0.7688442211061692 0 0
+26 -0.7587939698499158 0 0
+27 -0.7487437185936622 0 0
+28 -0.7386934673374087 0 0
+29 -0.7286432160811551 0 0
+30 -0.7185929648249016 0 0
+31 -0.7085427135686482 0 0
+32 -0.6984924623123947 0 0
+33 -0.6884422110561411 0 0
+34 -0.6783919597998875 0 0
+35 -0.6683417085436341 0 0
+36 -0.6582914572873806 0 0
+37 -0.6482412060311271 0 0
+38 -0.6381909547748736 0 0
+39 -0.6281407035186201 0 0
+40 -0.6180904522623665 0 0
+41 -0.608040201006113 0 0
+42 -0.5979899497498595 0 0
+43 -0.5879396984936061 0 0
+44 -0.5778894472373525 0 0
+45 -0.567839195981099 0 0
+46 -0.5577889447248454 0 0
+47 -0.5477386934685919 0 0
+48 -0.5376884422123385 0 0
+49 -0.527638190956085 0 0
+50 -0.5175879396998314 0 0
+51 -0.5075376884435778 0 0
+52 -0.4974874371873209 0 0
+53 -0.4874371859310535 0 0
+54 -0.477386934674786 0 0
+55 -0.4673366834185185 0 0
+56 -0.4572864321622511 0 0
+57 -0.4472361809059836 0 0
+58 -0.4371859296497161 0 0
+59 -0.4271356783934487 0 0
+60 -0.4170854271371812 0 0
+61 -0.4070351758809138 0 0
+62 -0.3969849246246463 0 0
+63 -0.3869346733683788 0 0
+64 -0.3768844221121114 0 0
+65 -0.3668341708558439 0 0
+66 -0.3567839195995764 0 0
+67 -0.346733668343309 0 0
+68 -0.3366834170870415 0 0
+69 -0.3266331658307741 0 0
+70 -0.3165829145745066 0 0
+71 -0.3065326633182391 0 0
+72 -0.2964824120619717 0 0
+73 -0.2864321608057042 0 0
+74 -0.2763819095494368 0 0
+75 -0.2663316582931693 0 0
+76 -0.2562814070369018 0 0
+77 -0.2462311557806344 0 0
+78 -0.2361809045243669 0 0
+79 -0.2261306532680994 0 0
+80 -0.216080402011832 0 0
+81 -0.2060301507555645 0 0
+82 -0.1959798994992971 0 0
+83 -0.1859296482430296 0 0
+84 -0.1758793969867621 0 0
+85 -0.1658291457304947 0 0
+86 -0.1557788944742272 0 0
+87 -0.1457286432179598 0 0
+88 -0.1356783919616923 0 0
+89 -0.1256281407054248 0 0
+90 -0.1155778894491574 0 0
+91 -0.1055276381928899 0 0
+92 -0.09547738693662244 0 0
+93 -0.08542713568035498 0 0
+94 -0.07537688442408752 0 0
+95 -0.06532663316782006 0 0
+96 -0.0552763819115526 0 0
+97 -0.04522613065528525 0 0
+98 -0.03517587939901778 0 0
+99 -0.02512562814275032 0 0
+100 -0.01507537688648286 0 0
+101 -0.005025125630215399 0 0
+102 0.005025125626066052 0 0
+103 0.01507537688236149 0 0
+104 0.02512562813865671 0 0
+105 0.03517587939495215 0 0
+106 0.04522613065124759 0 0
+107 0.0552763819075428 0 0
+108 0.06532663316383824 0 0
+109 0.07537688442013346 0 0
+110 0.0854271356764289 0 0
+111 0.09547738693272434 0 0
+112 0.1055276381890196 0 0
+113 0.115577889445315 0 0
+114 0.1256281407016102 0 0
+115 0.1356783919579057 0 0
+116 0.1457286432142011 0 0
+117 0.1557788944704963 0 0
+118 0.1658291457267917 0 0
+119 0.1758793969830872 0 0
+120 0.1859296482393824 0 0
+121 0.1959798994956778 0 0
+122 0.2060301507519731 0 0
+123 0.2160804020082685 0 0
+124 0.2261306532645639 0 0
+125 0.2361809045208592 0 0
+126 0.2462311557771546 0 0
+127 0.25628140703345 0 0
+128 0.2663316582897453 0 0
+129 0.2763819095460407 0 0
+130 0.2864321608023359 0 0
+131 0.2964824120586314 0 0
+132 0.3065326633149268 0 0
+133 0.316582914571222 0 0
+134 0.3266331658275174 0 0
+135 0.3366834170838127 0 0
+136 0.3467336683401081 0 0
+137 0.3567839195964035 0 0
+138 0.3668341708526988 0 0
+139 0.3768844221089942 0 0
+140 0.3869346733652896 0 0
+141 0.3969849246215849 0 0
+142 0.4070351758778803 0 0
+143 0.4170854271341755 0 0
+144 0.427135678390471 0 0
+145 0.4371859296467664 0 0
+146 0.4472361809030616 0 0
+147 0.457286432159357 0 0
+148 0.4673366834156525 0 0
+149 0.4773869346719477 0 0
+150 0.4874371859282431 0 0
+151 0.4974874371845384 0 0
+152 0.5075376884408442 0 0
+153 0.5175879396971537 0 0
+154 0.5276381909534629 0 0
+155 0.5376884422097721 0 0
+156 0.5477386934660815 0 0
+157 0.5577889447223907 0 0
+158 0.5678391959787001 0 0
+159 0.5778894472350093 0 0
+160 0.5879396984913188 0 0
+161 0.597989949747628 0 0
+162 0.6080402010039372 0 0
+163 0.6180904522602466 0 0
+164 0.6281407035165558 0 0
+165 0.6381909547728652 0 0
+166 0.6482412060291745 0 0
+167 0.6582914572854837 0 0
+168 0.6683417085417931 0 0
+169 0.6783919597981023 0 0
+170 0.6884422110544117 0 0
+171 0.6984924623107209 0 0
+172 0.7085427135670304 0 0
+173 0.7185929648233396 0 0
+174 0.7286432160796488 0 0
+175 0.7386934673359582 0 0
+176 0.7487437185922674 0 0
+177 0.7587939698485768 0 0
+178 0.768844221104886 0 0
+179 0.7788944723611955 0 0
+180 0.7889447236175047 0 0
+181 0.7989949748738139 0 0
+182 0.8090452261301233 0 0
+183 0.8190954773864325 0 0
+184 0.8291457286427419 0 0
+185 0.8391959798990511 0 0
+186 0.8492462311553606 0 0
+187 0.8592964824116698 0 0
+188 0.869346733667979 0 0
+189 0.8793969849242884 0 0
+190 0.8894472361805976 0 0
+191 0.8994974874369071 0 0
+192 0.9095477386932163 0 0
+193 0.9195979899495255 0 0
+194 0.9296482412058349 0 0
+195 0.9396984924621441 0 0
+196 0.9497487437184535 0 0
+197 0.9597989949747627 0 0
+198 0.9698492462310722 0 0
+199 0.9798994974873814 0 0
+200 0.9899497487436906 0 0
+$EndNodes
+$Elements
+201
+1 15 2 1 1 1
+2 15 2 2 2 2
+3 1 2 3 1 1 3
+4 1 2 3 1 3 4
+5 1 2 3 1 4 5
+6 1 2 3 1 5 6
+7 1 2 3 1 6 7
+8 1 2 3 1 7 8
+9 1 2 3 1 8 9
+10 1 2 3 1 9 10
+11 1 2 3 1 10 11
+12 1 2 3 1 11 12
+13 1 2 3 1 12 13
+14 1 2 3 1 13 14
+15 1 2 3 1 14 15
+16 1 2 3 1 15 16
+17 1 2 3 1 16 17
+18 1 2 3 1 17 18
+19 1 2 3 1 18 19
+20 1 2 3 1 19 20
+21 1 2 3 1 20 21
+22 1 2 3 1 21 22
+23 1 2 3 1 22 23
+24 1 2 3 1 23 24
+25 1 2 3 1 24 25
+26 1 2 3 1 25 26
+27 1 2 3 1 26 27
+28 1 2 3 1 27 28
+29 1 2 3 1 28 29
+30 1 2 3 1 29 30
+31 1 2 3 1 30 31
+32 1 2 3 1 31 32
+33 1 2 3 1 32 33
+34 1 2 3 1 33 34
+35 1 2 3 1 34 35
+36 1 2 3 1 35 36
+37 1 2 3 1 36 37
+38 1 2 3 1 37 38
+39 1 2 3 1 38 39
+40 1 2 3 1 39 40
+41 1 2 3 1 40 41
+42 1 2 3 1 41 42
+43 1 2 3 1 42 43
+44 1 2 3 1 43 44
+45 1 2 3 1 44 45
+46 1 2 3 1 45 46
+47 1 2 3 1 46 47
+48 1 2 3 1 47 48
+49 1 2 3 1 48 49
+50 1 2 3 1 49 50
+51 1 2 3 1 50 51
+52 1 2 3 1 51 52
+53 1 2 3 1 52 53
+54 1 2 3 1 53 54
+55 1 2 3 1 54 55
+56 1 2 3 1 55 56
+57 1 2 3 1 56 57
+58 1 2 3 1 57 58
+59 1 2 3 1 58 59
+60 1 2 3 1 59 60
+61 1 2 3 1 60 61
+62 1 2 3 1 61 62
+63 1 2 3 1 62 63
+64 1 2 3 1 63 64
+65 1 2 3 1 64 65
+66 1 2 3 1 65 66
+67 1 2 3 1 66 67
+68 1 2 3 1 67 68
+69 1 2 3 1 68 69
+70 1 2 3 1 69 70
+71 1 2 3 1 70 71
+72 1 2 3 1 71 72
+73 1 2 3 1 72 73
+74 1 2 3 1 73 74
+75 1 2 3 1 74 75
+76 1 2 3 1 75 76
+77 1 2 3 1 76 77
+78 1 2 3 1 77 78
+79 1 2 3 1 78 79
+80 1 2 3 1 79 80
+81 1 2 3 1 80 81
+82 1 2 3 1 81 82
+83 1 2 3 1 82 83
+84 1 2 3 1 83 84
+85 1 2 3 1 84 85
+86 1 2 3 1 85 86
+87 1 2 3 1 86 87
+88 1 2 3 1 87 88
+89 1 2 3 1 88 89
+90 1 2 3 1 89 90
+91 1 2 3 1 90 91
+92 1 2 3 1 91 92
+93 1 2 3 1 92 93
+94 1 2 3 1 93 94
+95 1 2 3 1 94 95
+96 1 2 3 1 95 96
+97 1 2 3 1 96 97
+98 1 2 3 1 97 98
+99 1 2 3 1 98 99
+100 1 2 3 1 99 100
+101 1 2 3 1 100 101
+102 1 2 3 1 101 102
+103 1 2 3 1 102 103
+104 1 2 3 1 103 104
+105 1 2 3 1 104 105
+106 1 2 3 1 105 106
+107 1 2 3 1 106 107
+108 1 2 3 1 107 108
+109 1 2 3 1 108 109
+110 1 2 3 1 109 110
+111 1 2 3 1 110 111
+112 1 2 3 1 111 112
+113 1 2 3 1 112 113
+114 1 2 3 1 113 114
+115 1 2 3 1 114 115
+116 1 2 3 1 115 116
+117 1 2 3 1 116 117
+118 1 2 3 1 117 118
+119 1 2 3 1 118 119
+120 1 2 3 1 119 120
+121 1 2 3 1 120 121
+122 1 2 3 1 121 122
+123 1 2 3 1 122 123
+124 1 2 3 1 123 124
+125 1 2 3 1 124 125
+126 1 2 3 1 125 126
+127 1 2 3 1 126 127
+128 1 2 3 1 127 128
+129 1 2 3 1 128 129
+130 1 2 3 1 129 130
+131 1 2 3 1 130 131
+132 1 2 3 1 131 132
+133 1 2 3 1 132 133
+134 1 2 3 1 133 134
+135 1 2 3 1 134 135
+136 1 2 3 1 135 136
+137 1 2 3 1 136 137
+138 1 2 3 1 137 138
+139 1 2 3 1 138 139
+140 1 2 3 1 139 140
+141 1 2 3 1 140 141
+142 1 2 3 1 141 142
+143 1 2 3 1 142 143
+144 1 2 3 1 143 144
+145 1 2 3 1 144 145
+146 1 2 3 1 145 146
+147 1 2 3 1 146 147
+148 1 2 3 1 147 148
+149 1 2 3 1 148 149
+150 1 2 3 1 149 150
+151 1 2 3 1 150 151
+152 1 2 3 1 151 152
+153 1 2 3 1 152 153
+154 1 2 3 1 153 154
+155 1 2 3 1 154 155
+156 1 2 3 1 155 156
+157 1 2 3 1 156 157
+158 1 2 3 1 157 158
+159 1 2 3 1 158 159
+160 1 2 3 1 159 160
+161 1 2 3 1 160 161
+162 1 2 3 1 161 162
+163 1 2 3 1 162 163
+164 1 2 3 1 163 164
+165 1 2 3 1 164 165
+166 1 2 3 1 165 166
+167 1 2 3 1 166 167
+168 1 2 3 1 167 168
+169 1 2 3 1 168 169
+170 1 2 3 1 169 170
+171 1 2 3 1 170 171
+172 1 2 3 1 171 172
+173 1 2 3 1 172 173
+174 1 2 3 1 173 174
+175 1 2 3 1 174 175
+176 1 2 3 1 175 176
+177 1 2 3 1 176 177
+178 1 2 3 1 177 178
+179 1 2 3 1 178 179
+180 1 2 3 1 179 180
+181 1 2 3 1 180 181
+182 1 2 3 1 181 182
+183 1 2 3 1 182 183
+184 1 2 3 1 183 184
+185 1 2 3 1 184 185
+186 1 2 3 1 185 186
+187 1 2 3 1 186 187
+188 1 2 3 1 187 188
+189 1 2 3 1 188 189
+190 1 2 3 1 189 190
+191 1 2 3 1 190 191
+192 1 2 3 1 191 192
+193 1 2 3 1 192 193
+194 1 2 3 1 193 194
+195 1 2 3 1 194 195
+196 1 2 3 1 195 196
+197 1 2 3 1 196 197
+198 1 2 3 1 197 198
+199 1 2 3 1 198 199
+200 1 2 3 1 199 200
+201 1 2 3 1 200 2
+$EndElements

Solver/dgConservationLawShallowWater2d.cpp

@@ -176,7 +176,7 @@ dgBoundaryCondition *dgConservationLawShallowWater2d::newBoundaryWall(){
 #include "Bindings.h"
 void dgConservationLawShallowWater2dRegisterBindings (binding *b){
   classBinding *cb = b->addClass<dgConservationLawShallowWater2d>("dgConservationLawShallowWater2d");
-  cb->setDescription("The conservative sallow water conservation law. (H, Hu, Hv)");
+  cb->setDescription("The non-conservative shallow water conservation law. (eta, u, v)");
   methodBinding *cm;
   cm = cb->addMethod("newBoundaryWall",&dgConservationLawShallowWater2d::newBoundaryWall);
   cm->setDescription("slip wall boundary");

Solver/dgConservationLawShallowWater2d.h

@@ -17,7 +17,7 @@ class dgConservationLawShallowWater2d : public dgConservationLaw {
   dataCacheDouble *newSourceTerm (dataCacheMap &cacheMap) const;
   dgConservationLawShallowWater2d() 
   {
-    _nbf = 3; // H U(=Hu) V(=Hv)
+    _nbf = 3; // eta u v
   }
   dgBoundaryCondition *newBoundaryWall();
 };

Solver/dgResidual.cpp

@@ -97,7 +97,6 @@ void dgResidualVolume::compute1Group(dgGroupOfElements &group, fullMatrix<double
       }
     }
   }
-
   // ----- 3 ---- do the redistribution at nodes using as many BLAS3 operations as there are local coordinates
   if(_convectiveFlux || _diffusiveFlux){
     for (int iUVW=0;iUVW<group.getDimUVW();iUVW++){

Solver/dgRungeKutta.cpp

 #include <stdlib.h>
 #include "dgRungeKutta.h"
+#include "dgAlgorithm.h"
 #include "dgConservationLaw.h"
 #include "dgDofContainer.h"
 #include "dgLimiter.h"
 #include "dgTransformNodalValue.h"
 #include "dgResidual.h"
 #include "dgGroupOfElements.h"
+#include <limits.h>
+#include <stdio.h>
 
 double dgRungeKutta::iterateEuler(const dgConservationLaw *claw, double dt, dgDofContainer *solution) {
   double A[] = {0};
@@ -211,6 +214,24 @@ double dgRungeKutta::nonDiagonalRK(const dgConservationLaw *claw,
   return solution->norm();
 }
 
+double dgRungeKutta::computeInvSpectralRadius(const dgConservationLaw *claw, 
+					      dgDofContainer *solution){   
+  double sr = 1.e22;
+  dgGroupCollection *groups=solution->getGroups();
+  for (int i=0;i<groups->getNbElementGroups();i++){
+    std::vector<double> DTS;
+    dgAlgorithm::computeElementaryTimeSteps(*claw, *(groups->getElementGroup(i)), solution->getGroupProxy(i), DTS);
+    for (int k=0;k<DTS.size();k++) sr = std::min(sr,DTS[k]);
+  }
+  #ifdef HAVE_MPI
+  double sr_min;
+  MPI_Allreduce((void *)&sr, &sr_min, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
+  return sr_min;
+  #else
+  return sr;
+  #endif
+}
+
 
 dgRungeKutta::dgRungeKutta():_limiter(NULL),_TransformNodalValue(NULL){}
 
@@ -247,7 +268,9 @@ void dgRungeKutta::registerBindings(binding *b) {
   cm = cb->addMethod("setLimiter",&dgRungeKutta::setLimiter);
   cm->setArgNames("limiter",NULL);
   cm->setDescription("if a limiter is set, it is applied after each RK stage");
-
+  cm = cb->addMethod("computeInvSpectralRadius",&dgRungeKutta::computeInvSpectralRadius);
+  cm->setArgNames("law","solution",NULL);
+  cm->setDescription("Returns the inverse of the spectral radius of L(U). Useful for computing stable explicit time step.");
  cm = cb->addMethod("setTransformNodalValue",&dgRungeKutta::setTransformNodalValue);
    cm->setArgNames("TransformNodalValue",NULL);
    cm->setDescription("if the Nodal values is transformed in first step of RK");

Solver/dgRungeKutta.h

@@ -5,6 +5,7 @@
 class dgConservationLaw;
 class dgDofContainer;
 class dgLimiter;
+class dgAlgorithm;
 class dgTransformNodalValue;
 class dgGroupCollection;
 class dgGroupOfElements;
@@ -21,6 +22,7 @@ class dgRungeKutta {
   dgTransformNodalValue *_TransformNodalValue;
   public:
   void setLimiter(dgLimiter *limiter) { _limiter = limiter; }
+  double computeInvSpectralRadius(const dgConservationLaw *claw, dgDofContainer *solution);
   void setTransformNodalValue(dgTransformNodalValue *TransformNodalValue) { _TransformNodalValue = TransformNodalValue; }
 
   double iterateEuler(const dgConservationLaw *claw, double dt, dgDofContainer *solution);

Solver/eigenSolver.cpp

@@ -52,7 +52,7 @@ bool eigenSolver::solve(int numEigenValues, std::string which)
 
   // set some default options
   _try(EPSSetDimensions(eps, numEigenValues, PETSC_DECIDE, PETSC_DECIDE));
-  _try(EPSSetTolerances(eps,1.e-7,100));
+  _try(EPSSetTolerances(eps,1.e-7,50));
   _try(EPSSetType(eps,EPSARNOLDI)); //EPSKRYLOVSCHUR is default
   //_try(EPSSetType(eps,EPSARPACK)); 
   //_try(EPSSetType(eps,EPSPOWER));