proteus  1.8.1
C/C++/Fortran libraries
VOF3P.h
Go to the documentation of this file.
1 #ifndef VOF3P_H
2 #define VOF3P_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 #define POWER_SMOOTHNESS_INDICATOR 2
12 #define IS_BETAij_ONE 0
13 #define GLOBAL_FCT 0
14 
15 // Cell based methods:
16 // * SUPG with BDF1 or BDF2 time integration
17 // * Explicit Taylor Galerkin with EV stabilization
18 // Edge based methods.
19 // Low order via D. Kuzmin's
20 // High order methods: Smoothness indicator with MC, EV commutator with MC, D.K with ML
21 // Zalesak's FCT
22 
23 namespace proteus
24 {
25  // Power entropy //
26  inline double ENTROPY(const double& phi, const double& phiL, const double& phiR){
27  return 1./2.*std::pow(fabs(phi),2.);
28  }
29  inline double DENTROPY(const double& phi, const double& phiL, const double& phiR){
30  return fabs(phi)*(phi>=0 ? 1 : -1);
31  }
32  // Log entropy // for level set from 0 to 1
33  inline double ENTROPY_LOG(const double& phi, const double& phiL, const double& phiR){
34  return std::log(fabs((phi-phiL)*(phiR-phi))+1E-14);
35  }
36  inline double DENTROPY_LOG(const double& phi, const double& phiL, const double& phiR){
37  return (phiL+phiR-2*phi)*((phi-phiL)*(phiR-phi)>=0 ? 1 : -1)/(fabs((phi-phiL)*(phiR-phi))+1E-14);
38  }
39 }
40 
41 namespace proteus
42 {
43  class cppVOF3P_base
44  {
45  //The base class defining the interface
46  public:
47  std::valarray<double> Rpos, Rneg;
48  std::valarray<double> FluxCorrectionMatrix;
49  std::valarray<double> TransportMatrix, TransposeTransportMatrix;
50  std::valarray<double> psi, eta, global_entropy_residual, boundary_integral;
51  virtual ~cppVOF3P_base(){}
52  virtual void calculateResidualElementBased(arguments_dict& args)=0;
53  virtual void calculateJacobian(arguments_dict& args)=0;
54  virtual void FCTStep(arguments_dict& args)=0;
55  virtual void calculateResidualEdgeBased(arguments_dict& args)=0;
56  };
57 
58  template<class CompKernelType,
59  int nSpace,
60  int nQuadraturePoints_element,
61  int nDOF_mesh_trial_element,
62  int nDOF_trial_element,
63  int nDOF_test_element,
64  int nQuadraturePoints_elementBoundary>
65  class cppVOF : public cppVOF3P_base
66  {
67  public:
68  const int nDOF_test_X_trial_element;
69  CompKernelType ck;
70  cppVOF():
71  nDOF_test_X_trial_element(nDOF_test_element*nDOF_trial_element),
72  ck()
73  {}
74 
75  inline
76  void calculateCFL(const double& elementDiameter,
77  const double df[nSpace],
78  double& cfl)
79  {
80  double h,nrm_v;
81  h = elementDiameter;
82  nrm_v=0.0;
83  for(int I=0;I<nSpace;I++)
84  nrm_v+=df[I]*df[I];
85  nrm_v = sqrt(nrm_v);
86  cfl = nrm_v/h;
87  }
88 
89  inline
90  void evaluateCoefficients(const double v[nSpace],
91  const double& u,
92  const double& porosity, //VRANS specific
93  double& m,
94  double& dm,
95  double f[nSpace],
96  double df[nSpace])
97  {
98  m = porosity*u;
99  dm= porosity;
100  for (int I=0; I < nSpace; I++)
101  {
102  f[I] = v[I]*porosity*u;
103  df[I] = v[I]*porosity;
104  }
105  }
106 
107  inline
108  void calculateSubgridError_tau(const double& elementDiameter,
109  const double& dmt,
110  const double dH[nSpace],
111  double& cfl,
112  double& tau)
113  {
114  double h,nrm_v,oneByAbsdt;
115  h = elementDiameter;
116  nrm_v=0.0;
117  for(int I=0;I<nSpace;I++)
118  nrm_v+=dH[I]*dH[I];
119  nrm_v = sqrt(nrm_v);
120  cfl = nrm_v/h;
121  oneByAbsdt = fabs(dmt);
122  tau = 1.0/(2.0*nrm_v/h + oneByAbsdt + 1.0e-8);
123  }
124 
125  inline
126  void calculateSubgridError_tau( const double& Ct_sge,
127  const double G[nSpace*nSpace],
128  const double& A0,
129  const double Ai[nSpace],
130  double& tau_v,
131  double& cfl)
132  {
133  double v_d_Gv=0.0;
134  for(int I=0;I<nSpace;I++)
135  for (int J=0;J<nSpace;J++)
136  v_d_Gv += Ai[I]*G[I*nSpace+J]*Ai[J];
137 
138  tau_v = 1.0/sqrt(Ct_sge*A0*A0 + v_d_Gv + 1.0e-8);
139  }
140 
141  inline
142  void calculateNumericalDiffusion(const double& shockCapturingDiffusion,
143  const double& elementDiameter,
144  const double& strong_residual,
145  const double grad_u[nSpace],
146  double& numDiff)
147  {
148  double h,
149  num,
150  den,
151  n_grad_u;
152  h = elementDiameter;
153  n_grad_u = 0.0;
154  for (int I=0;I<nSpace;I++)
155  n_grad_u += grad_u[I]*grad_u[I];
156  num = shockCapturingDiffusion*0.5*h*fabs(strong_residual);
157  den = sqrt(n_grad_u) + 1.0e-8;
158  numDiff = num/den;
159  }
160 
161  inline
162  void exteriorNumericalAdvectiveFlux(const int& isDOFBoundary_u,
163  const int& isFluxBoundary_u,
164  const double n[nSpace],
165  const double& bc_u,
166  const double& bc_flux_u,
167  const double& u,
168  const double velocity[nSpace],
169  double& flux)
170  {
171 
172  double flow=0.0;
173  for (int I=0; I < nSpace; I++)
174  flow += n[I]*velocity[I];
175  //std::cout<<" isDOFBoundary_u= "<<isDOFBoundary_u<<" flow= "<<flow<<std::endl;
176  if (isDOFBoundary_u == 1)
177  {
178  //std::cout<<"Dirichlet boundary u and bc_u "<<u<<'\t'<<bc_u<<std::endl;
179  if (flow >= 0.0)
180  {
181  flux = u*flow;
182  //flux = flow;
183  }
184  else
185  {
186  flux = bc_u*flow;
187  //flux = flow;
188  }
189  }
190  else if (isFluxBoundary_u == 1)
191  {
192  flux = bc_flux_u;
193  //std::cout<<"Flux boundary flux and flow"<<flux<<'\t'<<flow<<std::endl;
194  }
195  else
196  {
197  //std::cout<<"No BC boundary flux and flow"<<flux<<'\t'<<flow<<std::endl;
198  if (flow >= 0.0)
199  {
200  flux = u*flow;
201  }
202  else
203  {
204  std::cout<<"warning: VOF open boundary with no external trace, setting to zero for inflow"<<std::endl;
205  flux = 0.0;
206  }
207 
208  }
209  }
210 
211  inline
212  void exteriorNumericalAdvectiveFluxDerivative(const int& isDOFBoundary_u,
213  const int& isFluxBoundary_u,
214  const double n[nSpace],
215  const double velocity[nSpace],
216  double& dflux)
217  {
218  double flow=0.0;
219  for (int I=0; I < nSpace; I++)
220  flow += n[I]*velocity[I];
221  //double flow=n[0]*velocity[0]+n[1]*velocity[1]+n[2]*velocity[2];
222  dflux=0.0;//default to no flux
223  if (isDOFBoundary_u == 1)
224  {
225  if (flow >= 0.0)
226  {
227  dflux = flow;
228  }
229  else
230  {
231  dflux = 0.0;
232  }
233  }
234  else if (isFluxBoundary_u == 1)
235  {
236  dflux = 0.0;
237  }
238  else
239  {
240  if (flow >= 0.0)
241  {
242  dflux = flow;
243  }
244  }
245  }
246 
247  void calculateResidualElementBased(arguments_dict& args)
248  {
249  double dt = args.scalar<double>("dt");
250  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
251  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
252  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
253  xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
254  double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
255  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
256  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
257  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
258  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
259  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
260  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
261  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
262  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
263  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
264  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
265  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
266  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
267  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
268  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
269  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
270  int nElements_global = args.scalar<int>("nElements_global");
271  double useMetrics = args.scalar<double>("useMetrics");
272  double alphaBDF = args.scalar<double>("alphaBDF");
273  int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
274  double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
275  double sc_uref = args.scalar<double>("sc_uref");
276  double sc_alpha = args.scalar<double>("sc_alpha");
277  const xt::pyarray<double>& q_vos = args.array<double>("q_vos");
278  const xt::pyarray<double>& vos_dof = args.array<double>("vos_dof");
279  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
280  xt::pyarray<int>& r_l2g = args.array<int>("r_l2g");
281  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
282  double degree_polynomial = args.scalar<double>("degree_polynomial");
283  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
284  xt::pyarray<double>& u_dof_old = args.array<double>("u_dof_old");
285  xt::pyarray<double>& velocity = args.array<double>("velocity");
286  xt::pyarray<double>& q_m = args.array<double>("q_m");
287  xt::pyarray<double>& q_u = args.array<double>("q_u");
288  xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
289  xt::pyarray<double>& q_dV = args.array<double>("q_dV");
290  xt::pyarray<double>& q_dV_last = args.array<double>("q_dV_last");
291  xt::pyarray<double>& cfl = args.array<double>("cfl");
292  xt::pyarray<double>& edge_based_cfl = args.array<double>("edge_based_cfl");
293  xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
294  xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
295  int offset_u = args.scalar<int>("offset_u");
296  int stride_u = args.scalar<int>("stride_u");
297  xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
298  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
299  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
300  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
301  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
302  xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
303  const xt::pyarray<double>& ebqe_vos_ext = args.array<double>("ebqe_vos_ext");
304  xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
305  xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
306  xt::pyarray<int>& isFluxBoundary_u = args.array<int>("isFluxBoundary_u");
307  xt::pyarray<double>& ebqe_bc_flux_u_ext = args.array<double>("ebqe_bc_flux_u_ext");
308  xt::pyarray<double>& ebqe_phi = args.array<double>("ebqe_phi");
309  double epsFact = args.scalar<double>("epsFact");
310  xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
311  xt::pyarray<double>& ebqe_flux = args.array<double>("ebqe_flux");
312  int stage = args.scalar<int>("stage");
313  xt::pyarray<double>& uTilde_dof = args.array<double>("uTilde_dof");
314  double cE = args.scalar<double>("cE");
315  double cMax = args.scalar<double>("cMax");
316  double cK = args.scalar<double>("cK");
317  double uL = args.scalar<double>("uL");
318  double uR = args.scalar<double>("uR");
319  int numDOFs = args.scalar<int>("numDOFs");
320  int NNZ = args.scalar<int>("NNZ");
321  xt::pyarray<int>& csrRowIndeces_DofLoops = args.array<int>("csrRowIndeces_DofLoops");
322  xt::pyarray<int>& csrColumnOffsets_DofLoops = args.array<int>("csrColumnOffsets_DofLoops");
323  xt::pyarray<int>& csrRowIndeces_CellLoops = args.array<int>("csrRowIndeces_CellLoops");
324  xt::pyarray<int>& csrColumnOffsets_CellLoops = args.array<int>("csrColumnOffsets_CellLoops");
325  xt::pyarray<int>& csrColumnOffsets_eb_CellLoops = args.array<int>("csrColumnOffsets_eb_CellLoops");
326  xt::pyarray<double>& ML = args.array<double>("ML");
327  int LUMPED_MASS_MATRIX = args.scalar<int>("LUMPED_MASS_MATRIX");
328  int STABILIZATION_TYPE = args.scalar<int>("STABILIZATION_TYPE");
329  int ENTROPY_TYPE = args.scalar<int>("ENTROPY_TYPE");
330  xt::pyarray<double>& uLow = args.array<double>("uLow");
331  xt::pyarray<double>& dLow = args.array<double>("dLow");
332  xt::pyarray<double>& dt_times_dH_minus_dL = args.array<double>("dt_times_dH_minus_dL");
333  xt::pyarray<double>& min_u_bc = args.array<double>("min_u_bc");
334  xt::pyarray<double>& max_u_bc = args.array<double>("max_u_bc");
335  xt::pyarray<double>& quantDOFs = args.array<double>("quantDOFs");
336  double meanEntropy = 0., meanOmega = 0., maxEntropy = -1E10, minEntropy = 1E10;
337  register double maxVel[nElements_global], maxEntRes[nElements_global];
338  double Ct_sge = 4.0;
339  //
340  //loop over elements to compute volume integrals and load them into element and global residual
341  //
342  //eN is the element index
343  //eN_k is the quadrature point index for a scalar
344  //eN_k_nSpace is the quadrature point index for a vector
345  //eN_i is the element test function index
346  //eN_j is the element trial function index
347  //eN_k_j is the quadrature point index for a trial function
348  //eN_k_i is the quadrature point index for a trial function
349  for(int eN=0;eN<nElements_global;eN++)
350  {
351  // init maxVel and maxEntRes
352  maxVel[eN] = 0.;
353  maxEntRes[eN] = 0.;
354  //declare local storage for element residual and initialize
355  register double elementResidual_u[nDOF_test_element];
356  for (int i=0;i<nDOF_test_element;i++)
357  {
358  elementResidual_u[i]=0.0;
359  }//i
360  //loop over quadrature points and compute integrands
361  for (int k=0;k<nQuadraturePoints_element;k++)
362  {
363  //compute indeces and declare local storage
364  register int eN_k = eN*nQuadraturePoints_element+k,
365  eN_k_nSpace = eN_k*nSpace,
366  eN_nDOF_trial_element = eN*nDOF_trial_element;
367  register double
368  entVisc_minus_artComp,
369  u=0.0,un=0.0,
370  grad_u[nSpace],grad_u_old[nSpace],grad_uTilde[nSpace],
371  m=0.0,dm=0.0,
372  H=0.0,Hn=0.0,HTilde=0.0,
373  f[nSpace],fn[nSpace],df[nSpace],
374  m_t=0.0,dm_t=0.0,
375  pdeResidual_u=0.0,
376  Lstar_u[nDOF_test_element],
377  subgridError_u=0.0,
378  tau=0.0,tau0=0.0,tau1=0.0,
379  numDiff0=0.0,numDiff1=0.0,
380  jac[nSpace*nSpace],
381  jacDet,
382  jacInv[nSpace*nSpace],
383  u_grad_trial[nDOF_trial_element*nSpace],
384  u_test_dV[nDOF_trial_element],
385  u_grad_test_dV[nDOF_test_element*nSpace],
386  dV,x,y,z,xt,yt,zt,
387  //VRANS
388  porosity,
389  //
390  G[nSpace*nSpace],G_dd_G,tr_G;//norm_Rv;
391  // //
392  // //compute solution and gradients at quadrature points
393  // //
394  // u=0.0;
395  // for (int I=0;I<nSpace;I++)
396  // {
397  // grad_u[I]=0.0;
398  // }
399  // for (int j=0;j<nDOF_trial_element;j++)
400  // {
401  // int eN_j=eN*nDOF_trial_element+j;
402  // int eN_k_j=eN_k*nDOF_trial_element+j;
403  // int eN_k_j_nSpace = eN_k_j*nSpace;
404  // u += valFromDOF_c(u_dof[u_l2g[eN_j]],u_trial[eN_k_j]);
405  // for (int I=0;I<nSpace;I++)
406  // {
407  // grad_u[I] += gradFromDOF_c(u_dof[u_l2g[eN_j]],u_grad_trial[eN_k_j_nSpace+I]);
408  // }
409  // }
410  ck.calculateMapping_element(eN,
411  k,
412  mesh_dof.data(),
413  mesh_l2g.data(),
414  mesh_trial_ref.data(),
415  mesh_grad_trial_ref.data(),
416  jac,
417  jacDet,
418  jacInv,
419  x,y,z);
420  ck.calculateMappingVelocity_element(eN,
421  k,
422  mesh_velocity_dof.data(),
423  mesh_l2g.data(),
424  mesh_trial_ref.data(),
425  xt,yt,zt);
426  //get the physical integration weight
427  dV = fabs(jacDet)*dV_ref[k];
428  ck.calculateG(jacInv,G,G_dd_G,tr_G);
429  //get the trial function gradients
430  ck.gradTrialFromRef(&u_grad_trial_ref[k*nDOF_trial_element*nSpace],
431  jacInv,
432  u_grad_trial);
433  //get the solution
434  ck.valFromDOF(u_dof.data(),
435  &u_l2g[eN_nDOF_trial_element],
436  &u_trial_ref[k*nDOF_trial_element],
437  u);
438  ck.valFromDOF(u_dof_old.data(),
439  &u_l2g[eN_nDOF_trial_element],
440  &u_trial_ref[k*nDOF_trial_element],
441  un);
442  //get the solution gradients
443  ck.gradFromDOF(u_dof.data(),
444  &u_l2g[eN_nDOF_trial_element],
445  u_grad_trial,
446  grad_u);
447  ck.gradFromDOF(u_dof_old.data(),
448  &u_l2g[eN_nDOF_trial_element],
449  u_grad_trial,
450  grad_u_old);
451  ck.gradFromDOF(uTilde_dof.data(),
452  &u_l2g[eN_nDOF_trial_element],
453  u_grad_trial,
454  grad_uTilde);
455  //precalculate test function products with integration weights
456  for (int j=0;j<nDOF_trial_element;j++)
457  {
458  u_test_dV[j] = u_test_ref[k*nDOF_trial_element+j]*dV;
459  for (int I=0;I<nSpace;I++)
460  {
461  u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
462  }
463  }
464  //VRANS
465  porosity = 1.0 - q_vos[eN_k];
466  //
467  //
468  //calculate pde coefficients at quadrature points
469  //
470  evaluateCoefficients(&velocity[eN_k_nSpace],
471  u,
472  //VRANS
473  porosity,
474  //
475  m,
476  dm,
477  f,
478  df);
479  //
480  //moving mesh
481  //
482  double mesh_velocity[3];
483  mesh_velocity[0] = xt;
484  mesh_velocity[1] = yt;
485  mesh_velocity[2] = zt;
486 
487  for (int I=0;I<nSpace;I++)
488  {
489  f[I] -= MOVING_DOMAIN*m*mesh_velocity[I];
490  df[I] -= MOVING_DOMAIN*dm*mesh_velocity[I];
491  }
492  //
493  //calculate time derivative at quadrature points
494  //
495  if (q_dV_last[eN_k] <= -100)
496  q_dV_last[eN_k] = dV;
497  q_dV[eN_k] = dV;
498  ck.bdf(alphaBDF,
499  q_m_betaBDF[eN_k]*q_dV_last[eN_k]/dV,//ensure prior mass integral is correct for m_t with BDF1
500  m,
501  dm,
502  m_t,
503  dm_t);
504 
505  if (STABILIZATION_TYPE==1)
506  {
507  double normVel=0., norm_grad_un=0.;
508  for (int I=0;I<nSpace;I++)
509  {
510  Hn += df[I]*grad_u_old[I];
511  HTilde += df[I]*grad_uTilde[I];
512  fn[I] = porosity*df[I]*un-MOVING_DOMAIN*m*mesh_velocity[I];
513  H += df[I]*grad_u[I];
514  normVel += df[I]*df[I];
515  norm_grad_un += grad_u_old[I]*grad_u_old[I];
516  }
517  normVel = std::sqrt(normVel);
518  norm_grad_un = std::sqrt(norm_grad_un)+1E-10;
519 
520  // calculate CFL
521  calculateCFL(elementDiameter[eN]/degree_polynomial,df,cfl[eN_k]);
522 
523 
524  // compute max velocity at cell
525  maxVel[eN] = fmax(normVel,maxVel[eN]);
526 
527  // Strong entropy residual
528  double entRes = (ENTROPY(u,0,1)-ENTROPY(un,0,1))/dt + 0.5*(DENTROPY(u,0,1)*H +
529  DENTROPY(un,0,1)*Hn);
530  maxEntRes[eN] = fmax(maxEntRes[eN],fabs(entRes));
531 
532  // Quantities for normalization factor //
533  meanEntropy += ENTROPY(u,0,1)*dV;
534  meanOmega += dV;
535  maxEntropy = fmax(maxEntropy,ENTROPY(u,0,1));
536  minEntropy = fmin(minEntropy,ENTROPY(u,0,1));
537 
538  // artificial compression
539  double hK=elementDiameter[eN]/degree_polynomial;
540  entVisc_minus_artComp = fmax(1-cK*fmax(un*(1-un),0)/hK/norm_grad_un,0);
541  }
542  else
543  {
544  //
545  //calculate subgrid error (strong residual and adjoint)
546  //
547  //calculate strong residual
548  pdeResidual_u = ck.Mass_strong(m_t) + ck.Advection_strong(df,grad_u);
549  //calculate adjoint
550  for (int i=0;i<nDOF_test_element;i++)
551  {
552  // register int eN_k_i_nSpace = (eN_k*nDOF_trial_element+i)*nSpace;
553  // Lstar_u[i] = ck.Advection_adjoint(df,&u_grad_test_dV[eN_k_i_nSpace]);
554  register int i_nSpace = i*nSpace;
555  Lstar_u[i] = ck.Advection_adjoint(df,&u_grad_test_dV[i_nSpace]);
556  }
557  //calculate tau and tau*Res
558  calculateSubgridError_tau(elementDiameter[eN],dm_t,df,cfl[eN_k],tau0);
559  calculateSubgridError_tau(Ct_sge,
560  G,
561  dm_t,
562  df,
563  tau1,
564  cfl[eN_k]);
565  tau = useMetrics*tau1+(1.0-useMetrics)*tau0;
566 
567  subgridError_u = -tau*pdeResidual_u;
568  //
569  //calculate shock capturing diffusion
570  //
571 
572  ck.calculateNumericalDiffusion(shockCapturingDiffusion,
573  elementDiameter[eN],
574  pdeResidual_u,
575  grad_u,
576  numDiff0);
577  //ck.calculateNumericalDiffusion(shockCapturingDiffusion,G,pdeResidual_u,grad_u_old,numDiff1);
578  ck.calculateNumericalDiffusion(shockCapturingDiffusion,
579  sc_uref,
580  sc_alpha,
581  G,
582  G_dd_G,
583  pdeResidual_u,
584  grad_u,
585  numDiff1);
586  q_numDiff_u[eN_k] = useMetrics*numDiff1+(1.0-useMetrics)*numDiff0;
587  //std::cout<<tau<<" "<<q_numDiff_u[eN_k]<<'\t'<<numDiff0<<'\t'<<numDiff1<<'\t'<<pdeResidual_u<<std::endl;
588 
589  //
590  //update element residual
591  //
592 
593 
594  /* std::cout<<m_t<<'\t'
595  <<f[0]<<'\t'
596  <<f[1]<<'\t'
597  <<df[0]<<'\t'
598  <<df[1]<<'\t'
599  <<subgridError_u<<'\t'
600  <<q_numDiff_u_last[eN_k]<<std::endl;*/
601  }
602 
603  for(int i=0;i<nDOF_test_element;i++)
604  {
605  //register int eN_k_i=eN_k*nDOF_test_element+i,
606  //eN_k_i_nSpace = eN_k_i*nSpace,
607  register int i_nSpace=i*nSpace;
608  if (STABILIZATION_TYPE==1)
609  {
610  if (stage == 1)
611  elementResidual_u[i] +=
612  ck.Mass_weak(dt*m_t,u_test_dV[i]) + // time derivative
613  1./3*dt*ck.Advection_weak(fn,&u_grad_test_dV[i_nSpace]) +
614  1./9*dt*dt*ck.NumericalDiffusion(Hn,df,&u_grad_test_dV[i_nSpace]) +
615  1./3*dt*entVisc_minus_artComp*ck.NumericalDiffusion(q_numDiff_u_last[eN_k],
616  grad_u_old,
617  &u_grad_test_dV[i_nSpace]);
618  // TODO: Add part about moving mesh
619  else //stage 2
620  elementResidual_u[i] +=
621  ck.Mass_weak(dt*m_t,u_test_dV[i]) + // time derivative
622  dt*ck.Advection_weak(fn,&u_grad_test_dV[i_nSpace]) +
623  0.5*dt*dt*ck.NumericalDiffusion(HTilde,df,&u_grad_test_dV[i_nSpace]) +
624  dt*entVisc_minus_artComp*ck.NumericalDiffusion(q_numDiff_u_last[eN_k],
625  grad_u_old,
626  &u_grad_test_dV[i_nSpace]);
627  }
628  else //supg
629  {
630  elementResidual_u[i] +=
631  ck.Mass_weak(m_t,u_test_dV[i]) +
632  ck.Advection_weak(f,&u_grad_test_dV[i_nSpace]) +
633  ck.SubgridError(subgridError_u,Lstar_u[i]) +
634  ck.NumericalDiffusion(q_numDiff_u_last[eN_k],
635  grad_u,
636  &u_grad_test_dV[i_nSpace]);
637  }
638  }//i
639  //
640  //cek/ido todo, get rid of m, since u=m
641  //save momentum for time history and velocity for subgrid error
642  //save solution for other models
643  //
644  q_u[eN_k] = u;
645  q_m[eN_k] = m;
646  }
647  //
648  //load element into global residual and save element residual
649  //
650  for(int i=0;i<nDOF_test_element;i++)
651  {
652  register int eN_i=eN*nDOF_test_element+i;
653  globalResidual[offset_u+stride_u*r_l2g[eN_i]] += elementResidual_u[i];
654  }//i
655  }//elements
656  //
657  //loop over exterior element boundaries to calculate surface integrals and load into element and global residuals
658  //
659  //ebNE is the Exterior element boundary INdex
660  //ebN is the element boundary INdex
661  //eN is the element index
662  for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
663  {
664  register int ebN = exteriorElementBoundariesArray[ebNE],
665  eN = elementBoundaryElementsArray[ebN*2+0],
666  ebN_local = elementBoundaryLocalElementBoundariesArray[ebN*2+0],
667  eN_nDOF_trial_element = eN*nDOF_trial_element;
668  register double elementResidual_u[nDOF_test_element];
669  for (int i=0;i<nDOF_test_element;i++)
670  {
671  elementResidual_u[i]=0.0;
672  }
673  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
674  {
675  register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
676  ebNE_kb_nSpace = ebNE_kb*nSpace,
677  ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
678  ebN_local_kb_nSpace = ebN_local_kb*nSpace;
679  register double u_ext=0.0,
680  grad_u_ext[nSpace],
681  m_ext=0.0,
682  dm_ext=0.0,
683  f_ext[nSpace],
684  df_ext[nSpace],
685  flux_ext=0.0,
686  bc_u_ext=0.0,
687  //bc_grad_u_ext[nSpace],
688  bc_m_ext=0.0,
689  bc_dm_ext=0.0,
690  bc_f_ext[nSpace],
691  bc_df_ext[nSpace],
692  jac_ext[nSpace*nSpace],
693  jacDet_ext,
694  jacInv_ext[nSpace*nSpace],
695  boundaryJac[nSpace*(nSpace-1)],
696  metricTensor[(nSpace-1)*(nSpace-1)],
697  metricTensorDetSqrt,
698  dS,
699  u_test_dS[nDOF_test_element],
700  u_grad_trial_trace[nDOF_trial_element*nSpace],
701  normal[nSpace],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
702  //VRANS
703  porosity_ext,
704  //
705  G[nSpace*nSpace],G_dd_G,tr_G;
706  //
707  //calculate the solution and gradients at quadrature points
708  //
709  //compute information about mapping from reference element to physical element
710  ck.calculateMapping_elementBoundary(eN,
711  ebN_local,
712  kb,
713  ebN_local_kb,
714  mesh_dof.data(),
715  mesh_l2g.data(),
716  mesh_trial_trace_ref.data(),
717  mesh_grad_trial_trace_ref.data(),
718  boundaryJac_ref.data(),
719  jac_ext,
720  jacDet_ext,
721  jacInv_ext,
722  boundaryJac,
723  metricTensor,
724  metricTensorDetSqrt,
725  normal_ref.data(),
726  normal,
727  x_ext,y_ext,z_ext);
728  ck.calculateMappingVelocity_elementBoundary(eN,
729  ebN_local,
730  kb,
731  ebN_local_kb,
732  mesh_velocity_dof.data(),
733  mesh_l2g.data(),
734  mesh_trial_trace_ref.data(),
735  xt_ext,yt_ext,zt_ext,
736  normal,
737  boundaryJac,
738  metricTensor,
739  integralScaling);
740 
741  //std::cout<<"metricTensorDetSqrt "<<metricTensorDetSqrt<<" integralScaling "<<integralScaling<<std::endl;
742  dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref[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[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[eN_nDOF_trial_element],
756  &u_trial_trace_ref[ebN_local_kb*nDOF_test_element],
757  u_ext);
758  ck.gradFromDOF(u_dof_old.data(),
759  &u_l2g[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[eN_nDOF_trial_element],
767  &u_trial_trace_ref[ebN_local_kb*nDOF_test_element],
768  u_ext);
769  ck.gradFromDOF(u_dof.data(),
770  &u_l2g[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[ebN_local_kb*nDOF_test_element+j]*dS;
778  }
779  //
780  //load the boundary values
781  //
782  bc_u_ext = isDOFBoundary_u[ebNE_kb]*ebqe_bc_u_ext[ebNE_kb]+(1-isDOFBoundary_u[ebNE_kb])*u_ext;
783  //VRANS
784  porosity_ext = 1.0 - ebqe_vos_ext[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[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[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[ebNE_kb],
827  isFluxBoundary_u[ebNE_kb],
828  normal,
829  bc_u_ext,
830  ebqe_bc_flux_u_ext[ebNE_kb],
831  u_ext,
832  df_ext,//VRANS includes porosity
833  flux_ext);
834  ebqe_flux[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[ebNE_kb] = u_ext;
838  else
839  ebqe_u[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[offset_u+stride_u*r_l2g[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[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[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_vos = args.array<double>("q_vos");
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_vos_ext = args.array<double>("ebqe_vos_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<<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[u_l2g[eN_j]],u_trial[eN_k_j]);
991  // for (int I=0;I<nSpace;I++)
992  // {
993  // grad_u[I] += gradFromDOF_c(u_dof[u_l2g[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[k];
1015  ck.calculateG(jacInv,G,G_dd_G,tr_G);
1016  //get the trial function gradients
1017  ck.gradTrialFromRef(&u_grad_trial_ref[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
1018  //get the solution
1019  ck.valFromDOF(u_dof.data(),&u_l2g[eN_nDOF_trial_element],&u_trial_ref[k*nDOF_trial_element],u);
1020  //get the solution gradients
1021  ck.gradFromDOF(u_dof.data(),&u_l2g[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[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 = 1.0 - q_vos[eN_k];
1033  //
1034  //
1035  //calculate pde coefficients and derivatives at quadrature points
1036  //
1037  evaluateCoefficients(&velocity[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[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[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[eN],
1091  dm_t,
1092  df,
1093  cfl[eN_k],
1094  tau0);
1095 
1096  calculateSubgridError_tau(Ct_sge,
1097  G,
1098  dm_t,
1099  df,
1100  tau1,
1101  cfl[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[k*nDOF_trial_element+j],
1122  u_test_dV[i]) +
1123  ck.AdvectionJacobian_weak(df,
1124  u_trial_ref[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[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[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[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_u_u[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[ebNE];
1161  register int eN = elementBoundaryElementsArray[ebN*2+0],
1162  ebN_local = elementBoundaryLocalElementBoundariesArray[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[u_l2g[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[u_l2g[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[kb];
1251  //dS = metricTensorDetSqrt*dS_ref[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[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[eN_nDOF_trial_element],&u_trial_trace_ref[ebN_local_kb*nDOF_test_element],u_ext);
1258  ck.gradFromDOF(u_dof.data(),&u_l2g[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[ebN_local_kb*nDOF_test_element+j]*dS;
1263  }
1264  //
1265  //load the boundary values
1266  //
1267  bc_u_ext = isDOFBoundary_u[ebNE_kb]*ebqe_bc_u_ext[ebNE_kb]+(1-isDOFBoundary_u[ebNE_kb])*u_ext;
1268  //VRANS
1269  porosity_ext = 1.0 - ebqe_vos_ext[ebNE_kb];
1270  //
1271  //
1272  //calculate the internal and external trace of the pde coefficients
1273  //
1274  evaluateCoefficients(&ebqe_velocity_ext[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[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[ebNE_kb],
1312  isFluxBoundary_u[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[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[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[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), Rneg.resize(numDOFs,0.0);
1362  FluxCorrectionMatrix.resize(NNZ,0.0);
1364  // LOOP in DOFs //
1366  int ij=0;
1367  for (int i=0; i<numDOFs; i++)
1368  {
1369  //read some vectors
1370  double solHi = solH[i];
1371  double solni = soln[i];
1372  double mi = lumped_mass_matrix[i];
1373  double uLowi = uLow[i];
1374  double uDotLowi = (uLowi - solni)/dt;
1375  double mini=min_u_bc[i], maxi=max_u_bc[i]; // init min/max with value at BCs (NOTE: if no boundary then min=1E10, max=-1E10)
1376  if (GLOBAL_FCT==1)
1377  {
1378  mini = 0.;
1379  maxi = 1.;
1380  }
1381 
1382  double Pposi=0, Pnegi=0;
1383  // LOOP OVER THE SPARSITY PATTERN (j-LOOP)//
1384  for (int offset=csrRowIndeces_DofLoops[i]; offset<csrRowIndeces_DofLoops[i+1]; offset++)
1385  {
1386  int j = csrColumnOffsets_DofLoops[offset];
1387  double solnj = soln[j];
1389  // COMPUTE THE BOUNDS //
1391  if (GLOBAL_FCT == 0)
1392  {
1393  mini = fmin(mini,solnj);
1394  maxi = fmax(maxi,solnj);
1395  }
1396  double uLowj = uLow[j];
1397  double uDotLowj = (uLowj - solnj)/dt;
1398  // i-th row of flux correction matrix
1399  if (STABILIZATION_TYPE==4) // DK high-order, linearly stable anti-dif. flux
1400  {
1401  FluxCorrectionMatrix[ij] = dt*(MassMatrix[ij]*(uDotLowi-uDotLowj)
1402  + dLow[ij]*(uLowi-uLowj));
1403  }
1404  else
1405  {
1406  double ML_minus_MC =
1407  (LUMPED_MASS_MATRIX == 1 ? 0. : (i==j ? 1. : 0.)*mi - MassMatrix[ij]);
1408  FluxCorrectionMatrix[ij] = ML_minus_MC * (solH[j]-solnj - (solHi-solni))
1409  + dt_times_dH_minus_dL[ij]*(solnj-solni);
1410  }
1411 
1413  // COMPUTE P VECTORS //
1415  Pposi += FluxCorrectionMatrix[ij]*((FluxCorrectionMatrix[ij] > 0) ? 1. : 0.);
1416  Pnegi += FluxCorrectionMatrix[ij]*((FluxCorrectionMatrix[ij] < 0) ? 1. : 0.);
1417 
1418  //update ij
1419  ij+=1;
1420  }
1422  // COMPUTE Q VECTORS //
1424  double Qposi = mi*(maxi-uLow[i]);
1425  double Qnegi = mi*(mini-uLow[i]);
1426 
1428  // COMPUTE R VECTORS //
1430  Rpos[i] = ((Pposi==0) ? 1. : fmin(1.0,Qposi/Pposi));
1431  Rneg[i] = ((Pnegi==0) ? 1. : fmin(1.0,Qnegi/Pnegi));
1432  } // i DOFs
1433 
1435  // COMPUTE LIMITERS //
1437  ij=0;
1438  for (int i=0; i<numDOFs; i++)
1439  {
1440  double ith_Limiter_times_FluxCorrectionMatrix = 0.;
1441  double Rposi = Rpos[i], Rnegi = Rneg[i];
1442  // LOOP OVER THE SPARSITY PATTERN (j-LOOP)//
1443  for (int offset=csrRowIndeces_DofLoops[i]; offset<csrRowIndeces_DofLoops[i+1]; offset++)
1444  {
1445  int j = csrColumnOffsets_DofLoops[offset];
1446  double Lij = 1;
1447  Lij = ((FluxCorrectionMatrix[ij]>0) ? fmin(Rposi,Rneg[j]) : fmin(Rnegi,Rpos[j]));
1448  ith_Limiter_times_FluxCorrectionMatrix += Lij * FluxCorrectionMatrix[ij];
1449  //update ij
1450  ij+=1;
1451  }
1452  limited_solution[i] = uLow[i] + 1./lumped_mass_matrix[i]*ith_Limiter_times_FluxCorrectionMatrix;
1453  }
1454  }
1455 
1456  void calculateResidualEdgeBased(arguments_dict& args)
1457  {
1458  double dt = args.scalar<double>("dt");
1459  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1460  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1461  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1462  xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
1463  double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
1464  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1465  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1466  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1467  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1468  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1469  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
1470  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
1471  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
1472  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
1473  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
1474  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
1475  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
1476  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
1477  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
1478  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
1479  int nElements_global = args.scalar<int>("nElements_global");
1480  double useMetrics = args.scalar<double>("useMetrics");
1481  double alphaBDF = args.scalar<double>("alphaBDF");
1482  int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
1483  double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
1484  double sc_uref = args.scalar<double>("sc_uref");
1485  double sc_alpha = args.scalar<double>("sc_alpha");
1486  const xt::pyarray<double>& q_vos = args.array<double>("q_vos");
1487  const xt::pyarray<double>& vos_dof = args.array<double>("vos_dof");
1488  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1489  xt::pyarray<int>& r_l2g = args.array<int>("r_l2g");
1490  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1491  double degree_polynomial = args.scalar<double>("degree_polynomial");
1492  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1493  xt::pyarray<double>& u_dof_old = args.array<double>("u_dof_old");
1494  xt::pyarray<double>& velocity = args.array<double>("velocity");
1495  xt::pyarray<double>& q_m = args.array<double>("q_m");
1496  xt::pyarray<double>& q_u = args.array<double>("q_u");
1497  xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
1498  xt::pyarray<double>& q_dV = args.array<double>("q_dV");
1499  xt::pyarray<double>& q_dV_last = args.array<double>("q_dV_last");
1500  xt::pyarray<double>& cfl = args.array<double>("cfl");
1501  xt::pyarray<double>& edge_based_cfl = args.array<double>("edge_based_cfl");
1502  xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
1503  xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
1504  int offset_u = args.scalar<int>("offset_u");
1505  int stride_u = args.scalar<int>("stride_u");
1506  xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
1507  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
1508  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
1509  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
1510  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
1511  xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
1512  const xt::pyarray<double>& ebqe_vos_ext = args.array<double>("ebqe_vos_ext");
1513  xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
1514  xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
1515  xt::pyarray<int>& isFluxBoundary_u = args.array<int>("isFluxBoundary_u");
1516  xt::pyarray<double>& ebqe_bc_flux_u_ext = args.array<double>("ebqe_bc_flux_u_ext");
1517  xt::pyarray<double>& ebqe_phi = args.array<double>("ebqe_phi");
1518  double epsFact = args.scalar<double>("epsFact");
1519  xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
1520  xt::pyarray<double>& ebqe_flux = args.array<double>("ebqe_flux");
1521  int stage = args.scalar<int>("stage");
1522  xt::pyarray<double>& uTilde_dof = args.array<double>("uTilde_dof");
1523  double cE = args.scalar<double>("cE");
1524  double cMax = args.scalar<double>("cMax");
1525  double cK = args.scalar<double>("cK");
1526  double uL = args.scalar<double>("uL");
1527  double uR = args.scalar<double>("uR");
1528  int numDOFs = args.scalar<int>("numDOFs");
1529  int NNZ = args.scalar<int>("NNZ");
1530  xt::pyarray<int>& csrRowIndeces_DofLoops = args.array<int>("csrRowIndeces_DofLoops");
1531  xt::pyarray<int>& csrColumnOffsets_DofLoops = args.array<int>("csrColumnOffsets_DofLoops");
1532  xt::pyarray<int>& csrRowIndeces_CellLoops = args.array<int>("csrRowIndeces_CellLoops");
1533  xt::pyarray<int>& csrColumnOffsets_CellLoops = args.array<int>("csrColumnOffsets_CellLoops");
1534  xt::pyarray<int>& csrColumnOffsets_eb_CellLoops = args.array<int>("csrColumnOffsets_eb_CellLoops");
1535  xt::pyarray<double>& ML = args.array<double>("ML");
1536  int LUMPED_MASS_MATRIX = args.scalar<int>("LUMPED_MASS_MATRIX");
1537  int STABILIZATION_TYPE = args.scalar<int>("STABILIZATION_TYPE");
1538  int ENTROPY_TYPE = args.scalar<int>("ENTROPY_TYPE");
1539  xt::pyarray<double>& uLow = args.array<double>("uLow");
1540  xt::pyarray<double>& dLow = args.array<double>("dLow");
1541  xt::pyarray<double>& dt_times_dH_minus_dL = args.array<double>("dt_times_dH_minus_dL");
1542  xt::pyarray<double>& min_u_bc = args.array<double>("min_u_bc");
1543  xt::pyarray<double>& max_u_bc = args.array<double>("max_u_bc");
1544  xt::pyarray<double>& quantDOFs = args.array<double>("quantDOFs");
1545  // NOTE: This function follows a different (but equivalent) implementation of the smoothness based indicator than NCLS.h
1546  // Allocate space for the transport matrices
1547  // This is used for first order KUZMIN'S METHOD
1548  TransportMatrix.resize(NNZ,0.0), TransposeTransportMatrix.resize(NNZ,0.0);
1549  for (int i=0; i<NNZ; i++)
1550  {
1551  TransportMatrix[i] = 0.;
1552  TransposeTransportMatrix[i] = 0.;
1553  }
1554 
1555  // compute entropy and init global_entropy_residual and boundary_integral
1556  psi.resize(numDOFs,0.0);
1557  eta.resize(numDOFs,0.0);
1558  global_entropy_residual.resize(numDOFs,0.0);
1559  boundary_integral.resize(numDOFs,0.0);
1560  for (int i=0; i<numDOFs; i++)
1561  {
1562  // NODAL ENTROPY //
1563  if (STABILIZATION_TYPE==2) //EV stab
1564  {
1565  double porosity_times_solni = (1.0-vos_dof[i])*u_dof_old[i];
1566  eta[i] = ENTROPY_TYPE == 0 ? ENTROPY(porosity_times_solni,uL,uR) : ENTROPY_LOG(porosity_times_solni,uL,uR);
1567  global_entropy_residual[i]=0.;
1568  }
1569  boundary_integral[i]=0.;
1570  }
1571 
1573  // ** LOOP IN CELLS FOR CELL BASED TERMS ** //
1575  // HERE WE COMPUTE:
1576  // * Time derivative term. porosity*u_t
1577  // * cell based CFL (for reference)
1578  // * Entropy residual
1579  // * Transport matrices
1580  for(int eN=0;eN<nElements_global;eN++)
1581  {
1582  //declare local storage for local contributions and initialize
1583  register double
1584  elementResidual_u[nDOF_test_element],
1585  element_entropy_residual[nDOF_test_element];
1586  register double elementTransport[nDOF_test_element][nDOF_trial_element];
1587  register double elementTransposeTransport[nDOF_test_element][nDOF_trial_element];
1588  for (int i=0;i<nDOF_test_element;i++)
1589  {
1590  elementResidual_u[i]=0.0;
1591  element_entropy_residual[i]=0.0;
1592  for (int j=0;j<nDOF_trial_element;j++)
1593  {
1594  elementTransport[i][j]=0.0;
1595  elementTransposeTransport[i][j]=0.0;
1596  }
1597  }
1598  //loop over quadrature points and compute integrands
1599  for (int k=0;k<nQuadraturePoints_element;k++)
1600  {
1601  //compute indeces and declare local storage
1602  register int eN_k = eN*nQuadraturePoints_element+k,
1603  eN_k_nSpace = eN_k*nSpace,
1604  eN_nDOF_trial_element = eN*nDOF_trial_element;
1605  register double
1606  // for entropy residual
1607  aux_entropy_residual=0., DENTROPY_un, DENTROPY_uni,
1608  //for mass matrix contributions
1609  u=0.0, un=0.0, grad_un[nSpace], porosity_times_velocity[nSpace],
1610  u_test_dV[nDOF_trial_element],
1611  u_grad_trial[nDOF_trial_element*nSpace],
1612  u_grad_test_dV[nDOF_test_element*nSpace],
1613  //for general use
1614  jac[nSpace*nSpace], jacDet, jacInv[nSpace*nSpace],
1615  dV,x,y,z,xt,yt,zt,
1616  //VRANS
1617  porosity;
1618  //get the physical integration weight
1619  ck.calculateMapping_element(eN,
1620  k,
1621  mesh_dof.data(),
1622  mesh_l2g.data(),
1623  mesh_trial_ref.data(),
1624  mesh_grad_trial_ref.data(),
1625  jac,
1626  jacDet,
1627  jacInv,
1628  x,y,z);
1629  ck.calculateMappingVelocity_element(eN,
1630  k,
1631  mesh_velocity_dof.data(),
1632  mesh_l2g.data(),
1633  mesh_trial_ref.data(),
1634  xt,yt,zt);
1635  dV = fabs(jacDet)*dV_ref[k];
1636  //get the solution (of Newton's solver). To compute time derivative term
1637  ck.valFromDOF(u_dof.data(),&u_l2g[eN_nDOF_trial_element],&u_trial_ref[k*nDOF_trial_element],u);
1638  //get the solution at quad point at tn and tnm1 for entropy viscosity
1639  ck.valFromDOF(u_dof_old.data(),&u_l2g[eN_nDOF_trial_element],&u_trial_ref[k*nDOF_trial_element],un);
1640  //get the solution gradients at tn for entropy viscosity
1641  ck.gradTrialFromRef(&u_grad_trial_ref[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
1642  ck.gradFromDOF(u_dof_old.data(),&u_l2g[eN_nDOF_trial_element],u_grad_trial,grad_un);
1643 
1644  //precalculate test function products with integration weights for mass matrix terms
1645  for (int j=0;j<nDOF_trial_element;j++)
1646  {
1647  u_test_dV[j] = u_test_ref[k*nDOF_trial_element+j]*dV;
1648  for (int I=0;I<nSpace;I++)
1649  u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
1650  }
1651 
1652  //calculate time derivative at quadrature points
1653  if (q_dV_last[eN_k] <= -100)
1654  q_dV_last[eN_k] = dV;
1655  q_dV[eN_k] = dV;
1656  //VRANS
1657  porosity = 1.-q_vos[eN_k];
1658  //
1659  //moving mesh
1660  //
1661  double mesh_velocity[3];
1662  mesh_velocity[0] = xt;
1663  mesh_velocity[1] = yt;
1664  mesh_velocity[2] = zt;
1665  //relative velocity at tn
1666  for (int I=0;I<nSpace;I++)
1667  porosity_times_velocity[I] = porosity*(velocity[eN_k_nSpace+I]-MOVING_DOMAIN*mesh_velocity[I]);
1668 
1670  // CALCULATE CELL BASED CFL //
1672  calculateCFL(elementDiameter[eN]/degree_polynomial,porosity_times_velocity,cfl[eN_k]);
1673 
1675  // CALCULATE ENTROPY RESIDUAL AT QUAD POINT //
1677  if (STABILIZATION_TYPE==2) // EV stab
1678  {
1679  for (int I=0;I<nSpace;I++)
1680  aux_entropy_residual += porosity_times_velocity[I]*grad_un[I];
1681  DENTROPY_un = ENTROPY_TYPE==0 ? DENTROPY(porosity*un,uL,uR) : DENTROPY_LOG(porosity*un,uL,uR);
1682  }
1684  // ith-LOOP //
1686  for(int i=0;i<nDOF_test_element;i++)
1687  {
1688  // VECTOR OF ENTROPY RESIDUAL //
1689  int eN_i=eN*nDOF_test_element+i;
1690  if (STABILIZATION_TYPE==2) // EV stab
1691  {
1692  int gi = offset_u+stride_u*u_l2g[eN_i]; //global i-th index
1693  double porosity_times_uni = (1.-vos_dof[gi])*u_dof_old[gi];
1694  DENTROPY_uni = ENTROPY_TYPE == 0 ? DENTROPY(porosity_times_uni,uL,uR) : DENTROPY_LOG(porosity_times_uni,uL,uR);
1695  element_entropy_residual[i] += (DENTROPY_un - DENTROPY_uni)*aux_entropy_residual*u_test_dV[i];
1696  }
1697  elementResidual_u[i] += porosity*(u-un)*u_test_dV[i];
1699  // j-th LOOP // To construct transport matrices
1701  for(int j=0;j<nDOF_trial_element;j++)
1702  {
1703  int j_nSpace = j*nSpace;
1704  int i_nSpace = i*nSpace;
1705  elementTransport[i][j] += // -int[(vel.grad_wi)*wj*dx]
1706  ck.AdvectionJacobian_weak(porosity_times_velocity,
1707  u_trial_ref[k*nDOF_trial_element+j],&u_grad_test_dV[i_nSpace]);
1708  elementTransposeTransport[i][j] += // -int[(vel.grad_wj)*wi*dx]
1709  ck.AdvectionJacobian_weak(porosity_times_velocity,
1710  u_trial_ref[k*nDOF_trial_element+i],&u_grad_test_dV[j_nSpace]);
1711  }
1712  }//i
1713  //save solution for other models
1714  q_u[eN_k] = u;
1715  q_m[eN_k] = porosity*u;
1716  }
1718  // DISTRIBUTE // load cell based element into global residual
1720  for(int i=0;i<nDOF_test_element;i++)
1721  {
1722  int eN_i=eN*nDOF_test_element+i;
1723  int gi = offset_u+stride_u*u_l2g[eN_i]; //global i-th index
1724 
1725  // distribute global residual for (lumped) mass matrix
1726  globalResidual[gi] += elementResidual_u[i];
1727  // distribute entropy_residual
1728  if (STABILIZATION_TYPE==2) // EV Stab
1729  global_entropy_residual[gi] += element_entropy_residual[i];
1730 
1731  // distribute transport matrices
1732  for (int j=0;j<nDOF_trial_element;j++)
1733  {
1734  int eN_i_j = eN_i*nDOF_trial_element+j;
1735  TransportMatrix[csrRowIndeces_CellLoops[eN_i] +
1736  csrColumnOffsets_CellLoops[eN_i_j]] += elementTransport[i][j];
1737  TransposeTransportMatrix[csrRowIndeces_CellLoops[eN_i] +
1738  csrColumnOffsets_CellLoops[eN_i_j]]
1739  += elementTransposeTransport[i][j];
1740  }//j
1741  }//i
1742  }//elements
1743 
1745  // ADD OUTFLOW BOUNDARY TERM TO TRANSPORT MATRICES AND COMPUTE INFLOW BOUNDARY INTEGRAL //
1747  // * Compute outflow boundary integral as a matrix; i.e., int_B[ (vel.normal)*wi*wj*dx]
1748  for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
1749  {
1750  double min_u_bc_local = 1E10, max_u_bc_local = -1E10;
1751  register int ebN = exteriorElementBoundariesArray[ebNE];
1752  register int eN = elementBoundaryElementsArray[ebN*2+0],
1753  ebN_local = elementBoundaryLocalElementBoundariesArray[ebN*2+0],
1754  eN_nDOF_trial_element = eN*nDOF_trial_element;
1755  register double elementResidual_u[nDOF_test_element];
1756  for (int i=0;i<nDOF_test_element;i++)
1757  elementResidual_u[i]=0.0;
1758  // loop on quad points
1759  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
1760  {
1761  register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
1762  ebNE_kb_nSpace = ebNE_kb*nSpace,
1763  ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
1764  ebN_local_kb_nSpace = ebN_local_kb*nSpace;
1765  register double
1766  u_ext=0.0, bc_u_ext=0.0,
1767  porosity_times_velocity[nSpace],
1768  flux_ext=0.0, dflux_ext=0.0,
1769  fluxTransport[nDOF_trial_element],
1770  jac_ext[nSpace*nSpace],
1771  jacDet_ext,
1772  jacInv_ext[nSpace*nSpace],
1773  boundaryJac[nSpace*(nSpace-1)],
1774  metricTensor[(nSpace-1)*(nSpace-1)],
1775  metricTensorDetSqrt,
1776  dS,
1777  u_test_dS[nDOF_test_element],
1778  normal[nSpace],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,porosity_ext;
1779  // calculate mappings
1780  ck.calculateMapping_elementBoundary(eN,
1781  ebN_local,
1782  kb,
1783  ebN_local_kb,
1784  mesh_dof.data(),
1785  mesh_l2g.data(),
1786  mesh_trial_trace_ref.data(),
1787  mesh_grad_trial_trace_ref.data(),
1788  boundaryJac_ref.data(),
1789  jac_ext,
1790  jacDet_ext,
1791  jacInv_ext,
1792  boundaryJac,
1793  metricTensor,
1794  metricTensorDetSqrt,
1795  normal_ref.data(),
1796  normal,
1797  x_ext,y_ext,z_ext);
1798  ck.calculateMappingVelocity_elementBoundary(eN,
1799  ebN_local,
1800  kb,
1801  ebN_local_kb,
1802  mesh_velocity_dof.data(),
1803  mesh_l2g.data(),
1804  mesh_trial_trace_ref.data(),
1805  xt_ext,yt_ext,zt_ext,
1806  normal,
1807  boundaryJac,
1808  metricTensor,
1809  integralScaling);
1810  dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref[kb];
1811  //compute shape and solution information
1812  ck.valFromDOF(u_dof.data(),&u_l2g[eN_nDOF_trial_element],&u_trial_trace_ref[ebN_local_kb*nDOF_test_element],u_ext);
1813  //precalculate test function products with integration weights
1814  for (int j=0;j<nDOF_trial_element;j++)
1815  u_test_dS[j] = u_test_trace_ref[ebN_local_kb*nDOF_test_element+j]*dS;
1816 
1817  //VRANS
1818  porosity_ext = 1.0-ebqe_vos_ext[ebNE_kb];
1819  //
1820  //moving mesh
1821  //
1822  double mesh_velocity[3];
1823  mesh_velocity[0] = xt_ext;
1824  mesh_velocity[1] = yt_ext;
1825  mesh_velocity[2] = zt_ext;
1826  //std::cout<<"mesh_velocity ext"<<std::endl;
1827  for (int I=0;I<nSpace;I++)
1828  porosity_times_velocity[I] = porosity_ext*(ebqe_velocity_ext[ebNE_kb_nSpace+I] - MOVING_DOMAIN*mesh_velocity[I]);
1829  //
1830  //calculate the fluxes
1831  //
1832  double flow = 0.;
1833  for (int I=0; I < nSpace; I++)
1834  flow += normal[I]*porosity_times_velocity[I];
1835 
1836  if (flow >= 0) //outflow. This is handled via the transport matrices. Then flux_ext=0 and dflux_ext!=0
1837  {
1838  dflux_ext = flow;
1839  flux_ext = 0;
1840  // save external u
1841  ebqe_u[ebNE_kb] = u_ext;
1842  }
1843  else // inflow. This is handled via the boundary integral. Then flux_ext!=0 and dflux_ext=0
1844  {
1845  dflux_ext = 0;
1846  // save external u
1847  ebqe_u[ebNE_kb] = isDOFBoundary_u[ebNE_kb]*ebqe_bc_u_ext[ebNE_kb]+(1-isDOFBoundary_u[ebNE_kb])*u_ext;
1848  if (isDOFBoundary_u[ebNE_kb] == 1)
1849  flux_ext = ebqe_bc_u_ext[ebNE_kb]*flow;
1850  else if (isFluxBoundary_u[ebNE_kb] == 1)
1851  flux_ext = ebqe_bc_flux_u_ext[ebNE_kb];
1852  else
1853  {
1854  std::cout<<"warning: VOF open boundary with no external trace, setting to zero for inflow"<<std::endl;
1855  flux_ext = 0.0;
1856  }
1857  }
1858 
1859  for (int j=0;j<nDOF_trial_element;j++)
1860  {
1861  // elementResidual. This is to include the inflow boundary integral.
1862  // NOTE: here I assume that we use a Galerkin approach st nDOF_test_element = nDOF_trial_element
1863  elementResidual_u[j] += flux_ext*u_test_dS[j];
1864  register int ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j;
1865  fluxTransport[j] = dflux_ext*u_trial_trace_ref[ebN_local_kb_j];
1866  }//j
1868  // DISTRIBUTE OUTFLOW BOUNDARY TO TRANSPORT MATRICES //
1870  for (int i=0;i<nDOF_test_element;i++)
1871  {
1872  register int eN_i = eN*nDOF_test_element+i;
1873  for (int j=0;j<nDOF_trial_element;j++)
1874  {
1875  register int ebN_i_j = ebN*4*nDOF_test_X_trial_element + i*nDOF_trial_element + j;
1876  TransportMatrix[csrRowIndeces_CellLoops[eN_i] + csrColumnOffsets_eb_CellLoops[ebN_i_j]]
1877  += fluxTransport[j]*u_test_dS[i];
1878  TransposeTransportMatrix[csrRowIndeces_CellLoops[eN_i] + csrColumnOffsets_eb_CellLoops[ebN_i_j]]
1879  += fluxTransport[i]*u_test_dS[j];
1880  }//j
1881  }//i
1882  // local min/max at boundary
1883  min_u_bc_local = fmin(ebqe_u[ebNE_kb], min_u_bc_local);
1884  max_u_bc_local = fmax(ebqe_u[ebNE_kb], max_u_bc_local);
1885  }//kb
1886  // global min/max at boundary
1887  for (int i=0;i<nDOF_test_element;i++)
1888  {
1889  int eN_i = eN*nDOF_test_element+i;
1890  int gi = offset_u+stride_u*u_l2g[eN_i]; //global i-th index
1891  globalResidual[gi] += dt*elementResidual_u[i];
1892  boundary_integral[gi] += elementResidual_u[i];
1893  min_u_bc[gi] = fmin(min_u_bc_local,min_u_bc[gi]);
1894  max_u_bc[gi] = fmax(max_u_bc_local,max_u_bc[gi]);
1895  }
1896  }//ebNE
1897  // END OF ADDING BOUNDARY TERM TO TRANSPORT MATRICES and COMPUTING BOUNDARY INTEGRAL //
1898 
1900  // COMPUTE SMOOTHNESS INDICATOR and NORMALIZE ENTROPY RESIDUAL //
1902  // NOTE: see NCLS.h for a different but equivalent implementation of this.
1903  int ij = 0;
1904  for (int i=0; i<numDOFs; i++)
1905  {
1906  double etaMaxi, etaMini;
1907  if (STABILIZATION_TYPE==2) //EV
1908  {
1909  // For eta min and max
1910  etaMaxi = fabs(eta[i]);
1911  etaMini = fabs(eta[i]);
1912  }
1913  double porosity_times_solni = (1.-vos_dof[i])*u_dof_old[i];
1914  // for smoothness indicator //
1915  double alpha_numerator = 0., alpha_denominator = 0.;
1916  for (int offset=csrRowIndeces_DofLoops[i]; offset<csrRowIndeces_DofLoops[i+1]; offset++)
1917  { // First loop in j (sparsity pattern)
1918  int j = csrColumnOffsets_DofLoops[offset];
1919  if (STABILIZATION_TYPE==2) //EV Stabilization
1920  {
1921  // COMPUTE ETA MIN AND ETA MAX //
1922  etaMaxi = fmax(etaMaxi,fabs(eta[j]));
1923  etaMini = fmin(etaMini,fabs(eta[j]));
1924  }
1925  double porosity_times_solnj = (1.-vos_dof[i])*u_dof_old[j];
1926  alpha_numerator += porosity_times_solni - porosity_times_solnj;
1927  alpha_denominator += fabs(porosity_times_solni - porosity_times_solnj);
1928  //update ij
1929  ij+=1;
1930  }
1931  if (STABILIZATION_TYPE==2) //EV Stab
1932  {
1933  // Normalize entropy residual
1934  global_entropy_residual[i] *= etaMini == etaMaxi ? 0. : 2*cE/(etaMaxi-etaMini);
1935  quantDOFs[i] = fabs(global_entropy_residual[i]);
1936  }
1937 
1938  double alphai = alpha_numerator/(alpha_denominator+1E-15);
1939  quantDOFs[i] = alphai;
1940 
1941  if (POWER_SMOOTHNESS_INDICATOR==0)
1942  psi[i] = 1.0;
1943  else
1944  psi[i] = std::pow(alphai,POWER_SMOOTHNESS_INDICATOR); //NOTE: they use alpha^2 in the paper
1945  }
1947  // ** LOOP IN DOFs FOR EDGE BASED TERMS ** //
1949  ij=0;
1950  for (int i=0; i<numDOFs; i++)
1951  {
1952  // NOTE: Transport matrices already have the porosity considered. ---> Dissipation matrices as well.
1953  double solni = u_dof_old[i]; // solution at time tn for the ith DOF
1954  double porosityi = 1.-vos_dof[i];
1955  double ith_dissipative_term = 0;
1956  double ith_low_order_dissipative_term = 0;
1957  double ith_flux_term = 0;
1958  double dLii = 0.;
1959 
1960  // loop over the sparsity pattern of the i-th DOF
1961  for (int offset=csrRowIndeces_DofLoops[i]; offset<csrRowIndeces_DofLoops[i+1]; offset++)
1962  {
1963  int j = csrColumnOffsets_DofLoops[offset];
1964  double solnj = u_dof_old[j]; // solution at time tn for the jth DOF
1965  double porosityj = 1.-vos_dof[j];
1966  double dLowij, dLij, dEVij, dHij;
1967 
1968  ith_flux_term += TransportMatrix[ij]*solnj;
1969  if (i != j) //NOTE: there is really no need to check for i!=j (see formula for ith_dissipative_term)
1970  {
1971  // artificial compression
1972  double solij = 0.5*(porosityi*solni+porosityj*solnj);
1973  double Compij = cK*fmax(solij*(1.0-solij),0.0)/(fabs(porosityi*solni-porosityj*solnj)+1E-14);
1974  // first-order dissipative operator
1975  dLowij = fmax(fabs(TransportMatrix[ij]),fabs(TransposeTransportMatrix[ij]));
1976  //dLij = fmax(0.,fmax(psi[i]*TransportMatrix[ij], // Approach by S. Badia
1977  // psi[j]*TransposeTransportMatrix[ij]));
1978  dLij = dLowij*fmax(psi[i],psi[j]); // Approach by JLG & BP
1979  if (STABILIZATION_TYPE==2) //EV Stab
1980  {
1981  // high-order (entropy viscosity) dissipative operator
1982  dEVij = fmax(fabs(global_entropy_residual[i]),fabs(global_entropy_residual[j]));
1983  dHij = fmin(dLowij,dEVij) * fmax(1.0-Compij,0.0); // artificial compression
1984  }
1985  else // smoothness based indicator
1986  {
1987  dHij = dLij * fmax(1.0-Compij,0.0); // artificial compression
1988  }
1989  //dissipative terms
1990  ith_dissipative_term += dHij*(solnj-solni);
1991  ith_low_order_dissipative_term += dLowij*(solnj-solni);
1992  //dHij - dLij. This matrix is needed during FCT step
1993  dt_times_dH_minus_dL[ij] = dt*(dHij - dLowij);
1994  dLii -= dLij;
1995  dLow[ij] = dLowij;
1996  }
1997  else //i==j
1998  {
1999  // NOTE: this is incorrect. Indeed, dLii = -sum_{j!=i}(dLij) and similarly for dCii.
2000  // However, it is irrelevant since during the FCT step we do (dL-dC)*(solnj-solni)
2001  dt_times_dH_minus_dL[ij]=0;
2002  dLow[ij]=0;
2003  }
2004  //update ij
2005  ij+=1;
2006  }
2007  double mi = ML[i];
2008  // compute edge_based_cfl
2009  edge_based_cfl[i] = 2.*fabs(dLii)/mi;
2010  uLow[i] = u_dof_old[i] - dt/mi*(ith_flux_term
2011  + boundary_integral[i]
2012  - ith_low_order_dissipative_term);
2013 
2014  // update residual
2015  if (LUMPED_MASS_MATRIX==1)
2016  globalResidual[i] = u_dof_old[i] - dt/mi*(ith_flux_term
2017  + boundary_integral[i]
2018  - ith_dissipative_term);
2019  else
2020  globalResidual[i] += dt*(ith_flux_term - ith_dissipative_term);
2021  }
2022  }
2023  };//VOF
2024 
2025  inline cppVOF3P_base* newVOF3P(int nSpaceIn,
2026  int nQuadraturePoints_elementIn,
2027  int nDOF_mesh_trial_elementIn,
2028  int nDOF_trial_elementIn,
2029  int nDOF_test_elementIn,
2030  int nQuadraturePoints_elementBoundaryIn,
2031  int CompKernelFlag)
2032  {
2033  if (nSpaceIn == 2)
2034  return proteus::chooseAndAllocateDiscretization2D<cppVOF3P_base,cppVOF,CompKernel>(nSpaceIn,
2035  nQuadraturePoints_elementIn,
2036  nDOF_mesh_trial_elementIn,
2037  nDOF_trial_elementIn,
2038  nDOF_test_elementIn,
2039  nQuadraturePoints_elementBoundaryIn,
2040  CompKernelFlag);
2041  else
2042  return proteus::chooseAndAllocateDiscretization<cppVOF3P_base,cppVOF,CompKernel>(nSpaceIn,
2043  nQuadraturePoints_elementIn,
2044  nDOF_mesh_trial_elementIn,
2045  nDOF_trial_elementIn,
2046  nDOF_test_elementIn,
2047  nQuadraturePoints_elementBoundaryIn,
2048  CompKernelFlag);
2049  }
2050 }//proteus
2051 #endif
proteus::cppVOF::calculateSubgridError_tau
void calculateSubgridError_tau(const double &elementDiameter, const double &dmt, const double dH[nSpace], double &cfl, double &tau)
Definition: VOF3P.h:108
proteus::cppVOF::cppVOF
cppVOF()
Definition: VOF3P.h:70
proteus::cppVOF3P_base::FluxCorrectionMatrix
std::valarray< double > FluxCorrectionMatrix
Definition: VOF3P.h:48
proteus::newVOF3P
cppVOF3P_base * newVOF3P(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition: VOF3P.h:2025
GLOBAL_FCT
#define GLOBAL_FCT
Definition: VOF3P.h:13
proteus::cppVOF3P_base::TransportMatrix
std::valarray< double > TransportMatrix
Definition: VOF3P.h:49
proteus::cppVOF3P_base::global_entropy_residual
std::valarray< double > global_entropy_residual
Definition: VOF3P.h:50
proteus::cppVOF3P_base::psi
std::valarray< double > psi
Definition: VOF3P.h:50
proteus::cppVOF::nDOF_test_X_trial_element
const int nDOF_test_X_trial_element
Definition: VOF3P.h:68
proteus::cppVOF3P_base::FCTStep
virtual void FCTStep(arguments_dict &args)=0
proteus::cppVOF3P_base::~cppVOF3P_base
virtual ~cppVOF3P_base()
Definition: VOF3P.h:51
proteus::cppVOF3P_base::calculateResidualEdgeBased
virtual void calculateResidualEdgeBased(arguments_dict &args)=0
proteus::DENTROPY
double DENTROPY(const double &phi, const double &dummyL, const double &dummyR)
Definition: NCLS.h:25
n
Int n
Definition: Headers.h:28
df
double df(double C, double b, double a, int q, int r)
Definition: analyticalSolutions.c:2209
ModelFactory.h
CompKernel.h
proteus::arguments_dict::scalar
T & scalar(const std::string &key)
proteus::arguments_dict::array
xt::pyarray< T > & array(const std::string &key)
proteus::cppVOF3P_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)=0
num
Int num
Definition: Headers.h:32
H
Double H
Definition: Headers.h:65
proteus::cppVOF::ck
CompKernelType ck
Definition: VOF3P.h:69
proteus::cppVOF::calculateJacobian
void calculateJacobian(arguments_dict &args)
Definition: VOF3P.h:884
proteus::cppVOF
Definition: VOF3P.h:66
proteus::cppVOF::calculateNumericalDiffusion
void calculateNumericalDiffusion(const double &shockCapturingDiffusion, const double &elementDiameter, const double &strong_residual, const double grad_u[nSpace], double &numDiff)
Definition: VOF3P.h:142
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
proteus::DENTROPY_LOG
double DENTROPY_LOG(const double &phi, const double &phiL, const double &phiR)
Definition: NCLS.h:32
proteus::cppVOF::evaluateCoefficients
void evaluateCoefficients(const double v[nSpace], const double &u, const double &porosity, double &m, double &dm, double f[nSpace], double df[nSpace])
Definition: VOF3P.h:90
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::cppVOF3P_base::Rneg
std::valarray< double > Rneg
Definition: VOF3P.h:47
proteus::cppVOF::FCTStep
void FCTStep(arguments_dict &args)
Definition: VOF3P.h:1342
proteus::ENTROPY_LOG
double ENTROPY_LOG(const double &phi, const double &phiL, const double &phiR)
Definition: NCLS.h:29
xt
Definition: AddedMass.cpp:7
cMax
#define cMax
Definition: NCLS3P.h:11
proteus::cppVOF3P_base::TransposeTransportMatrix
std::valarray< double > TransposeTransportMatrix
Definition: VOF3P.h:49
proteus::cppVOF::calculateResidualElementBased
void calculateResidualElementBased(arguments_dict &args)
Definition: VOF3P.h:247
proteus::cppVOF::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: VOF3P.h:126
proteus
Definition: ADR.h:17
proteus::cppVOF3P_base
Definition: VOF3P.h:44
proteus::f
double f(const double &g, const double &h, const double &hZ)
Definition: SW2DCV.h:58
proteus::cppVOF::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: VOF3P.h:162
proteus::cppVOF3P_base::boundary_integral
std::valarray< double > boundary_integral
Definition: VOF3P.h:50
proteus::cppVOF::exteriorNumericalAdvectiveFluxDerivative
void exteriorNumericalAdvectiveFluxDerivative(const int &isDOFBoundary_u, const int &isFluxBoundary_u, const double n[nSpace], const double velocity[nSpace], double &dflux)
Definition: VOF3P.h:212
proteus::arguments_dict
Definition: ArgumentsDict.h:70
proteus::cppVOF3P_base::calculateResidualElementBased
virtual void calculateResidualElementBased(arguments_dict &args)=0
proteus::cppVOF3P_base::eta
std::valarray< double > eta
Definition: VOF3P.h:50
proteus::cppVOF::calculateResidualEdgeBased
void calculateResidualEdgeBased(arguments_dict &args)
Definition: VOF3P.h:1456
ArgumentsDict.h
proteus::cppVOF::calculateCFL
void calculateCFL(const double &elementDiameter, const double df[nSpace], double &cfl)
Definition: VOF3P.h:76
proteus::cppVOF3P_base::Rpos
std::valarray< double > Rpos
Definition: VOF3P.h:47
POWER_SMOOTHNESS_INDICATOR
#define POWER_SMOOTHNESS_INDICATOR
Definition: VOF3P.h:11
cE
#define cE
Definition: NCLS3P.h:10
Generated on Fri Jul 1 2022 11:17:46 for proteus by doxygen 1.8.20