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