diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index 5521ba2b0118ae3509eb276a0b25a0bbd099eaef..9e14abe960f7972a376269446b0fa04d5816171f 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -16,65 +16,74 @@
\int \vec{f} \cdot \grad \phi dv
*/
-void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe useless here)
- const dgConservationLaw &claw, // the conservation law
- const dgGroupOfElements & group,
- fullMatrix<double> &solution,
- fullMatrix<double> &residual // the residual
- )
-{
+class dgResidualVolume {
+ dataCacheMap _cacheMap;
+ const dgConservationLaw &_claw;
+ int _nbFields;
+ dataCacheElement &_cacheElement;
+ dataCacheDouble &_UVW, &_solutionQPe, &_gradientSolutionQPe;
+ dataCacheDouble *_sourceTerm, *_convectiveFlux, *_diffusiveFlux;
+ public:
+ dgResidualVolume (const dgConservationLaw &claw);
+ void compute1Group(dgGroupOfElements &group, fullMatrix<double> &solution, fullMatrix<double> &residual);
+ void computeAndMap1Group(dgGroupOfElements &group, dgDofContainer &solution, dgDofContainer &residual);
+};
+
+dgResidualVolume::dgResidualVolume(const dgConservationLaw &claw):
+ _claw(claw),
+ _nbFields(_claw.getNbFields()),
+ _cacheElement(_cacheMap.getElement()),
+ _UVW(_cacheMap.provideData("UVW", 1, 3)),
+ _solutionQPe(_cacheMap.provideData("Solution", 1, _nbFields)),
+ _gradientSolutionQPe(_cacheMap.provideData("SolutionGradient", 3, _nbFields)),
+ _sourceTerm(_claw.newSourceTerm(_cacheMap)),
+ _convectiveFlux(_claw.newConvectiveFlux(_cacheMap)),
+ _diffusiveFlux(_claw.newDiffusiveFlux(_cacheMap))
+{
+}
+
+void dgResidualVolume::compute1Group(dgGroupOfElements &group, fullMatrix<double> &solution, fullMatrix<double> &residual)
+{
+ residual.scale(0);
+ _cacheMap.setNbEvaluationPoints(group.getNbIntegrationPoints());
+ _UVW.set(group.getIntegrationPointsMatrix());
// ----- 1 ---- get the solution at quadrature points
// ----- 1.1 --- allocate a matrix of size (nbFields * nbElements, nbQuadraturePoints)
- int nbFields = claw.getNbFields();
- fullMatrix<double> solutionQP (group.getNbIntegrationPoints(),group.getNbElements() * nbFields);
+ fullMatrix<double> solutionQP (group.getNbIntegrationPoints(),group.getNbElements() * _nbFields);
// ----- 1.2 --- multiply the solution by the collocation matrix
group.getCollocationMatrix().mult(solution , solutionQP);
-
// ----- 2 ---- compute all fluxes (diffusive and convective) at integration points
// ----- 2.1 --- allocate elementary fluxes (compued in XYZ coordinates)
- fullMatrix<double> fConv[3] = {fullMatrix<double>( group.getNbIntegrationPoints(), nbFields ),
- fullMatrix<double>( group.getNbIntegrationPoints(), nbFields ),
- fullMatrix<double>( group.getNbIntegrationPoints(), nbFields )};
- fullMatrix<double> fDiff[3] = {fullMatrix<double>( group.getNbIntegrationPoints(), nbFields ),
- fullMatrix<double>( group.getNbIntegrationPoints(), nbFields ),
- fullMatrix<double>( group.getNbIntegrationPoints(), nbFields )};
- fullMatrix<double> Source = fullMatrix<double> (group.getNbIntegrationPoints(),group.getNbElements() * nbFields);
+ fullMatrix<double> fConv[3] = {fullMatrix<double>( group.getNbIntegrationPoints(), _nbFields ),
+ fullMatrix<double>( group.getNbIntegrationPoints(), _nbFields ),
+ fullMatrix<double>( group.getNbIntegrationPoints(), _nbFields )};
+ fullMatrix<double> fDiff[3] = {fullMatrix<double>( group.getNbIntegrationPoints(), _nbFields ),
+ fullMatrix<double>( group.getNbIntegrationPoints(), _nbFields ),
+ fullMatrix<double>( group.getNbIntegrationPoints(), _nbFields )};
+ fullMatrix<double> Source = fullMatrix<double> (group.getNbIntegrationPoints(),group.getNbElements() * _nbFields);
// ----- 2.2 --- allocate parametric fluxes (computed in UVW coordinates) for all elements at all integration points
- fullMatrix<double> Fuvw[3] = {fullMatrix<double> ( group.getNbIntegrationPoints(), group.getNbElements() * nbFields),
- fullMatrix<double> (group.getNbIntegrationPoints(), group.getNbElements() * nbFields),
- fullMatrix<double> (group.getNbIntegrationPoints(), group.getNbElements() * nbFields)};
- dataCacheMap cacheMap;
- cacheMap.setNbEvaluationPoints(group.getNbIntegrationPoints());
- dataCacheElement &cacheElement = cacheMap.getElement();
- // provided dataCache
- cacheMap.provideData("UVW",1,3).set(group.getIntegrationPointsMatrix());
- dataCacheDouble &solutionQPe = cacheMap.provideData("Solution",1,nbFields);
- dataCacheDouble &gradientSolutionQPe = cacheMap.provideData("SolutionGradient",3,nbFields);
- // needed dataCache
- dataCacheDouble *sourceTerm = claw.newSourceTerm(cacheMap);
- // fluxes in XYZ coordinates;
- dataCacheDouble *convectiveFlux = claw.newConvectiveFlux(cacheMap);
- dataCacheDouble *diffusiveFlux = claw.newDiffusiveFlux(cacheMap);
-
+ fullMatrix<double> Fuvw[3] = {fullMatrix<double> ( group.getNbIntegrationPoints(), group.getNbElements() * _nbFields),
+ fullMatrix<double> (group.getNbIntegrationPoints(), group.getNbElements() * _nbFields),
+ fullMatrix<double> (group.getNbIntegrationPoints(), group.getNbElements() * _nbFields)};
fullMatrix<double> fuvwe;
fullMatrix<double> source;
fullMatrix<double> dPsiDx,dofs;
// ----- 2.3 --- iterate on elements
for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {
// ----- 2.3.1 --- build a small object that contains elementary solution, jacobians, gmsh element
- solutionQPe.setAsProxy(solutionQP, iElement*nbFields, nbFields );
-
- if(gradientSolutionQPe.somethingDependOnMe()){
+ _solutionQPe.setAsProxy(solutionQP, iElement*_nbFields, _nbFields );
+
+ if(_gradientSolutionQPe.somethingDependOnMe()){
dPsiDx.setAsProxy(group.getDPsiDx(),iElement*group.getNbNodes(),group.getNbNodes());
- dofs.setAsProxy(solution, nbFields*iElement, nbFields);
- dPsiDx.mult(dofs, gradientSolutionQPe.set());
+ dofs.setAsProxy(solution, _nbFields*iElement, _nbFields);
+ dPsiDx.mult(dofs, _gradientSolutionQPe.set());
}
- cacheElement.set(group.getElement(iElement));
- if(convectiveFlux || diffusiveFlux) {
+ _cacheElement.set(group.getElement(iElement));
+ if(_convectiveFlux || _diffusiveFlux) {
// ----- 2.3.3 --- compute fluxes in UVW coordinates
for (int iUVW=0;iUVW<group.getDimUVW();iUVW++) {
// ----- 2.3.3.1 --- get a proxy on the big local flux matrix
- fuvwe.setAsProxy(Fuvw[iUVW], iElement*nbFields,nbFields);
+ fuvwe.setAsProxy(Fuvw[iUVW], iElement*_nbFields,_nbFields);
// ----- 2.3.3.1 --- compute \vec{f}^{UVW} =\vec{f}^{XYZ} J^{-1} and multiply bt detJ
for (int iXYZ=0;iXYZ<group.getDimXYZ();iXYZ++) {
for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
@@ -84,40 +93,85 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
const double detJ = group.getDetJ (iElement, iPt);
const double factor = invJ * detJ;
// compute fluxes in the reference coordinate system at this integration point
- for (int k=0;k<nbFields;k++){
- if(convectiveFlux)
- fuvwe(iPt,k) += ( (*convectiveFlux)(iPt,iXYZ*nbFields+k)) * factor;
- if(diffusiveFlux){
- fuvwe(iPt,k) += ( (*diffusiveFlux)(iPt,iXYZ*nbFields+k)) * factor;
+ for (int k=0;k<_nbFields;k++){
+ if(_convectiveFlux)
+ fuvwe(iPt,k) += ( (*_convectiveFlux)(iPt,iXYZ*_nbFields+k)) * factor;
+ if(_diffusiveFlux){
+ fuvwe(iPt,k) += ( (*_diffusiveFlux)(iPt,iXYZ*_nbFields+k)) * factor;
}
}
}
}
}
}
- if (sourceTerm){
- source.setAsProxy(Source, iElement*claw.getNbFields(),claw.getNbFields());
+ if (_sourceTerm){
+ source.setAsProxy(Source, iElement*_nbFields,_nbFields);
for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
const double detJ = group.getDetJ (iElement, iPt);
- for (int k=0;k<nbFields;k++){
- source(iPt,k) = (*sourceTerm)(iPt,k)*detJ;
+ for (int k=0;k<_nbFields;k++){
+ source(iPt,k) = (*_sourceTerm)(iPt,k)*detJ;
}
}
}
}
+
// ----- 3 ---- do the redistribution at nodes using as many BLAS3 operations as there are local coordinates
- if(convectiveFlux || diffusiveFlux){
+ if(_convectiveFlux || _diffusiveFlux){
for (int iUVW=0;iUVW<group.getDimUVW();iUVW++){
residual.gemm(group.getFluxRedistributionMatrix(iUVW),Fuvw[iUVW]);
}
}
- if(sourceTerm){
+ if(_sourceTerm){
residual.gemm(group.getSourceRedistributionMatrix(),Source);
}
}
-void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (maybe useless here)
- const dgConservationLaw &claw, // the conservation law
+void dgResidualVolume::computeAndMap1Group(dgGroupOfElements &group, dgDofContainer &solution, dgDofContainer &residual)
+{
+ compute1Group (group, solution.getGroupProxy(&group), residual.getGroupProxy(&group));
+}
+
+class dgResidualInterface {
+ dataCacheMap _cacheMapLeft, _cacheMapRight;
+ const dgConservationLaw &_claw;
+ int _nbFields;
+ dataCacheElement &_cacheElementLeft, &_cacheElementRight;
+ dataCacheDouble &_uvwLeft, &_uvwRight, &_solutionQPLeft, &_solutionQPRight, &_gradientSolutionLeft, &_gradientSolutionRight;
+ dataCacheDouble &_normals;
+ dataCacheDouble *_riemannSolver, *_maximumDiffusivityLeft,*_maximumDiffusivityRight, *_diffusiveFluxLeft, *_diffusiveFluxRight;
+ public:
+ dgResidualInterface (const dgConservationLaw &claw);
+ void compute1Group ( //dofManager &dof, // the DOF manager (maybe useless here)
+ dgGroupOfFaces &group,
+ const fullMatrix<double> &solution, // solution !! at faces nodes
+ fullMatrix<double> &solutionLeft,
+ fullMatrix<double> &solutionRight,
+ fullMatrix<double> &residual // residual !! at faces nodes
+ );
+ void computeAndMap1Group (dgGroupOfFaces &faces, dgDofContainer &solution, dgDofContainer &residual);
+};
+
+dgResidualInterface::dgResidualInterface (const dgConservationLaw &claw) :
+ _claw(claw),
+ _nbFields(claw.getNbFields()),
+ _cacheElementLeft(_cacheMapLeft.getElement()),
+ _cacheElementRight(_cacheMapRight.getElement()),
+ _uvwLeft(_cacheMapLeft.provideData("UVW",1,3)),
+ _uvwRight(_cacheMapRight.provideData("UVW",1,3)),
+ _solutionQPLeft(_cacheMapLeft.provideData("Solution",1,_nbFields)),
+ _solutionQPRight(_cacheMapRight.provideData("Solution",1,_nbFields)),
+ _gradientSolutionLeft(_cacheMapLeft.provideData("SolutionGradient",3,_nbFields)),
+ _gradientSolutionRight(_cacheMapRight.provideData("SolutionGradient",3,_nbFields)),
+ _normals(_cacheMapLeft.provideData("Normals",1,1)),
+ _riemannSolver(claw.newRiemannSolver(_cacheMapLeft,_cacheMapRight)),
+ _diffusiveFluxLeft(claw.newDiffusiveFlux(_cacheMapLeft)),
+ _diffusiveFluxRight(claw.newDiffusiveFlux(_cacheMapRight)),
+ _maximumDiffusivityLeft(claw.newMaximumDiffusivity(_cacheMapLeft)),
+ _maximumDiffusivityRight(claw.newMaximumDiffusivity(_cacheMapRight))
+{
+}
+
+void dgResidualInterface::compute1Group ( //dofManager &dof, // the DOF manager (maybe useless here)
dgGroupOfFaces &group,
const fullMatrix<double> &solution, // solution !! at faces nodes
fullMatrix<double> &solutionLeft,
@@ -127,297 +181,120 @@ void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (may
{
// A) global operations before entering the loop over the faces
// A1 ) copy locally some references from the group for the sake of lisibility
- int nbFields = claw.getNbFields();
int nbNodesLeft = group.getGroupLeft().getNbNodes();
int nbNodesRight = group.getGroupRight().getNbNodes();
int nbFaces = group.getNbElements();
//get matrices needed to compute the gradient on both sides
fullMatrix<double> &DPsiLeftDx = group.getDPsiLeftDxMatrix();
fullMatrix<double> &DPsiRightDx = group.getDPsiRightDxMatrix();
-
// ----- 1 ---- get the solution at quadrature points
- fullMatrix<double> solutionQP (group.getNbIntegrationPoints(),nbFaces * nbFields*2);
+ fullMatrix<double> solutionQP (group.getNbIntegrationPoints(),nbFaces * _nbFields*2);
group.getCollocationMatrix().mult(solution, solutionQP);
// needed tocompute normal fluxes at integration points
- fullMatrix<double> NormalFluxQP ( group.getNbIntegrationPoints(), nbFaces*nbFields*2);
-
+ fullMatrix<double> NormalFluxQP ( group.getNbIntegrationPoints(), nbFaces*_nbFields*2);
// create one dataCache for each side
- dataCacheMap cacheMapLeft;
- cacheMapLeft.setNbEvaluationPoints(group.getNbIntegrationPoints());
- dataCacheMap cacheMapRight;
- cacheMapRight.setNbEvaluationPoints(group.getNbIntegrationPoints());
-
- // data this algorithm provide to the cache map (we can maybe move each data to a separate function but
- // I think It's easier like this)
- dataCacheDouble &normals = cacheMapLeft.provideData("Normals",1,1);
- dataCacheElement &cacheElementLeft = cacheMapLeft.getElement();
- dataCacheElement &cacheElementRight = cacheMapRight.getElement();
- dataCacheDouble &uvwLeft = cacheMapLeft.provideData("UVW",1,3);
- dataCacheDouble &uvwRight = cacheMapRight.provideData("UVW",1,3);
- dataCacheDouble &solutionQPLeft = cacheMapLeft.provideData("Solution",1,nbFields);
- dataCacheDouble &solutionQPRight = cacheMapRight.provideData("Solution",1,nbFields);
- dataCacheDouble &gradientSolutionLeft = cacheMapLeft.provideData("SolutionGradient",3,nbFields);
- dataCacheDouble &gradientSolutionRight = cacheMapRight.provideData("SolutionGradient",3,nbFields);
-
- // A2 ) ask the law to provide the fluxes (possibly NULL)
- dataCacheDouble *riemannSolver = claw.newRiemannSolver(cacheMapLeft,cacheMapRight);
- dataCacheDouble *diffusiveFluxLeft = claw.newDiffusiveFlux(cacheMapLeft);
- dataCacheDouble *maximumDiffusivityLeft = claw.newMaximumDiffusivity(cacheMapLeft);
- dataCacheDouble *maximumDiffusivityRight = claw.newMaximumDiffusivity(cacheMapRight);
- dataCacheDouble *diffusiveFluxRight = claw.newDiffusiveFlux(cacheMapRight);
+ _cacheMapLeft.setNbEvaluationPoints(group.getNbIntegrationPoints());
+ _cacheMapRight.setNbEvaluationPoints(group.getNbIntegrationPoints());
fullMatrix<double> dPsiLeftDx,dPsiRightDx,dofsLeft,dofsRight,normalFluxQP;
-
int p = group.getGroupLeft().getOrder();
int dim = group.getGroupLeft().getElement(0)->getDim();
-
for (int iFace=0 ; iFace < nbFaces ; ++iFace) {
// B1 ) adjust the proxies for this element
// NB : the B1 section will almost completely disapear
// if we write a function (from the function class) for moving proxies across big matrices
// give the current elements to the dataCacheMap
- cacheElementLeft.set(group.getElementLeft(iFace));
- cacheElementRight.set(group.getElementRight(iFace));
+ _cacheElementLeft.set(group.getElementLeft(iFace));
+ _cacheElementRight.set(group.getElementRight(iFace));
// proxies for the solution
- solutionQPLeft.setAsProxy(solutionQP, iFace*2*nbFields, nbFields);
- solutionQPRight.setAsProxy(solutionQP, (iFace*2+1)*nbFields, nbFields);
- normals.setAsProxy(group.getNormals(), iFace*group.getNbIntegrationPoints(),group.getNbIntegrationPoints());
+ _solutionQPLeft.setAsProxy(solutionQP, iFace*2*_nbFields, _nbFields);
+ _solutionQPRight.setAsProxy(solutionQP, (iFace*2+1)*_nbFields, _nbFields);
+ _normals.setAsProxy(group.getNormals(), iFace*group.getNbIntegrationPoints(),group.getNbIntegrationPoints());
// proxies needed to compute the gradient of the solution and the IP term
dPsiLeftDx.setAsProxy(DPsiLeftDx,iFace*nbNodesLeft,nbNodesLeft);
dPsiRightDx.setAsProxy(DPsiRightDx,iFace*nbNodesRight,nbNodesRight);
- dofsLeft.setAsProxy(solutionLeft, nbFields*group.getElementLeftId(iFace), nbFields);
- dofsRight.setAsProxy(solutionRight, nbFields*group.getElementRightId(iFace), nbFields);
- uvwLeft.setAsProxy(group.getIntegrationOnElementLeft(iFace));
- uvwRight.setAsProxy(group.getIntegrationOnElementRight(iFace));
+ dofsLeft.setAsProxy(solutionLeft, _nbFields*group.getElementLeftId(iFace), _nbFields);
+ dofsRight.setAsProxy(solutionRight, _nbFields*group.getElementRightId(iFace), _nbFields);
+ _uvwLeft.setAsProxy(group.getIntegrationOnElementLeft(iFace));
+ _uvwRight.setAsProxy(group.getIntegrationOnElementRight(iFace));
// proxies for the flux
- normalFluxQP.setAsProxy(NormalFluxQP, iFace*nbFields*2, nbFields*2);
-
+ normalFluxQP.setAsProxy(NormalFluxQP, iFace*_nbFields*2, _nbFields*2);
// B2 ) compute the gradient of the solution
- if(gradientSolutionLeft.somethingDependOnMe()){
- dPsiLeftDx.mult(dofsLeft, gradientSolutionLeft.set());
- dPsiRightDx.mult(dofsRight, gradientSolutionRight.set());
+ if(_gradientSolutionLeft.somethingDependOnMe()){
+ dPsiLeftDx.mult(dofsLeft, _gradientSolutionLeft.set());
+ dPsiRightDx.mult(dofsRight, _gradientSolutionRight.set());
}
-
// B3 ) compute fluxes
- if (diffusiveFluxLeft) {
+ if (_diffusiveFluxLeft) {
for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
const double detJ = group.getDetJ (iFace, iPt);
//just for the lisibility :
- const fullMatrix<double> &dfl = (*diffusiveFluxLeft)();
- const fullMatrix<double> &dfr = (*diffusiveFluxRight)();
- for (int k=0;k<nbFields;k++) {
+ const fullMatrix<double> &dfl = (*_diffusiveFluxLeft)();
+ const fullMatrix<double> &dfr = (*_diffusiveFluxRight)();
+ for (int k=0;k<_nbFields;k++) {
double meanNormalFlux =
- ((dfl(iPt,k+nbFields*0 )+dfr(iPt,k+nbFields*0)) * normals(0,iPt)
- +(dfl(iPt,k+nbFields*1 )+dfr(iPt,k+nbFields*1)) * normals(1,iPt)
- +(dfl(iPt,k+nbFields*2 )+dfr(iPt,k+nbFields*2)) * normals(2,iPt))/2;
+ ((dfl(iPt,k+_nbFields*0 )+dfr(iPt,k+_nbFields*0)) * _normals(0,iPt)
+ +(dfl(iPt,k+_nbFields*1 )+dfr(iPt,k+_nbFields*1)) * _normals(1,iPt)
+ +(dfl(iPt,k+_nbFields*2 )+dfr(iPt,k+_nbFields*2)) * _normals(2,iPt))/2;
double minl = std::min(group.getElementVolumeLeft(iFace),
group.getElementVolumeRight(iFace)
)/group.getInterfaceSurface(iFace);
- double nu = std::max((*maximumDiffusivityRight)(iPt,0),(*maximumDiffusivityLeft)(iPt,0));
+ double nu = std::max((*_maximumDiffusivityRight)(iPt,0),(*_maximumDiffusivityLeft)(iPt,0));
double mu = (p+1)*(p+dim)/dim*nu/minl;
- double solutionJumpPenalty = (solutionQPLeft(iPt,k)-solutionQPRight(iPt,k))/2*mu;
+ double solutionJumpPenalty = (_solutionQPLeft(iPt,k)-_solutionQPRight(iPt,k))/2*mu;
normalFluxQP(iPt,k) -= (meanNormalFlux+solutionJumpPenalty)*detJ;
- normalFluxQP(iPt,k+nbFields) += (meanNormalFlux+solutionJumpPenalty)*detJ;
+ normalFluxQP(iPt,k+_nbFields) += (meanNormalFlux+solutionJumpPenalty)*detJ;
}
}
}
- if (riemannSolver) {
+ if (_riemannSolver) {
for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
const double detJ = group.getDetJ (iFace, iPt);
- for (int k=0;k<nbFields*2;k++)
- normalFluxQP(iPt,k) += (*riemannSolver)(iPt,k)*detJ;
+ for (int k=0;k<_nbFields*2;k++)
+ normalFluxQP(iPt,k) += (*_riemannSolver)(iPt,k)*detJ;
}
}
}
-
// C ) redistribute the flux to the residual (at Faces nodes)
-
- if(riemannSolver || diffusiveFluxLeft)
+ if(_riemannSolver || _diffusiveFluxLeft)
residual.gemm(group.getRedistributionMatrix(),NormalFluxQP);
}
-void dgAlgorithm::multAddInverseMassMatrix ( /*dofManager &dof,*/
- const dgGroupOfElements & group,
- fullMatrix<double> &residu,
- fullMatrix<double> &sol)
+
+void dgResidualInterface::computeAndMap1Group (dgGroupOfFaces &faces, dgDofContainer &solution, dgDofContainer &residual)
{
- fullMatrix<double> residuEl;
- fullMatrix<double> solEl;
- fullMatrix<double> iMassEl;
- int nFields=sol.size2()/group.getNbElements();
- for(int i=0;i<group.getNbElements();i++) {
- residuEl.setAsProxy(residu,i*nFields,nFields);
- solEl.setAsProxy(sol,i*nFields,nFields);
- iMassEl.setAsProxy(group.getInverseMassMatrix(),i*group.getNbNodes(),group.getNbNodes());
- solEl.gemm(iMassEl,residuEl);
- }
+ fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*2*_nbFields);
+ fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*2*_nbFields);
+ faces.mapToInterface(_nbFields, solution.getGroupProxy(&faces.getGroupLeft()), solution.getGroupProxy(&faces.getGroupRight()), solInterface);
+ compute1Group(faces,solInterface,solution.getGroupProxy(&faces.getGroupLeft()), solution.getGroupProxy(&faces.getGroupRight()),residuInterface);
+ faces.mapFromInterface(_nbFields, residuInterface,residual.getGroupProxy(&faces.getGroupLeft()), residual.getGroupProxy(&faces.getGroupRight()));
}
+class dgResidualBoundary {
+ const dgConservationLaw &_claw;
+ public :
+ void compute1Group ( //dofManager &dof, // the DOF manager (maybe useless here)
+ dgGroupOfFaces &group,
+ const fullMatrix<double> &solution, // solution !! at faces nodes
+ fullMatrix<double> &solutionLeft,
+ fullMatrix<double> &residual // residual !! at faces nodes
+ );
+ void computeAndMap1Group (dgGroupOfFaces &faces, dgDofContainer &solution, dgDofContainer &residual);
+ dgResidualBoundary (const dgConservationLaw &claw);
+};
+dgResidualBoundary::dgResidualBoundary(const dgConservationLaw &claw):_claw(claw){
+}
-/*void dgAlgorithm::multirateRungeKutta (const dgConservationLaw &claw, // conservation law
- dgGroupCollection &groups,
- double h, // time-step
- dgDofContainer &sol,
- dgDofContainer &resd,
- int orderRK) // order of RK integrator
- {
-
- // U_{n+1}=U_n+h*(SUM_i a_i*K_i)
- // K_i=M^(-1)R(U_n+b_i*K_{i-1})
-
- int nStages=10;
-// classical RK44
-// double A[4][4]={
-// {0, 0, 0, 0},
-// {1.0/2.0, 0, 0 ,0},
-// {0, 1.0/2.0, 0, 0},
-// {0, 0, 1, 0}
-// };
-// double b[4]={1.0/6.0, 1.0/3.0, 1.0/3.0, 1.0/6.0};
-// double c[4]={0, 1.0/2.0, 1.0/2.0, 1};
-
-// 3/8 RK44
-// double A[4][4]={
-// {0, 0, 0, 0},
-// {1.0/3.0, 0, 0 ,0},
-// {-1.0/3.0, 1.0, 0, 0},
-// {1, -1, 1, 0}
-// };
-// double b[4]={1.0/8.0, 3.0/8.0, 3.0/8.0, 1.0/8.0};
-// double c[4]={0, 1.0/3.0, 2.0/3.0, 1};
-
-// RK43 from Schlegel et al. JCAM 2009
-// double A[4][4]={
-// {0, 0, 0, 0},
-// {1.0/2.0, 0, 0 ,0},
-// {-1.0/6.0, 2.0/3.0, 0, 0},
-// {1.0/3.0, -1.0/3.0, 1, 0}
-// };
-// double b[4]={1.0/6.0, 1.0/3.0, 1.0/3.0, 1.0/6.0};
-// double c[4]={0, 1.0/2.0, 1.0/2.0, 1};
-
- double A[10][10];
- double *b;
- double *c;
-// Small step RK43
- double A1[10][10]={
-{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-{1.0/4.0 ,0, 0, 0, 0, 0, 0, 0, 0, 0},
-{-1.0/12.0, 1.0/3.0, 0, 0, 0, 0, 0, 0, 0, 0},
-{1.0/6.0, -1.0/6.0, 1.0/2.0, 0, 0, 0, 0, 0, 0, 0},
-{1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 0, 0, 0, 0, 0},
-{1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 0, 0, 0, 0, 0},
-{1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 1.0/4.0, 0, 0, 0, 0},
-{1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, -1.0/12.0, 1.0/3.0, 0, 0, 0},
-{1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 1.0/6.0, -1.0/6.0, 1.0/2.0, 0, 0},
-{1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0}
- };
- double b1[10]={1.0/12.0, 1.0/6.0, 1.0/6.0,1.0/12.0,0,1.0/12.0, 1.0/6.0, 1.0/6.0,1.0/12.0,0};
- double c1[10]={0, 1.0/4.0, 1.0/4.0, 1.0/2.0, 1.0/2.0, 1.0/2.0, 3.0/4.0, 3.0/4.0, 1, 1 };
-
-// Big step RK43
- double A2[10][10]={
-{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-{1.0/4.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
-{1.0/4.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
-{1.0/2.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
-{1.0/2.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
-{-1.0/6.0, 0, 0, 0, 2.0/3.0, 0, 0, 0, 0, 0},
-{1.0/12.0, 0, 0, 0, 1.0/6.0, 1.0/2.0, 0, 0, 0, 0},
-{1.0/12.0, 0, 0, 0, 1.0/6.0, 1.0/2.0, 0, 0, 0, 0},
-{1.0/3.0 , 0, 0, 0, -1.0/3.0, 1, 0, 0, 0, 0},
-{1.0/3.0 , 0, 0, 0, -1.0/3.0, 1, 0, 0, 0, 0},
- };
- double b2[10]={1.0/6.0, 0 , 0 , 0 , 1.0/3.0,1.0/3.0, 0 , 0 , 0 , 1.0/6.0};
- double c2[10]={0, 1.0/4.0, 1.0/4.0, 1.0/2.0, 1.0/2.0, 1.0/2.0, 3.0/4.0, 3.0/4.0, 1, 1 };
-
-
- // fullMatrix<double> K(sol);
- // Current updated solution
- // fullMatrix<double> Unp(sol);
- fullMatrix<double> residuEl, KEl;
- fullMatrix<double> iMassEl;
-
- int nbFields = claw.getNbFields();
-
- dgDofContainer **K;
- K=new dgDofContainer*[nStages];
- for(int i=0;i<nStages;i++){
- K[i]=new dgDofContainer(&groups,nbFields);
- K[i]->scale(0.);
- }
- dgDofContainer Unp (&groups,nbFields);
- dgDofContainer tmp (&groups,nbFields);
-
- Unp.scale(0.0);
- Unp.axpy(sol);
- // Case of 2 groups: first with big steps, second with small steps
- for(int j=0; j<nStages;j++){
- tmp.scale(0.0);
- tmp.axpy(sol);
- for(int k=0;k < groups.getNbElementGroups();k++) {
- if (k==0) {
- for (int ii=0; ii<10; ii++) {
- for (int jj=0; jj<10; jj++) {
- A[ii][jj]=A2[ii][jj];
- }
- }
- }
- else{
- for (int ii=0; ii<10; ii++) {
- for (int jj=0; jj<10; jj++) {
- A[ii][jj]=A1[ii][jj];
- }
- }
- }
- for(int i=0;i<j;i++){
- if(fabs(A[j][i])>1e-12){
- tmp.getGroupProxy(k).add(K[i]->getGroupProxy(k),h*A[j][i]);
- }
- }
- }
- dgAlgorithm::residual(claw,groups,tmp,resd);
- for(int k=0;k < groups.getNbElementGroups(); k++) {
- if (k==0) {
- b=b2;
- c=c2;
- }
- else{
- b=b1;
- c=c1;
- }
-
- dgGroupOfElements *group = groups.getElementGroup(k);
- int nbNodes = group->getNbNodes();
- for(int i=0;i<group->getNbElements();i++) {
- residuEl.setAsProxy(resd.getGroupProxy(k),i*nbFields,nbFields);
- KEl.setAsProxy(K[j]->getGroupProxy(k),i*nbFields,nbFields);
- iMassEl.setAsProxy(group->getInverseMassMatrix(),i*nbNodes,nbNodes);
- iMassEl.mult(residuEl,KEl);
- }
- Unp.getGroupProxy(k).add(K[j]->getGroupProxy(k),h*b[j]);
- }
- }
- sol.scale(0.);
- sol.axpy(Unp);
-
- for(int i=0;i<nStages;i++){
- delete K[i];
- }
- delete []K;
- }*/
-
-void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (maybe useless here)
- const dgConservationLaw &claw, // the conservation law
- dgGroupOfFaces &group,
- const fullMatrix<double> &solution, // solution !! at faces nodes
- fullMatrix<double> &solutionLeft,
- fullMatrix<double> &residual // residual !! at faces nodes
- )
-{
- int nbFields = claw.getNbFields();
+void dgResidualBoundary::compute1Group(
+ dgGroupOfFaces &group,
+ const fullMatrix<double> &solution, // solution !! at faces nodes
+ fullMatrix<double> &solutionLeft,
+ fullMatrix<double> &residual // residual !! at faces nodes
+ )
+{
+ int nbFields = _claw.getNbFields();
int nbNodesLeft = group.getGroupLeft().getNbNodes();
- const dgBoundaryCondition *boundaryCondition = claw.getBoundaryCondition(group.getBoundaryTag());
+ const dgBoundaryCondition *boundaryCondition = _claw.getBoundaryCondition(group.getBoundaryTag());
// ----- 1 ---- get the solution at quadrature points
fullMatrix<double> solutionQP (group.getNbIntegrationPoints(),group.getNbElements() * nbFields);
group.getCollocationMatrix().mult(solution, solutionQP);
@@ -437,8 +314,8 @@ void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (mayb
// inviscid
dataCacheDouble *boundaryTerm = boundaryCondition->newBoundaryTerm(cacheMapLeft);
// viscous
- dataCacheDouble *diffusiveFluxLeft = claw.newDiffusiveFlux(cacheMapLeft);
- dataCacheDouble *maximumDiffusivityLeft = claw.newMaximumDiffusivity(cacheMapLeft);
+ dataCacheDouble *diffusiveFluxLeft = _claw.newDiffusiveFlux(cacheMapLeft);
+ dataCacheDouble *maximumDiffusivityLeft = _claw.newMaximumDiffusivity(cacheMapLeft);
dataCacheDouble *maximumDiffusivityOut = boundaryCondition->newMaximumDiffusivity(cacheMapLeft);
dataCacheDouble *neumann = boundaryCondition->newDiffusiveNeumannBC(cacheMapLeft);
dataCacheDouble *dirichlet = boundaryCondition->newDiffusiveDirichletBC(cacheMapLeft);
@@ -493,56 +370,32 @@ void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (mayb
residual.gemm(group.getRedistributionMatrix(),NormalFluxQP);
}
+void dgResidualBoundary::computeAndMap1Group(dgGroupOfFaces &faces, dgDofContainer &solution, dgDofContainer &residual)
+{
+ int _nbFields = _claw.getNbFields();
+ fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*_nbFields);
+ fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*_nbFields);
+ faces.mapToInterface(_nbFields, solution.getGroupProxy(&faces.getGroupLeft()), solution.getGroupProxy(&faces.getGroupRight()), solInterface);
+ compute1Group(faces,solInterface,solution.getGroupProxy(&faces.getGroupLeft()),residuInterface);
+ faces.mapFromInterface(_nbFields, residuInterface, residual.getGroupProxy(&faces.getGroupLeft()), residual.getGroupProxy(&faces.getGroupRight()));
+}
+
// works for any number of groups
-void dgAlgorithm::residual( const dgConservationLaw &claw,
- dgGroupCollection &groups,
- dgDofContainer &solution,
- dgDofContainer &residu)
+void dgAlgorithm::residual( const dgConservationLaw &claw, dgGroupCollection &groups, dgDofContainer &solution, dgDofContainer &residual)
{
solution.scatter();
- int nbFields=claw.getNbFields();
- //volume term
- for(size_t i=0;i<groups.getNbElementGroups() ; i++) {
- residu.getGroupProxy(i).scale(0);
- residualVolume(claw,*groups.getElementGroup(i),solution.getGroupProxy(i),residu.getGroupProxy(i));
- }
- // residu[0]->print("Volume");
- //interface term
- for(size_t i=0;i<groups.getNbFaceGroups() ; i++) {
- dgGroupOfFaces &faces = *groups.getFaceGroup(i);
- int iGroupLeft = -1, iGroupRight = -1;
- for(size_t j=0;j<groups.getNbElementGroups(); j++) {
- dgGroupOfElements *groupj = groups.getElementGroup(j);
- if (groupj == &faces.getGroupLeft()) iGroupLeft = j;
- if (groupj == &faces.getGroupRight()) iGroupRight= j;
- }
- for(size_t j=0;j<groups.getNbGhostGroups(); j++) {
- dgGroupOfElements *groupj = groups.getGhostGroup(j);
- if (groupj == &faces.getGroupLeft()) iGroupLeft = j;
- if (groupj == &faces.getGroupLeft())iGroupLeft = j + groups.getNbElementGroups();
- if (groupj == &faces.getGroupRight())iGroupRight= j + groups.getNbElementGroups();
- }
- fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*2*nbFields);
- fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*2*nbFields);
- faces.mapToInterface(nbFields, solution.getGroupProxy(iGroupLeft), solution.getGroupProxy(iGroupRight), solInterface);
- residualInterface(claw,faces,solInterface,solution.getGroupProxy(iGroupLeft), solution.getGroupProxy(iGroupRight),residuInterface);
- faces.mapFromInterface(nbFields, residuInterface,residu.getGroupProxy(iGroupLeft), residu.getGroupProxy(iGroupRight));
- }
- //boundaries
- for(size_t i=0;i<groups.getNbBoundaryGroups() ; i++) {
- dgGroupOfFaces &faces = *groups.getBoundaryGroup(i);
- int iGroupLeft = -1, iGroupRight = -1;
- for(size_t j=0;j<groups.getNbElementGroups(); j++) {
- dgGroupOfElements *groupj = groups.getElementGroup(j);
- if (groupj == &faces.getGroupLeft())iGroupLeft = j;
- if (groupj == &faces.getGroupRight())iGroupRight= j;
- }
- fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*nbFields);
- fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*nbFields);
- faces.mapToInterface(nbFields, solution.getGroupProxy(iGroupLeft), solution.getGroupProxy(iGroupRight), solInterface);
- residualBoundary(claw,faces,solInterface,solution.getGroupProxy(iGroupLeft),residuInterface);
- faces.mapFromInterface(nbFields, residuInterface, residu.getGroupProxy(iGroupLeft), residu.getGroupProxy(iGroupRight));
- }
+
+ dgResidualVolume residualVolume(claw);
+ for (int i=0; i<groups.getNbElementGroups(); i++)
+ residualVolume.computeAndMap1Group(*groups.getElementGroup(i), solution, residual);
+
+ dgResidualInterface residualInterface(claw);
+ for(size_t i=0;i<groups.getNbFaceGroups() ; i++)
+ residualInterface.computeAndMap1Group(*groups.getFaceGroup(i), solution, residual);
+
+ dgResidualBoundary residualBoundary(claw);
+ for(size_t i=0;i<groups.getNbBoundaryGroups() ; i++)
+ residualBoundary.computeAndMap1Group(*groups.getBoundaryGroup(i), solution, residual);
}
void dgAlgorithm::computeElementaryTimeSteps ( //dofManager &dof, // the DOF manager (maybe useless here)
@@ -590,3 +443,19 @@ void dgAlgorithm::computeElementaryTimeSteps ( //dofManager &dof, // the DOF man
}
}
+void dgAlgorithm::multAddInverseMassMatrix ( /*dofManager &dof,*/
+ const dgGroupOfElements & group,
+ fullMatrix<double> &residu,
+ fullMatrix<double> &sol)
+{
+ fullMatrix<double> residuEl;
+ fullMatrix<double> solEl;
+ fullMatrix<double> iMassEl;
+ int nFields=sol.size2()/group.getNbElements();
+ for(int i=0;i<group.getNbElements();i++) {
+ residuEl.setAsProxy(residu,i*nFields,nFields);
+ solEl.setAsProxy(sol,i*nFields,nFields);
+ iMassEl.setAsProxy(group.getInverseMassMatrix(),i*group.getNbNodes(),group.getNbNodes());
+ solEl.gemm(iMassEl,residuEl);
+ }
+}
diff --git a/Solver/dgAlgorithm.h b/Solver/dgAlgorithm.h
index f970010eb839eda4a325399ef3e73ccd79e40aa1..6032f9a639f8c7323fa549fcd06ed3e5ff21ca6e 100644
--- a/Solver/dgAlgorithm.h
+++ b/Solver/dgAlgorithm.h
@@ -16,18 +16,6 @@ class dgSystemOfEquations;
class dgAlgorithm {
public :
static void registerBindings(binding *b);
- static void residualVolume ( //dofManager &dof, // the DOF manager (maybe useless here)
- const dgConservationLaw &claw, // the conservation law
- const dgGroupOfElements & group,
- fullMatrix<double> &solution,
- fullMatrix<double> &residual);
- static void residualInterface ( //dofManager &dof, // the DOF manager (maybe useless here)
- const dgConservationLaw &claw, // the conservation law
- dgGroupOfFaces & group,
- const fullMatrix<double> &solution, // ! at face nodes
- fullMatrix<double> &solutionRight,
- fullMatrix<double> &solutionLeft,
- fullMatrix<double> &residual); // ! at face nodes
static void residualBoundary ( //dofManager &dof, // the DOF manager (maybe useless here)
const dgConservationLaw &claw, // the conservation law
dgGroupOfFaces &group,
diff --git a/Solver/dgConservationLawPerfectGas.cpp b/Solver/dgConservationLawPerfectGas.cpp
index 1622dc7c81e69fbc787ba6e50d5a6ea44ca828a4..7dab96977826844e784de4c7e960d4854e1ca0d8 100644
--- a/Solver/dgConservationLawPerfectGas.cpp
+++ b/Solver/dgConservationLawPerfectGas.cpp
@@ -332,7 +332,7 @@ class dgPerfectGasLaw2d::advection : public dataCacheDouble {
sol(cacheMap.get("Solution",this))
{};
void _eval () {
- const int nQP = sol().size1();
+ const int nQP = _value.size1();
const double GM1 = GAMMA - 1.0;
for (size_t k = 0 ; k < nQP; k++ ){
const double invrho = 1./sol(k,0);
diff --git a/Solver/function.cpp b/Solver/function.cpp
index 62b88d983d56d8d9752ef5cb9fd1e96bb29a02d9..2fab642de27a89506e20cc5ed7477395b289b86b 100644
--- a/Solver/function.cpp
+++ b/Solver/function.cpp
@@ -282,8 +282,10 @@ public:
};
void dataCacheMap::setNbEvaluationPoints(int nbEvaluationPoints) {
_nbEvaluationPoints = nbEvaluationPoints;
- for(std::map<std::string, dataCacheDouble*>::iterator it = _cacheDoubleMap.begin(); it!= _cacheDoubleMap.end(); it++)
- it->second->resize();
+ for(std::set<dataCacheDouble*>::iterator it = _toResize.begin(); it!= _toResize.end(); it++){
+ (*it)->resize();
+ (*it)->_valid = false;
+ }
}
dataCacheDouble *functionMesh2Mesh::newDataCache(dataCacheMap *m)
@@ -295,6 +297,7 @@ dataCacheDouble *functionMesh2Mesh::newDataCache(dataCacheMap *m)
dataCacheDouble::dataCacheDouble(dataCacheMap &map,int nRowByPoint, int nCol):
dataCache(&map),_cacheMap(map),_value(nRowByPoint==0?1:nRowByPoint*map.getNbEvaluationPoints(),nCol){
_nRowByPoint=nRowByPoint;
+ map.addDataCacheDouble(this);
};
void dataCacheDouble::resize() {
_value = fullMatrix<double>(_nRowByPoint==0?1:_nRowByPoint*_cacheMap.getNbEvaluationPoints(),_value.size2());
diff --git a/Solver/function.h b/Solver/function.h
index 607135a53a78b00e529b3e555dd9a68db5488755..bdfd9c64da7a4447e6402492021ad302549fa5cd 100644
--- a/Solver/function.h
+++ b/Solver/function.h
@@ -159,6 +159,7 @@ class dataCacheElement : public dataCache {
// more explanation at the head of this file
class dataCacheMap {
friend class dataCache;
+ friend class dataCacheDouble;
private:
int _nbEvaluationPoints;
// keep track of the current element and all the dataCaches that
@@ -176,11 +177,15 @@ class dataCacheMap {
}
};
std::set<dataCache*> _toDelete;
+ std::set<dataCacheDouble*> _toResize;
protected:
void addDataCache(dataCache *data){
_toDelete.insert(data);
}
+ void addDataCacheDouble(dataCacheDouble *data){
+ _toResize.insert(data);
+ }
public:
// dg Part
//list of dgDofContainer that have to be evaluated, this not work as a function as it depends on the algorithm