proteus  1.8.1
C/C++/Fortran libraries
CLSVOF.h
Go to the documentation of this file.
1 #ifndef CLSVOF_H
2 #define CLSVOF_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 USE_SIGN_FUNCTION 1
14 #define IMPLICIT_BCs 0
15 #define LAMBDA_SCALING 0
16 
17 namespace proteus
18 {
19 // True characteristic functions
20  inline double heaviside(const double& z){
21  return (z>0 ? 1. : (z<0 ? 0. : 0.5));
22  }
23  inline double Sign(const double& z){
24  return (z>0 ? 1. : (z<0 ? -1. : 0.));
25  }
26 }
27 
28 namespace proteus
29 {
30  class CLSVOF_base
31  {
32  //The base class defining the interface
33  public:
34  std::valarray<double> Rpos, Rneg, FluxCorrectionMatrix;
35  virtual ~CLSVOF_base(){}
36  virtual void calculateResidual(arguments_dict& args)=0;
37  virtual void calculateJacobian(arguments_dict& args)=0;
38  virtual std::tuple<double, double, double, double, double, double, double, double, double, double, double>
39  calculateMetricsAtEOS(arguments_dict& args)=0;
40  virtual std::tuple<double, double, double, double, double> calculateMetricsAtETS(arguments_dict& args)=0;
41  virtual void normalReconstruction(arguments_dict& args)=0;
42  virtual void calculateRhsL2Proj(arguments_dict& args)=0;
43  virtual void calculateLumpedMassMatrix(arguments_dict& args)=0;
44  virtual void assembleSpinUpSystem(arguments_dict& args)=0;
45  virtual void FCTStep(arguments_dict& args)=0;
46  };
47 
48  template<class CompKernelType,
49  int nSpace,
50  int nQuadraturePoints_element,
51  int nDOF_mesh_trial_element,
52  int nDOF_trial_element,
53  int nDOF_test_element,
54  int nQuadraturePoints_elementBoundary>
55  class CLSVOF : public CLSVOF_base
56  {
57  public:
58  const int nDOF_test_X_trial_element;
59  CompKernelType ck;
60  CLSVOF():
61  nDOF_test_X_trial_element(nDOF_test_element*nDOF_trial_element),
62  ck()
63  {}
64 
65  inline
66  void calculateCFL(const double& elementDiameter,
67  const double df[nSpace],
68  double& cfl)
69  {
70  double h,nrm_v;
71  h = elementDiameter;
72  nrm_v=0.0;
73  for(int I=0;I<nSpace;I++)
74  nrm_v+=df[I]*df[I];
75  nrm_v = sqrt(nrm_v);
76  cfl = nrm_v/h;
77  }
78 
79  inline void calculateNonlinearCFL(const double& elementDiameter,
80  const double df[nSpace],
81  const double norm_factor_lagged,
82  const double epsFactHeaviside,
83  const double lambdaFact,
84  double& cfl)
85  {
86  double h,nrm2_v;
87  h = elementDiameter;
88  nrm2_v=0.0;
89  for(int I=0;I<nSpace;I++)
90  nrm2_v+=df[I]*df[I];
91  cfl = nrm2_v*norm_factor_lagged/(epsFactHeaviside*lambdaFact*h*h*h);
92  cfl = std::sqrt(cfl);
93  }
94 
95  inline void exteriorNumericalAdvectiveFlux(const int& isDOFBoundary_u,
96  const int& isFluxBoundary_u,
97  const double n[nSpace],
98  const double& bc_u,
99  const double& bc_flux_u,
100  const double& u,
101  const double velocity[nSpace],
102  double& flux)
103  {
104  double flow=0.0;
105  for (int I=0; I < nSpace; I++)
106  flow += n[I]*velocity[I];
107  if (isDOFBoundary_u == 1)
108  {
109  if (flow >= 0.0)
110  {
111  flux = u*flow;
112  }
113  else
114  {
115  flux = bc_u*flow;
116  }
117  }
118  else if (isFluxBoundary_u == 1)
119  {
120  flux = bc_flux_u;
121  }
122  else
123  {
124  if (flow >= 0.0)
125  {
126  flux = u*flow;
127  }
128  else
129  {
130  std::cout<<"warning: CLSVOF open boundary with no external trace, setting to zero for inflow"<<std::endl;
131  flux = 0.0;
132  }
133 
134  }
135  }
136 
137  inline void exteriorNumericalAdvectiveFluxDerivative(const int& isDOFBoundary_u,
138  const int& isFluxBoundary_u,
139  const double n[nSpace],
140  const double& du,
141  const double velocity[nSpace],
142  double& dflux)
143  {
144  double flow=0.0;
145  for (int I=0; I < nSpace; I++)
146  flow += n[I]*velocity[I];
147  dflux=0.0;//default to no flux
148  if (isDOFBoundary_u == 1)
149  {
150  if (flow >= 0.0)
151  {
152  dflux = du*flow;
153  }
154  else
155  {
156  dflux = 0.0; //zero since inflow BC is given by data (so independent on soln)
157  }
158  }
159  else if (isFluxBoundary_u == 1)
160  {
161  dflux = 0.0;
162  }
163  else
164  {
165  if (flow >= 0.0)
166  {
167  dflux = flow;
168  }
169  else
170  {
171  dflux = 0.;
172  }
173  }
174  }
175 
176  inline double smoothedHeaviside(double eps, double u)
177  {
178  double H;
179  if (u > eps)
180  H=1.0;
181  else if (u < -eps)
182  H=0.0;
183  else if (u==0.0)
184  H=0.5;
185  else
186  H = 0.5*(1.0 + u/eps + sin(M_PI*u/eps)/M_PI);
187  return H;
188  }
189 
190  inline double smoothedDirac(double eps, double u)
191  {
192  double d;
193  if (u > eps)
194  d=0.0;
195  else if (u < -eps)
196  d=0.0;
197  else
198  d = 0.5*(1.0 + cos(M_PI*u/eps))/eps;
199  return d;
200  }
201 
202  inline double smoothedNormalizedDirac(double eps, double u)
203  {
204  double d;
205  if (u > eps)
206  d=0.0;
207  else if (u < -eps)
208  d=0.0;
209  else
210  d = 0.5*(1.0 + cos(M_PI*u/eps));
211  return d;
212  }
213 
214  inline double smoothedSign(double eps, double u)
215  {
216  double H;
217  if (u > eps)
218  H=1.0;
219  else if (u < -eps)
220  H=0.0;
221  else if (u==0.0)
222  H=0.5;
223  else
224  H = 0.5*(1.0 + u/eps + sin(M_PI*u/eps)/M_PI);
225  if (USE_SIGN_FUNCTION==1)
226  return 2*H-1;
227  else
228  return H;
229  }
230 
231  inline double smoothedDerivativeSign(double eps, double u)
232  {
233  double d;
234  if (u > eps)
235  d=0.0;
236  else if (u < -eps)
237  d=0.0;
238  else
239  d = 0.5*(1.0 + cos(M_PI*u/eps))/eps;
240  if (USE_SIGN_FUNCTION==1)
241  return 2*d;
242  else
243  return d;
244  }
245 
246  void calculateResidual(arguments_dict& args)
247  {
248  double dt = args.scalar<double>("dt");
249  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
250  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
251  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
252  xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
253  double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
254  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
255  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
256  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
257  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
258  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
259  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
260  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
261  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
262  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
263  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
264  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
265  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
266  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
267  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
268  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
269  int nElements_global = args.scalar<int>("nElements_global");
270  int nElements_owned = args.scalar<int>("nElements_owned");
271  double useMetrics = args.scalar<double>("useMetrics");
272  const xt::pyarray<double>& q_vos = args.array<double>("q_vos");
273  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
274  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
275  xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
276  int degree_polynomial = args.scalar<int>("degree_polynomial");
277  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
278  xt::pyarray<double>& u_dof_old = args.array<double>("u_dof_old");
279  xt::pyarray<double>& velocity = args.array<double>("velocity");
280  xt::pyarray<double>& velocity_old = args.array<double>("velocity_old");
281  xt::pyarray<double>& q_m = args.array<double>("q_m");
282  xt::pyarray<double>& q_u = args.array<double>("q_u");
283  xt::pyarray<double>& q_n = args.array<double>("q_n");
284  xt::pyarray<double>& q_H = args.array<double>("q_H");
285  xt::pyarray<double>& q_mH = args.array<double>("q_mH");
286  xt::pyarray<double>& q_dV = args.array<double>("q_dV");
287  xt::pyarray<double>& q_dV_last = args.array<double>("q_dV_last");
288  xt::pyarray<double>& cfl = args.array<double>("cfl");
289  int offset_u = args.scalar<int>("offset_u");
290  int stride_u = args.scalar<int>("stride_u");
291  xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
292  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
293  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
294  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
295  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
296  xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
297  const xt::pyarray<double>& ebqe_vos_ext = args.array<double>("ebqe_vos_ext");
298  xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
299  xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
300  xt::pyarray<int>& isFluxBoundary_u = args.array<int>("isFluxBoundary_u");
301  xt::pyarray<double>& ebqe_bc_flux_u_ext = args.array<double>("ebqe_bc_flux_u_ext");
302  xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
303  xt::pyarray<double>& ebqe_n = args.array<double>("ebqe_n");
304  xt::pyarray<double>& ebqe_H = args.array<double>("ebqe_H");
305  xt::pyarray<double>& ebqe_flux = args.array<double>("ebqe_flux");
306  int timeOrder = args.scalar<int>("timeOrder");
307  int timeStage = args.scalar<int>("timeStage");
308  double epsFactHeaviside = args.scalar<double>("epsFactHeaviside");
309  double epsFactDirac = args.scalar<double>("epsFactDirac");
310  double epsFactRedist = args.scalar<double>("epsFactRedist");
311  double lambdaFact = args.scalar<double>("lambdaFact");
312  xt::pyarray<double>& min_distance = args.array<double>("min_distance");
313  xt::pyarray<double>& max_distance = args.array<double>("max_distance");
314  xt::pyarray<double>& mean_distance = args.array<double>("mean_distance");
315  xt::pyarray<double>& volume_domain = args.array<double>("volume_domain");
316  double norm_factor_lagged = args.scalar<double>("norm_factor_lagged");
317  double VelMax = args.scalar<double>("VelMax");
318  xt::pyarray<double>& projected_qx_tn = args.array<double>("projected_qx_tn");
319  xt::pyarray<double>& projected_qy_tn = args.array<double>("projected_qy_tn");
320  xt::pyarray<double>& projected_qz_tn = args.array<double>("projected_qz_tn");
321  xt::pyarray<double>& projected_qx_tStar = args.array<double>("projected_qx_tStar");
322  xt::pyarray<double>& projected_qy_tStar = args.array<double>("projected_qy_tStar");
323  xt::pyarray<double>& projected_qz_tStar = args.array<double>("projected_qz_tStar");
324  int numDOFs = args.scalar<int>("numDOFs");
325  xt::pyarray<double>& lumped_mass_matrix = args.array<double>("lumped_mass_matrix");
326  xt::pyarray<double>& H_dof = args.array<double>("H_dof");
327  int preRedistancingStage = args.scalar<int>("preRedistancingStage");
328  xt::pyarray<double>& interface_locator = args.array<double>("interface_locator");
329  double alpha = args.scalar<double>("alpha");
330  min_distance.data()[0] = 1E10;
331  max_distance.data()[0] = -1E10;
332  mean_distance.data()[0] = 0.;
333  volume_domain.data()[0] = 0.;
334 
335  for(int eN=0;eN<nElements_global;eN++)
336  {
337  //declare local storage for local contributions and initialize
338  register double elementResidual_u[nDOF_test_element], element_rhs_L2proj_H[nDOF_test_element];
339  for (int i=0;i<nDOF_test_element;i++)
340  {
341  elementResidual_u[i]=0.0;
342  element_rhs_L2proj_H[i]=0.0;
343  }
344  //loop over quadrature points and compute integrands
345  for (int k=0;k<nQuadraturePoints_element;k++)
346  {
347  //compute indeces and declare local storage
348  register int eN_k = eN*nQuadraturePoints_element+k,
349  eN_k_nSpace = eN_k*nSpace,
350  eN_nDOF_trial_element = eN*nDOF_trial_element;
351  register double
352  //for mass matrix contributions
353  u, un, grad_u[nSpace], grad_un[nSpace],
354  qxn, qyn, qzn, //qxnStar, qynStar, qznStar,
355  normalReconstruction[3], // assume 3D always
356  u_test_dV[nDOF_trial_element],
357  u_grad_trial[nDOF_trial_element*nSpace],
358  u_grad_test_dV[nDOF_test_element*nSpace],
359  //for general use
360  jac[nSpace*nSpace], jacDet, jacInv[nSpace*nSpace],
361  dV,x,y,z,xt,yt,zt,h_phi,
362  porosity;
363  //get the physical integration weight
364  ck.calculateMapping_element(eN,
365  k,
366  mesh_dof.data(),
367  mesh_l2g.data(),
368  mesh_trial_ref.data(),
369  mesh_grad_trial_ref.data(),
370  jac,
371  jacDet,
372  jacInv,
373  x,y,z);
374  ck.calculateH_element(eN,
375  k,
376  nodeDiametersArray.data(),
377  mesh_l2g.data(),
378  mesh_trial_ref.data(),
379  h_phi);
380  ck.calculateMappingVelocity_element(eN,
381  k,
382  mesh_velocity_dof.data(),
383  mesh_l2g.data(),
384  mesh_trial_ref.data(),
385  xt,yt,zt);
386  dV = fabs(jacDet)*dV_ref.data()[k];
387  ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],
388  jacInv,
389  u_grad_trial);
390  // get the components of the normal reconstruction //
391  ck.valFromDOF(projected_qx_tn.data(),
392  &u_l2g.data()[eN_nDOF_trial_element],
393  &u_trial_ref.data()[k*nDOF_trial_element],
394  qxn);
395  ck.valFromDOF(projected_qy_tn.data(),
396  &u_l2g.data()[eN_nDOF_trial_element],
397  &u_trial_ref.data()[k*nDOF_trial_element],
398  qyn);
399  ck.valFromDOF(projected_qz_tn.data(),
400  &u_l2g.data()[eN_nDOF_trial_element],
401  &u_trial_ref.data()[k*nDOF_trial_element],
402  qzn);
403  //ck.valFromDOF(projected_qx_tStar.data(),
404  // &u_l2g.data()[eN_nDOF_trial_element],
405  // &u_trial_ref.data()[k*nDOF_trial_element],
406  // qxnStar);
407  //ck.valFromDOF(projected_qy_tStar.data(),
408  // &u_l2g.data()[eN_nDOF_trial_element],
409  // &u_trial_ref.data()[k*nDOF_trial_element],
410  // qynStar);
411  //ck.valFromDOF(projected_qz_tStar.data(),
412  // &u_l2g.data()[eN_nDOF_trial_element],
413  // &u_trial_ref.data()[k*nDOF_trial_element],
414  // qznStar);
415  // get the solution (of Newton's solver)
416  ck.valFromDOF(u_dof.data(),
417  &u_l2g.data()[eN_nDOF_trial_element],
418  &u_trial_ref.data()[k*nDOF_trial_element],
419  u);
420  // get old solution
421  ck.valFromDOF(u_dof_old.data(),
422  &u_l2g.data()[eN_nDOF_trial_element],
423  &u_trial_ref.data()[k*nDOF_trial_element],
424  un);
425  //get the solution gradients at quad points
426  ck.gradFromDOF(u_dof.data(),
427  &u_l2g.data()[eN_nDOF_trial_element],
428  u_grad_trial,
429  grad_u);
430  ck.gradFromDOF(u_dof_old.data(),
431  &u_l2g.data()[eN_nDOF_trial_element],
432  u_grad_trial,
433  grad_un);
434  //precalculate test function products with integration weights for mass matrix terms
435  for (int j=0;j<nDOF_trial_element;j++)
436  {
437  u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
438  for (int I=0;I<nSpace;I++)
439  u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;
440  }
441  double hK=(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN])/degree_polynomial;
442  //VRANS
443  porosity=1.0-q_vos.data()[eN_k];
444 
446  // NORMAL RECONSTRUCTION //
448  //if (timeOrder == 2 && timeStage == 2) // this is never reached
449  //{
450  //normalReconstruction[0] = 0.5*(qxnStar+qxn);
451  //normalReconstruction[1] = 0.5*(qynStar+qyn);
452  //if (nSpace==3)
453  //normalReconstruction[2] = 0.5*(qznStar+qzn);
454  //else
455  //normalReconstruction[2] = 0.;
456  //}
457  //else //timeOrder == 1 or timeStage==1
458  //{
459  normalReconstruction[0] = qxn;
460  normalReconstruction[1] = qyn;
461  if (nSpace==3)
462  normalReconstruction[2] = qzn;
463  else
464  normalReconstruction[2] = 0.;
465  //}
466 
468  // ADJUSTMENT ON dV DUE TO MESH MOTION //
470  if (q_dV_last.data()[eN_k] <= -100)
471  q_dV_last.data()[eN_k] = dV;
472  q_dV.data()[eN_k] = dV;
473 
474  double delta, residualEikonal, tau, backgroundDissipation=0.1*hK;
475  double time_derivative_residual, fnHalf[nSpace], lambda, Hnp1;
476  double sign;
477  int same_sign=1;
478  if (preRedistancingStage==1)
479  {
481  // *************** EIKONAL EQUATION *************** //
483  double norm_grad_u=0, norm_grad_un=0;
484  for (int I=0;I<nSpace; I++)
485  {
486  norm_grad_u += grad_u[I]*grad_u[I];
487  norm_grad_un += grad_un[I]*grad_un[I];
488  }
489  norm_grad_u = std::sqrt(norm_grad_u)+1E-10;
490  norm_grad_un = std::sqrt(norm_grad_un)+1E-10;
491  // residual of Eikonal equation
492  double epsRedist = epsFactRedist*hK;
493  double Si = -1.0+2.0*smoothedHeaviside(epsRedist,un);
494  residualEikonal = Si*(norm_grad_u-1.0);
495  delta = smoothedDirac(epsRedist,un);
496 
497  // compute (lagged) velocity for redistancing //
498  double Un[nSpace];
499  double normUn=0;
500  for (int I=0; I < nSpace; I++)
501  {
502  Un[I] = Si*grad_un[I]/norm_grad_un;
503  normUn += Un[I]*Un[I];
504  }
505  normUn = sqrt(normUn)+1E-10;
506  // compute coefficient for stabilization
507  tau = 0.5*hK;
508  }
509  else //clsvof model
510  {
512  // *************** CLSVOF MODEL *************** //
515  // MOVING MESH //
517  double mesh_velocity[3];
518  mesh_velocity[0] = xt;
519  mesh_velocity[1] = yt;
520  mesh_velocity[2] = zt;
521 
523  // GENERAL STUFF //
525  double epsHeaviside = epsFactHeaviside*hK;
526  double Sn = smoothedSign(epsHeaviside,un);
527  double Snp1 = smoothedSign(epsHeaviside,u);
528  Hnp1 = smoothedHeaviside(epsHeaviside,u);
529 
531  // LAMBDA //
533  lambda = lambdaFact*hK/norm_factor_lagged;
534  if (LAMBDA_SCALING==1)
535  {
536  double deltaHat = fmax(smoothedNormalizedDirac(2*epsHeaviside,un),1E-6);
537  lambda = lambdaFact*deltaHat;
538  }
539 
541  // TIME DERIVATIVE //
543  time_derivative_residual = porosity*(Snp1-Sn)/dt;
544 
546  // ADVECTIVE TERM //
548  double relative_velocity[nSpace], relative_velocity_old[nSpace];
549  for (int I=0;I<nSpace;I++)
550  {
551  // compute relative velocity //
552  relative_velocity[I] = (velocity.data()[eN_k_nSpace+I]
553  -MOVING_DOMAIN*mesh_velocity[I]);
554  relative_velocity_old[I] = (velocity_old.data()[eN_k_nSpace+I]
555  -MOVING_DOMAIN*mesh_velocity[I]);
556  // compute advection term //
557  //fnp1[I] = relative_velocity[I]*Snp1; //implicit advection via BDF1
558  //fnHalf[I] = 0.5*(relative_velocity[I]*Snp1
559  // +relative_velocity_old[I]*Sn); //implicit advection via CN
560  fnHalf[I] = 0.5*relative_velocity[I]*porosity*(Snp1+Sn);
561  }
562 
564  // CALCULATE CELL BASED CFL //
566  calculateCFL(elementDiameter.data()[eN]/degree_polynomial,relative_velocity,cfl.data()[eN_k]);
567  //calculateNonlinearCFL(elementDiameter.data()[eN]/degree_polynomial,
568  // relative_velocity,
569  // norm_factor_lagged,
570  // epsFactHeaviside,
571  // lambdaFact,
572  // cfl.data()[eN_k]);
573 
575  // CALCULATE min, max and mean distance //
577  if (eN<nElements_owned) // locally owned?
578  {
579  min_distance.data()[0] = fmin(min_distance.data()[0],fabs(u));
580  max_distance.data()[0] = fmax(max_distance.data()[0],fabs(u));
581  mean_distance.data()[0] += fabs(u)*dV;
582  //min_distance.data()[0] = fmin(min_distance.data()[0],u);
583  //max_distance.data()[0] = fmax(max_distance.data()[0],u);
584  //mean_distance.data()[0] += u*dV;
585  volume_domain.data()[0] += dV;
586  }
588  // SAVE SOLUTION // for other models
590  q_u.data()[eN_k] = u;
591  q_m.data()[eN_k] = porosity*Hnp1;
592  q_H.data()[eN_k] = Hnp1;
593  q_mH.data()[eN_k] = porosity*Hnp1; //porosity*H(\phi)=(1-q_vos.data())*H(\phi)
594  // gradient //
595  for (int I=0;I<nSpace;I++)
596  q_n.data()[eN_k_nSpace+I] = grad_u[I];
597  }
598 
600  // *************** LOOP ON DOFs *************** //
602  for(int i=0;i<nDOF_test_element;i++)
603  {
604  int eN_i=eN*nDOF_test_element+i;
605  int gi = offset_u+stride_u*u_l2g.data()[eN_i]; //global i-th index
606  register int i_nSpace=i*nSpace;
607  if (preRedistancingStage==1)
608  {
610  // EIKONAL RESIDUAL //
612  elementResidual_u[i] +=
613  alpha*(u_dof.data()[gi]-u_dof_old.data()[gi])*delta*u_test_dV[i] // BCs
614  + residualEikonal*u_test_dV[i] // Eikonal eqn
615  // Dmitri's ~SUPG + background dissipation //
616  + ck.NumericalDiffusion(tau+backgroundDissipation,
617  grad_u,
618  &u_grad_test_dV[i_nSpace])
619  - ck.NumericalDiffusion(tau,
620  normalReconstruction,
621  &u_grad_test_dV[i_nSpace]);
622  }
623  else // clsvof
624  {
626  // LOCATE THE INTERFACE //
628  if (i==0)
629  sign = Sign(u_dof.data()[gi]);
630  else if (same_sign==1)
631  {
632  same_sign = sign == Sign(u_dof.data()[gi]) ? 1 : 0;
633  sign = Sign(u_dof.data()[gi]);
634  }
636  // for visualization of VOF field //
638  element_rhs_L2proj_H[i] += Hnp1*u_test_dV[i];
639  //if (timeOrder==1)
640  //{
642  // CLSVOF RESIDUAL //
644  elementResidual_u[i] +=
645  // TIME DERIVATIVE
646  time_derivative_residual*u_test_dV[i]
647  // ADVECTION TERM. This is IMPLICIT
648  + ck.Advection_weak(fnHalf,&u_grad_test_dV[i_nSpace])
649  // REGULARIZATION TERM. This is IMPLICIT
650  + lambda*(ck.NumericalDiffusion(1.0,
651  grad_u,
652  &u_grad_test_dV[i_nSpace])
653  // TARGET for PENALIZATION. This is EXPLICIT
654  - ck.NumericalDiffusion(1.0,
655  normalReconstruction,
656  &u_grad_test_dV[i_nSpace]));
657  //}
658  //else // timeOrder=2
659  //{
660  //elementResidual_u[i] +=
661  // TIME DERIVATIVE
662  //time_derivative_residual*u_test_dV[i]
663  // ADVECTION TERM. This is IMPLICIT
664  //+ ck.Advection_weak(fnHalf,&u_grad_test_dV[i_nSpace])
665  // REGULARIZATION TERM. This is IMPLICIT
666  //+ lambda*ck.NumericalDiffusion(1.0,
667  // grad_unHalf,
668  // &u_grad_test_dV[i_nSpace])
669  // TARGET for PENALIZATION. This is EXPLICIT
670  //- lambda*ck.NumericalDiffusion(1.0,
671  // normalReconstruction,
672  // &u_grad_test_dV[i_nSpace]);
673  //}
674  }
675  }//i
676  if (preRedistancingStage==0)
677  {
678  if (same_sign == 0) // This cell contains the interface
679  {
680  for(int i=0;i<nDOF_test_element;i++)
681  {
682  int gi = offset_u+stride_u*u_l2g.data()[eN*nDOF_test_element+i];
683  interface_locator.data()[gi] = 1.0;
684  }
685  }
686  }
687  } //k
689  // DISTRIBUTE // load cell based element into global residual
691  for(int i=0;i<nDOF_test_element;i++)
692  {
693  int eN_i=eN*nDOF_test_element+i;
694  int gi = offset_u+stride_u*u_l2g.data()[eN_i]; //global i-th index
695  // distribute global residual for (lumped) mass matrix
696  globalResidual.data()[gi] += elementResidual_u[i];
697  if (preRedistancingStage==0)
698  H_dof.data()[gi] += element_rhs_L2proj_H[i]; // int(H*wi*dx)
699  }//i
700  }//elements
701  if (preRedistancingStage==0)
702  {
703  // COMPUTE LUMPED L2 PROJECTION
704  for (int i=0; i<numDOFs; i++)
705  H_dof.data()[i] /= lumped_mass_matrix.data()[i];
706  }
707 
709  // BOUNDARY //
711  //ebNE is the Exterior element boundary INdex
712  //ebN is the element boundary INdex
713  //eN is the element index
714  if (preRedistancingStage==0)
715  for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
716  {
717  register int ebN = exteriorElementBoundariesArray.data()[ebNE],
718  eN = elementBoundaryElementsArray.data()[ebN*2+0],
719  ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
720  eN_nDOF_trial_element = eN*nDOF_trial_element;
721  register double elementResidual_u[nDOF_test_element];
722  for (int i=0;i<nDOF_test_element;i++)
723  {
724  elementResidual_u[i]=0.0;
725  }
726  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
727  {
728  register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
729  ebNE_kb_nSpace = ebNE_kb*nSpace,
730  ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
731  ebN_local_kb_nSpace = ebN_local_kb*nSpace;
732  register double
733  u_ext=0.0,un_ext=0.0,grad_u_ext[nSpace],
734  df_ext[nSpace],
735  flux_ext=0.0,
736  bc_u_ext=0.0,
737  jac_ext[nSpace*nSpace],
738  jacDet_ext,
739  jacInv_ext[nSpace*nSpace],
740  boundaryJac[nSpace*(nSpace-1)],
741  metricTensor[(nSpace-1)*(nSpace-1)],
742  metricTensorDetSqrt,
743  dS,
744  u_test_dS[nDOF_test_element],
745  u_grad_trial_trace[nDOF_trial_element*nSpace],
746  normal[nSpace],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
747  //VRANS
748  porosity_ext;
749  //
750  //calculate the solution and gradients at quadrature points
751  //
752  //compute information about mapping from reference element to physical element
753  ck.calculateMapping_elementBoundary(eN,
754  ebN_local,
755  kb,
756  ebN_local_kb,
757  mesh_dof.data(),
758  mesh_l2g.data(),
759  mesh_trial_trace_ref.data(),
760  mesh_grad_trial_trace_ref.data(),
761  boundaryJac_ref.data(),
762  jac_ext,
763  jacDet_ext,
764  jacInv_ext,
765  boundaryJac,
766  metricTensor,
767  metricTensorDetSqrt,
768  normal_ref.data(),
769  normal,
770  x_ext,y_ext,z_ext);
771  ck.calculateMappingVelocity_elementBoundary(eN,
772  ebN_local,
773  kb,
774  ebN_local_kb,
775  mesh_velocity_dof.data(),
776  mesh_l2g.data(),
777  mesh_trial_trace_ref.data(),
778  xt_ext,yt_ext,zt_ext,
779  normal,
780  boundaryJac,
781  metricTensor,
782  integralScaling);
783  dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt
784  + MOVING_DOMAIN*integralScaling)*dS_ref.data()[kb];
785  ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],
786  jacInv_ext,
787  u_grad_trial_trace);
788  //solution at quad points
789  ck.valFromDOF(u_dof.data(),
790  &u_l2g.data()[eN_nDOF_trial_element],
791  &u_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],
792  u_ext);
793  ck.valFromDOF(u_dof_old.data(),
794  &u_l2g.data()[eN_nDOF_trial_element],
795  &u_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],
796  un_ext);
797  ck.gradFromDOF(u_dof.data(),
798  &u_l2g.data()[eN_nDOF_trial_element],
799  u_grad_trial_trace,
800  grad_u_ext);
801  //precalculate test function products with integration weights
802  for (int j=0;j<nDOF_trial_element;j++)
803  u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
804  //
805  //load the boundary values
806  //
807  double hK = elementDiameter.data()[eN]/degree_polynomial;
808  double epsHeaviside = epsFactHeaviside*hK;
809  double Su_ext = smoothedSign(epsHeaviside,u_ext); //Sign(u_ext)
810  double Sun_ext = smoothedSign(epsHeaviside,un_ext); //Sign(un_ext)
811  // NOTE: ebqe_bc_u_ext is provided by the user as BCs for VOF (i.e., 0 or 1)
812  // SuDBC = Sign(uBC)
813  double SuBC = (USE_SIGN_FUNCTION == 0 ? ebqe_bc_u_ext.data()[ebNE_kb] :
814  2*ebqe_bc_u_ext.data()[ebNE_kb] - 1);
815  if (IMPLICIT_BCs==1)
816  bc_u_ext = (isDOFBoundary_u.data()[ebNE_kb]*SuBC
817  +(1-isDOFBoundary_u.data()[ebNE_kb])*Su_ext);
818  else
819  bc_u_ext = (isDOFBoundary_u.data()[ebNE_kb]*SuBC
820  +(1-isDOFBoundary_u.data()[ebNE_kb])*Sun_ext);
821  //VRANS
822  porosity_ext = 1.-ebqe_vos_ext.data()[ebNE_kb];
823 
824  //
825  //moving mesh
826  //
827  double mesh_velocity[3];
828  mesh_velocity[0] = xt_ext;
829  mesh_velocity[1] = yt_ext;
830  mesh_velocity[2] = zt_ext;
831 
832  for (int I=0;I<nSpace;I++)
833  df_ext[I] = porosity_ext*(ebqe_velocity_ext.data()[ebNE_kb_nSpace+I]
834  - MOVING_DOMAIN*mesh_velocity[I]);
835  //
836  //calculate the numerical fluxes
837  //
838  exteriorNumericalAdvectiveFlux(isDOFBoundary_u.data()[ebNE_kb],
839  isFluxBoundary_u.data()[ebNE_kb],
840  normal,
841  bc_u_ext, //{-1,1} or {0,1}
842  ebqe_bc_flux_u_ext.data()[ebNE_kb],
843  IMPLICIT_BCs == 1 ? Su_ext : Sun_ext,
844  df_ext, //VRANS includes porosity
845  flux_ext);
846  ebqe_flux.data()[ebNE_kb] = flux_ext;
847 
849  // save solution // for other models? cek need to be consistent with numerical flux
851  ebqe_u.data()[ebNE_kb] = u_ext;
852  // TODO: do I need ebqe_m?
853  ebqe_H.data()[ebNE_kb] = smoothedHeaviside(epsHeaviside,u_ext);
854  // gradient //
855  for (int I=0;I<nSpace;I++)
856  ebqe_n.data()[ebNE_kb_nSpace+I] = grad_u_ext[I];
857 
858  //
859  //update residuals
860  //
861  for (int i=0;i<nDOF_test_element;i++)
862  {
863  elementResidual_u[i] += ck.ExteriorElementBoundaryFlux(flux_ext,u_test_dS[i]);
864  }//i
865  }//kb
866  //
867  //update the element and global residual storage
868  //
869  for (int i=0;i<nDOF_test_element;i++)
870  {
871  int eN_i = eN*nDOF_test_element+i;
872  globalResidual.data()[offset_u+stride_u*u_l2g.data()[eN_i]] += elementResidual_u[i];
873  }//i
874  }//ebNE
875  // END OF BOUNDARY //
876  }
877 
878  void calculateJacobian(arguments_dict& args)
879  {
880  double dt = args.scalar<double>("dt");
881  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
882  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
883  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
884  xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
885  double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
886  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
887  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
888  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
889  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
890  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
891  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
892  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
893  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
894  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
895  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
896  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
897  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
898  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
899  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
900  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
901  int nElements_global = args.scalar<int>("nElements_global");
902  double useMetrics = args.scalar<double>("useMetrics");
903  const xt::pyarray<double>& q_vos = args.array<double>("q_vos");
904  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
905  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
906  xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
907  int degree_polynomial = args.scalar<int>("degree_polynomial");
908  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
909  xt::pyarray<double>& u_dof_old = args.array<double>("u_dof_old");
910  xt::pyarray<double>& velocity = args.array<double>("velocity");
911  xt::pyarray<double>& cfl = args.array<double>("cfl");
912  xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
913  xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
914  xt::pyarray<double>& globalJacobian = args.array<double>("globalJacobian");
915  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
916  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
917  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
918  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
919  xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
920  const xt::pyarray<double>& ebqe_vos_ext = args.array<double>("ebqe_vos_ext");
921  xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
922  xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
923  xt::pyarray<int>& isFluxBoundary_u = args.array<int>("isFluxBoundary_u");
924  xt::pyarray<double>& ebqe_bc_flux_u_ext = args.array<double>("ebqe_bc_flux_u_ext");
925  xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
926  int timeOrder = args.scalar<int>("timeOrder");
927  int timeStage = args.scalar<int>("timeStage");
928  double epsFactHeaviside = args.scalar<double>("epsFactHeaviside");
929  double epsFactDirac = args.scalar<double>("epsFactDirac");
930  double epsFactRedist = args.scalar<double>("epsFactRedist");
931  double lambdaFact = args.scalar<double>("lambdaFact");
932  int preRedistancingStage = args.scalar<int>("preRedistancingStage");
933  double norm_factor_lagged = args.scalar<double>("norm_factor_lagged");
934  double alpha = args.scalar<double>("alpha");
935  for(int eN=0;eN<nElements_global;eN++)
936  {
937  register double elementJacobian_u_u[nDOF_test_element][nDOF_trial_element];
938  for (int i=0;i<nDOF_test_element;i++)
939  for (int j=0;j<nDOF_trial_element;j++)
940  elementJacobian_u_u[i][j]=0.0;
941  for (int k=0;k<nQuadraturePoints_element;k++)
942  {
943  int eN_k = eN*nQuadraturePoints_element+k, //index to a scalar at a quadrature point
944  eN_k_nSpace = eN_k*nSpace,
945  eN_nDOF_trial_element = eN*nDOF_trial_element; //index to a vector at a quadrature point
946  //declare local storage
947  register double
948  u, un, u_grad_trial[nDOF_trial_element*nSpace],
949  grad_u[nSpace], grad_un[nSpace],
950  jac[nSpace*nSpace], jacDet, jacInv[nSpace*nSpace],
951  u_test_dV[nDOF_test_element], u_grad_test_dV[nDOF_test_element*nSpace],
952  dV, x,y,z,xt,yt,zt,h_phi,
953  porosity;
954  //get jacobian, etc for mapping reference element
955  ck.calculateMapping_element(eN,
956  k,
957  mesh_dof.data(),
958  mesh_l2g.data(),
959  mesh_trial_ref.data(),
960  mesh_grad_trial_ref.data(),
961  jac,
962  jacDet,
963  jacInv,
964  x,y,z);
965  ck.calculateH_element(eN,
966  k,
967  nodeDiametersArray.data(),
968  mesh_l2g.data(),
969  mesh_trial_ref.data(),
970  h_phi);
971  ck.calculateMappingVelocity_element(eN,
972  k,
973  mesh_velocity_dof.data(),
974  mesh_l2g.data(),
975  mesh_trial_ref.data(),
976  xt,yt,zt);
977  //get the physical integration weight
978  dV = fabs(jacDet)*dV_ref.data()[k];
979  //get the trial function gradients
980  ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],
981  jacInv,
982  u_grad_trial);
983  //get the solution
984  ck.valFromDOF(u_dof.data(),
985  &u_l2g.data()[eN_nDOF_trial_element],
986  &u_trial_ref.data()[k*nDOF_trial_element],
987  u);
988  ck.valFromDOF(u_dof_old.data(),
989  &u_l2g.data()[eN_nDOF_trial_element],
990  &u_trial_ref.data()[k*nDOF_trial_element],
991  un);
992  //get the solution gradients
993  ck.gradFromDOF(u_dof.data(),
994  &u_l2g.data()[eN_nDOF_trial_element],
995  u_grad_trial,
996  grad_u);
997  ck.gradFromDOF(u_dof_old.data(),
998  &u_l2g.data()[eN_nDOF_trial_element],
999  u_grad_trial,
1000  grad_un);
1001  //precalculate test function products with integration weights
1002  for (int j=0;j<nDOF_trial_element;j++)
1003  {
1004  u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
1005  for (int I=0;I<nSpace;I++)
1006  u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;
1007  }
1008  double hK=(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN])/degree_polynomial;
1009  //VRANS
1010  porosity=1.0-q_vos.data()[eN_k];
1011 
1012  double delta, dH[nSpace], tau, backgroundDissipation=0.1*hK;
1013  double lambda, time_derivative_jacobian, df[nSpace];
1014  if (preRedistancingStage==1)
1015  {
1016  // ************************************************ //
1017  // *************** EIKONAL EQUATION *************** //
1018  // ************************************************ //
1019  double norm_grad_u=0, norm_grad_un=0;
1020  for (int I=0;I<nSpace; I++)
1021  {
1022  norm_grad_u += grad_u[I]*grad_u[I];
1023  norm_grad_un += grad_un[I]*grad_un[I];
1024  }
1025  norm_grad_u = std::sqrt(norm_grad_u)+1E-10;
1026  norm_grad_un = std::sqrt(norm_grad_un)+1E-10;
1027  // derivative of residual of Eikonal equation //
1028  double epsRedist = epsFactRedist*hK;
1029  double Si = -1.0+2.0*smoothedHeaviside(epsRedist,un);
1030  for (int I=0; I<nSpace;I++)
1031  dH[I] = Si*grad_u[I]/norm_grad_u;
1032  delta = smoothedDirac(epsRedist,un);
1033 
1034  // compute lagged velocity of redistancing
1035  double Un[nSpace];
1036  double normUn = 0.;
1037  for (int I=0; I < nSpace; I++)
1038  {
1039  Un[I] = Si*grad_un[I]/norm_grad_un;
1040  normUn += Un[I]*Un[I];
1041  }
1042  normUn = sqrt(normUn)+1E-10;
1043  // compute tau coefficient
1044  tau = 0.5*hK;
1045  }
1046  else
1047  {
1048  // ******************************************** //
1049  // *************** CLSVOF MODEL *************** //
1050  // ******************************************** //
1052  // MOVING MESH //
1054  double mesh_velocity[3];
1055  mesh_velocity[0] = xt;
1056  mesh_velocity[1] = yt;
1057  mesh_velocity[2] = zt;
1058 
1060  // GENERAL STUFF //
1062  double epsDirac = epsFactDirac*hK;
1063  double dSnp1 = smoothedDerivativeSign(epsDirac,u); //derivative of smoothed sign
1064 
1066  // LAMBDA //
1068  lambda = lambdaFact*hK/norm_factor_lagged;
1069  if (LAMBDA_SCALING==1)
1070  {
1071  double deltaHat = fmax(smoothedNormalizedDirac(2*epsDirac,un),1E-6);
1072  lambda = lambdaFact*deltaHat;
1073  }
1074 
1076  // TIME DERIVATIVE //
1078  time_derivative_jacobian = porosity*dSnp1/dt;
1079 
1081  // ADVECTIVE TERM //
1083  double relative_velocity[nSpace];
1084  for (int I=0;I<nSpace;I++)
1085  {
1086  relative_velocity[I] = (velocity.data()[eN_k_nSpace+I]
1087  -MOVING_DOMAIN*mesh_velocity[I]);
1088  df[I] = relative_velocity[I]*porosity*dSnp1;
1089  }
1090  }
1091 
1093  // LOOP ON DOFs //
1095  for(int i=0;i<nDOF_test_element;i++)
1096  {
1097  for(int j=0;j<nDOF_trial_element;j++)
1098  {
1099  int j_nSpace = j*nSpace;
1100  int i_nSpace = i*nSpace;
1101 
1102  if (preRedistancingStage==1)
1103  {
1105  // EIKONAL JACOBIAN //
1107  elementJacobian_u_u[i][j] +=
1108  (i == j ? alpha*delta*u_test_dV[i] : 0.) // BCs
1109  + ck.HamiltonianJacobian_weak(dH, // Eikonal equation
1110  &u_grad_trial[j_nSpace],
1111  u_test_dV[i])
1112  // Dmitri's ~SUPG + background dissipation //
1113  +ck.NumericalDiffusionJacobian(tau+backgroundDissipation,
1114  &u_grad_trial[j_nSpace],
1115  &u_grad_test_dV[i_nSpace]);
1116  }
1117  else // clsvof model
1118  {
1120  // CLSVOF JACOBIAN //
1122  elementJacobian_u_u[i][j] +=
1123  // TIME DERIVATIVE
1124  time_derivative_jacobian*(u_trial_ref.data()[k*nDOF_trial_element+j]
1125  *u_test_dV[i])
1126  // IMPLICIT TERMS: ADVECTION, DIFFUSION
1127  + 0.5*ck.AdvectionJacobian_weak(df,
1128  u_trial_ref.data()[k*nDOF_trial_element+j],
1129  &u_grad_test_dV[i_nSpace])
1130  + lambda*ck.NumericalDiffusionJacobian(1.0,
1131  &u_grad_trial[j_nSpace],
1132  &u_grad_test_dV[i_nSpace]);
1133  }
1134  }//j
1135  }//i
1136  }//k
1137  //
1138  //load into element Jacobian into global Jacobian
1139  //
1140  for (int i=0;i<nDOF_test_element;i++)
1141  {
1142  int eN_i = eN*nDOF_test_element+i;
1143  for (int j=0;j<nDOF_trial_element;j++)
1144  {
1145  int eN_i_j = eN_i*nDOF_trial_element+j;
1146  globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_u_u.data()[eN_i_j]] +=
1147  elementJacobian_u_u[i][j];
1148  }//j
1149  }//i
1150  }//elements
1151 
1153  // BOUNDARY LOOP //
1155  if (IMPLICIT_BCs==1 && preRedistancingStage==0)
1156  for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
1157  {
1158  register int ebN = exteriorElementBoundariesArray.data()[ebNE];
1159  register int eN = elementBoundaryElementsArray.data()[ebN*2+0],
1160  ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
1161  eN_nDOF_trial_element = eN*nDOF_trial_element;
1162  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
1163  {
1164  register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
1165  ebNE_kb_nSpace = ebNE_kb*nSpace,
1166  ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
1167  ebN_local_kb_nSpace = ebN_local_kb*nSpace;
1168  register double u_ext=0.0,
1169  grad_u_ext[nSpace],
1170  m_ext=0.0,
1171  dm_ext=0.0,
1172  f_ext[nSpace],
1173  df_ext[nSpace],
1174  dflux_u_u_ext=0.0,
1175  //bc_grad_u_ext[nSpace],
1176  bc_m_ext=0.0,
1177  bc_dm_ext=0.0,
1178  bc_f_ext[nSpace],
1179  bc_df_ext[nSpace],
1180  fluxJacobian_u_u[nDOF_trial_element],
1181  jac_ext[nSpace*nSpace],
1182  jacDet_ext,
1183  jacInv_ext[nSpace*nSpace],
1184  boundaryJac[nSpace*(nSpace-1)],
1185  metricTensor[(nSpace-1)*(nSpace-1)],
1186  metricTensorDetSqrt,
1187  dS,
1188  u_test_dS[nDOF_test_element],
1189  u_grad_trial_trace[nDOF_trial_element*nSpace],
1190  normal[nSpace],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
1191  //VRANS
1192  porosity_ext;
1193  ck.calculateMapping_elementBoundary(eN,
1194  ebN_local,
1195  kb,
1196  ebN_local_kb,
1197  mesh_dof.data(),
1198  mesh_l2g.data(),
1199  mesh_trial_trace_ref.data(),
1200  mesh_grad_trial_trace_ref.data(),
1201  boundaryJac_ref.data(),
1202  jac_ext,
1203  jacDet_ext,
1204  jacInv_ext,
1205  boundaryJac,
1206  metricTensor,
1207  metricTensorDetSqrt,
1208  normal_ref.data(),
1209  normal,
1210  x_ext,y_ext,z_ext);
1211  ck.calculateMappingVelocity_elementBoundary(eN,
1212  ebN_local,
1213  kb,
1214  ebN_local_kb,
1215  mesh_velocity_dof.data(),
1216  mesh_l2g.data(),
1217  mesh_trial_trace_ref.data(),
1218  xt_ext,yt_ext,zt_ext,
1219  normal,
1220  boundaryJac,
1221  metricTensor,
1222  integralScaling);
1223  dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt
1224  + MOVING_DOMAIN*integralScaling)*dS_ref.data()[kb];
1225  ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],
1226  jacInv_ext,
1227  u_grad_trial_trace);
1228  ck.valFromDOF(u_dof.data(),
1229  &u_l2g.data()[eN_nDOF_trial_element],
1230  &u_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],
1231  u_ext);
1232  ck.gradFromDOF(u_dof.data(),
1233  &u_l2g.data()[eN_nDOF_trial_element],
1234  u_grad_trial_trace,
1235  grad_u_ext);
1236  //precalculate test function products with integration weights
1237  for (int j=0;j<nDOF_trial_element;j++)
1238  u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
1239  //
1240  //load the boundary values
1241  //
1242  double hK = elementDiameter.data()[eN]/degree_polynomial;
1243  double epsHeaviside = epsFactHeaviside*hK;
1244  double dSu_ext = smoothedDerivativeSign(epsHeaviside,u_ext);
1245  //VRANS
1246  porosity_ext = 1.-ebqe_vos_ext.data()[ebNE_kb];
1247  //
1248  //moving domain
1249  //
1250  double mesh_velocity[3];
1251  mesh_velocity[0] = xt_ext;
1252  mesh_velocity[1] = yt_ext;
1253  mesh_velocity[2] = zt_ext;
1254  //std::cout<<"ext J mesh_velocity"<<std::endl;
1255  for (int I=0;I<nSpace;I++)
1256  df_ext[I] = porosity_ext*(ebqe_velocity_ext.data()[ebNE_kb_nSpace+I]
1257  - MOVING_DOMAIN*mesh_velocity[I]);
1258  //
1259  //calculate the numerical fluxes
1260  //
1261  exteriorNumericalAdvectiveFluxDerivative(isDOFBoundary_u.data()[ebNE_kb],
1262  isFluxBoundary_u.data()[ebNE_kb],
1263  normal,
1264  dSu_ext,
1265  df_ext,//VRANS holds porosity
1266  dflux_u_u_ext);
1267  //
1268  //calculate the flux jacobian
1269  //
1270  for (int j=0;j<nDOF_trial_element;j++)
1271  {
1272  //register int ebNE_kb_j = ebNE_kb*nDOF_trial_element+j;
1273  register int ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j;
1274  fluxJacobian_u_u[j]=
1275  ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_u_u_ext,
1276  u_trial_trace_ref.data()[ebN_local_kb_j]);
1277  }//j
1278  //
1279  //update the global Jacobian from the flux Jacobian
1280  //
1281  for (int i=0;i<nDOF_test_element;i++)
1282  {
1283  register int eN_i = eN*nDOF_test_element+i;
1284  //register int ebNE_kb_i = ebNE_kb*nDOF_test_element+i;
1285  for (int j=0;j<nDOF_trial_element;j++)
1286  {
1287  register int ebN_i_j = ebN*4*nDOF_test_X_trial_element + i*nDOF_trial_element + j;
1288  globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j]]+=
1289  fluxJacobian_u_u[j]*u_test_dS[i];
1290  }//j
1291  }//i
1292  }//kb
1293  }//ebNE
1294  }//computeJacobian for MCorr with CLSVOF
1295 
1296  std::tuple<double, double, double, double, double, double, double, double, double, double, double>
1297  calculateMetricsAtEOS(arguments_dict& args)
1298  {
1299  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1300  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1301  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1302  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1303  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1304  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1305  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1306  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1307  int nElements_global = args.scalar<int>("nElements_global");
1308  int nElements_owned = args.scalar<int>("nElements_owned");
1309  int useMetrics = args.scalar<int>("useMetrics");
1310  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1311  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1312  xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
1313  double degree_polynomial = args.scalar<double>("degree_polynomial");
1314  double epsFactHeaviside = args.scalar<double>("epsFactHeaviside");
1315  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1316  xt::pyarray<double>& u0_dof = args.array<double>("u0_dof");
1317  xt::pyarray<double>& u_exact = args.array<double>("u_exact");
1318  int offset_u = args.scalar<int>("offset_u");
1319  int stride_u = args.scalar<int>("stride_u");
1320  double global_I_err = 0.0;
1321  double global_sI_err = 0.0;
1322  double global_V = 0.0;
1323  double global_V0 = 0.0;
1324  double global_sV = 0.0;
1325  double global_sV0 = 0.0;
1326  double global_D_err = 0.0;
1327  double global_L2_err = 0.0;
1328  double global_L2Banded_err = 0.0;
1329  double global_area_band = 0.0;
1330  double global_sH_L2_err = 0.0;
1332  // ** LOOP IN CELLS //
1334  for(int eN=0;eN<nElements_global;eN++)
1335  {
1336  if (eN<nElements_owned) // just consider the locally owned cells
1337  {
1338  //declare local storage for local contributions and initialize
1339  double
1340  cell_I_err = 0., cell_sI_err = 0.,
1341  cell_V = 0., cell_V0 = 0., cell_sV = 0., cell_sV0 = 0.,
1342  cell_D_err = 0.,
1343  cell_L2_err = 0.,
1344  cell_L2Banded_err = 0., cell_area_band = 0.,
1345  cell_sH_L2_err = 0.;
1346 
1347  //loop over quadrature points and compute integrands
1348  for (int k=0;k<nQuadraturePoints_element;k++)
1349  {
1350  //compute indeces and declare local storage
1351  register int eN_k = eN*nQuadraturePoints_element+k,
1352  eN_k_nSpace = eN_k*nSpace,
1353  eN_nDOF_trial_element = eN*nDOF_trial_element;
1354  register double
1355  u, u0, uh,
1356  u_grad_trial[nDOF_trial_element*nSpace],
1357  grad_uh[nSpace],
1358  //for general use
1359  jac[nSpace*nSpace], jacDet, jacInv[nSpace*nSpace],
1360  dV,x,y,z,h_phi;
1361  //get the physical integration weight
1362  ck.calculateMapping_element(eN,
1363  k,
1364  mesh_dof.data(),
1365  mesh_l2g.data(),
1366  mesh_trial_ref.data(),
1367  mesh_grad_trial_ref.data(),
1368  jac,
1369  jacDet,
1370  jacInv,
1371  x,y,z);
1372  ck.calculateH_element(eN,
1373  k,
1374  nodeDiametersArray.data(),
1375  mesh_l2g.data(),
1376  mesh_trial_ref.data(),
1377  h_phi);
1378  dV = fabs(jacDet)*dV_ref.data()[k];
1379  // get functions at quad points
1380  ck.valFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],uh);
1381  ck.valFromDOF(u0_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],u0);
1382  u = u_exact.data()[eN_k];
1383  // get gradients
1384  ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
1385  ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial,grad_uh);
1386 
1387  double epsHeaviside = epsFactHeaviside*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN])/degree_polynomial;
1388  // compute (smoothed) heaviside functions //
1389  double Hu0 = heaviside(u0);
1390  double Hu = heaviside(u);
1391  double Huh = heaviside(uh);
1392  double sHu0 = smoothedHeaviside(epsHeaviside,u0);
1393  double sHu = smoothedHeaviside(epsHeaviside,u);
1394  double sHuh = smoothedHeaviside(epsHeaviside,uh);
1396  // compute cell metrics //
1398  // metrics on the interface
1399  cell_I_err += fabs(Hu - Huh)*dV;
1400  cell_sI_err += fabs(sHu - sHuh)*dV;
1401  // L2 metrics on the level set
1402  cell_L2_err += std::pow(u-uh,2)*dV;
1403  if (fabs(uh) <= 2*epsHeaviside)
1404  {
1405  cell_L2Banded_err += std::pow(u-uh,2)*dV;
1406  cell_area_band += dV;
1407  }
1408  // L2 metrics on the Heviside of the level set
1409  cell_sH_L2_err += std::pow(sHu-sHuh,2)*dV;
1410  // volume conservation
1411  cell_V += Huh*dV;
1412  cell_V0 += Hu0*dV;
1413  cell_sV += sHuh*dV;
1414  cell_sV0 += sHu0*dV;
1415 
1416  double norm2_grad_uh = 0.;
1417  for (int I=0; I<nSpace; I++)
1418  norm2_grad_uh += grad_uh[I]*grad_uh[I];
1419  cell_D_err += std::pow(std::sqrt(norm2_grad_uh) - 1, 2.)*dV;
1420  }
1421  global_V += cell_V;
1422  global_V0 += cell_V0;
1423  global_sV += cell_sV;
1424  global_sV0 += cell_sV0;
1425  // metrics //
1426  global_I_err += cell_I_err;
1427  global_sI_err += cell_sI_err;
1428  global_D_err += cell_D_err;
1429  global_L2_err += cell_L2_err;
1430  global_L2Banded_err += cell_L2Banded_err;
1431  global_area_band += cell_area_band;
1432  global_sH_L2_err += cell_sH_L2_err;
1433  }//elements
1434  }
1435  global_D_err *= 0.5;
1436 
1437  return std::tuple<double, double, double, double, double, double, double, double, double, double, double>(global_I_err, global_sI_err, global_V, global_V0, global_sV, global_sV0, global_D_err, global_L2_err, global_L2Banded_err, global_area_band, global_sH_L2_err);
1438  }
1439 
1440  std::tuple<double, double, double, double, double> calculateMetricsAtETS(arguments_dict& args)
1441  {
1442  double dt = args.scalar<double>("dt");
1443  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1444  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1445  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1446  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1447  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1448  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1449  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1450  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1451  int nElements_global = args.scalar<int>("nElements_global");
1452  int nElements_owned = args.scalar<int>("nElements_owned");
1453  int useMetrics = args.scalar<int>("useMetrics");
1454  xt::pyarray<double>& q_vos = args.array<double>("q_vos");
1455  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1456  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1457  xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
1458  double degree_polynomial = args.scalar<double>("degree_polynomial");
1459  double epsFactHeaviside = args.scalar<double>("epsFactHeaviside");
1460  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1461  xt::pyarray<double>& u_dof_old = args.array<double>("u_dof_old");
1462  xt::pyarray<double>& u0_dof = args.array<double>("u0_dof");
1463  xt::pyarray<double>& velocity = args.array<double>("velocity");
1464  int offset_u = args.scalar<int>("offset_u");
1465  int stride_u = args.scalar<int>("stride_u");
1466  int numDOFs = args.scalar<int>("numDOFs");
1467  xt::pyarray<double>& R_vector = args.array<double>("R_vector");
1468  xt::pyarray<double>& sR_vector = args.array<double>("sR_vector");
1469  double global_V = 0.0;
1470  double global_V0 = 0.0;
1471  double global_sV = 0.0;
1472  double global_sV0 = 0.0;
1473  double global_D_err = 0.0;
1475  // ** LOOP IN CELLS FOR CELL BASED TERMS ** //
1477  for(int eN=0;eN<nElements_global;eN++)
1478  {
1479  //declare local storage for local contributions and initialize
1480  register double element_R[nDOF_test_element], element_sR[nDOF_test_element];
1481  for (int i=0;i<nDOF_test_element;i++)
1482  {
1483  element_R[i] = 0.;
1484  element_sR[i] = 0.;
1485  }
1486  double
1487  cell_V = 0., cell_V0 = 0., cell_sV = 0., cell_sV0 = 0.,
1488  cell_D_err = 0.;
1489  //loop over quadrature points and compute integrands
1490  for (int k=0;k<nQuadraturePoints_element;k++)
1491  {
1492  //compute indeces and declare local storage
1493  register int eN_k = eN*nQuadraturePoints_element+k,
1494  eN_k_nSpace = eN_k*nSpace,
1495  eN_nDOF_trial_element = eN*nDOF_trial_element;
1496  register double
1497  unp1, un, u0,
1498  grad_unp1[nSpace], sFlux_np1[nSpace], Flux_np1[nSpace],
1499  u_grad_trial[nDOF_trial_element*nSpace],
1500  u_grad_test_dV[nDOF_test_element*nSpace],
1501  u_test_dV[nDOF_trial_element],
1502  //for general use
1503  jac[nSpace*nSpace], jacDet, jacInv[nSpace*nSpace],
1504  dV,x,y,z,h_phi,
1505  porosity;
1506  //get the physical integration weight
1507  ck.calculateMapping_element(eN,
1508  k,
1509  mesh_dof.data(),
1510  mesh_l2g.data(),
1511  mesh_trial_ref.data(),
1512  mesh_grad_trial_ref.data(),
1513  jac,
1514  jacDet,
1515  jacInv,
1516  x,y,z);
1517  ck.calculateH_element(eN,
1518  k,
1519  nodeDiametersArray.data(),
1520  mesh_l2g.data(),
1521  mesh_trial_ref.data(),
1522  h_phi);
1523  dV = fabs(jacDet)*dV_ref.data()[k];
1524  // get functions at quad points
1525  ck.valFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],unp1);
1526  ck.valFromDOF(u_dof_old.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],un);
1527  ck.valFromDOF(u0_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],u0);
1528  // get gradients
1529  ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
1530  ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial,grad_unp1);
1531  //precalculate test function products with integration weights for mass matrix terms
1532  for (int j=0;j<nDOF_trial_element;j++)
1533  {
1534  u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
1535  for (int I=0;I<nSpace;I++)
1536  u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;
1537  }
1538 
1539  porosity = 1.0-q_vos.data()[eN_k];
1540 
1541  double epsHeaviside = epsFactHeaviside*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN])/degree_polynomial;
1542  // compute (smoothed) heaviside functions //
1543  double Hu0 = heaviside(u0);
1544  double Hunp1 = heaviside(unp1);
1545  double sHu0 = smoothedHeaviside(epsHeaviside,u0);
1546  double sHunp1 = smoothedHeaviside(epsHeaviside,unp1);
1547 
1548  // compute cell metrics //
1549  cell_V += porosity*Hunp1*dV;
1550  cell_V0 += porosity*Hu0*dV;
1551  cell_sV += porosity*sHunp1*dV;
1552  cell_sV0 += porosity*sHu0*dV;
1553 
1554  double norm2_grad_unp1 = 0.;
1555  for (int I=0; I<nSpace; I++)
1556  norm2_grad_unp1 += grad_unp1[I]*grad_unp1[I];
1557  cell_D_err += std::pow(std::sqrt(norm2_grad_unp1) - 1, 2.)*dV;
1558 
1559  double Sunp1 = Sign(unp1);
1560  double Sun = Sign(un);
1561  double sSunp1 = smoothedSign(epsHeaviside,unp1);
1562  double sSun = smoothedSign(epsHeaviside,un);
1563  for (int I=0; I<nSpace; I++)
1564  {
1565  Flux_np1[I] = velocity.data()[eN_k_nSpace+I]*Sunp1;
1566  sFlux_np1[I] = velocity.data()[eN_k_nSpace+I]*sSunp1;
1567  }
1568 
1569  for(int i=0;i<nDOF_test_element;i++)
1570  {
1571  register int i_nSpace=i*nSpace;
1572  element_R[i] += ((Sunp1-Sun)/dt*u_test_dV[i]
1573  + ck.Advection_weak(Flux_np1,&u_grad_test_dV[i_nSpace]));
1574  element_sR[i] += ((sSunp1-sSun)/dt*u_test_dV[i]
1575  + ck.Advection_weak(sFlux_np1,&u_grad_test_dV[i_nSpace]));
1576  }
1577  }
1578  // DISTRIBUTE //
1579  for(int i=0;i<nDOF_test_element;i++)
1580  {
1581  int eN_i=eN*nDOF_test_element+i;
1582  int gi = offset_u+stride_u*u_l2g.data()[eN_i]; //global i-th index
1583  R_vector.data()[gi] += element_R[i];
1584  sR_vector.data()[gi] += element_sR[i];
1585  }
1586  if (eN<nElements_owned) // just consider the locally owned cells
1587  {
1588  global_V += cell_V;
1589  global_V0 += cell_V0;
1590  global_sV += cell_sV;
1591  global_sV0 += cell_sV0;
1592  global_D_err += cell_D_err;
1593  }
1594  }//elements
1595  global_D_err *= 0.5;
1596 
1597  return std::tuple<double, double, double, double, double>(global_V, global_V0, global_sV, global_sV0, global_D_err);
1598  }
1599 
1600  void normalReconstruction(arguments_dict& args)
1601  {
1602  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1603  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1604  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1605  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1606  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1607  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1608  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1609  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1610  int nElements_global = args.scalar<int>("nElements_global");
1611  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1612  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1613  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1614  int offset_u = args.scalar<int>("offset_u");
1615  int stride_u = args.scalar<int>("stride_u");
1616  int numDOFs = args.scalar<int>("numDOFs");
1617  xt::pyarray<double>& weighted_lumped_mass_matrix = args.array<double>("weighted_lumped_mass_matrix");
1618  xt::pyarray<double>& rhs_qx = args.array<double>("rhs_qx");
1619  xt::pyarray<double>& rhs_qy = args.array<double>("rhs_qy");
1620  xt::pyarray<double>& rhs_qz = args.array<double>("rhs_qz");
1621  xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
1622  xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
1623  xt::pyarray<double>& weighted_mass_matrix = args.array<double>("weighted_mass_matrix");
1624  for (int i=0; i<numDOFs; i++)
1625  {
1626  weighted_lumped_mass_matrix.data()[i]=0.;
1627  rhs_qx.data()[i]=0.;
1628  rhs_qy.data()[i]=0.;
1629  rhs_qz.data()[i]=0.;
1630  }
1631  for(int eN=0;eN<nElements_global;eN++)
1632  {
1633  //declare local storage for local contributions and initialize
1634  register double
1635  element_weighted_lumped_mass_matrix[nDOF_test_element],
1636  element_rhsx_normal_reconstruction[nDOF_test_element],
1637  element_rhsy_normal_reconstruction[nDOF_test_element],
1638  element_rhsz_normal_reconstruction[nDOF_test_element];
1639  register double element_weighted_mass_matrix[nDOF_test_element][nDOF_trial_element];
1640  for (int i=0;i<nDOF_test_element;i++)
1641  {
1642  element_weighted_lumped_mass_matrix[i]=0.0;
1643  element_rhsx_normal_reconstruction[i]=0.0;
1644  element_rhsy_normal_reconstruction[i]=0.0;
1645  element_rhsz_normal_reconstruction[i]=0.0;
1646  for (int j=0;j<nDOF_trial_element;j++)
1647  element_weighted_mass_matrix[i][j]=0.0;
1648  }
1649  //loop over quadrature points and compute integrands
1650  for(int k=0;k<nQuadraturePoints_element;k++)
1651  {
1652  //compute indeces and declare local storage
1653  register int eN_k = eN*nQuadraturePoints_element+k,
1654  eN_k_nSpace = eN_k*nSpace,
1655  eN_nDOF_trial_element = eN*nDOF_trial_element;
1656  register double
1657  //for mass matrix contributions
1658  grad_u[nSpace],
1659  u_grad_trial[nDOF_trial_element*nSpace],
1660  u_test_dV[nDOF_trial_element],
1661  //for general use
1662  jac[nSpace*nSpace], jacDet, jacInv[nSpace*nSpace],
1663  dV,x,y,z;
1664  //get the physical integration weight
1665  ck.calculateMapping_element(eN,
1666  k,
1667  mesh_dof.data(),
1668  mesh_l2g.data(),
1669  mesh_trial_ref.data(),
1670  mesh_grad_trial_ref.data(),
1671  jac,
1672  jacDet,
1673  jacInv,
1674  x,y,z);
1675  dV = fabs(jacDet)*dV_ref.data()[k];
1676  ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
1677  ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial,grad_u);
1678  //precalculate test function products with integration weights for mass matrix terms
1679  for (int j=0;j<nDOF_trial_element;j++)
1680  u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
1681 
1682  double rhsx = grad_u[0];
1683  double rhsy = grad_u[1];
1684  double rhsz = 0;
1685  if (nSpace==3)
1686  rhsz = grad_u[2];
1687 
1688  double norm_grad_u = 0;
1689  for (int I=0;I<nSpace; I++)
1690  norm_grad_u += grad_u[I]*grad_u[I];
1691  norm_grad_u = std::sqrt(norm_grad_u)+1E-10;
1692 
1693  for(int i=0;i<nDOF_test_element;i++)
1694  {
1695  element_weighted_lumped_mass_matrix[i] += norm_grad_u*u_test_dV[i];
1696  element_rhsx_normal_reconstruction[i] += rhsx*u_test_dV[i];
1697  element_rhsy_normal_reconstruction[i] += rhsy*u_test_dV[i];
1698  element_rhsz_normal_reconstruction[i] += rhsz*u_test_dV[i];
1699  for(int j=0;j<nDOF_trial_element;j++)
1700  element_weighted_mass_matrix[i][j] +=
1701  norm_grad_u*u_trial_ref.data()[k*nDOF_trial_element+j]*u_test_dV[i];
1702  }
1703  } //k
1704  // DISTRIBUTE //
1705  for(int i=0;i<nDOF_test_element;i++)
1706  {
1707  int eN_i=eN*nDOF_test_element+i;
1708  int gi = offset_u+stride_u*u_l2g.data()[eN_i]; //global i-th index
1709 
1710  weighted_lumped_mass_matrix.data()[gi] += element_weighted_lumped_mass_matrix[i];
1711  // rhs for reconstruction via consistent mass matrix
1712  rhs_qx.data()[gi] += element_rhsx_normal_reconstruction[i];
1713  rhs_qy.data()[gi] += element_rhsy_normal_reconstruction[i];
1714  rhs_qz.data()[gi] += element_rhsz_normal_reconstruction[i];
1715  for (int j=0;j<nDOF_trial_element;j++)
1716  {
1717  int eN_i_j = eN_i*nDOF_trial_element+j;
1718  weighted_mass_matrix.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_u_u.data()[eN_i_j]]
1719  += element_weighted_mass_matrix[i][j];
1720  }
1721  }//i
1722  }//elements
1723  }
1724 
1725  void calculateRhsL2Proj(arguments_dict& args)
1726  {
1727  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1728  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1729  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1730  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1731  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1732  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1733  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1734  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1735  int nElements_global = args.scalar<int>("nElements_global");
1736  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1737  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1738  double he_for_disc_ICs = args.scalar<double>("he_for_disc_ICs");
1739  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1740  int offset_u = args.scalar<int>("offset_u");
1741  int stride_u = args.scalar<int>("stride_u");
1742  int numDOFs = args.scalar<int>("numDOFs");
1743  xt::pyarray<double>& rhs_l2_proj = args.array<double>("rhs_l2_proj");
1744  for (int i=0; i<numDOFs; i++)
1745  rhs_l2_proj.data()[i]=0.;
1746  for(int eN=0;eN<nElements_global;eN++)
1747  {
1748  //declare local storage for local contributions and initialize
1749  register double element_rhs_l2_proj[nDOF_test_element];
1750  for (int i=0;i<nDOF_test_element;i++)
1751  element_rhs_l2_proj[i]=0.0;
1752 
1753  //loop over quadrature points and compute integrands
1754  for(int k=0;k<nQuadraturePoints_element;k++)
1755  {
1756  //compute indeces and declare local storage
1757  register int eN_k = eN*nQuadraturePoints_element+k,
1758  eN_k_nSpace = eN_k*nSpace,
1759  eN_nDOF_trial_element = eN*nDOF_trial_element;
1760  register double
1761  u,u_test_dV[nDOF_trial_element],
1762  //for general use
1763  jac[nSpace*nSpace], jacDet, jacInv[nSpace*nSpace],
1764  dV,x,y,z;
1765  //get the physical integration weight
1766  ck.calculateMapping_element(eN,
1767  k,
1768  mesh_dof.data(),
1769  mesh_l2g.data(),
1770  mesh_trial_ref.data(),
1771  mesh_grad_trial_ref.data(),
1772  jac,
1773  jacDet,
1774  jacInv,
1775  x,y,z);
1776  dV = fabs(jacDet)*dV_ref.data()[k];
1777  ck.valFromDOF(u_dof.data(),
1778  &u_l2g.data()[eN_nDOF_trial_element],
1779  &u_trial_ref.data()[k*nDOF_trial_element],
1780  u);
1781  //precalculate test function products with integration weights for mass matrix terms
1782  for (int j=0;j<nDOF_trial_element;j++)
1783  u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
1784 
1785  for(int i=0;i<nDOF_test_element;i++)
1786  element_rhs_l2_proj[i] += he_for_disc_ICs*u*u_test_dV[i];
1787  //element_rhs_l2_proj[i] += u*u_test_dV[i];
1788  } //k
1789  // DISTRIBUTE //
1790  for(int i=0;i<nDOF_test_element;i++)
1791  {
1792  int eN_i=eN*nDOF_test_element+i;
1793  int gi = offset_u+stride_u*u_l2g.data()[eN_i]; //global i-th index
1794  rhs_l2_proj.data()[gi] += element_rhs_l2_proj[i];
1795  }//i
1796  }//elements
1797  }
1798 
1799  void calculateLumpedMassMatrix(arguments_dict& args)
1800  {
1801  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1802  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1803  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1804  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1805  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1806  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1807  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1808  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1809  int nElements_global = args.scalar<int>("nElements_global");
1810  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1811  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1812  xt::pyarray<double>& lumped_mass_matrix = args.array<double>("lumped_mass_matrix");
1813  int offset_u = args.scalar<int>("offset_u");
1814  int stride_u = args.scalar<int>("stride_u");
1815  for(int eN=0;eN<nElements_global;eN++)
1816  {
1817  //declare local storage for local contributions and initialize
1818  register double element_lumped_mass_matrix[nDOF_test_element];
1819  for (int i=0;i<nDOF_test_element;i++)
1820  element_lumped_mass_matrix[i]=0.0;
1821  //loop over quadrature points and compute integrands
1822  for(int k=0;k<nQuadraturePoints_element;k++)
1823  {
1824  //compute indeces and declare local storage
1825  register int eN_k = eN*nQuadraturePoints_element+k,
1826  eN_k_nSpace = eN_k*nSpace,
1827  eN_nDOF_trial_element = eN*nDOF_trial_element;
1828  register double
1829  //for mass matrix contributions
1830  u_test_dV[nDOF_trial_element],
1831  //for general use
1832  jac[nSpace*nSpace], jacDet, jacInv[nSpace*nSpace],
1833  dV,x,y,z;
1834  //get the physical integration weight
1835  ck.calculateMapping_element(eN,
1836  k,
1837  mesh_dof.data(),
1838  mesh_l2g.data(),
1839  mesh_trial_ref.data(),
1840  mesh_grad_trial_ref.data(),
1841  jac,
1842  jacDet,
1843  jacInv,
1844  x,y,z);
1845  dV = fabs(jacDet)*dV_ref.data()[k];
1846  //precalculate test function products with integration weights for mass matrix terms
1847  for (int j=0;j<nDOF_trial_element;j++)
1848  u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
1849 
1850  for(int i=0;i<nDOF_test_element;i++)
1851  element_lumped_mass_matrix[i] += u_test_dV[i];
1852  } //k
1853  // DISTRIBUTE //
1854  for(int i=0;i<nDOF_test_element;i++)
1855  {
1856  int eN_i=eN*nDOF_test_element+i;
1857  int gi = offset_u+stride_u*u_l2g.data()[eN_i]; //global i-th index
1858  lumped_mass_matrix.data()[gi] += element_lumped_mass_matrix[i];
1859  }//i
1860  }//elements
1861  }
1862 
1863  void assembleSpinUpSystem(arguments_dict& args)
1864  {
1865  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1866  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1867  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1868  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1869  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1870  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1871  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1872  int nElements_global = args.scalar<int>("nElements_global");
1873  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1874  xt::pyarray<double>& uInitial = args.array<double>("uInitial");
1875  int offset_u = args.scalar<int>("offset_u");
1876  int stride_u = args.scalar<int>("stride_u");
1877  xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
1878  xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
1879  xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
1880  xt::pyarray<double>& globalMassMatrix = args.array<double>("globalMassMatrix");
1881  for(int eN=0;eN<nElements_global;eN++)
1882  {
1883  register double
1884  elementResidual_u[nDOF_test_element],
1885  elementMassMatrix_u_u[nDOF_test_element][nDOF_trial_element];
1886  for (int i=0;i<nDOF_test_element;i++)
1887  {
1888  elementResidual_u[i]=0;
1889  for (int j=0;j<nDOF_trial_element;j++)
1890  elementMassMatrix_u_u[i][j]=0.0;
1891  }
1892  for (int k=0;k<nQuadraturePoints_element;k++)
1893  {
1894  int eN_k = eN*nQuadraturePoints_element+k, //index to a scalar at a quadrature point
1895  eN_nDOF_trial_element = eN*nDOF_trial_element; //index to a vector at a quadrature point
1896  //declare local storage
1897  register double
1898  jac[nSpace*nSpace], jacDet, jacInv[nSpace*nSpace],
1899  u_test_dV[nDOF_test_element],
1900  dV, x,y,z;
1901  //get jacobian, etc for mapping reference element
1902  ck.calculateMapping_element(eN,
1903  k,
1904  mesh_dof.data(),
1905  mesh_l2g.data(),
1906  mesh_trial_ref.data(),
1907  mesh_grad_trial_ref.data(),
1908  jac,
1909  jacDet,
1910  jacInv,
1911  x,y,z);
1912  //get the physical integration weight
1913  dV = fabs(jacDet)*dV_ref.data()[k];
1914  //precalculate test function products with integration weights
1915  for (int j=0;j<nDOF_trial_element;j++)
1916  u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
1917 
1919  // LOOP ON DOFs //
1921  for(int i=0;i<nDOF_test_element;i++)
1922  {
1923  elementResidual_u[i] += uInitial.data()[eN_k]*u_test_dV[i];
1924  for(int j=0;j<nDOF_trial_element;j++)
1925  {
1926  elementMassMatrix_u_u[i][j] +=
1927  u_trial_ref.data()[k*nDOF_trial_element+j]*u_test_dV[i];
1928  }//j
1929  }//i
1930  }//k
1931  //
1932  //load into element Jacobian into global Jacobian
1933  //
1934  for (int i=0;i<nDOF_test_element;i++)
1935  {
1936  int eN_i = eN*nDOF_test_element+i;
1937  int gi = offset_u+stride_u*u_l2g.data()[eN_i]; //global i-th index
1938  globalResidual.data()[gi] += elementResidual_u[i];
1939  for (int j=0;j<nDOF_trial_element;j++)
1940  {
1941  int eN_i_j = eN_i*nDOF_trial_element+j;
1942  globalMassMatrix.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_u_u.data()[eN_i_j]] +=
1943  elementMassMatrix_u_u[i][j];
1944  }//j
1945  }//i
1946  }//elements
1947  }
1948 
1949  void FCTStep(arguments_dict& args)
1950  {
1951  int NNZ = args.scalar<int>("NNZ");
1952  int numDOFs = args.scalar<int>("numDOFs");
1953  xt::pyarray<double>& lumped_mass_matrix = args.array<double>("lumped_mass_matrix");
1954  xt::pyarray<double>& soln = args.array<double>("soln");
1955  xt::pyarray<double>& solH = args.array<double>("solH");
1956  xt::pyarray<double>& solL = args.array<double>("solL");
1957  xt::pyarray<double>& limited_solution = args.array<double>("limited_solution");
1958  xt::pyarray<int>& csrRowIndeces_DofLoops = args.array<int>("csrRowIndeces_DofLoops");
1959  xt::pyarray<int>& csrColumnOffsets_DofLoops = args.array<int>("csrColumnOffsets_DofLoops");
1960  xt::pyarray<double>& MassMatrix = args.array<double>("matrix");
1961  Rpos.resize(numDOFs, 0.0);
1962  Rneg.resize(numDOFs, 0.0);
1963  FluxCorrectionMatrix.resize(NNZ, 0.0);
1965  // LOOP in DOFs //
1967  int ij=0;
1968  for (int i=0; i<numDOFs; i++)
1969  {
1970  //read some vectors
1971  double solHi = solH.data()[i];
1972  double solLi = solL.data()[i];
1973  double mi = lumped_mass_matrix.data()[i];
1974 
1975  double mini=-1.0, maxi=1.0; // global FCT
1976  //double mini=1.0E10, maxi=-1.0E10;
1977  double Pposi=0, Pnegi=0;
1978  // LOOP OVER THE SPARSITY PATTERN (j-LOOP)//
1979  for (int offset=csrRowIndeces_DofLoops.data()[i]; offset<csrRowIndeces_DofLoops.data()[i+1]; offset++)
1980  {
1981  int j = csrColumnOffsets_DofLoops.data()[offset];
1982  // i-th row of flux correction matrix
1983  FluxCorrectionMatrix[ij] = ((i==j ? 1. : 0.)*mi - MassMatrix.data()[ij]) * (solH.data()[j]-solHi);
1984 
1985  //mini = fmin(mini,limited_solution.data()[j]);
1986  //maxi = fmax(maxi,limited_solution.data()[j]);
1987 
1989  // COMPUTE P VECTORS //
1991  Pposi += FluxCorrectionMatrix[ij]*((FluxCorrectionMatrix[ij] > 0) ? 1. : 0.);
1992  Pnegi += FluxCorrectionMatrix[ij]*((FluxCorrectionMatrix[ij] < 0) ? 1. : 0.);
1993 
1994  //update ij
1995  ij+=1;
1996  }
1998  // COMPUTE Q VECTORS //
2000  double Qposi = mi*(maxi-solLi);
2001  double Qnegi = mi*(mini-solLi);
2002 
2004  // COMPUTE R VECTORS //
2006  Rpos[i] = ((Pposi==0) ? 1. : std::min(1.0,Qposi/Pposi));
2007  Rneg[i] = ((Pnegi==0) ? 1. : std::min(1.0,Qnegi/Pnegi));
2008  } // i DOFs
2009 
2011  // COMPUTE LIMITERS //
2013  ij=0;
2014  for (int i=0; i<numDOFs; i++)
2015  {
2016  double ith_Limiter_times_FluxCorrectionMatrix = 0.;
2017  double Rposi = Rpos[i], Rnegi = Rneg[i];
2018  // LOOP OVER THE SPARSITY PATTERN (j-LOOP)//
2019  for (int offset=csrRowIndeces_DofLoops.data()[i]; offset<csrRowIndeces_DofLoops.data()[i+1]; offset++)
2020  {
2021  int j = csrColumnOffsets_DofLoops.data()[offset];
2022  ith_Limiter_times_FluxCorrectionMatrix +=
2023  ((FluxCorrectionMatrix[ij]>0) ? std::min(Rposi,Rneg[j]) : std::min(Rnegi,Rpos[j]))
2024  * FluxCorrectionMatrix[ij];
2025  //ith_Limiter_times_FluxCorrectionMatrix += FluxCorrectionMatrix[ij];
2026  //update ij
2027  ij+=1;
2028  }
2029  limited_solution.data()[i] = solL.data()[i] + 1./lumped_mass_matrix.data()[i]*ith_Limiter_times_FluxCorrectionMatrix;
2030  }
2031  }
2032  };//CLSVOF
2033 
2034  inline CLSVOF_base* newCLSVOF(int nSpaceIn,
2035  int nQuadraturePoints_elementIn,
2036  int nDOF_mesh_trial_elementIn,
2037  int nDOF_trial_elementIn,
2038  int nDOF_test_elementIn,
2039  int nQuadraturePoints_elementBoundaryIn,
2040  int CompKernelFlag)
2041  {
2042  if (nSpaceIn == 2)
2043  return proteus::chooseAndAllocateDiscretization2D<CLSVOF_base,CLSVOF,CompKernel>(nSpaceIn,
2044  nQuadraturePoints_elementIn,
2045  nDOF_mesh_trial_elementIn,
2046  nDOF_trial_elementIn,
2047  nDOF_test_elementIn,
2048  nQuadraturePoints_elementBoundaryIn,
2049  CompKernelFlag);
2050  else
2051  return proteus::chooseAndAllocateDiscretization<CLSVOF_base,CLSVOF,CompKernel>(nSpaceIn,
2052  nQuadraturePoints_elementIn,
2053  nDOF_mesh_trial_elementIn,
2054  nDOF_trial_elementIn,
2055  nDOF_test_elementIn,
2056  nQuadraturePoints_elementBoundaryIn,
2057  CompKernelFlag);
2058  }
2059 }//proteus
2060 #endif
proteus::CLSVOF::smoothedNormalizedDirac
double smoothedNormalizedDirac(double eps, double u)
Definition: CLSVOF.h:202
proteus::CLSVOF_base::FluxCorrectionMatrix
std::valarray< double > FluxCorrectionMatrix
Definition: CLSVOF.h:34
proteus::CLSVOF_base::Rpos
std::valarray< double > Rpos
Definition: CLSVOF.h:34
sign
#define sign(x, y)
Definition: jf.h:44
proteus::CLSVOF::assembleSpinUpSystem
void assembleSpinUpSystem(arguments_dict &args)
Definition: CLSVOF.h:1863
proteus::CLSVOF_base::normalReconstruction
virtual void normalReconstruction(arguments_dict &args)=0
proteus::CLSVOF_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)=0
proteus::CLSVOF::CLSVOF
CLSVOF()
Definition: CLSVOF.h:60
n
Int n
Definition: Headers.h:28
proteus::CLSVOF_base::calculateMetricsAtEOS
virtual std::tuple< double, double, double, double, double, double, double, double, double, double, double > calculateMetricsAtEOS(xt::pyarray< double > &mesh_trial_ref, xt::pyarray< double > &mesh_grad_trial_ref, xt::pyarray< double > &mesh_dof, xt::pyarray< int > &mesh_l2g, xt::pyarray< double > &dV_ref, xt::pyarray< double > &u_trial_ref, xt::pyarray< double > &u_grad_trial_ref, xt::pyarray< double > &u_test_ref, int nElements_global, int nElements_owned, int useMetrics, xt::pyarray< int > &u_l2g, xt::pyarray< double > &elementDiameter, xt::pyarray< double > &nodeDiametersArray, double degree_polynomial, double epsFactHeaviside, xt::pyarray< double > &u_dof, xt::pyarray< double > &u0_dof, xt::pyarray< double > &u_exact, int offset_u, int stride_u)=0
df
double df(double C, double b, double a, int q, int r)
Definition: analyticalSolutions.c:2209
proteus::CLSVOF::calculateMetricsAtEOS
std::tuple< double, double, double, double, double, double, double, double, double, double, double > calculateMetricsAtEOS(arguments_dict &args)
Definition: CLSVOF.h:1297
proteus::CLSVOF_base::FCTStep
virtual void FCTStep(arguments_dict &args)=0
ModelFactory.h
CompKernel.h
proteus::CLSVOF_base::~CLSVOF_base
virtual ~CLSVOF_base()
Definition: CLSVOF.h:35
proteus::arguments_dict::scalar
T & scalar(const std::string &key)
proteus::CLSVOF::calculateJacobian
void calculateJacobian(arguments_dict &args)
Definition: CLSVOF.h:878
proteus::arguments_dict::array
xt::pyarray< T > & array(const std::string &key)
proteus::CLSVOF::normalReconstruction
void normalReconstruction(arguments_dict &args)
Definition: CLSVOF.h:1600
proteus::newCLSVOF
CLSVOF_base * newCLSVOF(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition: CLSVOF.h:2034
H
Double H
Definition: Headers.h:65
proteus::CLSVOF::calculateMetricsAtETS
std::tuple< double, double, double, double, double > calculateMetricsAtETS(arguments_dict &args)
Definition: CLSVOF.h:1440
min
#define min(a, b)
Definition: jf.h:71
proteus::CLSVOF::FCTStep
void FCTStep(arguments_dict &args)
Definition: CLSVOF.h:1949
proteus::CLSVOF::calculateLumpedMassMatrix
void calculateLumpedMassMatrix(arguments_dict &args)
Definition: CLSVOF.h:1799
proteus::CLSVOF_base::assembleSpinUpSystem
virtual void assembleSpinUpSystem(arguments_dict &args)=0
proteus::CLSVOF::smoothedHeaviside
double smoothedHeaviside(double eps, double u)
Definition: CLSVOF.h:176
proteus::CLSVOF::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: CLSVOF.h:95
proteus::CLSVOF::smoothedSign
double smoothedSign(double eps, double u)
Definition: CLSVOF.h:214
z
Double * z
Definition: Headers.h:49
proteus::CLSVOF_base::calculateMetricsAtEOS
virtual std::tuple< double, double, double, double, double, double, double, double, double, double, double > calculateMetricsAtEOS(arguments_dict &args)=0
proteus::heaviside
double heaviside(const double &z)
Definition: CLSVOF.h:20
proteus::CLSVOF::calculateResidual
void calculateResidual(arguments_dict &args)
Definition: CLSVOF.h:246
u
Double u
Definition: Headers.h:89
proteus::CLSVOF_base::calculateLumpedMassMatrix
virtual void calculateLumpedMassMatrix(arguments_dict &args)=0
proteus::CLSVOF::smoothedDirac
double smoothedDirac(double eps, double u)
Definition: CLSVOF.h:190
xt
Definition: AddedMass.cpp:7
proteus::CLSVOF::calculateCFL
void calculateCFL(const double &elementDiameter, const double df[nSpace], double &cfl)
Definition: CLSVOF.h:66
proteus::CLSVOF_base::calculateResidual
virtual void calculateResidual(arguments_dict &args)=0
USE_SIGN_FUNCTION
#define USE_SIGN_FUNCTION
Definition: CLSVOF.h:13
proteus::CLSVOF::ck
CompKernelType ck
Definition: CLSVOF.h:59
LAMBDA_SCALING
#define LAMBDA_SCALING
Definition: CLSVOF.h:15
proteus::CLSVOF_base::calculateMetricsAtETS
virtual std::tuple< double, double, double, double, double > calculateMetricsAtETS(arguments_dict &args)=0
proteus
Definition: ADR.h:17
proteus::Sign
double Sign(const double &z)
Definition: CLSVOF.h:23
proteus::CLSVOF_base
Definition: CLSVOF.h:31
proteus::CLSVOF::calculateNonlinearCFL
void calculateNonlinearCFL(const double &elementDiameter, const double df[nSpace], const double norm_factor_lagged, const double epsFactHeaviside, const double lambdaFact, double &cfl)
Definition: CLSVOF.h:79
proteus::CLSVOF_base::Rneg
std::valarray< double > Rneg
Definition: CLSVOF.h:34
proteus::arguments_dict
Definition: ArgumentsDict.h:70
proteus::CLSVOF::smoothedDerivativeSign
double smoothedDerivativeSign(double eps, double u)
Definition: CLSVOF.h:231
proteus::CLSVOF
Definition: CLSVOF.h:56
proteus::CLSVOF::calculateRhsL2Proj
void calculateRhsL2Proj(arguments_dict &args)
Definition: CLSVOF.h:1725
proteus::CLSVOF::nDOF_test_X_trial_element
const int nDOF_test_X_trial_element
Definition: CLSVOF.h:58
proteus::CLSVOF::exteriorNumericalAdvectiveFluxDerivative
void exteriorNumericalAdvectiveFluxDerivative(const int &isDOFBoundary_u, const int &isFluxBoundary_u, const double n[nSpace], const double &du, const double velocity[nSpace], double &dflux)
Definition: CLSVOF.h:137
ArgumentsDict.h
IMPLICIT_BCs
#define IMPLICIT_BCs
Definition: CLSVOF.h:14
proteus::CLSVOF_base::calculateRhsL2Proj
virtual void calculateRhsL2Proj(arguments_dict &args)=0
Generated on Fri Jul 1 2022 11:15:47 for proteus by doxygen 1.8.20