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

dg : bind the residual at the level of one group

parent b265aa84
No related branches found
No related tags found
No related merge requests found
Expand all Show whitespace changes
Inline Side-by-side
Showing
with 550 additions and 409 deletions
Common/LuaBindings.cpp
+ 2
− 0
View file @ c58fd0ba
......@@ -23,6 +23,7 @@
#include "dgSystemOfEquations.h"
#include "dgLimiter.h"
#include "Bindings.h"
#include "dgResidual.h"
extern "C" {
#include "lua.h"
......@@ -340,6 +341,7 @@ binding::binding(){
dgGroupCollection::registerBindings(this);
dgLimiter::registerBindings(this);
dgPerfectGasLaw2dRegisterBindings(this);
dgResidual::registerBindings(this);
dgRungeKutta::registerBindings(this);
dgSlopeLimiterRegisterBindings(this);
dgSystemOfEquations::registerBindings(this);
......
Common/LuaBindings.h
+ 47
− 0
View file @ c58fd0ba
......@@ -252,6 +252,33 @@ class luaStack<const type &>{
template <typename cb>
class argTypeNames;
template <typename tr, typename tObj, typename t0, typename t1, typename t2, typename t3, typename t4, typename t5>
class argTypeNames<tr (tObj::*)(t0,t1,t2,t3,t4,t5)>{
public:
static void get(std::vector<std::string> &names){
names.clear();
names.push_back(luaStack<tr>::getName());
names.push_back(luaStack<t0>::getName());
names.push_back(luaStack<t1>::getName());
names.push_back(luaStack<t2>::getName());
names.push_back(luaStack<t3>::getName());
names.push_back(luaStack<t4>::getName());
names.push_back(luaStack<t5>::getName());
}
};
template <typename tr, typename tObj, typename t0, typename t1, typename t2, typename t3, typename t4>
class argTypeNames<tr (tObj::*)(t0,t1,t2,t3,t4)>{
public:
static void get(std::vector<std::string> &names){
names.clear();
names.push_back(luaStack<tr>::getName());
names.push_back(luaStack<t0>::getName());
names.push_back(luaStack<t1>::getName());
names.push_back(luaStack<t2>::getName());
names.push_back(luaStack<t3>::getName());
names.push_back(luaStack<t4>::getName());
}
};
template <typename tr, typename tObj, typename t0, typename t1, typename t2, typename t3>
class argTypeNames<tr (tObj::*)(t0,t1,t2,t3)>{
public:
......@@ -401,6 +428,16 @@ static int luaCall(lua_State *L,tRet (tObj::*_f)() const) {
};
//non const, return
template <typename tObj, typename tRet, typename t0, typename t1, typename t2, typename t3, typename t4, typename t5>
static int luaCall(lua_State *L,tRet (tObj::*_f)(t0,t1,t2,t3,t4,t5)) {
luaStack<tRet>::push(L,(luaStack<tObj*>::get(L,1)->*(_f))(luaStack<t0>::get(L,2),luaStack<t1>::get(L,3),luaStack<t2>::get(L,4),luaStack<t3>::get(L,5),luaStack<t4>::get(L,6),luaStack<t5>::get(L,7)));
return 1;
};
template <typename tObj, typename tRet, typename t0, typename t1, typename t2, typename t3, typename t4>
static int luaCall(lua_State *L,tRet (tObj::*_f)(t0,t1,t2,t3,t4)) {
luaStack<tRet>::push(L,(luaStack<tObj*>::get(L,1)->*(_f))(luaStack<t0>::get(L,2),luaStack<t1>::get(L,3),luaStack<t2>::get(L,4),luaStack<t3>::get(L,5),luaStack<t4>::get(L,6)));
return 1;
};
template <typename tObj, typename tRet, typename t0, typename t1, typename t2, typename t3>
static int luaCall(lua_State *L,tRet (tObj::*_f)(t0,t1,t2,t3)) {
luaStack<tRet>::push(L,(luaStack<tObj*>::get(L,1)->*(_f))(luaStack<t0>::get(L,2),luaStack<t1>::get(L,3),luaStack<t2>::get(L,4),luaStack<t3>::get(L,5)));
......@@ -457,6 +494,16 @@ static int luaCall(lua_State *L,void (tObj::*_f)() const) {
//non const, no return
template <typename tObj, typename t0, typename t1, typename t2, typename t3, typename t4, typename t5>
static int luaCall(lua_State *L,void (tObj::*_f)(t0,t1,t2,t3,t4,t5)) {
(luaStack<tObj*>::get(L,1)->*(_f))(luaStack<t0>::get(L,2),luaStack<t1>::get(L,3),luaStack<t2>::get(L,4),luaStack<t3>::get(L,5),luaStack<t4>::get(L,6),luaStack<t5>::get(L,7));
return 1;
};
template <typename tObj, typename t0, typename t1, typename t2, typename t3, typename t4>
static int luaCall(lua_State *L,void (tObj::*_f)(t0,t1,t2,t3,t4)) {
(luaStack<tObj*>::get(L,1)->*(_f))(luaStack<t0>::get(L,2),luaStack<t1>::get(L,3),luaStack<t2>::get(L,4),luaStack<t3>::get(L,5),luaStack<t4>::get(L,6));
return 1;
};
template <typename tObj, typename t0, typename t1, typename t2, typename t3>
static int luaCall(lua_State *L,void (tObj::*_f)(t0,t1,t2,t3)) {
(luaStack<tObj*>::get(L,1)->*(_f))(luaStack<t0>::get(L,2),luaStack<t1>::get(L,3),luaStack<t2>::get(L,4),luaStack<t3>::get(L,5));
......
Solver/CMakeLists.txt
+ 1
− 0
View file @ c58fd0ba
......@@ -22,6 +22,7 @@ set(SRC
dgConservationLawMaxwell.cpp
dgConservationLawPerfectGas.cpp
dgLimiter.cpp
dgResidual.cpp
dgDofContainer.cpp
function.cpp
functionDerivator.cpp
......
Solver/dgAlgorithm.cpp
+ 0
− 379
View file @ c58fd0ba
......@@ -16,387 +16,8 @@
\int \vec{f} \cdot \grad \phi dv
*/
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)
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);
// ----- 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)};
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()){
dPsiDx.setAsProxy(group.getDPsiDx(),iElement*group.getNbNodes(),group.getNbNodes());
dofs.setAsProxy(solution, _nbFields*iElement, _nbFields);
dPsiDx.mult(dofs, _gradientSolutionQPe.set());
}
_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);
// ----- 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++) {
// get the right inv jacobian matrix dU/dX element
const double invJ = group.getInvJ (iElement, iPt, iUVW, iXYZ);
// get the detJ at this point
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;
}
}
}
}
}
}
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;
}
}
}
}
// ----- 3 ---- do the redistribution at nodes using as many BLAS3 operations as there are local coordinates
if(_convectiveFlux || _diffusiveFlux){
for (int iUVW=0;iUVW<group.getDimUVW();iUVW++){
residual.gemm(group.getFluxRedistributionMatrix(iUVW),Fuvw[iUVW]);
}
}
if(_sourceTerm){
residual.gemm(group.getSourceRedistributionMatrix(),Source);
}
}
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,
fullMatrix<double> &solutionRight,
fullMatrix<double> &residual // residual !! at faces nodes
)
{
// A) global operations before entering the loop over the faces
// A1 ) copy locally some references from the group for the sake of lisibility
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);
group.getCollocationMatrix().mult(solution, solutionQP);
// needed tocompute normal fluxes at integration points
fullMatrix<double> NormalFluxQP ( group.getNbIntegrationPoints(), nbFaces*_nbFields*2);
// create one dataCache for each side
_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));
// 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());
// 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));
// proxies for the flux
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());
}
// B3 ) compute fluxes
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++) {
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;
double minl = std::min(group.getElementVolumeLeft(iFace),
group.getElementVolumeRight(iFace)
)/group.getInterfaceSurface(iFace);
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;
normalFluxQP(iPt,k) -= (meanNormalFlux+solutionJumpPenalty)*detJ;
normalFluxQP(iPt,k+_nbFields) += (meanNormalFlux+solutionJumpPenalty)*detJ;
}
}
}
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;
}
}
}
// C ) redistribute the flux to the residual (at Faces nodes)
if(_riemannSolver || _diffusiveFluxLeft)
residual.gemm(group.getRedistributionMatrix(),NormalFluxQP);
}
void dgResidualInterface::computeAndMap1Group (dgGroupOfFaces &faces, dgDofContainer &solution, dgDofContainer &residual)
{
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 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());
// ----- 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());
fullMatrix<double> &DPsiLeftDx = group.getDPsiLeftDxMatrix();
// provided dataCache
dataCacheDouble &uvw=cacheMapLeft.provideData("UVW",1,3);
dataCacheDouble &solutionQPLeft = cacheMapLeft.provideData("Solution",1,nbFields);
dataCacheDouble &gradientSolutionLeft = cacheMapLeft.provideData("SolutionGradient",3,nbFields);
dataCacheDouble &normals = cacheMapLeft.provideData("Normals",1,1);
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 = group.getGroupLeft().getOrder();
int dim = group.getGroupLeft().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(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);
dofsLeft.setAsProxy(solutionLeft, nbFields*group.getElementLeftId(iFace), nbFields);
uvw.setAsProxy(group.getIntegrationOnElementLeft(iFace));
cacheElementLeft.set(group.getElementLeft(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 = group.getElementVolumeLeft(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();
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 &residual)
{
solution.scatter();
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)
const dgConservationLaw &claw, // the conservation law
......
Solver/dgAlgorithm.h
+ 0
− 11
View file @ c58fd0ba
......@@ -15,23 +15,12 @@ class dgSystemOfEquations;
class dgAlgorithm {
public :
static void registerBindings(binding *b);
static void 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
);
static void residual( const dgConservationLaw &claw, dgGroupCollection &groups,
dgDofContainer &solution, dgDofContainer &residual);
static void computeElementaryTimeSteps ( //dofManager &dof, // the DOF manager (maybe useless here)
const dgConservationLaw &claw, // the conservation law
const dgGroupOfElements & group, // the group
fullMatrix<double> &solution, // the solution
std::vector<double> & DT // elementary time steps
);
static void multAddInverseMassMatrix ( /*dofManager &dof,*/
const dgGroupOfElements & group,
fullMatrix<double> &residu,
......
Solver/dgConservationLawMaxwell.cpp
+ 4
− 12
View file @ c58fd0ba
......@@ -7,14 +7,12 @@ class dgConservationLawMaxwell::advection : public dataCacheDouble {
dataCacheDouble &sol, &mu_eps;
public:
advection(std::string mu_epsFunctionName,dataCacheMap &cacheMap):
dataCacheDouble(cacheMap),
dataCacheDouble(cacheMap,1,18),
sol(cacheMap.get("Solution",this)),
mu_eps(cacheMap.get(mu_epsFunctionName,this))
{};
void _eval () {
int nQP = sol().size1();
if(_value.size1() != nQP)
_value=fullMatrix<double>(nQP,18);
for(int i=0; i< nQP; i++) {
double Ex = sol(i,0);
double Ey = sol(i,1);
......@@ -54,14 +52,12 @@ class dgConservationLawMaxwell::source: public dataCacheDouble {
dataCacheDouble &xyz, &solution;
public :
source(dataCacheMap &cacheMap) :
dataCacheDouble(cacheMap),
dataCacheDouble(cacheMap,1,6),
xyz(cacheMap.get("XYZ",this)),
solution(cacheMap.get("Solution",this))
{}
void _eval () {
int nQP = xyz().size1();
if(_value.size1() != nQP)
_value = fullMatrix<double>(nQP,6);
for (int i=0; i<nQP; i++) {
_value (i,0) = 0;
......@@ -79,7 +75,7 @@ class dgConservationLawMaxwell::riemann:public dataCacheDouble {
dataCacheDouble &normals, &solL, &solR, &mu_eps;
public:
riemann(std::string mu_epsFunctionName,dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight):
dataCacheDouble(cacheMapLeft),
dataCacheDouble(cacheMapLeft,1,12),
normals(cacheMapLeft.get("Normals", this)),
solL(cacheMapLeft.get("Solution", this)),
solR(cacheMapRight.get("Solution", this)),
......@@ -87,8 +83,6 @@ class dgConservationLawMaxwell::riemann:public dataCacheDouble {
{};
void _eval () {
int nQP = solL().size1();
if(_value.size1() != nQP)
_value = fullMatrix<double>(nQP,12);
for(int i=0; i< nQP; i++) {
double nx= normals(0,i);
......@@ -191,14 +185,12 @@ class dgBoundaryConditionMaxwellWall : public dgBoundaryCondition {
dataCacheDouble &sol,&normals;
public:
term(dataCacheMap &cacheMap):
dataCacheDouble(cacheMap),
dataCacheDouble(cacheMap,1,6),
sol(cacheMap.get("Solution",this)),
normals(cacheMap.get("Normals",this))
{}
void _eval () {
int nQP = sol().size1();
if(_value.size1() != nQP)
_value = fullMatrix<double>(nQP,6);
for(int i=0; i< nQP; i++) {
double nx= normals(0,i);
......
Solver/dgGroupOfElements.cpp
+ 23
− 2
View file @ c58fd0ba
......@@ -929,9 +929,15 @@ void dgGroupCollection::find (MElement*e, int &ig, int &index){
#include "LuaBindings.h"
void dgGroupCollection::registerBindings(binding *b){
classBinding *cb = b->addClass<dgGroupCollection>("dgGroupCollection");
cb->setDescription("The GroupCollection class let you access to group partitioning functions");
classBinding *cb;
methodBinding *cm;
cb = b->addClass<dgGroupOfElements>("dgGroupOfElements");
cb->setDescription("a group of element sharing the same properties");
cb = b->addClass<dgGroupOfFaces>("dgGroupOfFaces");
cb->setDescription("a group of faces. All faces of this group are either interior of the same group of elements or at the boundary between two group of elements");
cb = b->addClass<dgGroupCollection>("dgGroupCollection");
cb->setDescription("The GroupCollection class let you access to group partitioning functions");
cm = cb->setConstructor<dgGroupCollection,GModel*,int,int>();
cm->setDescription("Build the group of elements, separated by element type and element order. Additional partitioning can be done using splitGroupsForXXX specific functions");
cm->setArgNames("model","dimension","order",NULL);
......@@ -940,4 +946,19 @@ void dgGroupCollection::registerBindings(binding *b){
cm = cb->addMethod("splitGroupsForMultirate",&dgGroupCollection::splitGroupsForMultirate);
cm->setDescription("Split the groups according to their own stable time step");
cm->setArgNames("refDt","claw",NULL);
cm = cb->addMethod("getNbElementGroups", &dgGroupCollection::getNbElementGroups);
cm->setDescription("return the number of dgGroupOfElements");
cm = cb->addMethod("getNbFaceGroups", &dgGroupCollection::getNbFaceGroups);
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->setDescription("get 1 group of elements");
cm->setArgNames("id", NULL);
cm = cb->addMethod("getFaceGroup", &dgGroupCollection::getFaceGroup);
cm->setDescription("get 1 group of faces");
cm->setArgNames("id", NULL);
cm = cb->addMethod("getBoundaryGroup", &dgGroupCollection::getBoundaryGroup);
cm->setDescription("get 1 group of boundary faces");
cm->setArgNames("id", NULL);
}
Solver/dgResidual.cpp 0 → 100644
+ 398
− 0
View file @ c58fd0ba
This diff is collapsed.
Solver/dgResidual.h 0 → 100644
+ 70
− 0
View file @ c58fd0ba
#ifndef _DG_RESIDUAL_H_
#define _DG_RESIDUAL_H_
class dgConservationLaw;
class dgGroupCollection;
class dataCacheDouble;
class dataCacheElement;
class dgDofContainer;
class dataCacheMap;
class dgGroupOfElements;
class dgGroupOfFaces;
class binding;
#include "fullMatrix.h"
class dgResidualVolume {
dataCacheMap *_cacheMap;
const dgConservationLaw &_claw;
int _nbFields;
dataCacheElement &_cacheElement;
dataCacheDouble &_UVW, &_solutionQPe, &_gradientSolutionQPe;
dataCacheDouble *_sourceTerm, *_convectiveFlux, *_diffusiveFlux;
public:
dgResidualVolume (const dgConservationLaw &claw);
~dgResidualVolume();
void compute1Group(dgGroupOfElements &group, fullMatrix<double> &solution, fullMatrix<double> &residual);
void computeAndMap1Group(dgGroupOfElements &group, dgDofContainer &solution, dgDofContainer &residual);
};
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();
};
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);
};
class dgResidual {
const dgConservationLaw &_claw;
public:
dgResidual (const dgConservationLaw *claw) : _claw(*claw) {}
void compute(dgGroupCollection &groups, dgDofContainer &solution, dgDofContainer &residual);
static void registerBindings (binding *b);
};
#endif
Solver/dgRungeKutta.cpp
+ 5
− 5
View file @ c58fd0ba
......@@ -2,7 +2,7 @@
#include "dgConservationLaw.h"
#include "dgDofContainer.h"
#include "dgLimiter.h"
#include "dgAlgorithm.h"
#include "dgResidual.h"
#include "dgGroupOfElements.h"
double dgRungeKutta::iterateEuler(const dgConservationLaw *claw, double dt, dgDofContainer *solution) {
......@@ -111,6 +111,7 @@ double dgRungeKutta::diagonalRK (const dgConservationLaw *claw,
K.axpy(*sol);
Unp.scale(0.);
Unp.axpy(*sol);
dgResidual residual(claw);
for(int j=0; j<nStages;j++){
if(j){
......@@ -119,7 +120,7 @@ double dgRungeKutta::diagonalRK (const dgConservationLaw *claw,
}
if (_limiter) _limiter->apply(&K);
dgAlgorithm::residual(*claw,*groups,K,resd);
residual.compute(*groups,K,resd);
K.scale(0.);
for(int k=0; k < groups->getNbElementGroups(); k++) {
dgGroupOfElements *group = groups->getElementGroup(k);
......@@ -139,8 +140,6 @@ double dgRungeKutta::diagonalRK (const dgConservationLaw *claw,
return sol->norm();
}
double dgRungeKutta::nonDiagonalRK(const dgConservationLaw *claw,
double dt,
dgDofContainer *solution,
......@@ -168,6 +167,7 @@ double dgRungeKutta::nonDiagonalRK(const dgConservationLaw *claw,
Unp.scale(0.0);
Unp.axpy(*solution);
dgResidual residual(claw);
for(int i=0; i<nStages;i++){
tmp.scale(0.0);
tmp.axpy(*solution);
......@@ -177,7 +177,7 @@ double dgRungeKutta::nonDiagonalRK(const dgConservationLaw *claw,
}
}
if (_limiter) _limiter->apply(&tmp);
dgAlgorithm::residual(*claw,*groups,tmp,resd);
residual.compute(*groups,tmp,resd);
// Multiply the spatial residual by the inverse of the mass matrix
for(int k=0; k < groups->getNbElementGroups(); k++) {
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment