proteus 1.9.0
C/C++/Fortran libraries
Loading...
Searching...
No Matches
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
11namespace py = pybind11;
12
13#define USE_SIGN_FUNCTION 1
14#define IMPLICIT_BCs 0
15#define LAMBDA_SCALING 0
16
17namespace 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
28namespace 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 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 int eN_k = eN*nQuadraturePoints_element+k,
349 eN_k_nSpace = eN_k*nSpace,
350 eN_nDOF_trial_element = eN*nDOF_trial_element;
351 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 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 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 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 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 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 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 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 int ebN = exteriorElementBoundariesArray.data()[ebNE];
1159 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 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 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 //int ebNE_kb_j = ebNE_kb*nDOF_trial_element+j;
1273 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 int eN_i = eN*nDOF_test_element+i;
1284 //int ebNE_kb_i = ebNE_kb*nDOF_test_element+i;
1285 for (int j=0;j<nDOF_trial_element;j++)
1286 {
1287 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 int eN_k = eN*nQuadraturePoints_element+k,
1352 eN_k_nSpace = eN_k*nSpace,
1353 eN_nDOF_trial_element = eN*nDOF_trial_element;
1354 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 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 int eN_k = eN*nQuadraturePoints_element+k,
1494 eN_k_nSpace = eN_k*nSpace,
1495 eN_nDOF_trial_element = eN*nDOF_trial_element;
1496 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 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 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 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 int eN_k = eN*nQuadraturePoints_element+k,
1654 eN_k_nSpace = eN_k*nSpace,
1655 eN_nDOF_trial_element = eN*nDOF_trial_element;
1656 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 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 int eN_k = eN*nQuadraturePoints_element+k,
1758 eN_k_nSpace = eN_k*nSpace,
1759 eN_nDOF_trial_element = eN*nDOF_trial_element;
1760 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 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 int eN_k = eN*nQuadraturePoints_element+k,
1826 eN_k_nSpace = eN_k*nSpace,
1827 eN_nDOF_trial_element = eN*nDOF_trial_element;
1828 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 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 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
USE_SIGN_FUNCTION
#define USE_SIGN_FUNCTION
Definition CLSVOF.h:13
IMPLICIT_BCs
#define IMPLICIT_BCs
Definition CLSVOF.h:14
LAMBDA_SCALING
#define LAMBDA_SCALING
Definition CLSVOF.h:15
n
Int n
Definition Headers.h:28
H
Double H
Definition Headers.h:65
u
Double u
Definition Headers.h:89
z
Double * z
Definition Headers.h:49
proteus::CLSVOF_base
Definition CLSVOF.h:31
proteus::CLSVOF_base::FCTStep
virtual void FCTStep(arguments_dict &args)=0
proteus::CLSVOF_base::Rneg
std::valarray< double > Rneg
Definition CLSVOF.h:34
proteus::CLSVOF_base::calculateMetricsAtETS
virtual std::tuple< double, double, double, double, double > calculateMetricsAtETS(arguments_dict &args)=0
proteus::CLSVOF_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)=0
proteus::CLSVOF_base::FluxCorrectionMatrix
std::valarray< double > FluxCorrectionMatrix
Definition CLSVOF.h:34
proteus::CLSVOF_base::calculateLumpedMassMatrix
virtual void calculateLumpedMassMatrix(arguments_dict &args)=0
proteus::CLSVOF_base::calculateMetricsAtEOS
virtual std::tuple< double, double, double, double, double, double, double, double, double, double, double > calculateMetricsAtEOS(arguments_dict &args)=0
proteus::CLSVOF_base::~CLSVOF_base
virtual ~CLSVOF_base()
Definition CLSVOF.h:35
proteus::CLSVOF_base::normalReconstruction
virtual void normalReconstruction(arguments_dict &args)=0
proteus::CLSVOF_base::calculateResidual
virtual void calculateResidual(arguments_dict &args)=0
proteus::CLSVOF_base::Rpos
std::valarray< double > Rpos
Definition CLSVOF.h:34
proteus::CLSVOF_base::assembleSpinUpSystem
virtual void assembleSpinUpSystem(arguments_dict &args)=0
proteus::CLSVOF_base::calculateRhsL2Proj
virtual void calculateRhsL2Proj(arguments_dict &args)=0
proteus::CLSVOF::calculateRhsL2Proj
void calculateRhsL2Proj(arguments_dict &args)
Definition CLSVOF.h:1725
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
proteus::CLSVOF::calculateResidual
void calculateResidual(arguments_dict &args)
Definition CLSVOF.h:246
proteus::CLSVOF::calculateCFL
void calculateCFL(const double &elementDiameter, const double df[nSpace], double &cfl)
Definition CLSVOF.h:66
proteus::CLSVOF::smoothedHeaviside
double smoothedHeaviside(double eps, double u)
Definition CLSVOF.h:176
proteus::CLSVOF::smoothedNormalizedDirac
double smoothedNormalizedDirac(double eps, double u)
Definition CLSVOF.h:202
proteus::CLSVOF::smoothedSign
double smoothedSign(double eps, double u)
Definition CLSVOF.h:214
proteus::CLSVOF::calculateLumpedMassMatrix
void calculateLumpedMassMatrix(arguments_dict &args)
Definition CLSVOF.h:1799
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::nDOF_test_X_trial_element
const int nDOF_test_X_trial_element
Definition CLSVOF.h:58
proteus::CLSVOF::ck
CompKernelType ck
Definition CLSVOF.h:59
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::normalReconstruction
void normalReconstruction(arguments_dict &args)
Definition CLSVOF.h:1600
proteus::CLSVOF::calculateMetricsAtETS
std::tuple< double, double, double, double, double > calculateMetricsAtETS(arguments_dict &args)
Definition CLSVOF.h:1440
proteus::CLSVOF::assembleSpinUpSystem
void assembleSpinUpSystem(arguments_dict &args)
Definition CLSVOF.h:1863
proteus::CLSVOF::calculateJacobian
void calculateJacobian(arguments_dict &args)
Definition CLSVOF.h:878
proteus::CLSVOF::smoothedDerivativeSign
double smoothedDerivativeSign(double eps, double u)
Definition CLSVOF.h:231
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::smoothedDirac
double smoothedDirac(double eps, double u)
Definition CLSVOF.h:190
proteus::CLSVOF::FCTStep
void FCTStep(arguments_dict &args)
Definition CLSVOF.h:1949
df
double df(double C, double b, double a, int q, int r)
Definition analyticalSolutions.c:2209
sign
#define sign(x, y)
Definition jf.h:44
proteus
Definition ADR.h:19
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
proteus::chooseAndAllocateDiscretization
Model_Base * chooseAndAllocateDiscretization(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nDOF_v_trial_elementIn, int nDOF_v_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition MixedModelFactory.h:48
proteus::Sign
double Sign(const double &z)
Definition CLSVOF.h:23
proteus::heaviside
double heaviside(const double &z)
Definition CLSVOF.h:20
proteus::chooseAndAllocateDiscretization2D
Model_Base * chooseAndAllocateDiscretization2D(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nDOF_v_trial_elementIn, int nDOF_v_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition MixedModelFactory.h:309
proteus::arguments_dict
Definition ArgumentsDict.h:70
proteus::arguments_dict::scalar
T & scalar(const std::string &key)
proteus::arguments_dict::array
xt::pyarray< T > & array(const std::string &key)
Generated on for proteus by doxygen 1.18.0