proteus 1.9.0
C/C++/Fortran libraries
Loading...
Searching...
No Matches
SW2DCV.h
Go to the documentation of this file.
1#ifndef SW2DCV_H
2#define SW2DCV_H
3#include "ArgumentsDict.h"
4#include "CompKernel.h"
5#include "ModelFactory.h"
6#include "xtensor-python/pyarray.hpp"
7#include <assert.h>
8#include <cmath>
9#include <iostream>
10#include <valarray>
11
12namespace py = pybind11;
13
14#define POWER_SMOOTHNESS_INDICATOR 2
15#define VEL_FIX_POWER 2.
16#define REESTIMATE_MAX_EDGE_BASED_CFL 0
17#define LIMITING_ITERATION 2
18
19/* inline functions */
20namespace proteus {
21// FOR CELL BASED ENTROPY VISCOSITY
22inline double ENTROPY(const double &g, const double &h, const double &hu,
23 const double &hv, const double &z,
24 const double &one_over_hReg) {
25 return 0.5 *
26 (g * h * h + one_over_hReg * (hu * hu + hv * hv) + 2. * g * h * z);
27}
28inline double DENTROPY_DH(const double &g, const double &h, const double &hu,
29 const double &hv, const double &z,
30 const double &one_over_hReg) {
31 return g * h - 0.5 * (hu * hu + hv * hv) * std::pow(one_over_hReg, 2) + g * z;
32}
33inline double DENTROPY_DHU(const double &g, const double &h, const double &hu,
34 const double &hv, const double &z,
35 const double &one_over_hReg) {
36 return hu * one_over_hReg;
37}
38inline double DENTROPY_DHV(const double &g, const double &h, const double &hu,
39 const double &hv, const double &z,
40 const double &one_over_hReg) {
41 return hv * one_over_hReg;
42}
43inline double ENTROPY_FLUX1(const double &g, const double &h, const double &hu,
44 const double &hv, const double &z,
45 const double &one_over_hReg) {
46 return (ENTROPY(g, h, hu, hv, z, one_over_hReg) + 0.5 * g * h * h +
47 g * h * z) *
48 hu * one_over_hReg;
49}
50inline double ENTROPY_FLUX2(const double &g, const double &h, const double &hu,
51 const double &hv, const double &z,
52 const double &one_over_hReg) {
53 return (ENTROPY(g, h, hu, hv, z, one_over_hReg) + 0.5 * g * h * h +
54 g * h * z) *
55 hv * one_over_hReg;
56}
57// FOR ESTIMATING MAX WAVE SPEEDS
58inline double f(const double &g, const double &h, const double &hZ) {
59 return ((h <= hZ) ? 2. * (sqrt(g * h) - sqrt(g * hZ))
60 : (h - hZ) * sqrt(0.5 * g * (h + hZ) / h / hZ));
61}
62inline double phi(const double &g, const double &h, const double &hL,
63 const double &hR, const double &uL, const double &uR) {
64 return (f(g, h, hL) + f(g, h, hR) + uR - uL);
65}
66inline double fp(const double &g, const double &h, const double &hZ) {
67 return ((h <= hZ)
68 ? sqrt(g / h)
69 : g * (2 * h * h + h * hZ + hZ * hZ) /
70 (2 * sqrt(2 * g) * h * h * hZ * sqrt(1 / h + 1 / hZ)));
71}
72inline double phip(const double &g, const double &h, const double &hL,
73 const double &hR) {
74 return (fp(g, h, hL) + fp(g, h, hR));
75}
76inline double nu1(const double &g, const double &hStar, const double &hL,
77 const double &uL, const double &one_over_hL) {
78 return (uL - sqrt(g * hL) *
79 sqrt((1. + fmax((hStar - hL) / 2. * one_over_hL, 0.0)) *
80 (1. + fmax((hStar - hL) * one_over_hL, 0.))));
81}
82inline double nu3(const double &g, const double &hStar, const double &hR,
83 const double &uR, const double &one_over_hR) {
84 return (uR + sqrt(g * hR) *
85 sqrt((1. + fmax((hStar - hR) / 2. * one_over_hR, 0.0)) *
86 (1. + fmax((hStar - hR) * one_over_hR, 0.))));
87}
88inline double phiDiff(const double &g, const double &h1k, const double &h2k,
89 const double &hL, const double &hR, const double &uL,
90 const double &uR) {
91 return ((phi(g, h2k, hL, hR, uL, uR) - phi(g, h1k, hL, hR, uL, uR)) /
92 (h2k - h1k));
93}
94inline double phiDDiff1(const double &g, const double &h1k, const double &h2k,
95 const double &hL, const double &hR, const double &uL,
96 const double &uR) {
97 return ((phiDiff(g, h1k, h2k, hL, hR, uL, uR) - phip(g, h1k, hL, hR)) /
98 (h2k - h1k));
99}
100inline double phiDDiff2(const double &g, const double &h1k, const double &h2k,
101 const double &hL, const double &hR, const double &uL,
102 const double &uR) {
103 return ((phip(g, h2k, hL, hR) - phiDiff(g, h1k, h2k, hL, hR, uL, uR)) /
104 (h2k - h1k));
105}
106inline double hStarLFromQuadPhiFromAbove(const double &g, const double &hStarL,
107 const double &hStarR, const double &hL,
108 const double &hR, const double &uL,
109 const double &uR) {
110 return (hStarL -
111 2 * phi(g, hStarL, hL, hR, uL, uR) /
112 (phip(g, hStarL, hL, hR) +
113 sqrt(std::pow(phip(g, hStarL, hL, hR), 2) -
114 4 * phi(g, hStarL, hL, hR, uL, uR) *
115 phiDDiff1(g, hStarL, hStarR, hL, hR, uL, uR))));
116}
117inline double hStarRFromQuadPhiFromBelow(const double &g, const double &hStarL,
118 const double &hStarR, const double &hL,
119 const double &hR, const double &uL,
120 const double &uR) {
121 return (hStarR -
122 2 * phi(g, hStarR, hL, hR, uL, uR) /
123 (phip(g, hStarR, hL, hR) +
124 sqrt(std::pow(phip(g, hStarR, hL, hR), 2) -
125 4 * phi(g, hStarR, hL, hR, uL, uR) *
126 phiDDiff2(g, hStarL, hStarR, hL, hR, uL, uR))));
127}
128} // namespace proteus
129
130namespace proteus {
131
132class SW2DCV_base {
133public:
134 virtual ~SW2DCV_base() {}
135 virtual void convexLimiting(arguments_dict &args) = 0;
136 virtual double calculateEdgeBasedCFL(arguments_dict &args) = 0;
137 virtual void calculateEV(arguments_dict &args) = 0;
138 virtual void calculateResidual(arguments_dict &args) = 0;
139 virtual void calculateMassMatrix(arguments_dict &args) = 0;
140 virtual void calculateLumpedMassMatrix(arguments_dict &args) = 0;
141};
142
143template <class CompKernelType, int nSpace, int nQuadraturePoints_element,
144 int nDOF_mesh_trial_element, int nDOF_trial_element,
145 int nDOF_test_element, int nQuadraturePoints_elementBoundary>
146class SW2DCV : public SW2DCV_base {
147public:
148 const int nDOF_test_X_trial_element;
149 CompKernelType ck;
150 SW2DCV()
151 : nDOF_test_X_trial_element(nDOF_test_element * nDOF_trial_element),
152 ck() {
153 std::cout << "Constructing SW2DCV<CompKernelTemplate<" << nSpace << ","
154 << nQuadraturePoints_element << "," << nDOF_mesh_trial_element
155 << "," << nDOF_trial_element << "," << nDOF_test_element << ","
156 << nQuadraturePoints_elementBoundary << ">());" << std::endl
157 << std::flush;
158 }
159
160 inline double maxWaveSpeedSharpInitialGuess(double g, double nx, double ny,
161 double hL, double huL, double hvL,
162 double hR, double huR, double hvR,
163 double hEps, bool debugging) {
164 double lambda1, lambda3;
165 // 1-eigenvalue: uL-sqrt(g*hL)
166 // 3-eigenvalue: uR+sqrt(g*hR)
167
168 // To avoid division by 0
169 double one_over_hL = 2.0 * hL / (hL * hL + std::pow(fmax(hL, hEps), 2.0));
170 double one_over_hR = 2.0 * hR / (hR * hR + std::pow(fmax(hR, hEps), 2.0));
171
172 double hVelL = nx * huL + ny * hvL;
173 double hVelR = nx * huR + ny * hvR;
174 double velL = one_over_hL * hVelL;
175 double velR = one_over_hR * hVelR;
176
177 double x0 = std::pow(2. * sqrt(2.) - 1., 2.);
178 double hMin = fmin(hL, hR);
179 double hMax = fmax(hL, hR);
180
181 double hStar;
182 double fMin = phi(g, x0 * hMin, hL, hR, velL, velR);
183 double fMax = phi(g, x0 * hMax, hL, hR, velL, velR);
184
185 double sqrMin = sqrt(hMin);
186 double sqrMax = sqrt(hMax);
187
188 if (0. <= fMin) {
189 hStar = fmin(
190 x0 * hMin,
191 std::pow(fmax(0., velL - velR + 2. * sqrt(g) * (sqrt(hL) + sqrt(hR))),
192 2) /
193 16. / g);
194 } else if (0. <= fMax) {
195 double a = 1.0 / (2.0 * sqrt(2.0));
196 double c = -hMin * a - sqrMin * sqrMax +
197 sqrMin * (velR - velL) / (2.0 * sqrt(g));
198 double delta = hMin - 4.0 * a * c;
199 if (delta < 0.0) {
200 std::cout << "Bug in computing lambda. Exiting." << std::endl;
201 abort();
202 }
203 hStar = fmin(x0 * hMax, std::pow((-sqrMin + sqrt(delta)) / (2. * a), 2));
204 } else {
205 hStar = sqrMin * sqrMax *
206 (1.0 + sqrt(2.0 / g) * (velL - velR) / (sqrMin + sqrMax));
207 }
208
209 // return lambda_max
210 lambda1 = nu1(g, hStar, hL, velL, one_over_hL);
211 lambda3 = nu3(g, hStar, hR, velR, one_over_hR);
212 return fmax(fabs(lambda1), fabs(lambda3));
213 }
214
215 inline void calculateCFL(const double &elementDiameter, const double &g,
216 const double &h, const double &hu, const double &hv,
217 const double hEps, double &cfl) {
218 double cflx, cfly, c = sqrt(fmax(g * hEps, g * h));
219 double u = 2 * h / (h * h + std::pow(fmax(h, hEps), 2)) * hu;
220 double v = 2 * h / (h * h + std::pow(fmax(h, hEps), 2)) * hv;
221
222 if (u > 0.0)
223 cflx = (u + c) / elementDiameter;
224 else
225 cflx = fabs(u - c) / elementDiameter;
226
227 if (v > 0.0)
228 cfly = (v + c) / elementDiameter;
229 else
230 cfly = fabs(v - c) / elementDiameter;
231 cfl = sqrt(cflx * cflx + cfly * cfly); // hack, conservative estimate
232 }
233
234 void convexLimiting(arguments_dict &args) {
235 double dt = args.scalar<double>("dt");
236 int NNZ = args.scalar<int>("NNZ");
237 int numDOFs = args.scalar<int>("numDOFs");
238 xt::pyarray<double> &lumped_mass_matrix =
239 args.array<double>("lumped_mass_matrix");
240 xt::pyarray<double> &h_old = args.array<double>("h_old");
241 xt::pyarray<double> &hu_old = args.array<double>("hu_old");
242 xt::pyarray<double> &hv_old = args.array<double>("hv_old");
243 xt::pyarray<double> &b_dof = args.array<double>("b_dof");
244 xt::pyarray<double> &high_order_hnp1 =
245 args.array<double>("high_order_hnp1");
246 xt::pyarray<double> &high_order_hunp1 =
247 args.array<double>("high_order_hunp1");
248 xt::pyarray<double> &high_order_hvnp1 =
249 args.array<double>("high_order_hvnp1");
250 xt::pyarray<double> &extendedSourceTerm_hu =
251 args.array<double>("extendedSourceTerm_hu");
252 xt::pyarray<double> &extendedSourceTerm_hv =
253 args.array<double>("extendedSourceTerm_hv");
254 xt::pyarray<double> &limited_hnp1 = args.array<double>("limited_hnp1");
255 xt::pyarray<double> &limited_hunp1 = args.array<double>("limited_hunp1");
256 xt::pyarray<double> &limited_hvnp1 = args.array<double>("limited_hvnp1");
257 xt::pyarray<int> &csrRowIndeces_DofLoops =
258 args.array<int>("csrRowIndeces_DofLoops");
259 xt::pyarray<int> &csrColumnOffsets_DofLoops =
260 args.array<int>("csrColumnOffsets_DofLoops");
261 xt::pyarray<double> &MassMatrix = args.array<double>("MassMatrix");
262 xt::pyarray<double> &dH_minus_dL = args.array<double>("dH_minus_dL");
263 xt::pyarray<double> &muH_minus_muL = args.array<double>("muH_minus_muL");
264 double hEps = args.scalar<double>("hEps");
265 int LUMPED_MASS_MATRIX = args.scalar<int>("LUMPED_MASS_MATRIX");
266 xt::pyarray<double> &dLow = args.array<double>("dLow");
267 xt::pyarray<double> &new_SourceTerm_hu =
268 args.array<double>("new_SourceTerm_hu");
269 xt::pyarray<double> &new_SourceTerm_hv =
270 args.array<double>("new_SourceTerm_hv");
271 xt::pyarray<double> &hLow = args.array<double>("hLow");
272 xt::pyarray<double> &huLow = args.array<double>("huLow");
273 xt::pyarray<double> &hvLow = args.array<double>("hvLow");
274 xt::pyarray<double> &h_min = args.array<double>("h_min");
275 xt::pyarray<double> &h_max = args.array<double>("h_max");
276 xt::pyarray<double> &kin_max = args.array<double>("kin_max");
277 double KE_tiny = args.scalar<double>("KE_tiny");
278
279 // Declare stuff for limiting on h and h*heta
280 std::valarray<double> Rneg(0.0, numDOFs), Rpos(0.0, numDOFs);
281
282 // Create FCT component matrices in vector form
283 std::valarray<double> FCT_h(0.0, dH_minus_dL.size()),
284 FCT_hu(0.0, dH_minus_dL.size()), FCT_hv(0.0, dH_minus_dL.size());
285
287 // Loop to define FCT matrices for each component //
289 int ij = 0;
290 for (int i = 0; i < numDOFs; i++) {
291 // Read some vectors
292 double high_order_hnp1i = high_order_hnp1[i];
293 double high_order_hunp1i = high_order_hunp1[i];
294 double high_order_hvnp1i = high_order_hvnp1[i];
295 double hi = h_old[i];
296 double huni = hu_old[i];
297 double hvni = hv_old[i];
298 double Zi = b_dof[i];
299 double mi = lumped_mass_matrix[i];
300 double one_over_hiReg =
301 2 * hi / (hi * hi + std::pow(fmax(hi, hEps), 2)); // hEps
302
303 // LOOP OVER THE SPARSITY PATTERN (j-LOOP)//
304 for (int offset = csrRowIndeces_DofLoops[i];
305 offset < csrRowIndeces_DofLoops[i + 1]; offset++) {
306
307 int j = csrColumnOffsets_DofLoops[offset];
308
309 if (i != j) {
310 // Read some vectors
311 double hj = h_old[j];
312 double hunj = hu_old[j];
313 double hvnj = hv_old[j];
314 double Zj = b_dof[j];
315 double one_over_hjReg =
316 2. * hj / (hj * hj + std::pow(fmax(hj, hEps), 2));
317
318 // Compute star states
319 double hStarij = fmax(0., hi + Zi - fmax(Zi, Zj));
320 double huStarij = huni * hStarij * one_over_hiReg;
321 double hvStarij = hvni * hStarij * one_over_hiReg;
322
323 double hStarji = fmax(0., hj + Zj - fmax(Zi, Zj));
324 double huStarji = hunj * hStarji * one_over_hjReg;
325 double hvStarji = hvnj * hStarji * one_over_hjReg;
326
327 // i-th row of flux correction matrix
328 double ML_minus_MC = (LUMPED_MASS_MATRIX == 1
329 ? 0.
330 : (i == j ? 1. : 0.) * mi - MassMatrix[ij]);
331
332 FCT_h[ij] =
333 ML_minus_MC *
334 (high_order_hnp1[j] - hj - (high_order_hnp1i - hi)) +
335 dt * (dH_minus_dL[ij] - muH_minus_muL[ij]) * (hStarji - hStarij) +
336 dt * muH_minus_muL[ij] * (hj - hi);
337
338 FCT_hu[ij] = ML_minus_MC * (high_order_hunp1[j] - hunj -
339 (high_order_hunp1i - huni)) +
340 dt * (dH_minus_dL[ij] - muH_minus_muL[ij]) *
341 (huStarji - huStarij) +
342 dt * muH_minus_muL[ij] * (hunj - huni);
343
344 FCT_hv[ij] = ML_minus_MC * (high_order_hvnp1[j] - hvnj -
345 (high_order_hvnp1i - hvni)) +
346 dt * (dH_minus_dL[ij] - muH_minus_muL[ij]) *
347 (hvStarji - hvStarij) +
348 dt * muH_minus_muL[ij] * (hvnj - hvni);
349
350 } else {
351 FCT_h[ij] = 0.0;
352 FCT_hu[ij] = 0.0;
353 FCT_hv[ij] = 0.0;
354 }
355
356 // UPDATE ij //
357 ij += 1;
358 } // j loop ends here
359 } // i loop ends here
360
362 // Main loop to define limiters and computed limited solution //////
364
365 // Create Lij_array and initialize with 1
366 std::valarray<double> Lij_array(1.0, dH_minus_dL.size());
367
368 /* Loop over limiting iterations */
369 for (int limit_iter = 0; limit_iter < LIMITING_ITERATION; limit_iter++) {
370
371 /* Loop to define FCT Rpos and Rneg values */
372 ij = 0;
373 for (int i = 0; i < numDOFs; i++) {
374
375 double hi = h_old[i];
376 double mi = lumped_mass_matrix[i];
377 // Initialize Pneg and Ppos quantities at ith node
378 double Pnegi = 0., Pposi = 0.;
379
380 // LOOP OVER THE SPARSITY PATTERN (j-LOOP)
381 for (int offset = csrRowIndeces_DofLoops[i];
382 offset < csrRowIndeces_DofLoops[i + 1]; offset++) {
383
384 int j = csrColumnOffsets_DofLoops[offset];
385
386 // COMPUTE P VECTORS
387 Pnegi += FCT_h[ij] * ((FCT_h[ij] < 0) ? 1. : 0.);
388 Pposi += FCT_h[ij] * ((FCT_h[ij] > 0) ? 1. : 0.);
389
390 // UPDATE ij
391 ij += 1;
392 } // j loop ends here
393
394 double psmall_h = 1E-14 * fmax(fabs(Pnegi), fabs(Pposi));
395
397 // COMPUTE Q VECTORS //
399 double Qnegi = std::min(mi * (h_min[i] - hLow[i]), 0.0);
400 double Qposi = std::max(mi * (h_max[i] - hLow[i]), 0.0);
401
403 // COMPUTE R VECTORS //
405 if (hi <= hEps) {
406 Rneg[i] = 0.;
407 Rpos[i] = 0.;
408 } else {
409 // for h
410 if (Pnegi >= -psmall_h) {
411 Rneg[i] = 1.0;
412 } else {
413 Rneg[i] = std::min(1.0, Qnegi / Pnegi);
414 }
415 if (Pposi <= psmall_h) {
416 Rpos[i] = 1.0;
417 } else {
418 Rpos[i] = std::min(1.0, Qposi / Pposi);
419 }
420 }
421 } // i loop ends here
422
423 /* Here we compute the limiters */
424 ij = 0;
425 for (int i = 0; i < numDOFs; i++) {
426
427 double mi = lumped_mass_matrix[i];
428 double ci =
429 kin_max[i] * hLow[i] -
430 0.5 * (huLow[i] * huLow[i] + hvLow[i] * hvLow[i]); // for KE lim.
431
432 // LOOP OVER THE SPARSITY PATTERN (j-LOOP)//
433 for (int offset = csrRowIndeces_DofLoops[i];
434 offset < csrRowIndeces_DofLoops[i + 1]; offset++) {
435
436 int j = csrColumnOffsets_DofLoops[offset];
437
438 if (j != i) {
439
440 // Compute limiter based on water height
441 if (FCT_h[ij] >= 0.) {
442 Lij_array[ij] = fmin(Lij_array[ij], std::min(Rneg[j], Rpos[i]));
443 } else {
444 Lij_array[ij] = fmin(Lij_array[ij], std::min(Rneg[i], Rpos[j]));
445 }
446
447 /*======================================================*/
448 /* Kinetic Energy limiting */
449 double lambdaj =
450 csrRowIndeces_DofLoops[i + 1] - csrRowIndeces_DofLoops[i] - 1;
451 double Ph_ij = FCT_h[ij] / mi / lambdaj;
452 double Phu_ij = FCT_hu[ij] / mi / lambdaj;
453 double Phv_ij = FCT_hv[ij] / mi / lambdaj;
454
455 // Here we initialize limiter based on kinetic energy
456 double KE_limiter = 1;
457 double neg_root_i = 1., neg_root_j = 1.;
458
459 // We first check if local kinetic energy > 0
460 if (kin_max[i] * hLow[i] <= KE_tiny)
461 KE_limiter = 0.;
462
463 // We then check if KE bound is already satisfied
464 double hi_with_lijPij = hLow[i] + Lij_array[ij] * Ph_ij;
465 double hui_with_lijPij = huLow[i] + Lij_array[ij] * Phu_ij;
466 double hvi_with_lijPij = hvLow[i] + Lij_array[ij] * Phv_ij;
467 double psi = kin_max[i] * hi_with_lijPij -
468 0.5 * (hui_with_lijPij * hui_with_lijPij +
469 hvi_with_lijPij * hvi_with_lijPij);
470 if (psi > -KE_tiny) {
471 KE_limiter = fmin(KE_limiter, Lij_array[ij]);
472 }
473
474 /*======================================================*/
475
476 double ai = -0.5 * (Phu_ij * Phu_ij + Phv_ij * Phv_ij);
477 double bi =
478 kin_max[i] * Ph_ij - (huLow[i] * Phu_ij + hvLow[i] * Phv_ij);
479 double delta_i = bi * bi - 4. * ai * ci;
480
481 if (delta_i < 0. || ai >= -0.) {
482 KE_limiter = fmin(KE_limiter, Lij_array[ij]);
483 } else {
484 neg_root_i = (-bi - std::sqrt(delta_i)) / 2. / ai;
485 }
486
487 // root of jth-DOF (To compute transpose component)
488 double lambdai =
489 csrRowIndeces_DofLoops[j + 1] - csrRowIndeces_DofLoops[j] - 1;
490 double mj = lumped_mass_matrix[j];
491 double cj = kin_max[j] * hLow[j] -
492 0.5 * (huLow[j] * huLow[j] + hvLow[j] * hvLow[j]);
493 double Ph_ji = -FCT_h[ij] / mj / lambdai; // Aij=-Aji
494 double Phu_ji = -FCT_hu[ij] / mj / lambdai;
495 double Phv_ji = -FCT_hv[ij] / mj / lambdai;
496 double aj = -0.5 * (Phu_ji * Phu_ji + Phv_ji * Phv_ji);
497 double bj =
498 kin_max[j] * Ph_ji - (huLow[j] * Phu_ji + hvLow[j] * Phv_ji);
499 double delta_j = bj * bj - 4. * aj * cj;
500
501 if (delta_j < 0. || aj >= -0.) {
502 KE_limiter = fmin(KE_limiter, Lij_array[ij]);
503 } else {
504 neg_root_j = (-bj - std::sqrt(delta_j)) / 2. / aj;
505 }
506
507 // define final limiter based on KE
508 KE_limiter =
509 fmin(KE_limiter, fmin(fabs(neg_root_i), fabs(neg_root_j)));
510
511 // Here we set final limiter
512 Lij_array[ij] = fmin(KE_limiter, Lij_array[ij]);
513
514 } else {
515 // if i = j then lij = 0
516 Lij_array[ij] = 0.;
517 }
518
519 // update ij
520 ij += 1;
521 } // end j loop
522 } // end i loop
523
524 /* Final loop to apply limiting and then define limited solution */
525 ij = 0;
526 for (int i = 0; i < numDOFs; i++) {
527
528 double one_over_mi = 1.0 / lumped_mass_matrix[i];
529 double ith_Limiter_times_FluxCorrectionMatrix1 = 0.;
530 double ith_Limiter_times_FluxCorrectionMatrix2 = 0.;
531 double ith_Limiter_times_FluxCorrectionMatrix3 = 0.;
532
533 // LOOP OVER THE SPARSITY PATTERN (j-LOOP)//
534 for (int offset = csrRowIndeces_DofLoops[i];
535 offset < csrRowIndeces_DofLoops[i + 1]; offset++) {
536
537 int j = csrColumnOffsets_DofLoops[offset];
538
539 // COMPUTE LIMITED FLUX //
540 ith_Limiter_times_FluxCorrectionMatrix1 += Lij_array[ij] * FCT_h[ij];
541 ith_Limiter_times_FluxCorrectionMatrix2 += Lij_array[ij] * FCT_hu[ij];
542 ith_Limiter_times_FluxCorrectionMatrix3 += Lij_array[ij] * FCT_hv[ij];
543
544 // update ij
545 ij += 1;
546 } // end j loop
547
548 // then we add lij*Aij to uLow
549 hLow[i] += one_over_mi * ith_Limiter_times_FluxCorrectionMatrix1;
550 huLow[i] += one_over_mi * ith_Limiter_times_FluxCorrectionMatrix2;
551 hvLow[i] += one_over_mi * ith_Limiter_times_FluxCorrectionMatrix3;
552
553 // Finally define the limited solution
554 limited_hnp1[i] = hLow[i];
555 limited_hunp1[i] = huLow[i] + dt * one_over_mi * new_SourceTerm_hu[i];
556 limited_hvnp1[i] = hvLow[i] + dt * one_over_mi * new_SourceTerm_hv[i];
557
558 if (limited_hnp1[i] < -hEps && dt < 1.0) {
559 std::cout << "Limited water height is negative: \n "
560 << "hLow: " << hLow[i] << "\n"
561 << "hHigh: " << limited_hnp1[i] << "\n"
562 << "hEps: " << hEps << "\n"
563 << " ... aborting!" << std::endl;
564 abort();
565 } else {
566 // clean up uHigh from round off error
567 if (limited_hnp1[i] < hEps)
568 limited_hnp1[i] = 0.0;
569 double aux = fmax(limited_hnp1[i], hEps);
570 limited_hunp1[i] *= 2 * std::pow(limited_hnp1[i], VEL_FIX_POWER) /
571 (std::pow(limited_hnp1[i], VEL_FIX_POWER) +
572 std::pow(aux, VEL_FIX_POWER));
573 limited_hvnp1[i] *= 2 * std::pow(limited_hnp1[i], VEL_FIX_POWER) /
574 (std::pow(limited_hnp1[i], VEL_FIX_POWER) +
575 std::pow(aux, VEL_FIX_POWER));
576 }
577 } // end i loop
578
579 // update FCT matrices as Fct = (1 - Lij)*Fct
580 FCT_h = (1.0 - Lij_array) * FCT_h;
581 FCT_hu = (1.0 - Lij_array) * FCT_hu;
582 FCT_hv = (1.0 - Lij_array) * FCT_hv;
583 } // end loop for limiting iteration
584 } // end convex limiting function
585
586 double calculateEdgeBasedCFL(arguments_dict &args) {
587 double g = args.scalar<double>("g");
588 int numDOFsPerEqn = args.scalar<int>("numDOFsPerEqn");
589 xt::pyarray<double> &lumped_mass_matrix =
590 args.array<double>("lumped_mass_matrix");
591 xt::pyarray<double> &h_dof_old = args.array<double>("h_dof_old");
592 xt::pyarray<double> &hu_dof_old = args.array<double>("hu_dof_old");
593 xt::pyarray<double> &hv_dof_old = args.array<double>("hv_dof_old");
594 xt::pyarray<double> &b_dof = args.array<double>("b_dof");
595 xt::pyarray<int> &csrRowIndeces_DofLoops =
596 args.array<int>("csrRowIndeces_DofLoops");
597 xt::pyarray<int> &csrColumnOffsets_DofLoops =
598 args.array<int>("csrColumnOffsets_DofLoops");
599 double hEps = args.scalar<double>("hEps");
600 xt::pyarray<double> &Cx = args.array<double>("Cx");
601 xt::pyarray<double> &Cy = args.array<double>("Cy");
602 xt::pyarray<double> &CTx = args.array<double>("CTx");
603 xt::pyarray<double> &CTy = args.array<double>("CTy");
604 xt::pyarray<double> &dLow = args.array<double>("dLow");
605 double run_cfl = args.scalar<double>("run_cfl");
606 xt::pyarray<double> &edge_based_cfl = args.array<double>("edge_based_cfl");
607 int debug = args.scalar<int>("debug");
608
609 double max_edge_based_cfl = 0.;
610
611 int ij = 0;
612 for (int i = 0; i < numDOFsPerEqn; i++) {
613 // solution at time tn for the ith DOF
614 double hi = h_dof_old[i];
615 double hui = hu_dof_old[i];
616 double hvi = hv_dof_old[i];
617 double dLowii = 0.;
618
619 for (int offset = csrRowIndeces_DofLoops[i];
620 offset < csrRowIndeces_DofLoops[i + 1]; offset++) {
621
622 // loop in j (sparsity pattern)
623 int j = csrColumnOffsets_DofLoops[offset];
624
625 // solution at time tn for the jth DOF
626 double hj = h_dof_old[j];
627 double huj = hu_dof_old[j];
628 double hvj = hv_dof_old[j];
629
630 if (i != j) {
632 // DISSIPATIVE MATRIX //
634 double cij_norm = sqrt(Cx[ij] * Cx[ij] + Cy[ij] * Cy[ij]);
635 double cji_norm = sqrt(CTx[ij] * CTx[ij] + CTy[ij] * CTy[ij]);
636 double nxij = Cx[ij] / cij_norm, nyij = Cy[ij] / cij_norm;
637 double nxji = CTx[ij] / cji_norm, nyji = CTy[ij] / cji_norm;
638 dLow[ij] =
639 fmax(maxWaveSpeedSharpInitialGuess(g, nxij, nyij, hi, hui, hvi,
640 hj, huj, hvj, hEps, debug) *
641 cij_norm, // hEps
642 maxWaveSpeedSharpInitialGuess(g, nxji, nyji, hj, huj, hvj,
643 hi, hui, hvi, hEps, debug) *
644 cji_norm); // hEps
645 dLowii -= dLow[ij];
646 } else
647 dLow[ij] = 0.;
648 // update ij
649 ij += 1;
650 }
652 // CALCULATE EDGE BASED CFL //
654 double mi = lumped_mass_matrix[i];
655 edge_based_cfl[i] = 1.0 * fabs(dLowii) / mi;
656 max_edge_based_cfl = fmax(max_edge_based_cfl, edge_based_cfl[i]);
657 }
658 return max_edge_based_cfl;
659 } // End calculateEdgeBasedCFL
660
661 void calculateEV(arguments_dict &args) {
662 double g = args.scalar<double>("g");
663 xt::pyarray<double> &h_dof_old = args.array<double>("h_dof_old");
664 xt::pyarray<double> &hu_dof_old = args.array<double>("hu_dof_old");
665 xt::pyarray<double> &hv_dof_old = args.array<double>("hv_dof_old");
666 xt::pyarray<double> &b_dof = args.array<double>("b_dof");
667 xt::pyarray<double> &Cx = args.array<double>("Cx");
668 xt::pyarray<double> &Cy = args.array<double>("Cy");
669 xt::pyarray<double> &CTx = args.array<double>("CTx");
670 xt::pyarray<double> &CTy = args.array<double>("CTy");
671 int numDOFsPerEqn = args.scalar<int>("numDOFsPerEqn");
672 xt::pyarray<int> &csrRowIndeces_DofLoops =
673 args.array<int>("csrRowIndeces_DofLoops");
674 xt::pyarray<int> &csrColumnOffsets_DofLoops =
675 args.array<int>("csrColumnOffsets_DofLoops");
676 xt::pyarray<double> &lumped_mass_matrix =
677 args.array<double>("lumped_mass_matrix");
678 double eps = args.scalar<double>("eps");
679 double hEps = args.scalar<double>("hEps");
680 xt::pyarray<double> &global_entropy_residual =
681 args.array<double>("global_entropy_residual");
682 double &dij_small = args.scalar<double>("dij_small");
683
685 // ********** FIRST LOOP ON DOFs ********** //
687
688 // To compute:
689 // * Entropy at i-th node
690 std::valarray<double> eta(numDOFsPerEqn);
691 for (int i = 0; i < numDOFsPerEqn; i++) {
692 // COMPUTE ENTROPY. NOTE: WE CONSIDER A FLAT BOTTOM
693 double hi = h_dof_old[i];
694 double one_over_hiReg =
695 2 * hi / (hi * hi + std::pow(fmax(hi, hEps), 2)); // hEps
696 eta[i] = ENTROPY(g, hi, hu_dof_old[i], hv_dof_old[i], 0., one_over_hiReg);
697 }
698
699 // ********** END OF FIRST LOOP ON DOFs ********** //
700
702 // ********** SECOND LOOP ON DOFs ********** //
704 // To compute:
705 // * global entropy residual
706 // * dij_small to avoid division by 0
707
708 int ij = 0;
709 std::valarray<double> etaMax(numDOFsPerEqn), etaMin(numDOFsPerEqn);
710
711 // speed = sqrt(g max(h_0)), I divide by h_epsilon to get max(h_0)
712 double speed = std::sqrt(g * hEps / eps);
713 dij_small = 0.0;
714
715 for (int i = 0; i < numDOFsPerEqn; i++) {
716
717 // solution at time tn for the ith DOF
718 double hi = h_dof_old[i];
719 double hui = hu_dof_old[i];
720 double hvi = hv_dof_old[i];
721 double Zi = b_dof[i];
722
723 // Define some things using above
724 double one_over_hiReg =
725 2 * hi / (hi * hi + std::pow(fmax(hi, hEps), 2)); // hEps
726 double ui = hui * one_over_hiReg;
727 double vi = hvi * one_over_hiReg;
728 double mi = lumped_mass_matrix[i];
729
730 // initialize etaMax and etaMin
731 etaMax[i] = fabs(eta[i]);
732 etaMin[i] = fabs(eta[i]);
733
734 // FOR ENTROPY RESIDUAL, NOTE: FLAT BOTTOM //
735 double ith_flux_term1 = 0., ith_flux_term2 = 0., ith_flux_term3 = 0.;
736 double entropy_flux = 0.;
737 double sum_entprime_flux = 0.;
738 double eta_prime1 = DENTROPY_DH(g, hi, hui, hvi, 0., one_over_hiReg);
739 double eta_prime2 = DENTROPY_DHU(g, hi, hui, hvi, 0., one_over_hiReg);
740 double eta_prime3 = DENTROPY_DHV(g, hi, hui, hvi, 0., one_over_hiReg);
741
742 // loop in j (sparsity pattern)
743 for (int offset = csrRowIndeces_DofLoops[i];
744 offset < csrRowIndeces_DofLoops[i + 1]; offset++) {
745
746 int j = csrColumnOffsets_DofLoops[offset];
747
748 // solution at time tn for the jth DOF
749 double hj = h_dof_old[j];
750 double huj = hu_dof_old[j];
751 double hvj = hv_dof_old[j];
752 double Zj = b_dof[j];
753
754 // Then define some things here using above
755 double one_over_hjReg =
756 2.0 * hj / (hj * hj + std::pow(fmax(hj, hEps), 2));
757 double uj = huj * one_over_hjReg;
758 double vj = hvj * one_over_hjReg;
759
760 // auxiliary functions to compute fluxes
761 double aux_h =
762 (uj * hj - ui * hi) * Cx[ij] + (vj * hj - vi * hi) * Cy[ij];
763 double aux_hu =
764 (uj * huj - ui * hui) * Cx[ij] + (vj * huj - vi * hui) * Cy[ij];
765 double aux_hv =
766 (uj * hvj - ui * hvi) * Cx[ij] + (vj * hvj - vi * hvi) * Cy[ij];
767
768 // flux for entropy
769 ith_flux_term1 += aux_h;
770 ith_flux_term2 += aux_hu + 0.5 * g * hj * hj * Cx[ij];
771 ith_flux_term3 += aux_hv + 0.5 * g * hj * hj * Cy[ij];
772
773 // NOTE: WE CONSIDER FLAT BOTTOM
774 entropy_flux +=
775 (Cx[ij] * ENTROPY_FLUX1(g, hj, huj, hvj, 0., one_over_hjReg) +
776 Cy[ij] * ENTROPY_FLUX2(g, hj, huj, hvj, 0., one_over_hjReg));
777
778 // COMPUTE ETA MIN AND ETA MAX //
779 etaMax[i] = fmax(etaMax[i], fabs(eta[j]));
780 etaMin[i] = fmin(etaMin[i], fabs(eta[j]));
781
782 // define dij_small in j loop
783 double x = fabs(Cx[ij]) + fabs(Cy[ij]);
784 dij_small = fmax(dij_small, x * speed);
785
786 // update ij
787 ij += 1;
788 } // end j loop
789
790 // define sum of entprime*flux
791 sum_entprime_flux =
792 (ith_flux_term1 * eta_prime1 + ith_flux_term2 * eta_prime2 +
793 ith_flux_term3 * eta_prime3);
794
795 // define rescale for normalization
796 double small_rescale = g * hEps * hEps / eps;
797 double rescale = fmax(fabs(etaMax[i] - etaMin[i]) / 2., small_rescale);
798
799 // COMPUTE ENTROPY RESIDUAL //
800 double one_over_entNormFactori = 1.0 / rescale;
801 global_entropy_residual[i] =
802 one_over_entNormFactori * fabs(entropy_flux - sum_entprime_flux);
803 if (hi <= hEps) {
804 global_entropy_residual[i] = 1.0;
805 }
806 } // end i loop
807
808 // Finally dij_small here
809 dij_small = 1E-14 * dij_small;
810 // ********** END OF LOOP IN DOFs ********** //
811 } // end calculateEV
812
813 void calculateResidual(arguments_dict &args) {
814 xt::pyarray<double> &mesh_trial_ref = args.array<double>("mesh_trial_ref");
815 xt::pyarray<double> &mesh_grad_trial_ref =
816 args.array<double>("mesh_grad_trial_ref");
817 xt::pyarray<double> &mesh_dof = args.array<double>("mesh_dof");
818 xt::pyarray<int> &mesh_l2g = args.array<int>("mesh_l2g");
819 xt::pyarray<double> &dV_ref = args.array<double>("dV_ref");
820 xt::pyarray<double> &h_trial_ref = args.array<double>("h_trial_ref");
821 xt::pyarray<double> &h_grad_trial_ref =
822 args.array<double>("h_grad_trial_ref");
823 xt::pyarray<double> &h_test_ref = args.array<double>("h_test_ref");
824 xt::pyarray<double> &h_grad_test_ref =
825 args.array<double>("h_grad_test_ref");
826 xt::pyarray<double> &vel_trial_ref = args.array<double>("vel_trial_ref");
827 xt::pyarray<double> &vel_grad_trial_ref =
828 args.array<double>("vel_grad_trial_ref");
829 xt::pyarray<double> &vel_test_ref = args.array<double>("vel_test_ref");
830 xt::pyarray<double> &vel_grad_test_ref =
831 args.array<double>("vel_grad_test_ref");
832 xt::pyarray<double> &mesh_trial_trace_ref =
833 args.array<double>("mesh_trial_trace_ref");
834 xt::pyarray<double> &mesh_grad_trial_trace_ref =
835 args.array<double>("mesh_grad_trial_trace_ref");
836 xt::pyarray<double> &h_trial_trace_ref =
837 args.array<double>("h_trial_trace_ref");
838 xt::pyarray<double> &h_grad_trial_trace_ref =
839 args.array<double>("h_grad_trial_trace_ref");
840 xt::pyarray<double> &h_test_trace_ref =
841 args.array<double>("h_test_trace_ref");
842 xt::pyarray<double> &h_grad_test_trace_ref =
843 args.array<double>("h_grad_test_trace_ref");
844 xt::pyarray<double> &vel_trial_trace_ref =
845 args.array<double>("vel_trial_trace_ref");
846 xt::pyarray<double> &vel_grad_trial_trace_ref =
847 args.array<double>("vel_grad_trial_trace_ref");
848 xt::pyarray<double> &vel_test_trace_ref =
849 args.array<double>("vel_test_trace_ref");
850 xt::pyarray<double> &vel_grad_test_trace_ref =
851 args.array<double>("vel_grad_test_trace_ref");
852 xt::pyarray<double> &normal_ref = args.array<double>("normal_ref");
853 xt::pyarray<double> &boundaryJac_ref =
854 args.array<double>("boundaryJac_ref");
855 xt::pyarray<double> &elementDiameter =
856 args.array<double>("elementDiameter");
857 int nElements_global = args.scalar<int>("nElements_global");
858 double g = args.scalar<double>("g");
859 xt::pyarray<int> &h_l2g = args.array<int>("h_l2g");
860 xt::pyarray<int> &vel_l2g = args.array<int>("vel_l2g");
861 xt::pyarray<double> &h_dof_old = args.array<double>("h_dof_old");
862 xt::pyarray<double> &hu_dof_old = args.array<double>("hu_dof_old");
863 xt::pyarray<double> &hv_dof_old = args.array<double>("hv_dof_old");
864 xt::pyarray<double> &b_dof = args.array<double>("b_dof");
865 xt::pyarray<double> &h_dof = args.array<double>("h_dof");
866 xt::pyarray<double> &hu_dof = args.array<double>("hu_dof");
867 xt::pyarray<double> &hv_dof = args.array<double>("hv_dof");
868 xt::pyarray<double> &q_cfl = args.array<double>("q_cfl");
869 xt::pyarray<int> &sdInfo_hu_hu_rowptr =
870 args.array<int>("sdInfo_hu_hu_rowptr");
871 xt::pyarray<int> &sdInfo_hu_hu_colind =
872 args.array<int>("sdInfo_hu_hu_colind");
873 xt::pyarray<int> &sdInfo_hu_hv_rowptr =
874 args.array<int>("sdInfo_hu_hv_rowptr");
875 xt::pyarray<int> &sdInfo_hu_hv_colind =
876 args.array<int>("sdInfo_hu_hv_colind");
877 xt::pyarray<int> &sdInfo_hv_hv_rowptr =
878 args.array<int>("sdInfo_hv_hv_rowptr");
879 xt::pyarray<int> &sdInfo_hv_hv_colind =
880 args.array<int>("sdInfo_hv_hv_colind");
881 xt::pyarray<int> &sdInfo_hv_hu_rowptr =
882 args.array<int>("sdInfo_hv_hu_rowptr");
883 xt::pyarray<int> &sdInfo_hv_hu_colind =
884 args.array<int>("sdInfo_hv_hu_colind");
885 int offset_h = args.scalar<int>("offset_h");
886 int offset_hu = args.scalar<int>("offset_hu");
887 int offset_hv = args.scalar<int>("offset_hv");
888 int stride_h = args.scalar<int>("stride_h");
889 int stride_hu = args.scalar<int>("stride_hu");
890 int stride_hv = args.scalar<int>("stride_hv");
891 xt::pyarray<double> &globalResidual = args.array<double>("globalResidual");
892 int nExteriorElementBoundaries_global =
893 args.scalar<int>("nExteriorElementBoundaries_global");
894 xt::pyarray<int> &exteriorElementBoundariesArray =
895 args.array<int>("exteriorElementBoundariesArray");
896 xt::pyarray<int> &elementBoundaryElementsArray =
897 args.array<int>("elementBoundaryElementsArray");
898 xt::pyarray<int> &elementBoundaryLocalElementBoundariesArray =
899 args.array<int>("elementBoundaryLocalElementBoundariesArray");
900 xt::pyarray<int> &isDOFBoundary_h = args.array<int>("isDOFBoundary_h");
901 xt::pyarray<int> &isDOFBoundary_hu = args.array<int>("isDOFBoundary_hu");
902 xt::pyarray<int> &isDOFBoundary_hv = args.array<int>("isDOFBoundary_hv");
903 xt::pyarray<int> &isAdvectiveFluxBoundary_h =
904 args.array<int>("isAdvectiveFluxBoundary_h");
905 xt::pyarray<int> &isAdvectiveFluxBoundary_hu =
906 args.array<int>("isAdvectiveFluxBoundary_hu");
907 xt::pyarray<int> &isAdvectiveFluxBoundary_hv =
908 args.array<int>("isAdvectiveFluxBoundary_hv");
909 xt::pyarray<int> &isDiffusiveFluxBoundary_hu =
910 args.array<int>("isDiffusiveFluxBoundary_hu");
911 xt::pyarray<int> &isDiffusiveFluxBoundary_hv =
912 args.array<int>("isDiffusiveFluxBoundary_hv");
913 xt::pyarray<double> &ebqe_bc_h_ext = args.array<double>("ebqe_bc_h_ext");
914 xt::pyarray<double> &ebqe_bc_flux_mass_ext =
915 args.array<double>("ebqe_bc_flux_mass_ext");
916 xt::pyarray<double> &ebqe_bc_flux_mom_hu_adv_ext =
917 args.array<double>("ebqe_bc_flux_mom_hu_adv_ext");
918 xt::pyarray<double> &ebqe_bc_flux_mom_hv_adv_ext =
919 args.array<double>("ebqe_bc_flux_mom_hv_adv_ext");
920 xt::pyarray<double> &ebqe_bc_hu_ext = args.array<double>("ebqe_bc_hu_ext");
921 xt::pyarray<double> &ebqe_bc_flux_hu_diff_ext =
922 args.array<double>("ebqe_bc_flux_hu_diff_ext");
923 xt::pyarray<double> &ebqe_penalty_ext =
924 args.array<double>("ebqe_penalty_ext");
925 xt::pyarray<double> &ebqe_bc_hv_ext = args.array<double>("ebqe_bc_hv_ext");
926 xt::pyarray<double> &ebqe_bc_flux_hv_diff_ext =
927 args.array<double>("ebqe_bc_flux_hv_diff_ext");
928 xt::pyarray<double> &q_velocity = args.array<double>("q_velocity");
929 xt::pyarray<double> &ebqe_velocity = args.array<double>("ebqe_velocity");
930 xt::pyarray<double> &flux = args.array<double>("flux");
931 xt::pyarray<double> &elementResidual_h =
932 args.array<double>("elementResidual_h");
933 xt::pyarray<double> &Cx = args.array<double>("Cx");
934 xt::pyarray<double> &Cy = args.array<double>("Cy");
935 xt::pyarray<double> &CTx = args.array<double>("CTx");
936 xt::pyarray<double> &CTy = args.array<double>("CTy");
937 int numDOFsPerEqn = args.scalar<int>("numDOFsPerEqn");
938 int NNZ = args.scalar<int>("NNZ");
939 xt::pyarray<int> &csrRowIndeces_DofLoops =
940 args.array<int>("csrRowIndeces_DofLoops");
941 xt::pyarray<int> &csrColumnOffsets_DofLoops =
942 args.array<int>("csrColumnOffsets_DofLoops");
943 xt::pyarray<double> &lumped_mass_matrix =
944 args.array<double>("lumped_mass_matrix");
945 double cfl_run = args.scalar<double>("cfl_run");
946 double eps = args.scalar<double>("eps");
947 double hEps = args.scalar<double>("hEps");
948 xt::pyarray<double> &hnp1_at_quad_point =
949 args.array<double>("hnp1_at_quad_point");
950 xt::pyarray<double> &hunp1_at_quad_point =
951 args.array<double>("hunp1_at_quad_point");
952 xt::pyarray<double> &hvnp1_at_quad_point =
953 args.array<double>("hvnp1_at_quad_point");
954 xt::pyarray<double> &extendedSourceTerm_hu =
955 args.array<double>("extendedSourceTerm_hu");
956 xt::pyarray<double> &extendedSourceTerm_hv =
957 args.array<double>("extendedSourceTerm_hv");
958 xt::pyarray<double> &dH_minus_dL = args.array<double>("dH_minus_dL");
959 xt::pyarray<double> &muH_minus_muL = args.array<double>("muH_minus_muL");
960 double cE = args.scalar<double>("cE");
961 int LUMPED_MASS_MATRIX = args.scalar<int>("LUMPED_MASS_MATRIX");
962 double dt = args.scalar<double>("dt");
963 int LINEAR_FRICTION = args.scalar<int>("LINEAR_FRICTION");
964 double mannings = args.scalar<double>("mannings");
965 xt::pyarray<double> &quantDOFs = args.array<double>("quantDOFs");
966 int SECOND_CALL_CALCULATE_RESIDUAL =
967 args.scalar<int>("SECOND_CALL_CALCULATE_RESIDUAL");
968 int COMPUTE_NORMALS = args.scalar<int>("COMPUTE_NORMALS");
969 xt::pyarray<double> &normalx = args.array<double>("normalx");
970 xt::pyarray<double> &normaly = args.array<double>("normaly");
971 xt::pyarray<double> &dLow = args.array<double>("dLow");
972 int lstage = args.scalar<int>("lstage");
973 xt::pyarray<double> &new_SourceTerm_hu =
974 args.array<double>("new_SourceTerm_hu");
975 xt::pyarray<double> &new_SourceTerm_hv =
976 args.array<double>("new_SourceTerm_hv");
977 xt::pyarray<double> &global_entropy_residual =
978 args.array<double>("global_entropy_residual");
979 double dij_small = args.scalar<double>("dij_small");
980 xt::pyarray<double> &hLow = args.array<double>("hLow");
981 xt::pyarray<double> &huLow = args.array<double>("huLow");
982 xt::pyarray<double> &hvLow = args.array<double>("hvLow");
983 xt::pyarray<double> &h_min = args.array<double>("h_min");
984 xt::pyarray<double> &h_max = args.array<double>("h_max");
985 xt::pyarray<double> &kin_max = args.array<double>("kin_max");
986 xt::pyarray<double> &urelax = args.array<double>("urelax");
987 xt::pyarray<double> &drelax = args.array<double>("drelax");
988 // FOR FRICTION//
989 double n2 = std::pow(mannings, 2.);
990 double gamma = 4. / 3;
991 double xi = 10.;
992
994 // ********** CELL LOOPS ********** //
996 // To compute:
997 // * Time derivative term
998 // * Cell based CFL
999 // * Velocity and soln at quad points (for other models)
1000 for (int eN = 0; eN < nElements_global; eN++) {
1001 // declare local storage for element residual and initialize
1002 double elementResidual_h[nDOF_test_element],
1003 elementResidual_hu[nDOF_test_element],
1004 elementResidual_hv[nDOF_test_element];
1005
1006 for (int i = 0; i < nDOF_test_element; i++) {
1007 elementResidual_h[i] = 0.0;
1008 elementResidual_hu[i] = 0.0;
1009 elementResidual_hv[i] = 0.0;
1010 }
1011 //
1012 // loop over quadrature points and compute integrands
1013 //
1014 for (int k = 0; k < nQuadraturePoints_element; k++) {
1015 // compute indices and declare local storage
1016 int eN_k = eN * nQuadraturePoints_element + k,
1017 eN_k_nSpace = eN_k * nSpace,
1018 eN_nDOF_trial_element = eN * nDOF_trial_element;
1019 double h = 0.0, hu = 0.0,
1020 hv = 0.0, // solution at current time
1021 h_old = 0.0, hu_old = 0.0, hv_old = 0.0, // solution at lstage
1022 jac[nSpace * nSpace], jacDet, jacInv[nSpace * nSpace],
1023 h_test_dV[nDOF_trial_element], dV, x, y, xt, yt;
1024 // get jacobian, etc for mapping reference element
1025 ck.calculateMapping_element(
1026 eN, k, mesh_dof.data(), mesh_l2g.data(), mesh_trial_ref.data(),
1027 mesh_grad_trial_ref.data(), jac, jacDet, jacInv, x, y);
1028 // get the physical integration weight
1029 dV = fabs(jacDet) * dV_ref[k];
1030 // get the solution at current time
1031 ck.valFromDOF(h_dof.data(), &h_l2g.data()[eN_nDOF_trial_element],
1032 &h_trial_ref.data()[k * nDOF_trial_element], h);
1033 ck.valFromDOF(hu_dof.data(), &vel_l2g[eN_nDOF_trial_element],
1034 &vel_trial_ref.data()[k * nDOF_trial_element], hu);
1035 ck.valFromDOF(hv_dof.data(), &vel_l2g[eN_nDOF_trial_element],
1036 &vel_trial_ref.data()[k * nDOF_trial_element], hv);
1037 // get the solution at the lstage
1038 ck.valFromDOF(h_dof_old.data(), &h_l2g.data()[eN_nDOF_trial_element],
1039 &h_trial_ref.data()[k * nDOF_trial_element], h_old);
1040 ck.valFromDOF(hu_dof_old.data(), &vel_l2g.data()[eN_nDOF_trial_element],
1041 &vel_trial_ref.data()[k * nDOF_trial_element], hu_old);
1042 ck.valFromDOF(hv_dof_old.data(), &vel_l2g.data()[eN_nDOF_trial_element],
1043 &vel_trial_ref.data()[k * nDOF_trial_element], hv_old);
1044 // calculate cell based CFL to keep a reference
1045 calculateCFL(elementDiameter.data()[eN], g, h_old, hu_old, hv_old, hEps,
1046 q_cfl[eN_k]);
1047 // precalculate test function products with integration weights
1048 for (int j = 0; j < nDOF_trial_element; j++)
1049 h_test_dV[j] = h_test_ref[k * nDOF_trial_element + j] * dV;
1050
1051 // SAVE VELOCITY // at quadrature points for other models to use
1052 q_velocity[eN_k_nSpace + 0] =
1053 2 * h / (h * h + std::pow(fmax(h, hEps), 2)) * hu;
1054 q_velocity[eN_k_nSpace + 1] =
1055 2 * h / (h * h + std::pow(fmax(h, hEps), 2)) * hv;
1056 hnp1_at_quad_point[eN_k] = h;
1057 hunp1_at_quad_point[eN_k] = hu;
1058 hvnp1_at_quad_point[eN_k] = hv;
1059
1060 for (int i = 0; i < nDOF_test_element; i++) {
1061 // compute time derivative part of global residual. NOTE: no lumping
1062 elementResidual_h[i] += (h - h_old) * h_test_dV[i];
1063 elementResidual_hu[i] += (hu - hu_old) * h_test_dV[i];
1064 elementResidual_hv[i] += (hv - hv_old) * h_test_dV[i];
1065 }
1066 }
1067 // distribute
1068 for (int i = 0; i < nDOF_test_element; i++) {
1069 int eN_i = eN * nDOF_test_element + i;
1070
1071 // global i-th index for h (this is same for vel_l2g)
1072 int h_gi = h_l2g[eN_i];
1073
1074 // distribute time derivative to global residual
1075 globalResidual[offset_h + stride_h * h_gi] += elementResidual_h[i];
1076 globalResidual[offset_hu + stride_hu * h_gi] += elementResidual_hu[i];
1077 globalResidual[offset_hv + stride_hv * h_gi] += elementResidual_hv[i];
1078 }
1079 }
1080 // ********** END OF CELL LOOPS ********** //
1081
1082 if (SECOND_CALL_CALCULATE_RESIDUAL == 0) // This is to save some time
1083 {
1085 // ********** FIRST SET OF LOOP ON DOFs ********** //
1087 // To compute:
1088 // * Local bounds for limiting
1089 // * Low order solution (in terms of bar states)
1090
1091 // Here we declare some arrays for local bounds
1092 std::valarray<double> delta_Sqd_h(0.0, numDOFsPerEqn),
1093 bar_deltaSqd_h(0.0, numDOFsPerEqn), delta_Sqd_kin(0.0, numDOFsPerEqn),
1094 bar_deltaSqd_kin(0.0, numDOFsPerEqn);
1095 xt::pyarray<double> kin(numDOFsPerEqn), max_of_h_and_hEps(numDOFsPerEqn);
1096
1097 // Define kinetic energy, kin = 1/2 q^2 / h
1098 max_of_h_and_hEps = xt::where(h_dof_old > hEps, h_dof_old, hEps);
1099 kin = 0.5 * (hu_dof_old * hu_dof_old + hv_dof_old * hv_dof_old);
1100 kin = kin * (2.0 * h_dof_old /
1101 (h_dof_old * h_dof_old + max_of_h_and_hEps * max_of_h_and_hEps));
1102
1103 /* First loop to define: delta_Sqd_h, delta_Sqd_kin */
1104 for (int i = 0; i < numDOFsPerEqn; i++) {
1105
1106 for (int offset = csrRowIndeces_DofLoops[i];
1107 offset < csrRowIndeces_DofLoops[i + 1]; offset++) {
1108
1109 int j = csrColumnOffsets_DofLoops[offset];
1110
1111 if (i != j) {
1112 delta_Sqd_h[i] += h_dof_old[i] - h_dof_old[j];
1113 delta_Sqd_kin[i] += kin[i] - kin[j];
1114 }
1115 } // j loop ends here
1116 } // i loops ends here
1117
1118 // Stuff for bar states here (BT = BarTilde)
1119 std::valarray<double> hBT(0.0, dH_minus_dL.size()),
1120 huBT(0.0, dH_minus_dL.size()), hvBT(0.0, dH_minus_dL.size());
1121
1122 /* Second loop to compute bar states (variable)BT */
1123 int ij = 0;
1124 for (int i = 0; i < numDOFsPerEqn; i++) {
1125
1126 // define things at ith node
1127 double hi = h_dof_old[i];
1128 double hui = hu_dof_old[i];
1129 double hvi = hv_dof_old[i];
1130 double Zi = b_dof[i];
1131 double mi = lumped_mass_matrix[i];
1132 double one_over_hiReg =
1133 2.0 * hi / (hi * hi + std::pow(fmax(hi, hEps), 2));
1134 double ui = hui * one_over_hiReg;
1135 double vi = hvi * one_over_hiReg;
1136
1137 // define full pressure at ith node for definition of bar states
1138 double pressure_i = 0.5 * g * hi * hi;
1139
1140 // loop over the sparsity pattern of the i-th DOF
1141 for (int offset = csrRowIndeces_DofLoops[i];
1142 offset < csrRowIndeces_DofLoops[i + 1]; offset++) {
1143 int j = csrColumnOffsets_DofLoops[offset];
1144
1145 // define things at jth node
1146 double hj = h_dof_old[j];
1147 double huj = hu_dof_old[j];
1148 double hvj = hv_dof_old[j];
1149 double Zj = b_dof[j];
1150 double one_over_hjReg =
1151 2.0 * hj / (hj * hj + std::pow(fmax(hj, hEps), 2));
1152 double uj = huj * one_over_hjReg;
1153 double vj = hvj * one_over_hjReg;
1154
1155 // for bar_deltaSqd_h, bar_deltaSqd_heta, bar_deltaSqd_kin
1156 double muLowij = 0., muLij = 0., dLowij = 0., dLij = 0.;
1157
1158 if (i != j) {
1159 // Put these computations first before it gets messy
1160 bar_deltaSqd_h[i] += 0.5 * delta_Sqd_h[j] + 0.5 * delta_Sqd_h[i];
1161 bar_deltaSqd_kin[i] +=
1162 0.5 * delta_Sqd_kin[j] + 0.5 * delta_Sqd_kin[i];
1163
1164 if (lstage == 0)
1165 dLowij = dLow[ij];
1166 else {
1167 double cij_norm = sqrt(Cx[ij] * Cx[ij] + Cy[ij] * Cy[ij]);
1168 double cji_norm = sqrt(CTx[ij] * CTx[ij] + CTy[ij] * CTy[ij]);
1169 double nxij = Cx[ij] / cij_norm, nyij = Cy[ij] / cij_norm;
1170 double nxji = CTx[ij] / cji_norm, nyji = CTy[ij] / cji_norm;
1171 dLowij = fmax(
1172 maxWaveSpeedSharpInitialGuess(g, nxij, nyij, hi, hui, hvi, hj,
1173 huj, hvj, hEps, false) *
1174 cij_norm,
1175 maxWaveSpeedSharpInitialGuess(g, nxji, nyji, hj, huj, hvj, hi,
1176 hui, hvi, hEps, false) *
1177 cji_norm);
1178 }
1179 // save dLij
1180 dLij = dLowij;
1181
1182 // compute muij
1183 muLij = fmax(fmax(0., -(ui * Cx[ij] + vi * Cy[ij])),
1184 fmax(0., (uj * Cx[ij] + vj * Cy[ij])));
1185
1186 // Define dLij as max of dLij and muLij
1187 dLij = fmax(dLowij, muLij);
1188 dLow[ij] = fmax(dLij, muLij);
1189
1191 // COMPUTE BAR STATES //
1193
1194 // define pressure at jth node for bar states
1195 double pressure_j = 0.5 * g * hj * hj;
1196
1197 // Compute star states
1198 double hStarij = fmax(0., hi + Zi - fmax(Zi, Zj));
1199 double huStarij = hui * hStarij * one_over_hiReg;
1200 double hvStarij = hvi * hStarij * one_over_hiReg;
1201 double hStarji = fmax(0., hj + Zj - fmax(Zi, Zj));
1202 double huStarji = huj * hStarji * one_over_hjReg;
1203 double hvStarji = hvj * hStarji * one_over_hjReg;
1204
1205 double hBar_ij = 0., hTilde_ij = 0., huBar_ij = 0., huTilde_ij = 0.,
1206 hvBar_ij = 0., hvTilde_ij = 0.;
1207
1208 // h component
1209 hBar_ij = -1. / (2.0 * fmax(dLij, dij_small)) *
1210 ((uj * hj - ui * hi) * Cx[ij] +
1211 (vj * hj - vi * hi) * Cy[ij]) +
1212 0.5 * (hj + hi);
1213 hTilde_ij = (dLij - muLij) / (2.0 * fmax(dLij, dij_small)) *
1214 (hStarji - hj - (hStarij - hi));
1215 // hu component
1216 huBar_ij =
1217 -1. / (2.0 * fmax(dLij, dij_small)) *
1218 ((uj * huj - ui * hui + pressure_j - pressure_i) * Cx[ij] +
1219 (vj * huj - vi * hui) * Cy[ij]) +
1220 0.5 * (huj + hui);
1221 huTilde_ij = (dLij - muLij) / (2.0 * fmax(dLij, dij_small)) *
1222 (huStarji - huj - (huStarij - hui));
1223 // hv component
1224 hvBar_ij =
1225 -1. / (2.0 * fmax(dLij, dij_small)) *
1226 ((uj * hvj - ui * hvi) * Cx[ij] +
1227 (vj * hvj - vi * hvi + pressure_j - pressure_i) * Cy[ij]) +
1228 0.5 * (hvj + hvi);
1229 hvTilde_ij = (dLij - muLij) / (2.0 * fmax(dLij, dij_small)) *
1230 (hvStarji - hvj - (hvStarij - hvi));
1231
1232 // Here we define uBar + uTilde
1233 hBT[ij] = hBar_ij + hTilde_ij;
1234 huBT[ij] = huBar_ij + huTilde_ij;
1235 hvBT[ij] = hvBar_ij + hvTilde_ij;
1236 } else {
1237 // i==j
1238 // Bar states by definition satisfy Utilde_ii + Ubar_ii = U_i
1239 hBT[ij] = hi;
1240 huBT[ij] = hui;
1241 hvBT[ij] = hvi;
1242 }
1243
1244 // UPDATE ij //
1245 ij += 1;
1246 } // j loop ends here
1247
1248 // for bar_deltaSqd_h, bar_deltaSqd_heta, bar_deltaSqd_kin
1249 bar_deltaSqd_h[i] =
1250 bar_deltaSqd_h[i] /
1251 (csrRowIndeces_DofLoops[i + 1] - csrRowIndeces_DofLoops[i] - 1) /
1252 2.0;
1253 bar_deltaSqd_kin[i] =
1254 bar_deltaSqd_kin[i] /
1255 (csrRowIndeces_DofLoops[i + 1] - csrRowIndeces_DofLoops[i] - 1) /
1256 2.0;
1257
1258 } // i loops ends here
1259
1260 /* Then final loop of first set to get local bounds */
1261 ij = 0;
1262 for (int i = 0; i < numDOFsPerEqn; i++) {
1263
1264 // define m_i
1265 double mi = lumped_mass_matrix[i];
1266
1267 /* Initialize hmin, hmax */
1268 h_min[i] = h_dof_old[i];
1269 h_max[i] = h_dof_old[i];
1270
1271 /* Initialize low order solution */
1272 hLow[i] = h_dof_old[i];
1273 huLow[i] = hu_dof_old[i];
1274 hvLow[i] = hv_dof_old[i];
1275 kin_max[i] = kin[i];
1276
1277 // loop in j (sparsity pattern)
1278 for (int offset = csrRowIndeces_DofLoops[i];
1279 offset < csrRowIndeces_DofLoops[i + 1]; offset++) {
1280
1281 int j = csrColumnOffsets_DofLoops[offset];
1282
1283 double one_over_hBT =
1284 2.0 * hBT[ij] /
1285 (hBT[ij] * hBT[ij] + std::pow(fmax(hBT[ij], hEps), 2));
1286 double psi_ij = one_over_hBT *
1287 (huBT[ij] * huBT[ij] + hvBT[ij] * hvBT[ij]) /
1288 2.0; // Eqn (6.31)
1289
1290 // COMPUTE LOCAL BOUNDS //
1291 kin_max[i] = fmax(psi_ij, kin_max[i]);
1292 h_min[i] = std::min(h_min[i], hBT[ij]);
1293 h_max[i] = std::max(h_max[i], hBT[ij]);
1294
1295 /* COMPUTE LOW ORDER SOLUTION. See EQN 6.23 in SW friction paper */
1296 // This is low order solution WITHOUT sources
1297 if (i != j) {
1298 hLow[i] += h_dof_old[i] * (-dt / mi * 2 * dLow[ij]) +
1299 dt / mi * (2 * dLow[ij] * hBT[ij]);
1300 huLow[i] += hu_dof_old[i] * (-dt / mi * 2 * dLow[ij]) +
1301 dt / mi * (2 * dLow[ij] * huBT[ij]);
1302 hvLow[i] += hv_dof_old[i] * (-dt / mi * 2 * dLow[ij]) +
1303 dt / mi * (2 * dLow[ij] * hvBT[ij]);
1304 }
1305
1306 // UPDATE ij //
1307 ij += 1;
1308 } // j loop ends here
1309
1310 // Then do relaxation of bounds here. If confused, see convex
1311 // limiting paper
1312 kin_max[i] = std::min(urelax[i] * kin_max[i],
1313 kin_max[i] + std::abs(bar_deltaSqd_kin[i]));
1314 h_min[i] = std::max(drelax[i] * h_min[i],
1315 h_min[i] - std::abs(bar_deltaSqd_h[i]));
1316 h_max[i] = std::min(urelax[i] * h_max[i],
1317 h_max[i] + std::abs(bar_deltaSqd_h[i]));
1318
1319 // clean up hLow from round off error
1320 if (hLow[i] < hEps)
1321 hLow[i] = 0.0;
1322 } // i loop ends here
1323
1325 // ********** Second set of loops on dofs ********** //
1327 // To compute:
1328 // * Hyperbolic part of the flux
1329 // * Extended source terms
1330 // * Smoothness indicator
1331
1332 ij = 0;
1333 std::valarray<double> hyp_flux_h(numDOFsPerEqn),
1334 hyp_flux_hu(numDOFsPerEqn), hyp_flux_hv(numDOFsPerEqn),
1335 psi(numDOFsPerEqn), etaMax(numDOFsPerEqn), etaMin(numDOFsPerEqn);
1336
1337 for (int i = 0; i < numDOFsPerEqn; i++) {
1338 // solution at time tn for the ith DOF
1339 double hi = h_dof_old[i];
1340 double hui = hu_dof_old[i];
1341 double hvi = hv_dof_old[i];
1342 double Zi = b_dof[i];
1343 // Define some things using above
1344 double one_over_hiReg =
1345 2 * hi / (hi * hi + std::pow(fmax(hi, hEps), 2)); // hEps
1346 double ui = hui * one_over_hiReg;
1347 double vi = hvi * one_over_hiReg;
1348 double mi = lumped_mass_matrix[i];
1349
1350 /* COMPUTE EXTENDED SOURCE TERMS for all equations:
1351 * Friction terms
1352 * NOTE: Be careful with sign of source terms. Extended sources are on
1353 * left side of equations. "new_SourceTerm" variables are on right
1354 */
1355
1356 // FRICTION
1357 if (LINEAR_FRICTION == 1) {
1358 extendedSourceTerm_hu[i] = mannings * hui * mi;
1359 extendedSourceTerm_hv[i] = mannings * hvi * mi;
1360 // For use in the convex limiting function
1361 // actually didn't need to do this but it helps with signs
1362 new_SourceTerm_hu[i] = -mannings * hui * mi;
1363 new_SourceTerm_hv[i] = -mannings * hvi * mi;
1364 } else {
1365 double veli_norm = std::sqrt(ui * ui + vi * vi);
1366 double hi_to_the_gamma = std::pow(fmax(hi, hEps), gamma);
1367 double friction_aux =
1368 veli_norm == 0.
1369 ? 0.
1370 : (2 * g * n2 * veli_norm * mi /
1371 (hi_to_the_gamma +
1372 fmax(hi_to_the_gamma, xi * g * n2 * dt * veli_norm)));
1373 extendedSourceTerm_hu[i] = friction_aux * hui;
1374 extendedSourceTerm_hv[i] = friction_aux * hvi;
1375 // For use in the convex limiting function
1376 new_SourceTerm_hu[i] = -friction_aux * hui;
1377 new_SourceTerm_hv[i] = -friction_aux * hvi;
1378 double isVeg = (xi-5.173) > 0.0 ? 1.0 : 0.0;
1379 extendedSourceTerm_hu[i] += mi * isVeg * 0.5 * 1.0 * 0.006 * 194 * fmin(0.15,hi)*ui*veli_norm;
1380 extendedSourceTerm_hv[i] += mi * isVeg * 0.5 * 1.0 * 0.006 * 194 * fmin(0.15,hi)*vi*veli_norm;
1381 new_SourceTerm_hu[i] -= mi * isVeg * 0.5 * 1.0 * 0.006 * 194 * fmin(0.15,hi)*ui*veli_norm;
1382 new_SourceTerm_hv[i] -= mi * isVeg * 0.5 * 1.0 * 0.006 * 194 * fmin(0.15,hi)*vi*veli_norm;
1383 }
1384
1385 /* HYPERBOLIC FLUXES */
1386 hyp_flux_h[i] = 0;
1387 hyp_flux_hu[i] = 0;
1388 hyp_flux_hv[i] = 0;
1389
1390 // FOR SMOOTHNESS INDICATOR //
1391 double alphai;
1392 double alpha_numerator = 0;
1393 double alpha_denominator = 0;
1394 double alpha_zero = 0.5; // if only want smoothness
1395 double alpha_factor = 1.0 / (1.0 - alpha_zero);
1396
1397 // loop in j (sparsity pattern)
1398 for (int offset = csrRowIndeces_DofLoops[i];
1399 offset < csrRowIndeces_DofLoops[i + 1]; offset++) {
1400
1401 int j = csrColumnOffsets_DofLoops[offset];
1402
1403 // solution at time tn for the jth DOF
1404 double hj = h_dof_old[j];
1405 double huj = hu_dof_old[j];
1406 double hvj = hv_dof_old[j];
1407 double Zj = b_dof[j];
1408 // Then define some things here using above
1409 double one_over_hjReg =
1410 2.0 * hj / (hj * hj + std::pow(fmax(hj, hEps), 2));
1411 double uj = huj * one_over_hjReg;
1412 double vj = hvj * one_over_hjReg;
1413
1414 // auxiliary functions to compute fluxes
1415 double aux_h =
1416 (uj * hj - ui * hi) * Cx[ij] + (vj * hj - vi * hi) * Cy[ij];
1417 double aux_hu =
1418 (uj * huj - ui * hui) * Cx[ij] + (vj * huj - vi * hui) * Cy[ij];
1419 double aux_hv =
1420 (uj * hvj - ui * hvi) * Cx[ij] + (vj * hvj - vi * hvi) * Cy[ij];
1421
1422 /* HYPERBOLIC FLUX */
1423 hyp_flux_h[i] += aux_h;
1424 hyp_flux_hu[i] += aux_hu;
1425 hyp_flux_hv[i] += aux_hv;
1426
1427 // EXTENDED SOURCE, USING 6.13 //
1428 extendedSourceTerm_hu[i] += g * hi * (hj + Zj) * Cx[ij];
1429 extendedSourceTerm_hv[i] += g * hi * (hj + Zj) * Cy[ij];
1430
1431 new_SourceTerm_hu[i] +=
1432 g * (-hi * (Zj - Zi) + 0.5 * std::pow(hj - hi, 2)) * Cx[ij];
1433 new_SourceTerm_hv[i] +=
1434 g * (-hi * (Zj - Zi) + 0.5 * std::pow(hj - hi, 2)) * Cy[ij];
1435
1436 // FOR SMOOTHNESS INDICATOR //
1437 alpha_numerator += hj - hi;
1438 alpha_denominator += fabs(hj - hi);
1439
1440 // update ij
1441 ij += 1;
1442 } // end j loop
1443
1444 // COMPUTE SMOOTHNESS INDICATOR //
1445 if (hi <= hEps) {
1446 alphai = 1.0;
1447 psi[i] = 1.0;
1448 } else {
1449 // Force alphai=0 in constant states
1450 if (fabs(alpha_numerator) <= hEps) {
1451 alphai = 0.;
1452 } else {
1453 alphai =
1454 (fabs(alpha_numerator) - hEps) / fabs(alpha_denominator - hEps);
1455 }
1456 alphai = fmax(alphai - alpha_zero, 0.0) * alpha_factor;
1457 if (POWER_SMOOTHNESS_INDICATOR == 0)
1458 psi[i] = 1.0;
1459 else
1460 psi[i] = std::pow(alphai, POWER_SMOOTHNESS_INDICATOR);
1461 }
1462 }
1463 // ********** END OF 2nd LOOP ON DOFS ********** //
1464
1466 // ********** MAIN LOOP ON DOFs **********
1467 // To compute:
1468 // * dissipative terms
1469 // * bar states
1471
1472 ij = 0;
1473 for (int i = 0; i < numDOFsPerEqn; i++) {
1474 double hi = h_dof_old[i];
1475 double hui = hu_dof_old[i];
1476 double hvi = hv_dof_old[i];
1477 double Zi = b_dof[i];
1478 double mi = lumped_mass_matrix[i];
1479
1480 double one_over_hiReg =
1481 2.0 * hi / (hi * hi + std::pow(fmax(hi, hEps), 2));
1482 double ui = hui * one_over_hiReg;
1483 double vi = hvi * one_over_hiReg;
1484
1485 // Define full pressure at ith node for definition of bar states below
1486 double pressure_i = 0.5 * g * hi * hi;
1487
1488 // HIGH ORDER DISSIPATIVE TERMS, for Aij matrix
1489 double ith_dHij_minus_muHij_times_hStarStates = 0.,
1490 ith_dHij_minus_muHij_times_huStarStates = 0.,
1491 ith_dHij_minus_muHij_times_hvStarStates = 0.,
1492 ith_muHij_times_hStates = 0., ith_muHij_times_huStates = 0.,
1493 ith_muHij_times_hvStates = 0.;
1494
1495 // loop over the sparsity pattern of the i-th DOF
1496 for (int offset = csrRowIndeces_DofLoops[i];
1497 offset < csrRowIndeces_DofLoops[i + 1]; offset++) {
1498 int j = csrColumnOffsets_DofLoops[offset];
1499 double hj = h_dof_old[j];
1500 double huj = hu_dof_old[j];
1501 double hvj = hv_dof_old[j];
1502 double Zj = b_dof[j];
1503 double one_over_hjReg =
1504 2.0 * hj / (hj * hj + std::pow(fmax(hj, hEps), 2));
1505 double uj = huj * one_over_hjReg;
1506 double vj = hvj * one_over_hjReg;
1507
1508 // define pressure at jth node
1509 double pressure_j = 0.5 * g * hj * hj;
1510
1511 // COMPUTE STAR STATES
1512 double hStarij = fmax(0., hi + Zi - fmax(Zi, Zj));
1513 double huStarij = hui * hStarij * one_over_hiReg;
1514 double hvStarij = hvi * hStarij * one_over_hiReg;
1515
1516 double hStarji = fmax(0., hj + Zj - fmax(Zi, Zj));
1517 double huStarji = huj * hStarji * one_over_hjReg;
1518 double hvStarji = hvj * hStarji * one_over_hjReg;
1519
1520 // Dissipative well balancing term
1521 double muLowij = 0., muLij = 0., muHij = 0.;
1522 double dLowij = 0., dLij = 0., dHij = 0.;
1523 if (i != j) // This is not necessary. See formula for
1524 // ith_dissipative_terms
1525 {
1527 // DISSIPATIVE MATRIX //
1529 if (lstage == 0)
1530 dLowij = dLow[ij];
1531 else {
1532 double cij_norm = sqrt(Cx[ij] * Cx[ij] + Cy[ij] * Cy[ij]);
1533 double cji_norm = sqrt(CTx[ij] * CTx[ij] + CTy[ij] * CTy[ij]);
1534 double nxij = Cx[ij] / cij_norm, nyij = Cy[ij] / cij_norm;
1535 double nxji = CTx[ij] / cji_norm, nyji = CTy[ij] / cji_norm;
1536 dLowij = fmax(
1537 maxWaveSpeedSharpInitialGuess(g, nxij, nyij, hi, hui, hvi, hj,
1538 huj, hvj, hEps, false) *
1539 cij_norm,
1540 maxWaveSpeedSharpInitialGuess(g, nxji, nyji, hj, huj, hvj, hi,
1541 hui, hvi, hEps, false) *
1542 cji_norm);
1543 }
1544 // this is standard low-order dij, can you use dLij =
1545 // dLowij*fmax(psi[i],psi[j]) if want smoothness based as low order
1546 dLij = dLowij;
1547
1549 // WELL BALANCING DISSIPATIVE MATRIX //
1551 muLowij = fmax(fmax(0., -(ui * Cx[ij] + vi * Cy[ij])),
1552 fmax(0., (uj * Cx[ij] + vj * Cy[ij])));
1553 muLij = muLowij;
1554
1555 // Define dLij as low order dijs
1556 muLij = muLowij;
1557 dLij = fmax(dLowij, muLij);
1558
1559 // Then save dLow for limiting step, maybe a bit confusing
1560 dLow[ij] = fmax(dLij, muLij);
1561
1563 // ENTROPY VISCOSITY //
1565 double dEVij = cE * fmax(global_entropy_residual[i],
1566 global_entropy_residual[j]);
1567 dHij = fmin(dLowij, dEVij);
1568 muHij = fmin(muLowij, dEVij);
1569
1570 // compute dij_minus_muij times star solution terms
1571 // see: eqn (6.13)
1572 ith_dHij_minus_muHij_times_hStarStates +=
1573 (dHij - muHij) * (hStarji - hStarij);
1574 ith_dHij_minus_muHij_times_huStarStates +=
1575 (dHij - muHij) * (huStarji - huStarij);
1576 ith_dHij_minus_muHij_times_hvStarStates +=
1577 (dHij - muHij) * (hvStarji - hvStarij);
1578
1579 // compute muij times solution terms
1580 ith_muHij_times_hStates += muHij * (hj - hi);
1581 ith_muHij_times_huStates += muHij * (huj - hui);
1582 ith_muHij_times_hvStates += muHij * (hvj - hvi);
1583
1584 // compute dH_minus_dL
1585 dH_minus_dL[ij] = dHij - dLij;
1586 muH_minus_muL[ij] = muHij - muLij;
1587 } else // i==j
1588 {
1589 dH_minus_dL[ij] =
1590 0.; // Not true but the prod of this times Uj-Ui will be zero
1591 muH_minus_muL[ij] =
1592 0.; // Not true but the prod of this times Uj-Ui will be zero
1593 }
1594 // update ij
1595 ij += 1;
1596 } // j loop ends here
1597
1598 /* Define global residual */
1599 if (LUMPED_MASS_MATRIX == 1) {
1600 globalResidual[offset_h + stride_h * i] =
1601 hi - dt / mi *
1602 (hyp_flux_h[i] - ith_dHij_minus_muHij_times_hStarStates -
1603 ith_muHij_times_hStates);
1604 globalResidual[offset_hu + stride_hu * i] =
1605 hui - dt / mi *
1606 ((hyp_flux_hu[i] + extendedSourceTerm_hu[i]) -
1607 ith_dHij_minus_muHij_times_huStarStates -
1608 ith_muHij_times_huStates);
1609 globalResidual[offset_hv + stride_hv * i] =
1610 hvi - dt / mi *
1611 ((hyp_flux_hv[i] + extendedSourceTerm_hv[i]) -
1612 ith_dHij_minus_muHij_times_hvStarStates -
1613 ith_muHij_times_hvStates);
1614 // clean up potential negative water height due to machine precision
1615 if (globalResidual[offset_h + stride_h * i] >= -hEps &&
1616 globalResidual[offset_h + stride_h * i] < hEps)
1617 globalResidual[offset_h + stride_h * i] = 0;
1618 } else {
1619 // Distribute residual
1620 // NOTE: MASS MATRIX IS CONSISTENT
1621 globalResidual[offset_h + stride_h * i] +=
1622 dt * (hyp_flux_h[i] - ith_dHij_minus_muHij_times_hStarStates -
1623 ith_muHij_times_hStates);
1624 globalResidual[offset_hu + stride_hu * i] +=
1625 dt * (hyp_flux_hu[i] + extendedSourceTerm_hu[i] -
1626 ith_dHij_minus_muHij_times_huStarStates -
1627 ith_muHij_times_huStates);
1628 globalResidual[offset_hv + stride_hv * i] +=
1629 dt * (hyp_flux_hv[i] + extendedSourceTerm_hv[i] -
1630 ith_dHij_minus_muHij_times_hvStarStates -
1631 ith_muHij_times_hvStates);
1632 }
1633 }
1634 // ********** END OF LOOP IN DOFs ********** //
1635 } // end SECOND_CALL_CALCULATE_RESIDUAL
1636
1637 // ********** COMPUTE NORMALS ********** //
1638 if (COMPUTE_NORMALS == 1) {
1639 // This is to identify the normals and create a vector of normal
1640 // components
1641 for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++) {
1642 int
1643 ebN = exteriorElementBoundariesArray[ebNE],
1644 eN = elementBoundaryElementsArray[ebN * 2 + 0],
1645 ebN_local = elementBoundaryLocalElementBoundariesArray[ebN * 2 + 0];
1646 double normal[3];
1647 { // "Loop" in quad points
1648 int kb = 0; // NOTE: I need to consider just one quad point since
1649 // the element is not curved so the normal is constant
1650 // per element
1651 int ebN_local_kb =
1652 ebN_local * nQuadraturePoints_elementBoundary + kb;
1653 double jac_ext[nSpace * nSpace], jacDet_ext,
1654 jacInv_ext[nSpace * nSpace], boundaryJac[nSpace * (nSpace - 1)],
1655 metricTensor[(nSpace - 1) * (nSpace - 1)], metricTensorDetSqrt,
1656 x_ext, y_ext;
1657 /* compute information about mapping from reference element to
1658 * physical element */
1659 ck.calculateMapping_elementBoundary(
1660 eN, ebN_local, kb, ebN_local_kb, mesh_dof.data(), mesh_l2g.data(),
1661 mesh_trial_trace_ref.data(), mesh_grad_trial_trace_ref.data(),
1662 boundaryJac_ref.data(), jac_ext, jacDet_ext, jacInv_ext,
1663 boundaryJac, metricTensor, metricTensorDetSqrt, normal_ref.data(),
1664 normal, x_ext, y_ext);
1665 }
1666 // distribute the normal vectors
1667 for (int i = 0; i < nDOF_test_element; i++) {
1668 int eN_i = eN * nDOF_test_element + i;
1669 int gi = h_l2g[eN_i];
1670 normalx[gi] += 0.5 * normal[0] * (i == ebN_local ? 0. : 1.);
1671 normaly[gi] += 0.5 * normal[1] * (i == ebN_local ? 0. : 1.);
1672 }
1673 }
1674 // normalize
1675 for (int gi = 0; gi < numDOFsPerEqn; gi++) {
1676 double norm_factor =
1677 sqrt(std::pow(normalx[gi], 2) + std::pow(normaly[gi], 2));
1678 if (norm_factor != 0) {
1679 normalx[gi] /= norm_factor;
1680 normaly[gi] /= norm_factor;
1681 }
1682 }
1683 }
1684 // ********** END OF COMPUTING NORMALS ********** //
1685 } // end calculateResidual
1686
1687 void calculateMassMatrix(arguments_dict &args) {
1688 xt::pyarray<double> &mesh_trial_ref = args.array<double>("mesh_trial_ref");
1689 xt::pyarray<double> &mesh_grad_trial_ref =
1690 args.array<double>("mesh_grad_trial_ref");
1691 xt::pyarray<double> &mesh_dof = args.array<double>("mesh_dof");
1692 xt::pyarray<double> &mesh_velocity_dof =
1693 args.array<double>("mesh_velocity_dof");
1694 double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
1695 xt::pyarray<int> &mesh_l2g = args.array<int>("mesh_l2g");
1696 xt::pyarray<double> &dV_ref = args.array<double>("dV_ref");
1697 xt::pyarray<double> &h_trial_ref = args.array<double>("h_trial_ref");
1698 xt::pyarray<double> &h_grad_trial_ref =
1699 args.array<double>("h_grad_trial_ref");
1700 xt::pyarray<double> &h_test_ref = args.array<double>("h_test_ref");
1701 xt::pyarray<double> &h_grad_test_ref =
1702 args.array<double>("h_grad_test_ref");
1703 xt::pyarray<double> &vel_trial_ref = args.array<double>("vel_trial_ref");
1704 xt::pyarray<double> &vel_grad_trial_ref =
1705 args.array<double>("vel_grad_trial_ref");
1706 xt::pyarray<double> &vel_test_ref = args.array<double>("vel_test_ref");
1707 xt::pyarray<double> &vel_grad_test_ref =
1708 args.array<double>("vel_grad_test_ref");
1709 xt::pyarray<double> &mesh_trial_trace_ref =
1710 args.array<double>("mesh_trial_trace_ref");
1711 xt::pyarray<double> &mesh_grad_trial_trace_ref =
1712 args.array<double>("mesh_grad_trial_trace_ref");
1713 xt::pyarray<double> &dS_ref = args.array<double>("dS_ref");
1714 xt::pyarray<double> &h_trial_trace_ref =
1715 args.array<double>("h_trial_trace_ref");
1716 xt::pyarray<double> &h_grad_trial_trace_ref =
1717 args.array<double>("h_grad_trial_trace_ref");
1718 xt::pyarray<double> &h_test_trace_ref =
1719 args.array<double>("h_test_trace_ref");
1720 xt::pyarray<double> &h_grad_test_trace_ref =
1721 args.array<double>("h_grad_test_trace_ref");
1722 xt::pyarray<double> &vel_trial_trace_ref =
1723 args.array<double>("vel_trial_trace_ref");
1724 xt::pyarray<double> &vel_grad_trial_trace_ref =
1725 args.array<double>("vel_grad_trial_trace_ref");
1726 xt::pyarray<double> &vel_test_trace_ref =
1727 args.array<double>("vel_test_trace_ref");
1728 xt::pyarray<double> &vel_grad_test_trace_ref =
1729 args.array<double>("vel_grad_test_trace_ref");
1730 xt::pyarray<double> &normal_ref = args.array<double>("normal_ref");
1731 xt::pyarray<double> &boundaryJac_ref =
1732 args.array<double>("boundaryJac_ref");
1733 xt::pyarray<double> &elementDiameter =
1734 args.array<double>("elementDiameter");
1735 int nElements_global = args.scalar<int>("nElements_global");
1736 double g = args.scalar<double>("g");
1737 xt::pyarray<int> &h_l2g = args.array<int>("h_l2g");
1738 xt::pyarray<int> &vel_l2g = args.array<int>("vel_l2g");
1739 xt::pyarray<double> &b_dof = args.array<double>("b_dof");
1740 xt::pyarray<double> &h_dof = args.array<double>("h_dof");
1741 xt::pyarray<double> &hu_dof = args.array<double>("hu_dof");
1742 xt::pyarray<double> &hv_dof = args.array<double>("hv_dof");
1743 xt::pyarray<double> &q_cfl = args.array<double>("q_cfl");
1744 xt::pyarray<int> &sdInfo_hu_hu_rowptr =
1745 args.array<int>("sdInfo_hu_hu_rowptr");
1746 xt::pyarray<int> &sdInfo_hu_hu_colind =
1747 args.array<int>("sdInfo_hu_hu_colind");
1748 xt::pyarray<int> &sdInfo_hu_hv_rowptr =
1749 args.array<int>("sdInfo_hu_hv_rowptr");
1750 xt::pyarray<int> &sdInfo_hu_hv_colind =
1751 args.array<int>("sdInfo_hu_hv_colind");
1752 xt::pyarray<int> &sdInfo_hv_hv_rowptr =
1753 args.array<int>("sdInfo_hv_hv_rowptr");
1754 xt::pyarray<int> &sdInfo_hv_hv_colind =
1755 args.array<int>("sdInfo_hv_hv_colind");
1756 xt::pyarray<int> &sdInfo_hv_hu_rowptr =
1757 args.array<int>("sdInfo_hv_hu_rowptr");
1758 xt::pyarray<int> &sdInfo_hv_hu_colind =
1759 args.array<int>("sdInfo_hv_hu_colind");
1760 xt::pyarray<int> &csrRowIndeces_h_h = args.array<int>("csrRowIndeces_h_h");
1761 xt::pyarray<int> &csrColumnOffsets_h_h =
1762 args.array<int>("csrColumnOffsets_h_h");
1763 xt::pyarray<int> &csrRowIndeces_h_hu =
1764 args.array<int>("csrRowIndeces_h_hu");
1765 xt::pyarray<int> &csrColumnOffsets_h_hu =
1766 args.array<int>("csrColumnOffsets_h_hu");
1767 xt::pyarray<int> &csrRowIndeces_h_hv =
1768 args.array<int>("csrRowIndeces_h_hv");
1769 xt::pyarray<int> &csrColumnOffsets_h_hv =
1770 args.array<int>("csrColumnOffsets_h_hv");
1771 xt::pyarray<int> &csrRowIndeces_hu_h =
1772 args.array<int>("csrRowIndeces_hu_h");
1773 xt::pyarray<int> &csrColumnOffsets_hu_h =
1774 args.array<int>("csrColumnOffsets_hu_h");
1775 xt::pyarray<int> &csrRowIndeces_hu_hu =
1776 args.array<int>("csrRowIndeces_hu_hu");
1777 xt::pyarray<int> &csrColumnOffsets_hu_hu =
1778 args.array<int>("csrColumnOffsets_hu_hu");
1779 xt::pyarray<int> &csrRowIndeces_hu_hv =
1780 args.array<int>("csrRowIndeces_hu_hv");
1781 xt::pyarray<int> &csrColumnOffsets_hu_hv =
1782 args.array<int>("csrColumnOffsets_hu_hv");
1783 xt::pyarray<int> &csrRowIndeces_hv_h =
1784 args.array<int>("csrRowIndeces_hv_h");
1785 xt::pyarray<int> &csrColumnOffsets_hv_h =
1786 args.array<int>("csrColumnOffsets_hv_h");
1787 xt::pyarray<int> &csrRowIndeces_hv_hu =
1788 args.array<int>("csrRowIndeces_hv_hu");
1789 xt::pyarray<int> &csrColumnOffsets_hv_hu =
1790 args.array<int>("csrColumnOffsets_hv_hu");
1791 xt::pyarray<int> &csrRowIndeces_hv_hv =
1792 args.array<int>("csrRowIndeces_hv_hv");
1793 xt::pyarray<int> &csrColumnOffsets_hv_hv =
1794 args.array<int>("csrColumnOffsets_hv_hv");
1795 xt::pyarray<double> &globalJacobian = args.array<double>("globalJacobian");
1796 int nExteriorElementBoundaries_global =
1797 args.scalar<int>("nExteriorElementBoundaries_global");
1798 xt::pyarray<int> &exteriorElementBoundariesArray =
1799 args.array<int>("exteriorElementBoundariesArray");
1800 xt::pyarray<int> &elementBoundaryElementsArray =
1801 args.array<int>("elementBoundaryElementsArray");
1802 xt::pyarray<int> &elementBoundaryLocalElementBoundariesArray =
1803 args.array<int>("elementBoundaryLocalElementBoundariesArray");
1804 xt::pyarray<int> &isDOFBoundary_h = args.array<int>("isDOFBoundary_h");
1805 xt::pyarray<int> &isDOFBoundary_hu = args.array<int>("isDOFBoundary_hu");
1806 xt::pyarray<int> &isDOFBoundary_hv = args.array<int>("isDOFBoundary_hv");
1807 xt::pyarray<int> &isAdvectiveFluxBoundary_h =
1808 args.array<int>("isAdvectiveFluxBoundary_h");
1809 xt::pyarray<int> &isAdvectiveFluxBoundary_hu =
1810 args.array<int>("isAdvectiveFluxBoundary_hu");
1811 xt::pyarray<int> &isAdvectiveFluxBoundary_hv =
1812 args.array<int>("isAdvectiveFluxBoundary_hv");
1813 xt::pyarray<int> &isDiffusiveFluxBoundary_hu =
1814 args.array<int>("isDiffusiveFluxBoundary_hu");
1815 xt::pyarray<int> &isDiffusiveFluxBoundary_hv =
1816 args.array<int>("isDiffusiveFluxBoundary_hv");
1817 xt::pyarray<double> &ebqe_bc_h_ext = args.array<double>("ebqe_bc_h_ext");
1818 xt::pyarray<double> &ebqe_bc_flux_mass_ext =
1819 args.array<double>("ebqe_bc_flux_mass_ext");
1820 xt::pyarray<double> &ebqe_bc_flux_mom_hu_adv_ext =
1821 args.array<double>("ebqe_bc_flux_mom_hu_adv_ext");
1822 xt::pyarray<double> &ebqe_bc_flux_mom_hv_adv_ext =
1823 args.array<double>("ebqe_bc_flux_mom_hv_adv_ext");
1824 xt::pyarray<double> &ebqe_bc_hu_ext = args.array<double>("ebqe_bc_hu_ext");
1825 xt::pyarray<double> &ebqe_bc_flux_hu_diff_ext =
1826 args.array<double>("ebqe_bc_flux_hu_diff_ext");
1827 xt::pyarray<double> &ebqe_penalty_ext =
1828 args.array<double>("ebqe_penalty_ext");
1829 xt::pyarray<double> &ebqe_bc_hv_ext = args.array<double>("ebqe_bc_hv_ext");
1830 xt::pyarray<double> &ebqe_bc_flux_hv_diff_ext =
1831 args.array<double>("ebqe_bc_flux_hv_diff_ext");
1832 xt::pyarray<int> &csrColumnOffsets_eb_h_h =
1833 args.array<int>("csrColumnOffsets_eb_h_h");
1834 xt::pyarray<int> &csrColumnOffsets_eb_h_hu =
1835 args.array<int>("csrColumnOffsets_eb_h_hu");
1836 xt::pyarray<int> &csrColumnOffsets_eb_h_hv =
1837 args.array<int>("csrColumnOffsets_eb_h_hv");
1838 xt::pyarray<int> &csrColumnOffsets_eb_hu_h =
1839 args.array<int>("csrColumnOffsets_eb_hu_h");
1840 xt::pyarray<int> &csrColumnOffsets_eb_hu_hu =
1841 args.array<int>("csrColumnOffsets_eb_hu_hu");
1842 xt::pyarray<int> &csrColumnOffsets_eb_hu_hv =
1843 args.array<int>("csrColumnOffsets_eb_hu_hv");
1844 xt::pyarray<int> &csrColumnOffsets_eb_hv_h =
1845 args.array<int>("csrColumnOffsets_eb_hv_h");
1846 xt::pyarray<int> &csrColumnOffsets_eb_hv_hu =
1847 args.array<int>("csrColumnOffsets_eb_hv_hu");
1848 xt::pyarray<int> &csrColumnOffsets_eb_hv_hv =
1849 args.array<int>("csrColumnOffsets_eb_hv_hv");
1850 double dt = args.scalar<double>("dt");
1851 //
1852 // loop over elements to compute volume integrals and load them into the
1853 // element Jacobians and global Jacobian
1854 //
1855 for (int eN = 0; eN < nElements_global; eN++) {
1856 double elementJacobian_h_h[nDOF_test_element]
1857 [nDOF_trial_element],
1858 elementJacobian_hu_hu[nDOF_test_element][nDOF_trial_element],
1859 elementJacobian_hv_hv[nDOF_test_element][nDOF_trial_element];
1860 for (int i = 0; i < nDOF_test_element; i++)
1861 for (int j = 0; j < nDOF_trial_element; j++) {
1862 elementJacobian_h_h[i][j] = 0.0;
1863 elementJacobian_hu_hu[i][j] = 0.0;
1864 elementJacobian_hv_hv[i][j] = 0.0;
1865 }
1866 for (int k = 0; k < nQuadraturePoints_element; k++) {
1867 int eN_k = eN * nQuadraturePoints_element +
1868 k, // index to a scalar at a quadrature point
1869 eN_k_nSpace = eN_k * nSpace,
1870 eN_nDOF_trial_element =
1871 eN *
1872 nDOF_trial_element; // index to a vector at a quadrature point
1873
1874 // declare local storage
1875 double jac[nSpace * nSpace], jacDet, jacInv[nSpace * nSpace],
1876 dV, h_test_dV[nDOF_test_element], vel_test_dV[nDOF_test_element], x,
1877 y, xt, yt;
1878 // get jacobian, etc for mapping reference element
1879 ck.calculateMapping_element(
1880 eN, k, mesh_dof.data(), mesh_l2g.data(), mesh_trial_ref.data(),
1881 mesh_grad_trial_ref.data(), jac, jacDet, jacInv, x, y);
1882 // get the physical integration weight
1883 dV = fabs(jacDet) * dV_ref[k];
1884 // precalculate test function products with integration weights
1885 for (int j = 0; j < nDOF_trial_element; j++) {
1886 h_test_dV[j] = h_test_ref[k * nDOF_trial_element + j] * dV;
1887 vel_test_dV[j] = vel_test_ref[k * nDOF_trial_element + j] * dV;
1888 }
1889 for (int i = 0; i < nDOF_test_element; i++) {
1890 int i_nSpace = i * nSpace;
1891 for (int j = 0; j < nDOF_trial_element; j++) {
1892 int j_nSpace = j * nSpace;
1893 elementJacobian_h_h[i][j] +=
1894 h_trial_ref[k * nDOF_trial_element + j] * h_test_dV[i];
1895 elementJacobian_hu_hu[i][j] +=
1896 vel_trial_ref[k * nDOF_trial_element + j] * vel_test_dV[i];
1897 elementJacobian_hv_hv[i][j] +=
1898 vel_trial_ref[k * nDOF_trial_element + j] * vel_test_dV[i];
1899 } // j
1900 } // i
1901 } // k
1902 //
1903 // load into element Jacobian into global Jacobian
1904 //
1905 for (int i = 0; i < nDOF_test_element; i++) {
1906 int eN_i = eN * nDOF_test_element + i;
1907 for (int j = 0; j < nDOF_trial_element; j++) {
1908 int eN_i_j = eN_i * nDOF_trial_element + j;
1909 globalJacobian[csrRowIndeces_h_h[eN_i] +
1910 csrColumnOffsets_h_h[eN_i_j]] +=
1911 elementJacobian_h_h[i][j];
1912 globalJacobian[csrRowIndeces_hu_hu[eN_i] +
1913 csrColumnOffsets_hu_hu[eN_i_j]] +=
1914 elementJacobian_hu_hu[i][j];
1915 globalJacobian[csrRowIndeces_hv_hv[eN_i] +
1916 csrColumnOffsets_hv_hv[eN_i_j]] +=
1917 elementJacobian_hv_hv[i][j];
1918 } // j
1919 } // i
1920 } // elements
1921 }
1922
1923 void calculateLumpedMassMatrix(arguments_dict &args) {
1924 xt::pyarray<double> &mesh_trial_ref = args.array<double>("mesh_trial_ref");
1925 xt::pyarray<double> &mesh_grad_trial_ref =
1926 args.array<double>("mesh_grad_trial_ref");
1927 xt::pyarray<double> &mesh_dof = args.array<double>("mesh_dof");
1928 xt::pyarray<double> &mesh_velocity_dof =
1929 args.array<double>("mesh_velocity_dof");
1930 double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
1931 xt::pyarray<int> &mesh_l2g = args.array<int>("mesh_l2g");
1932 xt::pyarray<double> &dV_ref = args.array<double>("dV_ref");
1933 xt::pyarray<double> &h_trial_ref = args.array<double>("h_trial_ref");
1934 xt::pyarray<double> &h_grad_trial_ref =
1935 args.array<double>("h_grad_trial_ref");
1936 xt::pyarray<double> &h_test_ref = args.array<double>("h_test_ref");
1937 xt::pyarray<double> &h_grad_test_ref =
1938 args.array<double>("h_grad_test_ref");
1939 xt::pyarray<double> &vel_trial_ref = args.array<double>("vel_trial_ref");
1940 xt::pyarray<double> &vel_grad_trial_ref =
1941 args.array<double>("vel_grad_trial_ref");
1942 xt::pyarray<double> &vel_test_ref = args.array<double>("vel_test_ref");
1943 xt::pyarray<double> &vel_grad_test_ref =
1944 args.array<double>("vel_grad_test_ref");
1945 xt::pyarray<double> &mesh_trial_trace_ref =
1946 args.array<double>("mesh_trial_trace_ref");
1947 xt::pyarray<double> &mesh_grad_trial_trace_ref =
1948 args.array<double>("mesh_grad_trial_trace_ref");
1949 xt::pyarray<double> &dS_ref = args.array<double>("dS_ref");
1950 xt::pyarray<double> &h_trial_trace_ref =
1951 args.array<double>("h_trial_trace_ref");
1952 xt::pyarray<double> &h_grad_trial_trace_ref =
1953 args.array<double>("h_grad_trial_trace_ref");
1954 xt::pyarray<double> &h_test_trace_ref =
1955 args.array<double>("h_test_trace_ref");
1956 xt::pyarray<double> &h_grad_test_trace_ref =
1957 args.array<double>("h_grad_test_trace_ref");
1958 xt::pyarray<double> &vel_trial_trace_ref =
1959 args.array<double>("vel_trial_trace_ref");
1960 xt::pyarray<double> &vel_grad_trial_trace_ref =
1961 args.array<double>("vel_grad_trial_trace_ref");
1962 xt::pyarray<double> &vel_test_trace_ref =
1963 args.array<double>("vel_test_trace_ref");
1964 xt::pyarray<double> &vel_grad_test_trace_ref =
1965 args.array<double>("vel_grad_test_trace_ref");
1966 xt::pyarray<double> &normal_ref = args.array<double>("normal_ref");
1967 xt::pyarray<double> &boundaryJac_ref =
1968 args.array<double>("boundaryJac_ref");
1969 xt::pyarray<double> &elementDiameter =
1970 args.array<double>("elementDiameter");
1971 int nElements_global = args.scalar<int>("nElements_global");
1972 double g = args.scalar<double>("g");
1973 xt::pyarray<int> &h_l2g = args.array<int>("h_l2g");
1974 xt::pyarray<int> &vel_l2g = args.array<int>("vel_l2g");
1975 xt::pyarray<double> &b_dof = args.array<double>("b_dof");
1976 xt::pyarray<double> &h_dof = args.array<double>("h_dof");
1977 xt::pyarray<double> &hu_dof = args.array<double>("hu_dof");
1978 xt::pyarray<double> &hv_dof = args.array<double>("hv_dof");
1979 xt::pyarray<double> &q_cfl = args.array<double>("q_cfl");
1980 xt::pyarray<int> &sdInfo_hu_hu_rowptr =
1981 args.array<int>("sdInfo_hu_hu_rowptr");
1982 xt::pyarray<int> &sdInfo_hu_hu_colind =
1983 args.array<int>("sdInfo_hu_hu_colind");
1984 xt::pyarray<int> &sdInfo_hu_hv_rowptr =
1985 args.array<int>("sdInfo_hu_hv_rowptr");
1986 xt::pyarray<int> &sdInfo_hu_hv_colind =
1987 args.array<int>("sdInfo_hu_hv_colind");
1988 xt::pyarray<int> &sdInfo_hv_hv_rowptr =
1989 args.array<int>("sdInfo_hv_hv_rowptr");
1990 xt::pyarray<int> &sdInfo_hv_hv_colind =
1991 args.array<int>("sdInfo_hv_hv_colind");
1992 xt::pyarray<int> &sdInfo_hv_hu_rowptr =
1993 args.array<int>("sdInfo_hv_hu_rowptr");
1994 xt::pyarray<int> &sdInfo_hv_hu_colind =
1995 args.array<int>("sdInfo_hv_hu_colind");
1996 xt::pyarray<int> &csrRowIndeces_h_h = args.array<int>("csrRowIndeces_h_h");
1997 xt::pyarray<int> &csrColumnOffsets_h_h =
1998 args.array<int>("csrColumnOffsets_h_h");
1999 xt::pyarray<int> &csrRowIndeces_h_hu =
2000 args.array<int>("csrRowIndeces_h_hu");
2001 xt::pyarray<int> &csrColumnOffsets_h_hu =
2002 args.array<int>("csrColumnOffsets_h_hu");
2003 xt::pyarray<int> &csrRowIndeces_h_hv =
2004 args.array<int>("csrRowIndeces_h_hv");
2005 xt::pyarray<int> &csrColumnOffsets_h_hv =
2006 args.array<int>("csrColumnOffsets_h_hv");
2007 xt::pyarray<int> &csrRowIndeces_hu_h =
2008 args.array<int>("csrRowIndeces_hu_h");
2009 xt::pyarray<int> &csrColumnOffsets_hu_h =
2010 args.array<int>("csrColumnOffsets_hu_h");
2011 xt::pyarray<int> &csrRowIndeces_hu_hu =
2012 args.array<int>("csrRowIndeces_hu_hu");
2013 xt::pyarray<int> &csrColumnOffsets_hu_hu =
2014 args.array<int>("csrColumnOffsets_hu_hu");
2015 xt::pyarray<int> &csrRowIndeces_hu_hv =
2016 args.array<int>("csrRowIndeces_hu_hv");
2017 xt::pyarray<int> &csrColumnOffsets_hu_hv =
2018 args.array<int>("csrColumnOffsets_hu_hv");
2019 xt::pyarray<int> &csrRowIndeces_hv_h =
2020 args.array<int>("csrRowIndeces_hv_h");
2021 xt::pyarray<int> &csrColumnOffsets_hv_h =
2022 args.array<int>("csrColumnOffsets_hv_h");
2023 xt::pyarray<int> &csrRowIndeces_hv_hu =
2024 args.array<int>("csrRowIndeces_hv_hu");
2025 xt::pyarray<int> &csrColumnOffsets_hv_hu =
2026 args.array<int>("csrColumnOffsets_hv_hu");
2027 xt::pyarray<int> &csrRowIndeces_hv_hv =
2028 args.array<int>("csrRowIndeces_hv_hv");
2029 xt::pyarray<int> &csrColumnOffsets_hv_hv =
2030 args.array<int>("csrColumnOffsets_hv_hv");
2031 xt::pyarray<double> &globalJacobian = args.array<double>("globalJacobian");
2032 int nExteriorElementBoundaries_global =
2033 args.scalar<int>("nExteriorElementBoundaries_global");
2034 xt::pyarray<int> &exteriorElementBoundariesArray =
2035 args.array<int>("exteriorElementBoundariesArray");
2036 xt::pyarray<int> &elementBoundaryElementsArray =
2037 args.array<int>("elementBoundaryElementsArray");
2038 xt::pyarray<int> &elementBoundaryLocalElementBoundariesArray =
2039 args.array<int>("elementBoundaryLocalElementBoundariesArray");
2040 xt::pyarray<int> &isDOFBoundary_h = args.array<int>("isDOFBoundary_h");
2041 xt::pyarray<int> &isDOFBoundary_hu = args.array<int>("isDOFBoundary_hu");
2042 xt::pyarray<int> &isDOFBoundary_hv = args.array<int>("isDOFBoundary_hv");
2043 xt::pyarray<int> &isAdvectiveFluxBoundary_h =
2044 args.array<int>("isAdvectiveFluxBoundary_h");
2045 xt::pyarray<int> &isAdvectiveFluxBoundary_hu =
2046 args.array<int>("isAdvectiveFluxBoundary_hu");
2047 xt::pyarray<int> &isAdvectiveFluxBoundary_hv =
2048 args.array<int>("isAdvectiveFluxBoundary_hv");
2049 xt::pyarray<int> &isDiffusiveFluxBoundary_hu =
2050 args.array<int>("isDiffusiveFluxBoundary_hu");
2051 xt::pyarray<int> &isDiffusiveFluxBoundary_hv =
2052 args.array<int>("isDiffusiveFluxBoundary_hv");
2053 xt::pyarray<double> &ebqe_bc_h_ext = args.array<double>("ebqe_bc_h_ext");
2054 xt::pyarray<double> &ebqe_bc_flux_mass_ext =
2055 args.array<double>("ebqe_bc_flux_mass_ext");
2056 xt::pyarray<double> &ebqe_bc_flux_mom_hu_adv_ext =
2057 args.array<double>("ebqe_bc_flux_mom_hu_adv_ext");
2058 xt::pyarray<double> &ebqe_bc_flux_mom_hv_adv_ext =
2059 args.array<double>("ebqe_bc_flux_mom_hv_adv_ext");
2060 xt::pyarray<double> &ebqe_bc_hu_ext = args.array<double>("ebqe_bc_hu_ext");
2061 xt::pyarray<double> &ebqe_bc_flux_hu_diff_ext =
2062 args.array<double>("ebqe_bc_flux_hu_diff_ext");
2063 xt::pyarray<double> &ebqe_penalty_ext =
2064 args.array<double>("ebqe_penalty_ext");
2065 xt::pyarray<double> &ebqe_bc_hv_ext = args.array<double>("ebqe_bc_hv_ext");
2066 xt::pyarray<double> &ebqe_bc_flux_hv_diff_ext =
2067 args.array<double>("ebqe_bc_flux_hv_diff_ext");
2068 xt::pyarray<int> &csrColumnOffsets_eb_h_h =
2069 args.array<int>("csrColumnOffsets_eb_h_h");
2070 xt::pyarray<int> &csrColumnOffsets_eb_h_hu =
2071 args.array<int>("csrColumnOffsets_eb_h_hu");
2072 xt::pyarray<int> &csrColumnOffsets_eb_h_hv =
2073 args.array<int>("csrColumnOffsets_eb_h_hv");
2074 xt::pyarray<int> &csrColumnOffsets_eb_hu_h =
2075 args.array<int>("csrColumnOffsets_eb_hu_h");
2076 xt::pyarray<int> &csrColumnOffsets_eb_hu_hu =
2077 args.array<int>("csrColumnOffsets_eb_hu_hu");
2078 xt::pyarray<int> &csrColumnOffsets_eb_hu_hv =
2079 args.array<int>("csrColumnOffsets_eb_hu_hv");
2080 xt::pyarray<int> &csrColumnOffsets_eb_hv_h =
2081 args.array<int>("csrColumnOffsets_eb_hv_h");
2082 xt::pyarray<int> &csrColumnOffsets_eb_hv_hu =
2083 args.array<int>("csrColumnOffsets_eb_hv_hu");
2084 xt::pyarray<int> &csrColumnOffsets_eb_hv_hv =
2085 args.array<int>("csrColumnOffsets_eb_hv_hv");
2086 double dt = args.scalar<double>("dt");
2087 //
2088 // loop over elements to compute volume integrals and load them into the
2089 // element Jacobians and global Jacobian
2090 //
2091 for (int eN = 0; eN < nElements_global; eN++) {
2092 double elementJacobian_h_h[nDOF_test_element]
2093 [nDOF_trial_element],
2094 elementJacobian_hu_hu[nDOF_test_element][nDOF_trial_element],
2095 elementJacobian_hv_hv[nDOF_test_element][nDOF_trial_element];
2096 for (int i = 0; i < nDOF_test_element; i++)
2097 for (int j = 0; j < nDOF_trial_element; j++) {
2098 elementJacobian_h_h[i][j] = 0.0;
2099 elementJacobian_hu_hu[i][j] = 0.0;
2100 elementJacobian_hv_hv[i][j] = 0.0;
2101 }
2102 for (int k = 0; k < nQuadraturePoints_element; k++) {
2103 int eN_k = eN * nQuadraturePoints_element +
2104 k, // index to a scalar at a quadrature point
2105 eN_k_nSpace = eN_k * nSpace,
2106 eN_nDOF_trial_element =
2107 eN *
2108 nDOF_trial_element; // index to a vector at a quadrature point
2109
2110 // declare local storage
2111 double jac[nSpace * nSpace], jacDet, jacInv[nSpace * nSpace],
2112 dV, h_test_dV[nDOF_test_element], vel_test_dV[nDOF_test_element], x,
2113 y, xt, yt;
2114 // get jacobian, etc for mapping reference element
2115 ck.calculateMapping_element(
2116 eN, k, mesh_dof.data(), mesh_l2g.data(), mesh_trial_ref.data(),
2117 mesh_grad_trial_ref.data(), jac, jacDet, jacInv, x, y);
2118 // get the physical integration weight
2119 dV = fabs(jacDet) * dV_ref[k];
2120 // precalculate test function products with integration weights
2121 for (int j = 0; j < nDOF_trial_element; j++) {
2122 h_test_dV[j] = h_test_ref[k * nDOF_trial_element + j] * dV;
2123 vel_test_dV[j] = vel_test_ref[k * nDOF_trial_element + j] * dV;
2124 }
2125
2126 for (int i = 0; i < nDOF_test_element; i++) {
2127 int i_nSpace = i * nSpace;
2128 for (int j = 0; j < nDOF_trial_element; j++) {
2129 int j_nSpace = j * nSpace;
2130 elementJacobian_h_h[i][j] += (i == j ? 1.0 : 0.0) * h_test_dV[i];
2131 elementJacobian_hu_hu[i][j] +=
2132 (i == j ? 1.0 : 0.0) * vel_test_dV[i];
2133 elementJacobian_hv_hv[i][j] +=
2134 (i == j ? 1.0 : 0.0) * vel_test_dV[i];
2135 } // j
2136 } // i
2137 } // k
2138 //
2139 // load into element Jacobian into global Jacobian
2140 //
2141 for (int i = 0; i < nDOF_test_element; i++) {
2142 int eN_i = eN * nDOF_test_element + i;
2143 for (int j = 0; j < nDOF_trial_element; j++) {
2144 int eN_i_j = eN_i * nDOF_trial_element + j;
2145 globalJacobian[csrRowIndeces_h_h[eN_i] +
2146 csrColumnOffsets_h_h[eN_i_j]] +=
2147 elementJacobian_h_h[i][j];
2148 globalJacobian[csrRowIndeces_hu_hu[eN_i] +
2149 csrColumnOffsets_hu_hu[eN_i_j]] +=
2150 elementJacobian_hu_hu[i][j];
2151 globalJacobian[csrRowIndeces_hv_hv[eN_i] +
2152 csrColumnOffsets_hv_hv[eN_i_j]] +=
2153 elementJacobian_hv_hv[i][j];
2154 } // j
2155 } // i
2156 } // elements
2157 }
2158}; // namespace proteus
2159
2160inline SW2DCV_base *newSW2DCV(int nSpaceIn, int nQuadraturePoints_elementIn,
2161 int nDOF_mesh_trial_elementIn,
2162 int nDOF_trial_elementIn, int nDOF_test_elementIn,
2163 int nQuadraturePoints_elementBoundaryIn,
2164 int CompKernelFlag) {
2165 return proteus::chooseAndAllocateDiscretization2D<SW2DCV_base, SW2DCV,
2166 CompKernel>(
2167 nSpaceIn, nQuadraturePoints_elementIn, nDOF_mesh_trial_elementIn,
2168 nDOF_trial_elementIn, nDOF_test_elementIn,
2169 nQuadraturePoints_elementBoundaryIn, CompKernelFlag);
2170}
2171} // end namespace proteus
2172
2173#endif
LIMITING_ITERATION
#define LIMITING_ITERATION
Definition GN_SW2DCV.h:17
VEL_FIX_POWER
#define VEL_FIX_POWER
Definition GN_SW2DCV.h:14
u
Double u
Definition Headers.h:89
z
Double * z
Definition Headers.h:49
v
Double v
Definition Headers.h:95
psi
Double psi
Definition Headers.h:78
cE
#define cE
Definition NCLS3P.h:10
proteus::SW2DCV_base
Definition SW2DCV.h:132
proteus::SW2DCV_base::calculateResidual
virtual void calculateResidual(arguments_dict &args)=0
proteus::SW2DCV_base::calculateEV
virtual void calculateEV(arguments_dict &args)=0
proteus::SW2DCV_base::calculateLumpedMassMatrix
virtual void calculateLumpedMassMatrix(arguments_dict &args)=0
proteus::SW2DCV_base::calculateEdgeBasedCFL
virtual double calculateEdgeBasedCFL(arguments_dict &args)=0
proteus::SW2DCV_base::~SW2DCV_base
virtual ~SW2DCV_base()
Definition SW2DCV.h:134
proteus::SW2DCV_base::calculateMassMatrix
virtual void calculateMassMatrix(arguments_dict &args)=0
proteus::SW2DCV_base::convexLimiting
virtual void convexLimiting(arguments_dict &args)=0
proteus::SW2DCV::calculateEdgeBasedCFL
double calculateEdgeBasedCFL(arguments_dict &args)
Definition SW2DCV.h:586
proteus::SW2DCV::convexLimiting
void convexLimiting(arguments_dict &args)
Definition SW2DCV.h:234
proteus::SW2DCV::ck
CompKernelType ck
Definition SW2DCV.h:149
proteus::SW2DCV::nDOF_test_X_trial_element
const int nDOF_test_X_trial_element
Definition SW2DCV.h:148
proteus::SW2DCV::calculateMassMatrix
void calculateMassMatrix(arguments_dict &args)
Definition SW2DCV.h:1687
proteus::SW2DCV::maxWaveSpeedSharpInitialGuess
double maxWaveSpeedSharpInitialGuess(double g, double nx, double ny, double hL, double huL, double hvL, double hR, double huR, double hvR, double hEps, bool debugging)
Definition SW2DCV.h:160
proteus::SW2DCV::calculateResidual
void calculateResidual(arguments_dict &args)
Definition SW2DCV.h:813
proteus::SW2DCV::calculateCFL
void calculateCFL(const double &elementDiameter, const double &g, const double &h, const double &hu, const double &hv, const double hEps, double &cfl)
Definition SW2DCV.h:215
proteus::SW2DCV::calculateEV
void calculateEV(arguments_dict &args)
Definition SW2DCV.h:661
proteus::SW2DCV::calculateLumpedMassMatrix
void calculateLumpedMassMatrix(arguments_dict &args)
Definition SW2DCV.h:1923
c
#define c(i)
Definition jf.h:21
POWER_SMOOTHNESS_INDICATOR
#define POWER_SMOOTHNESS_INDICATOR
Definition m_comp_co2.h:22
proteus
Definition ADR.h:19
proteus::phi
double phi(const double &g, const double &h, const double &hL, const double &hR, const double &uL, const double &uR)
Definition SW2DCV.h:62
proteus::ENTROPY_FLUX1
double ENTROPY_FLUX1(const double &g, const double &h, const double &hu, const double &hv, const double &z, const double &one_over_hReg)
Definition GN_SW2DCV.h:71
proteus::phiDiff
double phiDiff(const double &g, const double &h1k, const double &h2k, const double &hL, const double &hR, const double &uL, const double &uR)
Definition SW2DCV.h:88
proteus::phiDDiff2
double phiDDiff2(const double &g, const double &h1k, const double &h2k, const double &hL, const double &hR, const double &uL, const double &uR)
Definition SW2DCV.h:100
proteus::phip
double phip(const double &g, const double &h, const double &hL, const double &hR)
Definition SW2DCV.h:72
proteus::nu1
double nu1(const double &g, const double &hStar, const double &hL, const double &uL, const double &one_over_hL)
Definition SW2DCV.h:76
proteus::ENTROPY
ENTROPY
Definition TADR.h:36
proteus::newSW2DCV
SW2DCV_base * newSW2DCV(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition SW2DCV.h:2160
proteus::DENTROPY_DHU
double DENTROPY_DHU(const double &g, const double &h, const double &hu, const double &hv, const double &z, const double &one_over_hReg)
Definition GN_SW2DCV.h:61
proteus::hStarLFromQuadPhiFromAbove
double hStarLFromQuadPhiFromAbove(const double &g, const double &hStarL, const double &hStarR, const double &hL, const double &hR, const double &uL, const double &uR)
Definition SW2DCV.h:106
proteus::phiDDiff1
double phiDDiff1(const double &g, const double &h1k, const double &h2k, const double &hL, const double &hR, const double &uL, const double &uR)
Definition SW2DCV.h:94
proteus::DENTROPY_DHV
double DENTROPY_DHV(const double &g, const double &h, const double &hu, const double &hv, const double &z, const double &one_over_hReg)
Definition GN_SW2DCV.h:66
proteus::hStarRFromQuadPhiFromBelow
double hStarRFromQuadPhiFromBelow(const double &g, const double &hStarL, const double &hStarR, const double &hL, const double &hR, const double &uL, const double &uR)
Definition SW2DCV.h:117
proteus::DENTROPY_DH
double DENTROPY_DH(const double &g, const double &h, const double &hu, const double &hv, const double &z, const double &one_over_hReg)
Definition GN_SW2DCV.h:56
proteus::ENTROPY_FLUX2
double ENTROPY_FLUX2(const double &g, const double &h, const double &hu, const double &hv, const double &z, const double &one_over_hReg)
Definition GN_SW2DCV.h:78
proteus::nu3
double nu3(const double &g, const double &hStar, const double &hR, const double &uR, const double &one_over_hR)
Definition SW2DCV.h:82
proteus::chooseAndAllocateDiscretization2D
Model_Base * chooseAndAllocateDiscretization2D(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nDOF_v_trial_elementIn, int nDOF_v_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition MixedModelFactory.h:309
proteus::fp
double fp(const double &g, const double &h, const double &hZ)
Definition SW2DCV.h:66
proteus::f
double f(const double &g, const double &h, const double &hZ)
Definition SW2DCV.h:58
proteus::arguments_dict
Definition ArgumentsDict.h:70
proteus::arguments_dict::scalar
T & scalar(const std::string &key)
proteus::arguments_dict::array
xt::pyarray< T > & array(const std::string &key)
Generated on for proteus by doxygen 1.18.0