proteus  1.8.1
C/C++/Fortran libraries
NCLS3P.h
Go to the documentation of this file.
1 #ifndef NCLS3P_H
2 #define NCLS3P_H
3 #include <cmath>
4 #include <iostream>
5 #include "CompKernel.h"
6 #include "ModelFactory.h"
7 #include "ArgumentsDict.h"
8 #include "xtensor-python/pyarray.hpp"
9 
10 #define cE 0.1
11 #define cMax 0.1
12 
13 namespace proteus
14 {
15  inline double ENTROPY(const double& u){
16  return u;
17  }
18  inline double DENTROPY(const double& u){
19  return 1.0;
20  }
21  inline double Sign(const double& z){
22  return (z >= 0.0 ? 1.0 : -1.0);
23  }
24 }
25 
26 namespace proteus
27 {
28  class cppNCLS3P_base
29  {
30  //The base class defining the interface
31  public:
32  virtual ~cppNCLS3P_base(){}
33  virtual void calculateResidual(arguments_dict& args)=0;
34  virtual void calculateJacobian(arguments_dict& args)=0;
35  virtual void calculateWaterline(arguments_dict& args)=0;
36  };
37 
38  template<class CompKernelType,
39  int nSpace,
40  int nQuadraturePoints_element,
41  int nDOF_mesh_trial_element,
42  int nDOF_trial_element,
43  int nDOF_test_element,
44  int nQuadraturePoints_elementBoundary>
45  class cppNCLS3P : public cppNCLS3P_base
46  {
47  public:
48  const int nDOF_test_X_trial_element;
49  CompKernelType ck;
50  cppNCLS3P():
51  nDOF_test_X_trial_element(nDOF_test_element*nDOF_trial_element),
52  ck()
53  {}
54 
55  inline
56  void calculateCFL(const double& elementDiameter,
57  const double df[nSpace],
58  double& cfl)
59  {
60  double h,nrm_v;
61  h = elementDiameter;
62  nrm_v=0.0;
63  for(int I=0;I<nSpace;I++)
64  nrm_v+=df[I]*df[I];
65  nrm_v = sqrt(nrm_v);
66  cfl = nrm_v/h;
67  }
68 
69  inline void evaluateCoefficients(const double v[nSpace],
70  const double& u,
71  const double grad_u[nSpace],
72  double& m,
73  double& dm,
74  double& H,
75  double dH[nSpace])
76  {
77  m = u;
78  dm=1.0;
79  H = 0.0;
80  for (int I=0; I < nSpace; I++)
81  {
82  H += v[I]*grad_u[I];
83  dH[I] = v[I];
84  }
85  }
86 
87  inline
88  void calculateSubgridError_tau(const double& elementDiameter,
89  const double& dmt,
90  const double dH[nSpace],
91  double& cfl,
92  double& tau)
93  {
94  double h,nrm_v,oneByAbsdt;
95  h = elementDiameter;
96  nrm_v=0.0;
97  for(int I=0;I<nSpace;I++)
98  nrm_v+=dH[I]*dH[I];
99  nrm_v = sqrt(nrm_v);
100  cfl = nrm_v/h;
101  oneByAbsdt = fabs(dmt);
102  tau = 1.0/(2.0*nrm_v/h + oneByAbsdt + 1.0e-8);
103  }
104 
105 
106  inline
107  void calculateSubgridError_tau( const double& Ct_sge,
108  const double G[nSpace*nSpace],
109  const double& A0,
110  const double Ai[nSpace],
111  double& tau_v,
112  double& cfl)
113  {
114  double v_d_Gv=0.0;
115  for(int I=0;I<nSpace;I++)
116  for (int J=0;J<nSpace;J++)
117  v_d_Gv += Ai[I]*G[I*nSpace+J]*Ai[J];
118 
119  tau_v = 1.0/sqrt(Ct_sge*A0*A0 + v_d_Gv + 1.0e-8);
120  }
121 
122  void exteriorNumericalFlux(const double n[nSpace],
123  const double& bc_u,
124  const double& u,
125  const double velocity[nSpace],
126  const double velocity_movingDomain[nSpace],
127  double& flux)
128  {
129  double flow_total=0.0,flow_fluid=0.0,flow_movingDomain=0.0;
130  for (int I=0; I < nSpace; I++)
131  {
132  flow_fluid += n[I]*velocity[I];
133  //flow_movingDomain -= n[I]*velocity_movingDomain[I];
134  }
135  flow_total = flow_fluid+flow_movingDomain;
136  if (flow_total > 0.0)
137  {
138  flux = u*flow_movingDomain;
139  }
140  else
141  {
142 
143  flux = bc_u*flow_movingDomain - flow_fluid*(u-bc_u);
144  //std::cout<<"bc_u "<<bc_u<<" flow_fluid "<<flow_fluid<<" u "<<u<<std::endl;
145  }
146  }
147 
148  inline
149  void exteriorNumericalFluxDerivative(const double n[nSpace],
150  const double velocity[nSpace],
151  const double velocity_movingDomain[nSpace],
152  double& dflux)
153  {
154  double flow_total=0.0,flow_fluid=0.0,flow_movingDomain=0.0;
155  for (int I=0; I < nSpace; I++)
156  {
157  flow_fluid += n[I]*velocity[I];
158  //flow_movingDomain -= n[I]*velocity_movingDomain[I];
159  }
160  flow_total=flow_fluid+flow_movingDomain;
161  if (flow_total > 0.0)
162  {
163  dflux = flow_movingDomain;
164  }
165  else
166  {
167  dflux = -flow_fluid;
168  }
169  }
170 
171  void calculateResidual(arguments_dict& args)
172  {
173  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
174  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
175  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
176  xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
177  double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
178  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
179  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
180  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
181  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
182  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
183  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
184  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
185  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
186  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
187  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
188  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
189  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
190  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
191  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
192  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
193  int nElements_global = args.scalar<int>("nElements_global");
194  double useMetrics = args.scalar<double>("useMetrics");
195  double alphaBDF = args.scalar<double>("alphaBDF");
196  int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
197  double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
198  double sc_uref = args.scalar<double>("sc_uref");
199  double sc_alpha = args.scalar<double>("sc_alpha");
200  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
201  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
202  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
203  xt::pyarray<double>& u_dof_old = args.array<double>("u_dof_old");
204  xt::pyarray<double>& velocity = args.array<double>("velocity");
205  xt::pyarray<double>& q_m = args.array<double>("q_m");
206  xt::pyarray<double>& q_u = args.array<double>("q_u");
207  xt::pyarray<double>& q_n = args.array<double>("q_n");
208  xt::pyarray<double>& q_dH = args.array<double>("q_dH");
209  xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
210  xt::pyarray<double>& q_dV = args.array<double>("q_dV");
211  xt::pyarray<double>& q_dV_last = args.array<double>("q_dV_last");
212  xt::pyarray<double>& cfl = args.array<double>("cfl");
213  xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
214  xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
215  int offset_u = args.scalar<int>("offset_u");
216  int stride_u = args.scalar<int>("stride_u");
217  xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
218  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
219  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
220  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
221  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
222  xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
223  xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
224  xt::pyarray<double>& ebqe_rd_u_ext = args.array<double>("ebqe_rd_u_ext");
225  xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
226  xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
227  xt::pyarray<double>& cell_interface_locator = args.array<double>("cell_interface_locator");
228  xt::pyarray<double>& interface_locator = args.array<double>("interface_locator");
229  int EXPLICIT_METHOD = args.scalar<int>("EXPLICIT_METHOD");
230  double degree_polynomial = args.scalar<double>("degree_polynomial");
231  int stage = args.scalar<int>("stage");
232  xt::pyarray<double>& uTilde_dof = args.array<double>("uTilde_dof");
233  double dt = args.scalar<double>("dt");
234  int PURE_BDF = args.scalar<int>("PURE_BDF");
235  double meanEntropy = 0., meanOmega = 0., maxEntropy = -1E10, minEntropy = 1E10;
236  register double maxVel[nElements_global], maxEntRes[nElements_global];
237  //cek should this be read in?
238  double Ct_sge = 4.0;
239 
240  //loop over elements to compute volume integrals and load them into element and global residual
241  //
242  //eN is the element index
243  //eN_k is the quadrature point index for a scalar
244  //eN_k_nSpace is the quadrature point index for a vector
245  //eN_i is the element test function index
246  //eN_j is the element trial function index
247  //eN_k_j is the quadrature point index for a trial function
248  //eN_k_i is the quadrature point index for a trial function
249  for(int eN=0;eN<nElements_global;eN++)
250  {
251  // init maxVel and maxEntRes
252  maxVel[eN] = 0.;
253  maxEntRes[eN] = 0.;
254  //declare local storage for element residual and initialize
255  register double elementResidual_u[nDOF_test_element];
256  for (int i=0;i<nDOF_test_element;i++)
257  {
258  elementResidual_u[i]=0.0;
259  }//i
260  //loop over quadrature points and compute integrands
261  for (int k=0;k<nQuadraturePoints_element;k++)
262  {
263  //compute indeces and declare local storage
264  register int eN_k = eN*nQuadraturePoints_element+k,
265  eN_k_nSpace = eN_k*nSpace,
266  eN_nDOF_trial_element = eN*nDOF_trial_element;
267  register double u=0.0,grad_u[nSpace],grad_u_old[nSpace],
268  un=0.0, Hn=0.0, HTilde=0.0, grad_uTilde[nSpace],
269  m=0.0,dm=0.0,
270  H=0.0,dH[nSpace],
271  f[nSpace],df[nSpace],//for MOVING_DOMAIN
272  m_t=0.0,dm_t=0.0,
273  pdeResidual_u=0.0,
274  Lstar_u[nDOF_test_element],
275  subgridError_u=0.0,
276  tau=0.0,tau0=0.0,tau1=0.0,
277  numDiff0=0.0,numDiff1=0.0,
278  jac[nSpace*nSpace],
279  jacDet,
280  jacInv[nSpace*nSpace],
281  u_grad_trial[nDOF_trial_element*nSpace],
282  u_test_dV[nDOF_trial_element],
283  u_grad_test_dV[nDOF_test_element*nSpace],
284  dV,x,y,z,xt,yt,zt,
285  G[nSpace*nSpace],G_dd_G,tr_G;//,norm_Rv;
286  //
287  //compute solution and gradients at quadrature points
288  //
289  ck.calculateMapping_element(eN,
290  k,
291  mesh_dof.data(),
292  mesh_l2g.data(),
293  mesh_trial_ref.data(),
294  mesh_grad_trial_ref.data(),
295  jac,
296  jacDet,
297  jacInv,
298  x,y,z);
299  ck.calculateMappingVelocity_element(eN,
300  k,
301  mesh_velocity_dof.data(),
302  mesh_l2g.data(),
303  mesh_trial_ref.data(),
304  xt,yt,zt);
305  //get the physical integration weight
306  dV = fabs(jacDet)*dV_ref[k];
307  ck.calculateG(jacInv,G,G_dd_G,tr_G);
308  //get the trial function gradients
309  ck.gradTrialFromRef(&u_grad_trial_ref[k*nDOF_trial_element*nSpace],
310  jacInv,
311  u_grad_trial);
312  //get the solution
313  ck.valFromDOF(u_dof.data(),
314  &u_l2g[eN_nDOF_trial_element],
315  &u_trial_ref[k*nDOF_trial_element],
316  u);
317  ck.valFromDOF(u_dof_old.data(),
318  &u_l2g[eN_nDOF_trial_element],
319  &u_trial_ref[k*nDOF_trial_element],
320  un);
321  //get the solution gradients
322  ck.gradFromDOF(u_dof.data(),
323  &u_l2g[eN_nDOF_trial_element],
324  u_grad_trial,
325  grad_u);
326  ck.gradFromDOF(u_dof_old.data(),
327  &u_l2g[eN_nDOF_trial_element],
328  u_grad_trial,
329  grad_u_old);
330  ck.gradFromDOF(uTilde_dof.data(),
331  &u_l2g[eN_nDOF_trial_element],
332  u_grad_trial,
333  grad_uTilde);
334  //precalculate test function products with integration weights
335  for (int j=0;j<nDOF_trial_element;j++)
336  {
337  u_test_dV[j] = u_test_ref[k*nDOF_trial_element+j]*dV;
338  for (int I=0;I<nSpace;I++)
339  {
340  u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
341  }
342  }
343  //save solution at quadrature points for other models to use
344  q_u[eN_k]=u;
345  for (int I=0;I<nSpace;I++)
346  q_n[eN_k_nSpace+I] = grad_u[I];
347  //
348  //calculate pde coefficients at quadrature points
349  //
350  evaluateCoefficients(&velocity[eN_k_nSpace],
351  u,
352  grad_u,
353  m,
354  dm,
355  H,
356  dH);
357  //
358  //moving mesh
359  //
360  double mesh_velocity[3];
361  mesh_velocity[0] = xt;
362  mesh_velocity[1] = yt;
363  mesh_velocity[2] = zt;
364  for (int I=0;I<nSpace;I++)
365  {
366  f[I] = -MOVING_DOMAIN*m*mesh_velocity[I];
367  df[I] = -MOVING_DOMAIN*dm*mesh_velocity[I];
368  }
369  //
370  //calculate time derivative at quadrature points
371  //
372  if (q_dV_last[eN_k] <= -100)
373  q_dV_last[eN_k] = dV;
374  q_dV[eN_k] = dV;
375  ck.bdf(alphaBDF,
376  q_m_betaBDF[eN_k],
377  m,
378  dm,
379  m_t,
380  dm_t);
381  if (EXPLICIT_METHOD==1)
382  {
383  double normVel=0.;
384  double relVelocity[nSpace];
385  for (int I=0;I<nSpace;I++)
386  {
387  Hn += velocity[eN_k_nSpace+I]*grad_u_old[I];
388  HTilde += velocity[eN_k_nSpace+I]*grad_uTilde[I];
389  H += velocity[eN_k_nSpace+I]*grad_u[I];
390  relVelocity[I] = dH[I] - MOVING_DOMAIN*df[I];
391  normVel += relVelocity[I]*relVelocity[I];
392  }
393  normVel = std::sqrt(normVel);
394 
395  // calculate CFL
396  calculateCFL(elementDiameter[eN]/degree_polynomial,relVelocity,cfl[eN_k]);
397 
398  // compute max velocity at cell
399  maxVel[eN] = fmax(normVel,maxVel[eN]);
400 
401  // entropy residual
402  double entRes = (ENTROPY(u)-ENTROPY(un))/dt + 0.5*(DENTROPY(u)*H +
403  DENTROPY(un)*Hn);
404  maxEntRes[eN] = fmax(maxEntRes[eN],fabs(entRes));
405 
406  // Quantities for normalization factor //
407  meanEntropy += ENTROPY(u)*dV;
408  meanOmega += dV;
409  maxEntropy = fmax(maxEntropy,ENTROPY(u));
410  minEntropy = fmin(minEntropy,ENTROPY(u));
411  }
412  else
413  {
414  //
415  //calculate subgrid error (strong residual and adjoint)
416  //
417  //calculate strong residual
418  pdeResidual_u = ck.Mass_strong(m_t) +
419  ck.Hamiltonian_strong(dH,grad_u)+
420  MOVING_DOMAIN*ck.Advection_strong(df,grad_u);//cek I don't think all mesh motion will be divergence free so we may need to go add the divergence
421 
422  //calculate adjoint
423  for (int i=0;i<nDOF_test_element;i++)
424  {
425  //register int eN_k_i_nSpace = (eN_k*nDOF_trial_element+i)*nSpace;
426  //Lstar_u[i] = ck.Hamiltonian_adjoint(dH,&u_grad_test_dV[eN_k_i_nSpace]);
427  register int i_nSpace = i*nSpace;
428  Lstar_u[i] = ck.Hamiltonian_adjoint(dH,&u_grad_test_dV[i_nSpace])
429  + MOVING_DOMAIN*ck.Advection_adjoint(df,&u_grad_test_dV[i_nSpace]);
430  }
431  //calculate tau and tau*Res
432  double subgridErrorVelocity[nSpace];
433  for (int I=0;I<nSpace;I++)
434  subgridErrorVelocity[I] = dH[I] - MOVING_DOMAIN*df[I];
435 
436  calculateSubgridError_tau(elementDiameter[eN],
437  dm_t,
438  subgridErrorVelocity,//dH,
439  cfl[eN_k],
440  tau0);
441 
442  calculateSubgridError_tau(Ct_sge,
443  G,
444  dm_t,
445  subgridErrorVelocity,//dH,
446  tau1,
447  cfl[eN_k]);
448 
449  tau = useMetrics*tau1+(1.0-useMetrics)*tau0;
450 
451  subgridError_u = -tau*pdeResidual_u;
452  //
453  //calculate shock capturing diffusion
454  //
455  ck.calculateNumericalDiffusion(shockCapturingDiffusion,
456  elementDiameter[eN],
457  pdeResidual_u,
458  grad_u,
459  numDiff0);
460  //ck.calculateNumericalDiffusion(shockCapturingDiffusion,G,pdeResidual_u,grad_u_old,numDiff1);
461  ck.calculateNumericalDiffusion(shockCapturingDiffusion,
462  sc_uref,
463  sc_alpha,
464  G,
465  G_dd_G,
466  pdeResidual_u,
467  grad_u,
468  numDiff1);
469  q_numDiff_u[eN_k] = useMetrics*numDiff1+(1.0-useMetrics)*numDiff0;
470  //std::cout<<tau<<" "<<q_numDiff_u[eN_k]<<std::endl;
471  }
472  //
473  //update element residual
474  //
475  double sign;
476  int same_sign=1;
477  for(int i=0;i<nDOF_test_element;i++)
478  {
479  int gi = offset_u+stride_u*u_l2g[eN*nDOF_test_element+i];
480  if (i==0)
481  sign = Sign(u_dof[gi]);
482  else if (same_sign==1)
483  {
484  same_sign = sign == Sign(u_dof[gi]) ? 1 : 0;
485  sign = Sign(u_dof[gi]);
486  }
487  //register int eN_k_i=eN_k*nDOF_test_element+i,
488  // eN_k_i_nSpace = eN_k_i*nSpace;
489  register int i_nSpace=i*nSpace;
490  if (EXPLICIT_METHOD==1)
491  {
492  if (stage == 1)
493  elementResidual_u[i] +=
494  //ck.Mass_weak(u-un,u_test_dV[i]) + // time derivative
495  ck.Mass_weak(dt*m_t,u_test_dV[i]) + // time derivative
496  1./3*dt*ck.Hamiltonian_weak(Hn,u_test_dV[i]) + // v*grad(phi)
497  1./9*dt*dt*ck.NumericalDiffusion(Hn,dH,&u_grad_test_dV[i_nSpace]) +
498  1./3*dt*ck.NumericalDiffusion(q_numDiff_u_last[eN_k],
499  grad_u_old,
500  &u_grad_test_dV[i_nSpace]);
501  // TODO: Add part about moving mesh
502  else
503  elementResidual_u[i] +=
504  //ck.Mass_weak(u-un,u_test_dV[i]) + // time derivative
505  ck.Mass_weak(dt*m_t,u_test_dV[i]) + // time derivative
506  dt*ck.Hamiltonian_weak(Hn,u_test_dV[i]) + // v*grad(phi)
507  0.5*dt*dt*ck.NumericalDiffusion(HTilde,dH,&u_grad_test_dV[i_nSpace]) +
508  dt*ck.NumericalDiffusion(q_numDiff_u_last[eN_k],
509  grad_u_old,
510  &u_grad_test_dV[i_nSpace]);
511  //elementResidual_u[i] += // semi-implicit Lax Wendroff
512  // ck.Mass_weak(u-un,u_test_dV[i]) +
513  // dt*ck.Hamiltonian_weak(Hn,u_test_dV[i]) +
514  // 0.5*dt*dt*ck.NumericalDiffusion(H,dH,&u_grad_test_dV[i_nSpace]);
515  }
516  else // Implicit SUPG with SHOCK CAPTURING
517  {
518  elementResidual_u[i] +=
519  ck.Mass_weak(m_t,u_test_dV[i]) +
520  ck.Hamiltonian_weak(H,u_test_dV[i]) +
521  MOVING_DOMAIN*ck.Advection_weak(f,&u_grad_test_dV[i_nSpace])+
522  (PURE_BDF == 1 ? 0. : 1.)*
523  (ck.SubgridError(subgridError_u,Lstar_u[i]) +
524  ck.NumericalDiffusion(q_numDiff_u_last[eN_k],
525  grad_u,
526  &u_grad_test_dV[i_nSpace]));
527  }
528  }//i
529  if (same_sign == 0) // This cell contains the interface
530  {
531  cell_interface_locator[eN] = 1.0;
532  for(int i=0;i<nDOF_test_element;i++)
533  {
534  int gi = offset_u+stride_u*u_l2g[eN*nDOF_test_element+i];
535  interface_locator[gi] = 1.0;
536  }
537  }
538  //
539  //cek/ido todo, get rid of m, since u=m
540  //save momentum for time history and velocity for subgrid error
541  //save solution for other models
542  //
543  q_m[eN_k] = m;
544  q_u[eN_k] = u;
545  for (int I=0;I<nSpace;I++)
546  {
547  int eN_k_I = eN_k*nSpace+I;
548  q_dH[eN_k_I] = dH[I];
549  }
550  }
551  //
552  //load element into global residual and save element residual
553  //
554  for(int i=0;i<nDOF_test_element;i++)
555  {
556  register int eN_i=eN*nDOF_test_element+i;
557  globalResidual[offset_u+stride_u*u_l2g[eN_i]] += elementResidual_u[i];
558  }//i
559  }//elements
560  if (EXPLICIT_METHOD==1)
561  {
562  meanEntropy /= meanOmega;
563  double norm_factor = fmax(fabs(maxEntropy - meanEntropy), fabs(meanEntropy-minEntropy));
564  for(int eN=0;eN<nElements_global;eN++)
565  {
566  double hK=elementDiameter[eN]/degree_polynomial;
567  double linear_viscosity = cMax*hK*maxVel[eN];
568  double entropy_viscosity = cE*hK*hK*maxEntRes[eN]/norm_factor;
569  for (int k=0;k<nQuadraturePoints_element;k++)
570  {
571  register int eN_k = eN*nQuadraturePoints_element+k;
572  q_numDiff_u[eN_k] = fmin(linear_viscosity,entropy_viscosity);
573  }
574  }
575  }
576  //
577  //loop over exterior element boundaries to calculate surface integrals and load into element and global residuals
578  //
579  //ebNE is the Exterior element boundary INdex
580  //ebN is the element boundary INdex
581  //eN is the element index
582  //std::cout <<nExteriorElementBoundaries_global<<std::endl;
583  for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
584  {
585  register int ebN = exteriorElementBoundariesArray[ebNE],
586  eN = elementBoundaryElementsArray[ebN*2+0],
587  ebN_local = elementBoundaryLocalElementBoundariesArray[ebN*2+0],
588  eN_nDOF_trial_element = eN*nDOF_trial_element;
589  register double elementResidual_u[nDOF_test_element];
590  for (int i=0;i<nDOF_test_element;i++)
591  {
592  elementResidual_u[i]=0.0;
593  }
594  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
595  {
596  register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
597  ebNE_kb_nSpace = ebNE_kb*nSpace,
598  ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
599  ebN_local_kb_nSpace = ebN_local_kb*nSpace;
600  register double u_ext=0.0,
601  grad_u_ext[nSpace],
602  m_ext=0.0,
603  dm_ext=0.0,
604  H_ext=0.0,
605  dH_ext[nSpace],
606  //f_ext[nSpace],//MOVING_DOMAIN
607  //df_ext[nSpace],//MOVING_DOMAIN
608  //flux_ext=0.0,
609  bc_u_ext=0.0,
610  bc_grad_u_ext[nSpace],
611  bc_m_ext=0.0,
612  bc_dm_ext=0.0,
613  bc_H_ext=0.0,
614  bc_dH_ext[nSpace],
615  jac_ext[nSpace*nSpace],
616  jacDet_ext,
617  jacInv_ext[nSpace*nSpace],
618  boundaryJac[nSpace*(nSpace-1)],
619  metricTensor[(nSpace-1)*(nSpace-1)],
620  metricTensorDetSqrt,
621  dS,
622  u_test_dS[nDOF_test_element],
623  u_grad_trial_trace[nDOF_trial_element*nSpace],
624  normal[nSpace],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
625  G[nSpace*nSpace],G_dd_G,tr_G,flux_ext;
626  //
627  //calculate the solution and gradients at quadrature points
628  //
629  //compute information about mapping from reference element to physical element
630  ck.calculateMapping_elementBoundary(eN,
631  ebN_local,
632  kb,
633  ebN_local_kb,
634  mesh_dof.data(),
635  mesh_l2g.data(),
636  mesh_trial_trace_ref.data(),
637  mesh_grad_trial_trace_ref.data(),
638  boundaryJac_ref.data(),
639  jac_ext,
640  jacDet_ext,
641  jacInv_ext,
642  boundaryJac,
643  metricTensor,
644  metricTensorDetSqrt,
645  normal_ref.data(),
646  normal,
647  x_ext,y_ext,z_ext);
648  ck.calculateMappingVelocity_elementBoundary(eN,
649  ebN_local,
650  kb,
651  ebN_local_kb,
652  mesh_velocity_dof.data(),
653  mesh_l2g.data(),
654  mesh_trial_trace_ref.data(),
655  xt_ext,yt_ext,zt_ext,
656  normal,
657  boundaryJac,
658  metricTensor,
659  integralScaling);
660  dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref[kb];
661  //get the metric tensor
662  //cek todo use symmetry
663  ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
664  //compute shape and solution information
665  //shape
666  ck.gradTrialFromRef(&u_grad_trial_trace_ref[ebN_local_kb_nSpace*nDOF_trial_element],
667  jacInv_ext,
668  u_grad_trial_trace);
669  //solution and gradients
670  if (EXPLICIT_METHOD==1) // explicit
671  {
672  ck.valFromDOF(u_dof_old.data(),
673  &u_l2g[eN_nDOF_trial_element],
674  &u_trial_trace_ref[ebN_local_kb*nDOF_test_element],
675  u_ext);
676  ck.gradFromDOF(u_dof_old.data(),
677  &u_l2g[eN_nDOF_trial_element],
678  u_grad_trial_trace,
679  grad_u_ext);
680  }
681  else
682  {
683  ck.valFromDOF(u_dof.data(),
684  &u_l2g[eN_nDOF_trial_element],
685  &u_trial_trace_ref[ebN_local_kb*nDOF_test_element],
686  u_ext);
687  ck.gradFromDOF(u_dof.data(),
688  &u_l2g[eN_nDOF_trial_element],
689  u_grad_trial_trace,
690  grad_u_ext);
691  }
692  //precalculate test function products with integration weights
693  for (int j=0;j<nDOF_trial_element;j++)
694  {
695  u_test_dS[j] = u_test_trace_ref[ebN_local_kb*nDOF_test_element+j]*dS;
696  }
697  //
698  //load the boundary values
699  //
700  bc_u_ext = (isDOFBoundary_u[ebNE_kb]*ebqe_bc_u_ext[ebNE_kb]
701  +(1-isDOFBoundary_u[ebNE_kb])*ebqe_rd_u_ext[ebNE_kb]);
702  //
703  //calculate the pde coefficients using the solution and the boundary values for the solution
704  //
705  evaluateCoefficients(&ebqe_velocity_ext[ebNE_kb_nSpace],
706  u_ext,
707  grad_u_ext,
708  m_ext,
709  dm_ext,
710  H_ext,
711  dH_ext);
712  evaluateCoefficients(&ebqe_velocity_ext[ebNE_kb_nSpace],
713  bc_u_ext,
714  bc_grad_u_ext,
715  bc_m_ext,
716  bc_dm_ext,
717  bc_H_ext,
718  bc_dH_ext);
719  //
720  //moving mesh
721  //
722  double velocity_ext[nSpace];
723  double mesh_velocity[3];
724  mesh_velocity[0] = xt_ext;
725  mesh_velocity[1] = yt_ext;
726  mesh_velocity[2] = zt_ext;
727  for (int I=0;I<nSpace;I++)
728  velocity_ext[I] = - MOVING_DOMAIN*mesh_velocity[I];
729  //
730  //calculate the numerical fluxes
731  //
732  exteriorNumericalFlux(normal,
733  bc_u_ext,
734  u_ext,
735  dH_ext,
736  velocity_ext,
737  flux_ext);
738  ebqe_u[ebNE_kb] = u_ext;
739 
740  if (EXPLICIT_METHOD==1)
741  if (stage==1)
742  flux_ext *= 1./3*dt;
743  else
744  flux_ext *= dt;
745  //std::cout<<u_ext<<ebqe_bc_u_ext
746  //
747  //update residuals
748  //
749  for (int i=0;i<nDOF_test_element;i++)
750  {
751  //int ebNE_kb_i = ebNE_kb*nDOF_test_element+i;
752  elementResidual_u[i] += ck.ExteriorElementBoundaryFlux(flux_ext,u_test_dS[i]);
753  }//i
754  }//kb
755  //
756  //update the element and global residual storage
757  //
758  for (int i=0;i<nDOF_test_element;i++)
759  {
760  int eN_i = eN*nDOF_test_element+i;
761  //globalResidual[offset_u+stride_u*u_l2g[eN_i]] += MOVING_DOMAIN*elementResidual_u[i];
762  globalResidual[offset_u+stride_u*u_l2g[eN_i]] += elementResidual_u[i];
763 
764  }//i
765  }//ebNE
766  }
767 
768  void calculateJacobian(arguments_dict& args)
769  {
770  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
771  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
772  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
773  xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
774  double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
775  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
776  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
777  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
778  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
779  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
780  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
781  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
782  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
783  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
784  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
785  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
786  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
787  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
788  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
789  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
790  int nElements_global = args.scalar<int>("nElements_global");
791  double useMetrics = args.scalar<double>("useMetrics");
792  double alphaBDF = args.scalar<double>("alphaBDF");
793  int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
794  double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
795  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
796  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
797  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
798  xt::pyarray<double>& velocity = args.array<double>("velocity");
799  xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
800  xt::pyarray<double>& cfl = args.array<double>("cfl");
801  xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
802  xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
803  xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
804  xt::pyarray<double>& globalJacobian = args.array<double>("globalJacobian");
805  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
806  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
807  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
808  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
809  xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
810  xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
811  xt::pyarray<double>& ebqe_rd_u_ext = args.array<double>("ebqe_rd_u_ext");
812  xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
813  xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
814  int EXPLICIT_METHOD = args.scalar<int>("EXPLICIT_METHOD");
815  int PURE_BDF = args.scalar<int>("PURE_BDF");
816  double Ct_sge = 4.0;
817 
818  //
819  //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
820  //
821  for(int eN=0;eN<nElements_global;eN++)
822  {
823  register double elementJacobian_u_u[nDOF_test_element][nDOF_trial_element];
824  for (int i=0;i<nDOF_test_element;i++)
825  for (int j=0;j<nDOF_trial_element;j++)
826  {
827  elementJacobian_u_u[i][j]=0.0;
828  }
829  for (int k=0;k<nQuadraturePoints_element;k++)
830  {
831  int eN_k = eN*nQuadraturePoints_element+k, //index to a scalar at a quadrature point
832  eN_k_nSpace = eN_k*nSpace,
833  eN_nDOF_trial_element = eN*nDOF_trial_element; //index to a vector at a quadrature point
834 
835  //declare local storage
836  register double u=0.0,
837  grad_u[nSpace],
838  m=0.0,dm=0.0,
839  H=0.0,dH[nSpace],
840  f[nSpace],df[nSpace],//MOVING_MESH
841  m_t=0.0,dm_t=0.0,
842  dpdeResidual_u_u[nDOF_trial_element],
843  Lstar_u[nDOF_test_element],
844  dsubgridError_u_u[nDOF_trial_element],
845  tau=0.0,tau0=0.0,tau1=0.0,
846  jac[nSpace*nSpace],
847  jacDet,
848  jacInv[nSpace*nSpace],
849  u_grad_trial[nDOF_trial_element*nSpace],
850  dV,
851  u_test_dV[nDOF_test_element],
852  u_grad_test_dV[nDOF_test_element*nSpace],
853  x,y,z,xt,yt,zt,
854  G[nSpace*nSpace],G_dd_G,tr_G;
855  //
856  //calculate solution and gradients at quadrature points
857  //
858  //get jacobian, etc for mapping reference element
859  ck.calculateMapping_element(eN,
860  k,
861  mesh_dof.data(),
862  mesh_l2g.data(),
863  mesh_trial_ref.data(),
864  mesh_grad_trial_ref.data(),
865  jac,
866  jacDet,
867  jacInv,
868  x,y,z);
869  ck.calculateMappingVelocity_element(eN,
870  k,
871  mesh_velocity_dof.data(),
872  mesh_l2g.data(),
873  mesh_trial_ref.data(),
874  xt,yt,zt);
875  //get the physical integration weight
876  dV = fabs(jacDet)*dV_ref[k];
877  ck.calculateG(jacInv,G,G_dd_G,tr_G);
878  //get the trial function gradients
879  ck.gradTrialFromRef(&u_grad_trial_ref[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
880  //get the solution
881  ck.valFromDOF(u_dof.data(),&u_l2g[eN_nDOF_trial_element],&u_trial_ref[k*nDOF_trial_element],u);
882  //get the solution gradients
883  ck.gradFromDOF(u_dof.data(),&u_l2g[eN_nDOF_trial_element],u_grad_trial,grad_u);
884  //precalculate test function products with integration weights
885  for (int j=0;j<nDOF_trial_element;j++)
886  {
887  u_test_dV[j] = u_test_ref[k*nDOF_trial_element+j]*dV;
888  for (int I=0;I<nSpace;I++)
889  {
890  u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
891  }
892  }
893  //
894  //calculate pde coefficients and derivatives at quadrature points
895  //
896  evaluateCoefficients(&velocity[eN_k_nSpace],
897  u,
898  grad_u,
899  m,
900  dm,
901  H,
902  dH);
903  //
904  //moving mesh
905  //
906  double mesh_velocity[3];
907  mesh_velocity[0] = xt;
908  mesh_velocity[1] = yt;
909  mesh_velocity[2] = zt;
910  for (int I=0;I<nSpace;I++)
911  {
912  f[I] = -MOVING_DOMAIN*m*mesh_velocity[I];
913  df[I] = -MOVING_DOMAIN*dm*mesh_velocity[I];
914  }
915  //
916  //calculate time derivatives
917  //
918  ck.bdf(alphaBDF,
919  q_m_betaBDF[eN_k],
920  m,
921  dm,
922  m_t,
923  dm_t);
924  //
925  //calculate subgrid error contribution to the Jacobian (strong residual, adjoint, jacobian of strong residual)
926  //
927  //calculate the adjoint times the test functions
928  for (int i=0;i<nDOF_test_element;i++)
929  {
930  //int eN_k_i_nSpace = (eN_k*nDOF_trial_element+i)*nSpace;
931  //Lstar_u[i]=ck.Hamiltonian_adjoint(dH,&u_grad_test_dV[eN_k_i_nSpace]);
932  register int i_nSpace = i*nSpace;
933  Lstar_u[i]=ck.Hamiltonian_adjoint(dH,&u_grad_test_dV[i_nSpace]) + MOVING_DOMAIN*ck.Advection_adjoint(df,&u_grad_test_dV[i_nSpace]);
934 
935  }
936  //calculate the Jacobian of strong residual
937  for (int j=0;j<nDOF_trial_element;j++)
938  {
939  //int eN_k_j=eN_k*nDOF_trial_element+j;
940  //int eN_k_j_nSpace = eN_k_j*nSpace;
941  int j_nSpace = j*nSpace;
942  dpdeResidual_u_u[j]=ck.MassJacobian_strong(dm_t,u_trial_ref[k*nDOF_trial_element+j]) +
943  ck.HamiltonianJacobian_strong(dH,&u_grad_trial[j_nSpace]) +
944  MOVING_DOMAIN*ck.AdvectionJacobian_strong(df,&u_grad_trial[j_nSpace]);
945 
946  }
947  //tau and tau*Res
948  double subgridErrorVelocity[nSpace];
949  for (int I=0;I<nSpace;I++)
950  subgridErrorVelocity[I] = dH[I] - MOVING_DOMAIN*df[I];
951 
952  calculateSubgridError_tau(elementDiameter[eN],
953  dm_t,
954  subgridErrorVelocity,//dH,
955  cfl[eN_k],
956  tau0);
957 
958  calculateSubgridError_tau(Ct_sge,
959  G,
960  dm_t,
961  subgridErrorVelocity,//dH,
962  tau1,
963  cfl[eN_k]);
964  tau = useMetrics*tau1+(1.0-useMetrics)*tau0;
965 
966  for(int j=0;j<nDOF_trial_element;j++)
967  dsubgridError_u_u[j] = -tau*dpdeResidual_u_u[j];
968  for(int i=0;i<nDOF_test_element;i++)
969  {
970  //int eN_k_i=eN_k*nDOF_test_element+i;
971  //int eN_k_i_nSpace=eN_k_i*nSpace;
972  for(int j=0;j<nDOF_trial_element;j++)
973  {
974  //int eN_k_j=eN_k*nDOF_trial_element+j;
975  //int eN_k_j_nSpace = eN_k_j*nSpace;
976  int j_nSpace = j*nSpace;
977  int i_nSpace = i*nSpace;
978  if (EXPLICIT_METHOD==1)
979  {
980  elementJacobian_u_u[i][j] +=
981  ck.MassJacobian_weak(1.0,
982  u_trial_ref[k*nDOF_trial_element+j],
983  u_test_dV[i]);
984  //ck.MassJacobian_weak(1.0, // semi-implicit Lax Wendroff
985  // u_trial_ref[k*nDOF_trial_element+j],
986  // u_test_dV[i]) +
987  //0.5*dt*dt*dV*
988  //ck.NumericalDiffusion(1.0,dH,&u_grad_trial[i_nSpace])*
989  //ck.NumericalDiffusion(1.0,dH,&u_grad_trial[j_nSpace]);
990  }
991  else // Implicit SUPG with SHOCK CAPTURING
992  elementJacobian_u_u[i][j] +=
993  ck.MassJacobian_weak(dm_t,
994  u_trial_ref[k*nDOF_trial_element+j],
995  u_test_dV[i]) +
996  ck.HamiltonianJacobian_weak(dH,&u_grad_trial[j_nSpace],u_test_dV[i]) +
997  MOVING_DOMAIN*ck.AdvectionJacobian_weak(df,
998  u_trial_ref[k*nDOF_trial_element+j],
999  &u_grad_test_dV[i_nSpace]) +
1000  (PURE_BDF == 1 ? 0. : 1.)*
1001  (ck.SubgridErrorJacobian(dsubgridError_u_u[j],Lstar_u[i]) +
1002  ck.NumericalDiffusionJacobian(q_numDiff_u_last[eN_k],
1003  &u_grad_trial[j_nSpace],
1004  &u_grad_test_dV[i_nSpace]));
1005  }//j
1006  }//i
1007  }//k
1008  //
1009  //load into element Jacobian into global Jacobian
1010  //
1011  for (int i=0;i<nDOF_test_element;i++)
1012  {
1013  int eN_i = eN*nDOF_test_element+i;
1014  for (int j=0;j<nDOF_trial_element;j++)
1015  {
1016  int eN_i_j = eN_i*nDOF_trial_element+j;
1017  globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_u_u[eN_i_j]] += elementJacobian_u_u[i][j];
1018  }//j
1019  }//i
1020  }//elements
1021  //
1022  //loop over exterior element boundaries to compute the surface integrals and load them into the global Jacobian
1023  //
1024  if (EXPLICIT_METHOD==0) //(semi)-implicit
1025  for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
1026  {
1027  register int ebN = exteriorElementBoundariesArray[ebNE];
1028  register int eN = elementBoundaryElementsArray[ebN*2+0],
1029  ebN_local = elementBoundaryLocalElementBoundariesArray[ebN*2+0],
1030  eN_nDOF_trial_element = eN*nDOF_trial_element;
1031  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
1032  {
1033  register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
1034  ebNE_kb_nSpace = ebNE_kb*nSpace,
1035  ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
1036  ebN_local_kb_nSpace = ebN_local_kb*nSpace;
1037 
1038  register double u_ext=0.0,
1039  grad_u_ext[nSpace],
1040  m_ext=0.0,
1041  dm_ext=0.0,
1042  H_ext=0.0,
1043  dH_ext[nSpace],
1044  bc_grad_u_ext[nSpace],
1045  bc_H_ext=0.0,
1046  bc_dH_ext[nSpace],
1047  //f_ext[nSpace],
1048  //df_ext[nSpace],
1049  dflux_u_u_ext=0.0,
1050  bc_u_ext=0.0,
1051  //bc_grad_u_ext[nSpace],
1052  bc_m_ext=0.0,
1053  bc_dm_ext=0.0,
1054  //bc_f_ext[nSpace],
1055  //bc_df_ext[nSpace],
1056  fluxJacobian_u_u[nDOF_trial_element],
1057  jac_ext[nSpace*nSpace],
1058  jacDet_ext,
1059  jacInv_ext[nSpace*nSpace],
1060  boundaryJac[nSpace*(nSpace-1)],
1061  metricTensor[(nSpace-1)*(nSpace-1)],
1062  metricTensorDetSqrt,
1063  dS,
1064  u_test_dS[nDOF_test_element],
1065  u_grad_trial_trace[nDOF_trial_element*nSpace],
1066  normal[nSpace],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
1067  G[nSpace*nSpace],G_dd_G,tr_G;
1068  //
1069  //calculate the solution and gradients at quadrature points
1070  //
1071  // u_ext=0.0;
1072  // for (int I=0;I<nSpace;I++)
1073  // {
1074  // grad_u_ext[I] = 0.0;
1075  // bc_grad_u_ext[I] = 0.0;
1076  // }
1077  // for (int j=0;j<nDOF_trial_element;j++)
1078  // {
1079  // register int eN_j = eN*nDOF_trial_element+j,
1080  // ebNE_kb_j = ebNE_kb*nDOF_trial_element+j,
1081  // ebNE_kb_j_nSpace= ebNE_kb_j*nSpace;
1082  // u_ext += valFromDOF_c(u_dof[u_l2g[eN_j]],u_trial_ext[ebNE_kb_j]);
1083 
1084  // for (int I=0;I<nSpace;I++)
1085  // {
1086  // grad_u_ext[I] += gradFromDOF_c(u_dof[u_l2g[eN_j]],u_grad_trial_ext[ebNE_kb_j_nSpace+I]);
1087  // }
1088  // }
1089  ck.calculateMapping_elementBoundary(eN,
1090  ebN_local,
1091  kb,
1092  ebN_local_kb,
1093  mesh_dof.data(),
1094  mesh_l2g.data(),
1095  mesh_trial_trace_ref.data(),
1096  mesh_grad_trial_trace_ref.data(),
1097  boundaryJac_ref.data(),
1098  jac_ext,
1099  jacDet_ext,
1100  jacInv_ext,
1101  boundaryJac,
1102  metricTensor,
1103  metricTensorDetSqrt,
1104  normal_ref.data(),
1105  normal,
1106  x_ext,y_ext,z_ext);
1107  ck.calculateMappingVelocity_elementBoundary(eN,
1108  ebN_local,
1109  kb,
1110  ebN_local_kb,
1111  mesh_velocity_dof.data(),
1112  mesh_l2g.data(),
1113  mesh_trial_trace_ref.data(),
1114  xt_ext,yt_ext,zt_ext,
1115  normal,
1116  boundaryJac,
1117  metricTensor,
1118  integralScaling);
1119  dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref[kb];
1120  //dS = metricTensorDetSqrt*dS_ref[kb];
1121  ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
1122  //compute shape and solution information
1123  //shape
1124  ck.gradTrialFromRef(&u_grad_trial_trace_ref[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,u_grad_trial_trace);
1125  //solution and gradients
1126  ck.valFromDOF(u_dof.data(),&u_l2g[eN_nDOF_trial_element],&u_trial_trace_ref[ebN_local_kb*nDOF_test_element],u_ext);
1127  ck.gradFromDOF(u_dof.data(),&u_l2g[eN_nDOF_trial_element],u_grad_trial_trace,grad_u_ext);
1128  //precalculate test function products with integration weights
1129  for (int j=0;j<nDOF_trial_element;j++)
1130  {
1131  u_test_dS[j] = u_test_trace_ref[ebN_local_kb*nDOF_test_element+j]*dS;
1132  }
1133  //
1134  //load the boundary values
1135  //
1136  bc_u_ext = isDOFBoundary_u[ebNE_kb]*ebqe_bc_u_ext[ebNE_kb]+(1-isDOFBoundary_u[ebNE_kb])*ebqe_rd_u_ext[ebNE_kb];
1137  //
1138  //calculate the internal and external trace of the pde coefficients
1139  //
1140  evaluateCoefficients(&ebqe_velocity_ext[ebNE_kb_nSpace],
1141  u_ext,
1142  grad_u_ext,
1143  m_ext,
1144  dm_ext,
1145  H_ext,
1146  dH_ext);
1147  evaluateCoefficients(&ebqe_velocity_ext[ebNE_kb_nSpace],
1148  bc_u_ext,
1149  bc_grad_u_ext,
1150  bc_m_ext,
1151  bc_dm_ext,
1152  bc_H_ext,
1153  bc_dH_ext);
1154  //
1155  //moving domain
1156  //
1157  double velocity_ext[nSpace];
1158  double mesh_velocity[3];
1159  mesh_velocity[0] = xt_ext;
1160  mesh_velocity[1] = yt_ext;
1161  mesh_velocity[2] = zt_ext;
1162  for (int I=0;I<nSpace;I++)
1163  {
1164  velocity_ext[I] = - MOVING_DOMAIN*mesh_velocity[I];
1165  }
1166  //
1167  //calculate the numerical fluxes
1168  //
1169  exteriorNumericalFluxDerivative(normal,
1170  dH_ext,
1171  velocity_ext,
1172  dflux_u_u_ext);
1173  //
1174  //calculate the flux jacobian
1175  //
1176  for (int j=0;j<nDOF_trial_element;j++)
1177  {
1178  //register int ebNE_kb_j = ebNE_kb*nDOF_trial_element+j;
1179  register int ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j;
1180 
1181  fluxJacobian_u_u[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_u_u_ext,u_trial_trace_ref[ebN_local_kb_j]);
1182  }//j
1183  //
1184  //update the global Jacobian from the flux Jacobian
1185  //
1186  for (int i=0;i<nDOF_test_element;i++)
1187  {
1188  register int eN_i = eN*nDOF_test_element+i;
1189  //register int ebNE_kb_i = ebNE_kb*nDOF_test_element+i;
1190  for (int j=0;j<nDOF_trial_element;j++)
1191  {
1192  register int ebN_i_j = ebN*4*nDOF_test_X_trial_element + i*nDOF_trial_element + j;
1193  globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] += fluxJacobian_u_u[j]*u_test_dS[i];
1194  //globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] += MOVING_DOMAIN*fluxJacobian_u_u[j]*u_test_dS[i];
1195  }//j
1196  }//i
1197  }//kb
1198  }//ebNE
1199  }//computeJacobian
1200 
1201  void calculateWaterline(arguments_dict& args)
1202  {
1203  xt::pyarray<int>& wlc = args.array<int>("wlc");
1204  xt::pyarray<double>& waterline = args.array<double>("waterline");
1205  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1206  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1207  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1208  xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
1209  double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
1210  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1211  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1212  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1213  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1214  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1215  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
1216  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
1217  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
1218  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
1219  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
1220  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
1221  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
1222  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
1223  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
1224  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
1225  int nElements_global = args.scalar<int>("nElements_global");
1226  double useMetrics = args.scalar<double>("useMetrics");
1227  double alphaBDF = args.scalar<double>("alphaBDF");
1228  int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
1229  double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
1230  double sc_uref = args.scalar<double>("sc_uref");
1231  double sc_alpha = args.scalar<double>("sc_alpha");
1232  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1233  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1234  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1235  xt::pyarray<double>& u_dof_old = args.array<double>("u_dof_old");
1236  xt::pyarray<double>& velocity = args.array<double>("velocity");
1237  xt::pyarray<double>& q_m = args.array<double>("q_m");
1238  xt::pyarray<double>& q_u = args.array<double>("q_u");
1239  xt::pyarray<double>& q_n = args.array<double>("q_n");
1240  xt::pyarray<double>& q_dH = args.array<double>("q_dH");
1241  xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
1242  xt::pyarray<double>& cfl = args.array<double>("cfl");
1243  xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
1244  xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
1245  int offset_u = args.scalar<int>("offset_u");
1246  int stride_u = args.scalar<int>("stride_u");
1247  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
1248  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
1249  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
1250  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
1251  xt::pyarray<int>& elementBoundaryMaterialTypes = args.array<int>("elementBoundaryMaterialTypes");
1252  xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
1253  xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
1254  xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
1255  xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
1256  // Tetrehedral elements specific extraction routine for waterline extraction
1257  // Loops over boundaries and checks if boundary is infact a hull (hardwired check if mattype > 6)
1258  // Extracts the nodal values of boundary triangle (4th point is dropped = hardwired assumption we are dealing with linear tet)
1259  // Then computes an average value and position for both negative and positive values
1260  // If both positive and negative values re found, and we are actually in a triangle containing the interface
1261  // a linear interpolation of negative and positive average is reported as interface location (exact in case of linear tets)
1262 
1263  for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
1264  {
1265  register int ebN = exteriorElementBoundariesArray[ebNE];
1266  register int eN = elementBoundaryElementsArray[ebN*2+0];
1267  register int bN = elementBoundaryLocalElementBoundariesArray[ebN*2+0];
1268 
1269  if (elementBoundaryMaterialTypes[ebN] >6)
1270  {
1271 
1272  double val,x,y,z;
1273  int pos=0, neg=0;
1274  double xn=0.0, yn=0.0, zn=0.0, vn=0.0;
1275  double xp=0.0, yp=0.0, zp=0.0, vp=0.0;
1276 
1277  for (int j=0;j<nDOF_trial_element;j++)
1278  {
1279  if (j != bN) {
1280  int eN_nDOF_trial_element_j = eN*nDOF_trial_element + j;
1281  int eN_nDOF_mesh_trial_element_j = eN*nDOF_mesh_trial_element + j;
1282  val = u_dof[u_l2g[eN_nDOF_trial_element_j]];
1283  x = mesh_dof[mesh_l2g[eN_nDOF_mesh_trial_element_j]*3+0];
1284  y = mesh_dof[mesh_l2g[eN_nDOF_mesh_trial_element_j]*3+1];
1285  z = mesh_dof[mesh_l2g[eN_nDOF_mesh_trial_element_j]*3+2];
1286 
1287  if (val < 0.0)
1288  {
1289  neg++;
1290  vn+=val;
1291  xn+=x;
1292  yn+=y;
1293  zn+=z;
1294  }
1295  else
1296  {
1297  pos++;
1298  vp+=val;
1299  xp+=x;
1300  yp+=y;
1301  zp+=z;
1302  }
1303  }
1304  } // trail for
1305 
1306 
1307  if ((pos > 0) && (neg > 0) )
1308  {
1309  vp /= pos;
1310  vn /= neg;
1311 
1312  double alpha = vp/(vp -vn);
1313 
1314  waterline[wlc[0]*3 + 0] = alpha*(xn/neg) + (1.0-alpha)*(xp/pos);
1315  waterline[wlc[0]*3 + 1] = alpha*(yn/neg) + (1.0-alpha)*(yp/pos);
1316  waterline[wlc[0]*3 + 2] = alpha*(zn/neg) + (1.0-alpha)*(zp/pos);
1317  wlc[0]++;
1318 
1319  } // end value if
1320 
1321  } // end bnd mat check
1322 
1323  }//ebNE
1324 
1325  //std::cout<<"CPP WLC "<<wlc[0]<<std::endl;
1326  } // calcWaterline
1327 
1328  };//cppNCLS3P
1329 
1330  inline cppNCLS3P_base* newNCLS3P(int nSpaceIn,
1331  int nQuadraturePoints_elementIn,
1332  int nDOF_mesh_trial_elementIn,
1333  int nDOF_trial_elementIn,
1334  int nDOF_test_elementIn,
1335  int nQuadraturePoints_elementBoundaryIn,
1336  int CompKernelFlag)
1337  {
1338  if (nSpaceIn == 2)
1339  return proteus::chooseAndAllocateDiscretization2D<cppNCLS3P_base,cppNCLS3P,CompKernel>(nSpaceIn,
1340  nQuadraturePoints_elementIn,
1341  nDOF_mesh_trial_elementIn,
1342  nDOF_trial_elementIn,
1343  nDOF_test_elementIn,
1344  nQuadraturePoints_elementBoundaryIn,
1345  CompKernelFlag);
1346  else
1347  return proteus::chooseAndAllocateDiscretization<cppNCLS3P_base,cppNCLS3P,CompKernel>(nSpaceIn,
1348  nQuadraturePoints_elementIn,
1349  nDOF_mesh_trial_elementIn,
1350  nDOF_trial_elementIn,
1351  nDOF_test_elementIn,
1352  nQuadraturePoints_elementBoundaryIn,
1353  CompKernelFlag);
1354  /* return proteus::chooseAndAllocateDiscretization<cppNCLS3P_base,cppNCLS3P>(nSpaceIn, */
1355  /* nQuadraturePoints_elementIn, */
1356  /* nDOF_mesh_trial_elementIn, */
1357  /* nDOF_trial_elementIn, */
1358  /* nDOF_test_elementIn, */
1359  /* nQuadraturePoints_elementBoundaryIn, */
1360  /* CompKernelFlag); */
1361  }
1362 }//proteus
1363 #endif
sign
#define sign(x, y)
Definition: jf.h:44
proteus::cppNCLS3P::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: NCLS3P.h:107
proteus::cppNCLS3P::calculateWaterline
void calculateWaterline(arguments_dict &args)
Definition: NCLS3P.h:1201
proteus::cppNCLS3P::ck
CompKernelType ck
Definition: NCLS3P.h:49
proteus::cppNCLS3P::calculateCFL
void calculateCFL(const double &elementDiameter, const double df[nSpace], double &cfl)
Definition: NCLS3P.h:56
proteus::cppNCLS3P::calculateResidual
void calculateResidual(arguments_dict &args)
Definition: NCLS3P.h:171
neg
double neg(double a)
Definition: testFMMandFSW.cpp:10
proteus::cppNCLS3P::exteriorNumericalFlux
void exteriorNumericalFlux(const double n[nSpace], const double &bc_u, const double &u, const double velocity[nSpace], const double velocity_movingDomain[nSpace], double &flux)
Definition: NCLS3P.h:122
proteus::cppNCLS3P::nDOF_test_X_trial_element
const int nDOF_test_X_trial_element
Definition: NCLS3P.h:48
proteus::DENTROPY
double DENTROPY(const double &phi, const double &dummyL, const double &dummyR)
Definition: NCLS.h:25
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::cppNCLS3P_base::~cppNCLS3P_base
virtual ~cppNCLS3P_base()
Definition: NCLS3P.h:32
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)
proteus::cppNCLS3P::calculateSubgridError_tau
void calculateSubgridError_tau(const double &elementDiameter, const double &dmt, const double dH[nSpace], double &cfl, double &tau)
Definition: NCLS3P.h:88
H
Double H
Definition: Headers.h:65
v
Double v
Definition: Headers.h:95
proteus::cppNCLS3P_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)=0
proteus::cppNCLS3P::evaluateCoefficients
void evaluateCoefficients(const double v[nSpace], const double &u, const double grad_u[nSpace], double &m, double &dm, double &H, double dH[nSpace])
Definition: NCLS3P.h:69
proteus::ENTROPY
double ENTROPY(const double &g, const double &h, const double &hu, const double &hv, const double &z, const double &one_over_hReg)
Definition: GN_SW2DCV.h:50
proteus::cppNCLS3P::calculateJacobian
void calculateJacobian(arguments_dict &args)
Definition: NCLS3P.h:768
proteus::newNCLS3P
cppNCLS3P_base * newNCLS3P(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition: NCLS3P.h:1330
z
Double * z
Definition: Headers.h:49
proteus::cppNCLS3P
Definition: NCLS3P.h:46
u
Double u
Definition: Headers.h:89
proteus::cppNCLS3P_base::calculateWaterline
virtual void calculateWaterline(arguments_dict &args)=0
xt
Definition: AddedMass.cpp:7
cMax
#define cMax
Definition: NCLS3P.h:11
proteus::cppNCLS3P_base
Definition: NCLS3P.h:29
proteus::cppNCLS3P_base::calculateResidual
virtual void calculateResidual(arguments_dict &args)=0
proteus
Definition: ADR.h:17
proteus::Sign
double Sign(const double &z)
Definition: CLSVOF.h:23
proteus::f
double f(const double &g, const double &h, const double &hZ)
Definition: SW2DCV.h:58
proteus::cppNCLS3P::cppNCLS3P
cppNCLS3P()
Definition: NCLS3P.h:50
proteus::cppNCLS3P::exteriorNumericalFluxDerivative
void exteriorNumericalFluxDerivative(const double n[nSpace], const double velocity[nSpace], const double velocity_movingDomain[nSpace], double &dflux)
Definition: NCLS3P.h:149
proteus::arguments_dict
Definition: ArgumentsDict.h:70
pos
double pos(double a)
Definition: testFMMandFSW.cpp:8
ArgumentsDict.h
cE
#define cE
Definition: NCLS3P.h:10
Generated on Fri Jul 1 2022 11:17:39 for proteus by doxygen 1.8.20