No commit message

No commit message

Solver/dgConservationLawPerfectGas.cpp

@@ -446,26 +446,30 @@ class dgPerfectGasLaw2d::clipToPhysics : public dataCacheDouble {
   double _presMin, _rhoMin;
 public:
   clipToPhysics(dataCacheMap &cacheMap, double presMin, double rhoMin):
-    sol(cacheMap.get("Solution",this))
+    sol(cacheMap.get("SolToClip",this))
   {
     _presMin=presMin;
     _rhoMin=rhoMin;
     _value=fullMatrix<double>(cacheMap.getNbEvaluationPoints(),4);
   };
   void _eval () { 
-    const int nQP = sol().size1();      
-    for (size_t k = 0 ; k < nQP; k++ ){
+    const int nDofs = sol().size1();
+    for (size_t k = 0 ; k < nDofs; k++ ){
       _value(k,0) = sol(k,0);
       _value(k,1) = sol(k,1);
       _value(k,2) = sol(k,1);
       _value(k,3) = sol(k,3);
-      if (sol(k,0) < _rhoMin)  
+      if (sol(k,0) < _rhoMin){
+	//printf("CL: clip rho min =%g \n", _rhoMin);
 	_value(k,0) = _rhoMin;
+      }
       double rhoV2 = sol(k,1)*sol(k,1)+sol(k,2)*sol(k,2);
       rhoV2 /= sol(k,0);
       const double p = (GAMMA-1)*(sol(k,3) - 0.5*rhoV2);
-      if (p < _presMin) 
+      if (p < _presMin) {
 	_value(k,3) = _presMin / (GAMMA-1) +  0.5 *rhoV2 ; 
+	//printf("CL: clip pres min =%g \n ", _value(k,3));
+      }
     }
   }
 };

Solver/dgSlopeLimiter.cpp

@@ -13,8 +13,8 @@ bool dgSlopeLimiter::apply ( dgDofContainer &solution,
   //TODO: make this more general   
 
   dgGroupOfFaces* group = fGroups[0];  
-  fullMatrix<double> &solleft = solution.getGroupProxy(0); //*(solution._dataProxys[0]);
-  fullMatrix<double> &solright = solution.getGroupProxy(0); //*(solution._dataProxys[0]);
+  fullMatrix<double> &solleft = solution.getGroupProxy(0);
+  fullMatrix<double> &solright = solution.getGroupProxy(0); 
   int nbFields =_claw->nbFields();    
   int totNbElems = solution.getNbElements();
 
@@ -92,48 +92,40 @@ bool dgSlopeLimiter::apply ( dgDofContainer &solution,
 
     }
   }  
-/*
+
   //  --- CLIPPING: check unphysical values
-  dataCacheMap cacheMap(group->getNbNodes());
-  dgGroupOfElements* group = eGroups[0];   
-  fullMatrix<double> &solGroup = solution.getGroupProxy(0);
-  dataCacheDouble &solutionE = cacheMap.provideData("Solution");
-  solutionE.set(fullMatrix<double>(group.getNbNodes(),nbFields));
-  dataCacheElement &cacheElement = cacheMap.getElement();
+  for (int iG = 0; iG < eGroups.size(); iG++){
+    dgGroupOfElements* egroup = eGroups[iG];  
+    fullMatrix<double> &solGroup = solution.getGroupProxy(iG);
+
+    dataCacheMap cacheMap(egroup->getNbNodes());//nbdofs for each element
+
+    dataCacheDouble &solutionE = cacheMap.provideData("SolToClip");
+    solutionE.set(fullMatrix<double>(egroup->getNbNodes(),nbFields));//e.g 3x4 for perfect gas
+
     dataCacheDouble *solutionEClipped = _claw->newClipToPhysics(cacheMap);
-  for (int iElement=0 ; iElement<group->getNbElements() ;++iElement) {
-    cacheElement.set(group.getElement(iElement));
+    dataCacheElement &cacheElement = cacheMap.getElement();
+    for (int iElement=0 ; iElement<egroup->getNbElements() ;++iElement) {
+
+      cacheElement.set(egroup->getElement(iElement));
       solutionE.setAsProxy(solGroup,iElement*nbFields,nbFields);
-    for (int iPt =0; iPt< group->getNbNodes(); iPt++) {
-      solutionE.set((*solutionEclipped)());
-    }
-  }
-  delete solutionEClipped;
+      solutionE.set((*solutionEClipped)());
 
-*/
-  //#if 0
-//   double rhomin = 1.e-3;
-//   double presmin= 1.e-3;
-//   for (int iElement=0;iElement<totNbElems;iElement++){
 //       fullMatrix<double> solElem;
-//     solElem.setAsProxy(solleft, nbFields*iElement, nbFields );   
+//       double min = 1.e-3;
+//       solElem.setAsProxy(solGroup, nbFields*iElement, nbFields );   
 //       for (int k=0;k< solElem.size1() ;k++){
-//        if (solElem(k,0) < rhomin) {
-//  	solElem(k,0) = rhomin;
-//        }
-//        double rhoV2 = 0;
-//        for (int j=0;j<2;j++) {
-// 	 double rhov = solElem(k,j+1);
-// 	 rhoV2 += rhov*rhov;
-//        }
+// 	if (solElem(k,0) < min) solElem(k,0) = min;
+// 	double rhoV2 = solElem(k,1)*solElem(k,1)+solElem(k,2)*solElem(k,2);
 // 	rhoV2 /= solElem(k,0);
-//        const double p = (1.4-1)*solElem(k,3) - 0.5*(1.4-1)* rhoV2;
-//        if (p < presmin) {
-// 	 solElem(k,3) = 0.5 *rhoV2 + presmin / (1.4-1);
-//        }
+//        const double p = (1.4-1)*(solElem(k,3) - 0.5* rhoV2);
+//        if (p < 1.e-3)  solElem(k,3) = 0.5 *rhoV2 + min / (1.4-1);
 //       }
-//   }
-  //#endif
+      
+    }
+    delete solutionEClipped;
+    
+  }
 
   return true;