Skip to content
Snippets Groups Projects
Commit b265aa84 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

dg :keep the same datacacheMap for all groups (except boundaries)

parent 369ed8b4
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Solver/dgAlgorithm.cpp
+ 210
− 341
View file @ b265aa84
......@@ -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 );
_solutionQPe.setAsProxy(solutionQP, iElement*_nbFields, _nbFields );
if(gradientSolutionQPe.somethingDependOnMe()){
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)
{
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);
}
}
/*void dgAlgorithm::multirateRungeKutta (const dgConservationLaw &claw, // conservation law
dgGroupCollection &groups,
double h, // time-step
dgDofContainer &sol,
dgDofContainer &resd,
int orderRK) // order of RK integrator
void dgResidualInterface::computeAndMap1Group (dgGroupOfFaces &faces, dgDofContainer &solution, dgDofContainer &residual)
{
// 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;
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()));
}
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);
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);
};
for(int i=0;i<nStages;i++){
delete K[i];
dgResidualBoundary::dgResidualBoundary(const dgConservationLaw &claw):_claw(claw){
}
delete []K;
}*/
void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (maybe useless here)
const dgConservationLaw &claw, // the conservation law
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 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);
}
}
Solver/dgAlgorithm.h
+ 0
− 12
View file @ b265aa84
......@@ -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,
......
Solver/dgConservationLawPerfectGas.cpp
+ 1
− 1
View file @ b265aa84
......@@ -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);
......
Solver/function.cpp
+ 5
− 2
View file @ b265aa84
......@@ -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());
......
Solver/function.h
+ 5
− 0
View file @ b265aa84
......@@ -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
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment