From 148885f72bf3b1ce8fde3f56cba6b43d50108719 Mon Sep 17 00:00:00 2001
From: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date: Sat, 6 Mar 2010 08:26:41 +0000
Subject: [PATCH] dg : jacobian implemented for volume term (not source) NOT
TESTED
---
Common/LuaBindings.cpp | 6 ++
Solver/dgGroupOfElements.cpp | 73 ++++++++++-------
Solver/dgGroupOfElements.h | 7 +-
Solver/dgResidual.cpp | 152 +++++++++++++++++++++++++++++++++++
Solver/dgResidual.h | 1 +
5 files changed, 209 insertions(+), 30 deletions(-)
diff --git a/Common/LuaBindings.cpp b/Common/LuaBindings.cpp
index 0f2fa8fa87..8eba43f94c 100644
--- a/Common/LuaBindings.cpp
+++ b/Common/LuaBindings.cpp
@@ -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);
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index 1836721e21..08c0e9c256 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -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;
+ }
}
}
}
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 987e86593a..dc22cfcfd2 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -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);}
diff --git a/Solver/dgResidual.cpp b/Solver/dgResidual.cpp
index c8f34f5647..0289ba87d3 100644
--- a/Solver/dgResidual.cpp
+++ b/Solver/dgResidual.cpp
@@ -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);
diff --git a/Solver/dgResidual.h b/Solver/dgResidual.h
index ed7a1c678d..47ab9226d2 100644
--- a/Solver/dgResidual.h
+++ b/Solver/dgResidual.h
@@ -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);
};
--
GitLab