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