Skip to content
Snippets Groups Projects
Commit 4be90823 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

dg : inviscid non conservative linear shallow water equation

parent 30fa96cc
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
Solver/dgConservationLaw.h
+ 1
− 0
View file @ 4be90823
...@@ -55,5 +55,6 @@ dgConservationLaw *dgNewConservationLawAdvection(const std::string vname); ...@@ -55,5 +55,6 @@ dgConservationLaw *dgNewConservationLawAdvection(const std::string vname);
dgConservationLaw *dgNewConservationLawShallowWater2d(); dgConservationLaw *dgNewConservationLawShallowWater2d();
dgConservationLaw *dgNewConservationLawWaveEquation(); dgConservationLaw *dgNewConservationLawWaveEquation();
dgBoundaryCondition *dgNewBoundaryConditionWaveEquationWall(); dgBoundaryCondition *dgNewBoundaryConditionWaveEquationWall();
dgBoundaryCondition *dgNewBoundaryConditionShallowWater2dWall();
#endif #endif
Solver/dgConservationLawShallowWater2d.cpp
+ 64
− 32
View file @ 4be90823
#include "dgConservationLaw.h" #include "dgConservationLaw.h"
#include "function.h" #include "function.h"
static double g = 9.81; static double g = 9.81;
static double h = 1000;
class dgConservationLawShallowWater2d : public dgConservationLaw { class dgConservationLawShallowWater2d : public dgConservationLaw {
class advection : public dataCacheDouble { class advection : public dataCacheDouble {
dataCacheDouble / dataCacheDouble /
...@@ -13,23 +15,22 @@ class dgConservationLawShallowWater2d : public dgConservationLaw { ...@@ -13,23 +15,22 @@ class dgConservationLawShallowWater2d : public dgConservationLaw {
if(_value.size1() != nQP) if(_value.size1() != nQP)
_value=fullMatrix<double>(nQP,9); _value=fullMatrix<double>(nQP,9);
for(int i=0; i< nQP; i++) { for(int i=0; i< nQP; i++) {
double H = sol(i,0); double eta = sol(i,0);
double U = sol(i,1); double u = sol(i,1);
double V = sol(i,2); double v = sol(i,2);
// flux_x // flux_x
_value(i,0) = U; _value(i,0) = h*u;
_value(i,1) = U*U/H + g*H*H/2; _value(i,1) = g*eta;
_value(i,2) = U*V/H; _value(i,2) = 0;
// flux_y // flux_y
_value(i,3) = V; _value(i,3) = h*v;
_value(i,4) = V*U/H; _value(i,4) = 0;
_value(i,5) = V*V/H + g*H*H/2; _value(i,5) = g*eta;
// flux_z // flux_z
_value(i,6) = 0; _value(i,6) = 0;
_value(i,7) = 0; _value(i,7) = 0;
_value(i,8) = 0; _value(i,8) = 0;
} }
_value.scale(0);
} }
}; };
class source : public dataCacheDouble { class source : public dataCacheDouble {
...@@ -44,16 +45,15 @@ class dgConservationLawShallowWater2d : public dgConservationLaw { ...@@ -44,16 +45,15 @@ class dgConservationLawShallowWater2d : public dgConservationLaw {
if(_value.size1() != nQP) if(_value.size1() != nQP)
_value = fullMatrix<double>(nQP,3); _value = fullMatrix<double>(nQP,3);
for (int i=0; i<nQP; i++) { for (int i=0; i<nQP; i++) {
double H = solution(i,0); double eta = solution(i,0);
double U = solution(i,1); double u = solution(i,1);
double V = solution(i,2); double v = solution(i,2);
double wind = 0.1*sin(xyz(i,1)/1e6*M_PI)/H/1000; double wind = 0.1*sin(xyz(i,1)/1e6*M_PI)/1000/h;
double f = 1e-4+xyz(i,1)*2e-11; double f = 1e-4+xyz(i,1)*2e-11;
double gamma = 1e-6; double gamma = 1e-6;
f=0;
_value (i,0) = 0; _value (i,0) = 0;
_value (i,1) = -gamma*U - f*V + wind; _value (i,1) = f*v + - gamma*u + wind;
_value (i,2) = -gamma*V + f*U; _value (i,2) = -f*u - gamma*v ;
} }
} }
}; };
...@@ -72,22 +72,21 @@ class dgConservationLawShallowWater2d : public dgConservationLaw { ...@@ -72,22 +72,21 @@ class dgConservationLawShallowWater2d : public dgConservationLaw {
for(int i=0; i< nQP; i++) { for(int i=0; i< nQP; i++) {
double nx = normals(0,i); double nx = normals(0,i);
double ny = normals(1,i); double ny = normals(1,i);
double UnL = nx*solL(i,1) + ny*solL(i,2); double unL = nx*solL(i,1) + ny*solL(i,2);
double UnR = nx*solR(i,1) + ny*solR(i,2); double unR = nx*solR(i,1) + ny*solR(i,2);
double UtL = ny*solL(i,1) - nx*solL(i,2); double utL = ny*solL(i,1) - nx*solL(i,2);
double UtR = ny*solR(i,1) - nx*solR(i,2); double utR = ny*solR(i,1) - nx*solR(i,2);
double HR = solR(i,0); double etaR = solR(i,0);
double HL = solL(i,0); double etaL = solL(i,0);
double c = sqrt(g*(HR+HL)/2); double sq_g_h = sqrt(g/h);
double FUn = -(UnL*UnL/HL + g*HL*HL/2 + UnR*UnR/HR + g*HR*HR/2)/2 + c*(UnL-UnR)/2; double etaRiemann = (etaL+etaR + (unL-unR)/sq_g_h)/2;
double FUt = -(UnL*UtL/HL + UnR*UtR/HR)/2 + c*(UtL-UtR)/2; double unRiemann = (unL+unR + (etaL-etaR)*sq_g_h)/2;
double FH = -(UnL+UnR)/2 + c*(HL-HR)/2; double Fun = -g*etaRiemann;
_value(i,0) = FH; double Feta = -h*unRiemann;
_value(i,1) = FUn*nx + FUt*ny; _value(i,0) = Feta;
_value(i,2) = FUn*ny - FUt*nx; _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,3) = -_value(i,0);
_value(i,4) = -_value(i,1); _value(i,4) = -_value(i,1);
_value(i,5) = -_value(i,2); _value(i,5) = -_value(i,2);
...@@ -113,6 +112,39 @@ class dgConservationLawShallowWater2d : public dgConservationLaw { ...@@ -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() { dgConservationLaw *dgNewConservationLawShallowWater2d() {
return new dgConservationLawShallowWater2d(); return new dgConservationLawShallowWater2d();
} }
Solver/dgConservationLawWaveEquation.cpp
+ 2
− 2
View file @ 4be90823
...@@ -55,8 +55,8 @@ class dgConservationLawWaveEquation : public dgConservationLaw { ...@@ -55,8 +55,8 @@ class dgConservationLawWaveEquation : public dgConservationLaw {
double pR = solR(i,0); double pR = solR(i,0);
double pL = solL(i,0); double pL = solL(i,0);
double pRiemann = 0.5 * (pL+pR - (unR-unL)*(rho*c)); double pRiemann = 0.5 * (pL+pR + (unL-unR)*(rho*c));
double unRiemann = 0.5 * (unL+unR - (pR-pL)/(rho*c)); double unRiemann = 0.5 * (unL+unR + (pL-pR)/(rho*c));
double Fp = -rho*c*c*unRiemann; double Fp = -rho*c*c*unRiemann;
double Fun = -pRiemann/rho; double Fun = -pRiemann/rho;
... ...
......
Solver/dgMainShallowWater2d.cpp
+ 7
− 119
View file @ 4be90823
...@@ -11,108 +11,6 @@ ...@@ -11,108 +11,6 @@
#include "MElement.h" #include "MElement.h"
void print (const char *filename,const dgGroupOfElements &els, double *v,int iField=1, int nbFields=1); 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) { int main(int argc, char **argv) {
GmshMergeFile("input/stommel.msh"); GmshMergeFile("input/stommel.msh");
//we probably need a class to group those three ones //we probably need a class to group those three ones
...@@ -129,32 +27,22 @@ int main(int argc, char **argv) { ...@@ -129,32 +27,22 @@ int main(int argc, char **argv) {
int nbNodes = elementGroups[0]->getNbNodes(); int nbNodes = elementGroups[0]->getNbNodes();
fullMatrix<double> sol(nbNodes,elementGroups[0]->getNbElements()*3); fullMatrix<double> sol(nbNodes,elementGroups[0]->getNbElements()*3);
fullMatrix<double> residu(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(); dgConservationLaw *law = dgNewConservationLawShallowWater2d();
law->addBoundaryCondition("Wall",dgNewBoundaryConditionShallowWater2dWall());
law->addBoundaryCondition("Top",dgBoundaryCondition::new0FluxCondition(*law)); for(int iT=0; iT<10000; iT++) {
law->addBoundaryCondition("Bottom",dgBoundaryCondition::new0FluxCondition(*law)); if(iT%100==0){
law->addBoundaryCondition("Left",dgBoundaryCondition::new0FluxCondition(*law));
law->addBoundaryCondition("Right",dgBoundaryCondition::new0FluxCondition(*law));
for(int iT=0; iT<100; iT++) {
if(iT%10==0){
printf("%i\n",iT); printf("%i\n",iT);
char name[100]; 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); 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); 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); 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; delete law;
} }
... ...
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment