diff --git a/Common/LuaBindings.cpp b/Common/LuaBindings.cpp
index 50e271064fbc6cc3b326833ce6dd046b24619502..0011274ee36142c3784c687412ef70059be0c106 100644
--- a/Common/LuaBindings.cpp
+++ b/Common/LuaBindings.cpp
@@ -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);
diff --git a/Common/LuaBindings.h b/Common/LuaBindings.h
index 22c476fd609ca531d8638f438bdae374e926c936..b625cef9aa8a64c6474bc33b5f7158480abed14a 100644
--- a/Common/LuaBindings.h
+++ b/Common/LuaBindings.h
@@ -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));
diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt
index 12e9ede52556e4a9c5104915bcd1f1f9bba77e30..e7a29aefc23ae52711181d5effb1e819aadf2b04 100644
--- a/Solver/CMakeLists.txt
+++ b/Solver/CMakeLists.txt
@@ -22,6 +22,7 @@ set(SRC
dgConservationLawMaxwell.cpp
dgConservationLawPerfectGas.cpp
dgLimiter.cpp
+ dgResidual.cpp
dgDofContainer.cpp
function.cpp
functionDerivator.cpp
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index 9e14abe960f7972a376269446b0fa04d5816171f..4246e56c1ca33ed786a51a7cf05d205795b088e0 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -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
diff --git a/Solver/dgAlgorithm.h b/Solver/dgAlgorithm.h
index 6032f9a639f8c7323fa549fcd06ed3e5ff21ca6e..5756cce915bcae4ddfd70b99e682788307d1b9d7 100644
--- a/Solver/dgAlgorithm.h
+++ b/Solver/dgAlgorithm.h
@@ -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,
diff --git a/Solver/dgConservationLawMaxwell.cpp b/Solver/dgConservationLawMaxwell.cpp
index a4f1b1b1d955b37314ba86bea5b0c4ee9dd697c2..9ca9c3846e83e65f5e592790bcb3ab25d7069bc9 100644
--- a/Solver/dgConservationLawMaxwell.cpp
+++ b/Solver/dgConservationLawMaxwell.cpp
@@ -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);
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index 824d280832b6449eb209ca0e86ce9cf534bcb94e..c21ccaf325f7af74111c244e99859be4fb4a60ad 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -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);
}
diff --git a/Solver/dgResidual.cpp b/Solver/dgResidual.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ade2dbbf2883f302cacb7bc67ec7c99e9d36f13a
--- /dev/null
+++ b/Solver/dgResidual.cpp
@@ -0,0 +1,398 @@
+#include "dgResidual.h"
+#include "dgConservationLaw.h"
+#include "dgGroupOfElements.h"
+#include "function.h"
+#include "dgDofContainer.h"
+#include "MElement.h"
+
+dgResidualVolume::dgResidualVolume(const dgConservationLaw &claw):
+ _cacheMap(new dataCacheMap),
+ _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))
+{
+}
+
+dgResidualVolume::~dgResidualVolume()
+{
+ delete _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));
+}
+
+
+dgResidualInterface::dgResidualInterface (const dgConservationLaw &claw) :
+ _cacheMapLeft (new dataCacheMap),
+ _cacheMapRight (new dataCacheMap),
+ _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()));
+}
+
+dgResidualInterface::~dgResidualInterface ()
+{
+ delete _cacheMapLeft;
+ delete _cacheMapRight;
+}
+
+
+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
+ )
+{
+ //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
+
+ 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()));
+}
+
+void dgResidual::compute(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);
+}
+
+#include "Bindings.h"
+void dgResidual::registerBindings (binding *b)
+{
+ classBinding *cb = b->addClass<dgResidual>("dgResidual");
+ cb->setDescription("compute the residual, given a conservation law and a solution dof container");
+ methodBinding *mb;
+ mb = cb->addMethod("compute", &dgResidual::compute);
+ mb->setDescription("compute the residual");
+ mb->setArgNames("groups", "solution", "residual",NULL);
+ mb = cb->setConstructor<dgResidual, dgConservationLaw*>();
+ mb->setDescription("a new object used to compute the residual");
+ mb->setArgNames("law",NULL);
+
+ cb = b->addClass<dgResidualVolume>("dgResidualVolume");
+ cb->setDescription("compute the residual for one dgGroupOfFaces");
+ mb = cb->addMethod("computeAndMap1Group",&dgResidualVolume::computeAndMap1Group);
+ mb->setDescription("compute the residual for one group given a dgDofContainer solution");
+ mb->setArgNames("group", "solution", "residual", NULL);
+ mb = cb->addMethod("compute1Group", &dgResidualVolume::compute1Group);
+ mb->setDescription("compute the residual for one group given a fullMatrix proxy to the solution relative to this group");
+ mb->setArgNames("group", "solution", "residual", NULL);
+
+ cb = b->addClass<dgResidualInterface>("dgResidualInterface");
+ cb->setDescription("compute the residual for one dgGroupOfElements");
+ mb = cb->addMethod("computeAndMap1Group",&dgResidualInterface::computeAndMap1Group);
+ mb->setDescription("compute the residual for one group given a dgDofContainer solution");
+ mb->setArgNames("group", "solution", "residual", NULL);
+ 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->setArgNames("group", "solutionFaces", "solutionGroupLeft", "solutionGroupRight", "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);
+}
diff --git a/Solver/dgResidual.h b/Solver/dgResidual.h
new file mode 100644
index 0000000000000000000000000000000000000000..ed7a1c678d60aeb55120aad7dda1c655dc5eb84e
--- /dev/null
+++ b/Solver/dgResidual.h
@@ -0,0 +1,70 @@
+#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
diff --git a/Solver/dgRungeKutta.cpp b/Solver/dgRungeKutta.cpp
index 3480c5f5cc31307e90475d9a70a925255b772476..2e6898738a1b3b466359ad04369c93af5b45a56d 100644
--- a/Solver/dgRungeKutta.cpp
+++ b/Solver/dgRungeKutta.cpp
@@ -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++) {