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