proteus 1.9.0
C/C++/Fortran libraries
Loading...
Searching...
No Matches
Kappa.h
Go to the documentation of this file.
1#ifndef Kappa_H
2#define Kappa_H
3#include <cmath>
4#include <iostream>
5#include "CompKernel.h"
6#include "ModelFactory.h"
7#include "SedClosure.h"
8#include "ArgumentsDict.h"
9#include "xtensor-python/pyarray.hpp"
10
11namespace py = pybind11;
12
13namespace proteus
14{
15 class Kappa_base
16 {
17 //The base class defining the interface
18 public:
19 virtual ~Kappa_base(){}
20 virtual void setSedClosure(double aDarcy,
21 double betaForch,
22 double grain,
23 double packFraction,
24 double packMargin,
25 double maxFraction,
26 double frFraction,
27 double sigmaC,
28 double C3e,
29 double C4e,
30 double eR,
31 double fContact,
32 double mContact,
33 double nContact,
34 double angFriction,
35 double vos_limiter,
36 double mu_fr_limiter){}
37 virtual void calculateResidual(arguments_dict& args)=0;
38 virtual void calculateJacobian(arguments_dict& args)=0; //VRANS
39 };
40
41 template<class CompKernelType,
42 int nSpace,
43 int nQuadraturePoints_element,
44 int nDOF_mesh_trial_element,
45 int nDOF_trial_element,
46 int nDOF_test_element,
47 int nQuadraturePoints_elementBoundary>
48 class Kappa : public Kappa_base
49 {
50 public:
51 cppHsuSedStress<3> closure;
52 const int nDOF_test_X_trial_element;
53 CompKernelType ck;
54 Kappa():
55 closure(150.0,
56 0.0,
57 0.0102,
58 0.2,
59 0.01,
60 0.635,
61 0.57,
62 1.1,
63 1.2,
64 1.0,
65 0.8,
66 0.02,
67 2.0,
68 5.0,
69 M_PI/6.,
70 0.05,
71 1.00),
72 nDOF_test_X_trial_element(nDOF_test_element*nDOF_trial_element),
73 ck()
74 {
75 }
76
77 void setSedClosure(double aDarcy,
78 double betaForch,
79 double grain,
80 double packFraction,
81 double packMargin,
82 double maxFraction,
83 double frFraction,
84 double sigmaC,
85 double C3e,
86 double C4e,
87 double eR,
88 double fContact,
89 double mContact,
90 double nContact,
91 double angFriction,
92 double vos_limiter,
93 double mu_fr_limiter)
94 {
95 closure = cppHsuSedStress<3>(aDarcy,
96 betaForch,
97 grain,
98 packFraction,
99 packMargin,
100 maxFraction,
101 frFraction,
102 sigmaC,
103 C3e,
104 C4e,
105 eR,
106 fContact,
107 mContact,
108 nContact,
109 angFriction,
110 vos_limiter,
111 mu_fr_limiter);
112 }
113
114 inline
115 void computeK_OmegaCoefficients(const double& div_eps,
116 const double& k,
117 const double& omega,
118 const double grad_k[nSpace],
119 const double grad_omega[nSpace],
120 const double grad_vx[nSpace], //gradient of x component of velocity
121 const double grad_vy[nSpace], //gradient of x component of velocity
122 double& inverse_sigma_k,
123 double& inverse_sigma_omega,
124 double& beta_star,
125 double& beta,
126 double& gamma)
127 {
128 //take these from NASA Langley Turbulence Model page
129 //brute force just to see if I can figure it out
130 //use inverse of sigma_k to match standard k-epsilon form
131 inverse_sigma_k = 2.0; inverse_sigma_omega=2.0; gamma = 13.0/25.0;
132 const double beta0_star = 0.09; const double beta0 = 9.0/125.0;
133 double Omega[nSpace][nSpace] = {{0.,0.},
134 {0.,0.}};
135 double S[nSpace][nSpace] = {{0.,0.},
136 {0.,0.}};
137
138
139 //Omega_ij = (\pd{u_i}{x_j} - \pd{u_j}{x_i})/2
140 Omega[0][1] = 0.5*(grad_vx[1]-grad_vy[0]);
141 Omega[1][0] =-Omega[0][1];
142
143 //S_ij = (\pd{u_i}{x_j} + \pd{u_j}{x_i})/2
144 S[0][0] = grad_vx[0]; S[0][1] = 0.5*(grad_vx[1]+grad_vy[0]);
145 S[1][0] = S[0][1]; S[1][1] = grad_vy[1];
146
147 double chi_omega = 0.0;
148 for (int i=0; i < nSpace; i++)
149 for (int k=0; k < nSpace; k++)
150 for (int j=0; j < nSpace; j++)
151 chi_omega += Omega[i][j]*Omega[j][k]*S[k][i];
152
153 if (fabs(omega) > div_eps)
154 {
155 chi_omega = fabs(chi_omega/(beta0_star*omega*beta0_star*omega*beta0_star*omega));
156
157 const double f_beta = (1.0+70.0*chi_omega)/(1.0 + 80.0*chi_omega);
158 beta = beta0*f_beta;
159 }
160 else
161 {
162 beta = beta0;
163 }
164
165 double chi_k = grad_k[0]*grad_omega[0] + grad_k[1]*grad_omega[1];
166 double f_beta_star = 1.0;
167
168 const double omega3 = omega*omega*omega;
169 if (fabs(omega3) > div_eps)
170 {
171 chi_k = chi_k/omega3;
172 f_beta_star = (1.0 + 680.0*chi_k*chi_k)/(1.0 + 400.0*chi_k*chi_k);
173 }
174 else if (chi_k > 0.0)
175 f_beta_star = 680.0/400.0;
176
177 beta_star = beta0_star*f_beta_star;
178 //if (beta < 0.875*beta0 || beta > beta0)
179 // {
180 // std::cout<<"Kappa2D K-Omega coef problem k= "<<k<<" omega= "<<omega<<" beta= "<<beta<<" beta0= "<<beta0 <<" chi_omega= "<<chi_omega<<std::endl;
181 // }
182 beta = fmax(0.875*beta0,fmin(beta,beta0));
183 //if (beta_star < beta0_star || beta_star > (680.0+1.0e-4)/400.0*beta0_star)
184 //{
185 // std::cout<<"Kappa2D K-Omega coef problem k= "<<k<<" omega= "<<omega<<" beta_star= "<<beta_star<<" beta0_star= "<<beta0_star <<" chi_k= "<<chi_k<<std::endl;
186 //}
187 beta_star = fmax(beta0_star,fmin(beta_star,(680.0/400.0)*beta0_star));
188 //mwf hack
189 //beta = beta0; beta_star = beta0_star; //0.0;
190
191 }
192 //Try Lew, Buscaglia approximation
193 inline
194 void evaluateCoefficients( double v[nSpace],
195 const double eps_mu,
196 const double phi,
197 const double nu_0,
198 const double nu_1,
199 const double sigma_k,
200 const double c_mu,
201 const double grad_vx[nSpace], //gradient of x component of velocity
202 const double grad_vy[nSpace], //gradient of x component of velocity
203 const double grad_vz[nSpace], //gradient of x component of velocity
204 const double& k,
205 double& k_old,
206 double& dissipation,
207 const double& porosity,
208// Argumentlist for sediment
209 int sedFlag,
210 double q_vos,
211 double q_vos_gradc[nSpace],
212 double rho_f,
213 double rho_s,
214 double vs[nSpace],
215 double g[nSpace],
216 //end sediment
217 int dissipation_model_flag,
218 const double grad_k_old[nSpace],
219 const double grad_dissipation[nSpace],
220 double& m,
221 double& dm,
222 double f[nSpace],
223 double df[nSpace],
224 double& a,
225 double& da_dk,
226 double& r,
227 double& dr_dk)
228 {
229
230 double nu_t=0.0,dnu_t_dk=0.0,PiD4=0.0;
231 double gamma_k=0.0,F_k=0.0,sigma_a=sigma_k,kSed=0.,dkSed=0.;
232 //either K-Epsilon or K-Omega
233 const double isKEpsilon = (dissipation_model_flag>=2) ? 0.0 : 1.0;
234 m = k*porosity;
235 dm = porosity;
236
237 for (int I=0; I < nSpace; I++)
238 {
239 f[I] = v[I]*k*porosity;
240 df[I] = v[I]*porosity;
241 }
242 const double H_mu = smoothedHeaviside(eps_mu,phi);
243 double nu = (1.0-H_mu)*nu_0 + H_mu*nu_1;
244 const double div_eps = 1.0e-2*fmin(nu_0,nu_1);
245 //eddy viscosity
246 nu_t = isKEpsilon*c_mu*k_old*k_old/(fabs(dissipation)+div_eps)
247 + (1.0-isKEpsilon)*k_old/(fabs(dissipation)+div_eps);
248
249 //if (nu_t > 1.e6*nu)
250 //{
251 // std::cout<<"Kappa WARNING isKEpsilon = "<<isKEpsilon<<" nu_t = " <<nu_t<<" nu= "<<nu<<" k= "<<k<<" k_old= "<<k_old<<" dissipation= "<<dissipation<<std::endl;
252 //}
253 nu_t = fmax(nu_t,1.e-4*nu);
254 //mwf hack
255 nu_t = fmin(nu_t,1.0e6*nu);
256
257 dnu_t_dk = 0.0;
258
259 //Production term
260 PiD4 = 2.0*(grad_vx[0]*grad_vx[0] +
261 grad_vy[1]*grad_vy[1] +
262 grad_vz[2]*grad_vz[2])
263 +
264 (grad_vx[1] + grad_vy[0])*(grad_vx[1] + grad_vy[0])
265 +
266 (grad_vx[2] + grad_vz[0])*(grad_vx[2] + grad_vz[0])
267 +
268 (grad_vy[2] + grad_vz[1])*(grad_vy[2] + grad_vz[1]);
269
270 //Sediment terms
271 double theta = 1e-10; //Granural temperature- currently set to (almost) zero.
272 //Response time only controled by drag, not collisions
273 //Switch on when collision stress model is on.
274
275 if (sedFlag == 1 && isKEpsilon > 0)
276 {
277
278 double kp = k_old;
279 kSed = closure.kappa_sed1(
280 q_vos,
281 rho_f,
282 rho_s,
283 v,
284 vs,
285 q_vos_gradc,
286 nu,
287 theta,
288 kp,
289 dissipation,
290 nu_t,
291 g);
292 dkSed = closure.dkappa_sed1_dk(
293 q_vos,
294 rho_f,
295 rho_s,
296 v,
297 vs,
298 q_vos_gradc,
299 nu,
300 theta,
301 kp,
302 dissipation,
303 nu_t);
304
305 }
306
307
308
309
310 //K-Omega, 1998
311 if (dissipation_model_flag==2)
312 {
313 //temporaries
314 double sigma_omega=1.,beta_star=1.,beta=1.,gamma=1.;
315 computeK_OmegaCoefficients(div_eps,
316 k_old,
317 dissipation,
318 grad_k_old,
319 grad_dissipation,
320 grad_vx,
321 grad_vy,
322 grad_vz,
323 sigma_a,
324 sigma_omega,
325 beta_star,
326 beta,
327 gamma);
328 //straight forward lagging
329 gamma_k=fmax(beta_star*dissipation,0.0);
330 //gamma_k = fmax(beta_star*k_old/nu_t,0.0);
331 }
332 else if (dissipation_model_flag==3) //K-Omega 1988
333 {
334 sigma_a=2.0; //1/sigma_k,
335 double beta_star = 0.09;
336 gamma_k=fmax(beta_star*dissipation,0.0);
337 }
338 else
339 {
340 //K-Epsilon
341 gamma_k = fmax(c_mu*k_old/nu_t,0.0);
342 sigma_a = sigma_k;
343 }
344 //note sigma_a is calculated as inverse of what is in Wilcox's standard K-Omega formulation to
345 //match standard K-Epsilon formulation where divide by sigma
346 a = porosity*(nu_t/sigma_a + nu);
347 da_dk = porosity*dnu_t_dk/sigma_a;
348
349 F_k = nu_t*PiD4;
350 r = -porosity*F_k + porosity*gamma_k*k +porosity*kSed;
351 dr_dk = porosity*(gamma_k+dkSed);
352
353 }
354 inline
355 void computeK_OmegaCoefficients(const double& div_eps,
356 const double& k,
357 const double& omega,
358 const double grad_k[nSpace],
359 const double grad_omega[nSpace],
360 const double grad_vx[nSpace], //gradient of x component of velocity
361 const double grad_vy[nSpace], //gradient of x component of velocity
362 const double grad_vz[nSpace], //gradient of x component of velocity
363 double& inverse_sigma_k,
364 double& inverse_sigma_omega,
365 double& beta_star,
366 double& beta,
367 double& gamma)
368 {
369 //take these from NASA Langley Turbulence Model page
370 //brute force just to see if I can figure it out
371 //use inverse of sigma_k to match standard k-epsilon form
372 inverse_sigma_k = 2.0; inverse_sigma_omega=2.0; gamma = 13.0/25.0;
373 const double beta0_star = 0.09; const double beta0 = 9.0/125.0;
374 double Omega[nSpace][nSpace] = {{0.,0.,0,},
375 {0.,0.,0.},
376 {0.,0.,0.}};
377 double S[nSpace][nSpace] = {{0.0,0.,0},
378 {0.,0.,0.},
379 {0.,0.,0.}};
380
381 //Omega_ij = (\pd{u_i}{x_j} - \pd{u_j}{x_i})/2
382 Omega[0][1] = 0.5*(grad_vx[1]-grad_vy[0]); Omega[0][2] = 0.5*(grad_vx[2]-grad_vz[0]);
383 Omega[1][0] =-Omega[0][1]; Omega[1][2] = 0.5*(grad_vy[2]-grad_vz[1]);
384 Omega[2][0] =-Omega[0][2]; Omega[2][1] =-Omega[1][2];
385 //S_ij = (\pd{u_i}{x_j} + \pd{u_j}{x_i})/2
386 S[0][0] = grad_vx[0]; S[0][1] = 0.5*(grad_vx[1]+grad_vy[0]); S[0][2] = 0.5*(grad_vx[2]+grad_vz[0]);
387 S[1][0] = S[0][1]; S[1][1] = grad_vy[1]; S[1][2] = 0.5*(grad_vy[2] + grad_vz[1]);
388 S[2][0] = S[0][2]; S[2][1] = S[1][2]; S[2][2] = grad_vz[2];
389
390 double chi_omega = 0.0;
391 for (int i=0; i < nSpace; i++)
392 for (int k=0; k < nSpace; k++)
393 for (int j=0; j < nSpace; j++)
394 chi_omega += Omega[i][j]*Omega[j][k]*S[k][i];
395
396 if (fabs(omega) > div_eps)
397 {
398 chi_omega = fabs(chi_omega/(beta0_star*omega*beta0_star*omega*beta0_star*omega));
399
400 const double f_beta = (1.0+70.0*chi_omega)/(1.0 + 80.0*chi_omega);
401 beta = beta0*f_beta;
402 }
403 else
404 {
405 beta = beta0;
406 }
407
408 double chi_k = grad_k[0]*grad_omega[0] + grad_k[1]*grad_omega[1] + grad_k[2]*grad_omega[2];
409 double f_beta_star = 1.0;
410
411 const double omega3 = omega*omega*omega;
412 if (fabs(omega3) > div_eps)
413 {
414 chi_k = chi_k/omega3;
415 f_beta_star = (1.0 + 680.0*chi_k*chi_k)/(1.0 + 400.0*chi_k*chi_k);
416 }
417 else if (chi_k > 0.0)
418 f_beta_star = 680.0/400.0;
419
420 beta_star = beta0_star*f_beta_star;
421 //if (beta < 0.875*beta0 || beta > beta0)
422 // {
423 // std::cout<<"Kappa K-Omega coef problem k= "<<k<<" omega= "<<omega<<" beta= "<<beta<<" beta0= "<<beta0 <<" chi_omega= "<<chi_omega<<std::endl;
424 // }
425 beta = fmax(0.875*beta0,fmin(beta,beta0));
426 //if (beta_star < beta0_star || beta_star > (680.0+1.0e-4)/400.0*beta0_star)
427 //{
428 // std::cout<<"Kappa K-Omega coef problem k= "<<k<<" omega= "<<omega<<" beta_star= "<<beta_star<<" beta0_star= "<<beta0_star <<" chi_k= "<<chi_k<<std::endl;
429 //}
430 beta_star = fmax(beta0_star,fmin(beta_star,(680.0/400.0)*beta0_star));
431 //mwf hack
432 //beta = beta0; beta_star = beta0_star; //0.0;
433
434 }
435
436
437 inline
438 void calculateSubgridError_tau(const double& elementDiameter,
439 const double& dmt,
440 const double dH[nSpace],
441 double& cfl,
442 double& tau)
443 {
444 double h,nrm_v,oneByAbsdt;
445 h = elementDiameter;
446 nrm_v=0.0;
447 for(int I=0;I<nSpace;I++)
448 nrm_v+=dH[I]*dH[I];
449 nrm_v = sqrt(nrm_v);
450 cfl = nrm_v/h;
451 oneByAbsdt = fabs(dmt);
452 tau = 1.0/(2.0*nrm_v/h + oneByAbsdt + 1.0e-8);
453 }
454
455
456 inline
457 void calculateSubgridError_tau( const double& Ct_sge,
458 const double G[nSpace*nSpace],
459 const double& A0,
460 const double Ai[nSpace],
461 double& tau_v,
462 double& cfl)
463 {
464 double v_d_Gv=0.0;
465 for(int I=0;I<nSpace;I++)
466 for (int J=0;J<nSpace;J++)
467 v_d_Gv += Ai[I]*G[I*nSpace+J]*Ai[J];
468
469 tau_v = 1.0/sqrt(Ct_sge*A0*A0 + v_d_Gv);
470 }
471
472
473
474 inline
475 void calculateNumericalDiffusion(const double& shockCapturingDiffusion,
476 const double& elementDiameter,
477 const double& strong_residual,
478 const double grad_u[nSpace],
479 double& numDiff)
480 {
481 double h,
482 num,
483 den,
484 n_grad_u;
485 h = elementDiameter;
486 n_grad_u = 0.0;
487 for (int I=0;I<nSpace;I++)
488 n_grad_u += grad_u[I]*grad_u[I];
489 num = shockCapturingDiffusion*0.5*h*fabs(strong_residual);
490 den = sqrt(n_grad_u) + 1.0e-8;
491 numDiff = num/den;
492 }
493
494 inline
495 void exteriorNumericalAdvectiveFlux(const int& isDOFBoundary_u,
496 const int& isAdvectiveFluxBoundary_u,
497 const double n[nSpace],
498 const double& bc_u,
499 const double& bc_flux_u,
500 const double& u,
501 const double velocity[nSpace],
502 double& flux)
503 {
504
505 double flow=0.0;
506 for (int I=0; I < nSpace; I++)
507 flow += n[I]*velocity[I];
508 //std::cout<<" isDOFBoundary_u= "<<isDOFBoundary_u<<" flow= "<<flow<<std::endl;
509 if (isDOFBoundary_u == 1)
510 {
511 //std::cout<<"Dirichlet boundary u and bc_u "<<u<<'\t'<<bc_u<<std::endl;
512 if (flow >= 0.0)
513 {
514 flux = u*flow;
515 //flux = flow;
516 }
517 else
518 {
519 flux = bc_u*flow;
520 //flux = flow;
521 }
522 }
523 else if (isAdvectiveFluxBoundary_u == 1)
524 {
525 flux = bc_flux_u;
526 //std::cout<<"Flux boundary flux and flow"<<flux<<'\t'<<flow<<std::endl;
527 }
528 else
529 {
530 if (flow >= 0.0)
531 {
532 flux = u*flow;
533 }
534 else
535 {
536 //std::cout<<"warning: Kappa open boundary with no external trace, setting to zero for inflow n= ["<<n[0]<<","<<n[1]<<","<<n[2]<<"]"<<std::endl;
537 flux = 0.0;
538 }
539 //std::cout<<"No BC boundary flux and flow "<<flux<<'\t'<<flow<<std::endl;
540
541 }
542 //flux = flow;
543 //std::cout<<"flux error "<<flux-flow<<std::endl;
544 //std::cout<<"flux in computationa"<<flux<<std::endl;
545 }
546 inline double smoothedHeaviside(double eps, double phi)
547 {
548 double H;
549 if (phi > eps)
550 H=1.0;
551 else if (phi < -eps)
552 H=0.0;
553 else if (phi==0.0)
554 H=0.5;
555 else
556 H = 0.5*(1.0 + phi/eps + sin(M_PI*phi/eps)/M_PI);
557 return H;
558 }
559
560 inline
561 void exteriorNumericalAdvectiveFluxDerivative(const int& isDOFBoundary_u,
562 const int& isAdvectiveFluxBoundary,
563 const double n[nSpace],
564 const double velocity[nSpace],
565 double& dflux)
566 {
567 double flow=0.0;
568 for (int I=0; I < nSpace; I++)
569 flow += n[I]*velocity[I];
570 //double flow=n[0]*velocity[0]+n[1]*velocity[1]+n[2]*velocity[2];
571 dflux=0.0;//default to no flux
572 if (isDOFBoundary_u == 1)
573 {
574 if (flow >= 0.0)
575 {
576 dflux = flow;
577 }
578 else
579 {
580 dflux = 0.0;
581 }
582 }
583 else if (isAdvectiveFluxBoundary == 1)
584 {
585 dflux = 0.0;
586 }
587 else
588 {
589 if (flow >= 0.0)
590 {
591 dflux = flow;
592 }
593 }
594 }
595 inline
596 void exteriorNumericalDiffusiveFlux(const double& bc_flux,
597 const int& isDOFBoundary,
598 const int& isDiffusiveFluxBoundary,
599 double n[nSpace],
600 double bc_u,
601 double a,
602 double grad_psi[nSpace],
603 double u,
604 double penalty,
605 double& flux)
606 {
607 double v_I;
608 flux = 0.0;
609 if (isDiffusiveFluxBoundary)
610 {
611 flux = bc_flux;
612 }
613 else if (isDOFBoundary)
614 {
615 flux = 0.0;
616 for(int I=0;I<nSpace;I++)
617 {
618 v_I = -a*grad_psi[I];
619 flux += v_I*n[I];
620 }
621 flux += penalty*(u-bc_u);
622 }
623 else
624 {
625 //std::cerr<<"warning, Kappa diffusion term with no boundary condition set, setting diffusive flux to 0.0"<<std::endl;
626 flux = 0.0;
627 }
628 }
629 inline
630 void exteriorNumericalDiffusiveFluxDerivative(const int& isDOFBoundary,
631 const int& isDiffusiveFluxBoundary,
632 double n[nSpace],
633 double a,
634 double da,
635 double grad_psi[nSpace],
636 const double grad_v[nSpace],
637 double v,
638 double penalty,
639 double& fluxJacobian)
640 {
641 if (isDiffusiveFluxBoundary == 0 && isDOFBoundary == 1)
642 {
643 fluxJacobian = 0.0;
644 for(int I=0;I<nSpace;I++)
645 {
646 fluxJacobian -= (a*grad_v[I] + da*v*grad_psi[I])*n[I];
647 }
648 fluxJacobian += penalty*v;
649 }
650 else
651 fluxJacobian = 0.0;
652 }
653
654 void calculateResidual(arguments_dict& args)
655 {
656 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
657 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
658 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
659 xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
660 double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
661 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
662 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
663 xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
664 xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
665 xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
666 xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
667 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
668 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
669 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
670 xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
671 xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
672 xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
673 xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
674 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
675 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
676 int nElements_global = args.scalar<int>("nElements_global");
677 double nu_0 = args.scalar<double>("nu_0");
678 double nu_1 = args.scalar<double>("nu_1");
679 double sigma_k = args.scalar<double>("sigma_k");
680 double c_mu = args.scalar<double>("c_mu");
681 double rho_0 = args.scalar<double>("rho_0");
682 double rho_1 = args.scalar<double>("rho_1");
683 double sedFlag = args.scalar<int>("sedFlag");
684 xt::pyarray<double>& q_vos = args.array<double>("q_vos");
685 xt::pyarray<double>& q_vos_gradc = args.array<double>("q_vos_gradc");
686 xt::pyarray<double>& ebqe_q_vos = args.array<double>("ebqe_q_vos");
687 xt::pyarray<double>& ebqe_q_vos_gradc = args.array<double>("ebqe_q_vos_gradc");
688 double rho_f = args.scalar<double>("rho_f");
689 double rho_s = args.scalar<double>("rho_s");
690 xt::pyarray<double>& vs = args.array<double>("vs");
691 xt::pyarray<double>& ebqe_vs = args.array<double>("ebqe_vs");
692 xt::pyarray<double>& g = args.array<double>("g");
693 int dissipation_model_flag = args.scalar<int>("dissipation_model_flag");
694 double useMetrics = args.scalar<double>("useMetrics");
695 double alphaBDF = args.scalar<double>("alphaBDF");
696 int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
697 double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
698 double sc_uref = args.scalar<double>("sc_uref");
699 double sc_alpha = args.scalar<double>("sc_alpha");
700 xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
701 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
702 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
703 xt::pyarray<double>& u_dof_old = args.array<double>("u_dof_old");
704 xt::pyarray<double>& velocity = args.array<double>("velocity");
705 xt::pyarray<double>& phi_ls = args.array<double>("phi_ls");
706 xt::pyarray<double>& q_dissipation = args.array<double>("q_dissipation");
707 xt::pyarray<double>& q_grad_dissipation = args.array<double>("q_grad_dissipation");
708 xt::pyarray<double>& q_porosity = args.array<double>("q_porosity");
709 xt::pyarray<double>& velocity_dof_u = args.array<double>("velocity_dof_u");
710 xt::pyarray<double>& velocity_dof_v = args.array<double>("velocity_dof_v");
711 xt::pyarray<double>& velocity_dof_w = args.array<double>("velocity_dof_w");
712 xt::pyarray<double>& q_m = args.array<double>("q_m");
713 xt::pyarray<double>& q_u = args.array<double>("q_u");
714 xt::pyarray<double>& q_grad_u = args.array<double>("q_grad_u");
715 xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
716 xt::pyarray<double>& cfl = args.array<double>("cfl");
717 xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
718 xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
719 xt::pyarray<double>& ebqe_penalty_ext = args.array<double>("ebqe_penalty_ext");
720 int offset_u = args.scalar<int>("offset_u");
721 int stride_u = args.scalar<int>("stride_u");
722 xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
723 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
724 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
725 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
726 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
727 xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
728 xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
729 xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
730 xt::pyarray<int>& isAdvectiveFluxBoundary_u = args.array<int>("isAdvectiveFluxBoundary_u");
731 xt::pyarray<double>& ebqe_bc_advectiveFlux_u_ext = args.array<double>("ebqe_bc_advectiveFlux_u_ext");
732 xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.array<int>("isDiffusiveFluxBoundary_u");
733 xt::pyarray<double>& ebqe_bc_diffusiveFlux_u_ext = args.array<double>("ebqe_bc_diffusiveFlux_u_ext");
734 xt::pyarray<double>& ebqe_phi = args.array<double>("ebqe_phi");
735 double epsFact = args.scalar<double>("epsFact");
736 xt::pyarray<double>& ebqe_dissipation = args.array<double>("ebqe_dissipation");
737 xt::pyarray<double>& ebqe_porosity = args.array<double>("ebqe_porosity");
738 xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
739 xt::pyarray<double>& ebqe_flux = args.array<double>("ebqe_flux");
740 //cek should this be read in?
741 double Ct_sge = 4.0;
742 //loop over elements to compute volume integrals and load them into element and global residual
743 //
744 //eN is the element index
745 //eN_k is the quadrature point index for a scalar
746 //eN_k_nSpace is the quadrature point index for a vector
747 //eN_i is the element test function index
748 //eN_j is the element trial function index
749 //eN_k_j is the quadrature point index for a trial function
750 //eN_k_i is the quadrature point index for a trial function
751 for(int eN=0;eN<nElements_global;eN++)
752 {
753 //declare local storage for element residual and initialize
754 double elementResidual_u[nDOF_test_element];
755 for (int i=0;i<nDOF_test_element;i++)
756 {
757 elementResidual_u[i]=0.0;
758 }//i
759 //loop over quadrature points and compute integrands
760 for (int k=0;k<nQuadraturePoints_element;k++)
761 {
762 //compute indeces and declare local storage
763 int eN_k = eN*nQuadraturePoints_element+k,
764 eN_k_nSpace = eN_k*nSpace,
765 eN_nDOF_trial_element = eN*nDOF_trial_element;
766 double u=0.0,u_old=0.0,grad_u[nSpace],grad_u_old[nSpace],
767 m=0.0,dm=0.0,
768 f[nSpace],df[nSpace],df_minus_da_grad_u[nSpace],
769 m_t=0.0,dm_t=0.0,
770 a=0.0,da=0.0,
771 r=0.0,dr=0.0,
772 grad_vx[nSpace],grad_vy[nSpace],grad_vz[nSpace],
773 pdeResidual_u=0.0,
774 Lstar_u[nDOF_test_element],
775 subgridError_u=0.0,
776 tau=0.0,tau0=0.0,tau1=0.0,
777 numDiff0=0.0,numDiff1=0.0,
778 jac[nSpace*nSpace],
779 jacDet,
780 jacInv[nSpace*nSpace],
781 u_grad_trial[nDOF_trial_element*nSpace],
782 u_test_dV[nDOF_trial_element],
783 u_grad_test_dV[nDOF_test_element*nSpace],
784 dV,x,y,z,xt,yt,zt,
785 G[nSpace*nSpace],G_dd_G,tr_G,norm_Rv;
786 // //
787 // //compute solution and gradients at quadrature points
788 // //
789 // u=0.0;
790 // for (int I=0;I<nSpace;I++)
791 // {
792 // grad_u[I]=0.0;
793 // }
794 // for (int j=0;j<nDOF_trial_element;j++)
795 // {
796 // int eN_j=eN*nDOF_trial_element+j;
797 // int eN_k_j=eN_k*nDOF_trial_element+j;
798 // int eN_k_j_nSpace = eN_k_j*nSpace;
799 // u += valFromDOF_c(u_dof[u_l2g[eN_j]],u_trial[eN_k_j]);
800 // for (int I=0;I<nSpace;I++)
801 // {
802 // grad_u[I] += gradFromDOF_c(u_dof[u_l2g[eN_j]],u_grad_trial[eN_k_j_nSpace+I]);
803 // }
804 // }
805 ck.calculateMapping_element(eN,
806 k,
807 mesh_dof.data(),
808 mesh_l2g.data(),
809 mesh_trial_ref.data(),
810 mesh_grad_trial_ref.data(),
811 jac,
812 jacDet,
813 jacInv,
814 x,y,z);
815 ck.calculateMappingVelocity_element(eN,
816 k,
817 mesh_velocity_dof.data(),
818 mesh_l2g.data(),
819 mesh_trial_ref.data(),
820 xt,yt,zt);
821 //get the physical integration weight
822 dV = fabs(jacDet)*dV_ref[k];
823 ck.calculateG(jacInv,G,G_dd_G,tr_G);
824 //get the trial function gradients
825 ck.gradTrialFromRef(&u_grad_trial_ref[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
826 //get the solution and previous solution
827 ck.valFromDOF(u_dof.data(),&u_l2g[eN_nDOF_trial_element],&u_trial_ref[k*nDOF_trial_element],u);
828 ck.valFromDOF(u_dof_old.data(),&u_l2g[eN_nDOF_trial_element],&u_trial_ref[k*nDOF_trial_element],u_old);
829
830 //get the solution gradients
831 ck.gradFromDOF(u_dof.data(),&u_l2g[eN_nDOF_trial_element],u_grad_trial,grad_u);
832 ck.gradFromDOF(u_dof_old.data(),&u_l2g[eN_nDOF_trial_element],u_grad_trial,grad_u_old);
833 //
834 //compute velocity production terms, ***assumes same spaces for velocity dofs and kappa!***
835 ck.gradFromDOF(velocity_dof_u.data(),&u_l2g[eN_nDOF_trial_element],u_grad_trial,grad_vx);
836 ck.gradFromDOF(velocity_dof_v.data(),&u_l2g[eN_nDOF_trial_element],u_grad_trial,grad_vy);
837 ck.gradFromDOF(velocity_dof_w.data(),&u_l2g[eN_nDOF_trial_element],u_grad_trial,grad_vz);
838 //
839
840 //
841 //precalculate test function products with integration weights
842 for (int j=0;j<nDOF_trial_element;j++)
843 {
844 u_test_dV[j] = u_test_ref[k*nDOF_trial_element+j]*dV;
845 for (int I=0;I<nSpace;I++)
846 {
847 u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
848 }
849 }
850 //
851 //calculate pde coefficients at quadrature points
852 //
853 evaluateCoefficients(&velocity[eN_k_nSpace],
854 epsFact,
855 phi_ls[eN_k],
856 nu_0,
857 nu_1,
858 sigma_k,
859 c_mu,
860 grad_vx,
861 grad_vy,
862 grad_vz,
863 u,
864 u_old,
865 q_dissipation[eN_k],
866 q_porosity[eN_k],
867// Argumentlist for sediment
868 sedFlag,
869 q_vos[eN_k],
870 &q_vos_gradc[eN_k_nSpace],
871 rho_f,
872 rho_s,
873 &vs[eN_k_nSpace],
874 &g[0],
875 //end sediment
876 dissipation_model_flag,
877 grad_u_old,
878 &q_grad_dissipation[eN_k_nSpace],
879 m,
880 dm,
881 f,
882 df,
883 a,
884 da,
885 r,
886 dr);
887 //
888 //moving mesh
889 //
890 f[0] -= MOVING_DOMAIN*m*xt;
891 f[1] -= MOVING_DOMAIN*m*yt;
892 f[2] -= MOVING_DOMAIN*m*zt;
893 df[0] -= MOVING_DOMAIN*dm*xt;
894 df[1] -= MOVING_DOMAIN*dm*yt;
895 df[2] -= MOVING_DOMAIN*dm*zt;
896
897 //combine df and da/du \grad u term for stabilization and jacobian calculations
898 df_minus_da_grad_u[0] = df[0] - da*grad_u[0];
899 df_minus_da_grad_u[1] = df[1] - da*grad_u[1];
900 df_minus_da_grad_u[2] = df[2] - da*grad_u[2];
901 //
902 //calculate time derivative at quadrature points
903 //
904 ck.bdf(alphaBDF,
905 q_m_betaBDF[eN_k],
906 m,
907 dm,
908 m_t,
909 dm_t);
910 //
911 //calculate subgrid error (strong residual and adjoint)
912 //
913 //calculate strong residual
914 pdeResidual_u = ck.Mass_strong(m_t) + ck.Advection_strong(df_minus_da_grad_u,grad_u)
915 + ck.Reaction_strong(r);
916 //calculate adjoint
917 for (int i=0;i<nDOF_test_element;i++)
918 {
919 // int eN_k_i_nSpace = (eN_k*nDOF_trial_element+i)*nSpace;
920 // Lstar_u[i] = ck.Advection_adjoint(df,&u_grad_test_dV[eN_k_i_nSpace]);
921 int i_nSpace = i*nSpace;
922 Lstar_u[i] = ck.Advection_adjoint(df_minus_da_grad_u,&u_grad_test_dV[i_nSpace]) +
923 ck.Reaction_adjoint(dr,u_test_dV[i]);
924 }
925 //calculate tau and tau*Res
926 calculateSubgridError_tau(elementDiameter[eN],dm_t + dr,df_minus_da_grad_u,cfl[eN_k],tau0);
927 calculateSubgridError_tau(Ct_sge,
928 G,
929 dm_t + dr,
930 df_minus_da_grad_u,
931 tau1,
932 cfl[eN_k]);
933
934 tau = useMetrics*tau1+(1.0-useMetrics)*tau0;
935
936 subgridError_u = -tau*pdeResidual_u;
937 //
938 //calculate shock capturing diffusion
939 //
940
941
942 ck.calculateNumericalDiffusion(shockCapturingDiffusion,elementDiameter[eN],pdeResidual_u,grad_u,numDiff0);
943 ck.calculateNumericalDiffusion(shockCapturingDiffusion,sc_uref, sc_alpha,G,G_dd_G,pdeResidual_u,grad_u,numDiff1);
944 q_numDiff_u[eN_k] = useMetrics*numDiff1+(1.0-useMetrics)*numDiff0;
945 //std::cout<<tau<<" "<<q_numDiff_u[eN_k]<<std::endl;
946 //
947 //update element residual
948 //
949 for(int i=0;i<nDOF_test_element;i++)
950 {
951 int eN_k_i=eN_k*nDOF_test_element+i,
952 eN_k_i_nSpace = eN_k_i*nSpace,
953 i_nSpace=i*nSpace;
954
955 elementResidual_u[i] += ck.Mass_weak(m_t,u_test_dV[i]) +
956 ck.Advection_weak(f,&u_grad_test_dV[i_nSpace]) +
957 ck.SubgridError(subgridError_u,Lstar_u[i]) +
958 ck.NumericalDiffusion(a,grad_u,&u_grad_test_dV[i_nSpace]) + //scalar diffusion so steal numericalDiffusion approximation
959 ck.NumericalDiffusion(q_numDiff_u_last[eN_k],grad_u,&u_grad_test_dV[i_nSpace]) +
960 ck.Reaction_weak(r,u_test_dV[i]);
961
962 }//i
963 //
964 //cek/ido todo, get rid of m, since u=m
965 //save momentum for time history and velocity for subgrid error
966 //save solution for other models
967 //
968 q_u[eN_k] = u;
969 q_m[eN_k] = m;
970 for (int I=0; I < nSpace; I++)
971 {
972 q_grad_u[eN_k_nSpace + I] = grad_u[I];
973 }
974 }
975 //
976 //load element into global residual and save element residual
977 //
978 for(int i=0;i<nDOF_test_element;i++)
979 {
980 int eN_i=eN*nDOF_test_element+i;
981
982 globalResidual[offset_u+stride_u*u_l2g[eN_i]] += elementResidual_u[i];
983 }//i
984 }//elements
985 //
986 //loop over exterior element boundaries to calculate surface integrals and load into element and global residuals
987 //
988 //ebNE is the Exterior element boundary INdex
989 //ebN is the element boundary INdex
990 //eN is the element index
991 for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
992 {
993 int ebN = exteriorElementBoundariesArray[ebNE],
994 eN = elementBoundaryElementsArray[ebN*2+0],
995 ebN_local = elementBoundaryLocalElementBoundariesArray[ebN*2+0],
996 eN_nDOF_trial_element = eN*nDOF_trial_element;
997 double elementResidual_u[nDOF_test_element];
998 for (int i=0;i<nDOF_test_element;i++)
999 {
1000 elementResidual_u[i]=0.0;
1001 }
1002 for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
1003 {
1004 int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
1005 ebNE_kb_nSpace = ebNE_kb*nSpace,
1006 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
1007 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
1008 double u_ext=0.0,u_old_ext,
1009 grad_u_ext[nSpace],grad_vx_ext[nSpace],grad_vy_ext[nSpace],grad_vz_ext[nSpace],
1010 grad_u_old_ext[nSpace],grad_dissipation_ext_dummy[nSpace],
1011 m_ext=0.0,
1012 dm_ext=0.0,
1013 f_ext[nSpace],
1014 df_ext[nSpace],
1015 a_ext=0.0,da_ext=0.0,
1016 r_ext=0.0,dr_ext=0.0,
1017 flux_ext=0.0,
1018 diffusive_flux_ext=0.0,
1019 bc_u_ext=0.0,
1020 bc_grad_u_ext[nSpace],
1021 bc_m_ext=0.0,
1022 bc_dm_ext=0.0,
1023 bc_f_ext[nSpace],
1024 bc_df_ext[nSpace],
1025 jac_ext[nSpace*nSpace],
1026 jacDet_ext,
1027 jacInv_ext[nSpace*nSpace],
1028 boundaryJac[nSpace*(nSpace-1)],
1029 metricTensor[(nSpace-1)*(nSpace-1)],
1030 metricTensorDetSqrt,
1031 dS,
1032 u_test_dS[nDOF_test_element],
1033 u_grad_trial_trace[nDOF_trial_element*nSpace],
1034 normal[3],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
1035 G[nSpace*nSpace],G_dd_G,tr_G;
1036 //
1037 //calculate the solution and gradients at quadrature points
1038 //
1039 //compute information about mapping from reference element to physical element
1040 ck.calculateMapping_elementBoundary(eN,
1041 ebN_local,
1042 kb,
1043 ebN_local_kb,
1044 mesh_dof.data(),
1045 mesh_l2g.data(),
1046 mesh_trial_trace_ref.data(),
1047 mesh_grad_trial_trace_ref.data(),
1048 boundaryJac_ref.data(),
1049 jac_ext,
1050 jacDet_ext,
1051 jacInv_ext,
1052 boundaryJac,
1053 metricTensor,
1054 metricTensorDetSqrt,
1055 normal_ref.data(),
1056 normal,
1057 x_ext,y_ext,z_ext);
1058 ck.calculateMappingVelocity_elementBoundary(eN,
1059 ebN_local,
1060 kb,
1061 ebN_local_kb,
1062 mesh_velocity_dof.data(),
1063 mesh_l2g.data(),
1064 mesh_trial_trace_ref.data(),
1065 xt_ext,yt_ext,zt_ext,
1066 normal,
1067 boundaryJac,
1068 metricTensor,
1069 integralScaling);
1070 dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref[kb];
1071 //get the metric tensor
1072 //cek todo use symmetry
1073 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
1074 //compute shape and solution information
1075 //shape
1076 ck.gradTrialFromRef(&u_grad_trial_trace_ref[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,u_grad_trial_trace);
1077 //solution and gradients
1078 ck.valFromDOF(u_dof.data(),&u_l2g[eN_nDOF_trial_element],&u_trial_trace_ref[ebN_local_kb*nDOF_test_element],u_ext);
1079 ck.valFromDOF(u_dof_old.data(),&u_l2g[eN_nDOF_trial_element],&u_trial_trace_ref[ebN_local_kb*nDOF_test_element],u_old_ext);
1080 ck.gradFromDOF(u_dof.data(),&u_l2g[eN_nDOF_trial_element],u_grad_trial_trace,grad_u_ext);
1081 ck.gradFromDOF(u_dof_old.data(),&u_l2g[eN_nDOF_trial_element],u_grad_trial_trace,grad_u_old_ext);
1082
1083 //mwf hack, skip on boundary for now
1084 grad_dissipation_ext_dummy[0] = 0.0; grad_dissipation_ext_dummy[1] = 0.0; grad_dissipation_ext_dummy[2] = 0.0;
1085 //
1086 //compute velocity production terms, ***assumes same spaces for velocity dofs and kappa!***
1087 ck.gradFromDOF(velocity_dof_u.data(),&u_l2g[eN_nDOF_trial_element],u_grad_trial_trace,grad_vx_ext);
1088 ck.gradFromDOF(velocity_dof_v.data(),&u_l2g[eN_nDOF_trial_element],u_grad_trial_trace,grad_vy_ext);
1089 ck.gradFromDOF(velocity_dof_w.data(),&u_l2g[eN_nDOF_trial_element],u_grad_trial_trace,grad_vz_ext);
1090 //
1091
1092 //precalculate test function products with integration weights
1093 for (int j=0;j<nDOF_trial_element;j++)
1094 {
1095 u_test_dS[j] = u_test_trace_ref[ebN_local_kb*nDOF_test_element+j]*dS;
1096 }
1097 //
1098 //load the boundary values
1099 //
1100 bc_u_ext = isDOFBoundary_u[ebNE_kb]*ebqe_bc_u_ext[ebNE_kb]+(1-isDOFBoundary_u[ebNE_kb])*u_ext;
1101 //
1102 //calculate the pde coefficients using the solution and the boundary values for the solution
1103 //
1104 evaluateCoefficients(&ebqe_velocity_ext[ebNE_kb_nSpace],
1105 epsFact,
1106 ebqe_phi[ebNE_kb],
1107 nu_0,
1108 nu_1,
1109 sigma_k,
1110 c_mu,
1111 grad_vx_ext,
1112 grad_vy_ext,
1113 grad_vz_ext,
1114 u_ext,
1115 u_old_ext,
1116 ebqe_dissipation[ebNE_kb],
1117 ebqe_porosity[ebNE_kb],
1118// Argumentlist for sediment
1119 sedFlag,
1120 ebqe_q_vos[ebNE_kb],
1121 &ebqe_q_vos_gradc[ebNE_kb_nSpace],
1122 rho_f,
1123 rho_s,
1124 &ebqe_vs[ebNE_kb_nSpace],
1125 &g[0],
1126 //end sediment
1127 dissipation_model_flag,
1128 grad_u_old_ext,
1129 grad_dissipation_ext_dummy,
1130 m_ext,
1131 dm_ext,
1132 f_ext,
1133 df_ext,
1134 a_ext,
1135 da_ext,
1136 r_ext,
1137 dr_ext);
1138 evaluateCoefficients(&ebqe_velocity_ext[ebNE_kb_nSpace],
1139 epsFact,
1140 ebqe_phi[ebNE_kb],
1141 nu_0,
1142 nu_1,
1143 sigma_k,
1144 c_mu,
1145 grad_vx_ext,
1146 grad_vy_ext,
1147 grad_vz_ext,
1148 bc_u_ext,
1149 bc_u_ext,
1150 ebqe_dissipation[ebNE_kb],
1151 ebqe_porosity[ebNE_kb],
1152// Argumentlist for sediment
1153 sedFlag,
1154 ebqe_q_vos[ebNE_kb],
1155 &ebqe_q_vos_gradc[ebNE_kb_nSpace],
1156 rho_f,
1157 rho_s,
1158 &ebqe_vs[ebNE_kb_nSpace],
1159 &g[0],
1160 //end sediment
1161 dissipation_model_flag,
1162 grad_u_old_ext,
1163 grad_dissipation_ext_dummy,
1164 bc_m_ext,
1165 bc_dm_ext,
1166 bc_f_ext,
1167 bc_df_ext,
1168 a_ext,
1169 da_ext,
1170 r_ext,
1171 dr_ext);
1172 //
1173 //moving mesh
1174 //
1175 double velocity_ext[nSpace];
1176 velocity_ext[0] = ebqe_velocity_ext[ebNE_kb_nSpace+0] - MOVING_DOMAIN*xt_ext;
1177 velocity_ext[1] = ebqe_velocity_ext[ebNE_kb_nSpace+1] - MOVING_DOMAIN*yt_ext;
1178 velocity_ext[2] = ebqe_velocity_ext[ebNE_kb_nSpace+2] - MOVING_DOMAIN*zt_ext;
1179 //
1180 //calculate the numerical fluxes
1181 //
1182 exteriorNumericalAdvectiveFlux(isDOFBoundary_u[ebNE_kb],
1183 isAdvectiveFluxBoundary_u[ebNE_kb],
1184 normal,
1185 bc_u_ext,
1186 ebqe_bc_advectiveFlux_u_ext[ebNE_kb],
1187 u_ext,//smoothedHeaviside(eps,ebqe_phi[ebNE_kb]),
1188 velocity_ext,
1189 flux_ext);
1190 //diffusive flux now as well
1191 //for now just apply flux boundary through advection term
1192 const double bc_diffusive_flux = ebqe_bc_diffusiveFlux_u_ext[ebNE_kb];
1193 exteriorNumericalDiffusiveFlux(bc_diffusive_flux,
1194 isDOFBoundary_u[ebNE_kb],
1195 isDiffusiveFluxBoundary_u[ebNE_kb],
1196 normal,
1197 bc_u_ext,
1198 a_ext,
1199 grad_u_ext,
1200 u_ext,
1201 ebqe_penalty_ext[ebNE_kb],//penalty,
1202 diffusive_flux_ext);
1203 //mwf debug
1204 //std::cout<<"Residual ebNE= "<<ebNE<<" kb= "<<kb <<" penalty= "<<ebqe_penalty_ext[ebNE_kb] <<std::endl;
1205 flux_ext += diffusive_flux_ext;
1206 ebqe_flux[ebNE_kb] = flux_ext;
1207 //save for other models? cek need to be consistent with numerical flux
1208 if(flux_ext >=0.0)
1209 ebqe_u[ebNE_kb] = u_ext;
1210 else
1211 ebqe_u[ebNE_kb] = bc_u_ext;
1212 //
1213 //update residuals
1214 //
1215 for (int i=0;i<nDOF_test_element;i++)
1216 {
1217 int ebNE_kb_i = ebNE_kb*nDOF_test_element+i;
1218
1219 elementResidual_u[i] += ck.ExteriorElementBoundaryFlux(flux_ext,u_test_dS[i]);
1220 }//i
1221 }//kb
1222 //
1223 //update the element and global residual storage
1224 //
1225 for (int i=0;i<nDOF_test_element;i++)
1226 {
1227 int eN_i = eN*nDOF_test_element+i;
1228
1229 globalResidual[offset_u+stride_u*u_l2g[eN_i]] += elementResidual_u[i];
1230 }//i
1231 }//ebNE
1232 }
1233
1234 void calculateJacobian(arguments_dict& args) //VRANS
1235 {
1236 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1237 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1238 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1239 xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
1240 double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
1241 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1242 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1243 xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1244 xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1245 xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1246 xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
1247 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
1248 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
1249 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
1250 xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
1251 xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
1252 xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
1253 xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
1254 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
1255 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
1256 int nElements_global = args.scalar<int>("nElements_global");
1257 double nu_0 = args.scalar<double>("nu_0");
1258 double nu_1 = args.scalar<double>("nu_1");
1259 double sigma_k = args.scalar<double>("sigma_k");
1260 double c_mu = args.scalar<double>("c_mu");
1261 double rho_0 = args.scalar<double>("rho_0");
1262 double rho_1 = args.scalar<double>("rho_1");
1263 int dissipation_model_flag = args.scalar<int>("dissipation_model_flag");
1264 double useMetrics = args.scalar<double>("useMetrics");
1265 double alphaBDF = args.scalar<double>("alphaBDF");
1266 int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
1267 double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
1268 xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1269 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1270 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1271 xt::pyarray<double>& u_dof_old = args.array<double>("u_dof_old");
1272 xt::pyarray<double>& velocity = args.array<double>("velocity");
1273 xt::pyarray<double>& phi_ls = args.array<double>("phi_ls");
1274 xt::pyarray<double>& q_dissipation = args.array<double>("q_dissipation");
1275 xt::pyarray<double>& q_grad_dissipation = args.array<double>("q_grad_dissipation");
1276 xt::pyarray<double>& q_porosity = args.array<double>("q_porosity");
1277 double sedFlag = args.scalar<int>("sedFlag");
1278 xt::pyarray<double>& q_vos = args.array<double>("q_vos");
1279 xt::pyarray<double>& q_vos_gradc = args.array<double>("q_vos_gradc");
1280 xt::pyarray<double>& ebqe_q_vos = args.array<double>("ebqe_q_vos");
1281 xt::pyarray<double>& ebqe_q_vos_gradc = args.array<double>("ebqe_q_vos_gradc");
1282 double rho_f = args.scalar<double>("rho_f");
1283 double rho_s = args.scalar<double>("rho_s");
1284 xt::pyarray<double>& vs = args.array<double>("vs");
1285 xt::pyarray<double>& ebqe_vs = args.array<double>("ebqe_vs");
1286 xt::pyarray<double>& g = args.array<double>("g");
1287 xt::pyarray<double>& velocity_dof_u = args.array<double>("velocity_dof_u");
1288 xt::pyarray<double>& velocity_dof_v = args.array<double>("velocity_dof_v");
1289 xt::pyarray<double>& velocity_dof_w = args.array<double>("velocity_dof_w");
1290 xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
1291 xt::pyarray<double>& cfl = args.array<double>("cfl");
1292 xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
1293 xt::pyarray<double>& ebqe_penalty_ext = args.array<double>("ebqe_penalty_ext");
1294 xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
1295 xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
1296 xt::pyarray<double>& globalJacobian = args.array<double>("globalJacobian");
1297 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
1298 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
1299 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
1300 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
1301 xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
1302 xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
1303 xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
1304 xt::pyarray<int>& isAdvectiveFluxBoundary_u = args.array<int>("isAdvectiveFluxBoundary_u");
1305 xt::pyarray<double>& ebqe_bc_advectiveFlux_u_ext = args.array<double>("ebqe_bc_advectiveFlux_u_ext");
1306 xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.array<int>("isDiffusiveFluxBoundary_u");
1307 xt::pyarray<double>& ebqe_bc_diffusiveFlux_u_ext = args.array<double>("ebqe_bc_diffusiveFlux_u_ext");
1308 xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
1309 xt::pyarray<double>& ebqe_phi = args.array<double>("ebqe_phi");
1310 double epsFact = args.scalar<double>("epsFact");
1311 xt::pyarray<double>& ebqe_dissipation = args.array<double>("ebqe_dissipation");
1312 xt::pyarray<double>& ebqe_porosity = args.array<double>("ebqe_porosity");
1313 double Ct_sge = 4.0;
1314
1315 //
1316 //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
1317 //
1318 for(int eN=0;eN<nElements_global;eN++)
1319 {
1320 double elementJacobian_u_u[nDOF_test_element][nDOF_trial_element];
1321 for (int i=0;i<nDOF_test_element;i++)
1322 for (int j=0;j<nDOF_trial_element;j++)
1323 {
1324 elementJacobian_u_u[i][j]=0.0;
1325 }
1326 for (int k=0;k<nQuadraturePoints_element;k++)
1327 {
1328 int eN_k = eN*nQuadraturePoints_element+k, //index to a scalar at a quadrature point
1329 eN_k_nSpace = eN_k*nSpace,
1330 eN_nDOF_trial_element = eN*nDOF_trial_element; //index to a vector at a quadrature point
1331
1332 //declare local storage
1333 double u=0.0,u_old,
1334 grad_u[nSpace],grad_vx[nSpace],grad_vy[nSpace],grad_vz[nSpace],grad_u_old[nSpace],
1335 m=0.0,dm=0.0,a=0.0,da=0.0,r=0.0,dr=0.0,
1336 f[nSpace],df[nSpace],df_minus_da_grad_u[nSpace],
1337 m_t=0.0,dm_t=0.0,
1338 dpdeResidual_u_u[nDOF_trial_element],
1339 Lstar_u[nDOF_test_element],
1340 dsubgridError_u_u[nDOF_trial_element],
1341 tau=0.0,tau0=0.0,tau1=0.0,
1342 jac[nSpace*nSpace],
1343 jacDet,
1344 jacInv[nSpace*nSpace],
1345 u_grad_trial[nDOF_trial_element*nSpace],
1346 dV,
1347 u_test_dV[nDOF_test_element],
1348 u_grad_test_dV[nDOF_test_element*nSpace],
1349 x,y,z,xt,yt,zt,
1350 G[nSpace*nSpace],G_dd_G,tr_G;
1351 //
1352 //calculate solution and gradients at quadrature points
1353 //
1354 // u=0.0;
1355 // for (int I=0;I<nSpace;I++)
1356 // {
1357 // grad_u[I]=0.0;
1358 // }
1359 // for (int j=0;j<nDOF_trial_element;j++)
1360 // {
1361 // int eN_j=eN*nDOF_trial_element+j;
1362 // int eN_k_j=eN_k*nDOF_trial_element+j;
1363 // int eN_k_j_nSpace = eN_k_j*nSpace;
1364
1365 // u += valFromDOF_c(u_dof[u_l2g[eN_j]],u_trial[eN_k_j]);
1366 // for (int I=0;I<nSpace;I++)
1367 // {
1368 // grad_u[I] += gradFromDOF_c(u_dof[u_l2g[eN_j]],u_grad_trial[eN_k_j_nSpace+I]);
1369 // }
1370 // }
1371 //get jacobian, etc for mapping reference element
1372 ck.calculateMapping_element(eN,
1373 k,
1374 mesh_dof.data(),
1375 mesh_l2g.data(),
1376 mesh_trial_ref.data(),
1377 mesh_grad_trial_ref.data(),
1378 jac,
1379 jacDet,
1380 jacInv,
1381 x,y,z);
1382 ck.calculateMappingVelocity_element(eN,
1383 k,
1384 mesh_velocity_dof.data(),
1385 mesh_l2g.data(),
1386 mesh_trial_ref.data(),
1387 xt,yt,zt);
1388 //get the physical integration weight
1389 dV = fabs(jacDet)*dV_ref[k];
1390 ck.calculateG(jacInv,G,G_dd_G,tr_G);
1391 //get the trial function gradients
1392 ck.gradTrialFromRef(&u_grad_trial_ref[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
1393 //get the solution at old and new time step
1394 ck.valFromDOF(u_dof.data(),&u_l2g[eN_nDOF_trial_element],&u_trial_ref[k*nDOF_trial_element],u);
1395 ck.valFromDOF(u_dof_old.data(),&u_l2g[eN_nDOF_trial_element],&u_trial_ref[k*nDOF_trial_element],u_old);
1396 //get the solution gradients
1397 ck.gradFromDOF(u_dof.data(),&u_l2g[eN_nDOF_trial_element],u_grad_trial,grad_u);
1398 //
1399 //compute velocity production terms, ***assumes same spaces for velocity dofs and kappa!***
1400 ck.gradFromDOF(velocity_dof_u.data(),&u_l2g[eN_nDOF_trial_element],u_grad_trial,grad_vx);
1401 ck.gradFromDOF(velocity_dof_v.data(),&u_l2g[eN_nDOF_trial_element],u_grad_trial,grad_vy);
1402 ck.gradFromDOF(velocity_dof_w.data(),&u_l2g[eN_nDOF_trial_element],u_grad_trial,grad_vz);
1403 //
1404
1405 //precalculate test function products with integration weights
1406 for (int j=0;j<nDOF_trial_element;j++)
1407 {
1408 u_test_dV[j] = u_test_ref[k*nDOF_trial_element+j]*dV;
1409 for (int I=0;I<nSpace;I++)
1410 {
1411 u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
1412 }
1413 }
1414 //
1415 //calculate pde coefficients and derivatives at quadrature points
1416 //
1417 evaluateCoefficients(&velocity[eN_k_nSpace],
1418 epsFact,
1419 phi_ls[eN_k],
1420 nu_0,
1421 nu_1,
1422 sigma_k,
1423 c_mu,
1424 grad_vx,
1425 grad_vy,
1426 grad_vz,
1427 u,
1428 u_old,
1429 q_dissipation[eN_k],
1430 q_porosity[eN_k],
1431// Argumentlist for sediment
1432 sedFlag,
1433 q_vos[eN_k],
1434 &q_vos_gradc[eN_k_nSpace],
1435 rho_f,
1436 rho_s,
1437 &vs[eN_k_nSpace],
1438 g.data(),
1439 //end sediment
1440 dissipation_model_flag,
1441 grad_u_old,
1442 &q_grad_dissipation[eN_k_nSpace],
1443 m,
1444 dm,
1445 f,
1446 df,
1447 a,
1448 da,
1449 r,
1450 dr);
1451 //
1452 //moving mesh
1453 //
1454 f[0] -= MOVING_DOMAIN*m*xt;
1455 f[1] -= MOVING_DOMAIN*m*yt;
1456 f[2] -= MOVING_DOMAIN*m*zt;
1457 df[0] -= MOVING_DOMAIN*dm*xt;
1458 df[1] -= MOVING_DOMAIN*dm*yt;
1459 df[2] -= MOVING_DOMAIN*dm*zt;
1460
1461 //
1462 //account for nonlinearity in diffusion coefficient by piggy backing on advection routines
1463 //
1464 df_minus_da_grad_u[0] = df[0] - da*grad_u[0];
1465 df_minus_da_grad_u[1] = df[1] - da*grad_u[1];
1466 df_minus_da_grad_u[2] = df[2] - da*grad_u[2];
1467
1468 //
1469 //calculate time derivatives
1470 //
1471 ck.bdf(alphaBDF,
1472 q_m_betaBDF[eN_k],
1473 m,
1474 dm,
1475 m_t,
1476 dm_t);
1477 //
1478 //calculate subgrid error contribution to the Jacobian (strong residual, adjoint, jacobian of strong residual)
1479 //
1480 //calculate the adjoint times the test functions
1481 for (int i=0;i<nDOF_test_element;i++)
1482 {
1483 // int eN_k_i_nSpace = (eN_k*nDOF_trial_element+i)*nSpace;
1484 // Lstar_u[i]=ck.Advection_adjoint(df,&u_grad_test_dV[eN_k_i_nSpace]);
1485 int i_nSpace = i*nSpace;
1486 Lstar_u[i]=ck.Advection_adjoint(df_minus_da_grad_u,&u_grad_test_dV[i_nSpace])
1487 + ck.Reaction_adjoint(dr,u_test_dV[i]);
1488 }
1489 //calculate the Jacobian of strong residual
1490 for (int j=0;j<nDOF_trial_element;j++)
1491 {
1492 //int eN_k_j=eN_k*nDOF_trial_element+j;
1493 //int eN_k_j_nSpace = eN_k_j*nSpace;
1494 int j_nSpace = j*nSpace;
1495 dpdeResidual_u_u[j]= ck.MassJacobian_strong(dm_t,u_trial_ref[k*nDOF_trial_element+j]) +
1496 ck.AdvectionJacobian_strong(df_minus_da_grad_u,&u_grad_trial[j_nSpace]) +
1497 ck.ReactionJacobian_strong(dr,u_trial_ref[k*nDOF_trial_element+j]);
1498 }
1499 //tau and tau*Res
1500 calculateSubgridError_tau(elementDiameter[eN],
1501 dm_t + dr,
1502 df_minus_da_grad_u,
1503 cfl[eN_k],
1504 tau0);
1505
1506 calculateSubgridError_tau(Ct_sge,
1507 G,
1508 dm_t + dr,
1509 df_minus_da_grad_u,
1510 tau1,
1511 cfl[eN_k]);
1512 tau = useMetrics*tau1+(1.0-useMetrics)*tau0;
1513
1514 for(int j=0;j<nDOF_trial_element;j++)
1515 dsubgridError_u_u[j] = -tau*dpdeResidual_u_u[j];
1516 for(int i=0;i<nDOF_test_element;i++)
1517 {
1518 //int eN_k_i=eN_k*nDOF_test_element+i;
1519 //int eN_k_i_nSpace=eN_k_i*nSpace;
1520 for(int j=0;j<nDOF_trial_element;j++)
1521 {
1522 //int eN_k_j=eN_k*nDOF_trial_element+j;
1523 //int eN_k_j_nSpace = eN_k_j*nSpace;
1524 int j_nSpace = j*nSpace;
1525 int i_nSpace = i*nSpace;
1526 elementJacobian_u_u[i][j] += ck.MassJacobian_weak(dm_t,u_trial_ref[k*nDOF_trial_element+j],u_test_dV[i]) +
1527 ck.AdvectionJacobian_weak(df_minus_da_grad_u,u_trial_ref[k*nDOF_trial_element+j],&u_grad_test_dV[i_nSpace]) +
1528 ck.SubgridErrorJacobian(dsubgridError_u_u[j],Lstar_u[i]) +
1529 ck.NumericalDiffusionJacobian(a,&u_grad_trial[j_nSpace],&u_grad_test_dV[i_nSpace]) + //steal numericalDiffusion for scalar term
1530 ck.NumericalDiffusionJacobian(q_numDiff_u_last[eN_k],&u_grad_trial[j_nSpace],&u_grad_test_dV[i_nSpace]) +
1531 ck.ReactionJacobian_weak(dr,u_trial_ref[k*nDOF_trial_element+j],u_test_dV[i]);
1532 }//j
1533 }//i
1534 }//k
1535 //
1536 //load into element Jacobian into global Jacobian
1537 //
1538 for (int i=0;i<nDOF_test_element;i++)
1539 {
1540 int eN_i = eN*nDOF_test_element+i;
1541 for (int j=0;j<nDOF_trial_element;j++)
1542 {
1543 int eN_i_j = eN_i*nDOF_trial_element+j;
1544 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_u_u[eN_i_j]] += elementJacobian_u_u[i][j];
1545 }//j
1546 }//i
1547 }//elements
1548 //
1549 //loop over exterior element boundaries to compute the surface integrals and load them into the global Jacobian
1550 //
1551 for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
1552 {
1553 int ebN = exteriorElementBoundariesArray[ebNE];
1554 int eN = elementBoundaryElementsArray[ebN*2+0],
1555 ebN_local = elementBoundaryLocalElementBoundariesArray[ebN*2+0],
1556 eN_nDOF_trial_element = eN*nDOF_trial_element;
1557 for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
1558 {
1559 int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
1560 ebNE_kb_nSpace = ebNE_kb*nSpace,
1561 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
1562 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
1563
1564 double u_ext=0.0,u_old_ext=0.0,
1565 grad_u_ext[nSpace],
1566 grad_vx_ext[nSpace],grad_vy_ext[nSpace],grad_vz_ext[nSpace],
1567 grad_u_old_ext[nSpace],grad_dissipation_ext_dummy[nSpace],
1568 m_ext=0.0,
1569 dm_ext=0.0,
1570 f_ext[nSpace],
1571 df_ext[nSpace],
1572 dflux_u_u_ext=0.0,
1573 a_ext=0.0,da_ext=0.0,r_ext=0.0,dr_ext=0.0,
1574 bc_u_ext=0.0,
1575 //bc_grad_u_ext[nSpace],
1576 bc_m_ext=0.0,
1577 bc_dm_ext=0.0,
1578 bc_f_ext[nSpace],
1579 bc_df_ext[nSpace],
1580 fluxJacobian_u_u[nDOF_trial_element],
1581 diffusiveFluxJacobian_u_u[nDOF_trial_element],
1582 jac_ext[nSpace*nSpace],
1583 jacDet_ext,
1584 jacInv_ext[nSpace*nSpace],
1585 boundaryJac[nSpace*(nSpace-1)],
1586 metricTensor[(nSpace-1)*(nSpace-1)],
1587 metricTensorDetSqrt,
1588 dS,
1589 u_test_dS[nDOF_test_element],
1590 u_grad_trial_trace[nDOF_trial_element*nSpace],
1591 normal[3],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
1592 G[nSpace*nSpace],G_dd_G,tr_G;
1593 //
1594 //calculate the solution and gradients at quadrature points
1595 //
1596 // u_ext=0.0;
1597 // for (int I=0;I<nSpace;I++)
1598 // {
1599 // grad_u_ext[I] = 0.0;
1600 // bc_grad_u_ext[I] = 0.0;
1601 // }
1602 // for (int j=0;j<nDOF_trial_element;j++)
1603 // {
1604 // int eN_j = eN*nDOF_trial_element+j,
1605 // ebNE_kb_j = ebNE_kb*nDOF_trial_element+j,
1606 // ebNE_kb_j_nSpace= ebNE_kb_j*nSpace;
1607 // u_ext += valFromDOF_c(u_dof[u_l2g[eN_j]],u_trial_ext[ebNE_kb_j]);
1608
1609 // for (int I=0;I<nSpace;I++)
1610 // {
1611 // grad_u_ext[I] += gradFromDOF_c(u_dof[u_l2g[eN_j]],u_grad_trial_ext[ebNE_kb_j_nSpace+I]);
1612 // }
1613 // }
1614 ck.calculateMapping_elementBoundary(eN,
1615 ebN_local,
1616 kb,
1617 ebN_local_kb,
1618 mesh_dof.data(),
1619 mesh_l2g.data(),
1620 mesh_trial_trace_ref.data(),
1621 mesh_grad_trial_trace_ref.data(),
1622 boundaryJac_ref.data(),
1623 jac_ext,
1624 jacDet_ext,
1625 jacInv_ext,
1626 boundaryJac,
1627 metricTensor,
1628 metricTensorDetSqrt,
1629 normal_ref.data(),
1630 normal,
1631 x_ext,y_ext,z_ext);
1632 ck.calculateMappingVelocity_elementBoundary(eN,
1633 ebN_local,
1634 kb,
1635 ebN_local_kb,
1636 mesh_velocity_dof.data(),
1637 mesh_l2g.data(),
1638 mesh_trial_trace_ref.data(),
1639 xt_ext,yt_ext,zt_ext,
1640 normal,
1641 boundaryJac,
1642 metricTensor,
1643 integralScaling);
1644 dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref[kb];
1645 //dS = metricTensorDetSqrt*dS_ref[kb];
1646 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
1647 //compute shape and solution information
1648 //shape
1649 ck.gradTrialFromRef(&u_grad_trial_trace_ref[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,u_grad_trial_trace);
1650 //solution and gradients
1651 ck.valFromDOF(u_dof.data(),&u_l2g[eN_nDOF_trial_element],&u_trial_trace_ref[ebN_local_kb*nDOF_test_element],u_ext);
1652 ck.valFromDOF(u_dof_old.data(),&u_l2g[eN_nDOF_trial_element],&u_trial_trace_ref[ebN_local_kb*nDOF_test_element],u_old_ext);
1653 ck.gradFromDOF(u_dof.data(),&u_l2g[eN_nDOF_trial_element],u_grad_trial_trace,grad_u_ext);
1654 ck.gradFromDOF(u_dof_old.data(),&u_l2g[eN_nDOF_trial_element],u_grad_trial_trace,grad_u_old_ext);
1655
1656 //mwf hack, skip on boundary for now
1657 grad_dissipation_ext_dummy[0] = 0.0; grad_dissipation_ext_dummy[1] = 0.0; grad_dissipation_ext_dummy[2] = 0.0;
1658 //
1659 //compute velocity production terms, ***assumes same spaces for velocity dofs and kappa!***
1660 ck.gradFromDOF(velocity_dof_u.data(),&u_l2g[eN_nDOF_trial_element],u_grad_trial_trace,grad_vx_ext);
1661 ck.gradFromDOF(velocity_dof_v.data(),&u_l2g[eN_nDOF_trial_element],u_grad_trial_trace,grad_vy_ext);
1662 ck.gradFromDOF(velocity_dof_w.data(),&u_l2g[eN_nDOF_trial_element],u_grad_trial_trace,grad_vz_ext);
1663 //
1664
1665 //precalculate test function products with integration weights
1666 for (int j=0;j<nDOF_trial_element;j++)
1667 {
1668 u_test_dS[j] = u_test_trace_ref[ebN_local_kb*nDOF_test_element+j]*dS;
1669 }
1670 //
1671 //load the boundary values
1672 //
1673 bc_u_ext = isDOFBoundary_u[ebNE_kb]*ebqe_bc_u_ext[ebNE_kb]+(1-isDOFBoundary_u[ebNE_kb])*u_ext;
1674 //
1675 //calculate the internal and external trace of the pde coefficients
1676 //
1677 evaluateCoefficients(&ebqe_velocity_ext[ebNE_kb_nSpace],
1678 epsFact,
1679 ebqe_phi[ebNE_kb],
1680 nu_0,
1681 nu_1,
1682 sigma_k,
1683 c_mu,
1684 grad_vx_ext,
1685 grad_vy_ext,
1686 grad_vz_ext,
1687 u_ext,
1688 u_old_ext,
1689 ebqe_dissipation[ebNE_kb],
1690 ebqe_porosity[ebNE_kb],
1691// Argumentlist for sediment
1692 sedFlag,
1693 ebqe_q_vos[ebNE_kb],
1694 &ebqe_q_vos_gradc[ebNE_kb_nSpace],
1695 rho_f,
1696 rho_s,
1697 &ebqe_vs[ebNE_kb_nSpace],
1698 &g[0],
1699 //end sediment
1700 dissipation_model_flag,
1701 grad_u_old_ext,
1702 grad_dissipation_ext_dummy,
1703 m_ext,
1704 dm_ext,
1705 f_ext,
1706 df_ext,
1707 a_ext,
1708 da_ext,
1709 r_ext,
1710 dr_ext);
1711 evaluateCoefficients(&ebqe_velocity_ext[ebNE_kb_nSpace],
1712 epsFact,
1713 ebqe_phi[ebNE_kb],
1714 nu_0,
1715 nu_1,
1716 sigma_k,
1717 c_mu,
1718 grad_vx_ext,
1719 grad_vy_ext,
1720 grad_vz_ext,
1721 bc_u_ext,
1722 bc_u_ext,
1723 ebqe_dissipation[ebNE_kb],
1724 ebqe_porosity[ebNE_kb],
1725// Argumentlist for sediment
1726 sedFlag,
1727 ebqe_q_vos[ebNE_kb],
1728 &ebqe_q_vos_gradc[ebNE_kb_nSpace],
1729 rho_f,
1730 rho_s,
1731 &ebqe_vs[ebNE_kb_nSpace],
1732 &g[0],
1733 //end sediment
1734 dissipation_model_flag,
1735 grad_u_old_ext,
1736 grad_dissipation_ext_dummy,
1737 bc_m_ext,
1738 bc_dm_ext,
1739 bc_f_ext,
1740 bc_df_ext,
1741 a_ext,
1742 da_ext,
1743 r_ext,
1744 dr_ext);
1745 //
1746 //moving domain
1747 //
1748 double velocity_ext[nSpace];
1749 velocity_ext[0] = ebqe_velocity_ext[ebNE_kb_nSpace+0] - MOVING_DOMAIN*xt_ext;
1750 velocity_ext[1] = ebqe_velocity_ext[ebNE_kb_nSpace+1] - MOVING_DOMAIN*yt_ext;
1751 velocity_ext[2] = ebqe_velocity_ext[ebNE_kb_nSpace+2] - MOVING_DOMAIN*zt_ext;
1752 //
1753 //calculate the numerical fluxes
1754 //
1755 exteriorNumericalAdvectiveFluxDerivative(isDOFBoundary_u[ebNE_kb],
1756 isAdvectiveFluxBoundary_u[ebNE_kb],
1757 normal,
1758 velocity_ext,//ebqe_velocity_ext[ebNE_kb_nSpace],
1759 dflux_u_u_ext);
1760 //
1761 //calculate the flux jacobian
1762 //
1763 for (int j=0;j<nDOF_trial_element;j++)
1764 {
1765 //int ebNE_kb_j = ebNE_kb*nDOF_trial_element+j;
1766 int ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j;
1767 //diffusive flux
1768 exteriorNumericalDiffusiveFluxDerivative(isDOFBoundary_u[ebNE_kb],
1769 isDiffusiveFluxBoundary_u[ebNE_kb],
1770 normal,
1771 a_ext,
1772 da_ext,
1773 grad_u_ext,
1774 &u_grad_trial_trace[j*nSpace],
1775 u_trial_trace_ref[ebN_local_kb_j],
1776 ebqe_penalty_ext[ebNE_kb],//penalty,
1777 diffusiveFluxJacobian_u_u[j]);
1778 //mwf debug
1779 //std::cout<<"Jacobian ebNE= "<<ebNE<<" kb= "<<kb <<" penalty= "<<ebqe_penalty_ext[ebNE_kb] <<std::endl;
1780 fluxJacobian_u_u[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_u_u_ext,u_trial_trace_ref[ebN_local_kb_j]);
1781 }//j
1782 //
1783 //update the global Jacobian from the flux Jacobian
1784 //
1785 for (int i=0;i<nDOF_test_element;i++)
1786 {
1787 int eN_i = eN*nDOF_test_element+i;
1788 //int ebNE_kb_i = ebNE_kb*nDOF_test_element+i;
1789 for (int j=0;j<nDOF_trial_element;j++)
1790 {
1791 int ebN_i_j = ebN*4*nDOF_test_X_trial_element + i*nDOF_trial_element + j;
1792
1793 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] += fluxJacobian_u_u[j]*u_test_dS[i]
1794 + diffusiveFluxJacobian_u_u[j]*u_test_dS[i];
1795 }//j
1796 }//i
1797 }//kb
1798 }//ebNE
1799 }//computeJacobian
1800 };//Kappa
1801
1802 inline Kappa_base* newKappa(int nSpaceIn,
1803 int nQuadraturePoints_elementIn,
1804 int nDOF_mesh_trial_elementIn,
1805 int nDOF_trial_elementIn,
1806 int nDOF_test_elementIn,
1807 int nQuadraturePoints_elementBoundaryIn,
1808 int CompKernelFlag,
1809 double aDarcy,
1810 double betaForch,
1811 double grain,
1812 double packFraction,
1813 double packMargin,
1814 double maxFraction,
1815 double frFraction,
1816 double sigmaC,
1817 double C3e,
1818 double C4e,
1819 double eR,
1820 double fContact,
1821 double mContact,
1822 double nContact,
1823 double angFriction,
1824 double vos_limiter,
1825 double mu_fr_limiter)
1826
1827 {
1828 Kappa_base* rvalue =
1829 proteus::chooseAndAllocateDiscretization<Kappa_base,Kappa,CompKernel>
1830 (nSpaceIn,
1831 nQuadraturePoints_elementIn,
1832 nDOF_mesh_trial_elementIn,
1833 nDOF_trial_elementIn,
1834 nDOF_test_elementIn,
1835 nQuadraturePoints_elementBoundaryIn,
1836 CompKernelFlag);
1837
1838 rvalue->setSedClosure(aDarcy,
1839 betaForch,
1840 grain,
1841 packFraction,
1842 packMargin,
1843 maxFraction,
1844 frFraction,
1845 sigmaC,
1846 C3e,
1847 C4e,
1848 eR,
1849 fContact,
1850 mContact,
1851 nContact,
1852 angFriction,
1853 vos_limiter,
1854 mu_fr_limiter);
1855 return rvalue;
1856 }
1857}//proteus
1858#endif
nu_0
double nu_0
Definition ErrorResidualMethod.cpp:22
rho_1
double rho_1
Definition ErrorResidualMethod.cpp:22
nu_1
double nu_1
Definition ErrorResidualMethod.cpp:22
rho_0
double rho_0
Definition ErrorResidualMethod.cpp:22
n
Int n
Definition Headers.h:28
r
Double r
Definition Headers.h:83
H
Double H
Definition Headers.h:65
u
Double u
Definition Headers.h:89
num
Int num
Definition Headers.h:32
z
Double * z
Definition Headers.h:49
v
Double v
Definition Headers.h:95
proteus::Kappa_base
Definition Kappa.h:16
proteus::Kappa_base::setSedClosure
virtual void setSedClosure(double aDarcy, double betaForch, double grain, double packFraction, double packMargin, double maxFraction, double frFraction, double sigmaC, double C3e, double C4e, double eR, double fContact, double mContact, double nContact, double angFriction, double vos_limiter, double mu_fr_limiter)
Definition Kappa.h:20
proteus::Kappa_base::calculateResidual
virtual void calculateResidual(arguments_dict &args)=0
proteus::Kappa_base::~Kappa_base
virtual ~Kappa_base()
Definition Kappa.h:19
proteus::Kappa_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)=0
proteus::Kappa::calculateSubgridError_tau
void calculateSubgridError_tau(const double &elementDiameter, const double &dmt, const double dH[nSpace], double &cfl, double &tau)
Definition Kappa.h:438
proteus::Kappa::closure
cppHsuSedStress< 3 > closure
Definition Kappa.h:51
proteus::Kappa::calculateSubgridError_tau
void calculateSubgridError_tau(const double &Ct_sge, const double G[nSpace *nSpace], const double &A0, const double Ai[nSpace], double &tau_v, double &cfl)
Definition Kappa.h:457
proteus::Kappa::computeK_OmegaCoefficients
void computeK_OmegaCoefficients(const double &div_eps, const double &k, const double &omega, const double grad_k[nSpace], const double grad_omega[nSpace], const double grad_vx[nSpace], const double grad_vy[nSpace], double &inverse_sigma_k, double &inverse_sigma_omega, double &beta_star, double &beta, double &gamma)
Definition Kappa.h:115
proteus::Kappa::ck
CompKernelType ck
Definition Kappa.h:53
proteus::Kappa::calculateNumericalDiffusion
void calculateNumericalDiffusion(const double &shockCapturingDiffusion, const double &elementDiameter, const double &strong_residual, const double grad_u[nSpace], double &numDiff)
Definition Kappa.h:475
proteus::Kappa::smoothedHeaviside
double smoothedHeaviside(double eps, double phi)
Definition Kappa.h:546
proteus::Kappa::computeK_OmegaCoefficients
void computeK_OmegaCoefficients(const double &div_eps, const double &k, const double &omega, const double grad_k[nSpace], const double grad_omega[nSpace], const double grad_vx[nSpace], const double grad_vy[nSpace], const double grad_vz[nSpace], double &inverse_sigma_k, double &inverse_sigma_omega, double &beta_star, double &beta, double &gamma)
Definition Kappa.h:355
proteus::Kappa::exteriorNumericalDiffusiveFlux
void exteriorNumericalDiffusiveFlux(const double &bc_flux, const int &isDOFBoundary, const int &isDiffusiveFluxBoundary, double n[nSpace], double bc_u, double a, double grad_psi[nSpace], double u, double penalty, double &flux)
Definition Kappa.h:596
proteus::Kappa::evaluateCoefficients
void evaluateCoefficients(double v[nSpace], const double eps_mu, const double phi, const double nu_0, const double nu_1, const double sigma_k, const double c_mu, const double grad_vx[nSpace], const double grad_vy[nSpace], const double grad_vz[nSpace], const double &k, double &k_old, double &dissipation, const double &porosity, int sedFlag, double q_vos, double q_vos_gradc[nSpace], double rho_f, double rho_s, double vs[nSpace], double g[nSpace], int dissipation_model_flag, const double grad_k_old[nSpace], const double grad_dissipation[nSpace], double &m, double &dm, double f[nSpace], double df[nSpace], double &a, double &da_dk, double &r, double &dr_dk)
Definition Kappa.h:194
proteus::Kappa::exteriorNumericalDiffusiveFluxDerivative
void exteriorNumericalDiffusiveFluxDerivative(const int &isDOFBoundary, const int &isDiffusiveFluxBoundary, double n[nSpace], double a, double da, double grad_psi[nSpace], const double grad_v[nSpace], double v, double penalty, double &fluxJacobian)
Definition Kappa.h:630
proteus::Kappa::setSedClosure
void setSedClosure(double aDarcy, double betaForch, double grain, double packFraction, double packMargin, double maxFraction, double frFraction, double sigmaC, double C3e, double C4e, double eR, double fContact, double mContact, double nContact, double angFriction, double vos_limiter, double mu_fr_limiter)
Definition Kappa.h:77
proteus::Kappa::exteriorNumericalAdvectiveFlux
void exteriorNumericalAdvectiveFlux(const int &isDOFBoundary_u, const int &isAdvectiveFluxBoundary_u, const double n[nSpace], const double &bc_u, const double &bc_flux_u, const double &u, const double velocity[nSpace], double &flux)
Definition Kappa.h:495
proteus::Kappa::exteriorNumericalAdvectiveFluxDerivative
void exteriorNumericalAdvectiveFluxDerivative(const int &isDOFBoundary_u, const int &isAdvectiveFluxBoundary, const double n[nSpace], const double velocity[nSpace], double &dflux)
Definition Kappa.h:561
proteus::Kappa::calculateResidual
void calculateResidual(arguments_dict &args)
Definition Kappa.h:654
proteus::Kappa::calculateJacobian
void calculateJacobian(arguments_dict &args)
Definition Kappa.h:1234
proteus::Kappa::nDOF_test_X_trial_element
const int nDOF_test_X_trial_element
Definition Kappa.h:52
proteus::cppHsuSedStress::kappa_sed1
double kappa_sed1(double sedF, double rhoFluid, double rhoSolid, const double uFluid[nSpace], const double uSolid[nSpace], const double gradC[nSpace], double nu, double theta_n, double kappa_n, double epsilon_n, double nuT_n, const double g[nSpace])
Definition SedClosure.h:231
proteus::cppHsuSedStress::dkappa_sed1_dk
double dkappa_sed1_dk(double sedF, double rhoFluid, double rhoSolid, const double uFluid[nSpace], const double uSolid[nSpace], const double gradC[nSpace], double nu, double theta_n, double kappa_n, double epsilon_n, double nuT_n)
Definition SedClosure.h:296
df
double df(double C, double b, double a, int q, int r)
Definition analyticalSolutions.c:2209
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::newKappa
Kappa_base * newKappa(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag, double aDarcy, double betaForch, double grain, double packFraction, double packMargin, double maxFraction, double frFraction, double sigmaC, double C3e, double C4e, double eR, double fContact, double mContact, double nContact, double angFriction, double vos_limiter, double mu_fr_limiter)
Definition Kappa.h:1802
proteus::chooseAndAllocateDiscretization
Model_Base * chooseAndAllocateDiscretization(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:48
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