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