From a97df3a00d592c21893f6c072c9767238c554c53 Mon Sep 17 00:00:00 2001
From: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date: Thu, 4 Feb 2010 09:10:44 +0000
Subject: [PATCH] dg:non linear shallow water equations
---
Solver/TESTCASES/Stommel.lua | 8 +-
Solver/dgConservationLawShallowWater2d.cpp | 85 +++++++++++++++-------
Solver/dgGroupOfElements.cpp | 2 +-
Solver/dgSlopeLimiter.cpp | 1 +
Solver/dgSystemOfEquations.cpp | 2 +-
5 files changed, 65 insertions(+), 33 deletions(-)
diff --git a/Solver/TESTCASES/Stommel.lua b/Solver/TESTCASES/Stommel.lua
index b7b1ef93c7..4342a7e812 100644
--- a/Solver/TESTCASES/Stommel.lua
+++ b/Solver/TESTCASES/Stommel.lua
@@ -1,7 +1,7 @@
model = GModel()
model:load ('stommel_square.msh')
dg = dgSystemOfEquations (model)
-order=1
+order=3
dg:setOrder (order)
claw = dgConservationLawShallowWater2d()
claw:addBoundaryCondition('Wall',claw:newBoundaryWall())
@@ -10,12 +10,12 @@ dg:setup()
dg:exportSolution('output/init')
-for i=1,10000 do
- norm = dg:RK44(80*(3/(2.*order+1)))
+for i=1,100000 do
+ norm = dg:RK44(150*(3/(2.*order+1)))
if ( i%100 ==0 ) then
print ('iter ', i, norm)
end
if ( i%100 ==0 ) then
- dg:exportSolution(string.format('output/solution-%05d',i))
+ dg:exportSolution(string.format('output/solution-%06d',i))
end
end
diff --git a/Solver/dgConservationLawShallowWater2d.cpp b/Solver/dgConservationLawShallowWater2d.cpp
index 0671ba0217..4873c3b9d0 100644
--- a/Solver/dgConservationLawShallowWater2d.cpp
+++ b/Solver/dgConservationLawShallowWater2d.cpp
@@ -1,31 +1,30 @@
#include "dgConservationLawShallowWater2d.h"
#include "function.h"
static double g = 9.81;
-static double h = 1000;
+static double h = 1000.;
+static double linearDissipation = 2e-7;
class dgConservationLawShallowWater2d::advection: public dataCacheDouble {
dataCacheDouble /
public:
advection(dataCacheMap &cacheMap):
- dataCacheDouble(cacheMap),
+ dataCacheDouble(cacheMap,cacheMap.getNbEvaluationPoints(),9),
sol(cacheMap.get("Solution",this))
{};
void _eval () {
int nQP = sol().size1();
- if(_value.size1() != nQP)
- _value=fullMatrix<double>(nQP,9);
for(int i=0; i< nQP; i++) {
double eta = sol(i,0);
double u = sol(i,1);
double v = sol(i,2);
// flux_x
- _value(i,0) = h*u;
- _value(i,1) = g*eta;
- _value(i,2) = 0;
+ _value(i,0) = (h+eta)*u;
+ _value(i,1) = g*eta+u*u;
+ _value(i,2) = u*v;
// flux_y
- _value(i,3) = h*v;
- _value(i,4) = 0;
- _value(i,5) = g*eta;
+ _value(i,3) = (h+eta)*v;
+ _value(i,4) = v*u;
+ _value(i,5) = g*eta+v*v;
// flux_z
_value(i,6) = 0;
_value(i,7) = 0;
@@ -35,27 +34,29 @@ class dgConservationLawShallowWater2d::advection: public dataCacheDouble {
};
class dgConservationLawShallowWater2d::source: public dataCacheDouble {
- dataCacheDouble &xyz, &solution;
+ dataCacheDouble &xyz, &solution,&solutionGradient;
public :
source(dataCacheMap &cacheMap) :
- dataCacheDouble(cacheMap),
+ dataCacheDouble(cacheMap,cacheMap.getNbEvaluationPoints(),3),
xyz(cacheMap.get("XYZ",this)),
- solution(cacheMap.get("Solution",this))
+ solution(cacheMap.get("Solution",this)),
+ solutionGradient(cacheMap.get("SolutionGradient",this))
{}
void _eval () {
- int nQP = xyz().size1();
- if(_value.size1() != nQP)
- _value = fullMatrix<double>(nQP,3);
+ int nQP =_value.size1();
for (int i=0; i<nQP; i++) {
double eta = solution(i,0);
double u = solution(i,1);
double v = solution(i,2);
+ double dudx = solutionGradient(i*3,1);
+ double dvdx = solutionGradient(i*3,2);
+ double dudy = solutionGradient(i*3+1,1);
+ double dvdy = solutionGradient(i*3+1,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;
_value (i,0) = 0;
- _value (i,1) = f*v + - gamma*u + wind;
- _value (i,2) = -f*u - gamma*v ;
+ _value (i,1) = f*v + - linearDissipation*u + wind + u*(dudx+dvdy);
+ _value (i,2) = -f*u - linearDissipation*v + v*(dudx+dvdy);
}
}
};
@@ -63,33 +64,63 @@ class dgConservationLawShallowWater2d::riemann:public dataCacheDouble {
dataCacheDouble &normals, &solL, &solR;
public:
riemann(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight):
- dataCacheDouble(cacheMapLeft),
+ dataCacheDouble(cacheMapLeft,cacheMapLeft.getNbEvaluationPoints(),6),
normals(cacheMapLeft.get("Normals", this)),
solL(cacheMapLeft.get("Solution", this)),
solR(cacheMapRight.get("Solution", this))
{};
void _eval () {
int nQP = solL().size1();
- if(_value.size1() != nQP)
- _value = fullMatrix<double>(nQP,6);
for(int i=0; i< nQP; i++) {
double nx = normals(0,i);
double ny = normals(1,i);
+ double HR = solR(i,0) + h, HL = solL(i,0) + h;
+ double uL = nx*solL(i,1) + ny*solL(i,2), uR = nx*solR(i,1) + ny*solR(i,2);
+ double vL = -ny*solL(i,1) + nx*solL(i,2), vR = - ny*solR(i,1) + nx*solR(i,2);
+ double HuL = uL*HL, HuR = uR*HR;
+ double HvL = vL*HL, HvR = vR*HR;
+ double HM = (HL+HR)/2, HJ = (HL-HR)/2;
+ double HuM = (HuL+HuR)/2, HuJ = (HuL-HuR)/2;
+ double HvM = (HvL+HvR)/2, HvJ = (HvL-HvR)/2;
+ double sqHL = sqrt(HL), sqHR = sqrt(HR);
+ double u_roe = (sqHL*uL + sqHR*uR) / (sqHL + sqHR);
+ double v_roe = (sqHL*vL + sqHR*vR) / (sqHL + sqHR);
+ double c_roe = sqrt(g*HM);
+ double H = HM + (HuJ - u_roe *HJ) / c_roe;
+ double Hu = HuM + (c_roe - u_roe*u_roe/c_roe) *HJ + u_roe/c_roe *HuJ;
+ double Hv = -v_roe*u_roe/c_roe*HJ + v_roe/c_roe*HuJ + (u_roe>0 ? -v_roe*HJ+HvL : v_roe*HJ+HvR);
+ double u = Hu / H;
+ double v = Hv / H;
+ double eta = H-h;
+
+ /*double u = (uL+uR)/2; //
+ double v = (vL+vR)/2;
+ double eta = (HL+HR)/2-h;
+ double c = hypot(u,v)+sqrt(g*(h+eta));
+ double Fun = -g*eta - u*u - (uL-uR)/2 * c;
+ double Fut = -u*( u>0 ? vL:vR);
+ double Feta = -(h+eta)*u - (HL-HR)/2 * c;*/
+
+ double Fun = -g*eta - u*u;
+ double Fut = -u*v;
+ double Feta = -(h+eta)*u;
+
+ /* //linear equations
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 Fut = 0;
double Feta = -h*unRiemann;
- _value(i,0) = Feta;
- _value(i,1) = Fun*nx;
- _value(i,2) = Fun*ny;
+ */
+ _value(i,0) = Feta;
+ _value(i,1) = Fun*nx - Fut*ny;
+ _value(i,2) = Fun*ny + Fut*nx;
_value(i,3) = -_value(i,0);
_value(i,4) = -_value(i,1);
_value(i,5) = -_value(i,2);
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index c897f27887..96ada631f2 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -940,7 +940,7 @@ void dgGroupCollection::buildGroups(GModel *model, int dim, int order,std::strin
Msg::Info("%d groups of elements",localElements.size());
for (int n=0; n<localElements.size(); n++) {
- printf("\n **** Group %d: %d elements **** \n", n,localElements[n].size());
+ printf("\n **** Group %i: %lu elements **** \n", n,localElements[n].size());
}
// do a group of elements for every element type in the mesh
for (std::map<int, std::vector<MElement *> >::iterator it = localElements.begin(); it !=localElements.end() ; ++it){
diff --git a/Solver/dgSlopeLimiter.cpp b/Solver/dgSlopeLimiter.cpp
index 6f1009f33b..f7c0ae70d1 100644
--- a/Solver/dgSlopeLimiter.cpp
+++ b/Solver/dgSlopeLimiter.cpp
@@ -6,6 +6,7 @@
//----------------------------------------------------------------------------------
bool dgSlopeLimiter::apply ( dgDofContainer &solution, dgGroupCollection &groups)
{
+ printf("limit \n");
solution.scatter();
int nbFields =_claw->nbFields();
diff --git a/Solver/dgSystemOfEquations.cpp b/Solver/dgSystemOfEquations.cpp
index 29af76404b..89189337d9 100644
--- a/Solver/dgSystemOfEquations.cpp
+++ b/Solver/dgSystemOfEquations.cpp
@@ -215,7 +215,7 @@ void dgSystemOfEquations::export_solution_as_is (const std::string &name){
fullMatrix<double> sol = getSolutionProxy (i, iElement);
fprintf(f,"%d %d",num,sol.size1());
for (int k=0;k<sol.size1();++k) {
- fprintf(f," %12.5E ",sol(k,ICOMP));
+ fprintf(f," %.16E ",sol(k,ICOMP));
}
fprintf(f,"\n");
}
--
GitLab