slim3dEquations.cpp 46.2 KB
Newer Older
lambrechts's avatar
lambrechts committed
1
2
3
4
#include "slim3dEquations.h"
#include "slim3dSolver.h"
#include "functionGeneric.h"
#include "dgSW3dVerticalModel.h"
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
5
#include "dgSW3dTurbulenceGOTM.h"
lambrechts's avatar
lambrechts committed
6
7
8
9
#include "dgExtrusion.h"

slim3dSolverDofs::slim3dSolverDofs(slim3dSolver *solver): _solver(solver) {
  // init all dofs to NULL
delandmeter's avatar
delandmeter committed
10
  uDof = vDof = uvDof = uvDofCenter = uvSurfDof2d = uvAvDof2d = NULL;
lambrechts's avatar
lambrechts committed
11
  uvHorDivDof = NULL; 
delandmeter's avatar
delandmeter committed
12
  uvIntDof2d = uvCorrDof3d = uvBotDof2d = NULL;
lambrechts's avatar
lambrechts committed
13
  windStressDof2dCenter = NULL;
14
  uvTauBDof2d = uvTauBDof2dCenter = NULL;
delandmeter's avatar
delandmeter committed
15
16
  wDof3d = wSurfDof2d = NULL;
  //wMeshDof3d  = wMeshSurfDof2d = NULL;
delandmeter's avatar
delandmeter committed
17
  //dwMeshDzDof3d = NULL;
18
  etaDof2d = etaDof2dCG = etaDof3dCenter = NULL;
lambrechts's avatar
lambrechts committed
19
  rDof3d = rIntDof2d = rBotDof2d = NULL;
karna's avatar
karna committed
20
  rhoDof3d = rGradDof3d = rhoGradDof3d = rhoSurfDof2d = rGradIntDof2d = NULL;
delandmeter's avatar
delandmeter committed
21
  SDof = SDofCenter = TDof = TDofCenter = SBottomDof2d= NULL;
22
23
  //sedDof = sedBottomDof2d = sedGroundDof2d = sedBottomFluxDof2d = NULL;
  sedDof.resize(0); sedBottomDof2d.resize(0); sedGroundDof2d.resize(0); sedBottomFluxDof2d.resize(0);
delandmeter's avatar
delandmeter committed
24
  zDof = deltaZDof3d = deltaZDofSurf2d = zDofCenter = zBotDof2d = NULL;
lambrechts's avatar
lambrechts committed
25
  uvDeviationProdIntDof2d = NULL;
delandmeter's avatar
delandmeter committed
26
  uvDeviationProdDof3d = NULL;
lambrechts's avatar
lambrechts committed
27
28
  xyzOrigDof = nnDof = ssDof = nuvDof = nuvDofCenter = kappavDof = kappavDofCenter = NULL;
  tkeDof = tkeDofCenter = epsDof = epsDofCenter = lDof = NULL;
29
  uvGradDof  = uvGradDof2d= NULL;
30
  nbDof3d = rhoGradGMVelDof3d = GMVelDof3d = streamFuncDof3d = NULL;
31
  bathDof2d = bathCGMinDof2d = NULL;
32
  arbitraryMovingMeshDof3d = kappaZEqDof3d = NULL;
33
  vJumpDof3d = errorDof3d = zStarDof3d = oneOverErrorDof3d = oneOverErrorIntDof2d = NULL;
34
  testDof3d = testDof2d = NULL;
lambrechts's avatar
lambrechts committed
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
  _initialized = false;
}

void slim3dSolverDofs::initialize()
{
  if (_initialized)
    return;
  dgGroupCollection *groups3d = _solver->groups3d;
  dgGroupCollection *groups2d = _solver->groups2d;
  // TODO allocate according to solver options
  uvDof = new dgDofContainer(*groups3d, 2);
  uvHorDivDof = new dgDofContainer(*groups3d, 1);
  wDof3d = new dgDofContainer(*groups3d, 1);
  //wDof3dCopy = new dgDofContainer(*groups3d, 1);
  uvAvDof2d = new dgDofContainer(*groups2d, 2);
  uvIntDof2d = new dgDofContainer(*groups2d, 2);
delandmeter's avatar
delandmeter committed
51
52
  uvCorrDof3d = new dgDofContainer(*groups3d, 2);
  uvDeviationProdDof3d = new dgDofContainer(*groups3d, 3);
lambrechts's avatar
lambrechts committed
53
54
55
56
57
  uvDeviationProdIntDof2d = new dgDofContainer(*groups2d, 3); // this is needed
  etaDof2d = new dgDofContainer(*groups2d, 1);
  etaDof2dCG = new dgDofContainer(*groups2d, 1);
  // pressure related
  if (true){
58
    rGradDof3d = new dgDofContainer(*groups3d, 2); // rGradDof is not an appropriate name. It is actually the integral of the density deviation horizontal gradient
lambrechts's avatar
lambrechts committed
59
60
61
62
    rhoGradDof3d = new dgDofContainer(*groups3d, 2);
    rhoDof3d = new dgDofContainer(*groups3d, 1);
    rhoSurfDof2d = new dgDofContainer(*groups2d, 1);
    rGradIntDof2d = new dgDofContainer(*groups2d, 2);
63
64
65
    dRhodzDof3d = new dgDofContainer(*groups3d, 3);
    intDRhodzDof3d = new dgDofContainer(*groups3d, 1);
    intDRhodzDofSurf2d = new dgDofContainer(*groups2d, 1);
lambrechts's avatar
lambrechts committed
66
67
68
69
70
71
72
73
74
75
  }
  else {
    rDof3d = new dgDofContainer(*groups3d, 1);
    rIntDof2d = new dgDofContainer(*groups2d, 1);
    rBotDof2d = new dgDofContainer(*groups2d, 1);
  }
  // bottom friction related
  uvBotDof2d = new dgDofContainer(*groups2d, 2);
  uvTauBDof2d = new dgDofContainer(*groups2d, 2);
  zBotDof2d = new dgDofContainer(*groups2d, 1);
76
  if ( _solver->_solveUVDiff || _solver->_solveUVAdvVert) {
lambrechts's avatar
lambrechts committed
77
78
79
80
81
82
    vDof = new dgDofContainer(*groups3d, 1);
    uDof = new dgDofContainer(*groups3d, 1);
  }
  // moving mesh related
  uvSurfDof2d = new dgDofContainer(*groups2d, 2);
  wSurfDof2d = new dgDofContainer(*groups2d, 1);
delandmeter's avatar
delandmeter committed
83
  //wMeshSurfDof2d = new dgDofContainer(*groups2d, 1);
84
  //wMeshDof3d = new dgDofContainer(*groups3d, 1);
delandmeter's avatar
delandmeter committed
85
  //dwMeshDzDof3d = new dgDofContainer(*groups3d, 3);
delandmeter's avatar
delandmeter committed
86
  //wMeshSurfDof2d->setAll(0.);
87
  //wMeshDof3d->setAll(0.);
lambrechts's avatar
lambrechts committed
88
89
  xyzOrigDof = new dgDofContainer(*groups3d, 3);
  zDof = new dgDofContainer(*groups3d, 1);
90
91
  deltaZDof3d = new dgDofContainer(*groups3d, 1);
  deltaZDofSurf2d = new dgDofContainer(*groups2d, 1);
lambrechts's avatar
lambrechts committed
92
  if ( _solver->_solveS || _solver->getSolveSImplicitVerticalDiffusion() ||
93
      _solver->_solveTurbulence || _solver->getSolveGMVel() ) {
lambrechts's avatar
lambrechts committed
94
95
96
    SDof = new dgDofContainer(*groups3d, 1);
  }
  if ( _solver->_solveT || _solver->getSolveTImplicitVerticalDiffusion() ||
97
      _solver->_solveTurbulence || _solver->getSolveGMVel()) {
lambrechts's avatar
lambrechts committed
98
99
100
    TDof = new dgDofContainer(*groups3d, 1);
  }
  if ( _solver->_solveSed || _solver->getSolveSedImplicitVerticalDiffusion() ){
101
102
103
104
105
106
107
108
109
110
    sedDof.resize(_solver->getNumSedEq());
    sedBottomDof2d.resize(_solver->getNumSedEq());
    sedGroundDof2d.resize(_solver->getNumSedEq());
    sedBottomFluxDof2d.resize(_solver->getNumSedEq());
    for (int iSedEq = 0; iSedEq < _solver->getNumSedEq(); ++iSedEq){
      sedDof[iSedEq] = new dgDofContainer(*groups3d, 1);
      sedBottomDof2d[iSedEq] = new dgDofContainer(*groups2d, 1);
      sedGroundDof2d[iSedEq] =  new dgDofContainer(*groups2d, 1);
      sedBottomFluxDof2d[iSedEq] =  new dgDofContainer(*groups2d, 1);
    }
delandmeter's avatar
delandmeter committed
111
    SBottomDof2d = new dgDofContainer(*groups2d, 1);
lambrechts's avatar
lambrechts committed
112
113
  }
  if ( _solver->_solveTurbulence ) {
114
115
    nuvDof = new dgDofContainer(*groups3d, 1);
    kappavDof = new dgDofContainer(*groups3d, 1);
lambrechts's avatar
lambrechts committed
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
    TDofCenter = new dgDofContainer(*groups3d, 1);
    SDofCenter = new dgDofContainer(*groups3d, 1);
    uvDofCenter = new dgDofContainer(*groups3d, 2);
    etaDof3dCenter = new dgDofContainer(*groups3d, 1);
    windStressDof2dCenter = new dgDofContainer(*groups2d, 2);
    uvTauBDof2dCenter = new dgDofContainer(*groups2d, 2);
    zDofCenter = new dgDofContainer(*groups3d, 1);
    nnDof = new dgDofContainer(*groups3d, 1);
    ssDof = new dgDofContainer(*groups3d, 1);
    nuvDofCenter = new dgDofContainer(*groups3d, 1);
    kappavDofCenter = new dgDofContainer(*groups3d, 1);
    tkeDof = new dgDofContainer(*groups3d, 1);
    tkeDofCenter = new dgDofContainer(*groups3d, 1);
    epsDof = new dgDofContainer(*groups3d, 1);
    epsDofCenter = new dgDofContainer(*groups3d, 1);
    lDof = new dgDofContainer(*groups3d, 1);
  }
133
  if ( _solver->getSolveUVGrad()){
134
    uvGradDof = new dgDofContainer(*groups3d, 6);
135
136
    uvGradDof2d = new dgDofContainer(*groups2d, 6);
  }
lambrechts's avatar
lambrechts committed
137
138
139
140
141
142
  if ( _solver->getSolveGMVel() ) {
    rhoGradGMVelDof3d = new dgDofContainer(*groups3d, 3);
    nbDof3d = new dgDofContainer(*groups3d, 1);
    streamFuncDof3d = new dgDofContainer(*groups3d, 3);
    GMVelDof3d = new dgDofContainer(*groups3d, 3);
  }
143
  bathDof2d = new dgDofContainer(*groups2d, 1);
144
  bathCGMinDof2d = new dgDofContainer(*groups2d, 1);
145
  arbitraryMovingMeshDof3d = new dgDofContainer(*groups3d, 1);
146
  arbitraryMovingMeshDof3d->setAll(0.);
147
148
149
150
151
152
153
154
  if ( _solver->getUseAdaptiveVerticalGrid() ) {
    kappaZEqDof3d = new dgDofContainer(*groups3d, 1);
    vJumpDof3d = new dgDofContainer(*groups3d, 1);
    errorDof3d = new dgDofContainer(*groups3d, 1);
    zStarDof3d = new dgDofContainer(*groups3d, 1);
    oneOverErrorDof3d = new dgDofContainer(*groups3d, 1);
    oneOverErrorIntDof2d = new dgDofContainer(*groups2d, 1);
  }
155
  testDof3d = new dgDofContainer(*groups3d, 1);
156
  testDof2d = new dgDofContainer(*groups2d, 1);
lambrechts's avatar
lambrechts committed
157
158
159
160
161
162
163
164
165
166
  _initialized = true;
}

slim3dSolverDofs::~slim3dSolverDofs()
{
  if ( uvDof )             delete uvDof;
  if ( uvHorDivDof )       delete uvHorDivDof;
  if ( wDof3d )            delete wDof3d;
  if ( uvAvDof2d )         delete uvAvDof2d;
  if ( uvIntDof2d )        delete uvIntDof2d;
delandmeter's avatar
delandmeter committed
167
  if ( uvCorrDof3d )       delete uvCorrDof3d;
lambrechts's avatar
lambrechts committed
168
169
170
171
172
173
174
175
176
  if ( etaDof2d )          delete etaDof2d;
  if ( etaDof2dCG )        delete etaDof2dCG;
  if ( uvBotDof2d )        delete uvBotDof2d;
  if ( uvTauBDof2d )       delete uvTauBDof2d;
  if ( zBotDof2d )         delete zBotDof2d;
  if ( vDof )              delete vDof;
  if ( uDof )              delete uDof;
  if ( uvSurfDof2d )       delete uvSurfDof2d;
  if ( wSurfDof2d )        delete wSurfDof2d;
delandmeter's avatar
delandmeter committed
177
  //if ( wMeshSurfDof2d )    delete wMeshSurfDof2d;
178
  //if ( wMeshDof3d )        delete wMeshDof3d;
lambrechts's avatar
lambrechts committed
179
180
  if ( xyzOrigDof )        delete xyzOrigDof;
  if ( zDof )              delete zDof;
181
182
  if ( deltaZDof3d )       delete deltaZDof3d;
  if ( deltaZDofSurf2d )   delete deltaZDofSurf2d;
lambrechts's avatar
lambrechts committed
183
184
185
186
187
188
189
190
191
192
  if ( rDof3d )            delete rDof3d;
  if ( rIntDof2d )         delete rIntDof2d;
  if ( rBotDof2d )         delete rBotDof2d;
  if ( rGradDof3d)         delete rGradDof3d;
  if ( rhoGradDof3d)       delete rhoGradDof3d;
  if ( rhoDof3d )          delete rhoDof3d;
  if ( rhoSurfDof2d)       delete rhoSurfDof2d;
  if ( rGradIntDof2d )     delete rGradIntDof2d;
  if ( SDof )              delete SDof;
  if ( TDof )              delete TDof;
193
194
195
196
  if ( sedDof.size() !=0 )            {for (int iSedEq = 0; iSedEq < _solver->getNumSedEq(); ++iSedEq){delete sedDof[iSedEq];} sedDof.resize(0);}
  if ( sedBottomDof2d.size() !=0 )    {for (int iSedEq = 0; iSedEq < _solver->getNumSedEq(); ++iSedEq){delete sedBottomDof2d[iSedEq];} sedBottomDof2d.resize(0);}
  if ( sedGroundDof2d.size() !=0 )    {for (int iSedEq = 0; iSedEq < _solver->getNumSedEq(); ++iSedEq){delete sedGroundDof2d[iSedEq];} sedGroundDof2d.resize(0);}
  if ( sedBottomFluxDof2d.size() !=0 )    {for (int iSedEq = 0; iSedEq < _solver->getNumSedEq(); ++iSedEq){delete sedBottomFluxDof2d[iSedEq];} sedBottomFluxDof2d.resize(0);}
delandmeter's avatar
delandmeter committed
197
  if ( SBottomDof2d )      delete SBottomDof2d;
lambrechts's avatar
lambrechts committed
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
  if ( nuvDof )            delete nuvDof;
  if ( kappavDof )         delete kappavDof;
  if ( TDofCenter )        delete TDofCenter;
  if ( SDofCenter )        delete SDofCenter;
  if ( uvDofCenter )       delete uvDofCenter;
  if ( etaDof3dCenter )    delete etaDof3dCenter;
  if ( windStressDof2dCenter ) delete windStressDof2dCenter;
  if ( uvTauBDof2dCenter ) delete uvTauBDof2dCenter;
  if ( zDofCenter )        delete zDofCenter;
  if ( nnDof )             delete nnDof;
  if ( ssDof )             delete ssDof;
  if ( nuvDofCenter )      delete nuvDofCenter;
  if ( kappavDofCenter )   delete kappavDofCenter;
  if ( tkeDof )            delete tkeDof;
  if ( tkeDofCenter )      delete tkeDofCenter;
  if ( epsDof )            delete epsDof;
  if ( epsDofCenter )      delete epsDofCenter;
  if ( lDof )              delete lDof;
delandmeter's avatar
delandmeter committed
216
  if (uvDeviationProdDof3d)    delete uvDeviationProdDof3d;
lambrechts's avatar
lambrechts committed
217
218
  if (uvDeviationProdIntDof2d) delete uvDeviationProdIntDof2d;
  if ( nbDof3d )           delete nbDof3d;
219
  if ( uvGradDof )         delete uvGradDof;
220
  if ( uvGradDof2d )       delete uvGradDof2d;
lambrechts's avatar
lambrechts committed
221
222
223
  if ( rhoGradGMVelDof3d)  delete rhoGradGMVelDof3d;
  if ( streamFuncDof3d)    delete streamFuncDof3d;
  if ( GMVelDof3d)         delete GMVelDof3d;
224
  if ( bathDof2d)          delete bathDof2d;
225
  if ( bathCGMinDof2d)     delete bathCGMinDof2d;
226
  if ( arbitraryMovingMeshDof3d) delete arbitraryMovingMeshDof3d;
227
  if ( kappaZEqDof3d )     delete kappaZEqDof3d;
228
229
230
231
232
  if ( vJumpDof3d )        delete vJumpDof3d;
  if ( errorDof3d )        delete errorDof3d;
  if ( zStarDof3d )        delete zStarDof3d;
  if ( oneOverErrorDof3d )   delete oneOverErrorDof3d;
  if ( oneOverErrorIntDof2d )delete oneOverErrorIntDof2d;
233
  if ( testDof3d)          delete testDof3d;
234
  if ( testDof2d)          delete testDof2d;
235

lambrechts's avatar
lambrechts committed
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
}

slim3dSolverFunctions::slim3dSolverFunctions(slim3dSolver *solver) : _solver(solver)
{
  if (!solver->groups3d)
    Msg::Fatal("No 3d mesh loaded in solver");
  zeroFunc = new functionConstant(0.0);
  std::vector<double> zeroVec2(2,0.0);
  zeroFunc2 = new functionConstant(zeroVec2);
  std::vector<double> zeroVec(3,0.0);
  zeroFunc3 = new functionConstant(zeroVec);
  tinyFunc = new functionConstant(1e-10);
  oneFunc = new functionConstant(1.0);
  xyzFunc3d = solver->groups3d->getFunctionCoordinates();
  xyzFunc2d = solver->groups2d->getFunctionCoordinates();
  zFunc = functionExtractCompNew(xyzFunc3d,2);
  timeFunc = function::getTime();

  // functions depending on user input
  uvCorrectorFunc = NULL;
  pHydroStatFunc = rhoFunc = rhoFuncCenter = eosAlphaFunc = eosBetaFunc = NULL;
  SCenterGradientFunc = TCenterGradientFunc = SInitGradientFunc = TInitGradientFunc = NULL;
258
  bathFunc2d = bathCGMinFunc2d = NULL;
lambrechts's avatar
lambrechts committed
259
  coriolisFunc = NULL;
delandmeter's avatar
delandmeter committed
260
  TInitFunc = SInitFunc = TFunc = SFunc = sedInitFunc = sedGroundInitFunc = NULL;
lambrechts's avatar
lambrechts committed
261
262
263
  uvInitFunc = uvAvInitFunc = etaInitFunc = wInitFunc = NULL;
  uvDevFunc = uvDevProdFunc = uvDevAdvectionFunc = NULL;
  newZFunc = wMeshSurfFunc = wMeshFunc = wMeshDzFunc = zBotDistFunc = NULL;
264
  nuvFunc = nuh = nuh2d = kappavFunc = kappavS = kappavT = kappavJumpS = kappavJumpT = kappahTotal = kappahUser = kappahJumpS = kappahJumpT = kappahS = kappahT = NULL;
265
266
267
268
  maxDiffS = 35.;
  maxDiffT = 20.;
  jumpPercentS = 0.05;
  jumpPercentT = 0.05;
269
270
271
  approxDT = 1e-6; 
  coeffCFLh = 720;
  coeffCFLv = 300;
lambrechts's avatar
lambrechts committed
272
  _nuv0Func = _kapv0Func = NULL;
delandmeter's avatar
delandmeter committed
273
  z0BFunc = z0SFunc = uvTauBFunc2d = windStressFunc = windFunc = NULL;
lambrechts's avatar
lambrechts committed
274
  bottomFrictionFunc = NULL;
delandmeter's avatar
delandmeter committed
275
  bottomFrictionDragCoeff2d = bottomUVDeviation = NULL;
lambrechts's avatar
lambrechts committed
276
277
  SInitGradientFunc = TInitGradientFunc = zeroFunc3;
  slopeFunc = slopeTap = NULL;
278
  _bathDof3d = NULL;
lambrechts's avatar
lambrechts committed
279
  eadySpeedFunc = GMIndepTermFunc = GMVelFunc = NULL;
280
  wBottomFunc = NULL;
281
  eta2dFunc = uvInt2dFunc = uvAv2dFunc = rhoSurf2dFunc = rhoSurfGrad2dFunc = NULL;
282
  atmPressFunc = NULL;
283
284
  uvIntTarget2dFunc = bottomFriction2dPrecompFunc = zBotDist2dPrecompFunc = NULL;
  eta2dCGFunc = eta2dCGFuncOld = eta2dCGGradFunc = NULL;
delandmeter's avatar
delandmeter committed
285
  wMeshSurf2dPrecompFunc = NULL;
286
  wMeshPrecompFunc = wMeshDzPrecompFunc = NULL;
287
288
  wSedFunc.resize(0);
  sedBottomFluxFunc2d.resize(0);
delandmeter's avatar
delandmeter committed
289
  sedBottomFlux2dPrecompFunc = NULL;
290
  arbitraryMovingMeshFunc = kappaZEqFunc = NULL;
291
  jumpDiffTracerS = jumpDiffTracerT = jumpDiffTracerVertS = jumpDiffTracerVertT = NULL;
lambrechts's avatar
lambrechts committed
292
293
294
295
296
297
298
299
  _initialized = false;
}

void slim3dSolverFunctions::initialize()
{
  slim3dSolverDofs *dofs = _solver->getDofs();
  if (!dofs->isInitialized())
    Msg::Fatal("Dofs must be allocated before calling slim3dSolverFunctions::initialize()");
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
300
  if (! bathFunc2d) {
delandmeter's avatar
delandmeter committed
301
302
    initializeBath();
  }
lambrechts's avatar
lambrechts committed
303
304
305
306
307
308
309
  
  if (!SInitFunc && _solver->getSolveS())
    Msg::Fatal("SInitFunc must be set");
  if (!TInitFunc && _solver->getSolveT())
    Msg::Fatal("TInitFunc must be set");
  if (!sedInitFunc && _solver->getSolveSed())
    Msg::Fatal("sedInitFunc must be set");
310
311
  if (!rhoFunc && (!SInitFunc || !TInitFunc))
    Msg::Fatal("Both S and T InitFunc must be set to compute rhoFunc, or rhoFunc must be set");
lambrechts's avatar
lambrechts committed
312
313

  // set T and S functions either to the dofs (that are marched in time) or to the initial condition
314
  const functor *TCenterFunc=NULL, *SCenterFunc=NULL;
delandmeter's avatar
delandmeter committed
315
  if ( _solver->_solveS  || _solver->getSolveGMVel()) {
lambrechts's avatar
lambrechts committed
316
317
318
319
320
321
322
323
324
325
    dofs->SDof->interpolate(SInitFunc);    
    SFunc = dofs->SDof->getFunction();
    if ( _solver->_solveTurbulence ) {
      SCenterFunc = dofs->SDofCenter->getFunction();
      SCenterGradientFunc = dofs->SDofCenter->getFunctionGradient();
    }
  } else {
    SFunc = SCenterFunc = SInitFunc;
    SCenterGradientFunc = SInitGradientFunc;
  }
delandmeter's avatar
delandmeter committed
326
  if ( _solver->_solveT || _solver->getSolveGMVel()) {
lambrechts's avatar
lambrechts committed
327
328
329
330
331
332
333
334
335
336
    dofs->TDof->interpolate(TInitFunc);
    TFunc = dofs->TDof->getFunction();
    if ( _solver->_solveTurbulence ) {
      TCenterFunc = dofs->TDofCenter->getFunction();
      TCenterGradientFunc = dofs->TDofCenter->getFunctionGradient();
    }
  } else {
    TFunc = TCenterFunc = TInitFunc;
    TCenterGradientFunc = TInitGradientFunc;
  }
ylebars's avatar
ylebars committed
337
  // rho related functions
lambrechts's avatar
lambrechts committed
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
  pHydroStatFunc = zeroFunc;
  if (! rhoFunc) {
    rhoFunc = new functionStateEquation(SFunc,TFunc,pHydroStatFunc,slim3dParameters::rho0);
  }
  else if (_solver->_solveTurbulence) {
    Msg::Fatal("changing the state equation while solving the turbulence is not implemented");
  }

  if (_solver->_solveTurbulence) {
    rhoFuncCenter = new functionStateEquation( SCenterFunc , TCenterFunc,
                                           pHydroStatFunc,slim3dParameters::rho0);
    eosAlphaFunc = new functionStateEquationAlpha(SCenterFunc , TCenterFunc,
                                                 pHydroStatFunc, slim3dParameters::rho0);
    eosBetaFunc = new functionStateEquationBeta( SCenterFunc, TCenterFunc,
                                                 pHydroStatFunc, slim3dParameters::rho0);
  }
  //std::vector<int> xyComp(2,0);
  //xyComp[0] = 0;
  //xyComp[1] = 1;
  //rGradFunc = functionExtractCompNew(dofs->rGradDof3d->getFunction(), xyComp);
  // moving mesh
359
360
361
  const functor *wSurfFunc2d = dofs->wSurfDof2d->getFunction();
  const functor *uvSurfFunc2d= dofs->uvSurfDof2d->getFunction();
  //const functor *wMeshSurfDofFunc = dofs->wMeshSurfDof2d->getFunction();
delandmeter's avatar
delandmeter committed
362
  wMeshSurfFunc = new wMeshSurf(eta2dCGGradFunc, wSurfFunc2d,uvSurfFunc2d);
363
  wMeshSurf2dPrecompFunc = new functionPrecomputedExtrusion(_solver->extrusion(), _solver->groups2d->getElementGroup(0)->getOrder() * 2 + 1);//*wMeshSurfFunc);
karna's avatar
karna committed
364
  // NOTE using wMeshSurfFunc directly here fails in parallel
delandmeter's avatar
delandmeter committed
365

366
  //classical
367
368
369
  //newZFunc = new zNew(xyzOrigFunc3d,eta2dCGFunc,bathFunc2d);
  //wMeshFunc = new wMesh(xyzOrigFunc3d, eta2dCGFunc,bathFunc2d,wMeshSurf2dPrecompFunc);
  //wMeshDzFunc = new dwMeshDz(wMeshSurf2dPrecompFunc,eta2dCGFunc,bathFunc2d);
370

371
372
373
374
  if (!_solver->getUseAdaptiveVerticalGrid())
    arbitraryMovingMeshFunc = new functionConstant(0.);
  else if (_solver->getUseConservativeALE())
    Msg::Fatal("Conservative ALE does not work with adaptive vertical coordinates yet");
375
  _dt3d = new double(-1);
376
377
  //newZFunc = new zNewAdaptative(dofs->xyzOrigDof->getFunction(), eta2dCGFunc, eta2dCGFuncOld, bathFunc2d, dofs->arbitraryMovingMeshDof3d->getFunction(), _dt3d, timeFunc);
  newZFunc = new zNewAdaptative(dofs->xyzOrigDof->getFunction(), eta2dCGFunc, eta2dCGFuncOld, _solver->verticalBottomRemover->heightAboveVertBottom->getFunction(), dofs->arbitraryMovingMeshDof3d->getFunction(), _dt3d, timeFunc);
378
379
  wMeshFunc = new wMeshAdapt(wMeshSurf2dPrecompFunc, eta2dCGFunc, bathFunc2d, dofs->arbitraryMovingMeshDof3d->getFunction());
  wMeshDzFunc = new dwMeshDzAdapt(wMeshSurf2dPrecompFunc, eta2dCGFunc, bathFunc2d, dofs->arbitraryMovingMeshDof3d->getFunctionGradient());
380

lambrechts's avatar
lambrechts committed
381
382
383
384
385
  // diffusivity
  _nuv0Func = new functionConstant(slim3dParameters::nuv0);
  // TODO: differentiate between kapv0S and kapv0T ?
  _kapv0Func = new functionConstant(slim3dParameters::kapv0S);
  if (_solver->_solveTurbulence) {
386
387
388
389
390
391
392
393
    if (_solver->getUseTurbulenceKappavNuvP0()) {
      nuvFunc = functionSumNew(dofs->nuvDofCenter->getFunction(),_nuv0Func);
      kappavFunc = functionSumNew(dofs->kappavDofCenter->getFunction(),_kapv0Func);
    }
    else {
      nuvFunc = functionSumNew(dofs->nuvDof->getFunction(),_nuv0Func);
      kappavFunc = functionSumNew(dofs->kappavDof->getFunction(),_kapv0Func);
    }
lambrechts's avatar
lambrechts committed
394
395
396
397
398
399
400
401
402
403
404
405
406
  }
  else {
    if ( !nuvFunc && _solver->getSolveUVImplicitVerticalDiffusion() ) {
      Msg::Warning("nuvFunction not defined, using molecular diffusivity");
      nuvFunc = _nuv0Func;
    }
    if ( !kappavFunc && (_solver->getSolveSImplicitVerticalDiffusion() ||
                         _solver->getSolveTImplicitVerticalDiffusion() ||
                         _solver->getSolveSedImplicitVerticalDiffusion() ) ) {
      Msg::Warning("kappavFunction not defined, using molecular diffusivity");
      kappavFunc = _kapv0Func;
    }
  }
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
  if (_solver->getFlagJumpDiffusionTracerVert()){
    if (_solver->_solveS){
      if (_solver->_solveSDiff)
        jumpDiffTracerVertS = new dgJumpDiffusionTracer(*_solver->groups3d, 0, dofs->wDof3d, 0., false, coeffCFLv, jumpPercentS, maxDiffS);
      else 
        jumpDiffTracerVertS = new dgJumpDiffusionTracer(*_solver->groups3d, 0, dofs->wDof3d, approxDT, false, coeffCFLv, jumpPercentS, maxDiffS);
      jumpDiffTracerVertS->apply(dofs->SDof);
      kappavJumpS = jumpDiffTracerVertS->diffusivity().getFunction();
      if (kappavFunc)
        kappavS = functionSumNew(kappavJumpS, kappavFunc);
      else
        kappavS = kappavJumpS;
    }
    if (_solver->_solveT){ 
      if (_solver->_solveTDiff)
        jumpDiffTracerVertT = new dgJumpDiffusionTracer(*_solver->groups3d, 0, dofs->wDof3d, 0., false, coeffCFLv, jumpPercentT, maxDiffT);
      else
        jumpDiffTracerVertT = new dgJumpDiffusionTracer(*_solver->groups3d, 0, dofs->wDof3d, approxDT, false, coeffCFLv, jumpPercentT, maxDiffT);
      jumpDiffTracerVertT->apply(dofs->TDof);
      kappavJumpT = jumpDiffTracerVertT->diffusivity().getFunction();
      if (kappavFunc)
        kappavT = functionSumNew(kappavJumpT, kappavFunc);
      else
        kappavT = kappavJumpT;
    } 
  }
  else if (kappavFunc){
    kappavS = kappavFunc;
    kappavT = kappavFunc;
  }
437
438
439
440
441
442
443
444
445
  // horizontal diffusivity
  // TODO Remove kappahTotal 
  if (kappahTotal){
    kappahS = kappahTotal;
    kappahT = kappahTotal;
  }
  if (!kappahTotal){
    if (_solver->getFlagJumpDiffusionTracer()){
      if (_solver->getSolveS()){
446
        jumpDiffTracerS = new dgJumpDiffusionTracer(*_solver->groups3d, 0, dofs->uvDof, approxDT, true, coeffCFLh, jumpPercentS, maxDiffS); 
447
448
449
450
451
452
453
454
        jumpDiffTracerS->apply(dofs->SDof);
        kappahJumpS = jumpDiffTracerS->diffusivity().getFunction();
        if (kappahUser)
          kappahS = functionSumNew(kappahUser, kappahJumpS);
        else
          kappahS = kappahJumpS;
      }
      if (_solver->getSolveT()){
455
        jumpDiffTracerT = new dgJumpDiffusionTracer(*_solver->groups3d, 0, dofs->uvDof, approxDT, true, coeffCFLh, jumpPercentT, maxDiffT); 
456
457
458
459
460
461
462
463
        jumpDiffTracerT->apply(dofs->TDof);
        kappahJumpT = jumpDiffTracerT->diffusivity().getFunction();
        if (kappahUser)
          kappahT = functionSumNew(kappahUser, kappahJumpT);
        else
          kappahT = kappahJumpT;
      }
    }
464
    else {
465
466
      kappahS = kappahUser;
      kappahT = kappahUser;
467
    }
468
  }
lambrechts's avatar
lambrechts committed
469
470
471
472
473
474
475
476
  // bottom friction
  if ( _solver->getComputeBottomFriction() ) {
    // TODO rename to bottomRoughnessFunc
    if (!z0BFunc)
      Msg::Fatal("bottomRoughnessFunc must be set");
    if (!z0SFunc)
      Msg::Fatal("surfaceRoughnessFunc must be set");
    // TODO rename
477
478
479
    const functor *zBottomFunc2d = dofs->zBotDof2d->getFunction();
    const functor *uvBottomFunc2d = dofs->uvBotDof2d->getFunction();
    const functor *uvIntFunc2d = dofs->uvIntDof2d->getFunction();    
lambrechts's avatar
lambrechts committed
480
481
    // bottom friction velocity in 2d and 3d geometries

delandmeter's avatar
delandmeter committed
482
    zBotDistFunc = new zBotDist(bathFunc2d, zBottomFunc2d);
483
    zBotDist2dPrecompFunc = new functionPrecomputedExtrusion(_solver->extrusion(), _solver->groups2d->getElementGroup(0)->getOrder() * 2 + 1);
delandmeter's avatar
delandmeter committed
484
    if (!bottomFrictionDragCoeff2d ){
lambrechts's avatar
lambrechts committed
485
486
487
488
489
490
491
492
      bottomFrictionDragCoeff2d = new bottomFrictionDragCoeffVonKarman(bathFunc2d, zBottomFunc2d, z0BFunc);
    }
    else if (_solver->getSolveTurbulence())
      Msg::Info("Not sure that GOTM can be used without von Karman bottom friction parametrization\n");
    if (_solver->_solveTurbulence) {
      uvTauBFunc2d = new uvBotStar(bottomFrictionDragCoeff2d, uvBottomFunc2d);
    }
    
493
    bottomFrictionFunc = new dragNormUb(bottomFrictionDragCoeff2d, uvBottomFunc2d);
494
    bottomFriction2dPrecompFunc = new functionPrecomputedExtrusion(_solver->extrusion(), _solver->groups2d->getElementGroup(0)->getOrder() * 2 + 1);
delandmeter's avatar
delandmeter committed
495
    bottomUVDeviation = new bottomVelocityDeviation(bathFunc2d,eta2dCGFunc,uvBottomFunc2d,
lambrechts's avatar
lambrechts committed
496
497
498
                                                    uvIntFunc2d);
  }
  if(_solver->getSolveGMVel()){
499
    eadySpeedFunc = new eadySpeed(dofs->nbDof3d->getFunction(), bathFunc2d, dofs->etaDof2d->getFunction());
lambrechts's avatar
lambrechts committed
500
    GMIndepTermFunc = new GMIndepTerm(dofs->rhoGradDof3d->getFunction());
501
    const functor *streamGradFunc = dofs->streamFuncDof3d->getFunctionGradient();
lambrechts's avatar
lambrechts committed
502
503
504
505
506
507
508
509
510
    GMVelFunc = new GMVelocity(streamGradFunc);
  }
  if(_solver->getSolveIsoDiff()){
    slopeFunc = new slopeFunction(dofs->rhoGradDof3d->getFunction(), dofs->nbDof3d->getFunction());
    slopeTap = new slopeTapering(dofs->rhoGradDof3d->getFunction(), dofs->nbDof3d->getFunction());
    if(!_solver->getSolveTurbulence())
      kappavFunc = new kappavConvFunc(dofs->nbDof3d->getFunction()); 
  }
  if(_solver->getSolveSed()){
511
512
513
514
    if ((!windFunc && windStressFunc) || (windFunc && !windStressFunc))
      Msg::Fatal("Either both windFunc and windStressFunc have to be set, or none of them");
    else if (!windFunc)
      windFunc = zeroFunc2;
515
516
517
518
519
520
521
522
523
524
525
526
    if (wSedFunc.size() == 0){
      wSedFunc.resize(_solver->getNumSedEq());
      for (int iSedEq = 0; iSedEq < _solver->getNumSedEq(); ++iSedEq){
        wSedFunc[iSedEq] = new verticalVelocitySediment(dofs->sedDof[iSedEq]->getFunction(), dofs->wDof3d->getFunction(), dofs->SDof->getFunction(), windFunc);
      }
    }
    if (sedBottomFluxFunc2d.size() == 0) {
      sedBottomFluxFunc2d.resize(_solver->getNumSedEq());
      for (int iSedEq = 0; iSedEq < _solver->getNumSedEq(); ++iSedEq){
        sedBottomFluxFunc2d[iSedEq] = new sedimentBottomFlux(dofs->sedBottomDof2d[iSedEq]->getFunction(), dofs->uvBotDof2d->getFunction(), windFunc, bathFunc2d, dofs->sedGroundDof2d[iSedEq]->getFunction(), dofs->SBottomDof2d->getFunction(), 0.245);
      }
    }
527
    sedBottomFlux2dPrecompFunc = new functionPrecomputedExtrusion(_solver->extrusion(), _solver->groups2d->getElementGroup(0)->getOrder() * 2 + 1);
lambrechts's avatar
lambrechts committed
528
  }
529
  // PrecomputedFunctions
delandmeter's avatar
delandmeter committed
530
  if (_solver->getFlagMovingMesh()){
531
532
    wMeshPrecompFunc = new functionPrecomputed(*_solver->groups3d, _solver->groups3d->getElementGroup(0)->getOrder() * 2 + 1);
    wMeshDzPrecompFunc = new functionPrecomputed(*_solver->groups3d, _solver->groups3d->getElementGroup(0)->getOrder() * 2 + 1);
delandmeter's avatar
delandmeter committed
533
  }
delandmeter's avatar
delandmeter committed
534
535
536
537
538
539
540
541
542
543
544
  
  /*uvCorrectorFunc = new uvCorrector( dofs->uvIntDof2d->getFunction(),
                                     dofs->uvIntTargetDof3d->getFunction(),
                                     eta2dCGFunc,
                                     bathFunc2d );*/
  uvDevFunc = new uvDeviation( dofs->uvDof->getFunction(),
                               uvInt2dFunc,
                               eta2dCGFunc,
                               bathFunc2d );
  uvDevProdFunc = new uvDeviationProducts(uvDevFunc);
  uvDevAdvectionFunc = new uvDevAdvectionTerm(dofs->uvDeviationProdIntDof2d->getFunctionGradient(),eta2dCGFunc,bathFunc2d);
545

lambrechts's avatar
lambrechts committed
546
547
548
549
550
551
552
553
  _initialized= true;
}

void slim3dSolverFunctions::initializeBath()
{
  slim3dSolverDofs *dofs = _solver->getDofs();
  if (!dofs->isInitialized())
    Msg::Fatal("Dofs must be allocated before calling slim3dSolverFunctions::initialize()");
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
554
  if (! bathFunc2d) {
555
556
    dgBathymetryDofGenerator::extractFromMesh(_solver->columnInfo, dofs->bathDof2d);
    dofs->bathDof2d->scatter();
557
    _solver->dgCG2d->dgCGMin(dofs->bathDof2d, dofs->bathCGMinDof2d);
558
    bathFunc2d = new functionPrecomputedExtrusion(_solver->extrusion(), _solver->groups2d->getElementGroup(0)->getOrder() * 2 + 1, *dofs->bathDof2d->getFunction());
559
    bathCGMinFunc2d = new functionPrecomputedExtrusion(_solver->extrusion(), _solver->groups2d->getElementGroup(0)->getOrder() * 2 + 1, *dofs->bathCGMinDof2d->getFunction());
delandmeter's avatar
delandmeter committed
560
  }
delandmeter's avatar
delandmeter committed
561
  // PrecomputedFunctions
562
563
564
565
566
  eta2dFunc = new functionPrecomputedExtrusion(_solver->extrusion(), _solver->groups2d->getElementGroup(0)->getOrder() * 2 + 1);//*dofs->etaDof2d->getFunction());
  uvInt2dFunc = new functionPrecomputedExtrusion(_solver->extrusion(), _solver->groups2d->getElementGroup(0)->getOrder() * 2 + 1);//*dofs->uvIntDof2d->getFunction());
  uvAv2dFunc = new functionPrecomputedExtrusion(_solver->extrusion(), _solver->groups2d->getElementGroup(0)->getOrder() * 2 + 1);//*dofs->uvAvDof2d->getFunction());
  uvIntTarget2dFunc = new functionPrecomputedExtrusion(_solver->extrusion(), _solver->groups2d->getElementGroup(0)->getOrder() * 2 + 1);
  rhoSurf2dFunc = new functionPrecomputedExtrusion(_solver->extrusion(), _solver->groups2d->getElementGroup(0)->getOrder() * 2 + 1);//*dofs->rhoSurfDof2d->getFunction());
567
  rhoSurfGrad2dFunc = new functionPrecomputedExtrusion(_solver->extrusion(), _solver->groups2d->getElementGroup(0)->getOrder() * 2 + 1);//*dofs->rhoSurfDof2d->getFunction());
568
  intDRhodzPrecompFunc = new functionPrecomputed(*_solver->groups3d, _solver->groups3d->getElementGroup(0)->getOrder() * 2 + 1);
569
  //intDRhodzGradPrecompFunc = new functionPrecomputed(*_solver->groups3d, _solver->groups3d->getElementGroup(0)->getOrder() * 2 + 1, 3);
570
571
572
573
  intDRhodzSurf2dPrecompFunc = new functionPrecomputedExtrusion(_solver->extrusion(), _solver->groups2d->getElementGroup(0)->getOrder() * 2 + 1);
  eta2dCGFunc = new functionPrecomputedExtrusion(_solver->extrusion(), _solver->groups2d->getElementGroup(0)->getOrder() * 2 + 1);//*dofs->etaDof2dCG->getFunction());
  eta2dCGFuncOld = new functionPrecomputedExtrusion(_solver->extrusion(), _solver->groups2d->getElementGroup(0)->getOrder() * 2 + 1);//*dofs->etaDof2dCG->getFunction());
  eta2dCGGradFunc = new functionPrecomputedExtrusion(_solver->extrusion(), _solver->groups2d->getElementGroup(0)->getOrder() * 2 + 1);//*dofs->etaDof2dCG->getFunctionGradient());
574
  _solver->verticalBottomRemover->computeHeightAboveVertBottom(bathFunc2d);
lambrechts's avatar
lambrechts committed
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
}

slim3dSolverFunctions::~slim3dSolverFunctions()
{
  if (zeroFunc) delete zeroFunc;
  if (zeroFunc3) delete zeroFunc3;
  if (tinyFunc) delete tinyFunc;
  if (oneFunc) delete oneFunc;
  if (uvCorrectorFunc) delete uvCorrectorFunc;
  if (rhoFunc) delete rhoFunc;
  if (rhoFuncCenter) delete rhoFuncCenter;
  if (eosAlphaFunc) delete eosAlphaFunc;
  if (eosBetaFunc) delete eosBetaFunc;
  if (newZFunc) delete newZFunc;
  if (wMeshSurfFunc) delete wMeshSurfFunc;
  if (wMeshFunc) delete wMeshFunc;
  if (wMeshDzFunc) delete wMeshDzFunc;
  if (_nuv0Func) delete _nuv0Func;
  if (_kapv0Func) delete _kapv0Func;
  if (nuvFunc) delete nuvFunc;
  if (kappavFunc) delete kappavFunc;
  if (uvTauBFunc2d) delete uvTauBFunc2d;
  if (zBotDistFunc) delete zBotDistFunc;
  if (bottomFrictionFunc) delete bottomFrictionFunc;
  if (eadySpeedFunc) delete eadySpeedFunc;
  if (GMIndepTermFunc) delete GMIndepTermFunc;
  if (GMVelFunc) delete GMVelFunc;
  if (slopeFunc) delete slopeFunc;
  if (slopeTap) delete slopeTap;
  if (_bathDof3d) delete _bathDof3d;
605
606
  if (wSedFunc.size() !=0 ) {for (int iSedEq = 0; iSedEq < _solver->getNumSedEq(); ++iSedEq){ delete wSedFunc[iSedEq];} wSedFunc.resize(0);}
  if (sedBottomFluxFunc2d.size() !=0 ) {for (int iSedEq = 0; iSedEq < _solver->getNumSedEq(); ++iSedEq){ delete sedBottomFluxFunc2d[iSedEq];} wSedFunc.resize(0);}
607
608
  if (windStressFunc) delete windStressFunc;
  if (windFunc) delete windFunc;
609
  if (wBottomFunc) delete wBottomFunc;
610
  if (eta2dFunc) delete eta2dFunc;
delandmeter's avatar
delandmeter committed
611
  if (eta2dCGFunc) delete eta2dCGFunc;
612
  if (eta2dCGFuncOld) delete eta2dCGFuncOld;
delandmeter's avatar
delandmeter committed
613
  if (eta2dCGGradFunc) delete eta2dCGGradFunc;
614
615
616
  if (uvInt2dFunc) delete uvInt2dFunc;
  if (uvAv2dFunc) delete uvAv2dFunc;
  if (rhoSurf2dFunc) delete rhoSurf2dFunc;
617
  if (rhoSurfGrad2dFunc) delete rhoSurfGrad2dFunc;
618
  if (atmPressFunc) delete atmPressFunc;
delandmeter's avatar
delandmeter committed
619
  if (wMeshSurf2dPrecompFunc) delete wMeshSurf2dPrecompFunc;
620
621
  if (wMeshPrecompFunc) delete wMeshPrecompFunc;
  if (wMeshDzPrecompFunc) delete wMeshDzPrecompFunc;
622
  if (uvIntTarget2dFunc) delete uvIntTarget2dFunc;
623
624
  if (bottomFriction2dPrecompFunc) delete bottomFriction2dPrecompFunc;
  if (zBotDist2dPrecompFunc) delete zBotDist2dPrecompFunc;
delandmeter's avatar
delandmeter committed
625
  if (sedBottomFlux2dPrecompFunc) delete sedBottomFlux2dPrecompFunc;
626
627
  if (arbitraryMovingMeshFunc) delete arbitraryMovingMeshFunc;
  if (kappaZEqFunc) delete kappaZEqFunc;
628
629
  if (jumpDiffTracerS) delete jumpDiffTracerS;
  if (jumpDiffTracerT) delete jumpDiffTracerT;
630
631
632
633
634
635
636
637
638
639
640
641
642
643
  if (jumpDiffTracerVertS) delete jumpDiffTracerVertS;
  if (jumpDiffTracerVertT) delete jumpDiffTracerVertT;
  if (nuh) delete nuh;
  if (nuh2d) delete nuh2d;
  if (kappavS) delete kappavS;
  if (kappavT) delete kappavT;
  if (kappavJumpS) delete kappavJumpS;
  if (kappavJumpT) delete kappavJumpT;
  if (kappahTotal) delete kappahTotal;
  if (kappahUser) delete kappahUser;
  if (kappahJumpS) delete kappahJumpS;
  if (kappahJumpT) delete kappahJumpT;
  if (kappahS) delete kappahS;
  if (kappahT) delete kappahT;
lambrechts's avatar
lambrechts committed
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
}


slim3dSolverEquations::~slim3dSolverEquations()
{
  delete horMomEq;
  delete vertMomUEq;
  delete wEq;
  delete newRGradEq;
#ifdef HAVE_GOTM
  if ( turbMod ) delete turbMod;
#endif
  delete SEq;
  delete TEq;
  delete tkeAdvEq;
  delete epsAdvEq;
  delete SDiffEq;
  delete TDiffEq;
  delete sedDiffEq;
  delete depthIntegratorRho;
  delete eta2dEq;
  delete uv2dEq;
  delete uvHorDivEq;
667
668
  if ( _solver->getSolveUVGrad())
    delete uvGradEq;
lambrechts's avatar
lambrechts committed
669
670
671
672
673
674
675
  if ( _solver->getSolveGMVel() ) {
    delete rhoGradEq;
    delete nbMod;
    delete rhoGradContProj;
    delete streamFuncContProj;
    delete GMEq;
  }
676
  if ( vJumpEq ) delete vJumpEq;
lambrechts's avatar
lambrechts committed
677
678
679
680
681
682
683
684
}

slim3dSolverEquations::slim3dSolverEquations(slim3dSolver* solver) : _solver(solver)
{
  // create all equations with no arguments
  horMomEq = new dgConservationLawSW3dMomentum();
  vertMomUEq = new dgConservationLawSW3dMomWx();
  wEq = new dgConservationLawSW3dContinuity(0);
685
  newRGradEq = new dgConservationLawSW3dPressure(INTEGRATE_FROM_TOP);
lambrechts's avatar
lambrechts committed
686
687
688
689
690
691
692
693
694
695
696
697
698
699
#ifdef HAVE_GOTM
  turbMod = NULL;
#endif
  SEq = new dgConservationLawSW3dTracer();
  TEq = new dgConservationLawSW3dTracer();
  tkeAdvEq = new dgConservationLawSW3dTracer();
  epsAdvEq = new dgConservationLawSW3dTracer();
  SDiffEq = new dgConservationLawSW3dTracerVDiff();
  TDiffEq = new dgConservationLawSW3dTracerVDiff();
  sedDiffEq = new dgConservationLawSW3dTracerVDiff();
  depthIntegratorRho = new dgSW3dDepthIntegrate(1,INTEGRATE_FROM_TOP);
  eta2dEq = new dgConservationLawSW2dEta();
  uv2dEq = new dgConservationLawSW2dU();
  uvHorDivEq = new dgUVHorDivergence();
delandmeter's avatar
delandmeter committed
700
  uvGradEq = NULL;
delandmeter's avatar
delandmeter committed
701
  //wMeshGradEq = NULL;
702
703
704
705
706
707
708
  if (_solver->getSolveGMVel()){
    rhoGradEq = NULL;
    rhoGradContProj = NULL;
    streamFuncContProj = NULL;
    nbMod = NULL;
    GMEq = NULL;
  }
709
  vJumpEq = NULL;
lambrechts's avatar
lambrechts committed
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
  _initialized = false;
}

void slim3dSolverEquations::initialize()
{
  slim3dSolverDofs* dofs = _solver->getDofs();
  if (!dofs->isInitialized())
    Msg::Fatal("Dofs must be allocated before calling slim3dSolverEquations::initialize()");
  slim3dSolverFunctions* funcs = _solver->functions;
  if (!funcs->isInitialized())
    Msg::Fatal("Functions must be allocated before calling slim3dSolverEquations::initialize()");
  if ( _solver->getSolveTurbulence() && !_solver->getComputeBottomFriction() )
    Msg::Fatal("Turbulence closure model requires bottom friction");
  // turbulence
#ifdef HAVE_GOTM
  if ( _solver->_solveTurbulence ) {
    turbMod = new dgSW3dTurbulenceGOTM(_solver->columnInfo,_solver->turbulenceSetupFile);
    // TODO check which of these are not actually needed!
    turbMod->setZFunction(dofs->zDofCenter->getFunction());
    turbMod->setUVFunction(dofs->uvDofCenter->getFunction());
    turbMod->setUVGradientFunction(dofs->uvDofCenter->getFunctionGradient());
    turbMod->setEOSAlphaFunction(funcs->eosAlphaFunc);
    turbMod->setEOSBetaFunction(funcs->eosBetaFunc);
    turbMod->setSGradientFunction(funcs->SCenterGradientFunc);
    turbMod->setTGradientFunction(funcs->TCenterGradientFunc);
    turbMod->setZ0SFunction(funcs->z0SFunc);
    turbMod->setZ0BFunction(funcs->z0BFunc);
    if ( funcs->windStressFunc)
      turbMod->setWindStressFunction(dofs->windStressDof2dCenter->getFunction());
    else
      turbMod->setWindStressFunction(funcs->zeroFunc3);
    turbMod->setUVStarBFunction(dofs->uvTauBDof2dCenter->getFunction());
    turbMod->setNuvContainer(dofs->nuvDofCenter);
    turbMod->setKappavContainer(dofs->kappavDofCenter);
    turbMod->setTkeContainer(dofs->tkeDofCenter);
    turbMod->setEpsContainer(dofs->epsDofCenter);
    turbMod->setLContainer(dofs->lDof);
    turbMod->setSSContainer(dofs->ssDof);
    turbMod->setNNContainer(dofs->nnDof);
    turbMod->setSContainer(dofs->SDofCenter);
    turbMod->setTContainer(dofs->TDofCenter);
  }
#else
  if ( _solver->_solveTurbulence )
    Msg::Fatal("GOTM not compiled");
#endif
    
  if ( _solver->getSolveGMVel() ){
    nbMod = new dgSW3dVerticalModel(_solver->columnInfo);
    nbMod->setComputeNNWithEOS(true);
760
    nbMod->setZFunction(dofs->zDof->getFunction());
lambrechts's avatar
lambrechts committed
761
762
763
764
    nbMod->setSContainer(dofs->SDof);
    nbMod->setTContainer(dofs->TDof);
    
  }
765
766
  //const functor *bath3d = funcs->bathFunc3d;
  const functor *bath2d = funcs->bathFunc2d;
767
  const functor *bathCGMin2d = funcs->bathCGMinFunc2d;
768
769
770
771
772
773
774
775
776
777
778
779
780
  //const functor *eta2d = dofs->etaDof2d->getFunction();
  const functor *eta2d = funcs->eta2dFunc;
  const functor *uv3d = dofs->uvDof->getFunction();
  //const functor *uvInt2d = dofs->uvIntDof2d->getFunction();
  const functor *uvInt2d = funcs->uvInt2dFunc;
  //const functor *uvAv2d = dofs->uvAvDof2d->getFunction();
  const functor *uvAv2d = funcs->uvAv2dFunc;
  const functor *w3d = dofs->wDof3d->getFunction();
  //const functor *w3dCopy = dofs->wDof3dCopy->getFunction();
  //const functor *wMesh = funcs->wMeshFunc;
  //const functor *wMeshDz = funcs->wMeshDzFunc;
  const functor *wMesh = funcs->zeroFunc;
  const functor *wMeshDz = funcs->zeroFunc;
delandmeter's avatar
delandmeter committed
781
782
783
784
  if (_solver->getFlagMovingMesh()){
    wMesh = funcs->wMeshPrecompFunc;
    wMeshDz = funcs->wMeshDzPrecompFunc;
  }
785
786
787
788
789
  const functor *rGrad = dofs->rGradDof3d->getFunction();
  //const functor *r = dofs->rDof3d->getFunction();
  const functor *rGradInt2d = dofs->rGradIntDof2d->getFunction();
  //const functor *rhoSurf2d = dofs->rhoSurfDof2d->getFunction();
  const functor *rhoSurf2d = funcs->rhoSurf2dFunc;
790
  const functor *rhoSurfGrad2d = funcs->rhoSurfGrad2dFunc;
791
792
793
  //const functor *rBot2d = dofs->rBotDof2d->getFunction();
  const functor* bottomFrictionDragCoeff2d = funcs->bottomFrictionDragCoeff2d;
  const functor* bottomUVDeviation = funcs->bottomUVDeviation;
794
  horMomEq->setBathymetry(bath2d);
795
  horMomEq->setBathymetryMinCG(bathCGMin2d);
796
797
  horMomEq->setEta(eta2d);
  horMomEq->setUV2d(uvInt2d);
798
799
800
801
802
803
804
805
806
807
  if (_solver->_solveUVAdvVert){
    horMomEq->setW(funcs->zeroFunc);
    horMomEq->setWMesh(funcs->zeroFunc);
    horMomEq->setWMeshDz(funcs->zeroFunc);
  }
  else{
    horMomEq->setW(w3d);
    horMomEq->setWMesh(wMesh);
    horMomEq->setWMeshDz(wMeshDz);
  }
lambrechts's avatar
lambrechts committed
808
  horMomEq->setRGrad(rGrad);
karna's avatar
karna committed
809
  horMomEq->setRhoSurf(rhoSurf2d);
810
  horMomEq->setRhoSurfGrad(rhoSurfGrad2d);
811
812
  if (funcs->atmPressFunc)
    horMomEq->setAtmPress(funcs->atmPressFunc);
lambrechts's avatar
lambrechts committed
813
814
  if (funcs->coriolisFunc)
    horMomEq->setCoriolisFactor(funcs->coriolisFunc);
815
816
  if (funcs->nuh)
    horMomEq->setNuH(funcs->nuh);
lambrechts's avatar
lambrechts committed
817
818
819
820
821
822
823
  if (!_solver->_solveUVDiff && funcs->nuvFunc) {
    Msg::Warning("Vertical diffusion equation of UV disabled, assigning nuv to horizontal momentum equation");
    horMomEq->setNuV(funcs->nuvFunc);
  }
  if ( _solver->_solveUVDiff ) {
    vertMomUEq->setNuV(funcs->nuvFunc);
  }
824
825
826
827
828
829
830
  if (_solver->_solveUVAdvVert){
    vertMomUEq->setW(w3d);
    vertMomUEq->setWMesh(wMesh);
    vertMomUEq->setWMeshDz(wMeshDz);
    if (_solver->getUseConservativeALE())
      vertMomUEq->setUseConservativeALE(true);
  }
lambrechts's avatar
lambrechts committed
831
832

  wEq->setUV(uv3d);
833
  wEq->setEta(eta2d);
834
  wEq->setBathymetry(bath2d);
835
  wEq->setBathymetryMinCG(bathCGMin2d);
lambrechts's avatar
lambrechts committed
836
  //wEq->setW(w3dCopy);
837
838
839
  //wMeshGradEq = new dgFEGradient(dofs->wMeshDof3d->getFunction());
  //dgBoundaryCondition* wMeshGradSymmBnd = wMeshGradEq->newSymmetryBoundary();
  //wMeshGradEq->addBoundaryCondition(_solver->physicalTags2d,wMeshGradSymmBnd);
lambrechts's avatar
lambrechts committed
840
 
841
842
843
844
845
  //newRGradEq->setRhoFunction(funcs->rhoFunc);
  newRGradEq->setRhoFunction(dofs->rhoDof3d->getFunction());
  newRGradEq->setRhoGradFunction(dofs->rhoDof3d->getFunctionGradient());
  newRGradEq->setUseRGradByPart(_solver->getUseRGradByPart());

846
  //newRGradEq->setBoundarySymmetry(_solver->physicalTags2d);
lambrechts's avatar
lambrechts committed
847

848

lambrechts's avatar
lambrechts committed
849
850
  if ( _solver->_solveS ) {
    SEq->setUV(uv3d);
851
    SEq->setEta(eta2d);
852
853
    //SEq->setUV(funcs->zeroFunc2);
    //SEq->setEta(funcs->zeroFunc);
854
    SEq->setBathymetryMinCG(bathCGMin2d);
855
856
857
858
859
860
861
    if ( _solver->_solveSAdvVert) {
      SEq->setW(funcs->zeroFunc);
      SEq->setWMesh(funcs->zeroFunc);
      SEq->setWMeshDz(funcs->zeroFunc);
    }
    else{
      SEq->setW(w3d);
862
      //SEq->setW(funcs->zeroFunc);
863
864
865
      SEq->setWMesh(wMesh);
      SEq->setWMeshDz(wMeshDz);
    }
866
867
    if (funcs->kappahS)
      SEq->setKappaH(funcs->kappahS);
868
    if (!_solver->_solveSDiff && funcs->kappavS) {
lambrechts's avatar
lambrechts committed
869
      Msg::Warning("Vertical diffusion equation of S disabled, assigning kappav to S equation");
870
      SEq->setKappaV(funcs->kappavS);
lambrechts's avatar
lambrechts committed
871
872
873
    }
    if (_solver->getSolveGMVel()) {
      SEq->setGMVel(dofs->GMVelDof3d->getFunction());
874
      SEq->setKappaV(funcs->kappavS);
lambrechts's avatar
lambrechts committed
875
876
877
878
879
880
881
882
883
    }
    if (_solver->getSolveIsoDiff()){
      if (!_solver->getSolveGMVel())
	Msg::Fatal("to Solve Iso Diff, solving GMVel is required!");
      SEq->setIsoDiff(funcs->slopeFunc, funcs->slopeTap);
    }
  }

  if ( _solver->_solveSDiff ) {
884
    SDiffEq->setKappaV(funcs->kappavS);
885
    SDiffEq->setBoundary0Flux(_solver->physicalTags2d);
lambrechts's avatar
lambrechts committed
886
  }
887
888
889
890
891
892
893
  if (_solver->_solveSAdvVert ){
    SDiffEq->setW(w3d);
    SDiffEq->setWMesh(wMesh);
    SDiffEq->setWMeshDz(wMeshDz);
    if (_solver->getUseConservativeALE())
      SDiffEq->setUseConservativeALE(true);
  }
lambrechts's avatar
lambrechts committed
894
895
896

  if ( _solver->_solveT ) {
    TEq->setUV(uv3d);
897
    TEq->setEta(eta2d);
898
    TEq->setBathymetryMinCG(bathCGMin2d);
899
900
901
902
903
904
905
906
907
908
    if ( _solver->_solveTAdvVert) {
      TEq->setW(funcs->zeroFunc);
      TEq->setWMesh(funcs->zeroFunc);
      TEq->setWMeshDz(funcs->zeroFunc);
    }
    else{
      TEq->setW(w3d);
      TEq->setWMesh(wMesh);
      TEq->setWMeshDz(wMeshDz);
    }
909
910
    if (funcs->kappahT)
      TEq->setKappaH(funcs->kappahT);
911
    if (!_solver->_solveTDiff && funcs->kappavT) {
lambrechts's avatar
lambrechts committed
912
      Msg::Warning("Vertical diffusion equation of T disabled, assigning kappav to T equation");
913
      TEq->setKappaV(funcs->kappavT);
lambrechts's avatar
lambrechts committed
914
915
    }
    if (_solver->getSolveGMVel()){
delandmeter's avatar
delandmeter committed
916
      rhoGradEq = new dgFEGradient(funcs->rhoFunc);
917
      rhoGradEq->setBoundarySymmetry(_solver->physicalTags2d);
lambrechts's avatar
lambrechts committed
918
      TEq->setGMVel(dofs->GMVelDof3d->getFunction());
919
      TEq->setKappaV(funcs->kappavT);
lambrechts's avatar
lambrechts committed
920
921
922
    }
    if (_solver->getSolveIsoDiff()){
      if (!_solver->getSolveGMVel())
923
      	Msg::Fatal("to Solve Iso Diff, solving GMVel is required!");
lambrechts's avatar
lambrechts committed
924
925
926
927
928
      TEq->setIsoDiff(funcs->slopeFunc, funcs->slopeTap);
    }
  }

  if ( _solver->_solveTDiff ) {
929
    TDiffEq->setKappaV(funcs->kappavT);
930
    TDiffEq->setBoundary0Flux(_solver->physicalTags2d);
lambrechts's avatar
lambrechts committed
931
  }
932
933
934
935
936
937
938
  if (_solver->_solveTAdvVert ){
    TDiffEq->setW(w3d);
    TDiffEq->setWMesh(wMesh);
    TDiffEq->setWMeshDz(wMeshDz);
    if (_solver->getUseConservativeALE())
      TDiffEq->setUseConservativeALE(true);
  }
lambrechts's avatar
lambrechts committed
939
940
  
  if ( _solver->_solveSed ){
941
942
    sedEq.resize(_solver->getNumSedEq());
    for ( int iSedEq = 0; iSedEq < _solver->getNumSedEq(); ++iSedEq) {
943
944
      if (iSedEq == 1 && _solver->_solveSedAdvVert)
        Msg::Fatal("For implicit vertical advection, various sedDiffEq must be implemented. It is not done yet.");
945
946
947
948
      const functor* wSed = funcs->wSedFunc[iSedEq];
      sedEq[iSedEq] = new dgConservationLawSW3dTracer();
      sedEq[iSedEq]->setUV(uv3d);
      sedEq[iSedEq]->setEta(eta2d);
949
      sedEq[iSedEq]->setBathymetryMinCG(bathCGMin2d);
950
951
952
953
954
955
956
957
958
959
      if ( _solver->_solveSedAdvVert) {
        sedEq[iSedEq]->setW(funcs->zeroFunc);
        sedEq[iSedEq]->setWMesh(funcs->zeroFunc);
        sedEq[iSedEq]->setWMeshDz(funcs->zeroFunc);
      }
      else{
        sedEq[iSedEq]->setW(wSed);
        sedEq[iSedEq]->setWMesh(wMesh);
        sedEq[iSedEq]->setWMeshDz(wMeshDz);
      }
960
961
962
963
964
965
966
967
968
      if (funcs->kappahTotal)
        sedEq[iSedEq]->setKappaH(funcs->kappahTotal);
      if (!_solver->_solveSedDiff && funcs->kappavFunc) {
        Msg::Warning("Vertical diffusion equation of SedEq disabled, assigning kappav to Sediment equation");
        sedEq[iSedEq]->setKappaV(funcs->kappavFunc);
      }
      //dgBoundaryCondition* sedBottomFluxBnd = sedEq[iSedEq]->newInFluxBoundary(dofs->sedBottomFluxDof2d->getFunction()); // must be the dof and not directly the function. To not be updated between sedBottom Equation and sed Equation
      dgBoundaryCondition* sedBottomFluxBnd = sedEq[iSedEq]->newNeumannBoundary(funcs->sedBottomFlux2dPrecompFunc); // must be the dof and not directly the function. To not be updated between sedBottom Equation and sed Equation
      sedEq[iSedEq]->addBoundaryCondition(_solver->bottomTags , sedBottomFluxBnd);
lambrechts's avatar
lambrechts committed
969
970
971
972
973
    }
  }
  
  if ( _solver->_solveSedDiff ) {
    sedDiffEq->setKappaV(funcs->kappavFunc);
974
    sedDiffEq->setBoundary0Flux(_solver->physicalTags2d);
lambrechts's avatar
lambrechts committed
975
  }
976
977
978
979
980
981
982
983
984
  if (_solver->_solveSedAdvVert ){
    const functor* wSed = funcs->wSedFunc[0];
    sedDiffEq->setW(wSed);
    sedDiffEq->setWMesh(wMesh);
    sedDiffEq->setWMeshDz(wMeshDz);
    if (_solver->getUseConservativeALE())
      sedDiffEq->setUseConservativeALE(true);
  }

lambrechts's avatar
lambrechts committed
985
986
987
  if ( _solver->_advectTKE ) {
    tkeAdvEq->setUV(uv3d);
    tkeAdvEq->setW(w3d);
988
    tkeAdvEq->setEta(eta2d);
989
    tkeAdvEq->setBathymetryMinCG(bathCGMin2d);
lambrechts's avatar
lambrechts committed
990
991
992
993
    tkeAdvEq->setWMesh(wMesh);
    tkeAdvEq->setWMeshDz(wMeshDz);
    epsAdvEq->setUV(uv3d);
    epsAdvEq->setW(w3d);
994
    epsAdvEq->setEta(eta2d);
995
    epsAdvEq->setBathymetryMinCG(bathCGMin2d);
lambrechts's avatar
lambrechts committed
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
    epsAdvEq->setWMesh(wMesh);
    epsAdvEq->setWMeshDz(wMeshDz);
  }

  eta2dEq->setBathymetry(bath2d);
  if ( _solver->_useModeSplit )
    eta2dEq->setUV(uvAv2d);
  else {
    eta2dEq->setUsesDepthIntegratedUV(true);
    eta2dEq->setUV(uvInt2d);
  }
  uv2dEq->setBathymetry(bath2d);
  uv2dEq->setEta(eta2d);
  //uv2dEq->setRInt(rInt2d);
  //uv2dEq->setRInt(funcs->zeroFunc);
  //uv2dEq->setRBot(rBot2d);
  uv2dEq->setRGradInt(rGradInt2d);
  uv2dEq->setRhoSurf(rhoSurf2d);
1014
  uv2dEq->setRhoSurfGrad(rhoSurfGrad2d);
1015
1016
  if (funcs->atmPressFunc)
    uv2dEq->setAtmPress(funcs->atmPressFunc);
lambrechts's avatar
lambrechts committed
1017
1018
1019
1020
1021
1022
1023
1024
1025
  if (funcs->windStressFunc)
    uv2dEq->setWindStress(funcs->windStressFunc);
  
  if ( _solver->getComputeBottomFriction() ) {
    uv2dEq->setBottomFrictionCoefficient(bottomFrictionDragCoeff2d);
    uv2dEq->setBottomUV(bottomUVDeviation);
  }
  if (funcs->coriolisFunc)
    uv2dEq->setCoriolisFactor(funcs->coriolisFunc);
1026
1027
  if (funcs->nuh2d)
    uv2dEq->setNuH(funcs->nuh2d);
lambrechts's avatar
lambrechts committed
1028
  if (funcs->uvDevAdvectionFunc)
1029
    uv2dEq->setUVDevAdvection(funcs->uvDevAdvectionFunc);
1030
1031
1032
1033
1034
1035
1036

  bool linear = false;
  eta2dEq->setIsLinear(linear);
  uv2dEq->setIsLinear(linear);
  bool lf = false;
  eta2dEq->setLaxFriedrichs(lf);
  uv2dEq->setLaxFriedrichs(lf);
lambrechts's avatar
lambrechts committed
1037
  
1038
  uvHorDivEq->setBathymetryMinCG(bathCGMin2d);
lambrechts's avatar
lambrechts committed
1039
  uvHorDivEq->setUV(uv3d);
1040
  uvHorDivEq->setEta(eta2d);
1041
  uvHorDivEq->setBoundarySymmetry(_solver->physicalTags2d);
1042
1043

  if (_solver->getSolveUVGrad()){
delandmeter's avatar
delandmeter committed
1044
    uvGradEq = new dgFEVectorGradient(uv3d);
1045
    uvGradEq->setBathymetryMinCG(bathCGMin2d);
1046
    uvGradEq->setEta(eta2d);
1047
    uvGradEq->setBoundarySymmetry(_solver->physicalTags2d);
1048
  }
lambrechts's avatar
lambrechts committed
1049
1050
1051
1052
  
  depthIntegratorRho->setFunction(funcs->rhoFunc);
  for (size_t i = 0; i < _solver->physicalTags2d.size(); i++ ) {
    if ( _solver->physicalTags2d[i] != "top" )
1053
      depthIntegratorRho->setBoundary0Flux(_solver->physicalTags2d[i]);
lambrechts's avatar
lambrechts committed
1054
    else
1055
      depthIntegratorRho->setBoundarySymmetry(_solver->physicalTags2d[i], "");
lambrechts's avatar
lambrechts committed
1056
1057
1058
1059
  }
  
  if (_solver->getSolveGMVel()){
    GMEq = new dgEddyTransportFlux(funcs->eadySpeedFunc, funcs->GMIndepTermFunc, dofs->nbDof3d->getFunction());
1060
    GMEq->setBoundary0Flux(_solver->physicalTags2d);
lambrechts's avatar
lambrechts committed
1061
    for (size_t i = 0; i < _solver->physicalTags2d.size(); ++i)
1062
      GMEq->addStrongBoundaryCondition(2,_solver->physicalTags2d[i],funcs->zeroFunc3);
lambrechts's avatar
lambrechts committed
1063
  }
1064
1065
1066
1067
1068

  if ( _solver->getUseAdaptiveVerticalGrid() ) {
    vJumpEq =  new dgConservationLawSW3dTracer();
    vJumpEq->setUV(uv3d);
    vJumpEq->setEta(eta2d);
1069
    vJumpEq->setBathymetryMinCG(bathCGMin2d);
1070
1071
1072
1073
1074
1075
1076
1077
1078
    vJumpEq->setW(w3d);
    //vJumpEq->setW(funcs->zeroFunc);
    vJumpEq->setWMesh(wMesh);
    vJumpEq->setWMeshDz(wMeshDz);
    dgBoundaryCondition* vJumpSymmBnd = vJumpEq->newSymmetryBoundary();
    vJumpEq->addBoundaryCondition(_solver->physicalTags2d, vJumpSymmBnd);
  }


lambrechts's avatar
lambrechts committed
1079
1080
  _initialized = true;
}