From 4be90823c72c7068fe0cc48bec893551cf8cb437 Mon Sep 17 00:00:00 2001
From: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date: Wed, 2 Dec 2009 10:49:54 +0000
Subject: [PATCH] dg : inviscid non conservative linear shallow water equation
---
Solver/dgConservationLaw.h | 1 +
Solver/dgConservationLawShallowWater2d.cpp | 96 ++++++++++------
Solver/dgConservationLawWaveEquation.cpp | 4 +-
Solver/dgMainShallowWater2d.cpp | 126 ++-------------------
4 files changed, 74 insertions(+), 153 deletions(-)
diff --git a/Solver/dgConservationLaw.h b/Solver/dgConservationLaw.h
index f2370de706..93df46b575 100644
--- a/Solver/dgConservationLaw.h
+++ b/Solver/dgConservationLaw.h
@@ -55,5 +55,6 @@ dgConservationLaw *dgNewConservationLawAdvection(const std::string vname);
dgConservationLaw *dgNewConservationLawShallowWater2d();
dgConservationLaw *dgNewConservationLawWaveEquation();
dgBoundaryCondition *dgNewBoundaryConditionWaveEquationWall();
+dgBoundaryCondition *dgNewBoundaryConditionShallowWater2dWall();
#endif
diff --git a/Solver/dgConservationLawShallowWater2d.cpp b/Solver/dgConservationLawShallowWater2d.cpp
index cb2caa2052..dbf70dded9 100644
--- a/Solver/dgConservationLawShallowWater2d.cpp
+++ b/Solver/dgConservationLawShallowWater2d.cpp
@@ -1,6 +1,8 @@
#include "dgConservationLaw.h"
#include "function.h"
static double g = 9.81;
+static double h = 1000;
+
class dgConservationLawShallowWater2d : public dgConservationLaw {
class advection : public dataCacheDouble {
dataCacheDouble /
@@ -13,23 +15,22 @@ class dgConservationLawShallowWater2d : public dgConservationLaw {
if(_value.size1() != nQP)
_value=fullMatrix<double>(nQP,9);
for(int i=0; i< nQP; i++) {
- double H = sol(i,0);
- double U = sol(i,1);
- double V = sol(i,2);
+ double eta = sol(i,0);
+ double u = sol(i,1);
+ double v = sol(i,2);
// flux_x
- _value(i,0) = U;
- _value(i,1) = U*U/H + g*H*H/2;
- _value(i,2) = U*V/H;
+ _value(i,0) = h*u;
+ _value(i,1) = g*eta;
+ _value(i,2) = 0;
// flux_y
- _value(i,3) = V;
- _value(i,4) = V*U/H;
- _value(i,5) = V*V/H + g*H*H/2;
+ _value(i,3) = h*v;
+ _value(i,4) = 0;
+ _value(i,5) = g*eta;
// flux_z
_value(i,6) = 0;
_value(i,7) = 0;
_value(i,8) = 0;
}
- _value.scale(0);
}
};
class source : public dataCacheDouble {
@@ -44,16 +45,15 @@ class dgConservationLawShallowWater2d : public dgConservationLaw {
if(_value.size1() != nQP)
_value = fullMatrix<double>(nQP,3);
for (int i=0; i<nQP; i++) {
- double H = solution(i,0);
- double U = solution(i,1);
- double V = solution(i,2);
- double wind = 0.1*sin(xyz(i,1)/1e6*M_PI)/H/1000;
+ double eta = solution(i,0);
+ double u = solution(i,1);
+ double v = solution(i,2);
+ double wind = 0.1*sin(xyz(i,1)/1e6*M_PI)/1000/h;
double f = 1e-4+xyz(i,1)*2e-11;
double gamma = 1e-6;
- f=0;
_value (i,0) = 0;
- _value (i,1) = -gamma*U - f*V + wind;
- _value (i,2) = -gamma*V + f*U;
+ _value (i,1) = f*v + - gamma*u + wind;
+ _value (i,2) = -f*u - gamma*v ;
}
}
};
@@ -72,22 +72,21 @@ class dgConservationLawShallowWater2d : public dgConservationLaw {
for(int i=0; i< nQP; i++) {
double nx = normals(0,i);
double ny = normals(1,i);
- double UnL = nx*solL(i,1) + ny*solL(i,2);
- double UnR = nx*solR(i,1) + ny*solR(i,2);
- double UtL = ny*solL(i,1) - nx*solL(i,2);
- double UtR = ny*solR(i,1) - nx*solR(i,2);
- double HR = solR(i,0);
- double HL = solL(i,0);
- double c = sqrt(g*(HR+HL)/2);
- double FUn = -(UnL*UnL/HL + g*HL*HL/2 + UnR*UnR/HR + g*HR*HR/2)/2 + c*(UnL-UnR)/2;
- double FUt = -(UnL*UtL/HL + UnR*UtR/HR)/2 + c*(UtL-UtR)/2;
- double FH = -(UnL+UnR)/2 + c*(HL-HR)/2;
- _value(i,0) = FH;
- _value(i,1) = FUn*nx + FUt*ny;
- _value(i,2) = FUn*ny - FUt*nx;
+ double unL = nx*solL(i,1) + ny*solL(i,2);
+ double unR = nx*solR(i,1) + ny*solR(i,2);
+ double utL = ny*solL(i,1) - nx*solL(i,2);
+ double utR = ny*solR(i,1) - nx*solR(i,2);
+ double etaR = solR(i,0);
+ double etaL = solL(i,0);
+ double sq_g_h = sqrt(g/h);
+ double etaRiemann = (etaL+etaR + (unL-unR)/sq_g_h)/2;
+ double unRiemann = (unL+unR + (etaL-etaR)*sq_g_h)/2;
+ double Fun = -g*etaRiemann;
+ double Feta = -h*unRiemann;
+ _value(i,0) = Feta;
+ _value(i,1) = Fun*nx;
+ _value(i,2) = Fun*ny;
- _value(i,1)=0;
- _value(i,2)=0;
_value(i,3) = -_value(i,0);
_value(i,4) = -_value(i,1);
_value(i,5) = -_value(i,2);
@@ -113,6 +112,39 @@ class dgConservationLawShallowWater2d : public dgConservationLaw {
}
};
+class dgBoundaryConditionShallowWater2dWall : public dgBoundaryCondition {
+ class term : public dataCacheDouble {
+ dataCacheDouble &sol,&normals;
+ public:
+ term(dataCacheMap &cacheMap):
+ sol(cacheMap.get("Solution",this)),
+ normals(cacheMap.get("Normals",this)){}
+ void _eval () {
+ int nQP = sol().size1();
+ if(_value.size1() != nQP)
+ _value = fullMatrix<double>(nQP,3);
+ for(int i=0; i< nQP; i++) {
+ double nx = normals(0,i);
+ double ny = normals(1,i);
+ double eta = sol(i,0);
+
+ _value(i,0) = 0;
+ _value(i,1) = -g*eta*nx;
+ _value(i,2) = -g*eta*ny;
+ }
+ }
+ };
+ public:
+ dgBoundaryConditionShallowWater2dWall(){}
+ dataCacheDouble *newBoundaryTerm(dataCacheMap &cacheMapLeft) const {
+ return new term(cacheMapLeft);
+ }
+};
+
+dgBoundaryCondition *dgNewBoundaryConditionShallowWater2dWall() {
+ return new dgBoundaryConditionShallowWater2dWall();
+}
+
dgConservationLaw *dgNewConservationLawShallowWater2d() {
return new dgConservationLawShallowWater2d();
}
diff --git a/Solver/dgConservationLawWaveEquation.cpp b/Solver/dgConservationLawWaveEquation.cpp
index cb9223c232..b5d2e457fb 100644
--- a/Solver/dgConservationLawWaveEquation.cpp
+++ b/Solver/dgConservationLawWaveEquation.cpp
@@ -55,8 +55,8 @@ class dgConservationLawWaveEquation : public dgConservationLaw {
double pR = solR(i,0);
double pL = solL(i,0);
- double pRiemann = 0.5 * (pL+pR - (unR-unL)*(rho*c));
- double unRiemann = 0.5 * (unL+unR - (pR-pL)/(rho*c));
+ double pRiemann = 0.5 * (pL+pR + (unL-unR)*(rho*c));
+ double unRiemann = 0.5 * (unL+unR + (pL-pR)/(rho*c));
double Fp = -rho*c*c*unRiemann;
double Fun = -pRiemann/rho;
diff --git a/Solver/dgMainShallowWater2d.cpp b/Solver/dgMainShallowWater2d.cpp
index efb8462689..e3aed3f67d 100644
--- a/Solver/dgMainShallowWater2d.cpp
+++ b/Solver/dgMainShallowWater2d.cpp
@@ -11,108 +11,6 @@
#include "MElement.h"
void print (const char *filename,const dgGroupOfElements &els, double *v,int iField=1, int nbFields=1);
-class dgConservationLawInitialCondition : public dgConservationLaw {
- class gaussian : public dataCacheDouble {
- dataCacheDouble &xyz;
- double _xc,_yc;
- public:
- gaussian(dataCacheMap &cacheMap,double xc, double yc):xyz(cacheMap.get("XYZ",this)),_xc(xc),_yc(yc){};
- void _eval () {
- if(_value.size1() != xyz().size1())
- _value=fullMatrix<double>(xyz().size1(),1);
- for(int i=0; i< _value.size1(); i++) {
- _value(i,0)=exp(-(pow(xyz(i,0)-_xc,2)+pow(xyz(i,1)-_yc,2))*100);
- }
- }
- };
- public:
- dgConservationLawInitialCondition() {
- _nbf = 1;
- }
- dataCacheDouble *newSourceTerm(dataCacheMap &cacheMap)const {
- return new gaussian(cacheMap,0.2,0.3);
- }
-};
-
-/*int main(int argc, char **argv) {
- GmshMergeFile("input/square1.msh");
- //we probably need a class to group those three ones
- std::vector<dgGroupOfElements*> elementGroups;
- std::vector<dgGroupOfFaces*> faceGroups;
- std::vector<dgGroupOfFaces*> boundaryGroups;
- int order=1;
- int dimension=2;
- dgAlgorithm algo;
- function::registerDefaultFunctions();
- algo.buildGroups(GModel::current(),dimension,order,elementGroups,faceGroups,boundaryGroups);
-
- //for now, we suppose there is only one group of elements
- int nbNodes = elementGroups[0]->getNbNodes();
- fullMatrix<double> sol(nbNodes,elementGroups[0]->getNbElements());
- fullMatrix<double> residu(nbNodes,elementGroups[0]->getNbElements());
- // initial condition
- {
- dgConservationLawInitialCondition initLaw;
- algo.residualVolume(initLaw,*elementGroups[0],sol,residu);
- algo.multAddInverseMassMatrix(*elementGroups[0],residu,sol);
- }
-
- fullMatrix<double> advectionSpeed(3,1);
- advectionSpeed(0,0)=0.15;
- advectionSpeed(1,0)=0.05;
- advectionSpeed(2,0)=0.;
- function::add("advectionSpeed",function::newFunctionConstant(advectionSpeed));
-
- dgConservationLaw *law = dgNewConservationLawAdvection("advectionSpeed");
-
- fullMatrix<double> outsideValue(1,1);
- outsideValue(0,0)=0;
- function::add("outsideValue",function::newFunctionConstant(outsideValue));
- law->addBoundaryCondition("Left",dgBoundaryCondition::newOutsideValueCondition(*law,"outsideValue"));
- law->addBoundaryCondition("Right",dgBoundaryCondition::newOutsideValueCondition(*law,"outsideValue"));
-
- 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<100; iT++) {
- algo.rungeKutta(*law,elementGroups,faceGroups,boundaryGroups,0.01,residu,sol);
- if(iT%10==0){
- char name[100];
- sprintf(name,"output/test_%05i.pos",iT/10);
- printf("%i\n",iT);
- print(name,*elementGroups[0],&sol(0,0));
- }
- }
- delete law;
-}*/
-class dgConservationLawStommel : public dgConservationLaw {
- class initbath : public dataCacheDouble {
- dataCacheDouble &uvw;
- public:
- initbath(dataCacheMap &cacheMap):
- uvw(cacheMap.get("UVW"))
- {}
- void _eval () {
- if(_value.size1() != uvw().size1())
- _value=fullMatrix<double>(uvw().size1(),3);
- for(int i=0; i< _value.size1(); i++) {
- _value(i,0)=1000;
- _value(i,1)=0;
- _value(i,2)=0;
- }
- }
- };
- public:
- dgConservationLawStommel() {
- _nbf = 3;
- }
- dataCacheDouble *newSourceTerm(dataCacheMap &cacheMap)const {
- return new initbath(cacheMap);
- }
-};
int main(int argc, char **argv) {
GmshMergeFile("input/stommel.msh");
//we probably need a class to group those three ones
@@ -129,32 +27,22 @@ int main(int argc, char **argv) {
int nbNodes = elementGroups[0]->getNbNodes();
fullMatrix<double> sol(nbNodes,elementGroups[0]->getNbElements()*3);
fullMatrix<double> residu(nbNodes,elementGroups[0]->getNbElements()*3);
- // initial condition
- {
- dgConservationLawStommel initLaw;
- algo.residualVolume(initLaw,*elementGroups[0],sol,residu);
- algo.multAddInverseMassMatrix(*elementGroups[0],residu,sol);
- }
dgConservationLaw *law = dgNewConservationLawShallowWater2d();
-
- law->addBoundaryCondition("Top",dgBoundaryCondition::new0FluxCondition(*law));
- law->addBoundaryCondition("Bottom",dgBoundaryCondition::new0FluxCondition(*law));
- law->addBoundaryCondition("Left",dgBoundaryCondition::new0FluxCondition(*law));
- law->addBoundaryCondition("Right",dgBoundaryCondition::new0FluxCondition(*law));
+ law->addBoundaryCondition("Wall",dgNewBoundaryConditionShallowWater2dWall());
- for(int iT=0; iT<100; iT++) {
- if(iT%10==0){
+ for(int iT=0; iT<10000; iT++) {
+ if(iT%100==0){
printf("%i\n",iT);
char name[100];
- sprintf(name,"output/H_%05i.pos",iT/10);
+ sprintf(name,"output/eta_%05i.pos",iT/100);
print(name,*elementGroups[0],&sol(0,0),0,3);
- sprintf(name,"output/U_%05i.pos",iT/10);
+ sprintf(name,"output/u_%05i.pos",iT/100);
print(name,*elementGroups[0],&sol(0,0),1,3);
- sprintf(name,"output/V_%05i.pos",iT/10);
+ sprintf(name,"output/v_%05i.pos",iT/100);
print(name,*elementGroups[0],&sol(0,0),2,3);
}
- algo.rungeKutta(*law,elementGroups,faceGroups,boundaryGroups,1e-8,residu,sol);
+ algo.rungeKutta(*law,elementGroups,faceGroups,boundaryGroups,50,residu,sol);
}
delete law;
}
--
GitLab