proteus  1.8.1
C/C++/Fortran libraries
NCLS.h
Go to the documentation of this file.
1 #ifndef NCLS_H
2 #define NCLS_H
3 #include <cmath>
4 #include <iostream>
5 #include <valarray>
6 #include "CompKernel.h"
7 #include "ModelFactory.h"
8 #include "ArgumentsDict.h"
9 #include "xtensor-python/pyarray.hpp"
10 
11 namespace py = pybind11;
12 
13 #define POWER_SMOOTHNESS_INDICATOR 2
14 #define IS_BETAij_ONE 0
15 
16 namespace proteus
17 {
18  inline double Sign(const double& z){
19  return (z >= 0.0 ? 1.0 : -1.0);
20  }
21  // Power entropy //
22  inline double ENTROPY(const double& phi, const double& dummyL, const double& dummyR){
23  return 1./2.*std::pow(fabs(phi),2.);
24  }
25  inline double DENTROPY(const double& phi, const double& dummyL, const double& dummyR){
26  return fabs(phi)*(phi >= 0. ? 1. : -1.);
27  }
28  // Log entropy //
29  inline double ENTROPY_LOG(const double& phi, const double& phiL, const double& phiR){
30  return std::log(fabs((phi-phiL)*(phiR-phi))+1E-14);
31  }
32  inline double DENTROPY_LOG(const double& phi, const double& phiL, const double& phiR){
33  return (phiL+phiR-2*phi)*((phi-phiL)*(phiR-phi)>=0 ? 1 : -1)/(fabs((phi-phiL)*(phiR-phi))+1E-14);
34  }
35 }
36 
37 namespace proteus
38 {
39  class NCLS_base
40  {
41  //The base class defining the interface
42  public:
43  std::valarray<double> L2_norm_per_node;
44  std::valarray<double> TransportMatrix, TransposeTransportMatrix;
45  std::valarray<double> psi, etaMax, etaMin;
46  std::valarray<double> global_entropy_residual;
47  std::valarray<double> gx, gy, gz, eta,
48  alpha_numerator_pos, alpha_numerator_neg, alpha_denominator_pos, alpha_denominator_neg;
49  virtual ~NCLS_base(){}
50  virtual void calculateResidual(arguments_dict& args)=0;
51  virtual void calculateJacobian(arguments_dict& args)=0;
52  virtual void calculateWaterline(arguments_dict& args)=0;
53  virtual double calculateRedistancingResidual(arguments_dict& args)=0;
54  virtual double calculateRhsSmoothing(arguments_dict& args)=0;
55  virtual void calculateResidual_entropy_viscosity(arguments_dict& args) = 0;
56  virtual void calculateMassMatrix(arguments_dict& args)=0;
57  virtual void calculateSmoothingMatrix(arguments_dict& args) = 0;
58  };
59 
60  template<class CompKernelType,
61  int nSpace,
62  int nQuadraturePoints_element,
63  int nDOF_mesh_trial_element,
64  int nDOF_trial_element,
65  int nDOF_test_element,
66  int nQuadraturePoints_elementBoundary>
67  class NCLS : public NCLS_base
68  {
69  public:
70  const int nDOF_test_X_trial_element;
71  CompKernelType ck;
72  NCLS():
73  nDOF_test_X_trial_element(nDOF_test_element*nDOF_trial_element),
74  ck()
75  {}
76 
77  inline void evaluateCoefficients(const double v[nSpace],
78  const double& u,
79  const double grad_u[nSpace],
80  double& m,
81  double& dm,
82  double& H,
83  double dH[nSpace])
84  {
85  m = u;
86  dm=1.0;
87  H = 0.0;
88  for (int I=0; I < nSpace; I++)
89  {
90  H += v[I]*grad_u[I];
91  dH[I] = v[I];
92  }
93  }
94 
95  inline
96  void calculateCFL(const double& elementDiameter,
97  const double df[nSpace],
98  double& cfl)
99  {
100  double h,nrm_v;
101  h = elementDiameter;
102  nrm_v=0.0;
103  for(int I=0;I<nSpace;I++)
104  nrm_v+=df[I]*df[I];
105  nrm_v = sqrt(nrm_v);
106  cfl = nrm_v/h;
107  }
108 
109  inline
110  void calculateSubgridError_tau(const double& elementDiameter,
111  const double& dmt,
112  const double dH[nSpace],
113  double& cfl,
114  double& tau)
115  {
116  double h,nrm_v,oneByAbsdt;
117  h = elementDiameter;
118  nrm_v=0.0;
119  for(int I=0;I<nSpace;I++)
120  nrm_v+=dH[I]*dH[I];
121  nrm_v = sqrt(nrm_v);
122  cfl = nrm_v/h;
123  oneByAbsdt = fabs(dmt);
124  tau = 1.0/(2.0*nrm_v/h + oneByAbsdt + 1.0e-8);
125  }
126 
127 
128  inline
129  void calculateSubgridError_tau( const double& Ct_sge,
130  const double G[nSpace*nSpace],
131  const double& A0,
132  const double Ai[nSpace],
133  double& tau_v,
134  double& cfl)
135  {
136  double v_d_Gv=0.0;
137  for(int I=0;I<nSpace;I++)
138  for (int J=0;J<nSpace;J++)
139  v_d_Gv += Ai[I]*G[I*nSpace+J]*Ai[J];
140 
141  tau_v = 1.0/sqrt(Ct_sge*A0*A0 + v_d_Gv + 1.0e-8);
142  }
143 
144  void exteriorNumericalFlux(const double n[nSpace],
145  const double& bc_u,
146  const double& u,
147  const double velocity[nSpace],
148  const double velocity_movingDomain[nSpace],
149  double& flux)
150  {
151  double flow_total=0.0,flow_fluid=0.0,flow_movingDomain=0.0;
152  for (int I=0; I < nSpace; I++)
153  {
154  flow_fluid += n[I]*velocity[I];
155  //flow_movingDomain -= n[I]*velocity_movingDomain[I];
156  }
157  flow_total = flow_fluid+flow_movingDomain;
158  if (flow_total > 0.0)
159  {
160  flux = u*flow_movingDomain;
161  }
162  else
163  {
164  flux = bc_u*flow_movingDomain - flow_fluid*(u-bc_u);
165  }
166  }
167 
168  inline
169  void exteriorNumericalFluxDerivative(const double n[nSpace],
170  const double velocity[nSpace],
171  const double velocity_movingDomain[nSpace],
172  double& dflux)
173  {
174  double flow_total=0.0,flow_fluid=0.0,flow_movingDomain=0.0;
175  for (int I=0; I < nSpace; I++)
176  {
177  flow_fluid += n[I]*velocity[I];
178  //flow_movingDomain -= n[I]*velocity_movingDomain[I];
179  }
180  flow_total=flow_fluid+flow_movingDomain;
181  if (flow_total > 0.0)
182  {
183  dflux = flow_movingDomain;
184  }
185  else
186  {
187  dflux = -flow_fluid;
188  }
189  }
190 
191  inline
192  double sign(const double phi,
193  const double eps) //eps=epsFactRedistancing=0.33*he (for instance)
194  {
195  double H;
196  if (phi > eps)
197  H = 1.0;
198  else if (phi < -eps)
199  H = 0.0;
200  else if (phi == 0.0)
201  H = 0.5;
202  else
203  H = 0.5*(1.0 + phi/eps + std::sin(M_PI*phi/eps)/M_PI);
204  return -1.0 + 2.0*H;
205  //return (u > 0 ? 1. : -1.)*(u ==0 ? 0. : 1.);
206  //double tol_sign 0.1;
207  //return (u > tol_sign*epsFactRedistancing ? 1. : -1.)*((u > -tol_sign*epsFactRedistancing && u < tol_sign*epsFactRedistancing) ? 0. : 1.);
208  }
209 
210  void calculateResidual(arguments_dict& args)
211  {
212  double dt = args.scalar<double>("dt");
213  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
214  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
215  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
216  xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
217  double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
218  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
219  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
220  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
221  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
222  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
223  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
224  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
225  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
226  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
227  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
228  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
229  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
230  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
231  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
232  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
233  int nElements_global = args.scalar<int>("nElements_global");
234  double useMetrics = args.scalar<double>("useMetrics");
235  double alphaBDF = args.scalar<double>("alphaBDF");
236  int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
237  double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
238  double sc_uref = args.scalar<double>("sc_uref");
239  double sc_alpha = args.scalar<double>("sc_alpha");
240  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
241  xt::pyarray<int>& r_l2g = args.array<int>("r_l2g");
242  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
243  xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
244  int degree_polynomial = args.scalar<int>("degree_polynomial");
245  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
246  xt::pyarray<double>& u_dof_old = args.array<double>("u_dof_old");
247  xt::pyarray<double>& uStar_dof = args.array<double>("uStar_dof");
248  xt::pyarray<double>& velocity = args.array<double>("velocity");
249  xt::pyarray<double>& q_m = args.array<double>("q_m");
250  xt::pyarray<double>& q_u = args.array<double>("q_u");
251  xt::pyarray<double>& q_n = args.array<double>("q_n");
252  xt::pyarray<double>& q_dH = args.array<double>("q_dH");
253  xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
254  xt::pyarray<double>& q_dV = args.array<double>("q_dV");
255  xt::pyarray<double>& q_dV_last = args.array<double>("q_dV_last");
256  xt::pyarray<double>& cfl = args.array<double>("cfl");
257  xt::pyarray<double>& edge_based_cfl = args.array<double>("edge_based_cfl");
258  xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
259  xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
260  int offset_u = args.scalar<int>("offset_u");
261  int stride_u = args.scalar<int>("stride_u");
262  xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
263  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
264  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
265  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
266  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
267  xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
268  xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
269  xt::pyarray<double>& ebqe_rd_u_ext = args.array<double>("ebqe_rd_u_ext");
270  xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
271  xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
272  xt::pyarray<double>& ebqe_grad_u = args.array<double>("ebqe_grad_u");
273  xt::pyarray<double>& interface_locator = args.array<double>("interface_locator");
274  int PURE_BDF = args.scalar<int>("PURE_BDF");
275  int numDOFs = args.scalar<int>("numDOFs");
276  int NNZ = args.scalar<int>("NNZ");
277  xt::pyarray<int>& csrRowIndeces_DofLoops = args.array<int>("csrRowIndeces_DofLoops");
278  xt::pyarray<int>& csrColumnOffsets_DofLoops = args.array<int>("csrColumnOffsets_DofLoops");
279  xt::pyarray<int>& csrRowIndeces_CellLoops = args.array<int>("csrRowIndeces_CellLoops");
280  xt::pyarray<int>& csrColumnOffsets_CellLoops = args.array<int>("csrColumnOffsets_CellLoops");
281  xt::pyarray<int>& csrColumnOffsets_eb_CellLoops = args.array<int>("csrColumnOffsets_eb_CellLoops");
282  int LUMPED_MASS_MATRIX = args.scalar<int>("LUMPED_MASS_MATRIX");
283  xt::pyarray<double>& quantDOFs = args.array<double>("quantDOFs");
284  double lambda_coupez = args.scalar<double>("lambda_coupez");
285  double epsCoupez = args.scalar<double>("epsCoupez");
286  double epsFactRedistancing = args.scalar<double>("epsFactRedistancing");
287  int COUPEZ = args.scalar<int>("COUPEZ");
288  int SATURATED_LEVEL_SET = args.scalar<int>("SATURATED_LEVEL_SET");
289  xt::pyarray<double>& Cx = args.array<double>("Cx");
290  xt::pyarray<double>& Cy = args.array<double>("Cy");
291  xt::pyarray<double>& Cz = args.array<double>("Cz");
292  xt::pyarray<double>& ML = args.array<double>("ML");
293  int STABILIZATION_TYPE = args.scalar<int>("STABILIZATION_TYPE");
294  int ENTROPY_TYPE = args.scalar<int>("ENTROPY_TYPE");
295  double cE = args.scalar<double>("cE");
296  //cek should this be read in?
297  double Ct_sge = 4.0;
298 
299  //loop over elements to compute volume integrals and load them into element and global residual
300  //
301  //eN is the element index
302  //eN_k is the quadrature point index for a scalar
303  //eN_k_nSpace is the quadrature point index for a vector
304  //eN_i is the element test function index
305  //eN_j is the element trial function index
306  //eN_k_j is the quadrature point index for a trial function
307  //eN_k_i is the quadrature point index for a trial function
308  for(int eN=0;eN<nElements_global;eN++)
309  {
310  //declare local storage for element residual and initialize
311  register double elementResidual_u[nDOF_test_element];
312  for (int i=0;i<nDOF_test_element;i++)
313  {
314  elementResidual_u[i]=0.0;
315  }//i
316  //loop over quadrature points and compute integrands
317  for (int k=0;k<nQuadraturePoints_element;k++)
318  {
319  //compute indeces and declare local storage
320  register int eN_k = eN*nQuadraturePoints_element+k,
321  eN_k_nSpace = eN_k*nSpace,
322  eN_nDOF_trial_element = eN*nDOF_trial_element;
323  register double u=0.0, u_old=0.0, grad_u[nSpace], grad_u_old[nSpace],
324  m=0.0, dm=0.0,
325  H=0.0, dH[nSpace],
326  f[nSpace],df[nSpace],//for MOVING_DOMAIN
327  m_t=0.0,dm_t=0.0,
328  pdeResidual_u=0.0,
329  Lstar_u[nDOF_test_element],
330  subgridError_u=0.0,
331  tau=0.0,tau0=0.0,tau1=0.0,
332  numDiff0=0.0,numDiff1=0.0,
333  jac[nSpace*nSpace],
334  jacDet,
335  jacInv[nSpace*nSpace],
336  u_grad_trial[nDOF_trial_element*nSpace],
337  u_test_dV[nDOF_trial_element],
338  u_grad_test_dV[nDOF_test_element*nSpace],
339  dV,x,y,z,xt,yt,zt,
340  G[nSpace*nSpace],G_dd_G,tr_G;//,norm_Rv;
341  //
342  //compute solution and gradients at quadrature points
343  //
344  ck.calculateMapping_element(eN,
345  k,
346  mesh_dof.data(),
347  mesh_l2g.data(),
348  mesh_trial_ref.data(),
349  mesh_grad_trial_ref.data(),
350  jac,
351  jacDet,
352  jacInv,
353  x,y,z);
354  ck.calculateMappingVelocity_element(eN,
355  k,
356  mesh_velocity_dof.data(),
357  mesh_l2g.data(),
358  mesh_trial_ref.data(),
359  xt,yt,zt);
360  //get the physical integration weight
361  dV = fabs(jacDet)*dV_ref.data()[k];
362  ck.calculateG(jacInv,G,G_dd_G,tr_G);
363  //get the trial function gradients
364  ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
365  //get the solution
366  ck.valFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],u);
367  //get the solution gradients
368  ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial,grad_u);
369  //precalculate test function products with integration weights
370  for (int j=0;j<nDOF_trial_element;j++)
371  {
372  u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
373  for (int I=0;I<nSpace;I++)
374  {
375  u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
376  }
377  }
378  //save solution at quadrature points for other models to use
379  q_u.data()[eN_k]=u;
380  for (int I=0;I<nSpace;I++)
381  q_n.data()[eN_k_nSpace+I] = grad_u[I];
382  //
383  //calculate pde coefficients at quadrature points
384  //
385  evaluateCoefficients(&velocity.data()[eN_k_nSpace],
386  u,
387  grad_u,
388  m,
389  dm,
390  H,
391  dH);
392  //
393  //moving mesh
394  //
395  double mesh_velocity[3];
396  mesh_velocity[0] = xt;
397  mesh_velocity[1] = yt;
398  mesh_velocity[2] = zt;
399  for (int I=0;I<nSpace;I++)
400  {
401  f[I] = -MOVING_DOMAIN*m*mesh_velocity[I];
402  df[I] = -MOVING_DOMAIN*dm*mesh_velocity[I];
403  }
404  //
405  //calculate time derivative at quadrature points
406  //
407  if (q_dV_last.data()[eN_k] <= -100)
408  q_dV_last.data()[eN_k] = dV;
409  q_dV.data()[eN_k] = dV;
410  ck.bdf(alphaBDF,
411  q_m_betaBDF.data()[eN_k],
412  m,
413  dm,
414  m_t,
415  dm_t);
416  //
417  //calculate subgrid error (strong residual and adjoint)
418  //
419  //calculate strong residual
420  pdeResidual_u =
421  ck.Mass_strong(m_t) +
422  ck.Hamiltonian_strong(dH,grad_u)+
423  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
424 
425  //calculate adjoint
426  for (int i=0;i<nDOF_test_element;i++)
427  {
428  //register int eN_k_i_nSpace = (eN_k*nDOF_trial_element+i)*nSpace;
429  //Lstar_u[i] = ck.Hamiltonian_adjoint(dH,&u_grad_test_dV[eN_k_i_nSpace]);
430  register int i_nSpace = i*nSpace;
431  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]);
432  }
433  //calculate tau and tau*Res
434  double subgridErrorVelocity[nSpace];
435  for (int I=0;I<nSpace;I++)
436  subgridErrorVelocity[I] = dH[I] - MOVING_DOMAIN*df[I];
437 
438  calculateSubgridError_tau(elementDiameter.data()[eN],
439  dm_t,
440  subgridErrorVelocity,//dH,
441  cfl.data()[eN_k],
442  tau0);
443 
444  calculateSubgridError_tau(Ct_sge,
445  G,
446  dm_t,
447  subgridErrorVelocity,//dH,
448  tau1,
449  cfl.data()[eN_k]);
450 
451  tau = useMetrics*tau1+(1.0-useMetrics)*tau0;
452 
453  subgridError_u = -tau*pdeResidual_u;
454  //
455  //calculate shock capturing diffusion
456  //
457  ck.calculateNumericalDiffusion(shockCapturingDiffusion,elementDiameter.data()[eN],pdeResidual_u,grad_u,numDiff0);
458  ck.calculateNumericalDiffusion(shockCapturingDiffusion,sc_uref, sc_alpha,G,G_dd_G,pdeResidual_u,grad_u,numDiff1);
459  q_numDiff_u.data()[eN_k] = useMetrics*numDiff1+(1.0-useMetrics)*numDiff0;
460 
461  //
462  //update element residual
463  //
464  double sign;
465  int same_sign=1;
466  for(int i=0;i<nDOF_test_element;i++)
467  {
468  int gi = offset_u+stride_u*u_l2g.data()[eN*nDOF_test_element+i];
469  if (i==0)
470  sign = Sign(u_dof.data()[gi]);
471  else if (same_sign==1)
472  {
473  same_sign = sign == Sign(u_dof.data()[gi]) ? 1 : 0;
474  sign = Sign(u_dof.data()[gi]);
475  }
476  //register int eN_k_i=eN_k*nDOF_test_element+i,
477  // eN_k_i_nSpace = eN_k_i*nSpace;
478  register int i_nSpace=i*nSpace;
479 
480  elementResidual_u[i] +=
481  ck.Mass_weak(m_t,u_test_dV[i]) +
482  ck.Hamiltonian_weak(H,u_test_dV[i]) +
483  MOVING_DOMAIN*ck.Advection_weak(f,&u_grad_test_dV[i_nSpace])+
484  (PURE_BDF == 1 ? 0. : 1.)*
485  (ck.SubgridError(subgridError_u,Lstar_u[i]) +
486  ck.NumericalDiffusion(q_numDiff_u_last.data()[eN_k],
487  grad_u,
488  &u_grad_test_dV[i_nSpace]));
489  }//i
490  if (same_sign == 0) // This cell contains the interface
491  for(int i=0;i<nDOF_test_element;i++)
492  {
493  int gi = offset_u+stride_u*u_l2g.data()[eN*nDOF_test_element+i];
494  interface_locator.data()[gi] = 1.0;
495  }
496  //
497  //cek/ido todo, get rid of m, since u=m
498  //save momentum for time history and velocity for subgrid error
499  //save solution for other models
500  //
501  q_m.data()[eN_k] = m;
502  q_u.data()[eN_k] = u;
503  for (int I=0;I<nSpace;I++)
504  {
505  int eN_k_I = eN_k*nSpace+I;
506  q_dH.data()[eN_k_I] = dH[I];
507  }
508  }
509  //
510  //load element into global residual and save element residual
511  //
512  for(int i=0;i<nDOF_test_element;i++)
513  {
514  register int eN_i=eN*nDOF_test_element+i;
515  globalResidual.data()[offset_u+stride_u*r_l2g.data()[eN_i]] += elementResidual_u[i];
516  }//i
517  }//elements
518  //
519  //loop over exterior element boundaries to calculate surface integrals and load into element and global residuals
520  //
521  //ebNE is the Exterior element boundary INdex
522  //ebN is the element boundary INdex
523  //eN is the element index
524  //std::cout <<nExteriorElementBoundaries_global<<std::endl;
525  for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
526  {
527  register int ebN = exteriorElementBoundariesArray.data()[ebNE],
528  eN = elementBoundaryElementsArray.data()[ebN*2+0],
529  ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
530  eN_nDOF_trial_element = eN*nDOF_trial_element;
531  register double elementResidual_u[nDOF_test_element];
532  for (int i=0;i<nDOF_test_element;i++)
533  {
534  elementResidual_u[i]=0.0;
535  }
536  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
537  {
538  register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
539  ebNE_kb_nSpace = ebNE_kb*nSpace,
540  ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
541  ebN_local_kb_nSpace = ebN_local_kb*nSpace;
542  register double u_ext=0.0,
543  grad_u_ext[nSpace],
544  m_ext=0.0,
545  dm_ext=0.0,
546  H_ext=0.0,
547  dH_ext[nSpace],
548  //f_ext[nSpace],//MOVING_DOMAIN
549  //df_ext[nSpace],//MOVING_DOMAIN
550  //flux_ext=0.0,
551  bc_u_ext=0.0,
552  bc_grad_u_ext[nSpace],
553  bc_m_ext=0.0,
554  bc_dm_ext=0.0,
555  bc_H_ext=0.0,
556  bc_dH_ext[nSpace],
557  jac_ext[nSpace*nSpace],
558  jacDet_ext,
559  jacInv_ext[nSpace*nSpace],
560  boundaryJac[nSpace*(nSpace-1)],
561  metricTensor[(nSpace-1)*(nSpace-1)],
562  metricTensorDetSqrt,
563  dS,
564  u_test_dS[nDOF_test_element],
565  u_grad_trial_trace[nDOF_trial_element*nSpace],
566  normal[nSpace],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
567  G[nSpace*nSpace],G_dd_G,tr_G,flux_ext;
568  //
569  //calculate the solution and gradients at quadrature points
570  //
571  //compute information about mapping from reference element to physical element
572  ck.calculateMapping_elementBoundary(eN,
573  ebN_local,
574  kb,
575  ebN_local_kb,
576  mesh_dof.data(),
577  mesh_l2g.data(),
578  mesh_trial_trace_ref.data(),
579  mesh_grad_trial_trace_ref.data(),
580  boundaryJac_ref.data(),
581  jac_ext,
582  jacDet_ext,
583  jacInv_ext,
584  boundaryJac,
585  metricTensor,
586  metricTensorDetSqrt,
587  normal_ref.data(),
588  normal,
589  x_ext,y_ext,z_ext);
590  ck.calculateMappingVelocity_elementBoundary(eN,
591  ebN_local,
592  kb,
593  ebN_local_kb,
594  mesh_velocity_dof.data(),
595  mesh_l2g.data(),
596  mesh_trial_trace_ref.data(),
597  xt_ext,yt_ext,zt_ext,
598  normal,
599  boundaryJac,
600  metricTensor,
601  integralScaling);
602  dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref.data()[kb];
603  //get the metric tensor
604  //cek todo use symmetry
605  ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
606  //compute shape and solution information
607  //shape
608  ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,u_grad_trial_trace);
609  //solution and gradients
610  ck.valFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],u_ext);
611  ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial_trace,grad_u_ext);
612  //precalculate test function products with integration weights
613  for (int j=0;j<nDOF_trial_element;j++)
614  {
615  u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
616  }
617  //
618  //load the boundary values
619  //
620  bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+(1-isDOFBoundary_u.data()[ebNE_kb])*ebqe_rd_u_ext.data()[ebNE_kb];
621  //
622  //calculate the pde coefficients using the solution and the boundary values for the solution
623  //
624  evaluateCoefficients(&ebqe_velocity_ext.data()[ebNE_kb_nSpace],
625  u_ext,
626  grad_u_ext,
627  m_ext,
628  dm_ext,
629  H_ext,
630  dH_ext);
631  evaluateCoefficients(&ebqe_velocity_ext.data()[ebNE_kb_nSpace],
632  bc_u_ext,
633  bc_grad_u_ext,
634  bc_m_ext,
635  bc_dm_ext,
636  bc_H_ext,
637  bc_dH_ext);
638  //
639  //moving mesh
640  //
641  double velocity_ext[nSpace];
642  double mesh_velocity[3];
643  mesh_velocity[0] = xt_ext;
644  mesh_velocity[1] = yt_ext;
645  mesh_velocity[2] = zt_ext;
646  for (int I=0;I<nSpace;I++)
647  velocity_ext[I] = - MOVING_DOMAIN*mesh_velocity[I];
648  //
649  //calculate the numerical fluxes
650  //
651  exteriorNumericalFlux(normal,
652  bc_u_ext,
653  u_ext,
654  dH_ext,
655  velocity_ext,
656  flux_ext);
657  ebqe_u.data()[ebNE_kb] = u_ext;
658 
659  // gradient //
660 
661  register double norm = 1.0e-8;
662  for (int I=0;I<nSpace;I++)
663  norm += grad_u_ext[I]*grad_u_ext[I];
664  norm = sqrt(norm);
665  for (int I=0;I<nSpace;I++)
666  ebqe_grad_u.data()[ebNE_kb_nSpace+I] = grad_u_ext[I]/norm;
667 
668  //
669  //update residuals
670  //
671  for (int i=0;i<nDOF_test_element;i++)
672  {
673  //int ebNE_kb_i = ebNE_kb*nDOF_test_element+i;
674  elementResidual_u[i] += ck.ExteriorElementBoundaryFlux(flux_ext,u_test_dS[i]);
675  }//i
676  }//kb
677  //
678  //update the element and global residual storage
679  //
680  for (int i=0;i<nDOF_test_element;i++)
681  {
682  int eN_i = eN*nDOF_test_element+i;
683  //globalResidual.data()[offset_u+stride_u*u_l2g.data()[eN_i]] += MOVING_DOMAIN*elementResidual_u[i];
684  globalResidual.data()[offset_u+stride_u*r_l2g.data()[eN_i]] += elementResidual_u[i];
685 
686  }//i
687  }//ebNE
688  }
689 
690  void calculateJacobian(arguments_dict& args)
691  {
692  double dt = args.scalar<double>("dt");
693  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
694  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
695  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
696  xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
697  double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
698  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
699  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
700  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
701  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
702  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
703  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
704  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
705  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
706  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
707  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
708  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
709  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
710  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
711  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
712  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
713  int nElements_global = args.scalar<int>("nElements_global");
714  double useMetrics = args.scalar<double>("useMetrics");
715  double alphaBDF = args.scalar<double>("alphaBDF");
716  int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
717  double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
718  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
719  xt::pyarray<int>& r_l2g = args.array<int>("r_l2g");
720  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
721  int degree_polynomial = args.scalar<int>("degree_polynomial");
722  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
723  xt::pyarray<double>& velocity = args.array<double>("velocity");
724  xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
725  xt::pyarray<double>& cfl = args.array<double>("cfl");
726  xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
727  xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
728  xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
729  xt::pyarray<double>& globalJacobian = args.array<double>("globalJacobian");
730  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
731  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
732  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
733  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
734  xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
735  xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
736  xt::pyarray<double>& ebqe_rd_u_ext = args.array<double>("ebqe_rd_u_ext");
737  xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
738  xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
739  int PURE_BDF = args.scalar<int>("PURE_BDF");
740  int LUMPED_MASS_MATRIX = args.scalar<int>("LUMPED_MASS_MATRIX");
741  double Ct_sge = 4.0;
742 
743  //
744  //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
745  //
746  for(int eN=0;eN<nElements_global;eN++)
747  {
748  register double elementJacobian_u_u[nDOF_test_element][nDOF_trial_element];
749  for (int i=0;i<nDOF_test_element;i++)
750  for (int j=0;j<nDOF_trial_element;j++)
751  {
752  elementJacobian_u_u[i][j]=0.0;
753  }
754  for (int k=0;k<nQuadraturePoints_element;k++)
755  {
756  int eN_k = eN*nQuadraturePoints_element+k, //index to a scalar at a quadrature point
757  eN_k_nSpace = eN_k*nSpace,
758  eN_nDOF_trial_element = eN*nDOF_trial_element; //index to a vector at a quadrature point
759 
760  //declare local storage
761  register double u=0.0,
762  grad_u[nSpace],
763  m=0.0,dm=0.0,
764  H=0.0,dH[nSpace],
765  f[nSpace],df[nSpace],//MOVING_MESH
766  m_t=0.0,dm_t=0.0,
767  dpdeResidual_u_u[nDOF_trial_element],
768  Lstar_u[nDOF_test_element],
769  dsubgridError_u_u[nDOF_trial_element],
770  tau=0.0,tau0=0.0,tau1=0.0,
771  jac[nSpace*nSpace],
772  jacDet,
773  jacInv[nSpace*nSpace],
774  u_grad_trial[nDOF_trial_element*nSpace],
775  dV,
776  u_test_dV[nDOF_test_element],
777  u_grad_test_dV[nDOF_test_element*nSpace],
778  x,y,z,xt,yt,zt,
779  G[nSpace*nSpace],G_dd_G,tr_G;
780  //
781  //calculate solution and gradients at quadrature points
782  //
783  //get jacobian, etc for mapping reference element
784  ck.calculateMapping_element(eN,
785  k,
786  mesh_dof.data(),
787  mesh_l2g.data(),
788  mesh_trial_ref.data(),
789  mesh_grad_trial_ref.data(),
790  jac,
791  jacDet,
792  jacInv,
793  x,y,z);
794  ck.calculateMappingVelocity_element(eN,
795  k,
796  mesh_velocity_dof.data(),
797  mesh_l2g.data(),
798  mesh_trial_ref.data(),
799  xt,yt,zt);
800  //get the physical integration weight
801  dV = fabs(jacDet)*dV_ref.data()[k];
802  ck.calculateG(jacInv,G,G_dd_G,tr_G);
803  //get the trial function gradients
804  ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
805  //get the solution
806  ck.valFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],u);
807  //get the solution gradients
808  ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial,grad_u);
809  //precalculate test function products with integration weights
810  for (int j=0;j<nDOF_trial_element;j++)
811  {
812  u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
813  for (int I=0;I<nSpace;I++)
814  {
815  u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
816  }
817  }
818  //
819  //calculate pde coefficients and derivatives at quadrature points
820  //
821  evaluateCoefficients(&velocity.data()[eN_k_nSpace],
822  u,
823  grad_u,
824  m,
825  dm,
826  H,
827  dH);
828  //
829  //moving mesh
830  //
831  double mesh_velocity[3];
832  mesh_velocity[0] = xt;
833  mesh_velocity[1] = yt;
834  mesh_velocity[2] = zt;
835  for (int I=0;I<nSpace;I++)
836  {
837  f[I] = -MOVING_DOMAIN*m*mesh_velocity[I];
838  df[I] = -MOVING_DOMAIN*dm*mesh_velocity[I];
839  }
840  //
841  //calculate time derivatives
842  //
843  ck.bdf(alphaBDF,
844  q_m_betaBDF.data()[eN_k],
845  m,
846  dm,
847  m_t,
848  dm_t);
849  //
850  //calculate subgrid error contribution to the Jacobian (strong residual, adjoint, jacobian of strong residual)
851  //
852  //calculate the adjoint times the test functions
853  for (int i=0;i<nDOF_test_element;i++)
854  {
855  //int eN_k_i_nSpace = (eN_k*nDOF_trial_element+i)*nSpace;
856  //Lstar_u[i]=ck.Hamiltonian_adjoint(dH,&u_grad_test_dV[eN_k_i_nSpace]);
857  register int i_nSpace = i*nSpace;
858  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]);
859 
860  }
861  //calculate the Jacobian of strong residual
862  for (int j=0;j<nDOF_trial_element;j++)
863  {
864  //int eN_k_j=eN_k*nDOF_trial_element+j;
865  //int eN_k_j_nSpace = eN_k_j*nSpace;
866  int j_nSpace = j*nSpace;
867  dpdeResidual_u_u[j]=ck.MassJacobian_strong(dm_t,u_trial_ref.data()[k*nDOF_trial_element+j]) +
868  ck.HamiltonianJacobian_strong(dH,&u_grad_trial[j_nSpace]) +
869  MOVING_DOMAIN*ck.AdvectionJacobian_strong(df,&u_grad_trial[j_nSpace]);
870 
871  }
872  //tau and tau*Res
873  double subgridErrorVelocity[nSpace];
874  for (int I=0;I<nSpace;I++)
875  subgridErrorVelocity[I] = dH[I] - MOVING_DOMAIN*df[I];
876 
877  calculateSubgridError_tau(elementDiameter.data()[eN],
878  dm_t,
879  subgridErrorVelocity,//dH,
880  cfl.data()[eN_k],
881  tau0);
882 
883  calculateSubgridError_tau(Ct_sge,
884  G,
885  dm_t,
886  subgridErrorVelocity,//dH,
887  tau1,
888  cfl.data()[eN_k]);
889  tau = useMetrics*tau1+(1.0-useMetrics)*tau0;
890 
891  for(int j=0;j<nDOF_trial_element;j++)
892  dsubgridError_u_u[j] = -tau*dpdeResidual_u_u[j];
893  for(int i=0;i<nDOF_test_element;i++)
894  {
895  //int eN_k_i=eN_k*nDOF_test_element+i;
896  //int eN_k_i_nSpace=eN_k_i*nSpace;
897  for(int j=0;j<nDOF_trial_element;j++)
898  {
899  //int eN_k_j=eN_k*nDOF_trial_element+j;
900  //int eN_k_j_nSpace = eN_k_j*nSpace;
901  int j_nSpace = j*nSpace;
902  int i_nSpace = i*nSpace;
903  elementJacobian_u_u[i][j] +=
904  ck.MassJacobian_weak(dm_t,
905  u_trial_ref.data()[k*nDOF_trial_element+j],
906  u_test_dV[i]) +
907  ck.HamiltonianJacobian_weak(dH,&u_grad_trial[j_nSpace],u_test_dV[i]) +
908  MOVING_DOMAIN*ck.AdvectionJacobian_weak(df,
909  u_trial_ref.data()[k*nDOF_trial_element+j],
910  &u_grad_test_dV[i_nSpace]) +
911  (PURE_BDF == 1 ? 0. : 1.)*
912  (ck.SubgridErrorJacobian(dsubgridError_u_u[j],Lstar_u[i]) +
913  ck.NumericalDiffusionJacobian(q_numDiff_u_last.data()[eN_k],
914  &u_grad_trial[j_nSpace],
915  &u_grad_test_dV[i_nSpace]));
916  }//j
917  }//i
918  }//k
919  //
920  //load into element Jacobian into global Jacobian
921  //
922  for (int i=0;i<nDOF_test_element;i++)
923  {
924  int eN_i = eN*nDOF_test_element+i;
925  for (int j=0;j<nDOF_trial_element;j++)
926  {
927  int eN_i_j = eN_i*nDOF_trial_element+j;
928  globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_u_u.data()[eN_i_j]] += elementJacobian_u_u[i][j];
929  }//j
930  }//i
931  }//elements
932  //
933  //loop over exterior element boundaries to compute the surface integrals and load them into the global Jacobian
934  //
935  for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
936  {
937  register int ebN = exteriorElementBoundariesArray.data()[ebNE];
938  register int eN = elementBoundaryElementsArray.data()[ebN*2+0],
939  ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
940  eN_nDOF_trial_element = eN*nDOF_trial_element;
941  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
942  {
943  register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
944  ebNE_kb_nSpace = ebNE_kb*nSpace,
945  ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
946  ebN_local_kb_nSpace = ebN_local_kb*nSpace;
947 
948  register double u_ext=0.0,
949  grad_u_ext[nSpace],
950  m_ext=0.0,
951  dm_ext=0.0,
952  H_ext=0.0,
953  dH_ext[nSpace],
954  bc_grad_u_ext[nSpace],
955  bc_H_ext=0.0,
956  bc_dH_ext[nSpace],
957  //f_ext[nSpace],
958  //df_ext[nSpace],
959  dflux_u_u_ext=0.0,
960  bc_u_ext=0.0,
961  //bc_grad_u_ext[nSpace],
962  bc_m_ext=0.0,
963  bc_dm_ext=0.0,
964  //bc_f_ext[nSpace],
965  //bc_df_ext[nSpace],
966  fluxJacobian_u_u[nDOF_trial_element],
967  jac_ext[nSpace*nSpace],
968  jacDet_ext,
969  jacInv_ext[nSpace*nSpace],
970  boundaryJac[nSpace*(nSpace-1)],
971  metricTensor[(nSpace-1)*(nSpace-1)],
972  metricTensorDetSqrt,
973  dS,
974  u_test_dS[nDOF_test_element],
975  u_grad_trial_trace[nDOF_trial_element*nSpace],
976  normal[nSpace],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
977  G[nSpace*nSpace],G_dd_G,tr_G;
978  //
979  //calculate the solution and gradients at quadrature points
980  //
981  // u_ext=0.0;
982  // for (int I=0;I<nSpace;I++)
983  // {
984  // grad_u_ext[I] = 0.0;
985  // bc_grad_u_ext[I] = 0.0;
986  // }
987  // for (int j=0;j<nDOF_trial_element;j++)
988  // {
989  // register int eN_j = eN*nDOF_trial_element+j,
990  // ebNE_kb_j = ebNE_kb*nDOF_trial_element+j,
991  // ebNE_kb_j_nSpace= ebNE_kb_j*nSpace;
992  // u_ext += valFromDOF_c(u_dof.data()[u_l2g.data()[eN_j]],u_trial_ext[ebNE_kb_j]);
993 
994  // for (int I=0;I<nSpace;I++)
995  // {
996  // grad_u_ext[I] += gradFromDOF_c(u_dof.data()[u_l2g.data()[eN_j]],u_grad_trial_ext[ebNE_kb_j_nSpace+I]);
997  // }
998  // }
999  ck.calculateMapping_elementBoundary(eN,
1000  ebN_local,
1001  kb,
1002  ebN_local_kb,
1003  mesh_dof.data(),
1004  mesh_l2g.data(),
1005  mesh_trial_trace_ref.data(),
1006  mesh_grad_trial_trace_ref.data(),
1007  boundaryJac_ref.data(),
1008  jac_ext,
1009  jacDet_ext,
1010  jacInv_ext,
1011  boundaryJac,
1012  metricTensor,
1013  metricTensorDetSqrt,
1014  normal_ref.data(),
1015  normal,
1016  x_ext,y_ext,z_ext);
1017  ck.calculateMappingVelocity_elementBoundary(eN,
1018  ebN_local,
1019  kb,
1020  ebN_local_kb,
1021  mesh_velocity_dof.data(),
1022  mesh_l2g.data(),
1023  mesh_trial_trace_ref.data(),
1024  xt_ext,yt_ext,zt_ext,
1025  normal,
1026  boundaryJac,
1027  metricTensor,
1028  integralScaling);
1029  dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref.data()[kb];
1030  //dS = metricTensorDetSqrt*dS_ref.data()[kb];
1031  ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
1032  //compute shape and solution information
1033  //shape
1034  ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,u_grad_trial_trace);
1035  //solution and gradients
1036  ck.valFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],u_ext);
1037  ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial_trace,grad_u_ext);
1038  //precalculate test function products with integration weights
1039  for (int j=0;j<nDOF_trial_element;j++)
1040  {
1041  u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
1042  }
1043  //
1044  //load the boundary values
1045  //
1046  bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+(1-isDOFBoundary_u.data()[ebNE_kb])*ebqe_rd_u_ext.data()[ebNE_kb];
1047  //
1048  //calculate the internal and external trace of the pde coefficients
1049  //
1050  evaluateCoefficients(&ebqe_velocity_ext.data()[ebNE_kb_nSpace],
1051  u_ext,
1052  grad_u_ext,
1053  m_ext,
1054  dm_ext,
1055  H_ext,
1056  dH_ext);
1057  evaluateCoefficients(&ebqe_velocity_ext.data()[ebNE_kb_nSpace],
1058  bc_u_ext,
1059  bc_grad_u_ext,
1060  bc_m_ext,
1061  bc_dm_ext,
1062  bc_H_ext,
1063  bc_dH_ext);
1064  //
1065  //moving domain
1066  //
1067  double velocity_ext[nSpace];
1068  double mesh_velocity[3];
1069  mesh_velocity[0] = xt_ext;
1070  mesh_velocity[1] = yt_ext;
1071  mesh_velocity[2] = zt_ext;
1072  for (int I=0;I<nSpace;I++)
1073  {
1074  velocity_ext[I] = - MOVING_DOMAIN*mesh_velocity[I];
1075  }
1076  //
1077  //calculate the numerical fluxes
1078  //
1079  exteriorNumericalFluxDerivative(normal,
1080  dH_ext,
1081  velocity_ext,
1082  dflux_u_u_ext);
1083  //
1084  //calculate the flux jacobian
1085  //
1086  for (int j=0;j<nDOF_trial_element;j++)
1087  {
1088  //register int ebNE_kb_j = ebNE_kb*nDOF_trial_element+j;
1089  register int ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j;
1090  fluxJacobian_u_u[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_u_u_ext,u_trial_trace_ref.data()[ebN_local_kb_j]);
1091  }//j
1092  //
1093  //update the global Jacobian from the flux Jacobian
1094  //
1095  for (int i=0;i<nDOF_test_element;i++)
1096  {
1097  register int eN_i = eN*nDOF_test_element+i;
1098  //register int ebNE_kb_i = ebNE_kb*nDOF_test_element+i;
1099  for (int j=0;j<nDOF_trial_element;j++)
1100  {
1101  register int ebN_i_j = ebN*4*nDOF_test_X_trial_element + i*nDOF_trial_element + j;
1102  globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j]] += fluxJacobian_u_u[j]*u_test_dS[i];
1103  //globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j]] += MOVING_DOMAIN*fluxJacobian_u_u[j]*u_test_dS[i];
1104  }//j
1105  }//i
1106  }//kb
1107  }//ebNE
1108  }//computeJacobian
1109 
1110  void calculateWaterline(arguments_dict& args)
1111  {
1112  xt::pyarray<int>& wlc = args.array<int>("wlc");
1113  xt::pyarray<double>& waterline = args.array<double>("waterline");
1114  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1115  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1116  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1117  xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
1118  double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
1119  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1120  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1121  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1122  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1123  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1124  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
1125  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
1126  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
1127  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
1128  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
1129  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
1130  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
1131  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
1132  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
1133  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
1134  int nElements_global = args.scalar<int>("nElements_global");
1135  double useMetrics = args.scalar<double>("useMetrics");
1136  double alphaBDF = args.scalar<double>("alphaBDF");
1137  int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
1138  double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
1139  double sc_uref = args.scalar<double>("sc_uref");
1140  double sc_alpha = args.scalar<double>("sc_alpha");
1141  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1142  xt::pyarray<int>& r_l2g = args.array<int>("r_l2g");
1143  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1144  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1145  xt::pyarray<double>& u_dof_old = args.array<double>("u_dof_old");
1146  xt::pyarray<double>& velocity = args.array<double>("velocity");
1147  xt::pyarray<double>& q_m = args.array<double>("q_m");
1148  xt::pyarray<double>& q_u = args.array<double>("q_u");
1149  xt::pyarray<double>& q_n = args.array<double>("q_n");
1150  xt::pyarray<double>& q_dH = args.array<double>("q_dH");
1151  xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
1152  xt::pyarray<double>& cfl = args.array<double>("cfl");
1153  xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
1154  xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
1155  int offset_u = args.scalar<int>("offset_u");
1156  int stride_u = args.scalar<int>("stride_u");
1157  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
1158  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
1159  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
1160  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
1161  xt::pyarray<int>& elementBoundaryMaterialTypes = args.array<int>("elementBoundaryMaterialTypes");
1162  xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
1163  xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
1164  xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
1165  xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
1166  // Tetrehedral elements specific extraction routine for waterline extraction
1167  // Loops over boundaries and checks if boundary is infact a hull (hardwired check if mattype > 6)
1168  // Extracts the nodal values of boundary triangle (4th point is dropped = hardwired assumption we are dealing with linear tet)
1169  // Then computes an average value and position for both negative and positive values
1170  // If both positive and negative values re found, and we are actually in a triangle containing the interface
1171  // a linear interpolation of negative and positive average is reported as interface location (exact in case of linear tets)
1172 
1173  for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
1174  {
1175  register int ebN = exteriorElementBoundariesArray.data()[ebNE];
1176  register int eN = elementBoundaryElementsArray.data()[ebN*2+0];
1177  register int bN = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0];
1178 
1179  if (elementBoundaryMaterialTypes.data()[ebN] >6)
1180  {
1181 
1182  double val,x,y,z;
1183  int pos=0, neg=0;
1184  double xn=0.0, yn=0.0, zn=0.0, vn=0.0;
1185  double xp=0.0, yp=0.0, zp=0.0, vp=0.0;
1186 
1187  for (int j=0;j<nDOF_trial_element;j++)
1188  {
1189  if (j != bN) {
1190  int eN_nDOF_trial_element_j = eN*nDOF_trial_element + j;
1191  int eN_nDOF_mesh_trial_element_j = eN*nDOF_mesh_trial_element + j;
1192  val = u_dof.data()[u_l2g.data()[eN_nDOF_trial_element_j]];
1193  x = mesh_dof.data()[mesh_l2g.data()[eN_nDOF_mesh_trial_element_j]*3+0];
1194  y = mesh_dof.data()[mesh_l2g.data()[eN_nDOF_mesh_trial_element_j]*3+1];
1195  z = mesh_dof.data()[mesh_l2g.data()[eN_nDOF_mesh_trial_element_j]*3+2];
1196 
1197  if (val < 0.0)
1198  {
1199  neg++;
1200  vn+=val;
1201  xn+=x;
1202  yn+=y;
1203  zn+=z;
1204  }
1205  else
1206  {
1207  pos++;
1208  vp+=val;
1209  xp+=x;
1210  yp+=y;
1211  zp+=z;
1212  }
1213  }
1214  } // trail for
1215 
1216 
1217  if ((pos > 0) && (neg > 0) )
1218  {
1219  vp /= pos;
1220  vn /= neg;
1221 
1222  double alpha = vp/(vp -vn);
1223 
1224  waterline.data()[wlc.data()[0]*3 + 0] = alpha*(xn/neg) + (1.0-alpha)*(xp/pos);
1225  waterline.data()[wlc.data()[0]*3 + 1] = alpha*(yn/neg) + (1.0-alpha)*(yp/pos);
1226  waterline.data()[wlc.data()[0]*3 + 2] = alpha*(zn/neg) + (1.0-alpha)*(zp/pos);
1227  wlc.data()[0]++;
1228 
1229  } // end value if
1230 
1231  } // end bnd mat check
1232 
1233  }//ebNE
1234 
1235  //std::cout<<"CPP WLC "<<wlc.data()[0]<<std::endl;
1236  } // calcWaterline
1237 
1238  double calculateRedistancingResidual(arguments_dict& args)
1239  {
1240  double dt = args.scalar<double>("dt");
1241  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1242  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1243  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1244  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1245  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1246  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1247  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1248  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1249  int nElements_global = args.scalar<int>("nElements_global");
1250  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1251  xt::pyarray<int>& r_l2g = args.array<int>("r_l2g");
1252  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1253  xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
1254  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1255  xt::pyarray<double>& u_dof_old = args.array<double>("u_dof_old");
1256  xt::pyarray<double>& uStar_dof = args.array<double>("uStar_dof");
1257  int offset_u = args.scalar<int>("offset_u");
1258  int stride_u = args.scalar<int>("stride_u");
1259  xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
1260  int numDOFs = args.scalar<int>("numDOFs");
1261  int NNZ = args.scalar<int>("NNZ");
1262  xt::pyarray<int>& csrRowIndeces_DofLoops = args.array<int>("csrRowIndeces_DofLoops");
1263  xt::pyarray<int>& csrColumnOffsets_DofLoops = args.array<int>("csrColumnOffsets_DofLoops");
1264  xt::pyarray<int>& csrRowIndeces_CellLoops = args.array<int>("csrRowIndeces_CellLoops");
1265  xt::pyarray<int>& csrColumnOffsets_CellLoops = args.array<int>("csrColumnOffsets_CellLoops");
1266  double lambda_coupez = args.scalar<double>("lambda_coupez");
1267  double epsCoupez = args.scalar<double>("epsCoupez");
1268  double epsFactRedistancing = args.scalar<double>("epsFactRedistancing");
1269  xt::pyarray<double>& edge_based_cfl = args.array<double>("edge_based_cfl");
1270  int SATURATED_LEVEL_SET = args.scalar<int>("SATURATED_LEVEL_SET");
1271  xt::pyarray<double>& Cx = args.array<double>("Cx");
1272  xt::pyarray<double>& Cy = args.array<double>("Cy");
1273  xt::pyarray<double>& Cz = args.array<double>("Cz");
1274  xt::pyarray<double>& ML = args.array<double>("ML");
1275  L2_norm_per_node.resize(numDOFs,0.0);
1276  double L2_norm=0.;
1277  for (int i=0; i<numDOFs; i++)
1278  L2_norm_per_node[i] = 0.;
1279 
1280  // Allocate space for the transport matrices
1281  TransportMatrix.resize(NNZ,0.0);
1282  TransposeTransportMatrix.resize(NNZ,0.0);
1283 
1285  // ** LOOP IN CELLS FOR CELL BASED TERMS ** //
1287  // HERE WE COMPUTE:
1288  // * Time derivative term
1289  // * Transport matrices
1290  for(int eN=0;eN<nElements_global;eN++)
1291  {
1292  //declare local storage for local contributions and initialize
1293  register double
1294  elementResidual_u[nDOF_test_element], element_L2_norm_per_node[nDOF_test_element];
1295  register double elementTransport[nDOF_test_element][nDOF_trial_element];
1296  register double elementTransposeTransport[nDOF_test_element][nDOF_trial_element];
1297  for (int i=0;i<nDOF_test_element;i++)
1298  {
1299  elementResidual_u[i]=0.0;
1300  element_L2_norm_per_node[i]=0.;
1301  for (int j=0;j<nDOF_trial_element;j++)
1302  {
1303  elementTransport[i][j]=0.0;
1304  elementTransposeTransport[i][j]=0.0;
1305  }
1306  }
1307  //loop over quadrature points and compute integrands
1308  for (int k=0;k<nQuadraturePoints_element;k++)
1309  {
1310  //compute indeces and declare local storage
1311  register int eN_k = eN*nQuadraturePoints_element+k,
1312  eN_k_nSpace = eN_k*nSpace,
1313  eN_nDOF_trial_element = eN*nDOF_trial_element;
1314  register double
1315  hquad=0.0, u=0.0, u_test_dV[nDOF_trial_element], uStar=0.0,
1316  //for entropy viscosity
1317  un=0.0, grad_u[nSpace], grad_un[nSpace], vn[nSpace],
1318  u_grad_trial[nDOF_trial_element*nSpace],
1319  //for general use
1320  jac[nSpace*nSpace], jacDet, jacInv[nSpace*nSpace],
1321  dV,x,y,z;
1322  //get the physical integration weight
1323  ck.calculateMapping_element(eN,k,mesh_dof.data(),mesh_l2g.data(),mesh_trial_ref.data(),mesh_grad_trial_ref.data(),jac,jacDet,jacInv,x,y,z);
1324  dV = fabs(jacDet)*dV_ref.data()[k];
1325  // get h at quad points
1326  ck.valFromDOF(nodeDiametersArray.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],hquad);
1327  // get the star solution at quad points
1328  ck.valFromDOF(uStar_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],uStar);
1329  //get the solution at quad points
1330  ck.valFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],u);
1331  //get the solution at quad point at tn
1332  ck.valFromDOF(u_dof_old.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],un);
1333  // solution and grads at old times at quad points
1334  ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
1335  ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial,grad_u);
1336  ck.gradFromDOF(u_dof_old.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial,grad_un);
1337  //precalculate test function products with integration weights for mass matrix terms
1338  for (int j=0;j<nDOF_trial_element;j++)
1339  u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
1340 
1341  double norm_grad_un = std::sqrt(std::pow(grad_un[0],2) + std::pow(grad_un[1],2)) + 1E-10;
1342  double norm_grad_u = std::sqrt(std::pow(grad_u[0],2) + std::pow(grad_u[1],2));
1343  double dist_error = norm_grad_u - (1-SATURATED_LEVEL_SET*std::pow(u/epsCoupez,2));
1344 
1345  double sgn = sign(uStar,epsFactRedistancing);
1346  vn[0] = lambda_coupez*sgn*grad_un[0]/norm_grad_un;
1347  vn[1] = lambda_coupez*sgn*grad_un[1]/norm_grad_un;
1348 
1350  // ith-LOOP //
1352  double residual = (u-un)
1353  - dt*lambda_coupez*sgn*(1-SATURATED_LEVEL_SET*std::pow(un/epsCoupez,2));
1354  for(int i=0;i<nDOF_test_element;i++)
1355  {
1356  // lumped mass matrix
1357  elementResidual_u[i] += residual*u_test_dV[i];
1358  element_L2_norm_per_node[i] += dist_error*u_test_dV[i];
1360  // j-th LOOP // To construct transport matrices
1362  for(int j=0;j<nDOF_trial_element;j++)
1363  {
1364  int j_nSpace = j*nSpace;
1365  int i_nSpace = i*nSpace;
1366  // COMPUTE ELEMENT TRANSPORT MATRIX (MQL)
1367  elementTransport[i][j] += // int[(vel.grad_wj)*wi*dx]
1368  ck.HamiltonianJacobian_weak(vn,&u_grad_trial[j_nSpace],u_test_dV[i]);
1369  elementTransposeTransport[i][j] += // int[(vel.grad_wi)*wj*dx]
1370  ck.HamiltonianJacobian_weak(vn,&u_grad_trial[i_nSpace],u_test_dV[j]);
1371  }
1372  }//i
1373  }
1374 
1376  // DISTRIBUTE // load cell based element into global residual
1378  for(int i=0;i<nDOF_test_element;i++)
1379  {
1380  int eN_i=eN*nDOF_test_element+i;
1381  int gi = offset_u+stride_u*r_l2g.data()[eN_i]; //global i-th index
1382 
1383  // distribute global residual for (lumped) mass matrix
1384  globalResidual.data()[gi] += elementResidual_u[i];
1385  // distribute L2 norm per node
1386  L2_norm_per_node[gi] += element_L2_norm_per_node[i];
1387  // distribute transport matrices
1388  for (int j=0;j<nDOF_trial_element;j++)
1389  {
1390  int eN_i_j = eN_i*nDOF_trial_element+j;
1391  TransportMatrix[csrRowIndeces_CellLoops.data()[eN_i] + csrColumnOffsets_CellLoops.data()[eN_i_j]]
1392  += elementTransport[i][j];
1393  TransposeTransportMatrix[csrRowIndeces_CellLoops.data()[eN_i] + csrColumnOffsets_CellLoops.data()[eN_i_j]]
1394  += elementTransposeTransport[i][j];
1395  }//j
1396  }//i
1397  }//elements
1398 
1400  // COMPUTE SMOOTHNESS INDICATOR //
1402  // Smoothness indicator is based on the solution.
1403  // psi_i = psi_i(alpha_i); alpha_i = |sum(betaij*(uj-ui))|/sum(betaij*|uj-ui|)
1404  psi.resize(numDOFs,0.0);
1405  for (int i=0; i<numDOFs; i++)
1406  {
1407  double alphai;
1408  double solni = u_dof_old.data()[i]; // solution at time tn for the ith DOF
1409  // for smoothness indicator
1410  double alpha_numi=0;
1411  double alpha_deni=0;
1412  for (int offset=csrRowIndeces_DofLoops.data()[i]; offset<csrRowIndeces_DofLoops.data()[i+1]; offset++)
1413  { //loop in j (sparsity pattern)
1414  int j = csrColumnOffsets_DofLoops.data()[offset];
1415  double solnj = u_dof_old.data()[j]; // solution at time tn for the jth DOF
1416  alpha_numi += solnj-solni;
1417  alpha_deni += fabs(solnj-solni);
1418  }
1419  alphai = (fabs(alpha_numi) == alpha_deni ? 1. : fabs(alpha_numi)/(alpha_deni+1E-15));
1420  if (POWER_SMOOTHNESS_INDICATOR==0)
1421  psi[i] = 1.0;
1422  else
1423  psi[i] = std::pow(alphai,POWER_SMOOTHNESS_INDICATOR); //NOTE: they use alpha^2 in the paper
1424  // compute L2 norm
1425  L2_norm += std::pow(L2_norm_per_node[i],2);
1426  }
1427  L2_norm = std::sqrt(L2_norm);
1428  // END OF COMPUTING SMOOTHNESS INDICATOR
1429 
1431  // ** LOOP IN DOFs FOR EDGE BASED TERMS ** //
1433  int ij=0;
1434  for (int i=0; i<numDOFs; i++)
1435  {
1436  double solni = u_dof_old.data()[i]; // solution at time tn for the ith DOF
1437  double ith_dissipative_term = 0;
1438  double ith_flux_term = 0;
1439  double dLii = 0.;
1440 
1441  // loop over the sparsity pattern of the i-th DOF
1442  for (int offset=csrRowIndeces_DofLoops.data()[i]; offset<csrRowIndeces_DofLoops.data()[i+1]; offset++)
1443  {
1444  int j = csrColumnOffsets_DofLoops.data()[offset];
1445  double solnj = u_dof_old.data()[j]; // solution at time tn for the jth DOF
1446  double dLij=0.;
1447 
1448  ith_flux_term += TransportMatrix[ij]*solnj;
1449  if (i != j) //NOTE: there is really no need to check for i!=j (see formula for ith_dissipative_term)
1450  {
1451  // first-order dissipative operator
1452  dLij = fmax(fabs(TransportMatrix[ij]),fabs(TransposeTransportMatrix[ij]));
1453  dLij *= fmax(psi[i],psi[j]);
1454  //dissipative terms
1455  ith_dissipative_term += dLij*(solnj-solni);
1456  dLii -= dLij;
1457  }
1458  //update ij
1459  ij+=1;
1460  }
1461  // update residual
1462  double mi = ML.data()[i];
1463  // compute edge_based_cfl
1464  edge_based_cfl.data()[i] = 2*fabs(dLii)/mi;
1465  globalResidual.data()[i] += dt*(ith_flux_term - ith_dissipative_term);
1466  }
1467  return L2_norm;
1468  }
1469 
1470  double calculateRhsSmoothing(arguments_dict& args)
1471  {
1472  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1473  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1474  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1475  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1476  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1477  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1478  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1479  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1480  int nElements_global = args.scalar<int>("nElements_global");
1481  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1482  xt::pyarray<int>& r_l2g = args.array<int>("r_l2g");
1483  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1484  xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
1485  xt::pyarray<double>& u_dof_old = args.array<double>("u_dof_old");
1486  int offset_u = args.scalar<int>("offset_u");
1487  int stride_u = args.scalar<int>("stride_u");
1488  xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
1490  // ** LOOP IN CELLS FOR CELL BASED TERMS ** //
1492  for(int eN=0;eN<nElements_global;eN++)
1493  {
1494  //declare local storage for local contributions and initialize
1495  register double
1496  elementResidual_u[nDOF_test_element];
1497  for (int i=0;i<nDOF_test_element;i++)
1498  elementResidual_u[i]=0.0;
1499 
1500  //loop over quadrature points and compute integrands
1501  for (int k=0;k<nQuadraturePoints_element;k++)
1502  {
1503  //compute indeces and declare local storage
1504  register int eN_k = eN*nQuadraturePoints_element+k,
1505  eN_k_nSpace = eN_k*nSpace,
1506  eN_nDOF_trial_element = eN*nDOF_trial_element;
1507  register double
1508  un=0.0, u_test_dV[nDOF_trial_element], u_grad_trial[nDOF_trial_element*nSpace],
1509  //for general use
1510  jac[nSpace*nSpace], jacDet, jacInv[nSpace*nSpace],
1511  dV,x,y,z;
1512  //get the physical integration weight
1513  ck.calculateMapping_element(eN,k,mesh_dof.data(),mesh_l2g.data(),mesh_trial_ref.data(),mesh_grad_trial_ref.data(),jac,jacDet,jacInv,x,y,z);
1514  dV = fabs(jacDet)*dV_ref.data()[k];
1515  //get the solution at quad point at tn
1516  ck.valFromDOF(u_dof_old.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],un);
1517  //precalculate test function products with integration weights for mass matrix terms
1518  for (int j=0;j<nDOF_trial_element;j++)
1519  u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
1520 
1522  // ith-LOOP //
1524  for(int i=0;i<nDOF_test_element;i++)
1525  elementResidual_u[i] += un*u_test_dV[i];
1526  }
1527 
1529  // DISTRIBUTE // load cell based element into global residual
1531  for(int i=0;i<nDOF_test_element;i++)
1532  {
1533  int eN_i=eN*nDOF_test_element+i;
1534  int gi = offset_u+stride_u*r_l2g.data()[eN_i]; //global i-th index
1535 
1536  // distribute global residual
1537  globalResidual.data()[gi] += elementResidual_u[i];
1538  }//i
1539  }//elements
1540  }
1541 
1542  void calculateResidual_entropy_viscosity(arguments_dict& args)
1543  {
1544  double dt = args.scalar<double>("dt");
1545  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1546  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1547  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1548  xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
1549  double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
1550  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1551  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1552  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1553  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1554  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1555  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
1556  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
1557  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
1558  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
1559  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
1560  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
1561  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
1562  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
1563  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
1564  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
1565  int nElements_global = args.scalar<int>("nElements_global");
1566  double useMetrics = args.scalar<double>("useMetrics");
1567  double alphaBDF = args.scalar<double>("alphaBDF");
1568  int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
1569  double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
1570  double sc_uref = args.scalar<double>("sc_uref");
1571  double sc_alpha = args.scalar<double>("sc_alpha");
1572  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1573  xt::pyarray<int>& r_l2g = args.array<int>("r_l2g");
1574  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1575  xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
1576  int degree_polynomial = args.scalar<int>("degree_polynomial");
1577  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1578  xt::pyarray<double>& u_dof_old = args.array<double>("u_dof_old");
1579  xt::pyarray<double>& uStar_dof = args.array<double>("uStar_dof");
1580  xt::pyarray<double>& velocity = args.array<double>("velocity");
1581  xt::pyarray<double>& q_m = args.array<double>("q_m");
1582  xt::pyarray<double>& q_u = args.array<double>("q_u");
1583  xt::pyarray<double>& q_n = args.array<double>("q_n");
1584  xt::pyarray<double>& q_dH = args.array<double>("q_dH");
1585  xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
1586  xt::pyarray<double>& q_dV = args.array<double>("q_dV");
1587  xt::pyarray<double>& q_dV_last = args.array<double>("q_dV_last");
1588  xt::pyarray<double>& cfl = args.array<double>("cfl");
1589  xt::pyarray<double>& edge_based_cfl = args.array<double>("edge_based_cfl");
1590  xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
1591  xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
1592  int offset_u = args.scalar<int>("offset_u");
1593  int stride_u = args.scalar<int>("stride_u");
1594  xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
1595  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
1596  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
1597  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
1598  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
1599  xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
1600  xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
1601  xt::pyarray<double>& ebqe_rd_u_ext = args.array<double>("ebqe_rd_u_ext");
1602  xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
1603  xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
1604  xt::pyarray<double>& ebqe_grad_u = args.array<double>("ebqe_grad_u");
1605  xt::pyarray<double>& interface_locator = args.array<double>("interface_locator");
1606  int PURE_BDF = args.scalar<int>("PURE_BDF");
1607  int numDOFs = args.scalar<int>("numDOFs");
1608  int NNZ = args.scalar<int>("NNZ");
1609  xt::pyarray<int>& csrRowIndeces_DofLoops = args.array<int>("csrRowIndeces_DofLoops");
1610  xt::pyarray<int>& csrColumnOffsets_DofLoops = args.array<int>("csrColumnOffsets_DofLoops");
1611  xt::pyarray<int>& csrRowIndeces_CellLoops = args.array<int>("csrRowIndeces_CellLoops");
1612  xt::pyarray<int>& csrColumnOffsets_CellLoops = args.array<int>("csrColumnOffsets_CellLoops");
1613  xt::pyarray<int>& csrColumnOffsets_eb_CellLoops = args.array<int>("csrColumnOffsets_eb_CellLoops");
1614  int LUMPED_MASS_MATRIX = args.scalar<int>("LUMPED_MASS_MATRIX");
1615  xt::pyarray<double>& quantDOFs = args.array<double>("quantDOFs");
1616  double lambda_coupez = args.scalar<double>("lambda_coupez");
1617  double epsCoupez = args.scalar<double>("epsCoupez");
1618  double epsFactRedistancing = args.scalar<double>("epsFactRedistancing");
1619  int COUPEZ = args.scalar<int>("COUPEZ");
1620  int SATURATED_LEVEL_SET = args.scalar<int>("SATURATED_LEVEL_SET");
1621  xt::pyarray<double>& Cx = args.array<double>("Cx");
1622  xt::pyarray<double>& Cy = args.array<double>("Cy");
1623  xt::pyarray<double>& Cz = args.array<double>("Cz");
1624  xt::pyarray<double>& ML = args.array<double>("ML");
1625  int STABILIZATION_TYPE = args.scalar<int>("STABILIZATION_TYPE");
1626  int ENTROPY_TYPE = args.scalar<int>("ENTROPY_TYPE");
1627  double cE = args.scalar<double>("cE");
1628  // NOTE: This function follows a different (but equivalent) implementation of the smoothness based indicator than VOF.h
1629  // inverse of mass matrix in reference element
1630  //double inverse_elMassMatrix_degree2 [9][9]=
1631  //{
1632  //{81/4., 27/4., 9/4., 27/4., -81/4., -27/4., -27/4., -81/4., 81/4.},
1633  //{27/4., 81/4., 27/4., 9/4., -81/4., -81/4., -27/4., -27/4., 81/4.},
1634  //{9/4., 27/4., 81/4., 27/4., -27/4., -81/4., -81/4., -27/4., 81/4.},
1635  //{27/4., 9/4., 27/4., 81/4., -27/4., -27/4., -81/4., -81/4., 81/4.},
1636  //{-81/4., -81/4., -27/4., -27/4., 189/4., 81/4., 63/4., 81/4., -189/4.},
1637  //{-27/4., -81/4., -81/4., -27/4., 81/4., 189/4., 81/4., 63/4., -189/4.},
1638  //{-27/4., -27/4., -81/4., -81/4., 63/4., 81/4., 189/4., 81/4., -189/4.},
1639  //{-81/4., -27/4., -27/4., -81/4., 81/4., 63/4., 81/4., 189/4., -189/4.},
1640  //{81/4., 81/4., 81/4., 81/4., -189/4., -189/4., -189/4., -189/4., 441/4.}
1641  //};
1642  //double inverse_elMassMatrix_degree1 [4][4] = {{4, -2, 1, -2}, {-2, 4, -2, 1}, {1, -2, 4, -2}, {-2, 1, -2, 4}};
1643  double JacDet = 0;
1644 
1645  // Allocate space for the transport matrices
1646  // This is used for first order KUZMIN'S METHOD
1647  TransportMatrix.resize(NNZ,0.0);
1648  TransposeTransportMatrix.resize(NNZ,0.0);
1649  // Allocate and init to zero the Entropy residual vector
1650  global_entropy_residual.resize(numDOFs,0.0);
1651  if (STABILIZATION_TYPE==1) //EV stab
1652  for (int i=0; i<numDOFs; i++)
1653  global_entropy_residual[i] = 0.;
1654 
1656  // ** LOOP IN CELLS FOR CELL BASED TERMS ** //
1658  // HERE WE COMPUTE:
1659  // * Time derivative term
1660  // * Cell based CFL (for reference)
1661  // * Entropy residual
1662  // * Transport matrices
1663  for(int eN=0;eN<nElements_global;eN++)
1664  {
1665  //declare local storage for local contributions and initialize
1666  register double
1667  elementResidual_u[nDOF_test_element],
1668  element_entropy_residual[nDOF_test_element];
1669  register double elementTransport[nDOF_test_element][nDOF_trial_element];
1670  register double elementTransposeTransport[nDOF_test_element][nDOF_trial_element];
1671  //register double preconditioned_elementTransport[nDOF_test_element][nDOF_trial_element];
1672  for (int i=0;i<nDOF_test_element;i++)
1673  {
1674  elementResidual_u[i]=0.0;
1675  element_entropy_residual[i]=0.0;
1676  for (int j=0;j<nDOF_trial_element;j++)
1677  {
1678  elementTransport[i][j]=0.0;
1679  elementTransposeTransport[i][j]=0.0;
1680  // preconditioned elementTransport
1681  //preconditioned_elementTransport[i][j]=0.0;
1682  }
1683  }
1684  //loop over quadrature points and compute integrands
1685  for (int k=0;k<nQuadraturePoints_element;k++)
1686  {
1687  //compute indeces and declare local storage
1688  register int eN_k = eN*nQuadraturePoints_element+k,
1689  eN_k_nSpace = eN_k*nSpace,
1690  eN_nDOF_trial_element = eN*nDOF_trial_element;
1691  register double
1692  // for entropy residual
1693  aux_entropy_residual=0., DENTROPY_un, DENTROPY_uni,
1694  //for mass matrix contributions
1695  u=0.0, uStar=0.0, hquad=0.0,
1696  u_test_dV[nDOF_trial_element],
1697  //for entropy viscosity
1698  un=0.0, unm1=0.0, grad_u[nSpace], grad_un[nSpace], vn[nSpace],
1699  u_grad_trial[nDOF_trial_element*nSpace],
1700  //for general use
1701  jac[nSpace*nSpace], jacDet, jacInv[nSpace*nSpace],
1702  dV,x,y,z,xt,yt,zt;
1703  //get the physical integration weight
1704  ck.calculateMapping_element(eN,
1705  k,
1706  mesh_dof.data(),
1707  mesh_l2g.data(),
1708  mesh_trial_ref.data(),
1709  mesh_grad_trial_ref.data(),
1710  jac,
1711  jacDet,
1712  jacInv,x,y,z);
1713  ck.calculateMappingVelocity_element(eN,
1714  k,
1715  mesh_velocity_dof.data(),
1716  mesh_l2g.data(),
1717  mesh_trial_ref.data(),
1718  xt,yt,zt);
1719  dV = fabs(jacDet)*dV_ref.data()[k];
1720  JacDet = fabs(jacDet);
1721  // get h at quad points
1722  ck.valFromDOF(nodeDiametersArray.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],hquad);
1723  //get the solution (of Newton's solver). To compute time derivative term
1724  ck.valFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],u);
1725  ck.valFromDOF(uStar_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],uStar);
1726  //get the solution at quad point at tn and tnm1 for entropy viscosity
1727  // solution and grads at old times at quad points
1728  ck.valFromDOF(u_dof_old.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],un);
1729  //get the solution gradients at tn for entropy viscosity
1730  ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
1731  ck.gradFromDOF(u_dof_old.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial,grad_un);
1732  //precalculate test function products with integration weights for mass matrix terms
1733  for (int j=0;j<nDOF_trial_element;j++)
1734  u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
1735 
1736  //
1737  //moving mesh
1738  //
1739  double mesh_velocity[3];
1740  mesh_velocity[0] = xt;
1741  mesh_velocity[1] = yt;
1742  mesh_velocity[2] = zt;
1743  //epsCoupez = 3*hquad;
1744  double norm_grad_un = 0.;
1745  for (int I=0;I<nSpace;I++)
1746  norm_grad_un += std::pow(grad_un[I],2);
1747  norm_grad_un = std::sqrt(norm_grad_un) + 1E-10;
1748  double dist_error = fabs(norm_grad_un - (1-SATURATED_LEVEL_SET*std::pow(un/epsCoupez,2)));
1749  double sgn = sign(uStar,epsFactRedistancing);
1750 
1751  // get velocity
1752  for (int I=0;I<nSpace;I++)
1753  vn[I] = (velocity.data()[eN_k_nSpace+I] - MOVING_DOMAIN*mesh_velocity[I]
1754  + dist_error*COUPEZ*lambda_coupez*sgn*grad_un[I]/norm_grad_un);
1755 
1757  // CALCULATE CELL BASED CFL //
1759  calculateCFL(elementDiameter.data()[eN]/degree_polynomial,vn,cfl.data()[eN_k]);
1760 
1762  // CALCULATE ENTROPY RESIDUAL AT QUAD POINT //
1764  if (STABILIZATION_TYPE==1) // EV stab
1765  {
1766  //double entropy_residual = ((un - unm1)/dt + vn[0]*grad_un[0] + vn[1]*grad_un[1]
1767  // -dist_error*COUPEZ*lambda_coupez*sgn
1768  // *(1-SATURATED_LEVEL_SET*std::pow(un/epsCoupez,2)))*DENTROPY(un,-epsCoupez,epsCoupez);
1769  for (int I=0;I<nSpace;I++)
1770  aux_entropy_residual += vn[I]*grad_un[I];
1771  aux_entropy_residual -= dist_error*COUPEZ*lambda_coupez*sgn*(1-SATURATED_LEVEL_SET*std::pow(un/epsCoupez,2));
1772  DENTROPY_un = ENTROPY_TYPE==1 ? DENTROPY(un,-epsCoupez,epsCoupez) : DENTROPY_LOG(un,-epsCoupez,epsCoupez);
1773  }
1775  // ith-LOOP //
1777  double residual = (u-un) - dt*dist_error*COUPEZ*lambda_coupez*sgn*(1-SATURATED_LEVEL_SET*std::pow(un/epsCoupez,2));
1778 
1779  for(int i=0;i<nDOF_test_element;i++)
1780  {
1781  int eN_i=eN*nDOF_test_element+i;
1782  int gi = offset_u+stride_u*u_l2g.data()[eN_i]; //global i-th index
1783  // entropy residual
1784  if (STABILIZATION_TYPE==1) // EV Stab
1785  {
1786  DENTROPY_uni = ENTROPY_TYPE == 1 ? DENTROPY(u_dof_old.data()[gi],-epsCoupez,epsCoupez) : DENTROPY_LOG(u_dof_old.data()[gi],-epsCoupez,epsCoupez);
1787  element_entropy_residual[i] += (DENTROPY_un - DENTROPY_uni)*aux_entropy_residual*u_test_dV[i];
1788  }
1789  // element residual
1790  elementResidual_u[i] += residual*u_test_dV[i];
1792  // j-th LOOP // To construct transport matrices
1794  for(int j=0;j<nDOF_trial_element;j++)
1795  {
1796  int j_nSpace = j*nSpace;
1797  int i_nSpace = i*nSpace;
1798  // COMPUTE ELEMENT TRANSPORT MATRIX (MQL)
1799  elementTransport[i][j] += // int[(vel.grad_wj)*wi*dx]
1800  ck.HamiltonianJacobian_weak(vn,&u_grad_trial[j_nSpace],u_test_dV[i]);
1801  elementTransposeTransport[i][j] += // int[(vel.grad_wi)*wj*dx]
1802  ck.HamiltonianJacobian_weak(vn,&u_grad_trial[i_nSpace],u_test_dV[j]);
1803  }
1804  }//i
1805  //save solution at quadrature points for other models to use
1806  q_u.data()[eN_k]=u;
1807  q_m.data()[eN_k] = u;
1808  for (int I=0;I<nSpace;I++)
1809  q_n.data()[eN_k_nSpace+I] = grad_u[I];
1810  }
1811  // MULTIPLY Transport Matrix times preconditioner//
1812  //if (false)
1813  // {
1814  // for (int i=0; i < nDOF_test_element; i++)
1815  //for (int j=0; j < nDOF_test_element; j++)
1816  // for (int r=0; r < nDOF_test_element; r++)
1817  // if (degree_polynomial == 1)
1818  // preconditioned_elementTransport[i][j] +=
1819  // element_lumped_mass_matrix[i]*1./JacDet*inverse_elMassMatrix_degree1[i][r]*elementTransport[r][j];
1820  // else
1821  // preconditioned_elementTransport[i][j] +=
1822  // element_lumped_mass_matrix[i]*1./JacDet*inverse_elMassMatrix_degree2[i][r]*elementTransport[r][j];
1823  //for (int i=0; i < nDOF_test_element; i++)
1824  //for (int j=0; j < nDOF_test_element; j++)
1825  //{
1826  // elementTransport[i][j] = preconditioned_elementTransport[i][j];
1827  // elementTransposeTransport[i][j] = preconditioned_elementTransport[j][i];
1828  //}
1829  //}
1830 
1832  // DISTRIBUTE // load cell based element into global residual
1834  for(int i=0;i<nDOF_test_element;i++)
1835  {
1836  int eN_i=eN*nDOF_test_element+i;
1837  int gi = offset_u+stride_u*r_l2g.data()[eN_i]; //global i-th index
1838 
1839  // distribute global residual for (lumped) mass matrix
1840  globalResidual.data()[gi] += elementResidual_u[i];
1841  // distribute entropy_residual
1842  if (STABILIZATION_TYPE==1) // EV Stab
1843  global_entropy_residual[gi] += element_entropy_residual[i];
1844 
1845  // distribute transport matrices
1846  for (int j=0;j<nDOF_trial_element;j++)
1847  {
1848  int eN_i_j = eN_i*nDOF_trial_element+j;
1849  TransportMatrix[csrRowIndeces_CellLoops.data()[eN_i] + csrColumnOffsets_CellLoops.data()[eN_i_j]]
1850  += elementTransport[i][j];
1851  TransposeTransportMatrix[csrRowIndeces_CellLoops.data()[eN_i] + csrColumnOffsets_CellLoops.data()[eN_i_j]]
1852  += elementTransposeTransport[i][j];
1853  }//j
1854  }//i
1855  }//elements
1856 
1858  // COMPUTE g, ETA and others //
1860  // NOTE: see VOF.h for a different but equivalent implementation of this.
1861  gx.resize(numDOFs,0.0);
1862  gy.resize(numDOFs,0.0);
1863  gz.resize(numDOFs,0.0);
1864  eta.resize(numDOFs,0.0);
1865  alpha_numerator_pos.resize(numDOFs,0.0);
1866  alpha_numerator_neg.resize(numDOFs,0.0);
1867  alpha_denominator_pos.resize(numDOFs,0.0);
1868  alpha_denominator_neg.resize(numDOFs,0.0);
1869  int ij = 0;
1870  for (int i=0; i<numDOFs; i++)
1871  {
1872  double solni = u_dof_old.data()[i];
1873  if (STABILIZATION_TYPE==1) //EV Stab
1874  eta[i] = ENTROPY_TYPE == 1 ? ENTROPY(solni,-epsCoupez,epsCoupez) : ENTROPY_LOG(solni,-epsCoupez,epsCoupez);
1875  else // smoothness based stab
1876  {
1877  gx[i]=0.;
1878  gy[i]=0.;
1879  gz[i]=0.;
1880  // for smoothness indicator //
1881  alpha_numerator_pos[i] = 0.;
1882  alpha_numerator_neg[i] = 0.;
1883  alpha_denominator_pos[i] = 0.;
1884  alpha_denominator_neg[i] = 0.;
1885 
1886  for (int offset=csrRowIndeces_DofLoops.data()[i]; offset<csrRowIndeces_DofLoops.data()[i+1]; offset++)
1887  {
1888  int j = csrColumnOffsets_DofLoops.data()[offset];
1889  double solnj = u_dof_old.data()[j];
1890 
1891  // for gi vector
1892  gx[i] += Cx.data()[ij]*solnj;
1893  gy[i] += Cy.data()[ij]*solnj;
1894 #if nSpace==3
1895  gz[i] += Cz.data()[ij]*solnj;
1896 #endif
1897  double alpha_num = solni - solnj;
1898  if (alpha_num >= 0.)
1899  {
1900  alpha_numerator_pos[i] += alpha_num;
1901  alpha_denominator_pos[i] += alpha_num;
1902  }
1903  else
1904  {
1905  alpha_numerator_neg[i] += alpha_num;
1906  alpha_denominator_neg[i] += fabs(alpha_num);
1907  }
1908  //update ij
1909  ij+=1;
1910  }
1911  gx[i] /= ML.data()[i];
1912  gy[i] /= ML.data()[i];
1913 #if nSpace==3
1914  gz[i] /= ML.data()[i];
1915 #endif
1916  }
1917  }
1918 
1920  // COMPUTE SMOOTHNESS INDICATOR and ENTROPY MIN and MAX //
1922  // Smoothness indicator is based on the solution.
1923  // psi_i = psi_i(alpha_i); alpha_i = |sum(betaij*(uj-ui))|/sum(betaij*|uj-ui|)
1924  psi.resize(numDOFs,0.0);
1925  etaMax.resize(numDOFs,0.0);
1926  etaMin.resize(numDOFs,0.0);
1927  for (int i=0; i<numDOFs; i++)
1928  {
1929  double xi, yi, zi, SumPos=0., SumNeg=0.;
1930  if (STABILIZATION_TYPE==1) //EV Stabilization
1931  {
1932  // For eta min and max
1933  etaMax[i] = fabs(eta[i]);
1934  etaMin[i] = fabs(eta[i]);
1935  }
1936  else //smoothness based stab
1937  {
1938  xi = mesh_dof.data()[i*3+0];
1939  yi = mesh_dof.data()[i*3+1];
1940 #if nSpace==3
1941  zi = mesh_dof.data()[i*3+2];
1942 #endif
1943  }
1944  for (int offset=csrRowIndeces_DofLoops.data()[i]; offset<csrRowIndeces_DofLoops.data()[i+1]; offset++)
1945  { //loop in j (sparsity pattern)
1946  int j = csrColumnOffsets_DofLoops.data()[offset];
1947  if (STABILIZATION_TYPE==1) //EV stab
1948  {
1950  // COMPUTE ETA MIN AND ETA MAX //
1952  etaMax[i] = fmax(etaMax[i],fabs(eta[j]));
1953  etaMin[i] = fmin(etaMin[i],fabs(eta[j]));
1954  }
1955  else //smoothness based sta
1956  {
1957  double xj = mesh_dof.data()[j*3+0];
1958  double yj = mesh_dof.data()[j*3+1];
1959  double zj = 0;
1960 #if nSpace==3
1961  zj = mesh_dof.data()[j*3+2];
1962 #endif
1963  // FOR SMOOTHNESS INDICATOR //
1966  double gi_times_x = gx[i]*(xi-xj) + gy[i]*(yi-yj);
1967 #if nSpace==3
1968  gi_times_x += gz[i]*(zi-zj);
1969 #endif
1970  SumPos += gi_times_x > 0 ? gi_times_x : 0;
1971  SumNeg += gi_times_x < 0 ? gi_times_x : 0;
1972  }
1973  }
1974  if (STABILIZATION_TYPE==1) //EV stab
1975  {
1976  // Normalize entropy residual
1977  global_entropy_residual[i] *= etaMin[i] == etaMax[i] ? 0. : 2*cE/(etaMax[i]-etaMin[i]);
1978  quantDOFs.data()[i] = fabs(global_entropy_residual[i]);
1979  }
1980  else // smoothness based stab
1981  {
1982  // Compute sigmaPos and sigmaNeg
1983  double sigmaPosi = fmin(1.,(fabs(SumNeg)+1E-15)/(SumPos+1E-15));
1984  double sigmaNegi = fmin(1.,(SumPos+1E-15)/(fabs(SumNeg)+1E-15));
1985  double alpha_numi = fabs(sigmaPosi*alpha_numerator_pos[i] + sigmaNegi*alpha_numerator_neg[i]);
1986  double alpha_deni = sigmaPosi*alpha_denominator_pos[i] + sigmaNegi*alpha_denominator_neg[i];
1987  if (IS_BETAij_ONE == 1)
1988  {
1989  alpha_numi = fabs(alpha_numerator_pos[i] + alpha_numerator_neg[i]);
1990  alpha_deni = alpha_denominator_pos[i] + alpha_denominator_neg[i];
1991  }
1992  double alphai = alpha_numi/(alpha_deni+1E-15);
1993  quantDOFs.data()[i] = alphai;
1994 
1995  if (POWER_SMOOTHNESS_INDICATOR==0)
1996  psi[i] = 1.0;
1997  else
1998  psi[i] = std::pow(alphai,POWER_SMOOTHNESS_INDICATOR); //NOTE: they use alpha^2 in the paper
1999  }
2000  }
2001  // END OF COMPUTING SMOOTHNESS INDICATOR
2002 
2004  // ** LOOP IN DOFs FOR EDGE BASED TERMS ** //
2006  ij=0;
2007  for (int i=0; i<numDOFs; i++)
2008  {
2009  double solni = u_dof_old.data()[i]; // solution at time tn for the ith DOF
2010  double ith_dissipative_term = 0;
2011  double ith_flux_term = 0;
2012  double dLowii = 0.;
2013 
2014  // loop over the sparsity pattern of the i-th DOF
2015  for (int offset=csrRowIndeces_DofLoops.data()[i]; offset<csrRowIndeces_DofLoops.data()[i+1]; offset++)
2016  {
2017  int j = csrColumnOffsets_DofLoops.data()[offset];
2018  double solnj = u_dof_old.data()[j]; // solution at time tn for the jth DOF
2019  double dLowij, dHij, dEVij;
2020 
2021  ith_flux_term += TransportMatrix[ij]*solnj;
2022  if (i != j) //NOTE: there is really no need to check for i!=j (see formula for ith_dissipative_term)
2023  {
2024  // first-order dissipative operator
2025  dLowij = fmax(fabs(TransportMatrix[ij]),fabs(TransposeTransportMatrix[ij]));
2026  if (STABILIZATION_TYPE==1) //EV Stab
2027  {
2028  // high-order (entropy viscosity) dissipative operator
2029  dEVij = fmax(fabs(global_entropy_residual[i]),fabs(global_entropy_residual[j]));
2030  dHij = fmin(dLowij,dEVij);
2031  }
2032  else // smoothness based stab
2033  dHij = dLowij*fmax(psi[i],psi[j]); // enhance the order to 2nd order. No EV
2034 
2035  // compute dissipative term
2036  ith_dissipative_term += dHij*(solnj-solni);
2037  dLowii -= dLowij;
2038  }
2039  //update ij
2040  ij+=1;
2041  }
2042  // update residual
2043  double mi = ML.data()[i];
2044  // compute edge_based_cfl
2045  edge_based_cfl.data()[i] = 2.*fabs(dLowii)/mi;
2046 
2047  if (LUMPED_MASS_MATRIX==1)
2048  {
2049  //globalResidual.data()[i] = mi*(u_dof.data()[i] - u_dof_old.data()[i]) + dt*(ith_flux_term - ith_dissipative_term);
2050  globalResidual.data()[i] = u_dof_old.data()[i] - dt/mi*(ith_flux_term - ith_dissipative_term);
2051  }
2052  else
2053  globalResidual.data()[i] += dt*(ith_flux_term - ith_dissipative_term);
2054  }
2055 
2057  // BOUNDARY CONDITIONS //
2059  for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
2060  {
2061  register int ebN = exteriorElementBoundariesArray.data()[ebNE],
2062  eN = elementBoundaryElementsArray.data()[ebN*2+0],
2063  ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
2064  eN_nDOF_trial_element = eN*nDOF_trial_element;
2065  register double elementResidual_u[nDOF_test_element];
2066  for (int i=0;i<nDOF_test_element;i++)
2067  elementResidual_u[i]=0.0;
2068  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
2069  {
2070  register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
2071  ebNE_kb_nSpace = ebNE_kb*nSpace,
2072  ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
2073  ebN_local_kb_nSpace = ebN_local_kb*nSpace;
2074  register double u_ext=0.0,
2075  grad_u_ext[nSpace],
2076  m_ext=0.0,
2077  dm_ext=0.0,
2078  H_ext=0.0,
2079  dH_ext[nSpace],
2080  bc_u_ext=0.0,
2081  bc_grad_u_ext[nSpace],
2082  bc_m_ext=0.0,
2083  bc_dm_ext=0.0,
2084  bc_H_ext=0.0,
2085  bc_dH_ext[nSpace],
2086  jac_ext[nSpace*nSpace],
2087  jacDet_ext,
2088  jacInv_ext[nSpace*nSpace],
2089  boundaryJac[nSpace*(nSpace-1)],
2090  metricTensor[(nSpace-1)*(nSpace-1)],
2091  metricTensorDetSqrt,
2092  dS,
2093  u_test_dS[nDOF_test_element],
2094  u_grad_trial_trace[nDOF_trial_element*nSpace],
2095  normal[nSpace],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
2096  G[nSpace*nSpace],G_dd_G,tr_G,flux_ext;
2097  //
2098  //calculate the solution and gradients at quadrature points
2099  //
2100  //compute information about mapping from reference element to physical element
2101  ck.calculateMapping_elementBoundary(eN,
2102  ebN_local,
2103  kb,
2104  ebN_local_kb,
2105  mesh_dof.data(),
2106  mesh_l2g.data(),
2107  mesh_trial_trace_ref.data(),
2108  mesh_grad_trial_trace_ref.data(),
2109  boundaryJac_ref.data(),
2110  jac_ext,
2111  jacDet_ext,
2112  jacInv_ext,
2113  boundaryJac,
2114  metricTensor,
2115  metricTensorDetSqrt,
2116  normal_ref.data(),
2117  normal,
2118  x_ext,y_ext,z_ext);
2119  ck.calculateMappingVelocity_elementBoundary(eN,
2120  ebN_local,
2121  kb,
2122  ebN_local_kb,
2123  mesh_velocity_dof.data(),
2124  mesh_l2g.data(),
2125  mesh_trial_trace_ref.data(),
2126  xt_ext,yt_ext,zt_ext,
2127  normal,
2128  boundaryJac,
2129  metricTensor,
2130  integralScaling);
2131  dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref.data()[kb];
2132  //get the metric tensor
2133  //cek todo use symmetry
2134  ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
2135  //compute shape and solution information
2136  //shape
2137  ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,u_grad_trial_trace);
2138  //solution and gradients
2139  ck.valFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],u_ext);
2140  ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial_trace,grad_u_ext);
2141  //precalculate test function products with integration weights
2142  for (int j=0;j<nDOF_trial_element;j++)
2143  u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
2144  //
2145  //load the boundary values
2146  //
2147  bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+(1-isDOFBoundary_u.data()[ebNE_kb])*ebqe_rd_u_ext.data()[ebNE_kb];
2148  //
2149  //calculate the pde coefficients using the solution and the boundary values for the solution
2150  //
2151  evaluateCoefficients(&ebqe_velocity_ext.data()[ebNE_kb_nSpace],
2152  u_ext,
2153  grad_u_ext,
2154  m_ext,
2155  dm_ext,
2156  H_ext,
2157  dH_ext);
2158  evaluateCoefficients(&ebqe_velocity_ext.data()[ebNE_kb_nSpace],
2159  bc_u_ext,
2160  bc_grad_u_ext,
2161  bc_m_ext,
2162  bc_dm_ext,
2163  bc_H_ext,
2164  bc_dH_ext);
2165  //
2166  //moving mesh
2167  //
2168  double velocity_ext[nSpace];
2169  double mesh_velocity[3];
2170  mesh_velocity[0] = xt_ext;
2171  mesh_velocity[1] = yt_ext;
2172  mesh_velocity[2] = zt_ext;
2173  for (int I=0;I<nSpace;I++)
2174  velocity_ext[I] = - MOVING_DOMAIN*mesh_velocity[I];
2175  //
2176  //calculate the numerical fluxes
2177  //
2178  exteriorNumericalFlux(normal,
2179  bc_u_ext,
2180  u_ext,
2181  dH_ext,
2182  velocity_ext,
2183  flux_ext);
2184  ebqe_u.data()[ebNE_kb] = u_ext;
2185 
2186  // gradient //
2187 
2188  register double norm = 1.0e-8;
2189  for (int I=0;I<nSpace;I++)
2190  norm += grad_u_ext[I]*grad_u_ext[I];
2191  norm = sqrt(norm);
2192  for (int I=0;I<nSpace;I++)
2193  ebqe_grad_u.data()[ebNE_kb_nSpace+I] = grad_u_ext[I]/norm;
2194 
2195 
2196  //
2197  //update residuals
2198  //
2199  for (int i=0;i<nDOF_test_element;i++)
2200  {
2201  //int ebNE_kb_i = ebNE_kb*nDOF_test_element+i;
2202  elementResidual_u[i] += ck.ExteriorElementBoundaryFlux(flux_ext,u_test_dS[i]);
2203  }//i
2204  }//kb
2205  //
2206  //update the element and global residual storage
2207  //
2208  for (int i=0;i<nDOF_test_element;i++)
2209  {
2210  int eN_i = eN*nDOF_test_element+i;
2211  double mi = ML.data()[offset_u+stride_u*u_l2g.data()[eN_i]];
2212  if (LUMPED_MASS_MATRIX==1)
2213  globalResidual.data()[offset_u+stride_u*r_l2g.data()[eN_i]] -= dt/mi*elementResidual_u[i];
2214  else
2215  globalResidual.data()[offset_u+stride_u*r_l2g.data()[eN_i]] += dt*elementResidual_u[i];
2216  }//i
2217  }//ebNE
2219  }
2220 
2221  void calculateMassMatrix(arguments_dict& args)
2222  {
2223  double dt = args.scalar<double>("dt");
2224  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
2225  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
2226  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
2227  xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
2228  double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
2229  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
2230  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
2231  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
2232  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
2233  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
2234  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
2235  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
2236  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
2237  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
2238  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
2239  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
2240  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
2241  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
2242  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
2243  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
2244  int nElements_global = args.scalar<int>("nElements_global");
2245  double useMetrics = args.scalar<double>("useMetrics");
2246  double alphaBDF = args.scalar<double>("alphaBDF");
2247  int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
2248  double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
2249  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
2250  xt::pyarray<int>& r_l2g = args.array<int>("r_l2g");
2251  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
2252  int degree_polynomial = args.scalar<int>("degree_polynomial");
2253  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
2254  xt::pyarray<double>& velocity = args.array<double>("velocity");
2255  xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
2256  xt::pyarray<double>& cfl = args.array<double>("cfl");
2257  xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
2258  xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
2259  xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
2260  xt::pyarray<double>& globalJacobian = args.array<double>("globalJacobian");
2261  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
2262  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
2263  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
2264  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
2265  xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
2266  xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
2267  xt::pyarray<double>& ebqe_rd_u_ext = args.array<double>("ebqe_rd_u_ext");
2268  xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
2269  xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
2270  int PURE_BDF = args.scalar<int>("PURE_BDF");
2271  int LUMPED_MASS_MATRIX = args.scalar<int>("LUMPED_MASS_MATRIX");
2272  double Ct_sge = 4.0;
2273 
2274  //
2275  //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
2276  //
2277  for(int eN=0;eN<nElements_global;eN++)
2278  {
2279  register double elementJacobian_u_u[nDOF_test_element][nDOF_trial_element];
2280  for (int i=0;i<nDOF_test_element;i++)
2281  for (int j=0;j<nDOF_trial_element;j++)
2282  {
2283  elementJacobian_u_u[i][j]=0.0;
2284  }
2285  for (int k=0;k<nQuadraturePoints_element;k++)
2286  {
2287  int eN_k = eN*nQuadraturePoints_element+k, //index to a scalar at a quadrature point
2288  eN_k_nSpace = eN_k*nSpace,
2289  eN_nDOF_trial_element = eN*nDOF_trial_element; //index to a vector at a quadrature point
2290 
2291  //declare local storage
2292  register double u=0.0,
2293  grad_u[nSpace],
2294  m=0.0,dm=0.0,
2295  H=0.0,dH[nSpace],
2296  f[nSpace],df[nSpace],//MOVING_MESH
2297  m_t=0.0,dm_t=0.0,
2298  dpdeResidual_u_u[nDOF_trial_element],
2299  Lstar_u[nDOF_test_element],
2300  dsubgridError_u_u[nDOF_trial_element],
2301  tau=0.0,tau0=0.0,tau1=0.0,
2302  jac[nSpace*nSpace],
2303  jacDet,
2304  jacInv[nSpace*nSpace],
2305  u_grad_trial[nDOF_trial_element*nSpace],
2306  dV,
2307  u_test_dV[nDOF_test_element],
2308  u_grad_test_dV[nDOF_test_element*nSpace],
2309  x,y,z,xt,yt,zt,
2310  G[nSpace*nSpace],G_dd_G,tr_G;
2311  //
2312  //calculate solution and gradients at quadrature points
2313  //
2314  //get jacobian, etc for mapping reference element
2315  ck.calculateMapping_element(eN,
2316  k,
2317  mesh_dof.data(),
2318  mesh_l2g.data(),
2319  mesh_trial_ref.data(),
2320  mesh_grad_trial_ref.data(),
2321  jac,
2322  jacDet,
2323  jacInv,
2324  x,y,z);
2325  ck.calculateMappingVelocity_element(eN,
2326  k,
2327  mesh_velocity_dof.data(),
2328  mesh_l2g.data(),
2329  mesh_trial_ref.data(),
2330  xt,yt,zt);
2331  //get the physical integration weight
2332  dV = fabs(jacDet)*dV_ref.data()[k];
2333  ck.calculateG(jacInv,G,G_dd_G,tr_G);
2334  //get the trial function gradients
2335  ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
2336  //get the solution
2337  ck.valFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],u);
2338  //get the solution gradients
2339  ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial,grad_u);
2340  //precalculate test function products with integration weights
2341  for (int j=0;j<nDOF_trial_element;j++)
2342  {
2343  u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
2344  for (int I=0;I<nSpace;I++)
2345  {
2346  u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
2347  }
2348  }
2349  //
2350  //calculate pde coefficients and derivatives at quadrature points
2351  //
2352  evaluateCoefficients(&velocity.data()[eN_k_nSpace],
2353  u,
2354  grad_u,
2355  m,
2356  dm,
2357  H,
2358  dH);
2359  //
2360  //moving mesh
2361  //
2362  double mesh_velocity[3];
2363  mesh_velocity[0] = xt;
2364  mesh_velocity[1] = yt;
2365  mesh_velocity[2] = zt;
2366  for (int I=0;I<nSpace;I++)
2367  {
2368  f[I] = -MOVING_DOMAIN*m*mesh_velocity[I];
2369  df[I] = -MOVING_DOMAIN*dm*mesh_velocity[I];
2370  }
2371  //
2372  //calculate time derivatives
2373  //
2374  ck.bdf(alphaBDF,
2375  q_m_betaBDF.data()[eN_k],
2376  m,
2377  dm,
2378  m_t,
2379  dm_t);
2380  //
2381  //calculate subgrid error contribution to the Jacobian (strong residual, adjoint, jacobian of strong residual)
2382  //
2383  //calculate the adjoint times the test functions
2384  for (int i=0;i<nDOF_test_element;i++)
2385  {
2386  //int eN_k_i_nSpace = (eN_k*nDOF_trial_element+i)*nSpace;
2387  //Lstar_u[i]=ck.Hamiltonian_adjoint(dH,&u_grad_test_dV[eN_k_i_nSpace]);
2388  register int i_nSpace = i*nSpace;
2389  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]);
2390 
2391  }
2392  //calculate the Jacobian of strong residual
2393  for (int j=0;j<nDOF_trial_element;j++)
2394  {
2395  //int eN_k_j=eN_k*nDOF_trial_element+j;
2396  //int eN_k_j_nSpace = eN_k_j*nSpace;
2397  int j_nSpace = j*nSpace;
2398  dpdeResidual_u_u[j]=ck.MassJacobian_strong(dm_t,u_trial_ref.data()[k*nDOF_trial_element+j]) +
2399  ck.HamiltonianJacobian_strong(dH,&u_grad_trial[j_nSpace]) +
2400  MOVING_DOMAIN*ck.AdvectionJacobian_strong(df,&u_grad_trial[j_nSpace]);
2401 
2402  }
2403  //tau and tau*Res
2404  double subgridErrorVelocity[nSpace];
2405  for (int I=0;I<nSpace;I++)
2406  subgridErrorVelocity[I] = dH[I] - MOVING_DOMAIN*df[I];
2407 
2408  calculateSubgridError_tau(elementDiameter.data()[eN],
2409  dm_t,
2410  subgridErrorVelocity,//dH,
2411  cfl.data()[eN_k],
2412  tau0);
2413 
2414  calculateSubgridError_tau(Ct_sge,
2415  G,
2416  dm_t,
2417  subgridErrorVelocity,//dH,
2418  tau1,
2419  cfl.data()[eN_k]);
2420  tau = useMetrics*tau1+(1.0-useMetrics)*tau0;
2421 
2422  for(int j=0;j<nDOF_trial_element;j++)
2423  dsubgridError_u_u[j] = -tau*dpdeResidual_u_u[j];
2424 
2425  for(int i=0;i<nDOF_test_element;i++)
2426  {
2427  //int eN_k_i=eN_k*nDOF_test_element+i;
2428  //int eN_k_i_nSpace=eN_k_i*nSpace;
2429  for(int j=0;j<nDOF_trial_element;j++)
2430  {
2431  //int eN_k_j=eN_k*nDOF_trial_element+j;
2432  //int eN_k_j_nSpace = eN_k_j*nSpace;
2433  int j_nSpace = j*nSpace;
2434  int i_nSpace = i*nSpace;
2435  if (LUMPED_MASS_MATRIX==1)
2436  {
2437  if (i==j)
2438  elementJacobian_u_u[i][j] += u_test_dV[i];
2439  }
2440  else
2441  {
2442  elementJacobian_u_u[i][j] +=
2443  dt*ck.MassJacobian_weak(dm_t,u_trial_ref.data()[k*nDOF_trial_element+j],u_test_dV[i]);
2444  }
2445  }//j
2446  }//i
2447  }//k
2448  //
2449  //load into element Jacobian into global Jacobian
2450  //
2451  for (int i=0;i<nDOF_test_element;i++)
2452  {
2453  int eN_i = eN*nDOF_test_element+i;
2454  for (int j=0;j<nDOF_trial_element;j++)
2455  {
2456  int eN_i_j = eN_i*nDOF_trial_element+j;
2457  globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_u_u.data()[eN_i_j]] += elementJacobian_u_u[i][j];
2458  }//j
2459  }//i
2460  }//elements
2461  }//computeMassMatrix
2462 
2463  void calculateSmoothingMatrix(arguments_dict& args)
2464  {
2465  double dt = args.scalar<double>("dt");
2466  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
2467  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
2468  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
2469  xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
2470  double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
2471  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
2472  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
2473  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
2474  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
2475  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
2476  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
2477  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
2478  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
2479  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
2480  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
2481  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
2482  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
2483  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
2484  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
2485  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
2486  int nElements_global = args.scalar<int>("nElements_global");
2487  double useMetrics = args.scalar<double>("useMetrics");
2488  double alphaBDF = args.scalar<double>("alphaBDF");
2489  int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
2490  double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
2491  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
2492  xt::pyarray<int>& r_l2g = args.array<int>("r_l2g");
2493  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
2494  int degree_polynomial = args.scalar<int>("degree_polynomial");
2495  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
2496  xt::pyarray<double>& velocity = args.array<double>("velocity");
2497  xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
2498  xt::pyarray<double>& cfl = args.array<double>("cfl");
2499  xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
2500  xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
2501  xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
2502  xt::pyarray<double>& globalJacobian = args.array<double>("globalJacobian");
2503  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
2504  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
2505  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
2506  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
2507  xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
2508  xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
2509  xt::pyarray<double>& ebqe_rd_u_ext = args.array<double>("ebqe_rd_u_ext");
2510  xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
2511  xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
2512  double he = args.scalar<double>("he");
2513  double Ct_sge = 4.0;
2514 
2515  //
2516  //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
2517  //
2518  for(int eN=0;eN<nElements_global;eN++)
2519  {
2520  register double elementJacobian_u_u[nDOF_test_element][nDOF_trial_element];
2521  for (int i=0;i<nDOF_test_element;i++)
2522  for (int j=0;j<nDOF_trial_element;j++)
2523  {
2524  elementJacobian_u_u[i][j]=0.0;
2525  }
2526  for (int k=0;k<nQuadraturePoints_element;k++)
2527  {
2528  int eN_k = eN*nQuadraturePoints_element+k, //index to a scalar at a quadrature point
2529  eN_k_nSpace = eN_k*nSpace,
2530  eN_nDOF_trial_element = eN*nDOF_trial_element; //index to a vector at a quadrature point
2531 
2532  //declare local storage
2533  register double u=0.0,
2534  grad_u[nSpace],
2535  m=0.0,dm=0.0,
2536  H=0.0,dH[nSpace],
2537  f[nSpace],df[nSpace],//MOVING_MESH
2538  m_t=0.0,dm_t=0.0,
2539  dpdeResidual_u_u[nDOF_trial_element],
2540  Lstar_u[nDOF_test_element],
2541  dsubgridError_u_u[nDOF_trial_element],
2542  tau=0.0,tau0=0.0,tau1=0.0,
2543  jac[nSpace*nSpace],
2544  jacDet,
2545  jacInv[nSpace*nSpace],
2546  u_grad_trial[nDOF_trial_element*nSpace],
2547  dV,
2548  u_test_dV[nDOF_test_element],
2549  u_grad_test_dV[nDOF_test_element*nSpace],
2550  x,y,z,xt,yt,zt,
2551  G[nSpace*nSpace],G_dd_G,tr_G;
2552  //
2553  //calculate solution and gradients at quadrature points
2554  //
2555  //get jacobian, etc for mapping reference element
2556  ck.calculateMapping_element(eN,
2557  k,
2558  mesh_dof.data(),
2559  mesh_l2g.data(),
2560  mesh_trial_ref.data(),
2561  mesh_grad_trial_ref.data(),
2562  jac,
2563  jacDet,
2564  jacInv,
2565  x,y,z);
2566  ck.calculateMappingVelocity_element(eN,
2567  k,
2568  mesh_velocity_dof.data(),
2569  mesh_l2g.data(),
2570  mesh_trial_ref.data(),
2571  xt,yt,zt);
2572  //get the physical integration weight
2573  dV = fabs(jacDet)*dV_ref.data()[k];
2574  ck.calculateG(jacInv,G,G_dd_G,tr_G);
2575  //get the trial function gradients
2576  ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
2577  //get the solution
2578  ck.valFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],u);
2579  //get the solution gradients
2580  ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial,grad_u);
2581  //precalculate test function products with integration weights
2582  for (int j=0;j<nDOF_trial_element;j++)
2583  {
2584  u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
2585  for (int I=0;I<nSpace;I++)
2586  {
2587  u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
2588  }
2589  }
2590  //
2591  //calculate pde coefficients and derivatives at quadrature points
2592  //
2593  evaluateCoefficients(&velocity.data()[eN_k_nSpace],
2594  u,
2595  grad_u,
2596  m,
2597  dm,
2598  H,
2599  dH);
2600  //
2601  //moving mesh
2602  //
2603  double mesh_velocity[3];
2604  mesh_velocity[0] = xt;
2605  mesh_velocity[1] = yt;
2606  mesh_velocity[2] = zt;
2607  for (int I=0;I<nSpace;I++)
2608  {
2609  f[I] = -MOVING_DOMAIN*m*mesh_velocity[I];
2610  df[I] = -MOVING_DOMAIN*dm*mesh_velocity[I];
2611  }
2612  //
2613  //calculate time derivatives
2614  //
2615  ck.bdf(alphaBDF,
2616  q_m_betaBDF.data()[eN_k],
2617  m,
2618  dm,
2619  m_t,
2620  dm_t);
2621  //
2622  //calculate subgrid error contribution to the Jacobian (strong residual, adjoint, jacobian of strong residual)
2623  //
2624  //calculate the adjoint times the test functions
2625  for (int i=0;i<nDOF_test_element;i++)
2626  {
2627  //int eN_k_i_nSpace = (eN_k*nDOF_trial_element+i)*nSpace;
2628  //Lstar_u[i]=ck.Hamiltonian_adjoint(dH,&u_grad_test_dV[eN_k_i_nSpace]);
2629  register int i_nSpace = i*nSpace;
2630  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]);
2631 
2632  }
2633  //calculate the Jacobian of strong residual
2634  for (int j=0;j<nDOF_trial_element;j++)
2635  {
2636  //int eN_k_j=eN_k*nDOF_trial_element+j;
2637  //int eN_k_j_nSpace = eN_k_j*nSpace;
2638  int j_nSpace = j*nSpace;
2639  dpdeResidual_u_u[j]=ck.MassJacobian_strong(dm_t,u_trial_ref.data()[k*nDOF_trial_element+j]) +
2640  ck.HamiltonianJacobian_strong(dH,&u_grad_trial[j_nSpace]) +
2641  MOVING_DOMAIN*ck.AdvectionJacobian_strong(df,&u_grad_trial[j_nSpace]);
2642 
2643  }
2644  //tau and tau*Res
2645  double subgridErrorVelocity[nSpace];
2646  for (int I=0;I<nSpace;I++)
2647  subgridErrorVelocity[I] = dH[I] - MOVING_DOMAIN*df[I];
2648 
2649  calculateSubgridError_tau(elementDiameter.data()[eN],
2650  dm_t,
2651  subgridErrorVelocity,//dH,
2652  cfl.data()[eN_k],
2653  tau0);
2654 
2655  calculateSubgridError_tau(Ct_sge,
2656  G,
2657  dm_t,
2658  subgridErrorVelocity,//dH,
2659  tau1,
2660  cfl.data()[eN_k]);
2661  tau = useMetrics*tau1+(1.0-useMetrics)*tau0;
2662 
2663  for(int j=0;j<nDOF_trial_element;j++)
2664  dsubgridError_u_u[j] = -tau*dpdeResidual_u_u[j];
2665 
2666  for(int i=0;i<nDOF_test_element;i++)
2667  {
2668  //int eN_k_i=eN_k*nDOF_test_element+i;
2669  //int eN_k_i_nSpace=eN_k_i*nSpace;
2670  for(int j=0;j<nDOF_trial_element;j++)
2671  {
2672  //int eN_k_j=eN_k*nDOF_trial_element+j;
2673  //int eN_k_j_nSpace = eN_k_j*nSpace;
2674  int j_nSpace = j*nSpace;
2675  int i_nSpace = i*nSpace;
2676  elementJacobian_u_u[i][j] +=
2677  dt*ck.MassJacobian_weak(dm_t,u_trial_ref.data()[k*nDOF_trial_element+j],u_test_dV[i]) +
2678  ck.NumericalDiffusionJacobian(std::pow(he,2),&u_grad_trial[j_nSpace],&u_grad_test_dV[i_nSpace]);
2679  }//j
2680  }//i
2681  }//k
2682  //
2683  //load into element Jacobian into global Jacobian
2684  //
2685  for (int i=0;i<nDOF_test_element;i++)
2686  {
2687  int eN_i = eN*nDOF_test_element+i;
2688  for (int j=0;j<nDOF_trial_element;j++)
2689  {
2690  int eN_i_j = eN_i*nDOF_trial_element+j;
2691  globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_u_u.data()[eN_i_j]] += elementJacobian_u_u[i][j];
2692  }//j
2693  }//i
2694  }//elements
2695  }//computeMassMatrix
2696  };//NCLS
2697 
2698  inline NCLS_base* newNCLS(int nSpaceIn,
2699  int nQuadraturePoints_elementIn,
2700  int nDOF_mesh_trial_elementIn,
2701  int nDOF_trial_elementIn,
2702  int nDOF_test_elementIn,
2703  int nQuadraturePoints_elementBoundaryIn,
2704  int CompKernelFlag)
2705  {
2706  if (nSpaceIn == 2)
2707  return proteus::chooseAndAllocateDiscretization2D<NCLS_base,NCLS,CompKernel>(nSpaceIn,
2708  nQuadraturePoints_elementIn,
2709  nDOF_mesh_trial_elementIn,
2710  nDOF_trial_elementIn,
2711  nDOF_test_elementIn,
2712  nQuadraturePoints_elementBoundaryIn,
2713  CompKernelFlag);
2714  else
2715  return proteus::chooseAndAllocateDiscretization<NCLS_base,NCLS,CompKernel>(nSpaceIn,
2716  nQuadraturePoints_elementIn,
2717  nDOF_mesh_trial_elementIn,
2718  nDOF_trial_elementIn,
2719  nDOF_test_elementIn,
2720  nQuadraturePoints_elementBoundaryIn,
2721  CompKernelFlag);
2722  /* return proteus::chooseAndAllocateDiscretization<NCLS_base,NCLS>(nSpaceIn, */
2723  /* nQuadraturePoints_elementIn, */
2724  /* nDOF_mesh_trial_elementIn, */
2725  /* nDOF_trial_elementIn, */
2726  /* nDOF_test_elementIn, */
2727  /* nQuadraturePoints_elementBoundaryIn, */
2728  /* CompKernelFlag); */
2729  }
2730 }//proteus
2731 #endif
proteus::NCLS::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: NCLS.h:144
proteus::NCLS_base::L2_norm_per_node
std::valarray< double > L2_norm_per_node
Definition: NCLS.h:43
proteus::NCLS::calculateResidual
void calculateResidual(arguments_dict &args)
Definition: NCLS.h:210
proteus::NCLS::calculateMassMatrix
void calculateMassMatrix(arguments_dict &args)
Definition: NCLS.h:2221
IS_BETAij_ONE
#define IS_BETAij_ONE
Definition: NCLS.h:14
sgn
double sgn(double val)
Definition: RANS3PF2D.h:20
neg
double neg(double a)
Definition: testFMMandFSW.cpp:10
proteus::NCLS_base::gy
std::valarray< double > gy
Definition: NCLS.h:47
proteus::NCLS_base::alpha_numerator_pos
std::valarray< double > alpha_numerator_pos
Definition: NCLS.h:48
proteus::NCLS_base::etaMax
std::valarray< double > etaMax
Definition: NCLS.h:45
proteus::NCLS_base::alpha_denominator_pos
std::valarray< double > alpha_denominator_pos
Definition: NCLS.h:48
proteus::NCLS
Definition: NCLS.h:68
proteus::NCLS_base::global_entropy_residual
std::valarray< double > global_entropy_residual
Definition: NCLS.h:46
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
ModelFactory.h
CompKernel.h
proteus::arguments_dict::scalar
T & scalar(const std::string &key)
proteus::NCLS_base::calculateRhsSmoothing
virtual double calculateRhsSmoothing(arguments_dict &args)=0
proteus::NCLS_base::gx
std::valarray< double > gx
Definition: NCLS.h:47
proteus::NCLS::sign
double sign(const double phi, const double eps)
Definition: NCLS.h:192
proteus::arguments_dict::array
xt::pyarray< T > & array(const std::string &key)
proteus::NCLS::calculateJacobian
void calculateJacobian(arguments_dict &args)
Definition: NCLS.h:690
H
Double H
Definition: Headers.h:65
proteus::NCLS::ck
CompKernelType ck
Definition: NCLS.h:71
proteus::NCLS_base::eta
std::valarray< double > eta
Definition: NCLS.h:47
proteus::NCLS_base::calculateResidual_entropy_viscosity
virtual void calculateResidual_entropy_viscosity(arguments_dict &args)=0
v
Double v
Definition: Headers.h:95
proteus::NCLS::nDOF_test_X_trial_element
const int nDOF_test_X_trial_element
Definition: NCLS.h:70
proteus::NCLS_base
Definition: NCLS.h:40
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::DENTROPY_LOG
double DENTROPY_LOG(const double &phi, const double &phiL, const double &phiR)
Definition: NCLS.h:32
proteus::NCLS::calculateSubgridError_tau
void calculateSubgridError_tau(const double &elementDiameter, const double &dmt, const double dH[nSpace], double &cfl, double &tau)
Definition: NCLS.h:110
proteus::NCLS_base::calculateSmoothingMatrix
virtual void calculateSmoothingMatrix(arguments_dict &args)=0
z
Double * z
Definition: Headers.h:49
u
Double u
Definition: Headers.h:89
proteus::phi
double phi(const double &g, const double &h, const double &hL, const double &hR, const double &uL, const double &uR)
Definition: SW2DCV.h:62
proteus::NCLS_base::TransposeTransportMatrix
std::valarray< double > TransposeTransportMatrix
Definition: NCLS.h:44
proteus::NCLS::exteriorNumericalFluxDerivative
void exteriorNumericalFluxDerivative(const double n[nSpace], const double velocity[nSpace], const double velocity_movingDomain[nSpace], double &dflux)
Definition: NCLS.h:169
proteus::NCLS_base::alpha_denominator_neg
std::valarray< double > alpha_denominator_neg
Definition: NCLS.h:48
proteus::ENTROPY_LOG
double ENTROPY_LOG(const double &phi, const double &phiL, const double &phiR)
Definition: NCLS.h:29
xt
Definition: AddedMass.cpp:7
POWER_SMOOTHNESS_INDICATOR
#define POWER_SMOOTHNESS_INDICATOR
Definition: NCLS.h:13
proteus::NCLS::calculateResidual_entropy_viscosity
void calculateResidual_entropy_viscosity(arguments_dict &args)
Definition: NCLS.h:1542
proteus::NCLS::calculateRhsSmoothing
double calculateRhsSmoothing(arguments_dict &args)
Definition: NCLS.h:1470
proteus::NCLS_base::etaMin
std::valarray< double > etaMin
Definition: NCLS.h:45
proteus::NCLS::calculateRedistancingResidual
double calculateRedistancingResidual(arguments_dict &args)
Definition: NCLS.h:1238
proteus::NCLS_base::calculateResidual
virtual void calculateResidual(arguments_dict &args)=0
proteus::NCLS_base::psi
std::valarray< double > psi
Definition: NCLS.h:45
proteus::NCLS_base::calculateMassMatrix
virtual void calculateMassMatrix(arguments_dict &args)=0
proteus
Definition: ADR.h:17
proteus::NCLS::calculateWaterline
void calculateWaterline(arguments_dict &args)
Definition: NCLS.h:1110
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::NCLS_base::gz
std::valarray< double > gz
Definition: NCLS.h:47
proteus::NCLS_base::TransportMatrix
std::valarray< double > TransportMatrix
Definition: NCLS.h:44
proteus::NCLS_base::calculateWaterline
virtual void calculateWaterline(arguments_dict &args)=0
proteus::arguments_dict
Definition: ArgumentsDict.h:70
proteus::NCLS_base::alpha_numerator_neg
std::valarray< double > alpha_numerator_neg
Definition: NCLS.h:48
proteus::NCLS_base::~NCLS_base
virtual ~NCLS_base()
Definition: NCLS.h:49
pos
double pos(double a)
Definition: testFMMandFSW.cpp:8
proteus::NCLS::calculateSmoothingMatrix
void calculateSmoothingMatrix(arguments_dict &args)
Definition: NCLS.h:2463
proteus::NCLS_base::calculateRedistancingResidual
virtual double calculateRedistancingResidual(arguments_dict &args)=0
proteus::newNCLS
NCLS_base * newNCLS(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition: NCLS.h:2698
proteus::NCLS::calculateCFL
void calculateCFL(const double &elementDiameter, const double df[nSpace], double &cfl)
Definition: NCLS.h:96
proteus::NCLS::NCLS
NCLS()
Definition: NCLS.h:72
proteus::NCLS::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: NCLS.h:77
ArgumentsDict.h
proteus::NCLS::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: NCLS.h:129
proteus::NCLS_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)=0
cE
#define cE
Definition: NCLS3P.h:10
Generated on Fri Jul 1 2022 11:17:39 for proteus by doxygen 1.8.20