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