dg : implicit matrix ok for diffusion volume term

Solver/TESTCASES/implicit.lua

+model = GModel()
+model:load ('square2.msh')
+order = 1
+dimension = 2
+
+law = dgConservationLawAdvectionDiffusion:diffusionLaw(functionConstant({1}))
+law:addBoundaryCondition('Border',law:new0FluxBoundary())
+groups = dgGroupCollection(model, dimension, order)
+groups:buildGroupsOfInterfaces()
+solution = dgDofContainer(groups, law:getNbFields())
+residual = dgDofContainer(groups, law:getNbFields())
+r = dgResidualVolume (law)
+g = groups:getElementGroup(0)
+j = fullMatrix(3,3*g:getNbElements())
+r:compute1GroupWithJacobian(g,solution:getGroupProxy(0),residual:getGroupProxy(0),j)
+j:print("test")

Solver/dgGroupOfElements.cpp

@@ -121,7 +121,7 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e,
 
   _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);
+  _dPsiDXDPsiDXi = fullMatrix<double> (nbNodes*nbNodes, _dimUVW*_dimUVW*nbQP);
 
   for (int xi=0;xi<nbQP; xi++) {
     _fs.df((*_integration)(xi,0),
@@ -145,7 +145,7 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e,
     }
     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;
+        _dPsiDXDPsiDXi(i*nbNodes+j, (alpha*_dimUVW+beta)*nbQP+xi) = g[i][alpha]*g[j][beta]*weight;
       }
     }
     for (int alpha=0; alpha<_dimUVW; alpha++){
@@ -1212,6 +1212,8 @@ void dgGroupCollection::registerBindings(binding *b)
   methodBinding *cm;
   cb = b->addClass<dgGroupOfElements>("dgGroupOfElements");
   cb->setDescription("a group of element sharing the same properties");
+  cm = cb->addMethod("getNbElements",&dgGroupOfElements::getNbElements);
+  cm->setDescription("return the number of elements in the group");
   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");

Solver/dgResidual.cpp

@@ -135,10 +135,8 @@ void dgResidualVolume::compute1GroupWithJacobian(dgGroupOfElements &group, fullM
   fullMatrix<double> fuvwe;
   fullMatrix<double> source;
   fullMatrix<double> dPsiDx, dofs, gradSolProxy;
-  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 (group.getDimUVW()*group.getNbIntegrationPoints(), _nbFields*_nbFields*group.getNbElements());
+  fullMatrix<double> B (group.getDimUVW()*group.getDimUVW()*group.getNbIntegrationPoints(), _nbFields*_nbFields*group.getNbElements());
   fullMatrix<double> a, b;
   functionDerivator *dDiffusiveFluxdGradU,*dConvectiveFluxdU;
   if(_diffusiveFlux)
@@ -149,8 +147,8 @@ void dgResidualVolume::compute1GroupWithJacobian(dgGroupOfElements &group, fullM
   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);
+    a.setAsProxy(A, iElement*_nbFields*_nbFields, _nbFields*_nbFields);
+    b.setAsProxy(B, iElement*_nbFields*_nbFields, _nbFields*_nbFields);
 
     if(_gradientSolutionQPe.somethingDependOnMe()){
       dPsiDx.setAsProxy(group.getDPsiDx(),iElement*group.getNbNodes(),group.getNbNodes());
@@ -194,27 +192,26 @@ void dgResidualVolume::compute1GroupWithJacobian(dgGroupOfElements &group, fullM
       }
     }
     int nQP = group.getNbIntegrationPoints();
-    int dim = group.getDimUVW();
+    int dimUVW = group.getDimUVW();
+    int dimXYZ = group.getDimXYZ();
+    b.scale(0);
     if(_diffusiveFlux) {
-      printf("ok %i\n", __LINE__);
       dDiffusiveFluxdGradU->compute();
-      printf("ok %i\n", __LINE__);
-      exit(0);
-      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 alpha=0; alpha < dimUVW; alpha++) for (int beta=0; beta < dimUVW; beta++) {
+        for(int x=0;x <dimXYZ;x++) for(int y=0;y <dimXYZ;y++) {
           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;
+              b((alpha*dimUVW+beta)*nQP+xi,k*_nbFields+l) += dDiffusiveFluxdGradU->get(xi, k*dimXYZ+x, l*dimXYZ+y)*factor;
             }
           }
         }
       }
     }
-    if(_convectiveFlux) {
+    /*if(_convectiveFlux) {
       dConvectiveFluxdU->compute();
       for (int alpha=0; alpha < dim; alpha++) {
         for(int x=0;x <group.getDimXYZ();x++) {
@@ -228,33 +225,33 @@ void dgResidualVolume::compute1GroupWithJacobian(dgGroupOfElements &group, fullM
           }
         }
       }
-    }
+    }*/
   }
 
   // ----- 3 ---- do the redistribution at nodes using as many BLAS3 operations as there are local coordinates
-  if(_convectiveFlux || _diffusiveFlux){
+  /*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> jacobianK (nbNodes*nbNodes,_nbFields*_nbFields*group.getNbElements());
+  /*if (_convectiveFlux) {
+    jacobianK.gemm(group.getPsiDPsiDXi(),A);
+  }*/
+  if (_diffusiveFlux) {
+    jacobianK.gemm(group.getDPsiDXDPsiDXi(),B);
   }
   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);
+    jacobianE.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);
+        jacobianE(l*nbNodes+j, k*nbNodes+i) = jacobianKE(i*nbNodes+j, k*_nbFields+l);
       }
     }
   }

Solver/functionDerivator.cpp

@@ -4,7 +4,6 @@
 void functionDerivator::compute() {
   _xRef = _x();
   _fRef = _f();
-  printf("f.size = %i %i\n_x.size = %i %i\n",_f().size1(), _f().size2(), _x().size1(), _x().size2());
   _dfdx = fullMatrix<double>(_f().size1(),_f().size2()*_x().size2());
   for (int j=0;j<_xRef.size2();j++) {
     _xDx = _xRef;
@@ -16,5 +15,4 @@ void functionDerivator::compute() {
         _dfdx(i, k*_xRef.size2()+j) = (_f(i,k)-_fRef(i,k))/_epsilon;
   }
   _x.set(_xRef);
-  _dfdx.print();
 };