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