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

dg : merge residualBoundary into residualInterface : same code for boundaries,... 

dg : merge residualBoundary into residualInterface : same code for boundaries, regular 2-sided-faces and bifurcations
parent 8c7c06c1
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 63 additions and 207 deletions
Solver/TESTCASES/AdvectionDiffusion.lua
+ 1
− 1
View file @ e9522f01
...@@ -3,7 +3,7 @@ model:load ('square.geo') ...@@ -3,7 +3,7 @@ model:load ('square.geo')
model:load ('square.msh') model:load ('square.msh')
dg = dgSystemOfEquations (model) dg = dgSystemOfEquations (model)
dg:setOrder(4) dg:setOrder(1)
-- conservation law -- conservation law
law = dgConservationLawAdvectionDiffusion(functionConstant({0.15,0.05,0}),functionConstant({0.001})) law = dgConservationLawAdvectionDiffusion(functionConstant({0.15,0.05,0}),functionConstant({0.001}))
......
Solver/TESTCASES/Diffusion.lua
+ 7
− 10
View file @ e9522f01
...@@ -7,15 +7,12 @@ dg:setOrder(1) ...@@ -7,15 +7,12 @@ dg:setOrder(1)
-- conservation law -- conservation law
-- advection speed -- advection speed
nu=fullMatrix(1,1); law = dgConservationLawAdvectionDiffusion.diffusionLaw(functionConstant({0.01}))
nu:set(0,0,0.01)
law = dgConservationLawAdvectionDiffusion('',functionConstant(nu):getName())
dg:setConservationLaw(law) dg:setConservationLaw(law)
-- boundary condition -- boundary condition
outside=fullMatrix(1,1) law:addBoundaryCondition('Border',law:newOutsideValueBoundary(functionConstant({1})))
outside:set(0,0,0.) --law:addBoundaryCondition('Border',law:new0FluxBoundary())
law:addBoundaryCondition('Border',law:new0FluxBoundary())
dg:setup() dg:setup()
...@@ -28,14 +25,14 @@ function initial_condition( xyz , f ) ...@@ -28,14 +25,14 @@ function initial_condition( xyz , f )
f:set (i, 0, math.exp(-100*((x-0.2)^2 +(y-0.3)^2))) f:set (i, 0, math.exp(-100*((x-0.2)^2 +(y-0.3)^2)))
end end
end end
dg:L2Projection(functionLua(1,'initial_condition',{'XYZ'}):getName()) -- dg:L2Projection(functionLua(1,'initial_condition',{functionCoordinates.get()}))
dg:exportSolution('output/Diffusion_00000') dg:exportSolution('output/Diffusion_00000')
-- main loop -- main loop
for i=1,10000 do for i=1,100 do
--norm = dg:RK44(0.001) norm = dg:RK44(0.01)
norm = dg:RK44_TransformNodalValue(0.00001) --norm = dg:RK44_TransformNodalValue(0.00001)
if (i % 50 == 0) then if (i % 50 == 0) then
print('iter',i,norm) print('iter',i,norm)
dg:exportSolution(string.format("output/Diffusion-%05d", i)) dg:exportSolution(string.format("output/Diffusion-%05d", i))
......
Solver/TESTCASES/supra.lua
+ 5
− 7
View file @ e9522f01
...@@ -28,10 +28,8 @@ function valueLeft(f) ...@@ -28,10 +28,8 @@ function valueLeft(f)
end end
end end
law = dgConservationLawAdvectionDiffusion('',functionLua(1,'diffusivity',{'Solution'}):getName()) law = dgConservationLawAdvectionDiffusion.diffusionLaw(functionLua(1,'diffusivity',{functionSolution.get()}))
print(law)
dg:setConservationLaw(law) dg:setConservationLaw(law)
print(law)
-- boundary condition -- boundary condition
--[[ --[[
...@@ -42,8 +40,8 @@ outsideRight:set(0,0,1) ...@@ -42,8 +40,8 @@ outsideRight:set(0,0,1)
--]] --]]
law:addBoundaryCondition('Top',law:new0FluxBoundary()) law:addBoundaryCondition('Top',law:new0FluxBoundary())
law:addBoundaryCondition('Bottom',law:new0FluxBoundary()) law:addBoundaryCondition('Bottom',law:new0FluxBoundary())
leftFunction=functionLua(1,'valueLeft',{}):getName() leftFunction=functionLua(1,'valueLeft',{})
rightFunction=functionLua(1,'valueRight',{}):getName() rightFunction=functionLua(1,'valueRight',{})
--[[ --[[
law:addBoundaryCondition('Left',law:newNeumannBoundary(leftFunction)) law:addBoundaryCondition('Left',law:newNeumannBoundary(leftFunction))
law:addBoundaryCondition('Right',law:newNeumannBoundary(rightFunction)) law:addBoundaryCondition('Right',law:newNeumannBoundary(rightFunction))
...@@ -62,7 +60,7 @@ function initial_condition( xyz , f ) ...@@ -62,7 +60,7 @@ function initial_condition( xyz , f )
f:set (i, 0,0*math.exp(-100*((x-0.2)^2 +(y-0.3)^2))) f:set (i, 0,0*math.exp(-100*((x-0.2)^2 +(y-0.3)^2)))
end end
end end
dg:L2Projection(functionLua(1,'initial_condition',{'XYZ'}):getName()) dg:L2Projection(functionLua(1,'initial_condition',{functionCoordinates.get()}))
dg:exportSolution('output/supra-00000') dg:exportSolution('output/supra-00000')
...@@ -70,7 +68,7 @@ dg:exportSolution('output/supra-00000') ...@@ -70,7 +68,7 @@ dg:exportSolution('output/supra-00000')
for i=1,150000 do for i=1,150000 do
norm = dg:RK44_limiter(dt) norm = dg:RK44_limiter(dt)
t=t+dt t=t+dt
if (i % 100 == 0) then if (i % 1 == 0) then
print('iter',i,norm) print('iter',i,norm)
dg:exportSolution(string.format("output/supra-%05d", i)) dg:exportSolution(string.format("output/supra-%05d", i))
end end
......
Solver/dgGroupOfElements.cpp
+ 0
− 5
View file @ e9522f01
...@@ -1236,15 +1236,10 @@ void dgGroupCollection::registerBindings(binding *b) ...@@ -1236,15 +1236,10 @@ void dgGroupCollection::registerBindings(binding *b)
cm->setDescription("return the number of dgGroupOfElements"); cm->setDescription("return the number of dgGroupOfElements");
cm = cb->addMethod("getNbFaceGroups", &dgGroupCollection::getNbFaceGroups); cm = cb->addMethod("getNbFaceGroups", &dgGroupCollection::getNbFaceGroups);
cm->setDescription("return the number of dgGroupOfFaces (interior ones, not the domain boundaries)"); cm->setDescription("return the number of dgGroupOfFaces (interior ones, not the domain boundaries)");
cm = cb->addMethod("getNbBoundaryGroups", &dgGroupCollection::getNbBoundaryGroups);
cm->setDescription("return the number of group of boundary faces");
cm = cb->addMethod("getElementGroup", &dgGroupCollection::getElementGroup); cm = cb->addMethod("getElementGroup", &dgGroupCollection::getElementGroup);
cm->setDescription("get 1 group of elements"); cm->setDescription("get 1 group of elements");
cm->setArgNames("id", NULL); cm->setArgNames("id", NULL);
cm = cb->addMethod("getFaceGroup", &dgGroupCollection::getFaceGroup); cm = cb->addMethod("getFaceGroup", &dgGroupCollection::getFaceGroup);
cm->setDescription("get 1 group of faces"); cm->setDescription("get 1 group of faces");
cm->setArgNames("id", NULL); cm->setArgNames("id", NULL);
cm = cb->addMethod("getBoundaryGroup", &dgGroupCollection::getBoundaryGroup);
cm->setDescription("get 1 group of boundary faces");
cm->setArgNames("id", NULL);
} }
Solver/dgGroupOfElements.h
+ 0
− 2
View file @ e9522f01
...@@ -236,11 +236,9 @@ class dgGroupCollection { ...@@ -236,11 +236,9 @@ class dgGroupCollection {
inline GModel* getModel() {return _model;} inline GModel* getModel() {return _model;}
inline int getNbElementGroups() const {return _elementGroups.size();} inline int getNbElementGroups() const {return _elementGroups.size();}
inline int getNbFaceGroups() const {return _faceGroups.size();} inline int getNbFaceGroups() const {return _faceGroups.size();}
inline int getNbBoundaryGroups() const {return _boundaryGroups.size();}
inline int getNbGhostGroups() const {return _ghostGroups.size();} inline int getNbGhostGroups() const {return _ghostGroups.size();}
inline dgGroupOfElements *getElementGroup(int i) const {return i<getNbElementGroups()?_elementGroups[i]:_ghostGroups[i-getNbElementGroups()];} inline dgGroupOfElements *getElementGroup(int i) const {return i<getNbElementGroups()?_elementGroups[i]:_ghostGroups[i-getNbElementGroups()];}
inline dgGroupOfFaces *getFaceGroup(int i) const {return _faceGroups[i];} inline dgGroupOfFaces *getFaceGroup(int i) const {return _faceGroups[i];}
inline dgGroupOfFaces *getBoundaryGroup(int i) const {return _boundaryGroups[i];}
inline dgGroupOfElements *getGhostGroup(int i) const {return _ghostGroups[i];} inline dgGroupOfElements *getGhostGroup(int i) const {return _ghostGroups[i];}
inline int getNbImageElementsOnPartition(int partId) const {return _elementsToSend[partId].size();} inline int getNbImageElementsOnPartition(int partId) const {return _elementsToSend[partId].size();}
......
Solver/dgLimiter.cpp
+ 2
− 0
View file @ e9522f01
...@@ -22,6 +22,8 @@ int dgSlopeLimiter::apply ( dgDofContainer *solution) ...@@ -22,6 +22,8 @@ int dgSlopeLimiter::apply ( dgDofContainer *solution)
fullMatrix<double> TempL, TempR; fullMatrix<double> TempL, TempR;
for( int iGFace=0; iGFace<groups->getNbFaceGroups(); iGFace++) { for( int iGFace=0; iGFace<groups->getNbFaceGroups(); iGFace++) {
dgGroupOfFaces* group = groups->getFaceGroup(iGFace); dgGroupOfFaces* group = groups->getFaceGroup(iGFace);
if (group->getNbGroupOfConnections()!=2) // skip boundaries and bifurcations
continue;
const dgGroupOfConnections &left = group->getGroupOfConnections(0); const dgGroupOfConnections &left = group->getGroupOfConnections(0);
const dgGroupOfConnections &right = group->getGroupOfConnections(1); const dgGroupOfConnections &right = group->getGroupOfConnections(1);
const dgGroupOfElements *groupLeft = &left.getGroupOfElements(); const dgGroupOfElements *groupLeft = &left.getGroupOfElements();
......
Solver/dgResidual.cpp
+ 39
− 140
View file @ e9522f01
...@@ -292,7 +292,16 @@ void dgResidualInterface::compute1Group ( //dofManager &dof, // the DOF manager ...@@ -292,7 +292,16 @@ void dgResidualInterface::compute1Group ( //dofManager &dof, // the DOF manager
maximumDiffusivities[i] = _claw.newMaximumDiffusivity(caches[i]); maximumDiffusivities[i] = _claw.newMaximumDiffusivity(caches[i]);
caches[i].setNbEvaluationPoints(group.getNbIntegrationPoints()); caches[i].setNbEvaluationPoints(group.getNbIntegrationPoints());
} }
dataCacheDouble *riemannSolver = NULL; dataCacheDouble *riemannSolver = NULL, *neumann = NULL, *dirichlet = NULL, *maximumDiffusivityOut = NULL;
if (nbConnections == 1) {
const dgBoundaryCondition *boundaryCondition = _claw.getBoundaryCondition(group.getPhysicalTag());
riemannSolver = boundaryCondition->newBoundaryTerm(caches[0]);
neumann = boundaryCondition->newDiffusiveNeumannBC(caches[0]);
dirichlet = boundaryCondition->newDiffusiveDirichletBC(caches[0]);
maximumDiffusivityOut = boundaryCondition->newMaximumDiffusivity(caches[0]);
if (!maximumDiffusivityOut)
maximumDiffusivityOut = maximumDiffusivities[0];
}
if (nbConnections == 2) if (nbConnections == 2)
riemannSolver = _claw.newRiemannSolver(caches[0], caches[1]); riemannSolver = _claw.newRiemannSolver(caches[0], caches[1]);
if (nbConnections == 3) if (nbConnections == 3)
...@@ -324,13 +333,14 @@ void dgResidualInterface::compute1Group ( //dofManager &dof, // the DOF manager ...@@ -324,13 +333,14 @@ void dgResidualInterface::compute1Group ( //dofManager &dof, // the DOF manager
// proxies for the flux // proxies for the flux
normalFluxQP.setAsProxy(NormalFluxQP, iFace*_nbFields*nbConnections, _nbFields*nbConnections); normalFluxQP.setAsProxy(NormalFluxQP, iFace*_nbFields*nbConnections, _nbFields*nbConnections);
// B3 ) compute fluxes // B3 ) compute fluxes
if (diffusiveFluxes[0] && nbConnections==2) { // IP only implemented for 2-sides faces (not for bifurcations) if (diffusiveFluxes[0]) { // IP only implemented for 2-sides faces (not for bifurcations)
for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) { if(nbConnections == 2) {
const double detJ = group.getDetJ (iFace, iPt);
//just for the lisibility : //just for the lisibility :
const fullMatrix<double> &dfl = (*diffusiveFluxes[0])(); const fullMatrix<double> &dfl = (*diffusiveFluxes[0])();
const fullMatrix<double> &dfr = (*diffusiveFluxes[1])(); const fullMatrix<double> &dfr = (*diffusiveFluxes[1])();
const fullMatrix<double> &n = (caches[0].getNormals(NULL))(); const fullMatrix<double> &n = (caches[0].getNormals(NULL))();
for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
const double detJ = group.getDetJ (iFace, iPt);
for (int k=0;k<_nbFields;k++) { for (int k=0;k<_nbFields;k++) {
double meanNormalFlux = double meanNormalFlux =
((dfl(iPt,k+_nbFields*0)+dfr(iPt,k+_nbFields*0)) * n(0,iPt) ((dfl(iPt,k+_nbFields*0)+dfr(iPt,k+_nbFields*0)) * n(0,iPt)
...@@ -346,6 +356,18 @@ void dgResidualInterface::compute1Group ( //dofManager &dof, // the DOF manager ...@@ -346,6 +356,18 @@ void dgResidualInterface::compute1Group ( //dofManager &dof, // the DOF manager
normalFluxQP(iPt,k+_nbFields) += (meanNormalFlux+solutionJumpPenalty)*detJ; normalFluxQP(iPt,k+_nbFields) += (meanNormalFlux+solutionJumpPenalty)*detJ;
} }
} }
} else if(nbConnections == 1) {
for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
const double detJ = group.getDetJ (iFace, iPt);
for (int k=0;k<_nbFields;k++) {
double minl = elementGroups[0]->getElementVolume(connections[0]->getElementId(iFace))/group.getInterfaceSurface(iFace);
double nu = std::max((*maximumDiffusivities[0])(iPt,0),(*maximumDiffusivityOut)(iPt,0));
double mu = (p+1)*(p+dim)/dim*nu/minl;
double solutionJumpPenalty = (caches[0].getSolution(NULL)(iPt,k)-(*dirichlet)(iPt,k))/2*mu;
normalFluxQP(iPt,k) -= ((*neumann)(iPt,k)+solutionJumpPenalty)*detJ;
}
}
}
} }
if (riemannSolver) { if (riemannSolver) {
for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) { for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
...@@ -354,6 +376,7 @@ void dgResidualInterface::compute1Group ( //dofManager &dof, // the DOF manager ...@@ -354,6 +376,7 @@ void dgResidualInterface::compute1Group ( //dofManager &dof, // the DOF manager
normalFluxQP(iPt,k) += (*riemannSolver)(iPt,k)*detJ; normalFluxQP(iPt,k) += (*riemannSolver)(iPt,k)*detJ;
} }
} }
} }
// C ) redistribute the flux to the residual (at Faces nodes) // C ) redistribute the flux to the residual (at Faces nodes)
if(riemannSolver || diffusiveFluxes[0]) if(riemannSolver || diffusiveFluxes[0])
...@@ -384,118 +407,6 @@ dgResidualInterface::~dgResidualInterface () ...@@ -384,118 +407,6 @@ dgResidualInterface::~dgResidualInterface ()
{ {
} }
dgResidualBoundary::dgResidualBoundary(const dgConservationLaw &claw):_claw(claw){
}
void dgResidualBoundary::compute1Group(
dgGroupOfFaces &group,
const fullMatrix<double> &solution, // solution !! at faces nodes
const fullMatrix<double> &solutionLeft,
fullMatrix<double> &residual // residual !! at faces nodes
)
{
//should be splitted like dgResidualInterface and dgResidualVolume
//but i do not do it know because dgResidualBoundary will probably disapear when we will have list of elements on interfaces
const dgGroupOfConnections &left = group.getGroupOfConnections(0);
const dgGroupOfElements &groupLeft = left.getGroupOfElements();
int nbFields = _claw.getNbFields();
const dgBoundaryCondition *boundaryCondition = _claw.getBoundaryCondition(group.getPhysicalTag());
// ----- 1 ---- get the solution at quadrature points
fullMatrix<double> solutionQP (group.getNbIntegrationPoints(),group.getNbElements() * nbFields);
group.getCollocationMatrix().mult(solution, solutionQP);
// ----- 2 ---- compute normal fluxes at integration points
fullMatrix<double> NormalFluxQP ( group.getNbIntegrationPoints(), group.getNbElements()*nbFields);
dataCacheMap cacheMapLeft;
cacheMapLeft.setNbEvaluationPoints(group.getNbIntegrationPoints());
const fullMatrix<double> &DPsiLeftDx = left.getDPsiDx();
// provided dataCache
dataCacheDouble &uvw=cacheMapLeft.provideParametricCoordinates();
dataCacheDouble &solutionQPLeft = cacheMapLeft.provideSolution(nbFields);
dataCacheDouble &gradientSolutionLeft = cacheMapLeft.provideSolutionGradient(nbFields);
dataCacheDouble &normals = cacheMapLeft.provideNormals();
dataCacheElement &cacheElementLeft = cacheMapLeft.getElement();
// required data
// inviscid
dataCacheDouble *boundaryTerm = boundaryCondition->newBoundaryTerm(cacheMapLeft);
// viscous
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);
fullMatrix<double> normalFluxQP,dPsiLeftDx,dofsLeft;
int p = groupLeft.getOrder();
int dim = left.getElement(0)->getDim();
for (int iFace=0 ; iFace<group.getNbElements() ;++iFace) {
normalFluxQP.setAsProxy(NormalFluxQP, iFace*nbFields, nbFields);
// ----- 2.3.1 --- provide the data to the cacheMap
solutionQPLeft.setAsProxy(solutionQP, iFace*nbFields, nbFields);
normals.setAsProxy(left.getNormals(),iFace*group.getNbIntegrationPoints(),group.getNbIntegrationPoints());
// proxies needed to compute the gradient of the solution and the IP term
dPsiLeftDx.setAsProxy(DPsiLeftDx,iFace*groupLeft.getNbNodes(), groupLeft.getNbNodes());
dofsLeft.setAsProxy(solutionLeft, nbFields*left.getElementId(iFace), nbFields);
uvw.setAsProxy(left.getIntegrationPointsOnElement(iFace));
cacheElementLeft.set(left.getElement(iFace));
// compute the gradient of the solution
if(gradientSolutionLeft.somethingDependOnMe()){
dPsiLeftDx.mult(dofsLeft, gradientSolutionLeft.set());
}
// ----- 2.3.2 --- compute inviscid contribution
for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
const double detJ = group.getDetJ (iFace, iPt);
for (int k=0;k<nbFields;k++)
normalFluxQP(iPt,k) = (*boundaryTerm)(iPt,k)*detJ;
}
// ----- 2.3.3 --- compute viscous contribution
if (diffusiveFluxLeft) {
for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
const double detJ = group.getDetJ (iFace, iPt);
//just for the lisibility :
for (int k=0;k<nbFields;k++) {
double minl = groupLeft.getElementVolume(left.getElementId(iFace))/group.getInterfaceSurface(iFace);
double nu = (*maximumDiffusivityLeft)(iPt,0);
if(maximumDiffusivityOut)
nu = std::max(nu,(*maximumDiffusivityOut)(iPt,0));
double mu = (p+1)*(p+dim)/dim*nu/minl;
double solutionJumpPenalty = (solutionQPLeft(iPt,k)-(*dirichlet)(iPt,k))/2*mu;
normalFluxQP(iPt,k) -= ((*neumann)(iPt,k)+solutionJumpPenalty)*detJ;
}
}
}
}
// ----- 3 ---- do the redistribution at face nodes using BLAS3
residual.gemm(group.getRedistributionMatrix(),NormalFluxQP);
}
void dgResidualBoundary::computeAndMap1Group(dgGroupOfFaces &faces, dgDofContainer &solution, dgDofContainer &residual)
{
int _nbFields = _claw.getNbFields();
int nbConnections = faces.getNbGroupOfConnections();
fullMatrix<double> solInterface(faces.getNbNodes(), faces.getNbElements()*_nbFields*nbConnections);
fullMatrix<double> residuInterface(faces.getNbNodes(), faces.getNbElements()*_nbFields*nbConnections);
std::vector<const dgGroupOfElements *> elGroups(nbConnections);
std::vector<const fullMatrix<double> *> solutionProxies(nbConnections);
std::vector<fullMatrix<double> *> residualProxies(nbConnections);
for (int i=0; i<nbConnections; i++) {
elGroups[i] = &faces.getGroupOfConnections(i).getGroupOfElements();
solutionProxies[i] = &solution.getGroupProxy(elGroups[i]);
residualProxies[i] = &residual.getGroupProxy(elGroups[i]);
}
faces.mapToInterface(_nbFields, solutionProxies, solInterface);
compute1Group(faces, solInterface, *solutionProxies[0],residuInterface);
faces.mapFromInterface(_nbFields, residuInterface, residualProxies);
}
void dgResidual::compute(dgGroupCollection &groups, dgDofContainer &solution, dgDofContainer &residual) void dgResidual::compute(dgGroupCollection &groups, dgDofContainer &solution, dgDofContainer &residual)
{ {
solution.scatter(); solution.scatter();
...@@ -506,9 +417,6 @@ void dgResidual::compute(dgGroupCollection &groups, dgDofContainer &solution, dg ...@@ -506,9 +417,6 @@ void dgResidual::compute(dgGroupCollection &groups, dgDofContainer &solution, dg
for(size_t i=0;i<groups.getNbFaceGroups() ; i++){ for(size_t i=0;i<groups.getNbFaceGroups() ; i++){
residualInterface.computeAndMap1Group(*groups.getFaceGroup(i), solution, residual); 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);
} }
#include "Bindings.h" #include "Bindings.h"
...@@ -545,13 +453,4 @@ void dgResidual::registerBindings (binding *b) ...@@ -545,13 +453,4 @@ void dgResidual::registerBindings (binding *b)
// mb = cb->addMethod("compute1Group", &dgResidualInterface::compute1Group); // mb = cb->addMethod("compute1Group", &dgResidualInterface::compute1Group);
//mb->setDescription("compute the residual for one group given fullMatrices with the solution at faces nodes and at element nodes"); //mb->setDescription("compute the residual for one group given fullMatrices with the solution at faces nodes and at element nodes");
//mb->setArgNames("group", "solutionFaces", "solutionOnElements", "residual", NULL); //mb->setArgNames("group", "solutionFaces", "solutionOnElements", "residual", NULL);
cb = b->addClass<dgResidualBoundary>("dgResidualBoundary");
cb->setDescription("compute the residual for one boundary dgGroupOfFaces");
mb = cb->addMethod("computeAndMap1Group",&dgResidualBoundary::computeAndMap1Group);
mb->setDescription("compute the residual for one group given a dgDofContainer solution");
mb->setArgNames("group", "solution", "residual", NULL);
mb = cb->addMethod("compute1Group", &dgResidualBoundary::compute1Group);
mb->setDescription("compute the residual for one group given fullMatrices with the solution at faces nodes and at element nodes");
mb->setArgNames("group", "solutionFaces", "solutionGroupLeft", "residual", NULL);
} }
Solver/dgResidual.h
+ 0
− 13
View file @ e9522f01
...@@ -42,19 +42,6 @@ class dgResidualInterface { ...@@ -42,19 +42,6 @@ class dgResidualInterface {
~dgResidualInterface(); ~dgResidualInterface();
}; };
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
const fullMatrix<double> &solutionLeft,
fullMatrix<double> &residual // residual !! at faces nodes
);
void computeAndMap1Group (dgGroupOfFaces &faces, dgDofContainer &solution, dgDofContainer &residual);
dgResidualBoundary (const dgConservationLaw &claw);
};
class dgResidual { class dgResidual {
const dgConservationLaw &_claw; const dgConservationLaw &_claw;
public: public:
......
Solver/dgRungeKuttaMultirate.cpp
+ 9
− 28
View file @ e9522f01
...@@ -81,7 +81,6 @@ dgRungeKuttaMultirate::dgRungeKuttaMultirate(dgGroupCollection* gc,dgConservatio ...@@ -81,7 +81,6 @@ dgRungeKuttaMultirate::dgRungeKuttaMultirate(dgGroupCollection* gc,dgConservatio
_law=law; _law=law;
_residualVolume=new dgResidualVolume(*_law); _residualVolume=new dgResidualVolume(*_law);
_residualInterface=new dgResidualInterface(*_law); _residualInterface=new dgResidualInterface(*_law);
_residualBoundary=new dgResidualBoundary(*_law);
_gc=gc; _gc=gc;
_K=new dgDofContainer*[nStages]; _K=new dgDofContainer*[nStages];
for(int i=0;i<nStages;i++){ for(int i=0;i<nStages;i++){
...@@ -113,32 +112,15 @@ dgRungeKuttaMultirate::dgRungeKuttaMultirate(dgGroupCollection* gc,dgConservatio ...@@ -113,32 +112,15 @@ dgRungeKuttaMultirate::dgRungeKuttaMultirate(dgGroupCollection* gc,dgConservatio
} }
for(int iGroup=0;iGroup<gc->getNbFaceGroups();iGroup++){ for(int iGroup=0;iGroup<gc->getNbFaceGroups();iGroup++){
dgGroupOfFaces *gf=gc->getFaceGroup(iGroup); dgGroupOfFaces *gf=gc->getFaceGroup(iGroup);
const dgGroupOfElements *ge[2]; for(int i=0;i<gf->getNbGroupOfConnections();i++){
ge[0]=&gf->getGroupOfConnections(0).getGroupOfElements(); const dgGroupOfElements *ge = &gf->getGroupOfConnections(0).getGroupOfElements();
ge[1]=&gf->getGroupOfConnections(1).getGroupOfElements(); if(ge->getIsInnerMultirateBuffer()){
for(int i=0;i<2;i++){ _innerBufferGroupsOfElements[ge->getMultirateExponent()].second.push_back(gf);
if(ge[i]->getIsInnerMultirateBuffer()){
_innerBufferGroupsOfElements[ge[i]->getMultirateExponent()].second.push_back(gf);
}
else if(ge[i]->getIsOuterMultirateBuffer()){
_outerBufferGroupsOfElements[ge[i]->getMultirateExponent()].second.push_back(gf);
}else{
_bulkGroupsOfElements[ge[i]->getMultirateExponent()].second.push_back(gf);
}
}
}
for(int iGroup=0;iGroup<gc->getNbBoundaryGroups();iGroup++){
dgGroupOfFaces *gf=gc->getBoundaryGroup(iGroup);
const dgGroupOfElements *ge[1];
ge[0]=&gf->getGroupOfConnections(0).getGroupOfElements();
for(int i=0;i<1;i++){
if(ge[i]->getIsInnerMultirateBuffer()){
_innerBufferGroupsOfElements[ge[i]->getMultirateExponent()].second.push_back(gf);
} }
else if(ge[i]->getIsOuterMultirateBuffer()){ else if(ge->getIsOuterMultirateBuffer()){
_outerBufferGroupsOfElements[ge[i]->getMultirateExponent()].second.push_back(gf); _outerBufferGroupsOfElements[ge->getMultirateExponent()].second.push_back(gf);
}else{ }else{
_bulkGroupsOfElements[ge[i]->getMultirateExponent()].second.push_back(gf); _bulkGroupsOfElements[ge->getMultirateExponent()].second.push_back(gf);
} }
} }
} }
...@@ -163,7 +145,6 @@ dgRungeKuttaMultirate::~dgRungeKuttaMultirate(){ ...@@ -163,7 +145,6 @@ dgRungeKuttaMultirate::~dgRungeKuttaMultirate(){
delete _currentInput; delete _currentInput;
delete _residualVolume; delete _residualVolume;
delete _residualInterface; delete _residualInterface;
delete _residualBoundary;
} }
void dgRungeKuttaMultirate::computeK(int iK,int exponent,bool isBuffer){ void dgRungeKuttaMultirate::computeK(int iK,int exponent,bool isBuffer){
...@@ -191,7 +172,7 @@ void dgRungeKuttaMultirate::computeK(int iK,int exponent,bool isBuffer){ ...@@ -191,7 +172,7 @@ void dgRungeKuttaMultirate::computeK(int iK,int exponent,bool isBuffer){
for(std::set<dgGroupOfFaces *>::iterator it=gOFSet.begin();it!=gOFSet.end();it++){ for(std::set<dgGroupOfFaces *>::iterator it=gOFSet.begin();it!=gOFSet.end();it++){
dgGroupOfFaces *faces=*it; dgGroupOfFaces *faces=*it;
if(faces->getNbGroupOfConnections()==1){ if(faces->getNbGroupOfConnections()==1){
_residualBoundary->computeAndMap1Group(*faces,*_currentInput,*_residual); _residualInterface->computeAndMap1Group(*faces,*_currentInput,*_residual);
} }
else{ else{
const dgGroupOfElements *gL = &faces->getGroupOfConnections(0).getGroupOfElements(); const dgGroupOfElements *gL = &faces->getGroupOfConnections(0).getGroupOfElements();
...@@ -243,7 +224,7 @@ void dgRungeKuttaMultirate::computeK(int iK,int exponent,bool isBuffer){ ...@@ -243,7 +224,7 @@ void dgRungeKuttaMultirate::computeK(int iK,int exponent,bool isBuffer){
for(int i=0;i<vf.size();i++){ for(int i=0;i<vf.size();i++){
dgGroupOfFaces *faces=vf[i]; dgGroupOfFaces *faces=vf[i];
if(faces->getNbGroupOfConnections()==1){ if(faces->getNbGroupOfConnections()==1){
_residualBoundary->computeAndMap1Group(*faces,*_currentInput,*_residual); _residualInterface->computeAndMap1Group(*faces,*_currentInput,*_residual);
} }
else{ else{
const dgGroupOfElements *gL = &faces->getGroupOfConnections(0).getGroupOfElements(); const dgGroupOfElements *gL = &faces->getGroupOfConnections(0).getGroupOfElements();
......
Solver/dgRungeKuttaMultirate.h
+ 0
− 1
View file @ e9522f01
...@@ -12,7 +12,6 @@ class dgRungeKuttaMultirate: public dgRungeKutta{ ...@@ -12,7 +12,6 @@ class dgRungeKuttaMultirate: public dgRungeKutta{
dgDofContainer *_residual; dgDofContainer *_residual;
dgResidualVolume *_residualVolume; dgResidualVolume *_residualVolume;
dgResidualInterface *_residualInterface; dgResidualInterface *_residualInterface;
dgResidualBoundary *_residualBoundary;
dgGroupCollection *_gc; dgGroupCollection *_gc;
std::vector<std::pair<std::vector<dgGroupOfElements*>,std::vector<dgGroupOfFaces*> > >_bulkGroupsOfElements;// int is the multirateExponent std::vector<std::pair<std::vector<dgGroupOfElements*>,std::vector<dgGroupOfFaces*> > >_bulkGroupsOfElements;// int is the multirateExponent
std::vector<std::pair<std::vector<dgGroupOfElements*>,std::vector<dgGroupOfFaces*> > >_innerBufferGroupsOfElements;// int is the multirateExponent std::vector<std::pair<std::vector<dgGroupOfElements*>,std::vector<dgGroupOfFaces*> > >_innerBufferGroupsOfElements;// int is the multirateExponent
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment