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