dgConservationLawShallowWater2d.cpp 39.6 KB
Newer Older
lambrechts's avatar
lambrechts committed
1
2
3
4
5
6
7
8
#include "dgConservationLawShallowWater2d.h"
#include "dgGroupOfElements.h"
#include "dgDofContainer.h"
#include "functionGeneric.h"
#include "float.h"
#include "dgFunctionIntegratorInterface.h"
#include "dgMeshJacobian.h"
#include "slimMovingBathWettingDrying.h"
bertrand's avatar
bertrand committed
9
#include "functorMember.h"
10
11
#include "functor.h"

lambrechts's avatar
lambrechts committed
12
13
14
15
16
17
18
19
/*==============================================================================
 * DG-terms
 *============================================================================*/

class dgConservationLawShallowWater2d::clipToPhysics : public function {
  fullMatrix<double> sol, bathymetry;
  double _hThreshold;
 public:
20
  clipToPhysics(const functor *bathymetryF, double hMin) : function(3) {
lambrechts's avatar
lambrechts committed
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
    setArgument (sol, function::getSolution());
    setArgument (bathymetry, bathymetryF);
    _hThreshold = hMin;
  };
  void call (dataCacheMap *m, fullMatrix<double> &val) {
    const size_t nQP = val.size1();
    for (size_t k=0; k<nQP; k++){
      double u = sol(k,1);
      double v = sol(k,2);
      val(k,0) = sol(k,0);
      val(k,1) = u;
      val(k,2) = v;
      //val(k,1) = u * fu;
      //val(k,2) = v * fu;
      if (!((sol(k,0)+bathymetry(k,0))>0.05)) {
        Msg::Info("Solution Clipped, H_%d=%f",k, sol(k,0)+bathymetry(k,0));
        val(k,0) = -bathymetry(k,0)+0.05;
      }
    }
  }
  /*void call (dataCacheMap *m, fullMatrix<double> &val) {
    const size_t nQP = val.size1();
    double meanH=0;
    double minH = DBL_MAX;
    bool elementDry = false;
    int highestBathK = 0;
    int highestEtaK = 0;
    int highestHK = 0;
    for (size_t k = 0 ; k < nQP; k++ ){
      val(k,0) = sol(k,0) + bathymetry(k,0);
      if (bathymetry(k,0) < bathymetry(highestBathK,0))
        highestBathK = k;
      if (sol(k,0) > sol(highestEtaK,0))
        highestEtaK = k;
      if (sol(k,0) + bathymetry(k,0) > sol(highestHK,0)  + bathymetry(highestHK,0))
        highestHK = k;
      meanH += val (k,0);
      minH = std::min (minH, val (k,0));
    }
    meanH /= nQP;
    bool isDambreak = sol(highestHK,0) + bathymetry(highestBathK,0) - _hThreshold > 0;
    // from bunya paper
    elementDry = dryMap[m->getElement()] ?
        meanH < _hThreshold || !isDambreak :
        meanH < _hThreshold;
//     if(dryMap[cacheMap->getElement()] == 1 && elementDry == 0)
//       printf("element %p DRY->WET\n",cacheMap->getElement());
//     if(dryMap[cacheMap->getElement()] == 0 && elementDry == 1)
//       printf("element %p      WET->DRY\n",cacheMap->getElement());
    dryMap[m->getElement()] = elementDry;
    //     if(meanH<_hThreshold) Msg::Warning("meanH = %f\n",meanH);
    double _hLimiter = 0.5*_hThreshold;
    if ( minH < _hLimiter) {
      double alpha = std::max((_hLimiter - meanH) / (minH - meanH), 0.);
      for (size_t k = 0 ; k < nQP; k++) {
        double h = bathymetry(k,0);
        val (k,0) = meanH + (val(k,0) - meanH)*alpha - h;
        val (k,1) = isDambreak ? sol (k,1) : 0; // U needed for flooding
        val (k,2) = isDambreak ? sol (k,2) : 0;
      }
    } else {
      for (size_t k = 0 ; k < nQP; k++) {
        val (k,0) = sol(k,0);
        val(k,1) = sol(k,1);
        val(k,2) = sol(k,2);
      }
    }
  }*/
};

class dgConservationLawShallowWater2d::maxConvectiveSpeed: public function {
David Vincent's avatar
David Vincent committed
92
  fullMatrix<double> bath, sol, _xyz, _g;
lambrechts's avatar
lambrechts committed
93
94
95
  bool _linear, _filterLinear;
  double _R;
 public:
David Vincent's avatar
David Vincent committed
96
  maxConvectiveSpeed (const functor *bathymetry, bool linear, double R, bool filterLinear, const functor *g):function(1){
lambrechts's avatar
lambrechts committed
97
98
99
100
101
    setArgument(bath,bathymetry);
    setArgument(sol,function::getSolution());
    _linear = linear;
    _filterLinear = filterLinear;
    _R = R;
David Vincent's avatar
David Vincent committed
102
    setArgument(_g, g);
lambrechts's avatar
lambrechts committed
103
104
105
106
107
108
109
    if(_R>0) {
      setArgument(_xyz, function::getCoordinates());
    }
  }
  void call(dataCacheMap *m, fullMatrix<double> &val) {
    if (_linear){
      if (_filterLinear)
110
        val.setAll(0);
111
      else{
112
        for (int i=0; i<val.size1(); i++)
David Vincent's avatar
David Vincent committed
113
          val(i,0) = sqrt(bath(i,0)*_g(i,0));
114
          }
lambrechts's avatar
lambrechts committed
115
    }else
116
117
      if (!_filterLinear){
        for (int i=0; i<val.size1(); i++) {
David Vincent's avatar
David Vincent committed
118
          val(i,0) = sqrt((bath(i,0)+sol(i,0))*_g(i,0))+hypot(sol(i,1),sol(i,2));
119
120
121
        }
      }
      else{
122
        for (int i=0; i<val.size1(); i++)
David Vincent's avatar
David Vincent committed
123
          val(i,0) = sqrt(sol(i,0)*_g(i,0))+hypot(sol(i,1),sol(i,2));  }
124
125
      
      
lambrechts's avatar
lambrechts committed
126
127
    /*if (_R>0)
      for (int i=0; i<val.size1(); i++) {
128
129
130
131
      double x = _xyz(i,0);
      double y = _xyz(i,1);
      double J = 4.0*_R*_R/(4.0*_R*_R+x*x+y*y);
      val(i,0) /= J;
lambrechts's avatar
lambrechts committed
132
133
134
135
      }*/
  }
};

Valentin Vallaeys's avatar
Valentin Vallaeys committed
136
137
138
139
140
141
142
143
144
145
146
147
void dgConservationLawShallowWater2d::gradPsiTerm(functorCache &cm, fullMatrix<double> &val) const {
    const fullMatrix<double> &bath = cm.get(*_bathymetry);
    const fullMatrix<double> &sol = cm.get(*function::getSolution());
    const fullMatrix<double> *diff = _diffusiveFlux ? &cm.get(*_diffusiveFlux) : NULL;
    const fullMatrix<double> &g = cm.get(*_g);
    const fullMatrix<double> *xyz = _R ? &cm.get(*function::getCoordinates()) : NULL;
    const fullMatrix<double> &movingBathFactor = cm.get(*_movingBathFactor);
    const fullMatrix<double> &atmPress = cm.get(*_atmPress);
    const fullMatrix<double> &rhoSurf = cm.get(*_rhoSurf);
    const fullMatrix<double> *etaOld = _linearSolverFrom3D ? &cm.get(*_etaOld) : NULL;
    size_t nQP = cm.nEvaluationPoint();
    val.resize(nQP, 9, false);
lambrechts's avatar
lambrechts committed
148
149
150
    double J = 1.0;
    for(size_t i=0; i< nQP; i++) {
      if (_R > 0){
Valentin Vallaeys's avatar
Valentin Vallaeys committed
151
152
        double x = (*xyz)(i,0);
        double y = (*xyz)(i,1);
lambrechts's avatar
lambrechts committed
153
        J = 4.0*_R*_R/(4.0*_R*_R+x*x+y*y);
Valentin Vallaeys's avatar
Valentin Vallaeys committed
154
        J = 1.0;
lambrechts's avatar
lambrechts committed
155
156
157
158
159
160
161
162
      }
      double H = bath(i,0);
      double eta = sol(i,0);
      double u = sol(i,1);
      double v = sol(i,2);
      double A = movingBathFactor(i,0);//equals 1 if no WD
      if (!_linear)
        H+=eta;
Valentin Vallaeys's avatar
Valentin Vallaeys committed
163
164
      else if (_linearSolverFrom3D)
        H+=(*etaOld)(i,0);
lambrechts's avatar
lambrechts committed
165
166
      // flux_x
      val(i,0) = H*u * J/A;
Valentin Vallaeys's avatar
Valentin Vallaeys committed
167
      val(i,1) = ( g(i,0)*(1+rhoSurf(i,0)/_density)*eta + atmPress(i,0)/_density )*J;
lambrechts's avatar
lambrechts committed
168
169
170
171
      val(i,2) = 0;
      // flux_y
      val(i,3) = H*v * J/A;
      val(i,4) = 0;
Valentin Vallaeys's avatar
Valentin Vallaeys committed
172
173
      val(i,5) = ( g(i,0)*(1+rhoSurf(i,0)/_density) *eta + atmPress(i,0)/_density )*J;
      if (!_linear && !_from3D) {
lambrechts's avatar
lambrechts committed
174
175
176
177
178
        val(i,1) += u*u*J;
        val(i,2) += u*v*J;
        val(i,4) += v*u*J;
        val(i,5) += v*v*J;
      }
Valentin Vallaeys's avatar
Valentin Vallaeys committed
179
180
181
182
183
      if (diff) {
        val(i,1) += (*diff)(i,1);
        val(i,2) += (*diff)(i,2);
        val(i,4) += (*diff)(i,4);
        val(i,5) += (*diff)(i,5);
lambrechts's avatar
lambrechts committed
184
185
186
187
188
189
      }
      // flux_z
      val(i,6) = 0;
      val(i,7) = 0;
      val(i,8) = 0;
    }
Valentin Vallaeys's avatar
Valentin Vallaeys committed
190
}
lambrechts's avatar
lambrechts committed
191

Valentin Vallaeys's avatar
Valentin Vallaeys committed
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
void dgConservationLawShallowWater2d::source(functorCache &cm, fullMatrix<double> &val) const {
    const fullMatrix<double> &g = cm.get(*_g);
    const fullMatrix<double> &sol = cm.get(*function::getSolution());
    const fullMatrix<double> &coriolisFactor = cm.get(*_coriolisFactor);
    const fullMatrix<double> &linearDissipation = cm.get(*_linearDissipation);
    const fullMatrix<double> &quadraticDissipation = cm.get(*_quadraticDissipation);
    const fullMatrix<double> *nudgingCoeff = _nudgingCoeff ? &cm.get(*_nudgingCoeff) : NULL;
    const fullMatrix<double> *nudgingVel = _nudgingVel ? &cm.get(*_nudgingVel) : NULL;
    const fullMatrix<double> &source = cm.get(*_source);
    const fullMatrix<double> &windStress = cm.get(*_windStress);
    const fullMatrix<double> &bathymetry = cm.get(*_bathymetry);
    const fullMatrix<double> *solGradient = (_nu || !_linear) ? &cm.get(*function::getSolutionGradient()) : NULL;
    const fullMatrix<double> *nu = _nu ? &cm.get(*_nu) : NULL;
    const fullMatrix<double> *bathymetryGradient = _nu ? &cm.get(*_bathymetryGradient) : NULL;
    const fullMatrix<double> *xyz = _R ? &cm.get(*function::getCoordinates()) : NULL;
    const fullMatrix<double> *movingBathFactor = _useMovingBathWD ? &cm.get(*_movingBathFactor) : NULL;
    const fullMatrix<double> *movingBathFactorGradient = _useMovingBathWD ? &cm.get(*_movingBathFactorGradient) : NULL;
    const fullMatrix<double> &rhoSurfGrad = cm.get(*_rhoSurfGrad);
    size_t nQP = cm.nEvaluationPoint();
    val.resize(nQP, 3, false);
lambrechts's avatar
lambrechts committed
212
213
214
    double J = 1.0;
    for(size_t i=0; i< nQP; i++) {
      if(_R > 0){
Valentin Vallaeys's avatar
Valentin Vallaeys committed
215
216
        double x = (*xyz)(i,0);
        double y = (*xyz)(i,1);
lambrechts's avatar
lambrechts committed
217
        J = 4.0*_R*_R/(4.0*_R*_R+x*x+y*y);
Valentin Vallaeys's avatar
Valentin Vallaeys committed
218
        J = 1.0;
lambrechts's avatar
lambrechts committed
219
220
221
222
      }
      double eta = sol(i,0);
      double u = sol(i,1);
      double v = sol(i,2);
Valentin Vallaeys's avatar
Valentin Vallaeys committed
223
      double H = eta + bathymetry(i,0);
lambrechts's avatar
lambrechts committed
224
      val (i,0) = 0;//1/_R/_R*(_xyz(i,0)*u + _xyz(i,1)*v);//0.;
Valentin Vallaeys's avatar
Valentin Vallaeys committed
225
226
      val (i,1) =  coriolisFactor(i,0)*v - linearDissipation(i,0) * u + source(i,0) + windStress(i,0)/(H*_density) + g(i,0)*eta*rhoSurfGrad(i,0)/_density;
      val (i,2) = -coriolisFactor(i,0)*u - linearDissipation(i,0) * v + source(i,1) + windStress(i,1)/(H*_density) + g(i,0)*eta*rhoSurfGrad(i,1)/_density;
lambrechts's avatar
lambrechts committed
227
      if (_R>0) {
Valentin Vallaeys's avatar
Valentin Vallaeys committed
228
229
        val (i,1) -= g(i,0)*(eta)*(*xyz)(i,0)/2.0/_R/_R;
        val (i,2) -= g(i,0)*(eta)*(*xyz)(i,1)/2.0/_R/_R;
lambrechts's avatar
lambrechts committed
230
      }
Valentin Vallaeys's avatar
Valentin Vallaeys committed
231
232
233
234
      if(_nu){
        double H = bathymetry(i,0);
        double dHdx = (*bathymetryGradient)(i, 0);
        double dHdy = (*bathymetryGradient)(i, 1);
lambrechts's avatar
lambrechts committed
235
236
        if (!_linear) {
          H += sol(i,0);
Valentin Vallaeys's avatar
Valentin Vallaeys committed
237
238
          dHdx += (*solGradient)(i,0);
          dHdy += (*solGradient)(i,1);
lambrechts's avatar
lambrechts committed
239
        }
Valentin Vallaeys's avatar
Valentin Vallaeys committed
240
241
        val(i,1) += (*nu)(i, 0)/H*(dHdx*(*solGradient)(i, 3)+dHdy*(*solGradient)(i, 4));
        val(i,2) += (*nu)(i, 0)/H*(dHdx*(*solGradient)(i, 6)+dHdy*(*solGradient)(i, 7));
lambrechts's avatar
lambrechts committed
242
243
      }
      if(!_linear){
Valentin Vallaeys's avatar
Valentin Vallaeys committed
244
        double div = (*solGradient)(i,3) + (*solGradient)(i,7);
lambrechts's avatar
lambrechts committed
245
        double bottomFriction = hypot(u,v)*quadraticDissipation(i,0);
Valentin Vallaeys's avatar
Valentin Vallaeys committed
246
247
        val(i,1) += u*(- bottomFriction + div/J);
        val(i,2) += v*(- bottomFriction + div/J);
lambrechts's avatar
lambrechts committed
248
        if (_R > 0) {
Valentin Vallaeys's avatar
Valentin Vallaeys committed
249
250
          val(i,1)+= (-(*xyz)(i,0)*(u*u+v*v))/(2.0*_R*_R);
          val(i,2)+= (-(*xyz)(i,1)*(u*u+v*v))/(2.0*_R*_R);   
lambrechts's avatar
lambrechts committed
251
252
253
        }
      }
      if(_useMovingBathWD){
Valentin Vallaeys's avatar
Valentin Vallaeys committed
254
255
256
257
        double H = bathymetry(i,0) + eta;
        double A = (*movingBathFactor)(i,0);
        double dAdx = (*movingBathFactorGradient)(i,0);
        double dAdy = (*movingBathFactorGradient)(i,1);
lambrechts's avatar
lambrechts committed
258
259
260
261
262
        val(i,0) += - H / ( A*A ) * ( dAdx * u + dAdy * v);
      }
      val(i,0) *= J*J;
      val(i,1) *= J*J;
      val(i,2) *= J*J;
Valentin Vallaeys's avatar
Valentin Vallaeys committed
263
264
265
      if (_nudgingCoeff) {
        double uref = (*nudgingVel)(i,1);
        double vref = (*nudgingVel)(i,2);
Valentin Vallaeys's avatar
Valentin Vallaeys committed
266
        if (_nudgingVelIsTransport){
Valentin Vallaeys's avatar
Valentin Vallaeys committed
267
268
          uref = uref/bathymetry(i,0);
          vref = vref/bathymetry(i,0);
Valentin Vallaeys's avatar
Valentin Vallaeys committed
269
        }
270
        //double sigma[2] = {_absLayerCoef(i,0), _absLayerCoef(i,1)};
Valentin Vallaeys's avatar
Valentin Vallaeys committed
271
        double sigma = (*nudgingCoeff)(i,0);
Valentin Vallaeys's avatar
Valentin Vallaeys committed
272
273
        val(i,1) -= sigma * (u-uref);
        val(i,2) -= sigma * (v-vref);
274
275
276
277
278
279
        // Absorbing Layer of Lavelle & co (Pretty Good Sponge)
        /*double sigma[2] = {_absLayerCoef(i,0), _absLayerCoef(i,1)};
        val(i,0) -= (sigma[0]+sigma[1]) * (eta-_absLayerExtForc(i,0));
        val(i,1) -= sigma[0] * (u-_absLayerExtForc(i,1));
        val(i,2) -= sigma[1] * (v-_absLayerExtForc(i,2));*/
        // Absorbing Layer of Martinsen & co (Sponge Layer - FRS)
Valentin Vallaeys's avatar
Valentin Vallaeys committed
280
281
282
283
        //double sigma = _nudgingCoeff(i,0);
        //val(i,0) -= sigma * (eta-_nudgingVel(i,0));
        //val(i,1) -= sigma * (u-uref);
        //val(i,2) -= sigma * (v-vref);
amodave's avatar
amodave committed
284
      }
lambrechts's avatar
lambrechts committed
285
    }
Valentin Vallaeys's avatar
Valentin Vallaeys committed
286
}
lambrechts's avatar
lambrechts committed
287
288
289
290
291

class massFactor: public function {
  fullMatrix<double> _xyz;
  double _R;
 public :
292
  massFactor(double R, const functor *xyz) : function(3) {
lambrechts's avatar
lambrechts committed
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
    _R = R;
    setArgument(_xyz, xyz ? xyz : function::getCoordinates());
  }
  void call(dataCacheMap *m, fullMatrix<double> &val) {
    int nQP=val.size1();
    for (int i=0; i<nQP; i++) {
      double x = _xyz(i,0);
      double y = _xyz(i,1);
      double J = 4.0*_R*_R/(4.0*_R*_R+x*x+y*y);
      val(i,0) = J*J;
      val(i,1) = J*J;
      val(i,2) = J*J;
    }
  }
};

Valentin Vallaeys's avatar
Valentin Vallaeys committed
309
310
311
312
313
void dgConservationLawShallowWater2d::diffusiveFlux(functorCache &cm,fullMatrix<double> &val) const {
    const fullMatrix<double> &solGrad = cm.get(*function::getSolutionGradient());
    const fullMatrix<double> &nu = cm.get(*_nu);
    size_t nQP = cm.nEvaluationPoint();
    val.resize(nQP, 6, false);
lambrechts's avatar
lambrechts committed
314
315
316
317
318
319
320
321
    for (size_t i=0; i<nQP; i++) {
      val(i,0) = 0;
      val(i,1) = -nu(i,0)*solGrad(i,3);
      val(i,2) = -nu(i,0)*solGrad(i,6);
      val(i,3) = 0;
      val(i,4) = -nu(i,0)*solGrad(i,4);
      val(i,5) = -nu(i,0)*solGrad(i,7);
    }
Valentin Vallaeys's avatar
Valentin Vallaeys committed
322
}
lambrechts's avatar
lambrechts committed
323

Valentin Vallaeys's avatar
Valentin Vallaeys committed
324
325
326
327
void dgConservationLawShallowWater2d::diffusivity(functorCache &cm, fullMatrix<double> &val) const {
    const fullMatrix<double> &nu = cm.get(*_nu);
    size_t nQP = cm.nEvaluationPoint();
    val.resize(nQP, 3, false);
lambrechts's avatar
lambrechts committed
328
329
330
331
332
    for(size_t i=0; i<nQP; i++) {
      val(i,0) = 0;
      val(i,1) = nu(i,0);
      val(i,2) = nu(i,0);
    }
Valentin Vallaeys's avatar
Valentin Vallaeys committed
333
}
lambrechts's avatar
lambrechts committed
334

bertrand's avatar
bertrand committed
335
336
void dgConservationLawShallowWater2d::riemann(functorCache &cache,fullMatrix<double> &val) const {
  const size_t nQP = cache.nEvaluationPoint();
337
  val.resize(nQP,6);
lambrechts's avatar
lambrechts committed
338
339
  const fullMatrix<double> &solL = cache.get(cache.solution(), 0);
  const fullMatrix<double> &solR = cache.get(cache.solution(), 1);
340
341
  const fullMatrix<double> &bathL = cache.get(*_bathymetry,0);
  const fullMatrix<double> &bathR = cache.get(*_bathymetry,1);
David Vincent's avatar
David Vincent committed
342
  const fullMatrix<double> &g = cache.get(*_g,0);
343
344
  const fullMatrix<double> &normalsL = cache.get(*function::getNormals(),0);
  const fullMatrix<double> &normalsR = cache.get(*function::getNormals(),1);
bertrand's avatar
bertrand committed
345
346
  const fullMatrix<double> &movingBathFactorL = cache.get(*_movingBathFactor,0);
  const fullMatrix<double> &movingBathFactorR = cache.get(*_movingBathFactor,1);
Valentin Vallaeys's avatar
Valentin Vallaeys committed
347
348
349
350
351
352
  const fullMatrix<double> &atmPressL = cache.get(*_atmPress,0);
  const fullMatrix<double> &atmPressR = cache.get(*_atmPress,1);
  const fullMatrix<double> &rhoSurfL = cache.get(*_rhoSurf,0);
  const fullMatrix<double> &rhoSurfR = cache.get(*_rhoSurf,1);
  const fullMatrix<double> *etaOldL = (_linearSolverFrom3D) ? &cache.get(*_etaOld,0) : NULL;
  const fullMatrix<double> *etaOldR = (_linearSolverFrom3D) ? &cache.get(*_etaOld,1) : NULL;
bertrand's avatar
bertrand committed
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
  bool _isDiffusive = _ipTerm;
  double _upwindFactor = _upwindFactorRiemann;
  const fullMatrix<double> *ip = _isDiffusive ? &cache.get(*_ipTerm) : NULL;
  const fullMatrix<double> *xyz = _R>0 ? &cache.get(*_xyz) : NULL;

  int rev = 1;
  double J = 1.0;
  const dgGroupOfFaces &iGroup = *cache.interfaceGroup();
  int connectionR = iGroup.nConnection() == 2 ? 1 : 0;
  const dgGroupOfElements *groupR = &iGroup.elementGroup(connectionR);
  const dgGroupOfElements *groupL = &iGroup.elementGroup(0);
  for(size_t i=0; i< nQP; i++) {
    rev = 1;
    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);
      J = 1.0;
      if(groupR->getPhysicalTag() != groupL->getPhysicalTag())
        rev = -1;
lambrechts's avatar
lambrechts committed
374
    }
Valentin Vallaeys's avatar
Valentin Vallaeys committed
375
376
    double pa = 0.5*(atmPressL(i,0) + atmPressR(i,0));  
    double rhoSurf = 0.5*(rhoSurfL(i,0) + rhoSurfR(i,0));
bertrand's avatar
bertrand committed
377
    double Fun, Fut, Feta;
378
379
380
    double nxL = normalsL(i,0);
    double nyL = normalsL(i,1);
    double nxR,nyR;
381
    // the basis for nR is different from that of nL so we must multiply nR by -1 in order to express nR in the same basis (exept for the boundary nodes)
382
383
384
385
386
387
388
389
390
    if (cache.interfaceGroup()->nConnection() ==1){
        nxR = normalsR(i,0);
        nyR = normalsR(i,1);  
    }
    else {
        nxR = -1.*normalsR(i,0);
        nyR = -1.*normalsR(i,1);
    }    

lambrechts's avatar
lambrechts committed
391
/*      if(rev<0)
bertrand's avatar
bertrand committed
392
393
394
    {
      double RR = sqrt(_xyz.get(i,0)*_xyz.get(i,0) + _xyz.get(i,1)*_xyz.get(i,1));
      printf("nx = %.16g nx2 = %.16g, ny = %.16g ny2 = %.16g, 2R = %.16g , RR = %.16g\n ", nx,_xyz.get(i,0)/RR,ny,_xyz.get(i,1)/RR, 2*_R, RR);
lambrechts's avatar
lambrechts committed
395
}*/
396
397
398
    double uL = nxL*solL(i,1) + nyL*solL(i,2);
    double uR = nxR*solR(i,1) + nyR*solR(i,2);
    double vL = -nyL*solL(i,1) + nxL*solL(i,2);
399
    double vR = -nyR*solR(i,1) + nxR*solR(i,2); 
400
401
    double hL = bathL(i,0);
    double hR = bathR(i,0);
bertrand's avatar
bertrand committed
402
    uR *= rev;
403
    if (_laxFriedrichs) {
Valentin Vallaeys's avatar
Valentin Vallaeys committed
404
405
406
407
408
409
410
411
412
      double etaR = solR(i,0);
      double etaL = solL(i,0);
      const double etaM = (etaR + etaL) /2;
      const double hMin = std::min(hL, hR);
      double H = hMin;
      uR = uR * hR/hMin;
      uL = uL * hL/hMin;
      const double unM = (uR + uL) /2;
      double c, utM = 0.;
413
      if (_linear) {
Valentin Vallaeys's avatar
Valentin Vallaeys committed
414
415
416
        c = sqrt(g(i,0) * H);
        if (_linearSolverFrom3D)
          H += ((*etaOldL)(i,0) + (*etaOldR)(i,0))/2;
417
418
      }
      else {
Valentin Vallaeys's avatar
Valentin Vallaeys committed
419
        H += etaM;
420
421
        vR = vR * hR/hMin;
        vL = vL * hL/hMin;
Valentin Vallaeys's avatar
Valentin Vallaeys committed
422
423
424
425
426
427
428
429
430
        utM = (vR + vL) /2;
        c = sqrt(g(i,0) * H + fabs(unM));
      }
      Fun = J * (- g(i,0) * (1+rhoSurf/_density) * etaM - pa/_density ) + c * (uR - uL) / 2;
      Feta = J * (- H * unM)  + c * (etaR - etaL) / 2;
      Fut = 0.;
      if (!_linear && !_from3D){
        Fun += J * (- unM * unM);
        Fut += J * (- unM * utM) + c * (vR - vL) / 2;
431
432
433
      }
    }
    else {
Valentin Vallaeys's avatar
Valentin Vallaeys committed
434
435
      if (_from3D)
        Msg::Fatal("2D from 3D only implemented in LF");
436
437
438
439
      if (_linear) {
        //linear equations 
        double etaR = solR(i,0);
        double etaL = solL(i,0);
440
        const double hMin = std::min(hL, hR);
David Vincent's avatar
David Vincent committed
441
        double sq_g_h = sqrt(g(i,0)/hMin);
442
443
        double etaRiemann = (etaL+etaR + (uL-uR)/sq_g_h)/2;
        double unRiemann = (uL+uR + (etaL-etaR)*sq_g_h)/2;
Valentin Vallaeys's avatar
Valentin Vallaeys committed
444
        Fun = -(g(i,0)*(1+rhoSurf/_density)*etaRiemann + pa/_density)*J;
445
        Fut = 0;
446
        Feta = -(hMin)*unRiemann*J;
447
      }else{
448
        const double hMin = std::min(hL, hR);
449
450
        double HR = solR(i,0) + hMin, HL = solL(i,0) + hMin;
        if ( fabs(hL-hR) > 1e-3 )
451
          Msg::Error("Discontinuous Bathymetry is not implemented for non-linear equation with Roe Solver. Please use the linear equations and/or the LaxFriedrichs Flux.(left %g/%g right)", hL, hR);
452
453
        if(HR<0 || HL<0) printf(" HR = %e HL =%e\n", HR,HL);
        double u,v,H,Amb;
David Vincent's avatar
David Vincent committed
454
        roeSolver(uL, uR, vL, vR, HL, HR, u, v, H, Amb, g(i,0), movingBathFactorL(i,0), movingBathFactorR(i,0), _upwindFactor);
455
        double eta  = H-hMin;
Valentin Vallaeys's avatar
Valentin Vallaeys committed
456
        Fun = -(g(i,0)*(1+rhoSurf/_density)*eta + pa/_density)*J - u*u*J;
457
458
459
        Fut = -u*v*J;
        Feta = -H*u*J/Amb;
      }
bertrand's avatar
bertrand committed
460
461
    }
    val(i,0) = Feta;
462
463
    val(i,1) = Fun * nxL - Fut * nyL;
    val(i,2) = Fun * nyL + Fut * nxL;
464
    val(i,3) = -val(i,0);
465
466
    val(i,4) = -(Fun * nxR - Fut * nyR);
    val(i,5) = -(Fun * nyR + Fut * nxR);
bertrand's avatar
bertrand committed
467
468
469
    if (_isDiffusive) {
      val(i,1)+=(*ip)(i,1);
      val(i,2)+=(*ip)(i,2);
470
471
      val(i,4)+=(*ip)(i,4)*rev;
      val(i,5)+=(*ip)(i,5)*rev;
lambrechts's avatar
lambrechts committed
472
473
    }
  }
bertrand's avatar
bertrand committed
474
}
lambrechts's avatar
lambrechts committed
475
476
477

void roeSolver(double uL, double uR, double vL, double vR, double HL, double HR,
               double &uStar, double &vStar, double &HStar, double &AStar,
478
479
               double g, double mbFactL, double mbFactR , double upwindFact){
  double _g = g;
lambrechts's avatar
lambrechts committed
480
481
482
483
484
485
486
  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 sqHL = sqrt(HL), sqHR = sqrt(HR);
  double u_roe = (sqHL*uL + sqHR*uR) / (sqHL + sqHR);
  double v_roe = (sqHL*vL + sqHR*vR) / (sqHL + sqHR);
487
  double c_roe = sqrt(_g*HM);
lambrechts's avatar
lambrechts committed
488
489
490
491
492
493
494
495
496
497
498
499
500
501
  double Hu = HuM + (c_roe - u_roe*u_roe/c_roe) *HJ + u_roe/c_roe *HuJ;
  double Hv;
  if (upwindFact>0) {
    double upwindFactor=atan(u_roe/100*upwindFact)/M_PI+0.5;
    Hv = -v_roe*u_roe/c_roe*HJ + v_roe/c_roe*HuJ + upwindFactor *( -v_roe*HJ+HvL) + (1-upwindFactor)*(v_roe*HJ+HvR);
  }
  else
    Hv= -v_roe*u_roe/c_roe*HJ + v_roe/c_roe*HuJ + (u_roe>0 ? -v_roe*HJ+HvL : v_roe*HJ+HvR);
  HStar = HM + (HuJ - u_roe *HJ) / c_roe;
  uStar = Hu / HStar;
  vStar = Hv / HStar;
  AStar = (sqHL * mbFactL + sqHR * mbFactL)/(sqHL + sqHR); 
}
	
502
void dgConservationLawShallowWater2d::setCoordinatesFunction(functor *xyz){
lambrechts's avatar
lambrechts committed
503
504
505
506
507
508
509
510
511
512
513
514
515
516
  _xyz=xyz;
}


/*==============================================================================
 * Conservation law : constructor - setup - destructor
 *============================================================================*/

dgConservationLawShallowWater2d::dgConservationLawShallowWater2d() : dgConservationLawFunction(3) {
  // eta u v
  _fzero = new functionConstant(0.);
  _fzerov = new functionConstant(std::vector<double>(2,0));
  _bathymetry = NULL;
  _bathymetryGradient = NULL;
517
518
  _hydro = NULL;
  _hydroGradient = NULL;
lambrechts's avatar
lambrechts committed
519
520
521
522
523
524
  _originalBathymetry = NULL;
  _originalBathymetryGradient = NULL;
  _linearDissipation = _fzero;
  _coriolisFactor = _fzero;
  _quadraticDissipation = _fzero;
  _source = _fzerov;
525
  _windStress = _fzerov;
lambrechts's avatar
lambrechts committed
526
527
528
  _linear = false;
  _constantJac = false;
  _R = -1.0;
529
530
  static const functor *gdefault = new functionConstant(9.80616);
  _g = gdefault;
531
  _density = 1025;
lambrechts's avatar
lambrechts committed
532
533
534
535
536
537
538
539
  _nu = NULL;
  _useMovingBathWD = false;
  _movingBathFactor = new functionConstant(1.);
  _movingBathFactorGradient = _fzerov;
  _upwindFactorRiemann = -1;
  _ipTerm = NULL;
  _diffusiveFlux = NULL;
  _diffusion = NULL;
Valentin Vallaeys's avatar
Valentin Vallaeys committed
540
541
542
543
 // _imexMode = IMEX_ALL;
 // _linearIMEX=false;
 // _constantJacIMEX=false;
 // _linearFilter=false;
lambrechts's avatar
lambrechts committed
544
  _laxFriedrichs = false;
545
546
  _nudgingCoeff = NULL;
  _nudgingVel = NULL;
Valentin Vallaeys's avatar
Valentin Vallaeys committed
547
  _nudgingVelIsTransport = false;
Valentin Vallaeys's avatar
Valentin Vallaeys committed
548
549
550
551
552
553
554
  _atmPress = _fzero;
  _rhoSurf = _fzero;
  _rhoSurfGrad = _fzerov;
  _slip = false;
  _from3D = false;
  _linearSolverFrom3D = false;
  _etaOld = NULL;
lambrechts's avatar
lambrechts committed
555
556
}

Valentin Vallaeys's avatar
Valentin Vallaeys committed
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
//void dgConservationLawShallowWater2d::setImexMode(imexMode mode) {
//  if (mode == IMEX_IMPLICIT){
//    _volumeTerm0[""] = _sourceTermLin;
//    _volumeTerm1[""]=_gradPsiTermLin;
//    _interfaceTerm0[""] = _riemannTermLin;
//    _linearIMEX=true;
//    _constantJacIMEX=true;
//  }else{
//    _volumeTerm0[""] =_sourceTerm;
//    _volumeTerm1[""]=_gradPsiTerm;
//    _interfaceTerm0[""] = _riemannTerm;
//    _linearIMEX=false;
//    _constantJacIMEX=false;
//  }
//  _imexMode = mode;
//}
lambrechts's avatar
lambrechts committed
573

Valentin Vallaeys's avatar
Valentin Vallaeys committed
574
575
576
577
578
579
580
//void dgConservationLawShallowWater2d::setLinearFilterMode(bool set) {
//  if (set)
//    _maximumConvectiveSpeed[""] = _maxSpeedFilt;
//  else
//    _maximumConvectiveSpeed[""] = _maxSpeed;
//  _linearFilter=set;
//}
lambrechts's avatar
lambrechts committed
581
582

void dgConservationLawShallowWater2d::setup() {
Valentin Vallaeys's avatar
Valentin Vallaeys committed
583
  if(_nu && _bathymetry && !_bathymetryGradient) Msg::Fatal("The gradient of the bathymetry is missing since we want to use diffusion with a bathymetry");
584
  _diffusivity[""] = _nu;
Valentin Vallaeys's avatar
Valentin Vallaeys committed
585
586
  _diffusiveFlux = _nu ? newFunctorMember(*this, &dgConservationLawShallowWater2d::diffusiveFlux) : NULL;
  _diffusion = _nu ? newFunctorMember(*this, &dgConservationLawShallowWater2d::diffusivity) : NULL;
587
  _ipTerm = _nu ? dgNewIpTermIsotropicOnSphere(3, _diffusiveFlux, _diffusion) : NULL;
lambrechts's avatar
lambrechts committed
588
  //_volumeTerm0[""] = new source(_linearDissipation, _quadraticDissipation, _coriolisFactor, _source, _linear);
Valentin Vallaeys's avatar
Valentin Vallaeys committed
589
590
  _sourceTerm = newFunctorMember(*this,&dgConservationLawShallowWater2d::source);
  _gradPsiTerm = newFunctorMember(*this,&dgConservationLawShallowWater2d::gradPsiTerm);
bertrand's avatar
bertrand committed
591
  _riemannTerm = newFunctorMember(*this,&dgConservationLawShallowWater2d::riemann);
592
  _maxSpeed = new maxConvectiveSpeed(_bathymetry, _linear, _R, false, _g);
Valentin Vallaeys's avatar
Valentin Vallaeys committed
593
594
595
596
597
598
599
600
601
602
  _volumeTerm0[""] =_sourceTerm;
  _volumeTerm1[""]=_gradPsiTerm;
  _interfaceTerm0[""] = _riemannTerm;
  _maximumConvectiveSpeed[""] = _maxSpeed;
//  _sourceTermLin = new source(_linearDissipation, _quadraticDissipation, _coriolisFactor, _source, _windStress, _bathymetry, _bathymetryGradient, _nu, true, _R, _g, _density, _movingBathFactor, _movingBathFactorGradient, _xyz, _nudgingCoeff, _nudgingVel, _nudgingVelIsTransport);
//  _gradPsiTermLin = new gradPsiTerm(_bathymetry, _diffusiveFlux, true,_R, _movingBathFactor, _xyz, _g);
//  _riemannTermLin = newFunctorMember(*this,&dgConservationLawShallowWater2d::riemann);
//  _maxSpeedFilt = new maxConvectiveSpeed(_bathymetry, _linear, _R, true, _g);
//  setImexMode(_imexMode);
//  setLinearFilterMode(_linearFilter);
lambrechts's avatar
lambrechts committed
603
604
605
606
607
608
609
610
  /*if (_R>0)
    _massFactor[""] = std::make_pair(new massFactor(_R, _xyz), true);*/
  //_interfaceTerm0[""] = _nu ? : new riemann(_bathymetry,NULL, _linear);


  _clipToPhysics[""] = new clipToPhysics(_bathymetry, 0.01);
}

611
void dgConservationLawShallowWater2d::setMovingBathWettingDrying(const functor *hydro, const functor *hydroGradient, double alphaMovingBathWD){
612
613
614
615
  _useMovingBathWD = true;
  _alphaMovingBathWD = alphaMovingBathWD;
  _originalBathymetry = _bathymetry;
  _originalBathymetryGradient= _bathymetryGradient;
616
617
618
619
620
621
  _hydro = hydro; 
  _hydroGradient = hydroGradient; 
  _bathymetry = new movingBath(_originalBathymetry, _hydro, _alphaMovingBathWD);
  _bathymetryGradient = new movingBathGradient(_originalBathymetry, _originalBathymetryGradient, _hydro, _hydroGradient, _alphaMovingBathWD);
  _movingBathFactor = new movingBathFactor(_originalBathymetry, _hydro, _alphaMovingBathWD);
  _movingBathFactorGradient = new movingBathFactorGradient(_originalBathymetry,_originalBathymetryGradient, _hydro, _hydroGradient, _alphaMovingBathWD);
622
623
}

lambrechts's avatar
lambrechts committed
624
625
626
627
628
629
630
631
632
633
634
635
636
637
dgConservationLawShallowWater2d::~dgConservationLawShallowWater2d() {
  delete _fzero;
  delete _fzerov;
  if (_volumeTerm0[""]) delete _volumeTerm0[""];
  if (_volumeTerm1[""]) delete _volumeTerm1[""];
  if (_interfaceTerm0[""]) delete _interfaceTerm0[""];
  if (_ipTerm) delete _ipTerm;
  if (_diffusiveFlux) delete _diffusiveFlux;
  if (_diffusion) delete _diffusion;
  if (_maximumConvectiveSpeed[""]) delete _maximumConvectiveSpeed[""];
  if (_clipToPhysics[""]) delete _clipToPhysics[""];
  if (_massFactor.count("") != 0) delete _massFactor[""].first;
  if (_useMovingBathWD) {
    delete _bathymetry;
Valentin Vallaeys's avatar
Valentin Vallaeys committed
638
  if (_bathymetryGradient)  delete _bathymetryGradient;
lambrechts's avatar
lambrechts committed
639
640
641
642
643
644
645
646
647
  }
}

/*==============================================================================
 * Specific boundary conditions
 *============================================================================*/

// BC : Wall

648
/*class dgConservationLawShallowWater2d::boundaryWall : public dgBoundaryCondition {
lambrechts's avatar
lambrechts committed
649
650
651
652
653
  class term : public function {
    fullMatrix<double> normals, sol, _xyz, bath;
    bool _linear;
    double _R;
   public:
654
    term(bool linear, double R, const functor *bathF) : function(3) {
lambrechts's avatar
lambrechts committed
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
      setArgument(sol,function::getSolution());
      setArgument(_xyz,function::getCoordinates());
      setArgument(normals,function::getNormals());
      setArgument(bath,bathF);
      _linear = linear;
      _R = R;
    }
    void call(dataCacheMap *m, fullMatrix<double> &val) {
      size_t nQP = sol.size1();
      for (size_t i=0; i<nQP; i++) {
        double nx = normals(i,0);
        double ny = normals(i,1);
        double eta = sol(i,0);
        double J = 1.0;
        if (_R>0) {
          double x = _xyz(i,0);
          double y = _xyz(i,1);
          J = 4*_R*_R/(4*_R*_R+x*x+y*y);
	  J = 1.0;
        }
        double u = sol(i, 1);
        double v = sol(i, 2);
        double un = u * nx + v * ny;
        // val(i, 0) = 0;
679
680
        // val(i, 1) = -_g*eta*nx/J;
        // val(i, 2) = -_g*eta*ny/J;
lambrechts's avatar
lambrechts committed
681
682
683
684
685
686
        // if(!_linear && un > 0){
        //   val(i, 1)-=un*u/J;
        //   val(i, 2)-=un*v/J;
        // }
        if ( _linear) {
          val(i, 0) = 0;//bath(i, 0) * un * J;
687
688
          val(i, 1) = -_g*eta*nx;
          val(i, 2) = -_g*eta*ny;
lambrechts's avatar
lambrechts committed
689
690
691
692
693
694
        } else {
          double uL = un;
          double vL = u * (-ny) + v * nx;
          double h = bath(i,0);
          double HL = eta + h;
          double HStar,uStar,vStar,Amb;
695
696
          //roeSolver(uL, 0, vL, 0, HL, HL, uStar, vStar, HStar, Amb, _g);
          roeSolver(uL, -uL, vL, vL, HL, HL, uStar, vStar, HStar, Amb, _g);
lambrechts's avatar
lambrechts committed
697
          val(i,0) = 0;//(-HStar * uStar + (HL * un)) * J;
698
          double FunL = - _g * ((HStar -h)) * J - (uStar*uStar) * J;
lambrechts's avatar
lambrechts committed
699
700
701
          double FutL = - (uStar*vStar) * J;
          val(i,1) = FunL * nx - FutL * ny;
          val(i,2) = FunL * ny + FutL * nx;
702
703
704
705
706
	
	 // if (uStar>0) {
   //         val(i, 1) = (uStar * uStar - uL*uL) / J;
   //         val(i, 2) = (uStar * vStar - uL*vL) / J;
	 //
lambrechts's avatar
lambrechts committed
707
708
709
710
711
712
713
714
715
	}
      }
    }
  };
  term _term;
public:
  boundaryWall(dgConservationLawShallowWater2d *claw) : _term(claw->isLinear(), claw->getRadius(), claw->getBathymetry()) {
    _term0 = &_term;
  }
716
};*/
Valentin Vallaeys's avatar
Valentin Vallaeys committed
717
718
719
720
721
722
void dgConservationLawShallowWater2d::dgConservationLawShallowWater2dWallExtValue(functorCache &cm, fullMatrix<double> &val) const {
    const fullMatrix<double> &normals = cm.get(*function::getNormals(), 0);
    const fullMatrix<double> &solIn = cm.get(*function::getSolution(), 0);
    size_t nQP = cm.nEvaluationPoint();
    val.resize(nQP, 3, false);
    //cm.setDefaultToSide0(true);
723
724
725
    for (size_t i=0; i< nQP; i++) {
      const double nx = normals (i,0);
      const double ny = normals (i,1);
726
727
      val(i, 0) = solIn(i, 0);
      double un = solIn (i,1) * nx + solIn (i,2) * ny;
728
729
730
731
732
733
734
735
      if (_slip) {
        val(i,1) = solIn (i,1) - 2 * un * nx; // with utan different from 0
        val(i,2) = solIn (i,2) - 2 * un * ny;
      }
      else {
        val (i,1) = -un * nx; // with utan = 0
        val (i,2) = -un * ny;
      }
736
    }
Valentin Vallaeys's avatar
Valentin Vallaeys committed
737
}
738

Valentin Vallaeys's avatar
Valentin Vallaeys committed
739
740
741
742
743
744
void dgConservationLawShallowWater2d::dgConservationLawShallowWater2dWallExtValueGrad(functorCache &cm, fullMatrix<double> &val) const {
    const fullMatrix<double> &normals = cm.get(*function::getNormals(), 0);
    const fullMatrix<double> &solGradIn = cm.get(*function::getSolutionGradient(), 0);
    size_t nQP = cm.nEvaluationPoint();
    val.resize(nQP, 9, false);
    //cm.setDefaultToSide0(true);
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
    for (size_t i=0; i< nQP; i++) {
      const double nx = normals (i,0);
      const double ny = normals (i,1);
      val(i, 0) = solGradIn(i, 0);
      val(i, 1) = solGradIn(i, 1);
      val(i, 2) = solGradIn(i, 2);
      val(i, 5) = solGradIn(i, 5);
      val(i, 8) = solGradIn(i, 8);

      double unxIn = solGradIn(i,3)*nx + solGradIn(i,6)*ny;// dudx*nx+dvdx*ny
      double unyIn = solGradIn(i,4)*nx + solGradIn(i,7)*ny;// dudy*nx+dvdy*ny
      double utxIn = -solGradIn(i,3)*ny + solGradIn(i,6)*nx;// dudx*-ny+dvdx*nx
      double utyIn = -solGradIn(i,4)*ny + solGradIn(i,7)*nx;// dudy*-ny+dvdy*nx
      
      double unnIn = unxIn*nx + unyIn*ny; // comp normal de grad un
      //double untIn = -unxIn*ny + unyIn*nx;// comp tan de grad un
      double utnIn = utxIn*nx + utyIn*ny; // comp normal de grad ut
      //double uttIn = -utxIn*ny + utyIn*nx;// comp tan de grad ut
      
      double unnOut = unnIn;
      double utnOut = -utnIn;
      double untOut = 0;
      double uttOut = 0;
      
      double unxOut = unnOut*nx - untOut*ny;
      double unyOut = unnOut*ny + untOut*nx;
      double utxOut = utnOut*nx - uttOut*ny;
      double utyOut = utnOut*ny + uttOut*nx;
      if (_slip) {
        val(i,3) = nx*unxOut - ny*utxOut; 
        val(i,4) = nx*unyOut - ny*utyOut;
        val(i,6) = ny*unxOut + nx*utxOut;
        val(i,7) = ny*unyOut + nx*utyOut;
      }
      else {
        val (i,3) = solGradIn(i,3); 
        val (i,4) = solGradIn(i,4);
        val(i, 6) = solGradIn(i, 6);
        val(i, 7) = solGradIn(i, 7);
      }
    }
Valentin Vallaeys's avatar
Valentin Vallaeys committed
786
}
lambrechts's avatar
lambrechts committed
787

788
dgBoundaryCondition *dgConservationLawShallowWater2d::newBoundaryWall(bool slip){
789
  //return new boundaryWall(this);
Valentin Vallaeys's avatar
Valentin Vallaeys committed
790
791
792
  _slip = slip;
  functor *boundaryWallUV = newFunctorMember(*this, &dgConservationLawShallowWater2d::dgConservationLawShallowWater2dWallExtValue);
  functor *boundaryWallGradUV = newFunctorMember(*this, &dgConservationLawShallowWater2d::dgConservationLawShallowWater2dWallExtValueGrad);
793
794
  std::vector<const functor*> *toReplace = new std::vector<const functor*>();
  std::vector<const functor*> *replaceBy = new std::vector<const functor*>();
795
796
  toReplace->push_back(function::getSolutionGradient());
  replaceBy->push_back(boundaryWallGradUV);
797
798
799
  toReplace->push_back(function::getSolution());
  replaceBy->push_back(boundaryWallUV);
  return newOutsideValueBoundaryGeneric("",*toReplace,*replaceBy);
lambrechts's avatar
lambrechts committed
800
801
802
803
804
805
}

// BC : Forced Discharge

class dgConservationLawShallowWater2d::boundaryForcedDischarge : public dgBoundaryCondition {
  class term : public function {
David Vincent's avatar
David Vincent committed
806
    fullMatrix<double> normals, sol, _xyz, bath, discharge, mbFact, _g;
lambrechts's avatar
lambrechts committed
807
    bool _linear;
David Vincent's avatar
David Vincent committed
808
    double _R, S0;
lambrechts's avatar
lambrechts committed
809
  public:
David Vincent's avatar
David Vincent committed
810
    term(bool linear, double R, const functor *bathF, dgDofContainer *solution, const functor *dischargeF, std::string tag, const functor *movingBathF, const functor *g):function(3){
lambrechts's avatar
lambrechts committed
811
812
813
814
815
816
817
818
      setArgument(sol,function::getSolution());
      setArgument(_xyz,function::getCoordinates());
      setArgument(normals,function::getNormals());
      setArgument(bath,bathF);
      setArgument(discharge,dischargeF);
      setArgument(mbFact,movingBathF);
      _linear = linear;
      _R = R;
David Vincent's avatar
David Vincent committed
819
      setArgument(_g, g);
lambrechts's avatar
lambrechts committed
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
      dgFunctionIntegratorInterface integrator(solution->getGroups(), bathF, solution); 
      fullMatrix<double> Smat;
      integrator.compute(tag, Smat);
      S0 = Smat(0,0);
    }
    void call(dataCacheMap *m, fullMatrix<double> &val) {
      double J = 1;
      double Fun, Fut, Feta;
      for (int i=0; i<sol.size1(); i++) {
        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);
	  J = 1.0;
        }
        double nx = normals(i,0);
        double ny = normals(i,1);
        double uL = nx*sol(i,1) + ny*sol(i,2), uR =  - discharge(i,0) / S0;
        double vL = -ny*sol(i,1) + nx*sol(i,2), vR = vL;
        double h = bath(i,0);
        double etaL = sol(i,0);
        double etaR = etaL;
        if (_linear) {
David Vincent's avatar
David Vincent committed
843
          double sq_g_h = sqrt(_g(i,0)/h);
lambrechts's avatar
lambrechts committed
844
845
          double etaRiemann = (etaL+etaR + (uL-uR)/sq_g_h)/2;
          double unRiemann = (uL+uR + (etaL-etaR)*sq_g_h)/2;
David Vincent's avatar
David Vincent committed
846
          Fun = -_g(i,0)*etaRiemann*J;
lambrechts's avatar
lambrechts committed
847
848
849
850
851
852
          Fut = 0;
          Feta = -(h)*unRiemann*J;
        } else {
          double HR = etaR + h, HL = etaL + h;
          if(HR<0 || HL<0) printf("forced Discharge HR = %e HL =%e\n", HR,HL);
          double u,v,H,Amb;
David Vincent's avatar
David Vincent committed
853
          roeSolver(uL, uR, vL, vR, HL, HR, u, v, H, Amb, _g(i,0), mbFact(i,0), mbFact(i,0), -1);
lambrechts's avatar
lambrechts committed
854
          double eta  = H-h;
David Vincent's avatar
David Vincent committed
855
          Fun = -_g(i,0)*eta*J - u*u*J;
lambrechts's avatar
lambrechts committed
856
857
858
859
860
861
862
863
864
865
866
          Fut = -u*v*J;
          Feta = -(h+eta)*u*J/Amb;
        }
        val(i,0) = Feta;
        val(i,1) = Fun * nx - Fut * ny;
        val(i,2) = Fun * ny + Fut * nx;
      }
    }
  };
  term _term;
public:
867
  boundaryForcedDischarge(dgConservationLawShallowWater2d *claw, dgDofContainer *solution, const functor *discharge, std::string tag) : _term(claw->isLinear(), claw->getRadius(), claw->getBathymetry(), solution, discharge, tag, claw->_movingBathFactor, claw->_g) {
lambrechts's avatar
lambrechts committed
868
869
870
871
    _term0 = &_term;
  }
};

872
dgBoundaryCondition *dgConservationLawShallowWater2d::newForcedDischarge(dgDofContainer *solution, const functor *discharge, std::string tag){
lambrechts's avatar
lambrechts committed
873
874
875
876
  return new boundaryForcedDischarge(this, solution, discharge, tag);
}

// BC : Coupling SW1D
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
//
//class dgConservationLawShallowWater2d::boundaryCouplingSW1D:public dgBoundaryCondition {
//  class term:public function {
//    fullMatrix<double> normals, solExt, sol, bath, _xyz;
//    double _R, _g;
//  public:
//    term(const functor *solExtF, const functor *bathF,const double R, double g):function(3){
//      setArgument(solExt, solExtF);
//      setArgument(sol, function::getSolution());
//      setArgument(bath, bathF);
//      setArgument(normals, function::getNormals());
//      _R = R;
//      _g = g;
//      if (_R>0)
//        setArgument(_xyz, function::getCoordinates());
//    }
//    void call (dataCacheMap *map, fullMatrix<double> &val) { 
//      double J = 1;
//      for (int i=0; i<sol.size1(); i++){
//        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);
//	  J = 1.0;
//        }
//        double nx = normals(i,0), ny = normals(i,1);
//        double uIn = nx * sol(i,1) + ny * sol(i,2);
//        double vIn = -ny * sol(i,1) + nx * sol(i,2); 
//        double uExt = solExt(i,1);
//        double vExt = 0.0;
//        double etaIn = sol(i,0), etaExt = solExt(i,0);
//        double h = bath(i,0);
//        double HIn = h + etaIn;
//        double HExt = h + etaExt;
//        double uStar,vStar,HStar, AStar;
//        //printf("2D u = %.5e, %.5e  H = %.5e %.5e eta %.5e %.5e \n",uIn, uExt, HIn, HExt, etaIn, etaExt);
//        roeSolver(uIn, uExt, vIn, vExt, HIn, HExt, uStar, vStar, HStar, AStar, _g);
//        double etaStar = HStar - h;
//        double Feta = - HStar * uStar * J;
//        double Fun = - _g * etaStar * J - uStar * uStar * J;
//        double Fut = - uStar * vStar * J;
//        val(i,0) = Feta;
//        val(i,1) = Fun * nx - Fut * ny;
//        val(i,2) = Fun * ny + Fut * nx;
//      }
//    }
//  };
//public:
//  term  _term;
//  boundaryCouplingSW1D (dgConservationLawShallowWater2d *claw, const functor *solExtF) : _term(solExtF, claw->getBathymetry(), claw->getRadius(), claw->_g) {
//    _term0 = &_term;
//  }
//};
//
//dgBoundaryCondition *dgConservationLawShallowWater2d::newBoundaryCouplingSW1D(const functor *solExtF){
//  return new boundaryCouplingSW1D(this, solExtF);
//}
934
935
936
937
938
939
940
941
942
/*==============================================================================
 * IP Term
 *============================================================================*/

class dgIPTermOnSphere : public functor {
  const functor *_diffusiveFlux, *_diffusivity;
  bool _isotropic;
  public:
  void operator()(functorCache &cm, fullMatrix<double> &val) const {
943
    const fullMatrix<double> &diffusiveFluxL = cm.get(*_diffusiveFlux, 0); //derivatives nu*(deta_dx,du_dx,dv_dx,deta_dy,du_dy,dv_dy)
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
    const fullMatrix<double> &diffusiveFluxR = cm.get(*_diffusiveFlux, 1);
    const fullMatrix<double> &diffusivityL = cm.get(*_diffusivity, 0);
    const fullMatrix<double> &diffusivityR = cm.get(*_diffusivity, 1);
    const fullMatrix<double> &solutionL  = cm.get(cm.solution(), 0);
    const fullMatrix<double> &solutionR  = cm.get(cm.solution(), 1);
    const fullMatrix<double> &normalsL = cm.get(*function::getNormals(),0);
    const fullMatrix<double> &normalsR = cm.get(*function::getNormals(),1);
    int nbFields = solutionL.size2();

    fullMatrix<double> ipf = dgConservationLawFunction::muFactor(&cm);
    val.resize(cm.nEvaluationPoint(), nbFields * 2);
    double nxR, nyR, nxL, nyL;
    for (int iPt = 0; iPt < cm.nEvaluationPoint(); ++iPt) {
      // the sign of the right normals changes because the basis are different (except for element on the boundary) 
      if (cm.interfaceGroup()->nConnection() ==1){
        nxR = normalsR(iPt,0);
        nyR = normalsR(iPt,1);  
      }
      else {
        nxR = -1.*normalsR(iPt,0);
        nyR = -1.*normalsR(iPt,1);
      }    
      nxL = normalsL(iPt,0);
      nyL = normalsL(iPt,1);
      
      double mufactor = ipf(iPt, 0);
      // viscosities for u and v
      double nuU = std::max(diffusivityR(iPt, 1), diffusivityL(iPt, 1));
      double nuV = std::max(diffusivityR(iPt, 2), diffusivityL(iPt, 2));
      // left velocities
      double unL = nuU * solutionL(iPt,1)*nxL + nuV * solutionL(iPt,2)*nyL;
      double utL = nuV * solutionL(iPt,2)*nxL - nuU * solutionL(iPt,1)*nyL;
976
977
978
979
      double unxL = diffusiveFluxL(iPt,1)*nxL + diffusiveFluxL(iPt,2)*nyL;
      double unyL = diffusiveFluxL(iPt,4)*nxL + diffusiveFluxL(iPt,5)*nyL;
      double utxL = diffusiveFluxL(iPt,2)*nxL - diffusiveFluxL(iPt,1)*nyL;
      double utyL = diffusiveFluxL(iPt,5)*nxL - diffusiveFluxL(iPt,4)*nyL;
980
981
982
983
984
985
986
      
      double unnL = unxL*nxL + unyL*nyL;
      double utnL = utxL*nxL + utyL*nyL;
      
      // right velocities
      double unR = nuU * solutionR(iPt,1)*nxR + nuV * solutionR(iPt,2)*nyR;
      double utR = nuV * solutionR(iPt,2)*nxR - nuU * solutionR(iPt,1)*nyR;
987
988
989
990
      double unxR = diffusiveFluxR(iPt,1)*nxR + diffusiveFluxR(iPt,2)*nyR;
      double unyR = diffusiveFluxR(iPt,4)*nxR + diffusiveFluxR(iPt,5)*nyR;
      double utxR = diffusiveFluxR(iPt,2)*nxR - diffusiveFluxR(iPt,1)*nyR;
      double utyR = diffusiveFluxR(iPt,5)*nxR - diffusiveFluxR(iPt,4)*nyR;
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
      
      double unnR = unxR*nxR + unyR*nyR;
      double utnR = utxR*nxR + utyR*nyR;
      
      //mean velocities and jump
      double meanNormalFlux = 0.5*(unnR+unnL);
      double meanTangentFlux = 0.5*(utnR+utnL);
      double jumpNormalVel = 0.5*(unL-unR)* mufactor;// nu is in un and ut
      double jumpTangentVel = 0.5*(utL-utR)* mufactor;
      
      double Fun = meanNormalFlux + jumpNormalVel;
      double Fut = meanTangentFlux + jumpTangentVel;
      
1004
1005
      double FuL = Fun * nxL - Fut * nyL;
      double FvL = Fun * nyL + Fut * nxL;
1006
      
1007
1008
1009
      double FuR = Fun * nxR - Fut * nyR;
      double FvR = Fun * nyR + Fut * nxR;

1010
1011
      //left values
      val(iPt,0) = 0.;
1012
1013
1014
      val(iPt,1) = -FuL;
      val(iPt,2) = -FvL;
      
1015
1016
      //right values
      val(iPt,3) = 0.;
1017
1018
1019
      val(iPt,4) = FuR;
      val(iPt,5) = FvR;
          
1020
1021
1022
1023
1024
1025
1026
1027
    }
  }
  dgIPTermOnSphere (int nbFields, const functor *diffusiveFlux, const functor *diffusivity, bool isotropic){
    _diffusiveFlux = diffusiveFlux;
    _diffusivity = diffusivity;
    _isotropic = isotropic;
  }
};
Valentin Vallaeys's avatar
Valentin Vallaeys committed
1028

1029
1030
1031
1032
functor *dgNewIpTermIsotropicOnSphere (int nbFields, const functor *diffusiveFlux, const functor *diffusivity) {
  return new dgIPTermOnSphere (nbFields, diffusiveFlux, diffusivity, true);
}