dgConservationLawShallowWaterTracer2d.cpp 10 KB
Newer Older
lambrechts's avatar
lambrechts committed
1
2
3
4
5
6
7
8
9
#include "dgConservationLawShallowWaterTracer2d.h"
#include "dgConservationLawShallowWater2d.h"
#include "dgNeighbourMap.h"
#include "function.h"
#include "Numeric.h" 
#include "assert.h" 
#include "stdlib.h"
#include <string.h>
#include <algorithm>
10
#include "functorMember.h"
Valentin Vallaeys's avatar
Valentin Vallaeys committed
11
#include "slimMovingBathWettingDrying.h"
lambrechts's avatar
lambrechts committed
12
13
14
15
16


//*************************************
// Terms of the equation
//***************************************
17
18
19
20
21
22
23

void dgConservationLawShallowWaterTracer2d::massFactorCB(functorCache &m, fullMatrix<double> &val) const {
  const fullMatrix<double> &hydroSol = m.get(*_hydroSolution);
  const fullMatrix<double> &bath = m.get(*_bathymetry);
  val.resize(m.nEvaluationPoint(), 1);
  for(int i=0; i< m.nEvaluationPoint(); i++) {
    val(i,0) = bath(i,0) + hydroSol(i,0);
lambrechts's avatar
lambrechts committed
24
  }
25
26
27
28
29
30
}

void dgConservationLawShallowWaterTracer2d::gradPsiTermCB(functorCache &m, fullMatrix<double> &val)const{
  const fullMatrix<double> &hydroSol = m.get(*_hydroSolution);
  const fullMatrix<double> &bath = m.get(*_bathymetry);
  const fullMatrix<double> &sol = m.get(m.solution());
31
  const fullMatrix<double> *diffusiveFlux = _diffusiveFlux ? &m.get(*_diffusiveFlux) : NULL;
32
33
34
35
36
37
  val.resize(m.nEvaluationPoint(), 3);
  for(int i=0; i< m.nEvaluationPoint(); i++) {
    double eta=hydroSol(i,0);
    double H = bath(i,0) + eta;
    double u=hydroSol(i,1);
    double v=hydroSol(i,2);
38
39
    val(i,0) = H * u * sol(i,0);
    val(i,1) = H * v * sol(i,0);
40
    val(i,2) = 0.;
41
42
43
44
    if (diffusiveFlux){
      val(i,0) += (*diffusiveFlux)(i,0);
      val(i,1) += (*diffusiveFlux)(i,1);
    }
lambrechts's avatar
lambrechts committed
45
  }
46
47
48
49
50
51
52
53
54
}

void dgConservationLawShallowWaterTracer2d::diffusiveFluxCB(functorCache &m, fullMatrix<double> &val) const{
  const fullMatrix<double> &kappa = m.get(*_kappa);
  const fullMatrix<double> &bath = m.get(*_bathymetry);
  const fullMatrix<double> &hydroSol = m.get(*_hydroSolution);
  const fullMatrix<double> &solGrad = m.get(m.solutionGradient());
  val.resize(m.nEvaluationPoint(), 2);
  if (m.interfaceGroup() && m.interfaceGroup()->nConnection() == 1 && m.side() == 1) {
lambrechts's avatar
lambrechts committed
55
  }
56
57
58
59
  for(int i = 0; i < m.nEvaluationPoint(); i++) {
    double H = hydroSol(i,0) + bath(i,0);
    val(i,0) = -kappa(i,0) * H * solGrad(i,0);
    val(i,1) = -kappa(i,0) * H * solGrad(i,1);
lambrechts's avatar
lambrechts committed
60
  }
61
62
63
64
65
66
67
68
69
70
}

void dgConservationLawShallowWaterTracer2d::psiTermCB(functorCache &m, fullMatrix<double> &val) const{
  const fullMatrix<double> &hydroSol = m.get(*_hydroSolution);
  const fullMatrix<double> &source = m.get(*_source);
  const fullMatrix<double> &bath = m.get(*_bathymetry);
  val.resize(m.nEvaluationPoint(), 1);
  for (int i=0; i<m.nEvaluationPoint(); i++) {
    double H=hydroSol(i,0) + bath(i,0);
    val (i,0)  =  source(i,0) * H;
lambrechts's avatar
lambrechts committed
71
  }
72
73
74
75
76
77
78
}

void dgConservationLawShallowWaterTracer2d::riemannRoeCB(functorCache &m, fullMatrix<double>&val) const {
  const fullMatrix<double> &solL = m.get(m.solution(), 0);
  const fullMatrix<double> &solR = m.get(m.solution(), 1);
  const fullMatrix<double> &hydroSolL = m.get(*_hydroSolution, 0);
  const fullMatrix<double> &hydroSolR = m.get(*_hydroSolution, 1);
79
80
  const fullMatrix<double> &bathL = m.get(*_bathymetry, 0);
  const fullMatrix<double> &bathR = m.get(*_bathymetry, 1);
David Vincent's avatar
David Vincent committed
81
  const fullMatrix<double> &g = m.get(*_g, 0);
82
83
  const fullMatrix<double> &normalsL = m.get(*function::getNormals(), 0);
  const fullMatrix<double> &normalsR = m.get(*function::getNormals(),1);
84
85
86
87
88
89
  const fullMatrix<double> &movingBathFactorL = m.get(*_movingBathFactor, 0);
  const fullMatrix<double> &movingBathFactorR = m.get(*_movingBathFactor, 1);
  const fullMatrix<double> *ip = _ipTerm ? &m.get(*_ipTerm) : NULL;
  const fullMatrix<double> *xyz = _R > 0 ? &m.get(m.coordinates(), 0) : NULL;
  val.resize(m.nEvaluationPoint(), 2);
  for(int i=0; i< m.nEvaluationPoint(); i++) {
lambrechts's avatar
lambrechts committed
90
    double J = 1.0;
91
92
93
94
95
    if(_R > 0) {
      double x = (*xyz)(i,0);
      double y = (*xyz)(i,1);
      J = 4.0*_R*_R/(4.0*_R*_R+x*x+y*y);
    }
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
    double nxL = normalsL(i,0);
    double nyL = normalsL(i,1);
    double nxR,nyR;
    if (m.interfaceGroup()->nConnection() ==1){
        nxR = normalsR(i,0);
        nyR = normalsR(i,1);  
        
    }
    else {
        nxR = -normalsR(i,0);
        nyR = -normalsR(i,1);
    }
    //printf("nl=(%g,%g), nR=(%g,%g)\n",nxL,nyL,nxR,nyR);
    double uL = nxL*hydroSolL(i,1) + nyL*hydroSolL(i,2);
    double vL = -nyL*hydroSolL(i,1) + nxL*hydroSolL(i,2);
    double uR = nxR*hydroSolR(i,1) + nyR*hydroSolR(i,2);   
    double vR = - nyR*hydroSolR(i,1) + nxR*hydroSolR(i,2);
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
    double Feta, u;
    if (_laxFriedrichs) {
      double hL = bathL(i,0);
      double hR = bathR(i,0);
      double etaR = hydroSolR(i,0);
      double etaL = hydroSolL(i,0);
      const double etaM = (etaR + etaL) /2;
      const double hMin = std::min(hL, hR);
      uR = uR * hR/hMin;
      uL = uL * hL/hMin;
      vR = vR * hR/hMin;
      vL = vL * hL/hMin;
      const double unM = (uR + uL) /2;
      const double c = sqrt(g(i,0) * (hMin + etaM) + fabs(unM));
      Feta =  J * (- (hMin + etaM) * unM)  + c * (etaR - etaL) / 2;
      u =unM;
    }
    else{
      double h = (bathL(i,0)+bathR(i,0))/2.;
      double HR = hydroSolR(i,0) + h;
      double HL = hydroSolL(i,0) + h;
      if(HR<0 || HL<0) printf("in sw2d tracer Riemann:  HR = %e HL = %e\n", HR,HL);
      double v,H,Amb;
      roeSolver(uL, uR, vL, vR, HL, HR, u, v, H, Amb, g(i,0), movingBathFactorL(i,0), movingBathFactorR(i,0), _upwindFactorRiemann);
      Feta = -H*u/J/Amb;
    }    
    double cL = solL(i,0), cR = solR(i,0);
    val(i,0) = Feta  * (u > 0 ? cL : cR) ;
141
142
143
144
    val(i,1) = -val(i,0);
    if(ip){
      val(i,0) += (*ip)(i,0);
      val(i,1) += (*ip)(i,1);
lambrechts's avatar
lambrechts committed
145
146
    }
  }
147
}
lambrechts's avatar
lambrechts committed
148
149
150
151
152
153


//**********************************************
// Boundary conditions
//**********************************************

154
155
156
157
void dgConservationLawShallowWaterTracer2d::wallBoundaryCB(functorCache &m, fullMatrix<double> &val)const{
  val.resize(m.nEvaluationPoint(), 1);
  val.setAll(0.);
}
lambrechts's avatar
lambrechts committed
158
159

dgBoundaryCondition *dgConservationLawShallowWaterTracer2d::newBoundaryWall(){
160
  return newBoundaryConditionMember(*this, &dgConservationLawShallowWaterTracer2d::wallBoundaryCB);
lambrechts's avatar
lambrechts committed
161
162
}

163
164
165
166
167
168
169
void dgConservationLawShallowWaterTracer2d::forcedBoundaryCB(functorCache &m, fullMatrix<double>&val, const functor *hydroSolExtF, const functor *tracerExtF, const functor *tracerGradientExtF) const {
  const fullMatrix<double> &normals = m.get(*function::getNormals());
  const fullMatrix<double> &hydroIn = m.get(*_hydroSolution, 0);
  const fullMatrix<double> &hydroExt = m.get(*hydroSolExtF);
  const fullMatrix<double> &tracerIn = m.get(*function::getSolution(), 0);
  const fullMatrix<double> &tracerExt = m.get(*tracerExtF);
  const fullMatrix<double> &bath = m.get(*_bathymetry);
David Vincent's avatar
David Vincent committed
170
  const fullMatrix<double> &g = m.get(*_g);
171
172
173
174
175
176
177
178
179
180
181
182
183
184
  const fullMatrix<double> &movingBathFactor = m.get(*_movingBathFactor);
  val.resize(m.nEvaluationPoint(), 1);
  for (int i=0; i<m.nEvaluationPoint(); i++){
    double nx = normals(i,0);
    double ny = normals(i,1);
    double etaIn = hydroIn(i,0), etaExt = hydroExt(i,0);
    double uIn = nx * hydroIn(i,1) + ny * hydroIn(i,2);
    double uExt = nx * hydroExt(i,1) + ny * hydroExt(i,2);
    double vIn = -ny * hydroIn(i,1) + nx * hydroIn(i,2);
    double vExt = -ny * hydroExt(i,1) + nx * hydroExt(i,2);
    double cIn = tracerIn(i,0), cExt = tracerExt(i,0);
    double HIn= etaIn + bath(i,0) , HExt = etaExt + bath(i,0);
    if(HIn<0 || HExt<0) printf(" HIn = %e HExt =%e\n", HIn, HExt);
    double u,v,H,Amb;
David Vincent's avatar
David Vincent committed
185
    roeSolver(uIn, uExt, vIn, vExt, HIn, HExt, u, v, H, Amb, g(i,0),  movingBathFactor(i,0), movingBathFactor(i,0), _upwindFactorRiemann);
186
    val(i,0) = -H*u * ( uIn > 0 ? cIn : cExt );   
lambrechts's avatar
lambrechts committed
187
  }
188
}
lambrechts's avatar
lambrechts committed
189

190
dgBoundaryCondition *dgConservationLawShallowWaterTracer2d::newForceTracer(const functor *hydroSolExtF, const functor *tracerExtF, const functor *tracerGradientExtF){
191
192
  Msg::Warning("I do not think the force tracer boundary is correct, for example the tracerGradientExt parameter is not taken into account");
  return newBoundaryConditionMember(*this, &dgConservationLawShallowWaterTracer2d::forcedBoundaryCB, hydroSolExtF, tracerExtF, tracerGradientExtF);
lambrechts's avatar
lambrechts committed
193
194
195
196
197
198
199
}

//**********************************************
// Setup
//**********************************************
void dgConservationLawShallowWaterTracer2d::setup () {
  if(_useMovingBathWD){
200
201
    _bathymetry = new movingBath(_originalBathymetry, _hydro, _alphaMovingBathWD);
    _movingBathFactor = new movingBathFactor(_originalBathymetry, _hydro, _alphaMovingBathWD);
lambrechts's avatar
lambrechts committed
202
203
204
205
  }

  _diffusivity[""] = _kappa;
  if(_kappa){
206
    _diffusiveFlux = newFunctorMember(*this, &dgConservationLawShallowWaterTracer2d::diffusiveFluxCB);
lambrechts's avatar
lambrechts committed
207
    _diffusion = _kappa;
208
    _ipTerm = dgNewIpTermIsotropic(1, _diffusiveFlux, _diffusion);
lambrechts's avatar
lambrechts committed
209
210
  }

211
212
213
214
  _volumeTerm0[""] = newFunctorMember(*this, &dgConservationLawShallowWaterTracer2d::psiTermCB);
  _volumeTerm1[""] = newFunctorMember(*this, &dgConservationLawShallowWaterTracer2d::gradPsiTermCB);
  _interfaceTerm0[""] = newFunctorMember(*this, &dgConservationLawShallowWaterTracer2d::riemannRoeCB);
  _massFactor[""] = std::make_pair(newFunctorMember(*this, &dgConservationLawShallowWaterTracer2d::massFactorCB), MASS_FACTOR_TIME_DEPENDENT);
lambrechts's avatar
lambrechts committed
215
216
217
218
219
220
221
}

dgConservationLawShallowWaterTracer2d::dgConservationLawShallowWaterTracer2d() : dgConservationLawFunction(1) {
  fzero = new functionConstant(0.);
  fzerov = new functionConstant(std::vector<double>(2, 0.));
  _source = fzero;
  _hydroSolution=NULL;
222
  _kappa = NULL;
lambrechts's avatar
lambrechts committed
223
224
225
  _ipTerm = NULL;
  _diffusiveFlux = NULL;
  _R = -1.0;
226
227
  static const functor *gdefault = new functionConstant(9.80616);
  _g = gdefault;
lambrechts's avatar
lambrechts committed
228
229
230
  _alphaMovingBathWD = 0.0;
  _useMovingBathWD = false;
  _bathymetry = NULL;
231
  _hydro = NULL;
lambrechts's avatar
lambrechts committed
232
233
234
  _originalBathymetry = NULL;
  _movingBathFactor = new functionConstant(1.);
  _upwindFactorRiemann = -1;
235
236
  _isLinear = false;
  _isConstantJac = false;
lambrechts's avatar
lambrechts committed
237
238
239
240
241
242
243
244
245
246
247
248
249
250
}

dgConservationLawShallowWaterTracer2d::~dgConservationLawShallowWaterTracer2d () {
  delete fzero;
  delete fzerov;
  delete _movingBathFactor;
  if (_volumeTerm0[""]) delete _volumeTerm0[""];
  if (_volumeTerm1[""]) delete _volumeTerm1[""];
  if (_interfaceTerm0[""]) delete _interfaceTerm0[""];
  if (_massFactor.count("") != 0)     delete _massFactor[""].first;
  if (_useMovingBathWD) {
    delete _bathymetry;
  }
}