From 097c6c1651956268ea0953f8a1b0e69635f78853 Mon Sep 17 00:00:00 2001
From: Bruno Seny <bruno.seny@student.uclouvain.be>
Date: Tue, 24 Nov 2009 10:44:24 +0000
Subject: [PATCH] Runge-Kutta
---
Solver/dgAlgorithm.cpp | 46 ++++++++++++++++++++++++++++++++++--
Solver/dgAlgorithm.h | 12 ++++++++++
Solver/dgGroupOfElements.cpp | 6 ++---
Solver/dgMainTest.cpp | 17 ++++++++-----
4 files changed, 70 insertions(+), 11 deletions(-)
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index eb97d9578a..9edb3a4ca5 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/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,
diff --git a/Solver/dgAlgorithm.h b/Solver/dgAlgorithm.h
index 9a0c3c96b1..91061c7ef3 100644
--- a/Solver/dgAlgorithm.h
+++ b/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,
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index abbb35e326..5b28f8acbe 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/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();
diff --git a/Solver/dgMainTest.cpp b/Solver/dgMainTest.cpp
index 40e33204e4..03f4ee006d 100644
--- a/Solver/dgMainTest.cpp
+++ b/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;
}
--
GitLab