dgConservationLawShallowWater2d.cpp 34.8 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 {
Valentin Vallaeys's avatar
Valentin Vallaeys committed
92
  fullMatrix<double> bath, sol, _g;
lambrechts's avatar
lambrechts committed
93
94
  bool _linear, _filterLinear;
 public:
Valentin Vallaeys's avatar
Valentin Vallaeys committed
95
  maxConvectiveSpeed (const functor *bathymetry, bool linear, bool filterLinear, const functor *g):function(1){
lambrechts's avatar
lambrechts committed
96
97
98
99
    setArgument(bath,bathymetry);
    setArgument(sol,function::getSolution());
    _linear = linear;
    _filterLinear = filterLinear;
David Vincent's avatar
David Vincent committed
100
    setArgument(_g, g);
lambrechts's avatar
lambrechts committed
101
102
103
104
  }
  void call(dataCacheMap *m, fullMatrix<double> &val) {
    if (_linear){
      if (_filterLinear)
105
        val.setAll(0);
106
      else{
107
        for (int i=0; i<val.size1(); i++)
David Vincent's avatar
David Vincent committed
108
          val(i,0) = sqrt(bath(i,0)*_g(i,0));
109
          }
lambrechts's avatar
lambrechts committed
110
    }else
111
112
      if (!_filterLinear){
        for (int i=0; i<val.size1(); i++) {
David Vincent's avatar
David Vincent committed
113
          val(i,0) = sqrt((bath(i,0)+sol(i,0))*_g(i,0))+hypot(sol(i,1),sol(i,2));
114
115
116
        }
      }
      else{
117
        for (int i=0; i<val.size1(); i++)
David Vincent's avatar
David Vincent committed
118
          val(i,0) = sqrt(sol(i,0)*_g(i,0))+hypot(sol(i,1),sol(i,2));  }
lambrechts's avatar
lambrechts committed
119
120
121
  }
};

Valentin Vallaeys's avatar
Valentin Vallaeys committed
122
123
124
125
126
127
128
129
130
131
132
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> &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
133
134
135
136
137
138
    for(size_t i=0; i< nQP; i++) {
      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
139
      if (!_linear){
lambrechts's avatar
lambrechts committed
140
        H+=eta;
141
142
      }
      else if (_linearSolverFrom3D){
Valentin Vallaeys's avatar
Valentin Vallaeys committed
143
        H+=(*etaOld)(i,0);
144
      }
lambrechts's avatar
lambrechts committed
145
      // flux_x
Valentin Vallaeys's avatar
Valentin Vallaeys committed
146
147
      val(i,0) = H*u /A;
      val(i,1) = g(i,0)*(1+rhoSurf(i,0)/_density)*eta + atmPress(i,0)/_density;
lambrechts's avatar
lambrechts committed
148
149
      val(i,2) = 0;
      // flux_y
Valentin Vallaeys's avatar
Valentin Vallaeys committed
150
      val(i,3) = H*v /A;
lambrechts's avatar
lambrechts committed
151
      val(i,4) = 0;
Valentin Vallaeys's avatar
Valentin Vallaeys committed
152
      val(i,5) = g(i,0)*(1+rhoSurf(i,0)/_density) *eta + atmPress(i,0)/_density;
Valentin Vallaeys's avatar
Valentin Vallaeys committed
153
      if (!_linear && !_from3D) {
Valentin Vallaeys's avatar
Valentin Vallaeys committed
154
155
156
157
        val(i,1) += u*u;
        val(i,2) += u*v;
        val(i,4) += v*u;
        val(i,5) += v*v;
lambrechts's avatar
lambrechts committed
158
      }
Valentin Vallaeys's avatar
Valentin Vallaeys committed
159
160
161
162
163
      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
164
165
166
167
168
169
      }
      // flux_z
      val(i,6) = 0;
      val(i,7) = 0;
      val(i,8) = 0;
    }
Valentin Vallaeys's avatar
Valentin Vallaeys committed
170
}
lambrechts's avatar
lambrechts committed
171

Valentin Vallaeys's avatar
Valentin Vallaeys committed
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
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> *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
191
192
193
194
    for(size_t i=0; i< nQP; i++) {
      double eta = sol(i,0);
      double u = sol(i,1);
      double v = sol(i,2);
Valentin Vallaeys's avatar
Valentin Vallaeys committed
195
      double H = eta + bathymetry(i,0);
Valentin Vallaeys's avatar
Valentin Vallaeys committed
196
      val (i,0) = 0;
Valentin Vallaeys's avatar
Valentin Vallaeys committed
197
198
199
200
201
202
      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;
      if(_nu){
        double H = bathymetry(i,0);
        double dHdx = (*bathymetryGradient)(i, 0);
        double dHdy = (*bathymetryGradient)(i, 1);
lambrechts's avatar
lambrechts committed
203
204
        if (!_linear) {
          H += sol(i,0);
Valentin Vallaeys's avatar
Valentin Vallaeys committed
205
206
          dHdx += (*solGradient)(i,0);
          dHdy += (*solGradient)(i,1);
lambrechts's avatar
lambrechts committed
207
        }
Valentin Vallaeys's avatar
Valentin Vallaeys committed
208
209
        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
210
      }
211
      if(!_linear && !_from3D){
Valentin Vallaeys's avatar
Valentin Vallaeys committed
212
        double div = (*solGradient)(i,3) + (*solGradient)(i,7);
lambrechts's avatar
lambrechts committed
213
        double bottomFriction = hypot(u,v)*quadraticDissipation(i,0);
Valentin Vallaeys's avatar
Valentin Vallaeys committed
214
215
        val(i,1) += u*(- bottomFriction + div);
        val(i,2) += v*(- bottomFriction + div);
lambrechts's avatar
lambrechts committed
216
217
      }
      if(_useMovingBathWD){
Valentin Vallaeys's avatar
Valentin Vallaeys committed
218
219
220
221
        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
222
223
        val(i,0) += - H / ( A*A ) * ( dAdx * u + dAdy * v);
      }
Valentin Vallaeys's avatar
Valentin Vallaeys committed
224
225
226
      if (_nudgingCoeff) {
        double uref = (*nudgingVel)(i,1);
        double vref = (*nudgingVel)(i,2);
Valentin Vallaeys's avatar
Valentin Vallaeys committed
227
        if (_nudgingVelIsTransport){
Valentin Vallaeys's avatar
Valentin Vallaeys committed
228
229
          uref = uref/bathymetry(i,0);
          vref = vref/bathymetry(i,0);
Valentin Vallaeys's avatar
Valentin Vallaeys committed
230
        }
231
        //double sigma[2] = {_absLayerCoef(i,0), _absLayerCoef(i,1)};
Valentin Vallaeys's avatar
Valentin Vallaeys committed
232
        double sigma = (*nudgingCoeff)(i,0);
Valentin Vallaeys's avatar
Valentin Vallaeys committed
233
234
        val(i,1) -= sigma * (u-uref);
        val(i,2) -= sigma * (v-vref);
235
236
237
238
239
240
        // 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
241
242
243
244
        //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
245
      }
lambrechts's avatar
lambrechts committed
246
    }
Valentin Vallaeys's avatar
Valentin Vallaeys committed
247
}
lambrechts's avatar
lambrechts committed
248

Valentin Vallaeys's avatar
Valentin Vallaeys committed
249
250
251
252
253
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
254
255
256
257
258
259
260
261
    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
262
}
lambrechts's avatar
lambrechts committed
263

Valentin Vallaeys's avatar
Valentin Vallaeys committed
264
265
266
267
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
268
269
270
271
272
    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
273
}
lambrechts's avatar
lambrechts committed
274

bertrand's avatar
bertrand committed
275
276
void dgConservationLawShallowWater2d::riemann(functorCache &cache,fullMatrix<double> &val) const {
  const size_t nQP = cache.nEvaluationPoint();
277
  val.resize(nQP,6);
lambrechts's avatar
lambrechts committed
278
279
  const fullMatrix<double> &solL = cache.get(cache.solution(), 0);
  const fullMatrix<double> &solR = cache.get(cache.solution(), 1);
280
281
  const fullMatrix<double> &bathL = cache.get(*_bathymetry,0);
  const fullMatrix<double> &bathR = cache.get(*_bathymetry,1);
David Vincent's avatar
David Vincent committed
282
  const fullMatrix<double> &g = cache.get(*_g,0);
283
284
  const fullMatrix<double> &normalsL = cache.get(*function::getNormals(),0);
  const fullMatrix<double> &normalsR = cache.get(*function::getNormals(),1);
bertrand's avatar
bertrand committed
285
286
  const fullMatrix<double> &movingBathFactorL = cache.get(*_movingBathFactor,0);
  const fullMatrix<double> &movingBathFactorR = cache.get(*_movingBathFactor,1);
Valentin Vallaeys's avatar
Valentin Vallaeys committed
287
288
289
290
291
292
  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
293
294
295
296
  bool _isDiffusive = _ipTerm;
  double _upwindFactor = _upwindFactorRiemann;
  const fullMatrix<double> *ip = _isDiffusive ? &cache.get(*_ipTerm) : NULL;
  for(size_t i=0; i< nQP; i++) {
Valentin Vallaeys's avatar
Valentin Vallaeys committed
297
298
    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
299
    double Fun, Fut, Feta;
300
301
302
    double nxL = normalsL(i,0);
    double nyL = normalsL(i,1);
    double nxR,nyR;
303
    // 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)
304
305
306
307
308
309
310
311
312
313
314
    if (cache.interfaceGroup()->nConnection() ==1){
        nxR = normalsR(i,0);
        nyR = normalsR(i,1);  
    }
    else {
        nxR = -1.*normalsR(i,0);
        nyR = -1.*normalsR(i,1);
    }    
    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);
315
    double vR = -nyR*solR(i,1) + nxR*solR(i,2); 
316
317
    double hL = bathL(i,0);
    double hR = bathR(i,0);
318
    if (_laxFriedrichs) {
Valentin Vallaeys's avatar
Valentin Vallaeys committed
319
      //double Amb = (movingBathFactorL(i,0) + movingBathFactorR(i,0))/2.;
Valentin Vallaeys's avatar
Valentin Vallaeys committed
320
321
322
323
324
325
326
327
328
      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.;
329
      if (_linear) {
Valentin Vallaeys's avatar
Valentin Vallaeys committed
330
        c = sqrt(g(i,0) * H);
331
        if (_linearSolverFrom3D){
Valentin Vallaeys's avatar
Valentin Vallaeys committed
332
          H += ((*etaOldL)(i,0) + (*etaOldR)(i,0))/2;
333
        }
334
335
      }
      else {
Valentin Vallaeys's avatar
Valentin Vallaeys committed
336
        H += etaM;
337
338
        vR = vR * hR/hMin;
        vL = vL * hL/hMin;
Valentin Vallaeys's avatar
Valentin Vallaeys committed
339
        utM = (vR + vL) /2;
340
        c = sqrt(g(i,0) * H/* + fabs(unM)*/);
Valentin Vallaeys's avatar
Valentin Vallaeys committed
341
      }
Valentin Vallaeys's avatar
Valentin Vallaeys committed
342
343
344
      Fun  = - g(i,0) * (1+rhoSurf/_density) * etaM - pa/_density + c * (uR - uL) / 2;
      Feta = - H * unM + c * (etaR - etaL) / 2;
      Fut  = 0.;
Valentin Vallaeys's avatar
Valentin Vallaeys committed
345
      if (!_linear && !_from3D){
Valentin Vallaeys's avatar
Valentin Vallaeys committed
346
347
348
        Fun += - unM * unM;
        //Feta /= Amb; TODO verify if Wetting and drying and Lax Friedrichs !
        Fut += - unM * utM + c * (vR - vL) / 2; // jump only if advection
349
      }
350
351
352
353
354
355
356
      /*
      if (cache.interfaceGroup()->nConnection() == 2) {
        int el0 = cache.interfaceGroup()->element(cache.interfaceId(i/cache.nPointByElement()),0).num();
        int el1 = cache.interfaceGroup()->element(cache.interfaceId(i/cache.nPointByElement()),1).num();
        if (el0 == 560 && el1 == 578)
          printf("flux * H in 2d is %.16g eta is %.16g %.16g H (%.16g %.16g) is %.16g u is %.16g %.16g (%.16g)\n", Feta, solL(i,0), solR(i,0), H, hL, hR, uL/hMin*hL, uR/hMin*hR, uR);
      }*/
357
358
    }
    else {
Valentin Vallaeys's avatar
Valentin Vallaeys committed
359
360
      if (_from3D)
        Msg::Fatal("2D from 3D only implemented in LF");
361
362
363
364
      if (_linear) {
        //linear equations 
        double etaR = solR(i,0);
        double etaL = solL(i,0);
365
        const double hMin = std::min(hL, hR);
David Vincent's avatar
David Vincent committed
366
        double sq_g_h = sqrt(g(i,0)/hMin);
367
368
        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
369
        Fun = -g(i,0)*(1+rhoSurf/_density)*etaRiemann + pa/_density;
370
        Fut = 0;
Valentin Vallaeys's avatar
Valentin Vallaeys committed
371
        Feta = -hMin*unRiemann;
372
      }else{
373
        const double hMin = std::min(hL, hR);
374
375
        double HR = solR(i,0) + hMin, HL = solL(i,0) + hMin;
        if ( fabs(hL-hR) > 1e-3 )
376
          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);
377
378
        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
379
        roeSolver(uL, uR, vL, vR, HL, HR, u, v, H, Amb, g(i,0), movingBathFactorL(i,0), movingBathFactorR(i,0), _upwindFactor);
380
        double eta  = H-hMin;
Valentin Vallaeys's avatar
Valentin Vallaeys committed
381
382
383
        Fun = -g(i,0)*(1+rhoSurf/_density)*eta + pa/_density - u*u;
        Fut = -u*v;
        Feta = -H*u/Amb;
384
      }
bertrand's avatar
bertrand committed
385
386
    }
    val(i,0) = Feta;
387
388
    val(i,1) = Fun * nxL - Fut * nyL;
    val(i,2) = Fun * nyL + Fut * nxL;
389
    val(i,3) = -val(i,0);
390
391
    val(i,4) = -(Fun * nxR - Fut * nyR);
    val(i,5) = -(Fun * nyR + Fut * nxR);
bertrand's avatar
bertrand committed
392
393
394
    if (_isDiffusive) {
      val(i,1)+=(*ip)(i,1);
      val(i,2)+=(*ip)(i,2);
Valentin Vallaeys's avatar
Valentin Vallaeys committed
395
396
      val(i,4)+=(*ip)(i,4);
      val(i,5)+=(*ip)(i,5);
lambrechts's avatar
lambrechts committed
397
398
    }
  }
bertrand's avatar
bertrand committed
399
}
lambrechts's avatar
lambrechts committed
400
401
402

void roeSolver(double uL, double uR, double vL, double vR, double HL, double HR,
               double &uStar, double &vStar, double &HStar, double &AStar,
403
404
               double g, double mbFactL, double mbFactR , double upwindFact){
  double _g = g;
lambrechts's avatar
lambrechts committed
405
406
407
408
409
410
411
  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);
412
  double c_roe = sqrt(_g*HM);
lambrechts's avatar
lambrechts committed
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
  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); 
}

/*==============================================================================
 * 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;
437
438
  _hydro = NULL;
  _hydroGradient = NULL;
lambrechts's avatar
lambrechts committed
439
440
441
442
443
444
  _originalBathymetry = NULL;
  _originalBathymetryGradient = NULL;
  _linearDissipation = _fzero;
  _coriolisFactor = _fzero;
  _quadraticDissipation = _fzero;
  _source = _fzerov;
445
  _windStress = _fzerov;
lambrechts's avatar
lambrechts committed
446
447
  _linear = false;
  _constantJac = false;
448
449
  static const functor *gdefault = new functionConstant(9.80616);
  _g = gdefault;
450
  _density = 1025;
lambrechts's avatar
lambrechts committed
451
452
453
454
455
456
457
458
  _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
459
460
461
462
 // _imexMode = IMEX_ALL;
 // _linearIMEX=false;
 // _constantJacIMEX=false;
 // _linearFilter=false;
lambrechts's avatar
lambrechts committed
463
  _laxFriedrichs = false;
464
465
  _nudgingCoeff = NULL;
  _nudgingVel = NULL;
Valentin Vallaeys's avatar
Valentin Vallaeys committed
466
  _nudgingVelIsTransport = false;
Valentin Vallaeys's avatar
Valentin Vallaeys committed
467
468
469
470
471
472
  _atmPress = _fzero;
  _rhoSurf = _fzero;
  _rhoSurfGrad = _fzerov;
  _from3D = false;
  _linearSolverFrom3D = false;
  _etaOld = NULL;
lambrechts's avatar
lambrechts committed
473
474
}

Valentin Vallaeys's avatar
Valentin Vallaeys committed
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
//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
491

Valentin Vallaeys's avatar
Valentin Vallaeys committed
492
493
494
495
496
497
498
//void dgConservationLawShallowWater2d::setLinearFilterMode(bool set) {
//  if (set)
//    _maximumConvectiveSpeed[""] = _maxSpeedFilt;
//  else
//    _maximumConvectiveSpeed[""] = _maxSpeed;
//  _linearFilter=set;
//}
lambrechts's avatar
lambrechts committed
499
500

void dgConservationLawShallowWater2d::setup() {
Valentin Vallaeys's avatar
Valentin Vallaeys committed
501
  if(_nu && _bathymetry && !_bathymetryGradient) Msg::Fatal("The gradient of the bathymetry is missing since we want to use diffusion with a bathymetry");
502
  _diffusivity[""] = _nu;
Valentin Vallaeys's avatar
Valentin Vallaeys committed
503
504
  _diffusiveFlux = _nu ? newFunctorMember(*this, &dgConservationLawShallowWater2d::diffusiveFlux) : NULL;
  _diffusion = _nu ? newFunctorMember(*this, &dgConservationLawShallowWater2d::diffusivity) : NULL;
505
  _ipTerm = _nu ? dgNewIpTermIsotropicOnSphere(3, _diffusiveFlux, _diffusion) : NULL;
lambrechts's avatar
lambrechts committed
506
  //_volumeTerm0[""] = new source(_linearDissipation, _quadraticDissipation, _coriolisFactor, _source, _linear);
Valentin Vallaeys's avatar
Valentin Vallaeys committed
507
508
  _sourceTerm = newFunctorMember(*this,&dgConservationLawShallowWater2d::source);
  _gradPsiTerm = newFunctorMember(*this,&dgConservationLawShallowWater2d::gradPsiTerm);
bertrand's avatar
bertrand committed
509
  _riemannTerm = newFunctorMember(*this,&dgConservationLawShallowWater2d::riemann);
Valentin Vallaeys's avatar
Valentin Vallaeys committed
510
  _maxSpeed = new maxConvectiveSpeed(_bathymetry, _linear, false, _g);
Valentin Vallaeys's avatar
Valentin Vallaeys committed
511
512
513
514
515
516
517
  _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);
Valentin Vallaeys's avatar
Valentin Vallaeys committed
518
//  _maxSpeedFilt = new maxConvectiveSpeed(_bathymetry, _linear, true, _g);
Valentin Vallaeys's avatar
Valentin Vallaeys committed
519
520
//  setImexMode(_imexMode);
//  setLinearFilterMode(_linearFilter);
lambrechts's avatar
lambrechts committed
521
522
523
524
525
526
  //_interfaceTerm0[""] = _nu ? : new riemann(_bathymetry,NULL, _linear);


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

527
void dgConservationLawShallowWater2d::setMovingBathWettingDrying(const functor *hydro, const functor *hydroGradient, double alphaMovingBathWD){
528
529
530
531
  _useMovingBathWD = true;
  _alphaMovingBathWD = alphaMovingBathWD;
  _originalBathymetry = _bathymetry;
  _originalBathymetryGradient= _bathymetryGradient;
532
533
534
535
536
537
  _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);
538
539
}

lambrechts's avatar
lambrechts committed
540
541
542
543
544
545
546
547
548
549
550
551
552
553
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
554
  if (_bathymetryGradient)  delete _bathymetryGradient;
lambrechts's avatar
lambrechts committed
555
556
557
558
559
560
561
562
563
  }
}

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

// BC : Wall

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
void dgConservationLawShallowWater2d::dgConservationLawShallowWater2dWallExtValue(functorCache &cm, fullMatrix<double> &val, double slip_factor) 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);
  for (size_t i=0; i< nQP; i++) {
    const double nx = normals (i,0);
    const double ny = normals (i,1);
    val(i, 0) = solIn(i, 0);
    double un = solIn (i,1) * nx + solIn (i,2) * ny;
    double ut = -solIn (i,1)*ny + solIn (i,2) * nx;
    val (i,1) = -un*nx-slip_factor*ut*ny;
    val (i,2) = -un*ny+slip_factor*ut*nx;
  }
Valentin Vallaeys's avatar
Valentin Vallaeys committed
579
}
580

Valentin Vallaeys's avatar
Valentin Vallaeys committed
581
void dgConservationLawShallowWater2d::dgConservationLawShallowWater2dWallExtValueGrad(functorCache &cm, fullMatrix<double> &val) const {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
582
   // const fullMatrix<double> &normals = cm.get(*function::getNormals(), 0);
Valentin Vallaeys's avatar
Valentin Vallaeys committed
583
584
585
586
    const fullMatrix<double> &solGradIn = cm.get(*function::getSolutionGradient(), 0);
    size_t nQP = cm.nEvaluationPoint();
    val.resize(nQP, 9, false);
    //cm.setDefaultToSide0(true);
587
    for (size_t i=0; i< nQP; i++) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
588
589
      //const double nx = normals (i,0);
      //const double ny = normals (i,1);
590
591
592
593
594
      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);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
595
/*
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
      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;
      }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
621
      else {*/
622
623
624
625
        val (i,3) = solGradIn(i,3); 
        val (i,4) = solGradIn(i,4);
        val(i, 6) = solGradIn(i, 6);
        val(i, 7) = solGradIn(i, 7);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
626
      //}
627
    }
Valentin Vallaeys's avatar
Valentin Vallaeys committed
628
}
lambrechts's avatar
lambrechts committed
629

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
630
dgBoundaryCondition *dgConservationLawShallowWater2d::newBoundaryWall(double slip_factor){
631
  //return new boundaryWall(this);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
632
  functor *boundaryWallUV = newFunctorMember(*this, &dgConservationLawShallowWater2d::dgConservationLawShallowWater2dWallExtValue, slip_factor);
Valentin Vallaeys's avatar
Valentin Vallaeys committed
633
  functor *boundaryWallGradUV = newFunctorMember(*this, &dgConservationLawShallowWater2d::dgConservationLawShallowWater2dWallExtValueGrad);
634
635
  std::vector<const functor*> *toReplace = new std::vector<const functor*>();
  std::vector<const functor*> *replaceBy = new std::vector<const functor*>();
636
637
  toReplace->push_back(function::getSolutionGradient());
  replaceBy->push_back(boundaryWallGradUV);
638
639
640
  toReplace->push_back(function::getSolution());
  replaceBy->push_back(boundaryWallUV);
  return newOutsideValueBoundaryGeneric("",*toReplace,*replaceBy);
lambrechts's avatar
lambrechts committed
641
642
643
644
645
646
}

// BC : Forced Discharge

class dgConservationLawShallowWater2d::boundaryForcedDischarge : public dgBoundaryCondition {
  class term : public function {
Valentin Vallaeys's avatar
Valentin Vallaeys committed
647
    fullMatrix<double> normals, sol, bath, discharge, mbFact, _g;
lambrechts's avatar
lambrechts committed
648
    bool _linear;
Valentin Vallaeys's avatar
Valentin Vallaeys committed
649
    double S0;
lambrechts's avatar
lambrechts committed
650
  public:
Valentin Vallaeys's avatar
Valentin Vallaeys committed
651
    term(bool linear, const functor *bathF, dgDofContainer *solution, const functor *dischargeF, std::string tag, const functor *movingBathF, const functor *g):function(3){
lambrechts's avatar
lambrechts committed
652
653
654
655
656
657
      setArgument(sol,function::getSolution());
      setArgument(normals,function::getNormals());
      setArgument(bath,bathF);
      setArgument(discharge,dischargeF);
      setArgument(mbFact,movingBathF);
      _linear = linear;
David Vincent's avatar
David Vincent committed
658
      setArgument(_g, g);
lambrechts's avatar
lambrechts committed
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
      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 Fun, Fut, Feta;
      for (int i=0; i<sol.size1(); i++) {
        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
675
          double sq_g_h = sqrt(_g(i,0)/h);
lambrechts's avatar
lambrechts committed
676
677
          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
678
          Fun = -_g(i,0)*etaRiemann;
lambrechts's avatar
lambrechts committed
679
          Fut = 0;
Valentin Vallaeys's avatar
Valentin Vallaeys committed
680
          Feta = -h*unRiemann;
lambrechts's avatar
lambrechts committed
681
682
683
684
        } 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
685
          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
686
          double eta  = H-h;
Valentin Vallaeys's avatar
Valentin Vallaeys committed
687
688
689
          Fun = -_g(i,0)*eta - u*u;
          Fut = -u*v;
          Feta = -(h+eta)*u/Amb;
lambrechts's avatar
lambrechts committed
690
691
692
693
694
695
696
697
698
        }
        val(i,0) = Feta;
        val(i,1) = Fun * nx - Fut * ny;
        val(i,2) = Fun * ny + Fut * nx;
      }
    }
  };
  term _term;
public:
Valentin Vallaeys's avatar
Valentin Vallaeys committed
699
  boundaryForcedDischarge(dgConservationLawShallowWater2d *claw, dgDofContainer *solution, const functor *discharge, std::string tag) : _term(claw->isLinear(), claw->getBathymetry(), solution, discharge, tag, claw->_movingBathFactor, claw->_g) {
lambrechts's avatar
lambrechts committed
700
701
702
703
    _term0 = &_term;
  }
};

704
dgBoundaryCondition *dgConservationLawShallowWater2d::newForcedDischarge(dgDofContainer *solution, const functor *discharge, std::string tag){
lambrechts's avatar
lambrechts committed
705
706
707
708
  return new boundaryForcedDischarge(this, solution, discharge, tag);
}

// BC : Coupling SW1D
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
//
//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);
//}
766
767
768
769
770
771
772
773
774
/*==============================================================================
 * IP Term
 *============================================================================*/

class dgIPTermOnSphere : public functor {
  const functor *_diffusiveFlux, *_diffusivity;
  bool _isotropic;
  public:
  void operator()(functorCache &cm, fullMatrix<double> &val) const {
775
    const fullMatrix<double> &diffusiveFluxL = cm.get(*_diffusiveFlux, 0); //derivatives nu*(deta_dx,du_dx,dv_dx,deta_dy,du_dy,dv_dy)
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
    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;
808
809
810
811
      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;
812
813
814
815
816
817
818
      
      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;
819
820
821
822
      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;
823
824
825
826
827
828
829
830
831
832
833
834
835
      
      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;
      
836
837
      double FuL = Fun * nxL - Fut * nyL;
      double FvL = Fun * nyL + Fut * nxL;
838
      
839
840
841
      double FuR = Fun * nxR - Fut * nyR;
      double FvR = Fun * nyR + Fut * nxR;

842
843
      //left values
      val(iPt,0) = 0.;
844
845
846
      val(iPt,1) = -FuL;
      val(iPt,2) = -FvL;
      
847
848
      //right values
      val(iPt,3) = 0.;
849
850
851
      val(iPt,4) = FuR;
      val(iPt,5) = FvR;
          
852
853
854
855
856
857
858
859
    }
  }
  dgIPTermOnSphere (int nbFields, const functor *diffusiveFlux, const functor *diffusivity, bool isotropic){
    _diffusiveFlux = diffusiveFlux;
    _diffusivity = diffusivity;
    _isotropic = isotropic;
  }
};
Valentin Vallaeys's avatar
Valentin Vallaeys committed
860

861
862
863
864
functor *dgNewIpTermIsotropicOnSphere (int nbFields, const functor *diffusiveFlux, const functor *diffusivity) {
  return new dgIPTermOnSphere (nbFields, diffusiveFlux, diffusivity, true);
}