diff --git a/Solver/dgConservationLawPerfectGas.cpp b/Solver/dgConservationLawPerfectGas.cpp index 11ac62e37e33a8be5c9a2710e3be6a435e2242bb..49a2aaafa52624b340309097b0a9f71cb51c4c5f 100644 --- a/Solver/dgConservationLawPerfectGas.cpp +++ b/Solver/dgConservationLawPerfectGas.cpp @@ -446,27 +446,31 @@ 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)); + } + } } }; diff --git a/Solver/dgSlopeLimiter.cpp b/Solver/dgSlopeLimiter.cpp index 2f475d9743953be93f38ff20b3196e0544f9556a..026b5e990e02274e8fa266bfbcc94ade6f51ae89 100644 --- a/Solver/dgSlopeLimiter.cpp +++ b/Solver/dgSlopeLimiter.cpp @@ -12,9 +12,9 @@ bool dgSlopeLimiter::apply ( dgDofContainer &solution, //WARNING: ONLY FOR 1 GROUP OF FACES //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]); + dgGroupOfFaces* group = fGroups[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(); - dataCacheDouble *solutionEClipped = _claw->newClipToPhysics(cacheMap); - for (int iElement=0 ; iElement<group->getNbElements() ;++iElement) { - cacheElement.set(group.getElement(iElement)); - solutionE.setAsProxy(solGroup,iElement*nbFields,nbFields); - for (int iPt =0; iPt< group->getNbNodes(); iPt++) { - solutionE.set((*solutionEclipped)()); + 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); + dataCacheElement &cacheElement = cacheMap.getElement(); + for (int iElement=0 ; iElement<egroup->getNbElements() ;++iElement) { + + cacheElement.set(egroup->getElement(iElement)); + solutionE.setAsProxy(solGroup,iElement*nbFields,nbFields); + solutionE.set((*solutionEClipped)()); + +// fullMatrix<double> solElem; +// double min = 1.e-3; +// solElem.setAsProxy(solGroup, nbFields*iElement, nbFields ); +// for (int k=0;k< solElem.size1() ;k++){ +// 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* rhoV2); +// if (p < 1.e-3) solElem(k,3) = 0.5 *rhoV2 + min / (1.4-1); +// } + } + delete solutionEClipped; + } - delete 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 ); -// 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; -// } -// 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); -// } -// } -// } - //#endif return true;