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