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