Runge-Kutta

Solver/dgAlgorithm.cpp

@@ -15,8 +15,8 @@
 void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe useless here)
 				   const dgConservationLaw &claw,   // the conservation law
 				   const dgGroupOfElements & group, 
-           const fullMatrix<double> &solution,
-           fullMatrix<double> &residual // the residual
+				   const fullMatrix<double> &solution,
+				   fullMatrix<double> &residual // the residual
            )
 { 
   // ----- 1 ----  get the solution at quadrature points
@@ -163,6 +163,7 @@ void dgAlgorithm::multAddInverseMassMatrix ( /*dofManager &dof,*/
           const dgGroupOfElements & group,
           fullMatrix<double> &residu,
           fullMatrix<double> &sol)
+		  
 {
   for(int i=0;i<group.getNbElements();i++) {
     fullMatrix<double> residuEl(residu,i,1);
@@ -171,6 +172,47 @@ void dgAlgorithm::multAddInverseMassMatrix ( /*dofManager &dof,*/
   }
 }
 
+ void dgAlgorithm::rungeKutta (const dgConservationLaw &claw, // conservation law
+		  std::vector<dgGroupOfElements*> &eGroups, //group of elements
+		  std::vector<dgGroupOfFaces*> &fGroups,  // group of interfacs
+		  std::vector<dgGroupOfFaces*> &bGroups, // group of boundaries
+		  double h,		
+          fullMatrix<double> &residu,
+          fullMatrix<double> &sol,
+		  int orderRK)
+{
+
+
+// U_{n+1}=U_n+h*(SUM_i a_i*K_i)
+// K_i=M^(-1)R(U_n+b_i*K_{i-1})
+  
+  double a[4] = {h/6.0,h/3.0,h/3.0,h/6.0};
+  double b[4] = {0.,h/2.0,h/2.0,h};
+  	
+  fullMatrix<double> K(sol);
+  // Current updated solution
+  fullMatrix<double> Unp(sol);
+  
+  for(int j=0; j<orderRK;j++){
+	if(j!=0){
+		K.scale(b[j]);
+		K.add(sol);
+		}
+	this->residual(claw,eGroups,fGroups,bGroups,K,residu);
+	K.scale(0.);
+	for(int i=0;i<eGroups[0]->getNbElements();i++) {
+		fullMatrix<double> residuEl(residu,i,1);
+		fullMatrix<double> KEl(K,i,1);
+		(eGroups[0]->getInverseMassMatrix(i)).mult(residuEl,KEl);
+    }
+	
+	Unp.add(K,a[j]);
+  }
+  
+  sol=Unp;
+}
+
+
 void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (maybe useless here)
 				   const dgConservationLaw &claw,   // the conservation law
 				   const dgGroupOfFaces &group, 

Solver/dgAlgorithm.h

@@ -35,10 +35,22 @@ class dgAlgorithm {
       const fullMatrix<double> &solution, // solution !! at faces nodes
       fullMatrix<double> &residual // residual !! at faces nodes
       );
+	  
   void multAddInverseMassMatrix ( /*dofManager &dof,*/
       const dgGroupOfElements & group,
       fullMatrix<double> &residu,
       fullMatrix<double> &sol);
+
+  void rungeKutta (
+	  const dgConservationLaw &claw,
+      std::vector<dgGroupOfElements*> &eGroups, //group of elements
+      std::vector<dgGroupOfFaces*> &fGroups,  // group of interfacs
+      std::vector<dgGroupOfFaces*> &bGroups, // group of boundaries
+	  double h,	
+      fullMatrix<double> &residu,
+      fullMatrix<double> &sol,
+	  int orderRK=4);
+	  	  
   void buildGroups(GModel *model, int dim, int order,
       std::vector<dgGroupOfElements*> &eGroups,
       std::vector<dgGroupOfFaces*> &fGroups,

Solver/dgGroupOfElements.cpp

@@ -43,10 +43,10 @@ static fullMatrix<double> * dgGetFaceIntegrationRuleOnElement (
 
 dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e, int polyOrder)
   : _elements(e), 
-    _fs(*_elements[0]->getFunctionSpace(polyOrder)),
-    _integration(dgGetIntegrationRule (_elements[0], polyOrder)
-    )
+   _fs(*_elements[0]->getFunctionSpace(polyOrder)),
+   _integration(dgGetIntegrationRule (_elements[0], polyOrder))
 {
+
   _dimUVW = _dimXYZ = e[0]->getDim();
   // this is the biggest piece of data ... the mappings
   int nbNodes = _fs.coefficients.size1();

Solver/dgMainTest.cpp

@@ -36,7 +36,7 @@ class dgConservationLawInitialCondition : public dgConservationLaw {
 };
 
 int main(int argc, char **argv){
-  GmshMergeFile("square1.msh");
+  GmshMergeFile("input/square1.msh");
   //we probably need a class to group those three ones
   std::vector<dgGroupOfElements*> elementGroups;
   std::vector<dgGroupOfFaces*> faceGroups;
@@ -57,7 +57,8 @@ int main(int argc, char **argv){
     algo.residualVolume(initLaw,*elementGroups[0],sol,residu);
     algo.multAddInverseMassMatrix(*elementGroups[0],residu,sol);
   }
-  print("init.pos",*elementGroups[0],&sol(0,0));
+ 
+ 
 
 
   fullMatrix<double> advectionSpeed(3,1);
@@ -77,17 +78,21 @@ int main(int argc, char **argv){
   law->addBoundaryCondition("Top",dgBoundaryCondition::new0FluxCondition(*law));
   law->addBoundaryCondition("Bottom",dgBoundaryCondition::new0FluxCondition(*law));
 
+  
+  
+ print("output/init.pos",*elementGroups[0],&sol(0,0));
   for(int iT=0; iT<1000; iT++) {
     algo.residual(*law,elementGroups,faceGroups,boundaryGroups,sol,residu);
-    residu.scale(0.01);
-    algo.multAddInverseMassMatrix(*elementGroups[0],residu,sol);
+	algo.rungeKutta(*law,elementGroups,faceGroups,boundaryGroups,0.01,residu,sol);
     if(iT%10==0){
-      char name[100];
-      sprintf(name,"test_%05i.pos",iT/10);
+      char name[10];
+      sprintf(name,"output/test_%05i.pos",iT/10);
       printf("%i\n",iT);
       print(name,*elementGroups[0],&sol(0,0));
+
     }
   }
+
   delete law;
 }