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

dg : jacobian implemented for volume term (not source) NOT TESTED

parent 01ca5b4c
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Common/LuaBindings.cpp
+ 6
− 0
View file @ 148885f7
......@@ -27,6 +27,7 @@
#include "dgFunctionIntegrator.h"
#include "Bindings.h"
#include "dgResidual.h"
#include "drawContext.h"
extern "C" {
#include "lua.h"
......@@ -196,6 +197,10 @@ static int luaClear (lua_State *L){
return 0;
}
#endif
static int luaRefresh (lua_State *L){
drawContext::global()->draw();
return 0;
}
int binding::readFile(const char *filename)
{
......@@ -319,6 +324,7 @@ binding::binding(){
#ifdef HAVE_READLINE
lua_register(L,"saveHistory",luaSave);
lua_register(L,"clearHistory",luaClear);
lua_register(L,"refreshGraphics",luaRefresh);
#endif
// lua_pushcfunction(L, luaopen_io);
......
Solver/dgGroupOfElements.cpp
+ 44
− 29
View file @ 148885f7
......@@ -64,15 +64,16 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e,
_order=polyOrder;
_dimUVW = _dimXYZ = e[0]->getDim();
// this is the biggest piece of data ... the mappings
int nbPsi = _fs.coefficients.size1();
_redistributionFluxes[0] = new fullMatrix<double> (nbPsi,_integration->size1());
_redistributionFluxes[1] = new fullMatrix<double> (nbPsi,_integration->size1());
_redistributionFluxes[2] = new fullMatrix<double> (nbPsi,_integration->size1());
_redistributionSource = new fullMatrix<double> (nbPsi,_integration->size1());
_collocation = new fullMatrix<double> (_integration->size1(),nbPsi);
_mapping = new fullMatrix<double> (e.size(), 10 * _integration->size1());
_imass = new fullMatrix<double> (nbPsi,nbPsi*e.size());
_dPsiDx = new fullMatrix<double> ( _integration->size1()*3,nbPsi*e.size());
int nbNodes = _fs.coefficients.size1();
int nbQP = _integration->size1();
_redistributionFluxes[0] = new fullMatrix<double> (nbNodes, nbQP);
_redistributionFluxes[1] = new fullMatrix<double> (nbNodes, nbQP);
_redistributionFluxes[2] = new fullMatrix<double> (nbNodes, nbQP);
_redistributionSource = new fullMatrix<double> (nbNodes, nbQP);
_collocation = new fullMatrix<double> (nbQP, nbNodes);
_mapping = new fullMatrix<double> (e.size(), 10 * nbQP);
_imass = new fullMatrix<double> (nbNodes,nbNodes*e.size());
_dPsiDx = new fullMatrix<double> ( nbQP*3,nbNodes*e.size());
_elementVolume = new fullMatrix<double> (e.size(),1);
_innerRadii = new fullMatrix<double> (e.size(),1);
double g[256][3],f[256];
......@@ -80,7 +81,7 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e,
for (int i=0;i<_elements.size();i++){
MElement *e = _elements[i];
element_to_index[e] = i;
fullMatrix<double> imass(*_imass,nbPsi*i,nbPsi);
fullMatrix<double> imass(*_imass,nbNodes*i,nbNodes);
(*_innerRadii)(i,0)=e->getInnerRadius();
for (int j=0;j< _integration->size1() ; j++ ){
_fs.f((*_integration)(j,0), (*_integration)(j,1), (*_integration)(j,2), f);
......@@ -100,14 +101,14 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e,
(*_mapping)(i,10*j + 6) = ijac[2][0];
(*_mapping)(i,10*j + 7) = ijac[2][1];
(*_mapping)(i,10*j + 8) = ijac[2][2];
for (int k=0;k<nbPsi;k++){
for (int l=0;l<nbPsi;l++) {
for (int k=0;k<nbNodes;k++){
for (int l=0;l<nbNodes;l++) {
imass(k,l) += f[k]*f[l]*weight*detjac;
}
// (iQP*3+iXYZ , iFace*NPsi+iPsi)
(*_dPsiDx)(j*3 ,i*nbPsi+k) = g[k][0]*ijac[0][0]+g[k][1]*ijac[0][1]+g[k][2]*ijac[0][2];
(*_dPsiDx)(j*3+1,i*nbPsi+k) = g[k][0]*ijac[1][0]+g[k][1]*ijac[1][1]+g[k][2]*ijac[1][2];
(*_dPsiDx)(j*3+2,i*nbPsi+k) = g[k][0]*ijac[2][0]+g[k][1]*ijac[2][1]+g[k][2]*ijac[2][2];
(*_dPsiDx)(j*3 ,i*nbNodes+k) = g[k][0]*ijac[0][0]+g[k][1]*ijac[0][1]+g[k][2]*ijac[0][2];
(*_dPsiDx)(j*3+1,i*nbNodes+k) = g[k][0]*ijac[1][0]+g[k][1]*ijac[1][1]+g[k][2]*ijac[1][2];
(*_dPsiDx)(j*3+2,i*nbNodes+k) = g[k][0]*ijac[2][0]+g[k][1]*ijac[2][1]+g[k][2]*ijac[2][2];
}
}
imass.invertInPlace();
......@@ -115,20 +116,34 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e,
// redistribution matrix
// quadrature weight x parametric gradients in quadrature points
for (int j=0;j<_integration->size1();j++) {
_fs.df((*_integration)(j,0),
(*_integration)(j,1),
(*_integration)(j,2), g);
_fs.f((*_integration)(j,0),
(*_integration)(j,1),
(*_integration)(j,2), f);
const double weight = (*_integration)(j,3);
for (int k=0;k<_fs.coefficients.size1();k++){
(*_redistributionFluxes[0])(k,j) = g[k][0] * weight;
(*_redistributionFluxes[1])(k,j) = g[k][1] * weight;
(*_redistributionFluxes[2])(k,j) = g[k][2] * weight;
(*_redistributionSource)(k,j) = f[k] * weight;
(*_collocation)(j,k) = f[k];
_PsiDPsiDXi = fullMatrix<double> (_dimUVW*nbQP, nbNodes*nbNodes);
//reditribution of the diffusive jacobian dimUVW*dimUVW*nbIntegrationPoints x nbNodes*nbNodes
_dPsiDXDPsiDXi = fullMatrix<double> (_dimUVW*_dimUVW*nbQP, nbNodes*nbNodes);
for (int xi=0;xi<nbQP; xi++) {
_fs.df((*_integration)(xi,0),
(*_integration)(xi,1),
(*_integration)(xi,2), g);
_fs.f((*_integration)(xi,0),
(*_integration)(xi,1),
(*_integration)(xi,2), f);
const double weight = (*_integration)(xi,3);
for (int k=0;k<nbNodes; k++){
(*_redistributionFluxes[0])(k,xi) = g[k][0] * weight;
(*_redistributionFluxes[1])(k,xi) = g[k][1] * weight;
(*_redistributionFluxes[2])(k,xi) = g[k][2] * weight;
(*_redistributionSource)(k,xi) = f[k] * weight;
(*_collocation)(xi,k) = f[k];
}
for (int alpha=0; alpha<_dimUVW; alpha++) for (int beta=0; beta<_dimUVW; beta++) {
for (int i=0; i<nbNodes; i++) for (int j=0; j<nbNodes; j++) {
_dPsiDXDPsiDXi((alpha*_dimUVW+beta)*nbQP+xi, i*nbNodes+j) = g[i][alpha]*g[j][beta]*weight;
}
}
for (int alpha=0; alpha<_dimUVW; alpha++){
for (int i=0; i<nbNodes; i++) for (int j=0; j<nbNodes; j++) {
_dPsiDXDPsiDXi(alpha*nbQP+xi, i*nbNodes+j) = g[i][alpha]*f[j]*weight;
}
}
}
}
......
Solver/dgGroupOfElements.h
+ 6
− 1
View file @ 148885f7
......@@ -73,6 +73,10 @@ class dgGroupOfElements {
fullMatrix<double> *_imass;
//
fullMatrix<double> *_dPsiDx;
//redistributions of the convective jacobian diumUVW*nbIntegrationPoints x nbNodes*nbNodes
fullMatrix<double> _PsiDPsiDXi;
//reditribution of the diffusive jacobian dimUVW*dimUVW*nbIntegrationPoints x nbNodes*nbNodes
fullMatrix<double> _dPsiDXDPsiDXi;
// dimension of the parametric space and of the real space
// may be different if the domain is a surface in 3D (manifold)
int _dimUVW, _dimXYZ;
......@@ -91,7 +95,6 @@ protected:
bool _multirateInnerBuffer;
bool _multirateOuterBuffer;
public:
inline int getMultirateExponent() const {return _multirateExponent;}
inline int getIsInnerMultirateBuffer() const {return _multirateInnerBuffer;}
inline int getIsOuterMultirateBuffer() const {return _multirateOuterBuffer;}
......@@ -109,6 +112,8 @@ public:
inline const fullMatrix<double> & getCollocationMatrix () const {return *_collocation;}
inline const fullMatrix<double> & getFluxRedistributionMatrix (int i) const {return *_redistributionFluxes[i];}
inline const fullMatrix<double> & getSourceRedistributionMatrix () const {return *_redistributionSource;}
inline const fullMatrix<double> & getDPsiDXDPsiDXi() const {return _dPsiDXDPsiDXi;}
inline const fullMatrix<double> & getPsiDPsiDXi() const {return _PsiDPsiDXi;}
inline double getElementVolume (int iElement)const {return (*_elementVolume)(iElement,0);}
inline double getInnerRadius(int iElement)const {return (*_innerRadii)(iElement,0);}
inline double getDetJ (int iElement, int iGaussPoint) const {return (*_mapping)(iElement, 10*iGaussPoint + 9);}
......
Solver/dgResidual.cpp
+ 152
− 0
View file @ 148885f7
......@@ -2,6 +2,7 @@
#include "dgConservationLaw.h"
#include "dgGroupOfElements.h"
#include "function.h"
#include "functionDerivator.h"
#include "dgDofContainer.h"
#include "MElement.h"
......@@ -108,6 +109,153 @@ void dgResidualVolume::compute1Group(dgGroupOfElements &group, fullMatrix<double
}
}
void dgResidualVolume::compute1GroupWithJacobian(dgGroupOfElements &group, fullMatrix<double> &solution, fullMatrix<double> &residual, fullMatrix<double> &jacobian)
{
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;
int nColA = group.getDimUVW()*group.getNbIntegrationPoints();
int nColB = group.getDimUVW()*group.getDimUVW()*group.getNbIntegrationPoints();
fullMatrix<double> A (_nbFields*_nbFields, nColA*group.getNbElements());
fullMatrix<double> B (_nbFields*_nbFields, nColB*group.getNbElements());
fullMatrix<double> a, b;
functionDerivator *dDiffusiveFluxdGradU,*dConvectiveFluxdU;
if(_diffusiveFlux)
dDiffusiveFluxdGradU=new functionDerivator(*_diffusiveFlux,_gradientSolutionQPe,1e-6);
if(_convectiveFlux)
dConvectiveFluxdU=new functionDerivator(*_convectiveFlux,_solutionQPe,1e-6);
// ----- 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 );
a.setAsProxy(A, iElement*nColA, nColA);
b.setAsProxy(B, iElement*nColB, nColB);
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;
}
}
}
int nQP = group.getNbIntegrationPoints();
int dim = group.getDimUVW();
if(_diffusiveFlux) {
dDiffusiveFluxdGradU->compute();
for (int alpha=0; alpha < dim; alpha++) for (int beta=0; beta < dim; beta++) {
for(int x=0;x <group.getDimXYZ();x++) for(int y=0;y <group.getDimXYZ();x++) {
for (int xi=0; xi <group.getNbIntegrationPoints(); xi++) {
const double invJx = group.getInvJ (iElement, xi, alpha, x);
const double invJy = group.getInvJ (iElement, xi, beta, y);
const double detJ = group.getDetJ (iElement, xi);
const double factor = invJx * invJy * detJ;
for (int k=0; k< _nbFields; k++) for (int l=0; l<_nbFields; l++) {
b(k*_nbFields+l,(alpha*dim+beta)*nQP+xi) = dDiffusiveFluxdGradU->get(k*dim+alpha,l*dim+beta,xi)*factor;
}
}
}
}
}
if(_convectiveFlux) {
dConvectiveFluxdU->compute();
for (int alpha=0; alpha < dim; alpha++) {
for(int x=0;x <group.getDimXYZ();x++) {
for (int xi=0; xi <group.getNbIntegrationPoints(); xi++){
const double invJ = group.getInvJ (iElement, xi, alpha, x);
const double detJ = group.getDetJ (iElement, xi);
const double factor = invJ * detJ;
for (int k=0; k< _nbFields; k++) for (int l=0; l<_nbFields; l++){
a(k*_nbFields+l,alpha*nQP+xi) = dConvectiveFluxdU->get(k*dim+alpha,l,xi)*factor;
}
}
}
}
}
}
// ----- 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);
}
int nbNodes = group.getNbNodes();
fullMatrix<double> jacobianK (_nbFields*_nbFields,nbNodes*nbNodes);
if (_convectiveFlux) {
jacobianK.gemm(group.getPsiDPsiDXi(),B);
}
if (_convectiveFlux) {
jacobianK.gemm(group.getDPsiDXDPsiDXi(),A);
}
fullMatrix<double> jacobianE, jacobianKE;
for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {
jacobianKE.setAsProxy(jacobianK, iElement*_nbFields*_nbFields, _nbFields*_nbFields);
jacobianKE.setAsProxy(jacobian, iElement*_nbFields*nbNodes, _nbFields*nbNodes);
for (int k=0; k<_nbFields;k++) for (int l=0;l<_nbFields;l++) {
for(int i=0; i<nbNodes; i++) for(int j=0; j<nbNodes; j++) {
jacobianE(l*nbNodes+j, k*nbNodes+i) = jacobianKE(k*_nbFields+l, i*nbNodes+j);
}
}
}
}
void dgResidualVolume::computeAndMap1Group(dgGroupOfElements &group, dgDofContainer &solution, dgDofContainer &residual)
{
compute1Group (group, solution.getGroupProxy(&group), residual.getGroupProxy(&group));
......@@ -384,6 +532,10 @@ void dgResidual::registerBindings (binding *b)
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);
mb = cb->addMethod("compute1GroupWithJacobian", &dgResidualVolume::compute1GroupWithJacobian);
mb->setDescription("compute the residual for one group given a fullMatrix proxy to the solution relative to this group");
mb->setArgNames("group", "solution", "residual", "jacobian", NULL);
cb = b->addClass<dgResidualInterface>("dgResidualInterface");
cb->setDescription("compute the residual for one dgGroupOfElements");
mb = cb->addMethod("computeAndMap1Group",&dgResidualInterface::computeAndMap1Group);
......
Solver/dgResidual.h
+ 1
− 0
View file @ 148885f7
......@@ -22,6 +22,7 @@ class dgResidualVolume {
dgResidualVolume (const dgConservationLaw &claw);
~dgResidualVolume();
void compute1Group(dgGroupOfElements &group, fullMatrix<double> &solution, fullMatrix<double> &residual);
void compute1GroupWithJacobian(dgGroupOfElements &group, fullMatrix<double> &solution, fullMatrix<double> &residual, fullMatrix<double> &jacobian);
void computeAndMap1Group(dgGroupOfElements &group, dgDofContainer &solution, dgDofContainer &residual);
};
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment