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);
}