Bindings for the limiter on the way ...

Solver/TESTCASES/RayleighTaylor.lua

@@ -72,7 +72,7 @@ law:addBoundaryCondition('Top',law:newOutsideValueBoundary(FS))
 GC=dgGroupCollection(myModel,2,order)
 solution=dgDofContainer(GC,4)
 solution:L2Projection(FS)
--- limiter=dgSlopeLimiter()
+limiter=dgSlopeLimiter(law)
 -- limiter:apply(solution)
 GC:buildGroupsOfInterfaces(myModel,2,order)
 
@@ -87,7 +87,7 @@ LC = 0.1*.1
 dt = .0003;
 print('DT=',dt)
 RK=dgRungeKutta()
--- RK:setLimiter(limiter)
+RK:setLimiter(limiter)
 for i=1,10000 do
 --    norm = DG:RK44_limiter(dt)
     norm = RK:iterate44(law,dt,solution)

Solver/dgLimiter.cpp

@@ -4,29 +4,29 @@
 #include "function.h"
 
 //----------------------------------------------------------------------------------   
-bool dgSlopeLimiter::apply ( dgDofContainer &solution)
+bool dgSlopeLimiter::apply ( dgDofContainer *solution)
 {    
-  dgGroupCollection &groups=*solution.getGroups();
-  solution.scatter();
+  dgGroupCollection *groups=solution->getGroups();
+  solution->scatter();
   int nbFields =_claw->getNbFields();    
 	
   // first compute max and min of all fields for all stencils    
   //----------------------------------------------------------   
-  dgDofContainer MIN(&groups, nbFields);
-  dgDofContainer MAX(&groups, nbFields);
+  dgDofContainer MIN(groups, nbFields);
+  dgDofContainer MAX(groups, nbFields);
 
   MIN.setAll ( 1.e22);  
   MAX.setAll (-1.e22);  
 
   int iElementL, iElementR, fSize; 	 
   fullMatrix<double> TempL, TempR;
-  for( int iGFace=0; iGFace<groups.getNbFaceGroups(); iGFace++) {
-    dgGroupOfFaces* group = groups.getFaceGroup(iGFace);  
+  for( int iGFace=0; iGFace<groups->getNbFaceGroups(); iGFace++) {
+    dgGroupOfFaces* group = groups->getFaceGroup(iGFace);  
     const dgGroupOfElements *groupLeft = &group->getGroupLeft();
     const dgGroupOfElements *groupRight = &group->getGroupRight();
 
-    fullMatrix<double> &solleft = solution.getGroupProxy(groupLeft);
-    fullMatrix<double> &solright = solution.getGroupProxy(groupRight); 
+    fullMatrix<double> &solleft = solution->getGroupProxy(groupLeft);
+    fullMatrix<double> &solright = solution->getGroupProxy(groupRight); 
     fullMatrix<double> &MINLeft = MIN.getGroupProxy(groupLeft);
     fullMatrix<double> &MAXLeft = MAX.getGroupProxy(groupLeft);
     fullMatrix<double> &MINRight = MIN.getGroupProxy(groupRight);
@@ -62,9 +62,9 @@ bool dgSlopeLimiter::apply ( dgDofContainer &solution)
   // then limit the solution  
   //----------------------------------------------------------   
 
-  for (int iGroup=0 ; iGroup<groups.getNbElementGroups() ; iGroup++) {
-    dgGroupOfElements &group = *groups.getElementGroup(iGroup);
-    fullMatrix<double> &sol = solution.getGroupProxy(iGroup);
+  for (int iGroup=0 ; iGroup<groups->getNbElementGroups() ; iGroup++) {
+    dgGroupOfElements &group = *groups->getElementGroup(iGroup);
+    fullMatrix<double> &sol = solution->getGroupProxy(iGroup);
     fullMatrix<double> &MAXG = MAX.getGroupProxy(iGroup);
     fullMatrix<double> &MING = MIN.getGroupProxy(iGroup);
     fullMatrix<double> Temp;  
@@ -103,9 +103,9 @@ bool dgSlopeLimiter::apply ( dgDofContainer &solution)
     }  
   }
   //  --- CLIPPING: check unphysical values
-  for (int iG = 0; iG < groups.getNbElementGroups(); iG++){
-    dgGroupOfElements* egroup = groups.getElementGroup(iG);  
-    fullMatrix<double> &solGroup = solution.getGroupProxy(iG);
+  for (int iG = 0; iG < groups->getNbElementGroups(); iG++){
+    dgGroupOfElements* egroup = groups->getElementGroup(iG);  
+    fullMatrix<double> &solGroup = solution->getGroupProxy(iG);
 
     dataCacheMap cacheMap(egroup->getNbNodes());//nbdofs for each element
     dataCacheDouble &solutionE = cacheMap.provideData("Solution");
@@ -121,9 +121,17 @@ bool dgSlopeLimiter::apply ( dgDofContainer &solution)
   return true; 
 }
 
-/*
 #include "Bindings.h"
 
+void dgLimiter::registerBindings(binding *b) {
+  classBinding *cb = b->addClass<dgLimiter>("dgLimiter");
+  cb->setDescription("Parent class for limiters");
+  methodBinding *cm;
+//  cm = cb->addMethod("apply",&dgLimiter::apply);
+//  cm->setArgNames("solution",NULL);
+//  cm->setDescription("apply the limiter on the solution");
+}
+
 void dgSlopeLimiter::registerBindings(binding *b) {
   classBinding *cb = b->addClass<dgSlopeLimiter>("dgSlopeLimiter");
   cb->setDescription("The usual DG slope limiter: keep the mean, reduces uniformly the slope until it does not overshoot the neighbors' mean");
@@ -131,9 +139,5 @@ void dgSlopeLimiter::registerBindings(binding *b) {
   cm = cb->setConstructor<dgSlopeLimiter,dgConservationLaw *>();
   cm->setDescription("A new explicit slope limiter");
   cm->setArgNames("law",NULL);
-//  cm = cb->addMethod("apply",&dgLimiter::apply);
-//  cm->setArgNames("solution",NULL);
-//  cm->setDescription("apply the limiter on the solution");
 }
 
-*/

Solver/dgLimiter.h

@@ -13,14 +13,15 @@ protected:
   dgConservationLaw *_claw;
 public:
   dgLimiter (dgConservationLaw *claw) : _claw(claw) {}
-  virtual bool apply ( dgDofContainer &sol)=0;
+  virtual bool apply ( dgDofContainer *sol)=0;
+  static void registerBindings(binding *b);
 };
 
 class dgSlopeLimiter : public dgLimiter{
 public :
   dgSlopeLimiter (dgConservationLaw *claw) : dgLimiter (claw) {}
-  virtual bool apply ( dgDofContainer &solution);
-  //static void registerBindings(binding *b);
+  virtual bool apply ( dgDofContainer *solution);
+  static void registerBindings(binding *b);
 };
 
 #endif

Solver/dgRungeKutta.cpp

@@ -49,7 +49,7 @@ double dgRungeKutta::diagonalRK (const dgConservationLaw *claw,
       K.axpy(*sol);
     }
 
-    if (_limiter) _limiter->apply(K);
+    if (_limiter) _limiter->apply(&K);
     dgAlgorithm::residual(*claw,*groups,K,resd);
     K.scale(0.);
     for(int k=0; k < groups->getNbElementGroups(); k++) {
@@ -64,7 +64,7 @@ double dgRungeKutta::diagonalRK (const dgConservationLaw *claw,
     }
     Unp.axpy(K,b[j]*dt);
   }
-  if (_limiter) _limiter->apply(Unp);
+  if (_limiter) _limiter->apply(&Unp);
   sol->scale(0.);
   sol->axpy(Unp);
   return sol->norm();
@@ -90,7 +90,7 @@ void dgRungeKutta::registerBindings(binding *b) {
   cm = cb->addMethod("iterate44",&dgRungeKutta::iterate44);
   cm->setArgNames("law","dt","solution",NULL);
   cm->setDescription("update solution by doing classical fourth order runge-kutta step of time step dt for the conservation law");
-//  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("setLimiter",&dgRungeKutta::setLimiter);
+  cm->setArgNames("limiter",NULL);
+  cm->setDescription("if a limiter is set, it is applied after each RK stage");
 }

Solver/dgSystemOfEquations.cpp

@@ -142,7 +142,7 @@ void dgSystemOfEquations::exportSolution(std::string outputFile){
 
 void dgSystemOfEquations::limitSolution(){
   dgLimiter *sl = new dgSlopeLimiter(_claw);
-  sl->apply(*_solution);
+  sl->apply(_solution);
 
   delete sl;
 }