proteus  1.8.1
C/C++/Fortran libraries
VOF.h
Go to the documentation of this file.
1 #ifndef VOF_H
2 #define VOF_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 // Cell based methods:
18 // * SUPG with BDF1 or BDF2 time integration
19 // * Explicit Taylor Galerkin with EV stabilization
20 // Edge based methods.
21 // Low order via D. Kuzmin's
22 // High order methods: Smoothness indicator with MC, EV commutator with MC, D.K with ML
23 // Zalesak's FCT
24 
25 namespace proteus
26 {
27  // Power entropy //
28  inline double ENTROPY(const double& phi, const double& phiL, const double& phiR){
29  return 1./2.*std::pow(fabs(phi),2.);
30  }
31  inline double DENTROPY(const double& phi, const double& phiL, const double& phiR){
32  return fabs(phi)*(phi>=0 ? 1 : -1);
33  }
34  // Log entropy // for level set from 0 to 1
35  inline double ENTROPY_LOG(const double& phi, const double& phiL, const double& phiR){
36  return std::log(fabs((phi-phiL)*(phiR-phi))+1E-14);
37  }
38  inline double DENTROPY_LOG(const double& phi, const double& phiL, const double& phiR){
39  return (phiL+phiR-2*phi)*((phi-phiL)*(phiR-phi)>=0 ? 1 : -1)/(fabs((phi-phiL)*(phiR-phi))+1E-14);
40  }
41 }
42 
43 namespace proteus
44 {
45  class VOF_base
46  {
47  //The base class defining the interface
48  public:
49  std::valarray<double> Rpos, Rneg;
50  std::valarray<double> FluxCorrectionMatrix;
51  std::valarray<double> TransportMatrix, TransposeTransportMatrix;
52  std::valarray<double> psi, eta, global_entropy_residual, boundary_integral;
53  std::valarray<double> maxVel,maxEntRes;
54  virtual ~VOF_base(){}
55  virtual void calculateResidualElementBased(arguments_dict& args)=0;
56  virtual void calculateJacobian(arguments_dict& args)=0;
57  virtual void FCTStep(arguments_dict& args)=0;
58  virtual void calculateResidualEdgeBased(arguments_dict& args)=0;
59  };
60 
61  template<class CompKernelType,
62  int nSpace,
63  int nQuadraturePoints_element,
64  int nDOF_mesh_trial_element,
65  int nDOF_trial_element,
66  int nDOF_test_element,
67  int nQuadraturePoints_elementBoundary>
68  class VOF : public VOF_base
69  {
70  public:
71  const int nDOF_test_X_trial_element;
72  CompKernelType ck;
73  VOF():
74  nDOF_test_X_trial_element(nDOF_test_element*nDOF_trial_element),
75  ck()
76  {}
77 
78  inline
79  void calculateCFL(const double& elementDiameter,
80  const double df[nSpace],
81  double& cfl)
82  {
83  double h,nrm_v;
84  h = elementDiameter;
85  nrm_v=0.0;
86  for(int I=0;I<nSpace;I++)
87  nrm_v+=df[I]*df[I];
88  nrm_v = sqrt(nrm_v);
89  cfl = nrm_v/h;
90  }
91 
92  inline
93  void evaluateCoefficients(const double v[nSpace],
94  const double& u,
95  const double& porosity, //VRANS specific
96  double& m,
97  double& dm,
98  double f[nSpace],
99  double df[nSpace])
100  {
101  m = porosity*u;
102  dm= porosity;
103  for (int I=0; I < nSpace; I++)
104  {
105  f[I] = v[I]*porosity*u;
106  df[I] = v[I]*porosity;
107  }
108  }
109 
110  inline
111  void calculateSubgridError_tau(const double& elementDiameter,
112  const double& dmt,
113  const double dH[nSpace],
114  double& cfl,
115  double& tau)
116  {
117  double h,nrm_v,oneByAbsdt;
118  h = elementDiameter;
119  nrm_v=0.0;
120  for(int I=0;I<nSpace;I++)
121  nrm_v+=dH[I]*dH[I];
122  nrm_v = sqrt(nrm_v);
123  cfl = nrm_v/h;
124  oneByAbsdt = fabs(dmt);
125  tau = 1.0/(2.0*nrm_v/h + oneByAbsdt + 1.0e-8);
126  }
127 
128  inline
129  void calculateSubgridError_tau( const double& Ct_sge,
130  const double G[nSpace*nSpace],
131  const double& A0,
132  const double Ai[nSpace],
133  double& tau_v,
134  double& cfl)
135  {
136  double v_d_Gv=0.0;
137  for(int I=0;I<nSpace;I++)
138  for (int J=0;J<nSpace;J++)
139  v_d_Gv += Ai[I]*G[I*nSpace+J]*Ai[J];
140 
141  tau_v = 1.0/sqrt(Ct_sge*A0*A0 + v_d_Gv + 1.0e-8);
142  }
143 
144  inline
145  void calculateNumericalDiffusion(const double& shockCapturingDiffusion,
146  const double& elementDiameter,
147  const double& strong_residual,
148  const double grad_u[nSpace],
149  double& numDiff)
150  {
151  double h,
152  num,
153  den,
154  n_grad_u;
155  h = elementDiameter;
156  n_grad_u = 0.0;
157  for (int I=0;I<nSpace;I++)
158  n_grad_u += grad_u[I]*grad_u[I];
159  num = shockCapturingDiffusion*0.5*h*fabs(strong_residual);
160  den = sqrt(n_grad_u) + 1.0e-8;
161  numDiff = num/den;
162  }
163 
164  inline
165  void exteriorNumericalAdvectiveFlux(const int& isDOFBoundary_u,
166  const int& isFluxBoundary_u,
167  const double n[nSpace],
168  const double& bc_u,
169  const double& bc_flux_u,
170  const double& u,
171  const double velocity[nSpace],
172  double& flux)
173  {
174 
175  double flow=0.0;
176  for (int I=0; I < nSpace; I++)
177  flow += n[I]*velocity[I];
178  //std::cout<<" isDOFBoundary_u= "<<isDOFBoundary_u<<" flow= "<<flow<<std::endl;
179  if (isDOFBoundary_u == 1)
180  {
181  //std::cout<<"Dirichlet boundary u and bc_u "<<u<<'\t'<<bc_u<<std::endl;
182  if (flow >= 0.0)
183  {
184  flux = u*flow;
185  //flux = flow;
186  }
187  else
188  {
189  flux = bc_u*flow;
190  //flux = flow;
191  }
192  }
193  else if (isFluxBoundary_u == 1)
194  {
195  flux = bc_flux_u;
196  //std::cout<<"Flux boundary flux and flow"<<flux<<'\t'<<flow<<std::endl;
197  }
198  else
199  {
200  //std::cout<<"No BC boundary flux and flow"<<flux<<'\t'<<flow<<std::endl;
201  if (flow >= 0.0)
202  {
203  flux = u*flow;
204  }
205  else
206  {
207  std::cout<<"warning: VOF open boundary with no external trace, setting to zero for inflow"<<std::endl;
208  flux = 0.0;
209  }
210 
211  }
212  }
213 
214  inline
215  void exteriorNumericalAdvectiveFluxDerivative(const int& isDOFBoundary_u,
216  const int& isFluxBoundary_u,
217  const double n[nSpace],
218  const double velocity[nSpace],
219  double& dflux)
220  {
221  double flow=0.0;
222  for (int I=0; I < nSpace; I++)
223  flow += n[I]*velocity[I];
224  //double flow=n[0]*velocity[0]+n[1]*velocity[1]+n[2]*velocity[2];
225  dflux=0.0;//default to no flux
226  if (isDOFBoundary_u == 1)
227  {
228  if (flow >= 0.0)
229  {
230  dflux = flow;
231  }
232  else
233  {
234  dflux = 0.0;
235  }
236  }
237  else if (isFluxBoundary_u == 1)
238  {
239  dflux = 0.0;
240  }
241  else
242  {
243  if (flow >= 0.0)
244  {
245  dflux = flow;
246  }
247  }
248  }
249 
250  void calculateResidualElementBased(arguments_dict& args)
251  {
252  double dt = args.scalar<double>("dt");
253  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
254  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
255  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
256  xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
257  double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
258  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
259  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
260  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
261  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
262  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
263  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
264  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
265  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
266  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
267  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
268  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
269  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
270  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
271  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
272  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
273  int nElements_global = args.scalar<int>("nElements_global");
274  double useMetrics = args.scalar<double>("useMetrics");
275  double alphaBDF = args.scalar<double>("alphaBDF");
276  int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
277  double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
278  double sc_uref = args.scalar<double>("sc_uref");
279  double sc_alpha = args.scalar<double>("sc_alpha");
280  const xt::pyarray<double>& q_porosity = args.array<double>("q_porosity");
281  const xt::pyarray<double>& porosity_dof = args.array<double>("porosity_dof");
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  double degree_polynomial = args.scalar<double>("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_m_betaBDF = args.array<double>("q_m_betaBDF");
292  xt::pyarray<double>& q_dV = args.array<double>("q_dV");
293  xt::pyarray<double>& q_dV_last = args.array<double>("q_dV_last");
294  xt::pyarray<double>& cfl = args.array<double>("cfl");
295  xt::pyarray<double>& edge_based_cfl = args.array<double>("edge_based_cfl");
296  xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
297  xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
298  int offset_u = args.scalar<int>("offset_u");
299  int stride_u = args.scalar<int>("stride_u");
300  xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
301  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
302  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
303  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
304  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
305  xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
306  const xt::pyarray<double>& ebqe_porosity_ext = args.array<double>("ebqe_porosity_ext");
307  xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
308  xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
309  xt::pyarray<int>& isFluxBoundary_u = args.array<int>("isFluxBoundary_u");
310  xt::pyarray<double>& ebqe_bc_flux_u_ext = args.array<double>("ebqe_bc_flux_u_ext");
311  xt::pyarray<double>& ebqe_phi = args.array<double>("ebqe_phi");
312  double epsFact = args.scalar<double>("epsFact");
313  xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
314  xt::pyarray<double>& ebqe_flux = args.array<double>("ebqe_flux");
315  int stage = args.scalar<int>("stage");
316  xt::pyarray<double>& uTilde_dof = args.array<double>("uTilde_dof");
317  double cE = args.scalar<double>("cE");
318  double cMax = args.scalar<double>("cMax");
319  double cK = args.scalar<double>("cK");
320  double uL = args.scalar<double>("uL");
321  double uR = args.scalar<double>("uR");
322  int numDOFs = args.scalar<int>("numDOFs");
323  int NNZ = args.scalar<int>("NNZ");
324  xt::pyarray<int>& csrRowIndeces_DofLoops = args.array<int>("csrRowIndeces_DofLoops");
325  xt::pyarray<int>& csrColumnOffsets_DofLoops = args.array<int>("csrColumnOffsets_DofLoops");
326  xt::pyarray<int>& csrRowIndeces_CellLoops = args.array<int>("csrRowIndeces_CellLoops");
327  xt::pyarray<int>& csrColumnOffsets_CellLoops = args.array<int>("csrColumnOffsets_CellLoops");
328  xt::pyarray<int>& csrColumnOffsets_eb_CellLoops = args.array<int>("csrColumnOffsets_eb_CellLoops");
329  xt::pyarray<double>& ML = args.array<double>("ML");
330  int LUMPED_MASS_MATRIX = args.scalar<int>("LUMPED_MASS_MATRIX");
331  int STABILIZATION_TYPE = args.scalar<int>("STABILIZATION_TYPE");
332  int ENTROPY_TYPE = args.scalar<int>("ENTROPY_TYPE");
333  xt::pyarray<double>& uLow = args.array<double>("uLow");
334  xt::pyarray<double>& dLow = args.array<double>("dLow");
335  xt::pyarray<double>& dt_times_dH_minus_dL = args.array<double>("dt_times_dH_minus_dL");
336  xt::pyarray<double>& min_u_bc = args.array<double>("min_u_bc");
337  xt::pyarray<double>& max_u_bc = args.array<double>("max_u_bc");
338  xt::pyarray<double>& quantDOFs = args.array<double>("quantDOFs");
339  double meanEntropy = 0., meanOmega = 0., maxEntropy = -1E10, minEntropy = 1E10;
340  maxVel.resize(nElements_global, 0.0);
341  maxEntRes.resize(nElements_global, 0.0);
342  double Ct_sge = 4.0;
343  //
344  //loop over elements to compute volume integrals and load them into element and global residual
345  //
346  //eN is the element index
347  //eN_k is the quadrature point index for a scalar
348  //eN_k_nSpace is the quadrature point index for a vector
349  //eN_i is the element test function index
350  //eN_j is the element trial function index
351  //eN_k_j is the quadrature point index for a trial function
352  //eN_k_i is the quadrature point index for a trial function
353  for(int eN=0;eN<nElements_global;eN++)
354  {
355  //declare local storage for element residual and initialize
356  register double elementResidual_u[nDOF_test_element];
357  for (int i=0;i<nDOF_test_element;i++)
358  {
359  elementResidual_u[i]=0.0;
360  }//i
361  //loop over quadrature points and compute integrands
362  for (int k=0;k<nQuadraturePoints_element;k++)
363  {
364  //compute indeces and declare local storage
365  register int eN_k = eN*nQuadraturePoints_element+k,
366  eN_k_nSpace = eN_k*nSpace,
367  eN_nDOF_trial_element = eN*nDOF_trial_element;
368  register double
369  entVisc_minus_artComp,
370  u=0.0,un=0.0,
371  grad_u[nSpace],grad_u_old[nSpace],grad_uTilde[nSpace],
372  m=0.0,dm=0.0,
373  H=0.0,Hn=0.0,HTilde=0.0,
374  f[nSpace],fn[nSpace],df[nSpace],
375  m_t=0.0,dm_t=0.0,
376  pdeResidual_u=0.0,
377  Lstar_u[nDOF_test_element],
378  subgridError_u=0.0,
379  tau=0.0,tau0=0.0,tau1=0.0,
380  numDiff0=0.0,numDiff1=0.0,
381  jac[nSpace*nSpace],
382  jacDet,
383  jacInv[nSpace*nSpace],
384  u_grad_trial[nDOF_trial_element*nSpace],
385  u_test_dV[nDOF_trial_element],
386  u_grad_test_dV[nDOF_test_element*nSpace],
387  dV,x,y,z,xt,yt,zt,
388  //VRANS
389  porosity,
390  //
391  G[nSpace*nSpace],G_dd_G,tr_G;//norm_Rv;
392  // //
393  // //compute solution and gradients at quadrature points
394  // //
395  // u=0.0;
396  // for (int I=0;I<nSpace;I++)
397  // {
398  // grad_u[I]=0.0;
399  // }
400  // for (int j=0;j<nDOF_trial_element;j++)
401  // {
402  // int eN_j=eN*nDOF_trial_element+j;
403  // int eN_k_j=eN_k*nDOF_trial_element+j;
404  // int eN_k_j_nSpace = eN_k_j*nSpace;
405  // u += valFromDOF_c(u_dof.data()[u_l2g.data()[eN_j]],u_trial[eN_k_j]);
406  // for (int I=0;I<nSpace;I++)
407  // {
408  // grad_u[I] += gradFromDOF_c(u_dof.data()[u_l2g.data()[eN_j]],u_grad_trial[eN_k_j_nSpace+I]);
409  // }
410  // }
411  ck.calculateMapping_element(eN,
412  k,
413  mesh_dof.data(),
414  mesh_l2g.data(),
415  mesh_trial_ref.data(),
416  mesh_grad_trial_ref.data(),
417  jac,
418  jacDet,
419  jacInv,
420  x,y,z);
421  ck.calculateMappingVelocity_element(eN,
422  k,
423  mesh_velocity_dof.data(),
424  mesh_l2g.data(),
425  mesh_trial_ref.data(),
426  xt,yt,zt);
427  //get the physical integration weight
428  dV = fabs(jacDet)*dV_ref.data()[k];
429  ck.calculateG(jacInv,G,G_dd_G,tr_G);
430  //get the trial function gradients
431  ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],
432  jacInv,
433  u_grad_trial);
434  //get the solution
435  ck.valFromDOF(u_dof.data(),
436  &u_l2g.data()[eN_nDOF_trial_element],
437  &u_trial_ref.data()[k*nDOF_trial_element],
438  u);
439  ck.valFromDOF(u_dof_old.data(),
440  &u_l2g.data()[eN_nDOF_trial_element],
441  &u_trial_ref.data()[k*nDOF_trial_element],
442  un);
443  //get the solution gradients
444  ck.gradFromDOF(u_dof.data(),
445  &u_l2g.data()[eN_nDOF_trial_element],
446  u_grad_trial,
447  grad_u);
448  ck.gradFromDOF(u_dof_old.data(),
449  &u_l2g.data()[eN_nDOF_trial_element],
450  u_grad_trial,
451  grad_u_old);
452  ck.gradFromDOF(uTilde_dof.data(),
453  &u_l2g.data()[eN_nDOF_trial_element],
454  u_grad_trial,
455  grad_uTilde);
456  //precalculate test function products with integration weights
457  for (int j=0;j<nDOF_trial_element;j++)
458  {
459  u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
460  for (int I=0;I<nSpace;I++)
461  {
462  u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
463  }
464  }
465  //VRANS
466  porosity = q_porosity.data()[eN_k];
467  //
468  //
469  //calculate pde coefficients at quadrature points
470  //
471  evaluateCoefficients(&velocity.data()[eN_k_nSpace],
472  u,
473  //VRANS
474  porosity,
475  //
476  m,
477  dm,
478  f,
479  df);
480  //
481  //moving mesh
482  //
483  double mesh_velocity[3];
484  mesh_velocity[0] = xt;
485  mesh_velocity[1] = yt;
486  mesh_velocity[2] = zt;
487 
488  for (int I=0;I<nSpace;I++)
489  {
490  f[I] -= MOVING_DOMAIN*m*mesh_velocity[I];
491  df[I] -= MOVING_DOMAIN*dm*mesh_velocity[I];
492  }
493  //
494  //calculate time derivative at quadrature points
495  //
496  if (q_dV_last.data()[eN_k] <= -100)
497  q_dV_last.data()[eN_k] = dV;
498  q_dV.data()[eN_k] = dV;
499  ck.bdf(alphaBDF,
500  q_m_betaBDF.data()[eN_k]*q_dV_last.data()[eN_k]/dV,//ensure prior mass integral is correct for m_t with BDF1
501  m,
502  dm,
503  m_t,
504  dm_t);
505 
506  if (STABILIZATION_TYPE==1)
507  {
508  double normVel=0., norm_grad_un=0.;
509  for (int I=0;I<nSpace;I++)
510  {
511  Hn += df[I]*grad_u_old[I];
512  HTilde += df[I]*grad_uTilde[I];
513  fn[I] = porosity*df[I]*un-MOVING_DOMAIN*m*mesh_velocity[I];
514  H += df[I]*grad_u[I];
515  normVel += df[I]*df[I];
516  norm_grad_un += grad_u_old[I]*grad_u_old[I];
517  }
518  normVel = std::sqrt(normVel);
519  norm_grad_un = std::sqrt(norm_grad_un)+1E-10;
520 
521  // calculate CFL
522  calculateCFL(elementDiameter.data()[eN]/degree_polynomial,df,cfl.data()[eN_k]);
523 
524 
525  // compute max velocity at cell
526  maxVel[eN] = fmax(normVel,maxVel[eN]);
527 
528  // Strong entropy residual
529  double entRes = (ENTROPY(u,0,1)-ENTROPY(un,0,1))/dt + 0.5*(DENTROPY(u,0,1)*H +
530  DENTROPY(un,0,1)*Hn);
531  maxEntRes[eN] = fmax(maxEntRes[eN],fabs(entRes));
532 
533  // Quantities for normalization factor //
534  meanEntropy += ENTROPY(u,0,1)*dV;
535  meanOmega += dV;
536  maxEntropy = fmax(maxEntropy,ENTROPY(u,0,1));
537  minEntropy = fmin(minEntropy,ENTROPY(u,0,1));
538 
539  // artificial compression
540  double hK=elementDiameter.data()[eN]/degree_polynomial;
541  entVisc_minus_artComp = fmax(1-cK*fmax(un*(1-un),0)/hK/norm_grad_un,0);
542  }
543  else
544  {
545  //
546  //calculate subgrid error (strong residual and adjoint)
547  //
548  //calculate strong residual
549  pdeResidual_u = ck.Mass_strong(m_t) + ck.Advection_strong(df,grad_u);
550  //calculate adjoint
551  for (int i=0;i<nDOF_test_element;i++)
552  {
553  // register int eN_k_i_nSpace = (eN_k*nDOF_trial_element+i)*nSpace;
554  // Lstar_u[i] = ck.Advection_adjoint(df,&u_grad_test_dV[eN_k_i_nSpace]);
555  register int i_nSpace = i*nSpace;
556  Lstar_u[i] = ck.Advection_adjoint(df,&u_grad_test_dV[i_nSpace]);
557  }
558  //calculate tau and tau*Res
559  calculateSubgridError_tau(elementDiameter.data()[eN],dm_t,df,cfl.data()[eN_k],tau0);
560  calculateSubgridError_tau(Ct_sge,
561  G,
562  dm_t,
563  df,
564  tau1,
565  cfl.data()[eN_k]);
566  tau = useMetrics*tau1+(1.0-useMetrics)*tau0;
567 
568  subgridError_u = -tau*pdeResidual_u;
569  //
570  //calculate shock capturing diffusion
571  //
572 
573  ck.calculateNumericalDiffusion(shockCapturingDiffusion,
574  elementDiameter.data()[eN],
575  pdeResidual_u,
576  grad_u,
577  numDiff0);
578  //ck.calculateNumericalDiffusion(shockCapturingDiffusion,G,pdeResidual_u,grad_u_old,numDiff1);
579  ck.calculateNumericalDiffusion(shockCapturingDiffusion,
580  sc_uref,
581  sc_alpha,
582  G,
583  G_dd_G,
584  pdeResidual_u,
585  grad_u,
586  numDiff1);
587  q_numDiff_u.data()[eN_k] = useMetrics*numDiff1+(1.0-useMetrics)*numDiff0;
588  //std::cout<<tau<<" "<<q_numDiff_u.data()[eN_k]<<'\t'<<numDiff0<<'\t'<<numDiff1<<'\t'<<pdeResidual_u<<std::endl;
589 
590  //
591  //update element residual
592  //
593 
594 
595  /* std::cout<<m_t<<'\t'
596  <<f[0]<<'\t'
597  <<f[1]<<'\t'
598  <<df[0]<<'\t'
599  <<df[1]<<'\t'
600  <<subgridError_u<<'\t'
601  <<q_numDiff_u_last.data()[eN_k]<<std::endl;*/
602  }
603 
604  for(int i=0;i<nDOF_test_element;i++)
605  {
606  //register int eN_k_i=eN_k*nDOF_test_element+i,
607  //eN_k_i_nSpace = eN_k_i*nSpace,
608  register int i_nSpace=i*nSpace;
609  if (STABILIZATION_TYPE==1)
610  {
611  if (stage == 1)
612  elementResidual_u[i] +=
613  ck.Mass_weak(dt*m_t,u_test_dV[i]) + // time derivative
614  1./3*dt*ck.Advection_weak(fn,&u_grad_test_dV[i_nSpace]) +
615  1./9*dt*dt*ck.NumericalDiffusion(Hn,df,&u_grad_test_dV[i_nSpace]) +
616  1./3*dt*entVisc_minus_artComp*ck.NumericalDiffusion(q_numDiff_u_last.data()[eN_k],
617  grad_u_old,
618  &u_grad_test_dV[i_nSpace]);
619  // TODO: Add part about moving mesh
620  else //stage 2
621  elementResidual_u[i] +=
622  ck.Mass_weak(dt*m_t,u_test_dV[i]) + // time derivative
623  dt*ck.Advection_weak(fn,&u_grad_test_dV[i_nSpace]) +
624  0.5*dt*dt*ck.NumericalDiffusion(HTilde,df,&u_grad_test_dV[i_nSpace]) +
625  dt*entVisc_minus_artComp*ck.NumericalDiffusion(q_numDiff_u_last.data()[eN_k],
626  grad_u_old,
627  &u_grad_test_dV[i_nSpace]);
628  }
629  else //supg
630  {
631  elementResidual_u[i] +=
632  ck.Mass_weak(m_t,u_test_dV[i]) +
633  ck.Advection_weak(f,&u_grad_test_dV[i_nSpace]) +
634  ck.SubgridError(subgridError_u,Lstar_u[i]) +
635  ck.NumericalDiffusion(q_numDiff_u_last.data()[eN_k],
636  grad_u,
637  &u_grad_test_dV[i_nSpace]);
638  }
639  }//i
640  //
641  //cek/ido todo, get rid of m, since u=m
642  //save momentum for time history and velocity for subgrid error
643  //save solution for other models
644  //
645  q_u.data()[eN_k] = u;
646  q_m.data()[eN_k] = m;
647  }
648  //
649  //load element into global residual and save element residual
650  //
651  for(int i=0;i<nDOF_test_element;i++)
652  {
653  register int eN_i=eN*nDOF_test_element+i;
654  globalResidual.data()[offset_u+stride_u*r_l2g.data()[eN_i]] += elementResidual_u[i];
655  }//i
656  }//elements
657  //
658  //loop over exterior element boundaries to calculate surface integrals and load into element and global residuals
659  //
660  //ebNE is the Exterior element boundary INdex
661  //ebN is the element boundary INdex
662  //eN is the element index
663  for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
664  {
665  register int ebN = exteriorElementBoundariesArray.data()[ebNE],
666  eN = elementBoundaryElementsArray.data()[ebN*2+0],
667  ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
668  eN_nDOF_trial_element = eN*nDOF_trial_element;
669  register double elementResidual_u[nDOF_test_element];
670  for (int i=0;i<nDOF_test_element;i++)
671  {
672  elementResidual_u[i]=0.0;
673  }
674  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
675  {
676  register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
677  ebNE_kb_nSpace = ebNE_kb*nSpace,
678  ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
679  ebN_local_kb_nSpace = ebN_local_kb*nSpace;
680  register double u_ext=0.0,
681  grad_u_ext[nSpace],
682  m_ext=0.0,
683  dm_ext=0.0,
684  f_ext[nSpace],
685  df_ext[nSpace],
686  flux_ext=0.0,
687  bc_u_ext=0.0,
688  //bc_grad_u_ext[nSpace],
689  bc_m_ext=0.0,
690  bc_dm_ext=0.0,
691  bc_f_ext[nSpace],
692  bc_df_ext[nSpace],
693  jac_ext[nSpace*nSpace],
694  jacDet_ext,
695  jacInv_ext[nSpace*nSpace],
696  boundaryJac[nSpace*(nSpace-1)],
697  metricTensor[(nSpace-1)*(nSpace-1)],
698  metricTensorDetSqrt,
699  dS,
700  u_test_dS[nDOF_test_element],
701  u_grad_trial_trace[nDOF_trial_element*nSpace],
702  normal[nSpace],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
703  //VRANS
704  porosity_ext,
705  //
706  G[nSpace*nSpace],G_dd_G,tr_G;
707  //
708  //calculate the solution and gradients at quadrature points
709  //
710  //compute information about mapping from reference element to physical element
711  ck.calculateMapping_elementBoundary(eN,
712  ebN_local,
713  kb,
714  ebN_local_kb,
715  mesh_dof.data(),
716  mesh_l2g.data(),
717  mesh_trial_trace_ref.data(),
718  mesh_grad_trial_trace_ref.data(),
719  boundaryJac_ref.data(),
720  jac_ext,
721  jacDet_ext,
722  jacInv_ext,
723  boundaryJac,
724  metricTensor,
725  metricTensorDetSqrt,
726  normal_ref.data(),
727  normal,
728  x_ext,y_ext,z_ext);
729  ck.calculateMappingVelocity_elementBoundary(eN,
730  ebN_local,
731  kb,
732  ebN_local_kb,
733  mesh_velocity_dof.data(),
734  mesh_l2g.data(),
735  mesh_trial_trace_ref.data(),
736  xt_ext,yt_ext,zt_ext,
737  normal,
738  boundaryJac,
739  metricTensor,
740  integralScaling);
741  //std::cout<<"metricTensorDetSqrt "<<metricTensorDetSqrt<<" integralScaling "<<integralScaling<<std::endl;
742  dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref.data()[kb];
743  //get the metric tensor
744  //cek todo use symmetry
745  ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
746  //compute shape and solution information
747  //shape
748  ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],
749  jacInv_ext,
750  u_grad_trial_trace);
751  //solution and gradients
752  if (STABILIZATION_TYPE==1) //explicit
753  {
754  ck.valFromDOF(u_dof_old.data(),
755  &u_l2g.data()[eN_nDOF_trial_element],
756  &u_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],
757  u_ext);
758  ck.gradFromDOF(u_dof_old.data(),
759  &u_l2g.data()[eN_nDOF_trial_element],
760  u_grad_trial_trace,
761  grad_u_ext);
762  }
763  else
764  {
765  ck.valFromDOF(u_dof.data(),
766  &u_l2g.data()[eN_nDOF_trial_element],
767  &u_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],
768  u_ext);
769  ck.gradFromDOF(u_dof.data(),
770  &u_l2g.data()[eN_nDOF_trial_element],
771  u_grad_trial_trace,
772  grad_u_ext);
773  }
774  //precalculate test function products with integration weights
775  for (int j=0;j<nDOF_trial_element;j++)
776  {
777  u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
778  }
779  //
780  //load the boundary values
781  //
782  bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+(1-isDOFBoundary_u.data()[ebNE_kb])*u_ext;
783  //VRANS
784  porosity_ext = ebqe_porosity_ext.data()[ebNE_kb];
785  //
786  //
787  //calculate the pde coefficients using the solution and the boundary values for the solution
788  //
789  evaluateCoefficients(&ebqe_velocity_ext.data()[ebNE_kb_nSpace],
790  u_ext,
791  //VRANS
792  porosity_ext,
793  //
794  m_ext,
795  dm_ext,
796  f_ext,
797  df_ext);
798  evaluateCoefficients(&ebqe_velocity_ext.data()[ebNE_kb_nSpace],
799  bc_u_ext,
800  //VRANS
801  porosity_ext,
802  //
803  bc_m_ext,
804  bc_dm_ext,
805  bc_f_ext,
806  bc_df_ext);
807  //
808  //moving mesh
809  //
810  double mesh_velocity[3];
811  mesh_velocity[0] = xt_ext;
812  mesh_velocity[1] = yt_ext;
813  mesh_velocity[2] = zt_ext;
814  //std::cout<<"mesh_velocity ext"<<std::endl;
815  for (int I=0;I<nSpace;I++)
816  {
817  //std::cout<<mesh_velocity[I]<<std::endl;
818  f_ext[I] -= MOVING_DOMAIN*m_ext*mesh_velocity[I];
819  df_ext[I] -= MOVING_DOMAIN*dm_ext*mesh_velocity[I];
820  bc_f_ext[I] -= MOVING_DOMAIN*bc_m_ext*mesh_velocity[I];
821  bc_df_ext[I] -= MOVING_DOMAIN*bc_dm_ext*mesh_velocity[I];
822  }
823  //
824  //calculate the numerical fluxes
825  //
826  exteriorNumericalAdvectiveFlux(isDOFBoundary_u.data()[ebNE_kb],
827  isFluxBoundary_u.data()[ebNE_kb],
828  normal,
829  bc_u_ext,
830  ebqe_bc_flux_u_ext.data()[ebNE_kb],
831  u_ext,
832  df_ext,//VRANS includes porosity
833  flux_ext);
834  ebqe_flux.data()[ebNE_kb] = flux_ext;
835  //save for other models? cek need to be consistent with numerical flux
836  if(flux_ext >=0.0)
837  ebqe_u.data()[ebNE_kb] = u_ext;
838  else
839  ebqe_u.data()[ebNE_kb] = bc_u_ext;
840 
841  if (STABILIZATION_TYPE==1)
842  if (stage==1)
843  flux_ext *= 1./3*dt;
844  else
845  flux_ext *= dt;
846 
847  //
848  //update residuals
849  //
850  for (int i=0;i<nDOF_test_element;i++)
851  {
852  //int ebNE_kb_i = ebNE_kb*nDOF_test_element+i;
853  elementResidual_u[i] += ck.ExteriorElementBoundaryFlux(flux_ext,u_test_dS[i]);
854  }//i
855  }//kb
856  //
857  //update the element and global residual storage
858  //
859  for (int i=0;i<nDOF_test_element;i++)
860  {
861  int eN_i = eN*nDOF_test_element+i;
862  globalResidual.data()[offset_u+stride_u*r_l2g.data()[eN_i]] += elementResidual_u[i];
863  }//i
864  }//ebNE
865  if (STABILIZATION_TYPE==1)
866  {
867  meanEntropy /= meanOmega;
868  double norm_factor = fmax(fabs(maxEntropy - meanEntropy), fabs(meanEntropy-minEntropy));
869  for(int eN=0;eN<nElements_global;eN++)
870  {
871  double hK=elementDiameter.data()[eN]/degree_polynomial;
872  double linear_viscosity = cMax*hK*maxVel[eN];
873  double entropy_viscosity = cE*hK*hK*maxEntRes[eN]/norm_factor;
874  for (int k=0;k<nQuadraturePoints_element;k++)
875  {
876  register int eN_k = eN*nQuadraturePoints_element+k;
877  q_numDiff_u.data()[eN_k] = fmin(linear_viscosity,entropy_viscosity);
878  }
879  }
880  }
881 
882  }
883 
884  void calculateJacobian(arguments_dict& args)
885  {
886  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
887  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
888  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
889  xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
890  double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
891  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
892  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
893  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
894  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
895  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
896  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
897  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
898  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
899  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
900  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
901  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
902  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
903  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
904  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
905  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
906  int nElements_global = args.scalar<int>("nElements_global");
907  double useMetrics = args.scalar<double>("useMetrics");
908  double alphaBDF = args.scalar<double>("alphaBDF");
909  int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
910  double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
911  const xt::pyarray<double>& q_porosity = args.array<double>("q_porosity");
912  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
913  xt::pyarray<int>& r_l2g = args.array<int>("r_l2g");
914  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
915  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
916  xt::pyarray<double>& velocity = args.array<double>("velocity");
917  xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
918  xt::pyarray<double>& cfl = args.array<double>("cfl");
919  xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
920  xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
921  xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
922  xt::pyarray<double>& globalJacobian = args.array<double>("globalJacobian");
923  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
924  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
925  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
926  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
927  xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
928  const xt::pyarray<double>& ebqe_porosity_ext = args.array<double>("ebqe_porosity_ext");
929  xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
930  xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
931  xt::pyarray<int>& isFluxBoundary_u = args.array<int>("isFluxBoundary_u");
932  xt::pyarray<double>& ebqe_bc_flux_u_ext = args.array<double>("ebqe_bc_flux_u_ext");
933  xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
934  int STABILIZATION_TYPE = args.scalar<int>("STABILIZATION_TYPE");
935  //std::cout<<"ndjaco address "<<q_numDiff_u_last.data()<<std::endl;
936  double Ct_sge = 4.0;
937  //
938  //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
939  //
940  for(int eN=0;eN<nElements_global;eN++)
941  {
942  register double elementJacobian_u_u[nDOF_test_element][nDOF_trial_element];
943  for (int i=0;i<nDOF_test_element;i++)
944  for (int j=0;j<nDOF_trial_element;j++)
945  {
946  elementJacobian_u_u[i][j]=0.0;
947  }
948  for (int k=0;k<nQuadraturePoints_element;k++)
949  {
950  int eN_k = eN*nQuadraturePoints_element+k, //index to a scalar at a quadrature point
951  eN_k_nSpace = eN_k*nSpace,
952  eN_nDOF_trial_element = eN*nDOF_trial_element; //index to a vector at a quadrature point
953 
954  //declare local storage
955  register double u=0.0,
956  grad_u[nSpace],
957  m=0.0,dm=0.0,
958  f[nSpace],df[nSpace],
959  m_t=0.0,dm_t=0.0,
960  dpdeResidual_u_u[nDOF_trial_element],
961  Lstar_u[nDOF_test_element],
962  dsubgridError_u_u[nDOF_trial_element],
963  tau=0.0,tau0=0.0,tau1=0.0,
964  jac[nSpace*nSpace],
965  jacDet,
966  jacInv[nSpace*nSpace],
967  u_grad_trial[nDOF_trial_element*nSpace],
968  dV,
969  u_test_dV[nDOF_test_element],
970  u_grad_test_dV[nDOF_test_element*nSpace],
971  x,y,z,xt,yt,zt,
972  //VRANS
973  porosity,
974  //
975  G[nSpace*nSpace],G_dd_G,tr_G;
976  //
977  //calculate solution and gradients at quadrature points
978  //
979  // u=0.0;
980  // for (int I=0;I<nSpace;I++)
981  // {
982  // grad_u[I]=0.0;
983  // }
984  // for (int j=0;j<nDOF_trial_element;j++)
985  // {
986  // int eN_j=eN*nDOF_trial_element+j;
987  // int eN_k_j=eN_k*nDOF_trial_element+j;
988  // int eN_k_j_nSpace = eN_k_j*nSpace;
989 
990  // u += valFromDOF_c(u_dof.data()[u_l2g.data()[eN_j]],u_trial[eN_k_j]);
991  // for (int I=0;I<nSpace;I++)
992  // {
993  // grad_u[I] += gradFromDOF_c(u_dof.data()[u_l2g.data()[eN_j]],u_grad_trial[eN_k_j_nSpace+I]);
994  // }
995  // }
996  //get jacobian, etc for mapping reference element
997  ck.calculateMapping_element(eN,
998  k,
999  mesh_dof.data(),
1000  mesh_l2g.data(),
1001  mesh_trial_ref.data(),
1002  mesh_grad_trial_ref.data(),
1003  jac,
1004  jacDet,
1005  jacInv,
1006  x,y,z);
1007  ck.calculateMappingVelocity_element(eN,
1008  k,
1009  mesh_velocity_dof.data(),
1010  mesh_l2g.data(),
1011  mesh_trial_ref.data(),
1012  xt,yt,zt);
1013  //get the physical integration weight
1014  dV = fabs(jacDet)*dV_ref.data()[k];
1015  ck.calculateG(jacInv,G,G_dd_G,tr_G);
1016  //get the trial function gradients
1017  ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
1018  //get the solution
1019  ck.valFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],u);
1020  //get the solution gradients
1021  ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial,grad_u);
1022  //precalculate test function products with integration weights
1023  for (int j=0;j<nDOF_trial_element;j++)
1024  {
1025  u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
1026  for (int I=0;I<nSpace;I++)
1027  {
1028  u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
1029  }
1030  }
1031  //VRANS
1032  porosity = q_porosity.data()[eN_k];
1033  //
1034  //
1035  //calculate pde coefficients and derivatives at quadrature points
1036  //
1037  evaluateCoefficients(&velocity.data()[eN_k_nSpace],
1038  u,
1039  //VRANS
1040  porosity,
1041  //
1042  m,
1043  dm,
1044  f,
1045  df);
1046  //
1047  //moving mesh
1048  //
1049  double mesh_velocity[3];
1050  mesh_velocity[0] = xt;
1051  mesh_velocity[1] = yt;
1052  mesh_velocity[2] = zt;
1053  //std::cout<<"qj mesh_velocity"<<std::endl;
1054  for(int I=0;I<nSpace;I++)
1055  {
1056  //std::cout<<mesh_velocity[I]<<std::endl;
1057  f[I] -= MOVING_DOMAIN*m*mesh_velocity[I];
1058  df[I] -= MOVING_DOMAIN*dm*mesh_velocity[I];
1059  }
1060  //
1061  //calculate time derivatives
1062  //
1063  ck.bdf(alphaBDF,
1064  q_m_betaBDF.data()[eN_k],//since m_t isn't used, we don't have to correct mass
1065  m,
1066  dm,
1067  m_t,
1068  dm_t);
1069  //
1070  //calculate subgrid error contribution to the Jacobian (strong residual, adjoint, jacobian of strong residual)
1071  //
1072  //calculate the adjoint times the test functions
1073  for (int i=0;i<nDOF_test_element;i++)
1074  {
1075  // int eN_k_i_nSpace = (eN_k*nDOF_trial_element+i)*nSpace;
1076  // Lstar_u[i]=ck.Advection_adjoint(df,&u_grad_test_dV[eN_k_i_nSpace]);
1077  register int i_nSpace = i*nSpace;
1078  Lstar_u[i]=ck.Advection_adjoint(df,&u_grad_test_dV[i_nSpace]);
1079  }
1080  //calculate the Jacobian of strong residual
1081  for (int j=0;j<nDOF_trial_element;j++)
1082  {
1083  //int eN_k_j=eN_k*nDOF_trial_element+j;
1084  //int eN_k_j_nSpace = eN_k_j*nSpace;
1085  int j_nSpace = j*nSpace;
1086  dpdeResidual_u_u[j]= ck.MassJacobian_strong(dm_t,u_trial_ref.data()[k*nDOF_trial_element+j]) +
1087  ck.AdvectionJacobian_strong(df,&u_grad_trial[j_nSpace]);
1088  }
1089  //tau and tau*Res
1090  calculateSubgridError_tau(elementDiameter.data()[eN],
1091  dm_t,
1092  df,
1093  cfl.data()[eN_k],
1094  tau0);
1095 
1096  calculateSubgridError_tau(Ct_sge,
1097  G,
1098  dm_t,
1099  df,
1100  tau1,
1101  cfl.data()[eN_k]);
1102  tau = useMetrics*tau1+(1.0-useMetrics)*tau0;
1103 
1104  for(int j=0;j<nDOF_trial_element;j++)
1105  dsubgridError_u_u[j] = -tau*dpdeResidual_u_u[j];
1106 
1107  for(int i=0;i<nDOF_test_element;i++)
1108  {
1109  //int eN_k_i=eN_k*nDOF_test_element+i;
1110  //int eN_k_i_nSpace=eN_k_i*nSpace;
1111  for(int j=0;j<nDOF_trial_element;j++)
1112  {
1113  //int eN_k_j=eN_k*nDOF_trial_element+j;
1114  //int eN_k_j_nSpace = eN_k_j*nSpace;
1115  int j_nSpace = j*nSpace;
1116  int i_nSpace = i*nSpace;
1117  if (STABILIZATION_TYPE==0)
1118  {
1119  elementJacobian_u_u[i][j] +=
1120  ck.MassJacobian_weak(dm_t,
1121  u_trial_ref.data()[k*nDOF_trial_element+j],
1122  u_test_dV[i]) +
1123  ck.AdvectionJacobian_weak(df,
1124  u_trial_ref.data()[k*nDOF_trial_element+j],
1125  &u_grad_test_dV[i_nSpace]) +
1126  ck.SubgridErrorJacobian(dsubgridError_u_u[j],Lstar_u[i]) +
1127  ck.NumericalDiffusionJacobian(q_numDiff_u_last.data()[eN_k],
1128  &u_grad_trial[j_nSpace],
1129  &u_grad_test_dV[i_nSpace]); //implicit
1130  }
1131  else
1132  {
1133  elementJacobian_u_u[i][j] +=
1134  ck.MassJacobian_weak(1.0,
1135  u_trial_ref.data()[k*nDOF_trial_element+j],
1136  u_test_dV[i]);
1137  }
1138  }//j
1139  }//i
1140  }//k
1141  //
1142  //load into element Jacobian into global Jacobian
1143  //
1144  for (int i=0;i<nDOF_test_element;i++)
1145  {
1146  int eN_i = eN*nDOF_test_element+i;
1147  for (int j=0;j<nDOF_trial_element;j++)
1148  {
1149  int eN_i_j = eN_i*nDOF_trial_element+j;
1150  globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_u_u.data()[eN_i_j]] += elementJacobian_u_u[i][j];
1151  }//j
1152  }//i
1153  }//elements
1154  //
1155  //loop over exterior element boundaries to compute the surface integrals and load them into the global Jacobian
1156  //
1157  if (STABILIZATION_TYPE==0)
1158  for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
1159  {
1160  register int ebN = exteriorElementBoundariesArray.data()[ebNE];
1161  register int eN = elementBoundaryElementsArray.data()[ebN*2+0],
1162  ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
1163  eN_nDOF_trial_element = eN*nDOF_trial_element;
1164  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
1165  {
1166  register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
1167  ebNE_kb_nSpace = ebNE_kb*nSpace,
1168  ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
1169  ebN_local_kb_nSpace = ebN_local_kb*nSpace;
1170  register double u_ext=0.0,
1171  grad_u_ext[nSpace],
1172  m_ext=0.0,
1173  dm_ext=0.0,
1174  f_ext[nSpace],
1175  df_ext[nSpace],
1176  dflux_u_u_ext=0.0,
1177  bc_u_ext=0.0,
1178  //bc_grad_u_ext[nSpace],
1179  bc_m_ext=0.0,
1180  bc_dm_ext=0.0,
1181  bc_f_ext[nSpace],
1182  bc_df_ext[nSpace],
1183  fluxJacobian_u_u[nDOF_trial_element],
1184  jac_ext[nSpace*nSpace],
1185  jacDet_ext,
1186  jacInv_ext[nSpace*nSpace],
1187  boundaryJac[nSpace*(nSpace-1)],
1188  metricTensor[(nSpace-1)*(nSpace-1)],
1189  metricTensorDetSqrt,
1190  dS,
1191  u_test_dS[nDOF_test_element],
1192  u_grad_trial_trace[nDOF_trial_element*nSpace],
1193  normal[nSpace],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
1194  //VRANS
1195  porosity_ext,
1196  //
1197  G[nSpace*nSpace],G_dd_G,tr_G;
1198  //
1199  //calculate the solution and gradients at quadrature points
1200  //
1201  // u_ext=0.0;
1202  // for (int I=0;I<nSpace;I++)
1203  // {
1204  // grad_u_ext[I] = 0.0;
1205  // bc_grad_u_ext[I] = 0.0;
1206  // }
1207  // for (int j=0;j<nDOF_trial_element;j++)
1208  // {
1209  // register int eN_j = eN*nDOF_trial_element+j,
1210  // ebNE_kb_j = ebNE_kb*nDOF_trial_element+j,
1211  // ebNE_kb_j_nSpace= ebNE_kb_j*nSpace;
1212  // u_ext += valFromDOF_c(u_dof.data()[u_l2g.data()[eN_j]],u_trial_ext[ebNE_kb_j]);
1213 
1214  // for (int I=0;I<nSpace;I++)
1215  // {
1216  // grad_u_ext[I] += gradFromDOF_c(u_dof.data()[u_l2g.data()[eN_j]],u_grad_trial_ext[ebNE_kb_j_nSpace+I]);
1217  // }
1218  // }
1219  ck.calculateMapping_elementBoundary(eN,
1220  ebN_local,
1221  kb,
1222  ebN_local_kb,
1223  mesh_dof.data(),
1224  mesh_l2g.data(),
1225  mesh_trial_trace_ref.data(),
1226  mesh_grad_trial_trace_ref.data(),
1227  boundaryJac_ref.data(),
1228  jac_ext,
1229  jacDet_ext,
1230  jacInv_ext,
1231  boundaryJac,
1232  metricTensor,
1233  metricTensorDetSqrt,
1234  normal_ref.data(),
1235  normal,
1236  x_ext,y_ext,z_ext);
1237  ck.calculateMappingVelocity_elementBoundary(eN,
1238  ebN_local,
1239  kb,
1240  ebN_local_kb,
1241  mesh_velocity_dof.data(),
1242  mesh_l2g.data(),
1243  mesh_trial_trace_ref.data(),
1244  xt_ext,yt_ext,zt_ext,
1245  normal,
1246  boundaryJac,
1247  metricTensor,
1248  integralScaling);
1249  //std::cout<<"J mtsqrdet "<<metricTensorDetSqrt<<" integralScaling "<<integralScaling<<std::endl;
1250  dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref.data()[kb];
1251  //dS = metricTensorDetSqrt*dS_ref.data()[kb];
1252  ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
1253  //compute shape and solution information
1254  //shape
1255  ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,u_grad_trial_trace);
1256  //solution and gradients
1257  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);
1258  ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial_trace,grad_u_ext);
1259  //precalculate test function products with integration weights
1260  for (int j=0;j<nDOF_trial_element;j++)
1261  {
1262  u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
1263  }
1264  //
1265  //load the boundary values
1266  //
1267  bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+(1-isDOFBoundary_u.data()[ebNE_kb])*u_ext;
1268  //VRANS
1269  porosity_ext = ebqe_porosity_ext.data()[ebNE_kb];
1270  //
1271  //
1272  //calculate the internal and external trace of the pde coefficients
1273  //
1274  evaluateCoefficients(&ebqe_velocity_ext.data()[ebNE_kb_nSpace],
1275  u_ext,
1276  //VRANS
1277  porosity_ext,
1278  //
1279  m_ext,
1280  dm_ext,
1281  f_ext,
1282  df_ext);
1283  evaluateCoefficients(&ebqe_velocity_ext.data()[ebNE_kb_nSpace],
1284  bc_u_ext,
1285  //VRANS
1286  porosity_ext,
1287  //
1288  bc_m_ext,
1289  bc_dm_ext,
1290  bc_f_ext,
1291  bc_df_ext);
1292  //
1293  //moving domain
1294  //
1295  double mesh_velocity[3];
1296  mesh_velocity[0] = xt_ext;
1297  mesh_velocity[1] = yt_ext;
1298  mesh_velocity[2] = zt_ext;
1299  //std::cout<<"ext J mesh_velocity"<<std::endl;
1300  for (int I=0;I<nSpace;I++)
1301  {
1302  //std::cout<<mesh_velocity[I]<<std::endl;
1303  f_ext[I] -= MOVING_DOMAIN*m_ext*mesh_velocity[I];
1304  df_ext[I] -= MOVING_DOMAIN*dm_ext*mesh_velocity[I];
1305  bc_f_ext[I] -= MOVING_DOMAIN*bc_m_ext*mesh_velocity[I];
1306  bc_df_ext[I] -= MOVING_DOMAIN*bc_dm_ext*mesh_velocity[I];
1307  }
1308  //
1309  //calculate the numerical fluxes
1310  //
1311  exteriorNumericalAdvectiveFluxDerivative(isDOFBoundary_u.data()[ebNE_kb],
1312  isFluxBoundary_u.data()[ebNE_kb],
1313  normal,
1314  df_ext,//VRANS holds porosity
1315  dflux_u_u_ext);
1316  //
1317  //calculate the flux jacobian
1318  //
1319  for (int j=0;j<nDOF_trial_element;j++)
1320  {
1321  //register int ebNE_kb_j = ebNE_kb*nDOF_trial_element+j;
1322  register int ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j;
1323  fluxJacobian_u_u[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_u_u_ext,u_trial_trace_ref.data()[ebN_local_kb_j]);
1324  }//j
1325  //
1326  //update the global Jacobian from the flux Jacobian
1327  //
1328  for (int i=0;i<nDOF_test_element;i++)
1329  {
1330  register int eN_i = eN*nDOF_test_element+i;
1331  //register int ebNE_kb_i = ebNE_kb*nDOF_test_element+i;
1332  for (int j=0;j<nDOF_trial_element;j++)
1333  {
1334  register int ebN_i_j = ebN*4*nDOF_test_X_trial_element + i*nDOF_trial_element + j;
1335  globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j]] += fluxJacobian_u_u[j]*u_test_dS[i];
1336  }//j
1337  }//i
1338  }//kb
1339  }//ebNE
1340  }//computeJacobian
1341 
1342  void FCTStep(arguments_dict& args)
1343  {
1344  double dt = args.scalar<double>("dt");
1345  int NNZ = args.scalar<int>("NNZ");
1346  int numDOFs = args.scalar<int>("numDOFs");
1347  xt::pyarray<double>& lumped_mass_matrix = args.array<double>("lumped_mass_matrix");
1348  xt::pyarray<double>& soln = args.array<double>("soln");
1349  xt::pyarray<double>& solH = args.array<double>("solH");
1350  xt::pyarray<double>& uLow = args.array<double>("uLow");
1351  xt::pyarray<double>& dLow = args.array<double>("dLow");
1352  xt::pyarray<double>& limited_solution = args.array<double>("limited_solution");
1353  xt::pyarray<int>& csrRowIndeces_DofLoops = args.array<int>("csrRowIndeces_DofLoops");
1354  xt::pyarray<int>& csrColumnOffsets_DofLoops = args.array<int>("csrColumnOffsets_DofLoops");
1355  xt::pyarray<double>& MassMatrix = args.array<double>("MassMatrix");
1356  xt::pyarray<double>& dt_times_dH_minus_dL = args.array<double>("dt_times_dH_minus_dL");
1357  xt::pyarray<double>& min_u_bc = args.array<double>("min_u_bc");
1358  xt::pyarray<double>& max_u_bc = args.array<double>("max_u_bc");
1359  int LUMPED_MASS_MATRIX = args.scalar<int>("LUMPED_MASS_MATRIX");
1360  int STABILIZATION_TYPE = args.scalar<int>("STABILIZATION_TYPE");
1361  Rpos.resize(numDOFs,0.0);
1362  Rneg.resize(numDOFs,0.0);
1363  FluxCorrectionMatrix.resize(NNZ,0.0);
1365  // LOOP in DOFs //
1367  int ij=0;
1368  for (int i=0; i<numDOFs; i++)
1369  {
1370  //read some vectors
1371  double solHi = solH.data()[i];
1372  double solni = soln.data()[i];
1373  double mi = lumped_mass_matrix.data()[i];
1374  double uLowi = uLow.data()[i];
1375  double uDotLowi = (uLowi - solni)/dt;
1376  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)
1377  if (GLOBAL_FCT==1)
1378  {
1379  mini = 0.;
1380  maxi = 1.;
1381  }
1382 
1383  double Pposi=0, Pnegi=0;
1384  // LOOP OVER THE SPARSITY PATTERN (j-LOOP)//
1385  for (int offset=csrRowIndeces_DofLoops.data()[i]; offset<csrRowIndeces_DofLoops.data()[i+1]; offset++)
1386  {
1387  int j = csrColumnOffsets_DofLoops.data()[offset];
1388  double solnj = soln.data()[j];
1390  // COMPUTE THE BOUNDS //
1392  if (GLOBAL_FCT == 0)
1393  {
1394  mini = fmin(mini,solnj);
1395  maxi = fmax(maxi,solnj);
1396  }
1397  double uLowj = uLow.data()[j];
1398  double uDotLowj = (uLowj - solnj)/dt;
1399  // i-th row of flux correction matrix
1400  if (STABILIZATION_TYPE==4) // DK high-order, linearly stable anti-dif. flux
1401  {
1402  FluxCorrectionMatrix[ij] = dt*(MassMatrix.data()[ij]*(uDotLowi-uDotLowj)
1403  + dLow.data()[ij]*(uLowi-uLowj));
1404  }
1405  else
1406  {
1407  double ML_minus_MC =
1408  (LUMPED_MASS_MATRIX == 1 ? 0. : (i==j ? 1. : 0.)*mi - MassMatrix.data()[ij]);
1409  FluxCorrectionMatrix[ij] = ML_minus_MC * (solH.data()[j]-solnj - (solHi-solni))
1410  + dt_times_dH_minus_dL.data()[ij]*(solnj-solni);
1411  }
1412 
1414  // COMPUTE P VECTORS //
1416  Pposi += FluxCorrectionMatrix[ij]*((FluxCorrectionMatrix[ij] > 0) ? 1. : 0.);
1417  Pnegi += FluxCorrectionMatrix[ij]*((FluxCorrectionMatrix[ij] < 0) ? 1. : 0.);
1418 
1419  //update ij
1420  ij+=1;
1421  }
1423  // COMPUTE Q VECTORS //
1425  double Qposi = mi*(maxi-uLow.data()[i]);
1426  double Qnegi = mi*(mini-uLow.data()[i]);
1427 
1429  // 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  } // i DOFs
1434 
1436  // COMPUTE LIMITERS //
1438  ij=0;
1439  for (int i=0; i<numDOFs; i++)
1440  {
1441  double ith_Limiter_times_FluxCorrectionMatrix = 0.;
1442  double Rposi = Rpos[i], Rnegi = Rneg[i];
1443  // LOOP OVER THE SPARSITY PATTERN (j-LOOP)//
1444  for (int offset=csrRowIndeces_DofLoops.data()[i]; offset<csrRowIndeces_DofLoops.data()[i+1]; offset++)
1445  {
1446  int j = csrColumnOffsets_DofLoops.data()[offset];
1447  double Lij = 1;
1448  Lij = ((FluxCorrectionMatrix[ij]>0) ? fmin(Rposi,Rneg[j]) : fmin(Rnegi,Rpos[j]));
1449  ith_Limiter_times_FluxCorrectionMatrix += Lij * FluxCorrectionMatrix[ij];
1450  //update ij
1451  ij+=1;
1452  }
1453  limited_solution.data()[i] = uLow.data()[i] + 1./lumped_mass_matrix.data()[i]*ith_Limiter_times_FluxCorrectionMatrix;
1454  }
1455  }
1456 
1457  void calculateResidualEdgeBased(arguments_dict& args)
1458  {
1459  double dt = args.scalar<double>("dt");
1460  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1461  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1462  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1463  xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
1464  double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
1465  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1466  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1467  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1468  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1469  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1470  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
1471  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
1472  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
1473  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
1474  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
1475  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
1476  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
1477  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
1478  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
1479  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
1480  int nElements_global = args.scalar<int>("nElements_global");
1481  double useMetrics = args.scalar<double>("useMetrics");
1482  double alphaBDF = args.scalar<double>("alphaBDF");
1483  int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
1484  double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
1485  double sc_uref = args.scalar<double>("sc_uref");
1486  double sc_alpha = args.scalar<double>("sc_alpha");
1487  const xt::pyarray<double>& q_porosity = args.array<double>("q_porosity");
1488  const xt::pyarray<double>& porosity_dof = args.array<double>("porosity_dof");
1489  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1490  xt::pyarray<int>& r_l2g = args.array<int>("r_l2g");
1491  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1492  double degree_polynomial = args.scalar<double>("degree_polynomial");
1493  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1494  xt::pyarray<double>& u_dof_old = args.array<double>("u_dof_old");
1495  xt::pyarray<double>& velocity = args.array<double>("velocity");
1496  xt::pyarray<double>& q_m = args.array<double>("q_m");
1497  xt::pyarray<double>& q_u = args.array<double>("q_u");
1498  xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
1499  xt::pyarray<double>& q_dV = args.array<double>("q_dV");
1500  xt::pyarray<double>& q_dV_last = args.array<double>("q_dV_last");
1501  xt::pyarray<double>& cfl = args.array<double>("cfl");
1502  xt::pyarray<double>& edge_based_cfl = args.array<double>("edge_based_cfl");
1503  xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
1504  xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
1505  int offset_u = args.scalar<int>("offset_u");
1506  int stride_u = args.scalar<int>("stride_u");
1507  xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
1508  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
1509  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
1510  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
1511  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
1512  xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
1513  const xt::pyarray<double>& ebqe_porosity_ext = args.array<double>("ebqe_porosity_ext");
1514  xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
1515  xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
1516  xt::pyarray<int>& isFluxBoundary_u = args.array<int>("isFluxBoundary_u");
1517  xt::pyarray<double>& ebqe_bc_flux_u_ext = args.array<double>("ebqe_bc_flux_u_ext");
1518  xt::pyarray<double>& ebqe_phi = args.array<double>("ebqe_phi");
1519  double epsFact = args.scalar<double>("epsFact");
1520  xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
1521  xt::pyarray<double>& ebqe_flux = args.array<double>("ebqe_flux");
1522  int stage = args.scalar<int>("stage");
1523  xt::pyarray<double>& uTilde_dof = args.array<double>("uTilde_dof");
1524  double cE = args.scalar<double>("cE");
1525  double cMax = args.scalar<double>("cMax");
1526  double cK = args.scalar<double>("cK");
1527  double uL = args.scalar<double>("uL");
1528  double uR = args.scalar<double>("uR");
1529  int numDOFs = args.scalar<int>("numDOFs");
1530  int NNZ = args.scalar<int>("NNZ");
1531  xt::pyarray<int>& csrRowIndeces_DofLoops = args.array<int>("csrRowIndeces_DofLoops");
1532  xt::pyarray<int>& csrColumnOffsets_DofLoops = args.array<int>("csrColumnOffsets_DofLoops");
1533  xt::pyarray<int>& csrRowIndeces_CellLoops = args.array<int>("csrRowIndeces_CellLoops");
1534  xt::pyarray<int>& csrColumnOffsets_CellLoops = args.array<int>("csrColumnOffsets_CellLoops");
1535  xt::pyarray<int>& csrColumnOffsets_eb_CellLoops = args.array<int>("csrColumnOffsets_eb_CellLoops");
1536  xt::pyarray<double>& ML = args.array<double>("ML");
1537  int LUMPED_MASS_MATRIX = args.scalar<int>("LUMPED_MASS_MATRIX");
1538  int STABILIZATION_TYPE = args.scalar<int>("STABILIZATION_TYPE");
1539  int ENTROPY_TYPE = args.scalar<int>("ENTROPY_TYPE");
1540  xt::pyarray<double>& uLow = args.array<double>("uLow");
1541  xt::pyarray<double>& dLow = args.array<double>("dLow");
1542  xt::pyarray<double>& dt_times_dH_minus_dL = args.array<double>("dt_times_dH_minus_dL");
1543  xt::pyarray<double>& min_u_bc = args.array<double>("min_u_bc");
1544  xt::pyarray<double>& max_u_bc = args.array<double>("max_u_bc");
1545  xt::pyarray<double>& quantDOFs = args.array<double>("quantDOFs");
1546  // NOTE: This function follows a different (but equivalent) implementation of the smoothness based indicator than NCLS.h
1547  // Allocate space for the transport matrices
1548  // This is used for first order KUZMIN'S METHOD
1549  TransportMatrix.resize(NNZ,0.0);
1550  TransposeTransportMatrix.resize(NNZ,0.0);
1551  // compute entropy and init global_entropy_residual and boundary_integral
1552  psi.resize(numDOFs,0.0);
1553  eta.resize(numDOFs,0.0);
1554  global_entropy_residual.resize(numDOFs,0.0);
1555  boundary_integral.resize(numDOFs,0.0);
1556 
1557  for (int i=0; i<numDOFs; i++)
1558  {
1559  // NODAL ENTROPY //
1560  if (STABILIZATION_TYPE==2) //EV stab
1561  {
1562  double porosity_times_solni = porosity_dof.data()[i]*u_dof_old.data()[i];
1563  eta[i] = ENTROPY_TYPE == 0 ? ENTROPY(porosity_times_solni,uL,uR) : ENTROPY_LOG(porosity_times_solni,uL,uR);
1564  global_entropy_residual[i]=0.;
1565  }
1566  boundary_integral[i]=0.;
1567  }
1568 
1570  // ** LOOP IN CELLS FOR CELL BASED TERMS ** //
1572  // HERE WE COMPUTE:
1573  // * Time derivative term. porosity*u_t
1574  // * cell based CFL (for reference)
1575  // * Entropy residual
1576  // * Transport matrices
1577  for(int eN=0;eN<nElements_global;eN++)
1578  {
1579  //declare local storage for local contributions and initialize
1580  register double
1581  elementResidual_u[nDOF_test_element],
1582  element_entropy_residual[nDOF_test_element];
1583  register double elementTransport[nDOF_test_element][nDOF_trial_element];
1584  register double elementTransposeTransport[nDOF_test_element][nDOF_trial_element];
1585  for (int i=0;i<nDOF_test_element;i++)
1586  {
1587  elementResidual_u[i]=0.0;
1588  element_entropy_residual[i]=0.0;
1589  for (int j=0;j<nDOF_trial_element;j++)
1590  {
1591  elementTransport[i][j]=0.0;
1592  elementTransposeTransport[i][j]=0.0;
1593  }
1594  }
1595  //loop over quadrature points and compute integrands
1596  for (int k=0;k<nQuadraturePoints_element;k++)
1597  {
1598  //compute indeces and declare local storage
1599  register int eN_k = eN*nQuadraturePoints_element+k,
1600  eN_k_nSpace = eN_k*nSpace,
1601  eN_nDOF_trial_element = eN*nDOF_trial_element;
1602  register double
1603  // for entropy residual
1604  aux_entropy_residual=0., DENTROPY_un, DENTROPY_uni,
1605  //for mass matrix contributions
1606  u=0.0, un=0.0, grad_un[nSpace], porosity_times_velocity[nSpace],
1607  u_test_dV[nDOF_trial_element],
1608  u_grad_trial[nDOF_trial_element*nSpace],
1609  u_grad_test_dV[nDOF_test_element*nSpace],
1610  //for general use
1611  jac[nSpace*nSpace], jacDet, jacInv[nSpace*nSpace],
1612  dV,x,y,z,xt,yt,zt,
1613  //VRANS
1614  porosity;
1615  //get the physical integration weight
1616  ck.calculateMapping_element(eN,
1617  k,
1618  mesh_dof.data(),
1619  mesh_l2g.data(),
1620  mesh_trial_ref.data(),
1621  mesh_grad_trial_ref.data(),
1622  jac,
1623  jacDet,
1624  jacInv,
1625  x,y,z);
1626  ck.calculateMappingVelocity_element(eN,
1627  k,
1628  mesh_velocity_dof.data(),
1629  mesh_l2g.data(),
1630  mesh_trial_ref.data(),
1631  xt,yt,zt);
1632  dV = fabs(jacDet)*dV_ref.data()[k];
1633  //get the solution (of Newton's solver). To compute time derivative term
1634  ck.valFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],u);
1635  //get the solution at quad point at tn and tnm1 for entropy viscosity
1636  ck.valFromDOF(u_dof_old.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],un);
1637  //get the solution gradients at tn for entropy viscosity
1638  ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
1639  ck.gradFromDOF(u_dof_old.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial,grad_un);
1640 
1641  //precalculate test function products with integration weights for mass matrix terms
1642  for (int j=0;j<nDOF_trial_element;j++)
1643  {
1644  u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
1645  for (int I=0;I<nSpace;I++)
1646  u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
1647  }
1648 
1649  //calculate time derivative at quadrature points
1650  if (q_dV_last.data()[eN_k] <= -100)
1651  q_dV_last.data()[eN_k] = dV;
1652  q_dV.data()[eN_k] = dV;
1653  //VRANS
1654  porosity = q_porosity.data()[eN_k];
1655  //
1656  //moving mesh
1657  //
1658  double mesh_velocity[3];
1659  mesh_velocity[0] = xt;
1660  mesh_velocity[1] = yt;
1661  mesh_velocity[2] = zt;
1662  //relative velocity at tn
1663  for (int I=0;I<nSpace;I++)
1664  porosity_times_velocity[I] = porosity*(velocity.data()[eN_k_nSpace+I]-MOVING_DOMAIN*mesh_velocity[I]);
1665 
1667  // CALCULATE CELL BASED CFL //
1669  calculateCFL(elementDiameter.data()[eN]/degree_polynomial,porosity_times_velocity,cfl.data()[eN_k]);
1670 
1672  // CALCULATE ENTROPY RESIDUAL AT QUAD POINT //
1674  if (STABILIZATION_TYPE==2) // EV stab
1675  {
1676  for (int I=0;I<nSpace;I++)
1677  aux_entropy_residual += porosity_times_velocity[I]*grad_un[I];
1678  DENTROPY_un = ENTROPY_TYPE==0 ? DENTROPY(porosity*un,uL,uR) : DENTROPY_LOG(porosity*un,uL,uR);
1679  }
1681  // ith-LOOP //
1683  for(int i=0;i<nDOF_test_element;i++)
1684  {
1685  // VECTOR OF ENTROPY RESIDUAL //
1686  int eN_i=eN*nDOF_test_element+i;
1687  if (STABILIZATION_TYPE==2) // EV stab
1688  {
1689  int gi = offset_u+stride_u*u_l2g.data()[eN_i]; //global i-th index
1690  double porosity_times_uni = porosity_dof.data()[gi]*u_dof_old.data()[gi];
1691  DENTROPY_uni = ENTROPY_TYPE == 0 ? DENTROPY(porosity_times_uni,uL,uR) : DENTROPY_LOG(porosity_times_uni,uL,uR);
1692  element_entropy_residual[i] += (DENTROPY_un - DENTROPY_uni)*aux_entropy_residual*u_test_dV[i];
1693  }
1694  elementResidual_u[i] += porosity*(u-un)*u_test_dV[i];
1696  // j-th LOOP // To construct transport matrices
1698  for(int j=0;j<nDOF_trial_element;j++)
1699  {
1700  int j_nSpace = j*nSpace;
1701  int i_nSpace = i*nSpace;
1702  elementTransport[i][j] += // -int[(vel.grad_wi)*wj*dx]
1703  ck.AdvectionJacobian_weak(porosity_times_velocity,
1704  u_trial_ref.data()[k*nDOF_trial_element+j],&u_grad_test_dV[i_nSpace]);
1705  elementTransposeTransport[i][j] += // -int[(vel.grad_wj)*wi*dx]
1706  ck.AdvectionJacobian_weak(porosity_times_velocity,
1707  u_trial_ref.data()[k*nDOF_trial_element+i],&u_grad_test_dV[j_nSpace]);
1708  }
1709  }//i
1710  //save solution for other models
1711  q_u.data()[eN_k] = u;
1712  q_m.data()[eN_k] = porosity*u;
1713  }
1715  // DISTRIBUTE // load cell based element into global residual
1717  for(int i=0;i<nDOF_test_element;i++)
1718  {
1719  int eN_i=eN*nDOF_test_element+i;
1720  int gi = offset_u+stride_u*u_l2g.data()[eN_i]; //global i-th index
1721 
1722  // distribute global residual for (lumped) mass matrix
1723  globalResidual.data()[gi] += elementResidual_u[i];
1724  // distribute entropy_residual
1725  if (STABILIZATION_TYPE==2) // EV Stab
1726  global_entropy_residual[gi] += element_entropy_residual[i];
1727 
1728  // distribute transport matrices
1729  for (int j=0;j<nDOF_trial_element;j++)
1730  {
1731  int eN_i_j = eN_i*nDOF_trial_element+j;
1732  TransportMatrix[csrRowIndeces_CellLoops.data()[eN_i] +
1733  csrColumnOffsets_CellLoops.data()[eN_i_j]] += elementTransport[i][j];
1734  TransposeTransportMatrix[csrRowIndeces_CellLoops.data()[eN_i] +
1735  csrColumnOffsets_CellLoops.data()[eN_i_j]]
1736  += elementTransposeTransport[i][j];
1737  }//j
1738  }//i
1739  }//elements
1740 
1742  // ADD OUTFLOW BOUNDARY TERM TO TRANSPORT MATRICES AND COMPUTE INFLOW BOUNDARY INTEGRAL //
1744  // * Compute outflow boundary integral as a matrix; i.e., int_B[ (vel.normal)*wi*wj*dx]
1745  for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
1746  {
1747  double min_u_bc_local = 1E10, max_u_bc_local = -1E10;
1748  register int ebN = exteriorElementBoundariesArray.data()[ebNE];
1749  register int eN = elementBoundaryElementsArray.data()[ebN*2+0],
1750  ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
1751  eN_nDOF_trial_element = eN*nDOF_trial_element;
1752  register double elementResidual_u[nDOF_test_element];
1753  for (int i=0;i<nDOF_test_element;i++)
1754  elementResidual_u[i]=0.0;
1755  // loop on quad points
1756  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
1757  {
1758  register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
1759  ebNE_kb_nSpace = ebNE_kb*nSpace,
1760  ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
1761  ebN_local_kb_nSpace = ebN_local_kb*nSpace;
1762  register double
1763  u_ext=0.0, bc_u_ext=0.0,
1764  porosity_times_velocity[nSpace],
1765  flux_ext=0.0, dflux_ext=0.0,
1766  fluxTransport[nDOF_trial_element],
1767  jac_ext[nSpace*nSpace],
1768  jacDet_ext,
1769  jacInv_ext[nSpace*nSpace],
1770  boundaryJac[nSpace*(nSpace-1)],
1771  metricTensor[(nSpace-1)*(nSpace-1)],
1772  metricTensorDetSqrt,
1773  dS,
1774  u_test_dS[nDOF_test_element],
1775  normal[nSpace],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,porosity_ext;
1776  // calculate mappings
1777  ck.calculateMapping_elementBoundary(eN,
1778  ebN_local,
1779  kb,
1780  ebN_local_kb,
1781  mesh_dof.data(),
1782  mesh_l2g.data(),
1783  mesh_trial_trace_ref.data(),
1784  mesh_grad_trial_trace_ref.data(),
1785  boundaryJac_ref.data(),
1786  jac_ext,
1787  jacDet_ext,
1788  jacInv_ext,
1789  boundaryJac,
1790  metricTensor,
1791  metricTensorDetSqrt,
1792  normal_ref.data(),
1793  normal,
1794  x_ext,y_ext,z_ext);
1795  ck.calculateMappingVelocity_elementBoundary(eN,
1796  ebN_local,
1797  kb,
1798  ebN_local_kb,
1799  mesh_velocity_dof.data(),
1800  mesh_l2g.data(),
1801  mesh_trial_trace_ref.data(),
1802  xt_ext,yt_ext,zt_ext,
1803  normal,
1804  boundaryJac,
1805  metricTensor,
1806  integralScaling);
1807  dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref.data()[kb];
1808  //compute shape and solution information
1809  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);
1810  //precalculate test function products with integration weights
1811  for (int j=0;j<nDOF_trial_element;j++)
1812  u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
1813 
1814  //VRANS
1815  porosity_ext = ebqe_porosity_ext.data()[ebNE_kb];
1816  //
1817  //moving mesh
1818  //
1819  double mesh_velocity[3];
1820  mesh_velocity[0] = xt_ext;
1821  mesh_velocity[1] = yt_ext;
1822  mesh_velocity[2] = zt_ext;
1823  //std::cout<<"mesh_velocity ext"<<std::endl;
1824  for (int I=0;I<nSpace;I++)
1825  porosity_times_velocity[I] = porosity_ext*(ebqe_velocity_ext.data()[ebNE_kb_nSpace+I] - MOVING_DOMAIN*mesh_velocity[I]);
1826  //
1827  //calculate the fluxes
1828  //
1829  double flow = 0.;
1830  for (int I=0; I < nSpace; I++)
1831  flow += normal[I]*porosity_times_velocity[I];
1832 
1833  if (flow >= 0) //outflow. This is handled via the transport matrices. Then flux_ext=0 and dflux_ext!=0
1834  {
1835  dflux_ext = flow;
1836  flux_ext = 0;
1837  // save external u
1838  ebqe_u.data()[ebNE_kb] = u_ext;
1839  }
1840  else // inflow. This is handled via the boundary integral. Then flux_ext!=0 and dflux_ext=0
1841  {
1842  dflux_ext = 0;
1843  // save external u
1844  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;
1845  if (isDOFBoundary_u.data()[ebNE_kb] == 1)
1846  flux_ext = ebqe_bc_u_ext.data()[ebNE_kb]*flow;
1847  else if (isFluxBoundary_u.data()[ebNE_kb] == 1)
1848  flux_ext = ebqe_bc_flux_u_ext.data()[ebNE_kb];
1849  else
1850  {
1851  std::cout<<"warning: VOF open boundary with no external trace, setting to zero for inflow"<<std::endl;
1852  flux_ext = 0.0;
1853  }
1854  }
1855 
1856  for (int j=0;j<nDOF_trial_element;j++)
1857  {
1858  // elementResidual. This is to include the inflow boundary integral.
1859  // NOTE: here I assume that we use a Galerkin approach st nDOF_test_element = nDOF_trial_element
1860  elementResidual_u[j] += flux_ext*u_test_dS[j];
1861  register int ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j;
1862  fluxTransport[j] = dflux_ext*u_trial_trace_ref.data()[ebN_local_kb_j];
1863  }//j
1865  // DISTRIBUTE OUTFLOW BOUNDARY TO TRANSPORT MATRICES //
1867  for (int i=0;i<nDOF_test_element;i++)
1868  {
1869  register int eN_i = eN*nDOF_test_element+i;
1870  for (int j=0;j<nDOF_trial_element;j++)
1871  {
1872  register int ebN_i_j = ebN*4*nDOF_test_X_trial_element + i*nDOF_trial_element + j;
1873  TransportMatrix[csrRowIndeces_CellLoops.data()[eN_i] + csrColumnOffsets_eb_CellLoops.data()[ebN_i_j]]
1874  += fluxTransport[j]*u_test_dS[i];
1875  TransposeTransportMatrix[csrRowIndeces_CellLoops.data()[eN_i] + csrColumnOffsets_eb_CellLoops.data()[ebN_i_j]]
1876  += fluxTransport[i]*u_test_dS[j];
1877  }//j
1878  }//i
1879  // local min/max at boundary
1880  min_u_bc_local = fmin(ebqe_u.data()[ebNE_kb], min_u_bc_local);
1881  max_u_bc_local = fmax(ebqe_u.data()[ebNE_kb], max_u_bc_local);
1882  }//kb
1883  // global min/max at boundary
1884  for (int i=0;i<nDOF_test_element;i++)
1885  {
1886  int eN_i = eN*nDOF_test_element+i;
1887  int gi = offset_u+stride_u*u_l2g.data()[eN_i]; //global i-th index
1888  globalResidual.data()[gi] += dt*elementResidual_u[i];
1889  boundary_integral[gi] += elementResidual_u[i];
1890  min_u_bc[gi] = fmin(min_u_bc_local,min_u_bc[gi]);
1891  max_u_bc[gi] = fmax(max_u_bc_local,max_u_bc[gi]);
1892  }
1893  }//ebNE
1894  // END OF ADDING BOUNDARY TERM TO TRANSPORT MATRICES and COMPUTING BOUNDARY INTEGRAL //
1895 
1897  // COMPUTE SMOOTHNESS INDICATOR and NORMALIZE ENTROPY RESIDUAL //
1899  // NOTE: see NCLS.h for a different but equivalent implementation of this.
1900  int ij = 0;
1901  for (int i=0; i<numDOFs; i++)
1902  {
1903  double etaMaxi, etaMini;
1904  if (STABILIZATION_TYPE==2) //EV
1905  {
1906  // For eta min and max
1907  etaMaxi = fabs(eta[i]);
1908  etaMini = fabs(eta[i]);
1909  }
1910  double porosity_times_solni = porosity_dof.data()[i]*u_dof_old.data()[i];
1911  // for smoothness indicator //
1912  double alpha_numerator = 0., alpha_denominator = 0.;
1913  for (int offset=csrRowIndeces_DofLoops.data()[i]; offset<csrRowIndeces_DofLoops.data()[i+1]; offset++)
1914  { // First loop in j (sparsity pattern)
1915  int j = csrColumnOffsets_DofLoops.data()[offset];
1916  if (STABILIZATION_TYPE==2) //EV Stabilization
1917  {
1918  // COMPUTE ETA MIN AND ETA MAX //
1919  etaMaxi = fmax(etaMaxi,fabs(eta[j]));
1920  etaMini = fmin(etaMini,fabs(eta[j]));
1921  }
1922  double porosity_times_solnj = porosity_dof.data()[j]*u_dof_old.data()[j];
1923  alpha_numerator += porosity_times_solni - porosity_times_solnj;
1924  alpha_denominator += fabs(porosity_times_solni - porosity_times_solnj);
1925  //update ij
1926  ij+=1;
1927  }
1928  if (STABILIZATION_TYPE==2) //EV Stab
1929  {
1930  // Normalize entropy residual
1931  global_entropy_residual[i] *= etaMini == etaMaxi ? 0. : 2*cE/(etaMaxi-etaMini);
1932  quantDOFs[i] = fabs(global_entropy_residual[i]);
1933  }
1934 
1935  double alphai = alpha_numerator/(alpha_denominator+1E-15);
1936  quantDOFs[i] = alphai;
1937 
1938  if (POWER_SMOOTHNESS_INDICATOR==0)
1939  psi[i] = 1.0;
1940  else
1941  psi[i] = std::pow(alphai,POWER_SMOOTHNESS_INDICATOR); //NOTE: they use alpha^2 in the paper
1942  }
1944  // ** LOOP IN DOFs FOR EDGE BASED TERMS ** //
1946  ij=0;
1947  for (int i=0; i<numDOFs; i++)
1948  {
1949  // NOTE: Transport matrices already have the porosity considered. ---> Dissipation matrices as well.
1950  double solni = u_dof_old.data()[i]; // solution at time tn for the ith DOF
1951  double porosityi = porosity_dof.data()[i];
1952  double ith_dissipative_term = 0;
1953  double ith_low_order_dissipative_term = 0;
1954  double ith_flux_term = 0;
1955  double dLii = 0.;
1956 
1957  // loop over the sparsity pattern of the i-th DOF
1958  for (int offset=csrRowIndeces_DofLoops.data()[i]; offset<csrRowIndeces_DofLoops.data()[i+1]; offset++)
1959  {
1960  int j = csrColumnOffsets_DofLoops.data()[offset];
1961  double solnj = u_dof_old.data()[j]; // solution at time tn for the jth DOF
1962  double porosityj = porosity_dof.data()[j];
1963  double dLowij, dLij, dEVij, dHij;
1964 
1965  ith_flux_term += TransportMatrix[ij]*solnj;
1966  if (i != j) //NOTE: there is really no need to check for i!=j (see formula for ith_dissipative_term)
1967  {
1968  // artificial compression
1969  double solij = 0.5*(porosityi*solni+porosityj*solnj);
1970  double Compij = cK*fmax(solij*(1.0-solij),0.0)/(fabs(porosityi*solni-porosityj*solnj)+1E-14);
1971  // first-order dissipative operator
1972  dLowij = fmax(fabs(TransportMatrix[ij]),fabs(TransposeTransportMatrix[ij]));
1973  //dLij = fmax(0.,fmax(psi[i]*TransportMatrix[ij], // Approach by S. Badia
1974  // psi[j]*TransposeTransportMatrix[ij]));
1975  dLij = dLowij*fmax(psi[i],psi[j]); // Approach by JLG & BP
1976  if (STABILIZATION_TYPE==2) //EV Stab
1977  {
1978  // high-order (entropy viscosity) dissipative operator
1979  dEVij = fmax(fabs(global_entropy_residual[i]),fabs(global_entropy_residual[j]));
1980  dHij = fmin(dLowij,dEVij) * fmax(1.0-Compij,0.0); // artificial compression
1981  }
1982  else // smoothness based indicator
1983  {
1984  dHij = dLij * fmax(1.0-Compij,0.0); // artificial compression
1985  }
1986  //dissipative terms
1987  ith_dissipative_term += dHij*(solnj-solni);
1988  ith_low_order_dissipative_term += dLowij*(solnj-solni);
1989  //dHij - dLij. This matrix is needed during FCT step
1990  dt_times_dH_minus_dL[ij] = dt*(dHij - dLowij);
1991  dLii -= dLij;
1992  dLow[ij] = dLowij;
1993  }
1994  else //i==j
1995  {
1996  // NOTE: this is incorrect. Indeed, dLii = -sum_{j!=i}(dLij) and similarly for dCii.
1997  // However, it is irrelevant since during the FCT step we do (dL-dC)*(solnj-solni)
1998  dt_times_dH_minus_dL[ij]=0;
1999  dLow[ij]=0;
2000  }
2001  //update ij
2002  ij+=1;
2003  }
2004  double mi = ML.data()[i];
2005  // compute edge_based_cfl
2006  edge_based_cfl.data()[i] = 2.*fabs(dLii)/mi;
2007  uLow[i] = u_dof_old.data()[i] - dt/mi*(ith_flux_term
2008  + boundary_integral[i]
2009  - ith_low_order_dissipative_term);
2010 
2011  // update residual
2012  if (LUMPED_MASS_MATRIX==1)
2013  globalResidual.data()[i] = u_dof_old.data()[i] - dt/mi*(ith_flux_term
2014  + boundary_integral[i]
2015  - ith_dissipative_term);
2016  else
2017  globalResidual.data()[i] += dt*(ith_flux_term - ith_dissipative_term);
2018  }
2019  }
2020  };//VOF
2021 
2022  inline VOF_base* newVOF(int nSpaceIn,
2023  int nQuadraturePoints_elementIn,
2024  int nDOF_mesh_trial_elementIn,
2025  int nDOF_trial_elementIn,
2026  int nDOF_test_elementIn,
2027  int nQuadraturePoints_elementBoundaryIn,
2028  int CompKernelFlag)
2029  {
2030  if (nSpaceIn == 2)
2031  return proteus::chooseAndAllocateDiscretization2D<VOF_base,VOF,CompKernel>(nSpaceIn,
2032  nQuadraturePoints_elementIn,
2033  nDOF_mesh_trial_elementIn,
2034  nDOF_trial_elementIn,
2035  nDOF_test_elementIn,
2036  nQuadraturePoints_elementBoundaryIn,
2037  CompKernelFlag);
2038  else
2039  return proteus::chooseAndAllocateDiscretization<VOF_base,VOF,CompKernel>(nSpaceIn,
2040  nQuadraturePoints_elementIn,
2041  nDOF_mesh_trial_elementIn,
2042  nDOF_trial_elementIn,
2043  nDOF_test_elementIn,
2044  nQuadraturePoints_elementBoundaryIn,
2045  CompKernelFlag);
2046  }
2047 }//proteus
2048 #endif
proteus::newVOF
VOF_base * newVOF(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition: VOF.h:2022
proteus::VOF_base::maxEntRes
std::valarray< double > maxEntRes
Definition: VOF.h:53
proteus::VOF::exteriorNumericalAdvectiveFluxDerivative
void exteriorNumericalAdvectiveFluxDerivative(const int &isDOFBoundary_u, const int &isFluxBoundary_u, const double n[nSpace], const double velocity[nSpace], double &dflux)
Definition: VOF.h:215
proteus::VOF_base::boundary_integral
std::valarray< double > boundary_integral
Definition: VOF.h:52
proteus::VOF::ck
CompKernelType ck
Definition: VOF.h:72
proteus::VOF_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)=0
proteus::VOF_base::maxVel
std::valarray< double > maxVel
Definition: VOF.h:53
proteus::VOF::nDOF_test_X_trial_element
const int nDOF_test_X_trial_element
Definition: VOF.h:71
proteus::VOF::calculateCFL
void calculateCFL(const double &elementDiameter, const double df[nSpace], double &cfl)
Definition: VOF.h:79
proteus::DENTROPY
double DENTROPY(const double &phi, const double &dummyL, const double &dummyR)
Definition: NCLS.h:25
proteus::VOF::calculateJacobian
void calculateJacobian(arguments_dict &args)
Definition: VOF.h:884
n
Int n
Definition: Headers.h:28
df
double df(double C, double b, double a, int q, int r)
Definition: analyticalSolutions.c:2209
proteus::VOF_base
Definition: VOF.h:46
ModelFactory.h
CompKernel.h
proteus::arguments_dict::scalar
T & scalar(const std::string &key)
proteus::VOF::calculateSubgridError_tau
void calculateSubgridError_tau(const double &elementDiameter, const double &dmt, const double dH[nSpace], double &cfl, double &tau)
Definition: VOF.h:111
proteus::VOF
Definition: VOF.h:69
proteus::arguments_dict::array
xt::pyarray< T > & array(const std::string &key)
proteus::VOF::calculateNumericalDiffusion
void calculateNumericalDiffusion(const double &shockCapturingDiffusion, const double &elementDiameter, const double &strong_residual, const double grad_u[nSpace], double &numDiff)
Definition: VOF.h:145
proteus::VOF::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: VOF.h:165
proteus::VOF::FCTStep
void FCTStep(arguments_dict &args)
Definition: VOF.h:1342
proteus::VOF::calculateResidualEdgeBased
void calculateResidualEdgeBased(arguments_dict &args)
Definition: VOF.h:1457
proteus::VOF::evaluateCoefficients
void evaluateCoefficients(const double v[nSpace], const double &u, const double &porosity, double &m, double &dm, double f[nSpace], double df[nSpace])
Definition: VOF.h:93
proteus::VOF_base::Rpos
std::valarray< double > Rpos
Definition: VOF.h:49
num
Int num
Definition: Headers.h:32
H
Double H
Definition: Headers.h:65
v
Double v
Definition: Headers.h:95
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
POWER_SMOOTHNESS_INDICATOR
#define POWER_SMOOTHNESS_INDICATOR
Definition: VOF.h:13
proteus::VOF_base::eta
std::valarray< double > eta
Definition: VOF.h:52
proteus::DENTROPY_LOG
double DENTROPY_LOG(const double &phi, const double &phiL, const double &phiR)
Definition: NCLS.h:32
proteus::VOF_base::Rneg
std::valarray< double > Rneg
Definition: VOF.h:49
proteus::VOF::VOF
VOF()
Definition: VOF.h:73
z
Double * z
Definition: Headers.h:49
u
Double u
Definition: Headers.h:89
proteus::phi
double phi(const double &g, const double &h, const double &hL, const double &hR, const double &uL, const double &uR)
Definition: SW2DCV.h:62
proteus::ENTROPY_LOG
double ENTROPY_LOG(const double &phi, const double &phiL, const double &phiR)
Definition: NCLS.h:29
proteus::VOF_base::psi
std::valarray< double > psi
Definition: VOF.h:52
xt
Definition: AddedMass.cpp:7
proteus::VOF_base::FluxCorrectionMatrix
std::valarray< double > FluxCorrectionMatrix
Definition: VOF.h:50
cMax
#define cMax
Definition: NCLS3P.h:11
proteus::VOF_base::TransposeTransportMatrix
std::valarray< double > TransposeTransportMatrix
Definition: VOF.h:51
proteus::VOF_base::calculateResidualElementBased
virtual void calculateResidualElementBased(arguments_dict &args)=0
proteus::VOF_base::~VOF_base
virtual ~VOF_base()
Definition: VOF.h:54
proteus::VOF_base::calculateResidualEdgeBased
virtual void calculateResidualEdgeBased(arguments_dict &args)=0
proteus::VOF_base::global_entropy_residual
std::valarray< double > global_entropy_residual
Definition: VOF.h:52
proteus::VOF_base::FCTStep
virtual void FCTStep(arguments_dict &args)=0
proteus
Definition: ADR.h:17
proteus::VOF_base::TransportMatrix
std::valarray< double > TransportMatrix
Definition: VOF.h:51
proteus::f
double f(const double &g, const double &h, const double &hZ)
Definition: SW2DCV.h:58
proteus::arguments_dict
Definition: ArgumentsDict.h:70
GLOBAL_FCT
#define GLOBAL_FCT
Definition: VOF.h:15
proteus::VOF::calculateResidualElementBased
void calculateResidualElementBased(arguments_dict &args)
Definition: VOF.h:250
proteus::VOF::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: VOF.h:129
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