proteus 1.9.0
C/C++/Fortran libraries
Loading...
Searching...
No Matches
VOF.h
Go to the documentation of this file.
1#ifndef VOF_H
2#define VOF_H
3#include <cmath>
4#include <iostream>
5#include <set>
6#include <map>
7#include <valarray>
8#include "CompKernel.h"
9#include "ModelFactory.h"
10#include "equivalent_polynomials.h"
11#include "ArgumentsDict.h"
12#include "xtensor-python/pyarray.hpp"
13
14namespace py = pybind11;
15
16#define POWER_SMOOTHNESS_INDICATOR 2
17#define IS_BETAij_ONE 0
18#define GLOBAL_FCT 0
19
20// Cell based methods:
21// * SUPG with BDF1 or BDF2 time integration
22// * Explicit Taylor Galerkin with EV stabilization
23// Edge based methods.
24// Low order via D. Kuzmin's
25// High order methods: Smoothness indicator with MC, EV commutator with MC, D.K with ML
26// Zalesak's FCT
27
28namespace proteus
29{
30 // Power entropy //
31 inline double ENTROPY(const double& phi, const double& phiL, const double& phiR){
32 return 1./2.*std::pow(fabs(phi),2.);
33 }
34 inline double DENTROPY(const double& phi, const double& phiL, const double& phiR){
35 return fabs(phi)*(phi>=0 ? 1 : -1);
36 }
37 // Log entropy // for level set from 0 to 1
38 inline double ENTROPY_LOG(const double& phi, const double& phiL, const double& phiR){
39 return std::log(fabs((phi-phiL)*(phiR-phi))+1E-14);
40 }
41 inline double DENTROPY_LOG(const double& phi, const double& phiL, const double& phiR){
42 return (phiL+phiR-2*phi)*((phi-phiL)*(phiR-phi)>=0 ? 1 : -1)/(fabs((phi-phiL)*(phiR-phi))+1E-14);
43 }
44}
45
46namespace proteus
47{
48 template<int nSpace, int nP, int nQ, int nEBQ>
49 // The trailing flag is the IFEM gate. It has always been false here -- by
50 // default rather than by statement -- and must stay false: with it set, an
51 // element whose interface passes through an edge or corner node takes the
52 // IFEM branch in Simplex::set_quad, which forces D to 0 and H/ImH to a hard
53 // 0/1 instead of the moment fit. That deletes the interface measure this
54 // model integrates over. Written out so the choice is visible at the call
55 // site and cannot change underneath us if the template default changes.
56 using GeneralizedFunctions = equivalent_polynomials::GeneralizedFunctions_mix<nSpace, nP, nP, nQ, nEBQ, false>;
57
58 class VOF_base
59 {
60 //The base class defining the interface
61 public:
62 std::valarray<double> Rpos, Rneg;
63 std::valarray<double> FluxCorrectionMatrix;
64 std::valarray<double> TransportMatrix, TransposeTransportMatrix;
65 std::valarray<double> psi, eta, global_entropy_residual, boundary_integral;
66 std::valarray<double> maxVel,maxEntRes;
67 virtual ~VOF_base(){}
68 virtual void calculateResidualElementBased(arguments_dict& args)=0;
69 virtual void calculateJacobian(arguments_dict& args)=0;
70 virtual void FCTStep(arguments_dict& args)=0;
71 virtual void calculateResidualEdgeBased(arguments_dict& args)=0;
72 };
73
74 template<class CompKernelType,
75 int nSpace,
76 int nQuadraturePoints_element,
77 int nDOF_mesh_trial_element,
78 int nDOF_trial_element,
79 int nDOF_test_element,
80 int nQuadraturePoints_elementBoundary>
81 class VOF : public VOF_base
82 {
83 public:
84 std::set<int> cutfem_boundaries;
85 std::map<int, int> cutfem_local_boundaries;
86 const int nDOF_test_X_trial_element;
87 CompKernelType ck;
88 GeneralizedFunctions<nSpace,2,nQuadraturePoints_element,nQuadraturePoints_elementBoundary> gf_s;
89 VOF():
90 nDOF_test_X_trial_element(nDOF_test_element*nDOF_trial_element),
91 ck()
92 {}
93
94 inline
95 void calculateCFL(const double& elementDiameter,
96 const double df[nSpace],
97 double& cfl)
98 {
99 double h,nrm_v;
100 h = elementDiameter;
101 nrm_v=0.0;
102 for(int I=0;I<nSpace;I++)
103 nrm_v+=df[I]*df[I];
104 nrm_v = sqrt(nrm_v);
105 cfl = nrm_v/h;
106 }
107
108 inline
109 void evaluateCoefficients(const double v[nSpace],
110 const double& u,
111 const double& porosity, //VRANS specific
112 double& m,
113 double& dm,
114 double f[nSpace],
115 double df[nSpace])
116 {
117 m = porosity*u;
118 dm= porosity;
119 for (int I=0; I < nSpace; I++)
120 {
121 f[I] = v[I]*porosity*u;
122 df[I] = v[I]*porosity;
123 }
124 }
125
126 inline
127 void calculateSubgridError_tau(const double& elementDiameter,
128 const double& dmt,
129 const double dH[nSpace],
130 double& cfl,
131 double& tau)
132 {
133 double h,nrm_v,oneByAbsdt;
134 h = elementDiameter;
135 nrm_v=0.0;
136 for(int I=0;I<nSpace;I++)
137 nrm_v+=dH[I]*dH[I];
138 nrm_v = sqrt(nrm_v);
139 cfl = nrm_v/h;
140 oneByAbsdt = fabs(dmt);
141 tau = 1.0/(2.0*nrm_v/h + oneByAbsdt + 1.0e-8);
142 }
143
144 inline
145 void calculateSubgridError_tau( const double& Ct_sge,
146 const double G[nSpace*nSpace],
147 const double& A0,
148 const double Ai[nSpace],
149 double& tau_v,
150 double& cfl)
151 {
152 double v_d_Gv=0.0;
153 for(int I=0;I<nSpace;I++)
154 for (int J=0;J<nSpace;J++)
155 v_d_Gv += Ai[I]*G[I*nSpace+J]*Ai[J];
156
157 tau_v = 1.0/sqrt(Ct_sge*A0*A0 + v_d_Gv + 1.0e-8);
158 }
159
160 inline
161 void calculateNumericalDiffusion(const double& shockCapturingDiffusion,
162 const double& elementDiameter,
163 const double& strong_residual,
164 const double grad_u[nSpace],
165 double& numDiff)
166 {
167 double h,
168 num,
169 den,
170 n_grad_u;
171 h = elementDiameter;
172 n_grad_u = 0.0;
173 for (int I=0;I<nSpace;I++)
174 n_grad_u += grad_u[I]*grad_u[I];
175 num = shockCapturingDiffusion*0.5*h*fabs(strong_residual);
176 den = sqrt(n_grad_u) + 1.0e-8;
177 numDiff = num/den;
178 }
179
180 inline
181 void exteriorNumericalAdvectiveFlux(const int& isDOFBoundary_u,
182 const int& isFluxBoundary_u,
183 const double n[nSpace],
184 const double& bc_u,
185 const double& bc_flux_u,
186 const double& u,
187 const double velocity[nSpace],
188 double& flux)
189 {
190
191 double flow=0.0;
192 for (int I=0; I < nSpace; I++)
193 flow += n[I]*velocity[I];
194 if (isDOFBoundary_u == 1)
195 {
196 if (flow >= 0.0)
197 {
198 flux = u*flow;
199 }
200 else
201 {
202 flux = bc_u*flow;
203 }
204 }
205 else if (isFluxBoundary_u == 1)
206 {
207 flux = bc_flux_u;
208 }
209 else
210 {
211 if (flow >= 0.0)
212 {
213 flux = u*flow;
214 }
215 else
216 {
217 std::cout<<"warning: VOF open boundary with no external trace, setting to zero for inflow"<<std::endl;
218 flux = 0.0;
219 }
220
221 }
222 }
223
224 inline
225 void exteriorNumericalAdvectiveFluxDerivative(const int& isDOFBoundary_u,
226 const int& isFluxBoundary_u,
227 const double n[nSpace],
228 const double velocity[nSpace],
229 double& dflux)
230 {
231 double flow=0.0;
232 for (int I=0; I < nSpace; I++)
233 flow += n[I]*velocity[I];
234 dflux=0.0;//default to no flux
235 if (isDOFBoundary_u == 1)
236 {
237 if (flow >= 0.0)
238 {
239 dflux = flow;
240 }
241 else
242 {
243 dflux = 0.0;
244 }
245 }
246 else if (isFluxBoundary_u == 1)
247 {
248 dflux = 0.0;
249 }
250 else
251 {
252 if (flow >= 0.0)
253 {
254 dflux = flow;
255 }
256 }
257 }
258
259 void calculateResidualElementBased(arguments_dict& args)
260 {
261 double dt = args.scalar<double>("dt");
262 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
263 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
264 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
265 xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
266 double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
267 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
268 xt::pyarray<double>& x_ref = args.array<double>("x_ref");
269 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
270 xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
271 xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
272 xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
273 xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
274 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
275 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
276 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
277 xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
278 xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
279 xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
280 xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
281 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
282 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
283 int nElements_global = args.scalar<int>("nElements_global");
284 double useMetrics = args.scalar<double>("useMetrics");
285 double alphaBDF = args.scalar<double>("alphaBDF");
286 int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
287 double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
288 double sc_uref = args.scalar<double>("sc_uref");
289 double sc_alpha = args.scalar<double>("sc_alpha");
290 const xt::pyarray<double>& q_porosity = args.array<double>("q_porosity");
291 const xt::pyarray<double>& porosity_dof = args.array<double>("porosity_dof");
292 xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
293 xt::pyarray<int>& r_l2g = args.array<int>("r_l2g");
294 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
295 xt::pyarray<double>& elementBoundaryDiameter = args.array<double>("elementBoundaryDiameter");
296 double degree_polynomial = args.scalar<double>("degree_polynomial");
297 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
298 xt::pyarray<double>& u_dof_old = args.array<double>("u_dof_old");
299 xt::pyarray<double>& velocity = args.array<double>("velocity");
300 xt::pyarray<double>& q_m = args.array<double>("q_m");
301 xt::pyarray<double>& q_u = args.array<double>("q_u");
302 xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
303 xt::pyarray<double>& q_dV = args.array<double>("q_dV");
304 xt::pyarray<double>& q_dV_last = args.array<double>("q_dV_last");
305 xt::pyarray<double>& cfl = args.array<double>("cfl");
306 xt::pyarray<double>& edge_based_cfl = args.array<double>("edge_based_cfl");
307 xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
308 xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
309 int offset_u = args.scalar<int>("offset_u");
310 int stride_u = args.scalar<int>("stride_u");
311 xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
312 xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
313 xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
314 xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
315 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
316 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
317 xt::pyarray<int>& elementBoundariesArray = args.array<int>("elementBoundariesArray");
318 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
319 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
320 xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
321 const xt::pyarray<double>& ebqe_porosity_ext = args.array<double>("ebqe_porosity_ext");
322 xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
323 xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
324 xt::pyarray<int>& isFluxBoundary_u = args.array<int>("isFluxBoundary_u");
325 xt::pyarray<double>& ebqe_bc_flux_u_ext = args.array<double>("ebqe_bc_flux_u_ext");
326 xt::pyarray<double>& ebqe_phi = args.array<double>("ebqe_phi");
327 double epsFact = args.scalar<double>("epsFact");
328 xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
329 xt::pyarray<double>& ebqe_flux = args.array<double>("ebqe_flux");
330 int stage = args.scalar<int>("stage");
331 xt::pyarray<double>& uTilde_dof = args.array<double>("uTilde_dof");
332 double cE = args.scalar<double>("cE");
333 double cMax = args.scalar<double>("cMax");
334 double cK = args.scalar<double>("cK");
335 double uL = args.scalar<double>("uL");
336 double uR = args.scalar<double>("uR");
337 int numDOFs = args.scalar<int>("numDOFs");
338 int NNZ = args.scalar<int>("NNZ");
339 xt::pyarray<int>& csrRowIndeces_DofLoops = args.array<int>("csrRowIndeces_DofLoops");
340 xt::pyarray<int>& csrColumnOffsets_DofLoops = args.array<int>("csrColumnOffsets_DofLoops");
341 xt::pyarray<int>& csrRowIndeces_CellLoops = args.array<int>("csrRowIndeces_CellLoops");
342 xt::pyarray<int>& csrColumnOffsets_CellLoops = args.array<int>("csrColumnOffsets_CellLoops");
343 xt::pyarray<int>& csrColumnOffsets_eb_CellLoops = args.array<int>("csrColumnOffsets_eb_CellLoops");
344 xt::pyarray<double>& ML = args.array<double>("ML");
345 int LUMPED_MASS_MATRIX = args.scalar<int>("LUMPED_MASS_MATRIX");
346 int STABILIZATION_TYPE = args.scalar<int>("STABILIZATION_TYPE");
347 int ENTROPY_TYPE = args.scalar<int>("ENTROPY_TYPE");
348 xt::pyarray<double>& uLow = args.array<double>("uLow");
349 xt::pyarray<double>& dLow = args.array<double>("dLow");
350 xt::pyarray<double>& dt_times_dH_minus_dL = args.array<double>("dt_times_dH_minus_dL");
351 xt::pyarray<double>& min_u_bc = args.array<double>("min_u_bc");
352 xt::pyarray<double>& max_u_bc = args.array<double>("max_u_bc");
353 xt::pyarray<double>& quantDOFs = args.array<double>("quantDOFs");
354 xt::pyarray<double>& ebqe_phi_s = args.array<double>("ebqe_phi_s");
355 double ghost_penalty_constant = args.scalar<double>("ghost_penalty_constant");
356 const xt::pyarray<double>& phi_solid = args.array<double>("phi_solid");
357 xt::pyarray<double>& phi_solid_nodes = args.array<double>("phi_solid_nodes");
358 bool useExact = args.scalar<int>("useExact");
359 xt::pyarray<double>& isActiveR = args.array<double>("isActiveR");
360 xt::pyarray<double>& isActiveDOF = args.array<double>("isActiveDOF");
361 xt::pyarray<int>& isActiveElement = args.array<int>("isActiveElement");
362 double meanEntropy = 0., meanOmega = 0., maxEntropy = -1E10, minEntropy = 1E10;
363 maxVel.resize(nElements_global, 0.0);
364 maxEntRes.resize(nElements_global, 0.0);
365 double Ct_sge = 4.0;
366 gf_s.useExact = useExact;
367 cutfem_boundaries.clear();
368 cutfem_local_boundaries.clear();
369 //
370 //loop over elements to compute volume integrals and load them into element and global residual
371 //
372 //eN is the element index
373 //eN_k is the quadrature point index for a scalar
374 //eN_k_nSpace is the quadrature point index for a vector
375 //eN_i is the element test function index
376 //eN_j is the element trial function index
377 //eN_k_j is the quadrature point index for a trial function
378 //eN_k_i is the quadrature point index for a trial function
379 for(int eN=0;eN<nElements_global;eN++)
380 {
381 //declare local storage for element residual and initialize
382 double elementResidual_u[nDOF_test_element];
383 bool element_active=false;
384 isActiveElement[eN]=0;
385 for (int i=0;i<nDOF_test_element;i++)
386 {
387 elementResidual_u[i]=0.0;
388 }//i
389 double element_phi_s[nDOF_mesh_trial_element];
390 for (int j=0;j<nDOF_mesh_trial_element;j++)
391 {
392 int eN_j = eN*nDOF_mesh_trial_element+j;
393 element_phi_s[j] = phi_solid_nodes.data()[u_l2g.data()[eN_j]];
394 }
395 double element_nodes[nDOF_mesh_trial_element*3];
396 for (int i=0;i<nDOF_mesh_trial_element;i++)
397 {
398 int eN_i=eN*nDOF_mesh_trial_element+i;
399 for(int I=0;I<3;I++)
400 element_nodes[i*3 + I] = mesh_dof.data()[mesh_l2g.data()[eN_i]*3 + I];
401 }//i
402 int icase_s = gf_s.calculate(element_phi_s, element_nodes, x_ref.data(),false);
403 if (icase_s == 0)
404 {
405 element_active=true;
406 isActiveElement[eN]=1;
407 //only works for simplices
408 for (int ebN_element=0;ebN_element < nDOF_mesh_trial_element; ebN_element++)
409 {
410 const int ebN = elementBoundariesArray.data()[eN*nDOF_mesh_trial_element+ebN_element];
411 //internal and actually a cut edge
412 //if (elementBoundaryElementsArray.data()[ebN*2+1] != -1 && (ebN < nElementBoundaries_owned) && element_phi_s[(ebN_element+1)%nDOF_mesh_trial_element]*element_phi_s[(ebN_element+2)%nDOF_mesh_trial_element] < 0.0)
413 if (elementBoundaryElementsArray[ebN*2+1] != -1 && element_phi_s[(ebN_element+1)%nDOF_mesh_trial_element]*element_phi_s[(ebN_element+2)%nDOF_mesh_trial_element] <= 0.0)
414 {
415 cutfem_boundaries.insert(ebN);
416 if (elementBoundaryElementsArray[ebN*2 + 0] == eN)
417 cutfem_local_boundaries[ebN] = ebN_element;
418 }
419 }
420 }
421 else if (icase_s == 1)
422 {
423 element_active=true;
424 isActiveElement[eN]=1;
425 }
426 //loop over quadrature points and compute integrands
427 for (int k=0;k<nQuadraturePoints_element;k++)
428 {
429 //compute indeces and declare local storage
430 int eN_k = eN*nQuadraturePoints_element+k,
431 eN_k_nSpace = eN_k*nSpace,
432 eN_nDOF_trial_element = eN*nDOF_trial_element;
433 double
434 entVisc_minus_artComp,
435 u=0.0,un=0.0,
436 grad_u[nSpace],grad_u_old[nSpace],grad_uTilde[nSpace],
437 m=0.0,dm=0.0,
438 H=0.0,Hn=0.0,HTilde=0.0,
439 f[nSpace],fn[nSpace],df[nSpace],
440 m_t=0.0,dm_t=0.0,
441 pdeResidual_u=0.0,
442 Lstar_u[nDOF_test_element],
443 subgridError_u=0.0,
444 tau=0.0,tau0=0.0,tau1=0.0,
445 numDiff0=0.0,numDiff1=0.0,
446 jac[nSpace*nSpace],
447 jacDet,
448 jacInv[nSpace*nSpace],
449 u_grad_trial[nDOF_trial_element*nSpace],
450 u_test_dV[nDOF_trial_element],
451 u_grad_test_dV[nDOF_test_element*nSpace],
452 dV,x,y,z,xt,yt,zt,
453 //VRANS
454 porosity,
455 //
456 G[nSpace*nSpace],G_dd_G,tr_G;//norm_Rv;
457 gf_s.set_quad(k);
458 const double H_s = gf_s.H(0.0,phi_solid.data()[eN_k]);
459
460 ck.calculateMapping_element(eN,
461 k,
462 mesh_dof.data(),
463 mesh_l2g.data(),
464 mesh_trial_ref.data(),
465 mesh_grad_trial_ref.data(),
466 jac,
467 jacDet,
468 jacInv,
469 x,y,z);
470 ck.calculateMappingVelocity_element(eN,
471 k,
472 mesh_velocity_dof.data(),
473 mesh_l2g.data(),
474 mesh_trial_ref.data(),
475 xt,yt,zt);
476 //get the physical integration weight
477 dV = fabs(jacDet)*dV_ref.data()[k];
478 ck.calculateG(jacInv,G,G_dd_G,tr_G);
479 //get the trial function gradients
480 ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],
481 jacInv,
482 u_grad_trial);
483 //get the solution
484 ck.valFromDOF(u_dof.data(),
485 &u_l2g.data()[eN_nDOF_trial_element],
486 &u_trial_ref.data()[k*nDOF_trial_element],
487 u);
488 ck.valFromDOF(u_dof_old.data(),
489 &u_l2g.data()[eN_nDOF_trial_element],
490 &u_trial_ref.data()[k*nDOF_trial_element],
491 un);
492 //get the solution gradients
493 ck.gradFromDOF(u_dof.data(),
494 &u_l2g.data()[eN_nDOF_trial_element],
495 u_grad_trial,
496 grad_u);
497 ck.gradFromDOF(u_dof_old.data(),
498 &u_l2g.data()[eN_nDOF_trial_element],
499 u_grad_trial,
500 grad_u_old);
501 ck.gradFromDOF(uTilde_dof.data(),
502 &u_l2g.data()[eN_nDOF_trial_element],
503 u_grad_trial,
504 grad_uTilde);
505 //precalculate test function products with integration weights
506 for (int j=0;j<nDOF_trial_element;j++)
507 {
508 u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
509 for (int I=0;I<nSpace;I++)
510 {
511 u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
512 }
513 }
514 //VRANS
515 porosity = q_porosity.data()[eN_k];
516 //
517 //
518 //calculate pde coefficients at quadrature points
519 //
520 evaluateCoefficients(&velocity.data()[eN_k_nSpace],
521 u,
522 //VRANS
523 porosity,
524 //
525 m,
526 dm,
527 f,
528 df);
529 //
530 //moving mesh
531 //
532 double mesh_velocity[3];
533 mesh_velocity[0] = xt;
534 mesh_velocity[1] = yt;
535 mesh_velocity[2] = zt;
536
537 for (int I=0;I<nSpace;I++)
538 {
539 f[I] -= MOVING_DOMAIN*m*mesh_velocity[I];
540 df[I] -= MOVING_DOMAIN*dm*mesh_velocity[I];
541 }
542 //
543 //calculate time derivative at quadrature points
544 //
545 if (q_dV_last.data()[eN_k] <= -100)
546 q_dV_last.data()[eN_k] = dV;
547 q_dV.data()[eN_k] = dV;
548 ck.bdf(alphaBDF,
549 q_m_betaBDF.data()[eN_k]*q_dV_last.data()[eN_k]/dV,//ensure prior mass integral is correct for m_t with BDF1
550 m,
551 dm,
552 m_t,
553 dm_t);
554
555 if (STABILIZATION_TYPE==1)
556 {
557 double normVel=0., norm_grad_un=0.;
558 for (int I=0;I<nSpace;I++)
559 {
560 Hn += df[I]*grad_u_old[I];
561 HTilde += df[I]*grad_uTilde[I];
562 fn[I] = porosity*df[I]*un-MOVING_DOMAIN*m*mesh_velocity[I];
563 H += df[I]*grad_u[I];
564 normVel += df[I]*df[I];
565 norm_grad_un += grad_u_old[I]*grad_u_old[I];
566 }
567 normVel = std::sqrt(normVel);
568 norm_grad_un = std::sqrt(norm_grad_un)+1E-10;
569
570 // calculate CFL
571 calculateCFL(elementDiameter.data()[eN]/degree_polynomial,df,cfl.data()[eN_k]);
572
573
574 // compute max velocity at cell
575 maxVel[eN] = fmax(normVel,maxVel[eN]);
576
577 // Strong entropy residual
578 double entRes = (ENTROPY(u,0,1)-ENTROPY(un,0,1))/dt + 0.5*(DENTROPY(u,0,1)*H +
579 DENTROPY(un,0,1)*Hn);
580 maxEntRes[eN] = fmax(maxEntRes[eN],fabs(entRes));
581
582 // Quantities for normalization factor //
583 meanEntropy += ENTROPY(u,0,1)*dV;
584 meanOmega += dV;
585 maxEntropy = fmax(maxEntropy,ENTROPY(u,0,1));
586 minEntropy = fmin(minEntropy,ENTROPY(u,0,1));
587
588 // artificial compression
589 double hK=elementDiameter.data()[eN]/degree_polynomial;
590 entVisc_minus_artComp = fmax(1-cK*fmax(un*(1-un),0)/hK/norm_grad_un,0);
591 }
592 else
593 {
594 //
595 //calculate subgrid error (strong residual and adjoint)
596 //
597 //calculate strong residual
598 pdeResidual_u = ck.Mass_strong(m_t) + ck.Advection_strong(df,grad_u);
599 //calculate adjoint
600 for (int i=0;i<nDOF_test_element;i++)
601 {
602 // int eN_k_i_nSpace = (eN_k*nDOF_trial_element+i)*nSpace;
603 // Lstar_u[i] = ck.Advection_adjoint(df,&u_grad_test_dV[eN_k_i_nSpace]);
604 int i_nSpace = i*nSpace;
605 Lstar_u[i] = ck.Advection_adjoint(df,&u_grad_test_dV[i_nSpace]);
606 }
607 //calculate tau and tau*Res
608 calculateSubgridError_tau(elementDiameter.data()[eN],dm_t,df,cfl.data()[eN_k],tau0);
609 calculateSubgridError_tau(Ct_sge,
610 G,
611 dm_t,
612 df,
613 tau1,
614 cfl.data()[eN_k]);
615 tau = useMetrics*tau1+(1.0-useMetrics)*tau0;
616
617 subgridError_u = -tau*pdeResidual_u;
618 //
619 //calculate shock capturing diffusion
620 //
621
622 ck.calculateNumericalDiffusion(shockCapturingDiffusion,
623 elementDiameter.data()[eN],
624 pdeResidual_u,
625 grad_u,
626 numDiff0);
627 //ck.calculateNumericalDiffusion(shockCapturingDiffusion,G,pdeResidual_u,grad_u_old,numDiff1);
628 ck.calculateNumericalDiffusion(shockCapturingDiffusion,
629 sc_uref,
630 sc_alpha,
631 G,
632 G_dd_G,
633 pdeResidual_u,
634 grad_u,
635 numDiff1);
636 q_numDiff_u.data()[eN_k] = useMetrics*numDiff1+(1.0-useMetrics)*numDiff0;
637 //std::cout<<tau<<" "<<q_numDiff_u.data()[eN_k]<<'\t'<<numDiff0<<'\t'<<numDiff1<<'\t'<<pdeResidual_u<<std::endl;
638
639 //
640 //update element residual
641 //
642
643
644 /* std::cout<<m_t<<'\t'
645 <<f[0]<<'\t'
646 <<f[1]<<'\t'
647 <<df[0]<<'\t'
648 <<df[1]<<'\t'
649 <<subgridError_u<<'\t'
650 <<q_numDiff_u_last.data()[eN_k]<<std::endl;*/
651 }
652
653 for(int i=0;i<nDOF_test_element;i++)
654 {
655 int eN_i=eN*nDOF_test_element+i;
656 //int eN_k_i=eN_k*nDOF_test_element+i,
657 //eN_k_i_nSpace = eN_k_i*nSpace,
658 int i_nSpace=i*nSpace;
659 if (STABILIZATION_TYPE==1)
660 {
661 if (stage == 1)
662 elementResidual_u[i] +=
663 ck.Mass_weak(dt*m_t,u_test_dV[i]) + // time derivative
664 1./3*dt*ck.Advection_weak(fn,&u_grad_test_dV[i_nSpace]) +
665 1./9*dt*dt*ck.NumericalDiffusion(Hn,df,&u_grad_test_dV[i_nSpace]) +
666 1./3*dt*entVisc_minus_artComp*ck.NumericalDiffusion(q_numDiff_u_last.data()[eN_k],
667 grad_u_old,
668 &u_grad_test_dV[i_nSpace]);
669 // TODO: Add part about moving mesh
670 else //stage 2
671 elementResidual_u[i] +=
672 ck.Mass_weak(dt*m_t,u_test_dV[i]) + // time derivative
673 dt*ck.Advection_weak(fn,&u_grad_test_dV[i_nSpace]) +
674 0.5*dt*dt*ck.NumericalDiffusion(HTilde,df,&u_grad_test_dV[i_nSpace]) +
675 dt*entVisc_minus_artComp*ck.NumericalDiffusion(q_numDiff_u_last.data()[eN_k],
676 grad_u_old,
677 &u_grad_test_dV[i_nSpace]);
678 }
679 else //supg
680 {
681 elementResidual_u[i] +=
682 H_s*(ck.Mass_weak(m_t,u_test_dV[i]) +
683 ck.Advection_weak(f,&u_grad_test_dV[i_nSpace]) +
684 ck.SubgridError(subgridError_u,Lstar_u[i]) +
685 ck.NumericalDiffusion(q_numDiff_u_last.data()[eN_k],
686 grad_u,
687 &u_grad_test_dV[i_nSpace]));
688 }
689 if (element_active)
690 {
691 isActiveR.data()[offset_u + stride_u*r_l2g.data()[eN_i]] = 1.0;
692 isActiveDOF.data()[u_l2g.data()[eN_i]] = 1.0;
693 }
694 }//i
695 //
696 //todo, get rid of m, since u=m
697 //save momentum for time history and velocity for subgrid error
698 //save solution for other models
699 //
700 q_u.data()[eN_k] = u;
701 q_m.data()[eN_k] = m;
702 }
703 //
704 //load element into global residual and save element residual
705 //
706 for(int i=0;i<nDOF_test_element;i++)
707 {
708 int eN_i=eN*nDOF_test_element+i;
709 globalResidual.data()[offset_u+stride_u*r_l2g.data()[eN_i]] += elementResidual_u[i];
710 }//i
711 }//elements
712 std::set<int>::iterator it=cutfem_boundaries.begin();
713 while(it!=cutfem_boundaries.end())
714 {
715 if(isActiveElement[elementBoundaryElementsArray[(*it)*2+0]] && isActiveElement[elementBoundaryElementsArray[(*it)*2+1]])
716 {
717 std::map<int,double> DW_Dn_jump;
718 double gamma_cutfem=ghost_penalty_constant,
719 h_cutfem=elementBoundaryDiameter.data()[*it];
720 int eN_nDOF_trial_element = elementBoundaryElementsArray.data()[(*it)*2+0]*nDOF_trial_element;
721 //See Massing Schott Wall 2018
722 //double norm_v=0.0;
723 //for (int i_offset=1;i_offset<nDOF_trial_element;i_offset++)//MSW18 is just on face, so trying to just use face dof
724 // {
725 // int i = (cutfem_local_boundaries[*it] + i_offset)%nDOF_trial_element;
726 // double u=u_old_dof.data()[vel_l2g.data()[eN_nDOF_trial_element+i]];
727 // v=v_old_dof.data()[vel_l2g.data()[eN_nDOF_v_trial_element+i]],
728 // w=w_old_dof.data()[vel_l2g.data()[eN_nDOF_v_trial_element+i]];
729 // norm_v=fmax(norm_v,sqrt(u*u+v*v+w*w));
730 // }
731 //double gamma_v_dim = rho_0*(nu_0 + norm_v*h_cutfem + alphaBDF*h_cutfem*h_cutfem);
732 //gamma_cutfem_p *= h_cutfem*h_cutfem/gamma_v_dim;
733 //if (NONCONSERVATIVE_FORM)
734 // gamma_cutfem*=gamma_v_dim;
735 //else
736 // gamma_cutfem*=(gamma_v_dim/rho_0);
737 for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
738 {
739 double Du_Dn_jump=0.0, dS;
740 for (int eN_side=0;eN_side < 2; eN_side++)
741 {
742 int ebN = *it,
743 eN = elementBoundaryElementsArray.data()[ebN*2+eN_side];
744 for (int i=0;i<nDOF_test_element;i++)
745 {
746 DW_Dn_jump[r_l2g.data()[eN*nDOF_test_element+i]] = 0.0;
747 }
748 }
749 for (int eN_side=0;eN_side < 2; eN_side++)
750 {
751 int ebN = *it,
752 eN = elementBoundaryElementsArray[ebN*2+eN_side],
753 ebN_local = elementBoundaryLocalElementBoundariesArray[ebN*2+eN_side],
754 eN_nDOF_trial_element = eN*nDOF_trial_element,
755 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
756 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
757 double u_int=0.0,
758 grad_u_int[nSpace],
759 jac_int[nSpace*nSpace],
760 jacDet_int,
761 jacInv_int[nSpace*nSpace],
762 boundaryJac[nSpace*(nSpace-1)],
763 metricTensor[(nSpace-1)*(nSpace-1)],
764 metricTensorDetSqrt,
765 u_test_dS[nDOF_test_element],
766 u_grad_trial_trace[nDOF_trial_element*nSpace],
767 u_grad_test_dS[nDOF_trial_element*nSpace],
768 normal[nSpace],x_int,y_int,z_int,xt_int,yt_int,zt_int,integralScaling,
769 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty;
770 for (int I=0; I<nSpace;I++)
771 grad_u_int[I] = 0.0;
772 //compute information about mapping from reference element to physical element
773 ck.calculateMapping_elementBoundary(eN,
774 ebN_local,
775 kb,
776 ebN_local_kb,
777 mesh_dof.data(),
778 mesh_l2g.data(),
779 mesh_trial_trace_ref.data(),
780 mesh_grad_trial_trace_ref.data(),
781 boundaryJac_ref.data(),
782 jac_int,
783 jacDet_int,
784 jacInv_int,
785 boundaryJac,
786 metricTensor,
787 metricTensorDetSqrt,
788 normal_ref.data(),
789 normal,
790 x_int,y_int,z_int);
791 //todo: check that physical coordinates match
792 ck.calculateMappingVelocity_elementBoundary(eN,
793 ebN_local,
794 kb,
795 ebN_local_kb,
796 mesh_velocity_dof.data(),
797 mesh_l2g.data(),
798 mesh_trial_trace_ref.data(),
799 xt_int,yt_int,zt_int,
800 normal,
801 boundaryJac,
802 metricTensor,
803 integralScaling);
804 dS = metricTensorDetSqrt*dS_ref.data()[kb];
805 //compute shape and solution information
806 //shape
807 ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_int,u_grad_trial_trace);
808 //solution and gradients
809 ck.valFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],u_int);
810 ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial_trace,grad_u_int);
811 for (int I=0;I<nSpace;I++)
812 {
813 Du_Dn_jump += grad_u_int[I]*normal[I];
814 }
815 for (int i=0;i<nDOF_test_element;i++)
816 {
817 int eN_i = eN*nDOF_test_element + i;
818 for (int I=0;I<nSpace;I++)
819 DW_Dn_jump[r_l2g[eN_i]] += u_grad_trial_trace[i*nSpace+I]*normal[I];
820 }
821 }//eN_side
822 for (std::map<int,double>::iterator W_it=DW_Dn_jump.begin(); W_it!=DW_Dn_jump.end(); ++W_it)
823 {
824 int i_global = W_it->first;
825 double DW_Dn_jump_i = W_it->second;
826 globalResidual.data()[offset_u+stride_u*i_global]+=gamma_cutfem*h_cutfem*Du_Dn_jump*DW_Dn_jump_i*dS;
827 }
828 }//kb
829 it++;
830 }
831 else
832 {
833 it = cutfem_boundaries.erase(it);
834 }
835 }//cutfem element boundaries
836 // std::set<int>::iterator it=cutfem_boundaries.begin();
837 // while(it!=cutfem_boundaries.end())
838 // {
839 // std::map<int,double> DW_Dn_jump;
840 // std::map<std::pair<int, int>, int> u_u_nz;
841 // double gamma_cutfem=ghost_penalty_constant,
842 // h_cutfem=elementBoundaryDiameter.data()[*it];
843 // int eN_nDOF_trial_element = elementBoundaryElementsArray.data()[(*it)*2+0]*nDOF_trial_element;
844 // //See Massing Schott Wall 2018
845 // //double norm_v=0.0;
846 // //for (int i_offset=1;i_offset<nDOF_v_trial_element;i_offset++)//MSW18 is just on face
847 // // {
848 // // int i = (cutfem_local_boundaries[*it] + i_offset)%nDOF_v_trial_element;//cek hack only works for P1
849 // // double u=u_old_dof.data()[vel_l2g.data()[eN_nDOF_v_trial_element+i]],
850 // // v=v_old_dof.data()[vel_l2g.data()[eN_nDOF_v_trial_element+i]],
851 // // w=w_old_dof.data()[vel_l2g.data()[eN_nDOF_v_trial_element+i]];
852 // // norm_v=fmax(norm_v,sqrt(u*u+v*v+w*w));
853 // // }
854 // ///double gamma_v_dim = rho_0*(nu_0 + norm_v*h_cutfem + alphaBDF*h_cutfem*h_cutfem);
855 // //gamma_cutfem_p *= h_cutfem*h_cutfem/gamma_v_dim;
856 // //if (NONCONSERVATIVE_FORM)
857 // // gamma_cutfem*=gamma_v_dim;
858 // //else
859 // // gamma_cutfem*=(gamma_v_dim/rho_0);
860 // for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
861 // {
862 // double Du_Dn_jump=0.0, dS;
863 // for (int eN_side=0;eN_side < 2; eN_side++)
864 // {
865 // int ebN = *it,
866 // eN = elementBoundaryElementsArray.data()[ebN*2+eN_side];
867 // for (int i=0;i<nDOF_test_element;i++)
868 // {
869 // DW_Dn_jump[r_l2g.data()[eN*nDOF_test_element+i]] = 0.0;
870 // }
871 // }
872 // for (int eN_side=0;eN_side < 2; eN_side++)
873 // {
874 // int ebN = *it,
875 // eN = elementBoundaryElementsArray.data()[ebN*2+eN_side],
876 // ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+eN_side],
877 // eN_nDOF_trial_element = eN*nDOF_trial_element,
878 // ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
879 // ebN_local_kb_nSpace = ebN_local_kb*nSpace;
880 // double u_int=0.0,
881 // grad_u_int[nSpace],
882 // jac_int[nSpace*nSpace],
883 // jacDet_int,
884 // jacInv_int[nSpace*nSpace],
885 // boundaryJac[nSpace*(nSpace-1)],
886 // metricTensor[(nSpace-1)*(nSpace-1)],
887 // metricTensorDetSqrt,
888 // u_test_dS[nDOF_test_element],
889 // u_grad_trial_trace[nDOF_trial_element*nSpace],
890 // u_grad_test_dS[nDOF_trial_element*nSpace],
891 // normal[nSpace],x_int,y_int,z_int,xt_int,yt_int,zt_int,integralScaling,
892 // G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty;
893 // for (int I=0; I<nSpace;I++)
894 // grad_u_int[I] = 0.0;
895 // //compute information about mapping from reference element to physical element
896
897 // ck.calculateMapping_elementBoundary(eN,
898 // ebN_local,
899 // kb,
900 // ebN_local_kb,
901 // mesh_dof.data(),
902 // mesh_l2g.data(),
903 // mesh_trial_trace_ref.data(),
904 // mesh_grad_trial_trace_ref.data(),
905 // boundaryJac_ref.data(),
906 // jac_int,
907 // jacDet_int,
908 // jacInv_int,
909 // boundaryJac,
910 // metricTensor,
911 // metricTensorDetSqrt,
912 // normal_ref.data(),
913 // normal,
914 // x_int,y_int,z_int);
915 // //todo: check that physical coordinates match
916 // ck.calculateMappingVelocity_elementBoundary(eN,
917 // ebN_local,
918 // kb,
919 // ebN_local_kb,
920 // mesh_velocity_dof.data(),
921 // mesh_l2g.data(),
922 // mesh_trial_trace_ref.data(),
923 // xt_int,yt_int,zt_int,
924 // normal,
925 // boundaryJac,
926 // metricTensor,
927 // integralScaling);
928 // dS = metricTensorDetSqrt*dS_ref.data()[kb];
929 // //compute shape and solution information
930 // //shape
931 // ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_int,u_grad_trial_trace);
932 // for (int i=0;i<nDOF_test_element;i++)
933 // {
934 // int eN_i = eN*nDOF_test_element + i;
935 // for (int I=0;I<nSpace;I++)
936 // DW_Dn_jump[r_l2g.data()[eN_i]] += u_grad_trial_trace[i*nSpace+I]*normal[I];
937 // }
938 // }//eN_side
939 // for (int eN_side=0;eN_side < 2; eN_side++)
940 // {
941 // int ebN = *it,
942 // eN = elementBoundaryElementsArray.data()[ebN*2+eN_side];
943 // for (int i=0;i<nDOF_test_element;i++)
944 // {
945 // int eN_i = eN*nDOF_test_element+i;
946 // for (int eN_side2=0;eN_side2 < 2; eN_side2++)
947 // {
948 // int eN2 = elementBoundaryElementsArray.data()[ebN*2+eN_side2];
949 // for (int j=0;j<nDOF_test_element;j++)
950 // {
951 // int eN_i_j = eN_i*nDOF_test_element + j;
952 // int eN2_j = eN2*nDOF_test_element + j;
953 // int ebN_i_j = ebN*4*nDOF_test_X_trial_element +
954 // eN_side*2*nDOF_test_X_trial_element +
955 // eN_side2*nDOF_test_X_trial_element +
956 // i*nDOF_trial_element +
957 // j;
958 // std::pair<int,int> ij = std::make_pair(u_l2g.data()[eN_i], u_l2g.data()[eN2_j]);
959 // if (u_u_nz.count(ij))
960 // {
961 // assert(u_u_nz[ij] == csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j]);
962 // }
963 // else
964 // u_u_nz[ij] = csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j];
965 // }
966 // }
967 // }
968 // }
969 // for (std::map<int,double>::iterator Wi_it=DW_Dn_jump.begin(); Wi_it!=DW_Dn_jump.end(); ++Wi_it)
970 // for (std::map<int,double>::iterator Wj_it=DW_Dn_jump.begin(); Wj_it!=DW_Dn_jump.end(); ++Wj_it)
971 // {
972 // int i_global = Wi_it->first,
973 // j_global = Wj_it->first;
974 // double DW_Dn_jump_i = Wi_it->second,
975 // DW_Dn_jump_j = Wj_it->second;
976 // std::pair<int,int> ij = std::make_pair(i_global, j_global);
977 // //globalJacobian.data()[u_u_nz.at(ij)] += gamma_cutfem*h_cutfem*DW_Dn_jump_j*DW_Dn_jump_i*dS;
978 // globalResidual.data()[offset_u+stride_u*i_global]+=u_dof.data()[j_global]*gamma_cutfem*h_cutfem*DW_Dn_jump_j*DW_Dn_jump_i*dS;
979 // }//i,j
980 // }//kb
981 // it++;
982 // }//cutfem element boundaries
983 //
984 //loop over exterior element boundaries to calculate surface integrals and load into element and global residuals
985 //
986 //ebNE is the Exterior element boundary INdex
987 //ebN is the element boundary INdex
988 //eN is the element index
989 for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
990 {
991 int ebN = exteriorElementBoundariesArray.data()[ebNE],
992 eN = elementBoundaryElementsArray.data()[ebN*2+0],
993 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
994 eN_nDOF_trial_element = eN*nDOF_trial_element;
995 double elementResidual_u[nDOF_test_element];
996 for (int i=0;i<nDOF_test_element;i++)
997 {
998 elementResidual_u[i]=0.0;
999 }
1000 double element_phi_s[nDOF_mesh_trial_element];
1001 for (int j=0;j<nDOF_mesh_trial_element;j++)
1002 {
1003 int eN_j = eN*nDOF_mesh_trial_element+j;
1004 element_phi_s[j] = phi_solid_nodes[u_l2g.data()[eN_j]];
1005 }
1006 double element_nodes[nDOF_mesh_trial_element*3];
1007 for (int i=0;i<nDOF_mesh_trial_element;i++)
1008 {
1009 int eN_i=eN*nDOF_mesh_trial_element+i;
1010 for(int I=0;I<3;I++)
1011 element_nodes[i*3 + I] = mesh_dof[mesh_l2g.data()[eN_i]*3 + I];
1012 }//i
1013 double mesh_dof_ref[12]={0.,0.,0.,1.,0.,0.,0.,1.,0.,0.,0.,1.};
1014 double xb_ref_calc[nQuadraturePoints_elementBoundary*3];
1015 for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
1016 {
1017 double x=0.0,y=0.0,z=0.0;
1018 for (int j=0;j<nDOF_mesh_trial_element;j++)
1019 {
1020 int ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb;
1021 int ebN_local_kb_j = ebN_local_kb*nDOF_mesh_trial_element+j;
1022 x += mesh_dof_ref[j*3+0]*mesh_trial_trace_ref.data()[ebN_local_kb_j];
1023 y += mesh_dof_ref[j*3+1]*mesh_trial_trace_ref.data()[ebN_local_kb_j];
1024 z += mesh_dof_ref[j*3+2]*mesh_trial_trace_ref.data()[ebN_local_kb_j];
1025 }
1026 xb_ref_calc[3*kb+0] = x;
1027 xb_ref_calc[3*kb+1] = y;
1028 xb_ref_calc[3*kb+2] = z;
1029 }
1030 int icase_s = gf_s.calculate(element_phi_s, element_nodes, xb_ref_calc, true);
1031 for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
1032 {
1033 int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
1034 ebNE_kb_nSpace = ebNE_kb*nSpace,
1035 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
1036 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
1037 double u_ext=0.0,
1038 grad_u_ext[nSpace],
1039 m_ext=0.0,
1040 dm_ext=0.0,
1041 f_ext[nSpace],
1042 df_ext[nSpace],
1043 flux_ext=0.0,
1044 bc_u_ext=0.0,
1045 //bc_grad_u_ext[nSpace],
1046 bc_m_ext=0.0,
1047 bc_dm_ext=0.0,
1048 bc_f_ext[nSpace],
1049 bc_df_ext[nSpace],
1050 jac_ext[nSpace*nSpace],
1051 jacDet_ext,
1052 jacInv_ext[nSpace*nSpace],
1053 boundaryJac[nSpace*(nSpace-1)],
1054 metricTensor[(nSpace-1)*(nSpace-1)],
1055 metricTensorDetSqrt,
1056 dS,
1057 u_test_dS[nDOF_test_element],
1058 u_grad_trial_trace[nDOF_trial_element*nSpace],
1059 normal[nSpace],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
1060 //VRANS
1061 porosity_ext,
1062 //
1063 G[nSpace*nSpace],G_dd_G,tr_G;
1064 //
1065 //calculate the solution and gradients at quadrature points
1066 //
1067 //compute information about mapping from reference element to physical element
1068 gf_s.set_boundary_quad(kb);
1069 ck.calculateMapping_elementBoundary(eN,
1070 ebN_local,
1071 kb,
1072 ebN_local_kb,
1073 mesh_dof.data(),
1074 mesh_l2g.data(),
1075 mesh_trial_trace_ref.data(),
1076 mesh_grad_trial_trace_ref.data(),
1077 boundaryJac_ref.data(),
1078 jac_ext,
1079 jacDet_ext,
1080 jacInv_ext,
1081 boundaryJac,
1082 metricTensor,
1083 metricTensorDetSqrt,
1084 normal_ref.data(),
1085 normal,
1086 x_ext,y_ext,z_ext);
1087 ck.calculateMappingVelocity_elementBoundary(eN,
1088 ebN_local,
1089 kb,
1090 ebN_local_kb,
1091 mesh_velocity_dof.data(),
1092 mesh_l2g.data(),
1093 mesh_trial_trace_ref.data(),
1094 xt_ext,yt_ext,zt_ext,
1095 normal,
1096 boundaryJac,
1097 metricTensor,
1098 integralScaling);
1099 //std::cout<<"metricTensorDetSqrt "<<metricTensorDetSqrt<<" integralScaling "<<integralScaling<<std::endl;
1100 dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref.data()[kb];
1101 //get the metric tensor
1102 //cek todo use symmetry
1103 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
1104 //compute shape and solution information
1105 //shape
1106 ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],
1107 jacInv_ext,
1108 u_grad_trial_trace);
1109 //solution and gradients
1110 if (STABILIZATION_TYPE==1) //explicit
1111 {
1112 ck.valFromDOF(u_dof_old.data(),
1113 &u_l2g.data()[eN_nDOF_trial_element],
1114 &u_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],
1115 u_ext);
1116 ck.gradFromDOF(u_dof_old.data(),
1117 &u_l2g.data()[eN_nDOF_trial_element],
1118 u_grad_trial_trace,
1119 grad_u_ext);
1120 }
1121 else
1122 {
1123 ck.valFromDOF(u_dof.data(),
1124 &u_l2g.data()[eN_nDOF_trial_element],
1125 &u_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],
1126 u_ext);
1127 ck.gradFromDOF(u_dof.data(),
1128 &u_l2g.data()[eN_nDOF_trial_element],
1129 u_grad_trial_trace,
1130 grad_u_ext);
1131 }
1132 //precalculate test function products with integration weights
1133 for (int j=0;j<nDOF_trial_element;j++)
1134 {
1135 u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
1136 }
1137 //
1138 //load the boundary values
1139 //
1140 bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+(1-isDOFBoundary_u.data()[ebNE_kb])*u_ext;
1141 //VRANS
1142 porosity_ext = ebqe_porosity_ext.data()[ebNE_kb];
1143 //
1144 //
1145 //calculate the pde coefficients using the solution and the boundary values for the solution
1146 //
1147 evaluateCoefficients(&ebqe_velocity_ext.data()[ebNE_kb_nSpace],
1148 u_ext,
1149 //VRANS
1150 porosity_ext,
1151 //
1152 m_ext,
1153 dm_ext,
1154 f_ext,
1155 df_ext);
1156 evaluateCoefficients(&ebqe_velocity_ext.data()[ebNE_kb_nSpace],
1157 bc_u_ext,
1158 //VRANS
1159 porosity_ext,
1160 //
1161 bc_m_ext,
1162 bc_dm_ext,
1163 bc_f_ext,
1164 bc_df_ext);
1165 //
1166 //moving mesh
1167 //
1168 double mesh_velocity[3];
1169 mesh_velocity[0] = xt_ext;
1170 mesh_velocity[1] = yt_ext;
1171 mesh_velocity[2] = zt_ext;
1172 //std::cout<<"mesh_velocity ext"<<std::endl;
1173 for (int I=0;I<nSpace;I++)
1174 {
1175 //std::cout<<mesh_velocity[I]<<std::endl;
1176 f_ext[I] -= MOVING_DOMAIN*m_ext*mesh_velocity[I];
1177 df_ext[I] -= MOVING_DOMAIN*dm_ext*mesh_velocity[I];
1178 bc_f_ext[I] -= MOVING_DOMAIN*bc_m_ext*mesh_velocity[I];
1179 bc_df_ext[I] -= MOVING_DOMAIN*bc_dm_ext*mesh_velocity[I];
1180 }
1181 //
1182 //calculate the numerical fluxes
1183 //
1184 exteriorNumericalAdvectiveFlux(isDOFBoundary_u.data()[ebNE_kb],
1185 isFluxBoundary_u.data()[ebNE_kb],
1186 normal,
1187 bc_u_ext,
1188 ebqe_bc_flux_u_ext.data()[ebNE_kb],
1189 u_ext,
1190 df_ext,//VRANS includes porosity
1191 flux_ext);
1192 ebqe_flux.data()[ebNE_kb] = flux_ext;
1193 //save for other models? cek need to be consistent with numerical flux
1194 if(flux_ext >=0.0)
1195 ebqe_u.data()[ebNE_kb] = u_ext;
1196 else
1197 ebqe_u.data()[ebNE_kb] = bc_u_ext;
1198
1199 if (STABILIZATION_TYPE==1)
1200 if (stage==1)
1201 flux_ext *= 1./3*dt;
1202 else
1203 flux_ext *= dt;
1204
1205 //
1206 //update residuals
1207 //
1208 const double H_s = gf_s.H(0.0, ebqe_phi_s.data()[ebNE_kb]);
1209 if (isActiveElement[eN])
1210 {
1211 for (int i=0;i<nDOF_test_element;i++)
1212 {
1213 //int ebNE_kb_i = ebNE_kb*nDOF_test_element+i;
1214 elementResidual_u[i] += H_s*ck.ExteriorElementBoundaryFlux(flux_ext,u_test_dS[i]);
1215 }//i
1216 }
1217 }//kb
1218 //
1219 //update the element and global residual storage
1220 //
1221 for (int i=0;i<nDOF_test_element;i++)
1222 {
1223 int eN_i = eN*nDOF_test_element+i;
1224 globalResidual.data()[offset_u+stride_u*r_l2g.data()[eN_i]] += elementResidual_u[i];
1225 }//i
1226 }//ebNE
1227 if (STABILIZATION_TYPE==1)
1228 {
1229 meanEntropy /= meanOmega;
1230 double norm_factor = fmax(fabs(maxEntropy - meanEntropy), fabs(meanEntropy-minEntropy));
1231 for(int eN=0;eN<nElements_global;eN++)
1232 {
1233 double hK=elementDiameter.data()[eN]/degree_polynomial;
1234 double linear_viscosity = cMax*hK*maxVel[eN];
1235 double entropy_viscosity = cE*hK*hK*maxEntRes[eN]/norm_factor;
1236 for (int k=0;k<nQuadraturePoints_element;k++)
1237 {
1238 int eN_k = eN*nQuadraturePoints_element+k;
1239 q_numDiff_u.data()[eN_k] = fmin(linear_viscosity,entropy_viscosity);
1240 }
1241 }
1242 }
1243
1244 }
1245
1246 void calculateJacobian(arguments_dict& args)
1247 {
1248 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1249 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1250 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1251 xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
1252 double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
1253 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1254 xt::pyarray<double>& x_ref = args.array<double>("x_ref");
1255 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1256 xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1257 xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1258 xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1259 xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
1260 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
1261 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
1262 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
1263 xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
1264 xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
1265 xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
1266 xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
1267 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
1268 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
1269 int nElements_global = args.scalar<int>("nElements_global");
1270 double useMetrics = args.scalar<double>("useMetrics");
1271 double alphaBDF = args.scalar<double>("alphaBDF");
1272 int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
1273 double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
1274 const xt::pyarray<double>& q_porosity = args.array<double>("q_porosity");
1275 xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1276 xt::pyarray<int>& r_l2g = args.array<int>("r_l2g");
1277 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1278 xt::pyarray<double>& elementBoundaryDiameter = args.array<double>("elementBoundaryDiameter");
1279 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1280 xt::pyarray<double>& velocity = args.array<double>("velocity");
1281 xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
1282 xt::pyarray<double>& cfl = args.array<double>("cfl");
1283 xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
1284 xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
1285 xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
1286 xt::pyarray<double>& globalJacobian = args.array<double>("globalJacobian");
1287 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
1288 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
1289 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
1290 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
1291 xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
1292 const xt::pyarray<double>& ebqe_porosity_ext = args.array<double>("ebqe_porosity_ext");
1293 xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
1294 xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
1295 xt::pyarray<int>& isFluxBoundary_u = args.array<int>("isFluxBoundary_u");
1296 xt::pyarray<double>& ebqe_bc_flux_u_ext = args.array<double>("ebqe_bc_flux_u_ext");
1297 xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
1298 int STABILIZATION_TYPE = args.scalar<int>("STABILIZATION_TYPE");
1299 xt::pyarray<double>& ebqe_phi_s = args.array<double>("ebqe_phi_s");
1300 const xt::pyarray<double>& phi_solid = args.array<double>("phi_solid");
1301 double ghost_penalty_constant = args.scalar<double>("ghost_penalty_constant");
1302 xt::pyarray<double>& phi_solid_nodes = args.array<double>("phi_solid_nodes");
1303 bool useExact = args.scalar<int>("useExact");
1304 xt::pyarray<double>& isActiveR = args.array<double>("isActiveR");
1305 xt::pyarray<double>& isActiveDOF = args.array<double>("isActiveDOF");
1306 xt::pyarray<int>& isActiveElement = args.array<int>("isActiveElement");
1307 //std::cout<<"ndjaco address "<<q_numDiff_u_last.data()<<std::endl;
1308 double Ct_sge = 4.0;
1309 gf_s.useExact = useExact;
1310 //
1311 //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
1312 //
1313 for(int eN=0;eN<nElements_global;eN++)
1314 {
1315 double elementJacobian_u_u[nDOF_test_element][nDOF_trial_element];
1316 for (int i=0;i<nDOF_test_element;i++)
1317 for (int j=0;j<nDOF_trial_element;j++)
1318 {
1319 elementJacobian_u_u[i][j]=0.0;
1320 }
1321 double element_phi_s[nDOF_mesh_trial_element];
1322 for (int j=0;j<nDOF_mesh_trial_element;j++)
1323 {
1324 int eN_j = eN*nDOF_mesh_trial_element+j;
1325 element_phi_s[j] = phi_solid_nodes.data()[u_l2g.data()[eN_j]];
1326 }
1327 double element_nodes[nDOF_mesh_trial_element*3];
1328 for (int i=0;i<nDOF_mesh_trial_element;i++)
1329 {
1330 int eN_i=eN*nDOF_mesh_trial_element+i;
1331 for(int I=0;I<3;I++)
1332 element_nodes[i*3 + I] = mesh_dof.data()[mesh_l2g.data()[eN_i]*3 + I];
1333 }//i
1334 int icase_s = gf_s.calculate(element_phi_s, element_nodes, x_ref.data(), false);
1335 for (int k=0;k<nQuadraturePoints_element;k++)
1336 {
1337 int eN_k = eN*nQuadraturePoints_element+k, //index to a scalar at a quadrature point
1338 eN_k_nSpace = eN_k*nSpace,
1339 eN_nDOF_trial_element = eN*nDOF_trial_element; //index to a vector at a quadrature point
1340
1341 //declare local storage
1342 double u=0.0,
1343 grad_u[nSpace],
1344 m=0.0,dm=0.0,
1345 f[nSpace],df[nSpace],
1346 m_t=0.0,dm_t=0.0,
1347 dpdeResidual_u_u[nDOF_trial_element],
1348 Lstar_u[nDOF_test_element],
1349 dsubgridError_u_u[nDOF_trial_element],
1350 tau=0.0,tau0=0.0,tau1=0.0,
1351 jac[nSpace*nSpace],
1352 jacDet,
1353 jacInv[nSpace*nSpace],
1354 u_grad_trial[nDOF_trial_element*nSpace],
1355 dV,
1356 u_test_dV[nDOF_test_element],
1357 u_grad_test_dV[nDOF_test_element*nSpace],
1358 x,y,z,xt,yt,zt,
1359 //VRANS
1360 porosity,
1361 //
1362 G[nSpace*nSpace],G_dd_G,tr_G;
1363 gf_s.set_quad(k);
1364 //
1365 //calculate solution and gradients at quadrature points
1366 //
1367 // u=0.0;
1368 // for (int I=0;I<nSpace;I++)
1369 // {
1370 // grad_u[I]=0.0;
1371 // }
1372 // for (int j=0;j<nDOF_trial_element;j++)
1373 // {
1374 // int eN_j=eN*nDOF_trial_element+j;
1375 // int eN_k_j=eN_k*nDOF_trial_element+j;
1376 // int eN_k_j_nSpace = eN_k_j*nSpace;
1377
1378 // u += valFromDOF_c(u_dof.data()[u_l2g.data()[eN_j]],u_trial[eN_k_j]);
1379 // for (int I=0;I<nSpace;I++)
1380 // {
1381 // grad_u[I] += gradFromDOF_c(u_dof.data()[u_l2g.data()[eN_j]],u_grad_trial[eN_k_j_nSpace+I]);
1382 // }
1383 // }
1384 //get jacobian, etc for mapping reference element
1385 ck.calculateMapping_element(eN,
1386 k,
1387 mesh_dof.data(),
1388 mesh_l2g.data(),
1389 mesh_trial_ref.data(),
1390 mesh_grad_trial_ref.data(),
1391 jac,
1392 jacDet,
1393 jacInv,
1394 x,y,z);
1395 ck.calculateMappingVelocity_element(eN,
1396 k,
1397 mesh_velocity_dof.data(),
1398 mesh_l2g.data(),
1399 mesh_trial_ref.data(),
1400 xt,yt,zt);
1401 //get the physical integration weight
1402 dV = fabs(jacDet)*dV_ref.data()[k];
1403 ck.calculateG(jacInv,G,G_dd_G,tr_G);
1404 //get the trial function gradients
1405 ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
1406 //get the solution
1407 ck.valFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],u);
1408 //get the solution gradients
1409 ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial,grad_u);
1410 //precalculate test function products with integration weights
1411 for (int j=0;j<nDOF_trial_element;j++)
1412 {
1413 u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
1414 for (int I=0;I<nSpace;I++)
1415 {
1416 u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
1417 }
1418 }
1419 //VRANS
1420 porosity = q_porosity.data()[eN_k];
1421 //
1422 //
1423 //calculate pde coefficients and derivatives at quadrature points
1424 //
1425 const double H_s = gf_s.H(0.0, phi_solid.data()[eN_k]);
1426 evaluateCoefficients(&velocity.data()[eN_k_nSpace],
1427 u,
1428 //VRANS
1429 porosity,
1430 //
1431 m,
1432 dm,
1433 f,
1434 df);
1435 //
1436 //moving mesh
1437 //
1438 double mesh_velocity[3];
1439 mesh_velocity[0] = xt;
1440 mesh_velocity[1] = yt;
1441 mesh_velocity[2] = zt;
1442 //std::cout<<"qj mesh_velocity"<<std::endl;
1443 for(int I=0;I<nSpace;I++)
1444 {
1445 //std::cout<<mesh_velocity[I]<<std::endl;
1446 f[I] -= MOVING_DOMAIN*m*mesh_velocity[I];
1447 df[I] -= MOVING_DOMAIN*dm*mesh_velocity[I];
1448 }
1449 //
1450 //calculate time derivatives
1451 //
1452 ck.bdf(alphaBDF,
1453 q_m_betaBDF.data()[eN_k],//since m_t isn't used, we don't have to correct mass
1454 m,
1455 dm,
1456 m_t,
1457 dm_t);
1458 //
1459 //calculate subgrid error contribution to the Jacobian (strong residual, adjoint, jacobian of strong residual)
1460 //
1461 //calculate the adjoint times the test functions
1462 for (int i=0;i<nDOF_test_element;i++)
1463 {
1464 // int eN_k_i_nSpace = (eN_k*nDOF_trial_element+i)*nSpace;
1465 // Lstar_u[i]=ck.Advection_adjoint(df,&u_grad_test_dV[eN_k_i_nSpace]);
1466 int i_nSpace = i*nSpace;
1467 Lstar_u[i]=ck.Advection_adjoint(df,&u_grad_test_dV[i_nSpace]);
1468 }
1469 //calculate the Jacobian of strong residual
1470 for (int j=0;j<nDOF_trial_element;j++)
1471 {
1472 //int eN_k_j=eN_k*nDOF_trial_element+j;
1473 //int eN_k_j_nSpace = eN_k_j*nSpace;
1474 int j_nSpace = j*nSpace;
1475 dpdeResidual_u_u[j]= ck.MassJacobian_strong(dm_t,u_trial_ref.data()[k*nDOF_trial_element+j]) +
1476 ck.AdvectionJacobian_strong(df,&u_grad_trial[j_nSpace]);
1477 }
1478 //tau and tau*Res
1479 calculateSubgridError_tau(elementDiameter.data()[eN],
1480 dm_t,
1481 df,
1482 cfl.data()[eN_k],
1483 tau0);
1484
1485 calculateSubgridError_tau(Ct_sge,
1486 G,
1487 dm_t,
1488 df,
1489 tau1,
1490 cfl.data()[eN_k]);
1491 tau = useMetrics*tau1+(1.0-useMetrics)*tau0;
1492
1493 for(int j=0;j<nDOF_trial_element;j++)
1494 dsubgridError_u_u[j] = -tau*dpdeResidual_u_u[j];
1495
1496 for(int i=0;i<nDOF_test_element;i++)
1497 {
1498 //int eN_k_i=eN_k*nDOF_test_element+i;
1499 //int eN_k_i_nSpace=eN_k_i*nSpace;
1500 for(int j=0;j<nDOF_trial_element;j++)
1501 {
1502 //int eN_k_j=eN_k*nDOF_trial_element+j;
1503 //int eN_k_j_nSpace = eN_k_j*nSpace;
1504 int j_nSpace = j*nSpace;
1505 int i_nSpace = i*nSpace;
1506 if (STABILIZATION_TYPE==0)
1507 {
1508 elementJacobian_u_u[i][j] +=
1509 H_s*(ck.MassJacobian_weak(dm_t,
1510 u_trial_ref.data()[k*nDOF_trial_element+j],
1511 u_test_dV[i]) +
1512 ck.AdvectionJacobian_weak(df,
1513 u_trial_ref.data()[k*nDOF_trial_element+j],
1514 &u_grad_test_dV[i_nSpace]) +
1515 ck.SubgridErrorJacobian(dsubgridError_u_u[j],Lstar_u[i]) +
1516 ck.NumericalDiffusionJacobian(q_numDiff_u_last.data()[eN_k],
1517 &u_grad_trial[j_nSpace],
1518 &u_grad_test_dV[i_nSpace])); //implicit
1519 }
1520 else
1521 {
1522 elementJacobian_u_u[i][j] +=
1523 ck.MassJacobian_weak(1.0,
1524 u_trial_ref.data()[k*nDOF_trial_element+j],
1525 u_test_dV[i]);
1526 }
1527 }//j
1528 }//i
1529 }//k
1530 //
1531 //load into element Jacobian into global Jacobian
1532 //
1533 for (int i=0;i<nDOF_test_element;i++)
1534 {
1535 int eN_i = eN*nDOF_test_element+i;
1536 for (int j=0;j<nDOF_trial_element;j++)
1537 {
1538 int eN_i_j = eN_i*nDOF_trial_element+j;
1539 globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_u_u.data()[eN_i_j]] += elementJacobian_u_u[i][j];
1540 }//j
1541 }//i
1542 }//elements
1543 std::set<int>::iterator it=cutfem_boundaries.begin();
1544 while(it!=cutfem_boundaries.end())
1545 {
1546 std::map<int,double> DW_Dn_jump;
1547 std::map<std::pair<int, int>, int> u_u_nz;
1548 double gamma_cutfem=ghost_penalty_constant,
1549 h_cutfem=elementBoundaryDiameter.data()[*it];
1550 int eN_nDOF_trial_element = elementBoundaryElementsArray.data()[(*it)*2+0]*nDOF_trial_element;
1551 //See Massing Schott Wall 2018
1552 //double norm_v=0.0;
1553 //for (int i_offset=1;i_offset<nDOF_v_trial_element;i_offset++)//MSW18 is just on face
1554 // {
1555 // int i = (cutfem_local_boundaries[*it] + i_offset)%nDOF_v_trial_element;//cek hack only works for P1
1556 // double u=u_old_dof.data()[vel_l2g.data()[eN_nDOF_v_trial_element+i]],
1557 // v=v_old_dof.data()[vel_l2g.data()[eN_nDOF_v_trial_element+i]],
1558 // w=w_old_dof.data()[vel_l2g.data()[eN_nDOF_v_trial_element+i]];
1559 // norm_v=fmax(norm_v,sqrt(u*u+v*v+w*w));
1560 // }
1562 //gamma_cutfem_p *= h_cutfem*h_cutfem/gamma_v_dim;
1563 //if (NONCONSERVATIVE_FORM)
1564 // gamma_cutfem*=gamma_v_dim;
1565 //else
1566 // gamma_cutfem*=(gamma_v_dim/rho_0);
1567 for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
1568 {
1569 double Du_Dn_jump=0.0, dS;
1570 for (int eN_side=0;eN_side < 2; eN_side++)
1571 {
1572 int ebN = *it,
1573 eN = elementBoundaryElementsArray.data()[ebN*2+eN_side];
1574 for (int i=0;i<nDOF_test_element;i++)
1575 {
1576 DW_Dn_jump[r_l2g.data()[eN*nDOF_test_element+i]] = 0.0;
1577 }
1578 }
1579 for (int eN_side=0;eN_side < 2; eN_side++)
1580 {
1581 int ebN = *it,
1582 eN = elementBoundaryElementsArray.data()[ebN*2+eN_side],
1583 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+eN_side],
1584 eN_nDOF_trial_element = eN*nDOF_trial_element,
1585 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
1586 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
1587 double u_int=0.0,
1588 grad_u_int[nSpace],
1589 jac_int[nSpace*nSpace],
1590 jacDet_int,
1591 jacInv_int[nSpace*nSpace],
1592 boundaryJac[nSpace*(nSpace-1)],
1593 metricTensor[(nSpace-1)*(nSpace-1)],
1594 metricTensorDetSqrt,
1595 u_test_dS[nDOF_test_element],
1596 u_grad_trial_trace[nDOF_trial_element*nSpace],
1597 u_grad_test_dS[nDOF_trial_element*nSpace],
1598 normal[nSpace],x_int,y_int,z_int,xt_int,yt_int,zt_int,integralScaling,
1599 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty;
1600 for (int I=0; I<nSpace;I++)
1601 grad_u_int[I] = 0.0;
1602 //compute information about mapping from reference element to physical element
1603
1604 ck.calculateMapping_elementBoundary(eN,
1605 ebN_local,
1606 kb,
1607 ebN_local_kb,
1608 mesh_dof.data(),
1609 mesh_l2g.data(),
1610 mesh_trial_trace_ref.data(),
1611 mesh_grad_trial_trace_ref.data(),
1612 boundaryJac_ref.data(),
1613 jac_int,
1614 jacDet_int,
1615 jacInv_int,
1616 boundaryJac,
1617 metricTensor,
1618 metricTensorDetSqrt,
1619 normal_ref.data(),
1620 normal,
1621 x_int,y_int,z_int);
1622 //todo: check that physical coordinates match
1623 ck.calculateMappingVelocity_elementBoundary(eN,
1624 ebN_local,
1625 kb,
1626 ebN_local_kb,
1627 mesh_velocity_dof.data(),
1628 mesh_l2g.data(),
1629 mesh_trial_trace_ref.data(),
1630 xt_int,yt_int,zt_int,
1631 normal,
1632 boundaryJac,
1633 metricTensor,
1634 integralScaling);
1635 dS = metricTensorDetSqrt*dS_ref.data()[kb];
1636 //compute shape and solution information
1637 //shape
1638 ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_int,u_grad_trial_trace);
1639 for (int i=0;i<nDOF_test_element;i++)
1640 {
1641 int eN_i = eN*nDOF_test_element + i;
1642 for (int I=0;I<nSpace;I++)
1643 DW_Dn_jump[r_l2g.data()[eN_i]] += u_grad_trial_trace[i*nSpace+I]*normal[I];
1644 }
1645 }//eN_side
1646 for (int eN_side=0;eN_side < 2; eN_side++)
1647 {
1648 int ebN = *it,
1649 eN = elementBoundaryElementsArray.data()[ebN*2+eN_side];
1650 for (int i=0;i<nDOF_test_element;i++)
1651 {
1652 int eN_i = eN*nDOF_test_element+i;
1653 for (int eN_side2=0;eN_side2 < 2; eN_side2++)
1654 {
1655 int eN2 = elementBoundaryElementsArray.data()[ebN*2+eN_side2];
1656 for (int j=0;j<nDOF_test_element;j++)
1657 {
1658 int eN_i_j = eN_i*nDOF_test_element + j;
1659 int eN2_j = eN2*nDOF_test_element + j;
1660 int ebN_i_j = ebN*4*nDOF_test_X_trial_element +
1661 eN_side*2*nDOF_test_X_trial_element +
1662 eN_side2*nDOF_test_X_trial_element +
1663 i*nDOF_trial_element +
1664 j;
1665 std::pair<int,int> ij = std::make_pair(u_l2g.data()[eN_i], u_l2g.data()[eN2_j]);
1666 if (u_u_nz.count(ij))
1667 {
1668 assert(u_u_nz[ij] == csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j]);
1669 }
1670 else
1671 u_u_nz[ij] = csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j];
1672 }
1673 }
1674 }
1675 }
1676 for (std::map<int,double>::iterator Wi_it=DW_Dn_jump.begin(); Wi_it!=DW_Dn_jump.end(); ++Wi_it)
1677 for (std::map<int,double>::iterator Wj_it=DW_Dn_jump.begin(); Wj_it!=DW_Dn_jump.end(); ++Wj_it)
1678 {
1679 int i_global = Wi_it->first,
1680 j_global = Wj_it->first;
1681 double DW_Dn_jump_i = Wi_it->second,
1682 DW_Dn_jump_j = Wj_it->second;
1683 std::pair<int,int> ij = std::make_pair(i_global, j_global);
1684 globalJacobian.data()[u_u_nz.at(ij)] += gamma_cutfem*h_cutfem*DW_Dn_jump_j*DW_Dn_jump_i*dS;
1685 }//i,j
1686 }//kb
1687 it++;
1688 }//cutfem element boundaries
1689 //
1690 //loop over exterior element boundaries to compute the surface integrals and load them into the global Jacobian
1691 //
1692 if (STABILIZATION_TYPE==0)
1693 for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
1694 {
1695 int ebN = exteriorElementBoundariesArray.data()[ebNE];
1696 int eN = elementBoundaryElementsArray.data()[ebN*2+0],
1697 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
1698 eN_nDOF_trial_element = eN*nDOF_trial_element;
1699 double element_phi_s[nDOF_mesh_trial_element];
1700 for (int j=0;j<nDOF_mesh_trial_element;j++)
1701 {
1702 int eN_j = eN*nDOF_mesh_trial_element+j;
1703 element_phi_s[j] = phi_solid_nodes.data()[u_l2g.data()[eN_j]];
1704 }
1705 double element_nodes[nDOF_mesh_trial_element*3];
1706 for (int i=0;i<nDOF_mesh_trial_element;i++)
1707 {
1708 int eN_i=eN*nDOF_mesh_trial_element+i;
1709 for(int I=0;I<3;I++)
1710 element_nodes[i*3 + I] = mesh_dof[mesh_l2g.data()[eN_i]*3 + I];
1711 }//i
1712 double mesh_dof_ref[12]={0.,0.,0.,1.,0.,0.,0.,1.,0.,0.,0.,1.};
1713 double xb_ref_calc[nQuadraturePoints_elementBoundary*3];
1714 for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
1715 {
1716 double x=0.0,y=0.0,z=0.0;
1717 for (int j=0;j<nDOF_mesh_trial_element;j++)
1718 {
1719 int ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb;
1720 int ebN_local_kb_j = ebN_local_kb*nDOF_mesh_trial_element+j;
1721 x += mesh_dof_ref[j*3+0]*mesh_trial_trace_ref.data()[ebN_local_kb_j];
1722 y += mesh_dof_ref[j*3+1]*mesh_trial_trace_ref.data()[ebN_local_kb_j];
1723 z += mesh_dof_ref[j*3+2]*mesh_trial_trace_ref.data()[ebN_local_kb_j];
1724 }
1725 xb_ref_calc[3*kb+0] = x;
1726 xb_ref_calc[3*kb+1] = y;
1727 xb_ref_calc[3*kb+2] = z;
1728 }
1729 int icase_s = gf_s.calculate(element_phi_s, element_nodes, xb_ref_calc,true);
1730 for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
1731 {
1732 int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
1733 ebNE_kb_nSpace = ebNE_kb*nSpace,
1734 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
1735 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
1736 double u_ext=0.0,
1737 grad_u_ext[nSpace],
1738 m_ext=0.0,
1739 dm_ext=0.0,
1740 f_ext[nSpace],
1741 df_ext[nSpace],
1742 dflux_u_u_ext=0.0,
1743 bc_u_ext=0.0,
1744 //bc_grad_u_ext[nSpace],
1745 bc_m_ext=0.0,
1746 bc_dm_ext=0.0,
1747 bc_f_ext[nSpace],
1748 bc_df_ext[nSpace],
1749 fluxJacobian_u_u[nDOF_trial_element],
1750 jac_ext[nSpace*nSpace],
1751 jacDet_ext,
1752 jacInv_ext[nSpace*nSpace],
1753 boundaryJac[nSpace*(nSpace-1)],
1754 metricTensor[(nSpace-1)*(nSpace-1)],
1755 metricTensorDetSqrt,
1756 dS,
1757 u_test_dS[nDOF_test_element],
1758 u_grad_trial_trace[nDOF_trial_element*nSpace],
1759 normal[nSpace],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
1760 //VRANS
1761 porosity_ext,
1762 //
1763 G[nSpace*nSpace],G_dd_G,tr_G;
1764 gf_s.set_boundary_quad(kb);
1765 //
1766 //calculate the solution and gradients at quadrature points
1767 //
1768 // u_ext=0.0;
1769 // for (int I=0;I<nSpace;I++)
1770 // {
1771 // grad_u_ext[I] = 0.0;
1772 // bc_grad_u_ext[I] = 0.0;
1773 // }
1774 // for (int j=0;j<nDOF_trial_element;j++)
1775 // {
1776 // int eN_j = eN*nDOF_trial_element+j,
1777 // ebNE_kb_j = ebNE_kb*nDOF_trial_element+j,
1778 // ebNE_kb_j_nSpace= ebNE_kb_j*nSpace;
1779 // u_ext += valFromDOF_c(u_dof.data()[u_l2g.data()[eN_j]],u_trial_ext[ebNE_kb_j]);
1780
1781 // for (int I=0;I<nSpace;I++)
1782 // {
1783 // grad_u_ext[I] += gradFromDOF_c(u_dof.data()[u_l2g.data()[eN_j]],u_grad_trial_ext[ebNE_kb_j_nSpace+I]);
1784 // }
1785 // }
1786 ck.calculateMapping_elementBoundary(eN,
1787 ebN_local,
1788 kb,
1789 ebN_local_kb,
1790 mesh_dof.data(),
1791 mesh_l2g.data(),
1792 mesh_trial_trace_ref.data(),
1793 mesh_grad_trial_trace_ref.data(),
1794 boundaryJac_ref.data(),
1795 jac_ext,
1796 jacDet_ext,
1797 jacInv_ext,
1798 boundaryJac,
1799 metricTensor,
1800 metricTensorDetSqrt,
1801 normal_ref.data(),
1802 normal,
1803 x_ext,y_ext,z_ext);
1804 ck.calculateMappingVelocity_elementBoundary(eN,
1805 ebN_local,
1806 kb,
1807 ebN_local_kb,
1808 mesh_velocity_dof.data(),
1809 mesh_l2g.data(),
1810 mesh_trial_trace_ref.data(),
1811 xt_ext,yt_ext,zt_ext,
1812 normal,
1813 boundaryJac,
1814 metricTensor,
1815 integralScaling);
1816 //std::cout<<"J mtsqrdet "<<metricTensorDetSqrt<<" integralScaling "<<integralScaling<<std::endl;
1817 dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref.data()[kb];
1818 //dS = metricTensorDetSqrt*dS_ref.data()[kb];
1819 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
1820 //compute shape and solution information
1821 //shape
1822 ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,u_grad_trial_trace);
1823 //solution and gradients
1824 ck.valFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],u_ext);
1825 ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial_trace,grad_u_ext);
1826 //precalculate test function products with integration weights
1827 for (int j=0;j<nDOF_trial_element;j++)
1828 {
1829 u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
1830 }
1831 //
1832 //load the boundary values
1833 //
1834 bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+(1-isDOFBoundary_u.data()[ebNE_kb])*u_ext;
1835 //VRANS
1836 porosity_ext = ebqe_porosity_ext.data()[ebNE_kb];
1837 //
1838 //
1839 //calculate the internal and external trace of the pde coefficients
1840 //
1841 evaluateCoefficients(&ebqe_velocity_ext.data()[ebNE_kb_nSpace],
1842 u_ext,
1843 //VRANS
1844 porosity_ext,
1845 //
1846 m_ext,
1847 dm_ext,
1848 f_ext,
1849 df_ext);
1850 evaluateCoefficients(&ebqe_velocity_ext.data()[ebNE_kb_nSpace],
1851 bc_u_ext,
1852 //VRANS
1853 porosity_ext,
1854 //
1855 bc_m_ext,
1856 bc_dm_ext,
1857 bc_f_ext,
1858 bc_df_ext);
1859 //
1860 //moving domain
1861 //
1862 double mesh_velocity[3];
1863 mesh_velocity[0] = xt_ext;
1864 mesh_velocity[1] = yt_ext;
1865 mesh_velocity[2] = zt_ext;
1866 //std::cout<<"ext J mesh_velocity"<<std::endl;
1867 for (int I=0;I<nSpace;I++)
1868 {
1869 //std::cout<<mesh_velocity[I]<<std::endl;
1870 f_ext[I] -= MOVING_DOMAIN*m_ext*mesh_velocity[I];
1871 df_ext[I] -= MOVING_DOMAIN*dm_ext*mesh_velocity[I];
1872 bc_f_ext[I] -= MOVING_DOMAIN*bc_m_ext*mesh_velocity[I];
1873 bc_df_ext[I] -= MOVING_DOMAIN*bc_dm_ext*mesh_velocity[I];
1874 }
1875 //
1876 //calculate the numerical fluxes
1877 //
1878 exteriorNumericalAdvectiveFluxDerivative(isDOFBoundary_u.data()[ebNE_kb],
1879 isFluxBoundary_u.data()[ebNE_kb],
1880 normal,
1881 df_ext,//VRANS holds porosity
1882 dflux_u_u_ext);
1883 //
1884 //calculate the flux jacobian
1885 //
1886 for (int j=0;j<nDOF_trial_element;j++)
1887 {
1888 //int ebNE_kb_j = ebNE_kb*nDOF_trial_element+j;
1889 int ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j;
1890 fluxJacobian_u_u[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_u_u_ext,u_trial_trace_ref.data()[ebN_local_kb_j]);
1891 }//j
1892 //
1893 //update the global Jacobian from the flux Jacobian
1894 //
1895 const double H_s = gf_s.H(0.0, ebqe_phi_s[ebNE_kb]);
1896 if (isActiveElement[eN])
1897 {
1898 for (int i=0;i<nDOF_test_element;i++)
1899 {
1900 int eN_i = eN*nDOF_test_element+i;
1901 //int ebNE_kb_i = ebNE_kb*nDOF_test_element+i;
1902 for (int j=0;j<nDOF_trial_element;j++)
1903 {
1904 int ebN_i_j = ebN*4*nDOF_test_X_trial_element + i*nDOF_trial_element + j;
1905 globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j]] += H_s*fluxJacobian_u_u[j]*u_test_dS[i];
1906 }//j
1907 }//i
1908 }
1909 }//kb
1910 }//ebNE
1911 }//computeJacobian
1912
1913 void FCTStep(arguments_dict& args)
1914 {
1915 double dt = args.scalar<double>("dt");
1916 int NNZ = args.scalar<int>("NNZ");
1917 int numDOFs = args.scalar<int>("numDOFs");
1918 xt::pyarray<double>& lumped_mass_matrix = args.array<double>("lumped_mass_matrix");
1919 xt::pyarray<double>& soln = args.array<double>("soln");
1920 xt::pyarray<double>& solH = args.array<double>("solH");
1921 xt::pyarray<double>& uLow = args.array<double>("uLow");
1922 xt::pyarray<double>& dLow = args.array<double>("dLow");
1923 xt::pyarray<double>& limited_solution = args.array<double>("limited_solution");
1924 xt::pyarray<int>& csrRowIndeces_DofLoops = args.array<int>("csrRowIndeces_DofLoops");
1925 xt::pyarray<int>& csrColumnOffsets_DofLoops = args.array<int>("csrColumnOffsets_DofLoops");
1926 xt::pyarray<double>& MassMatrix = args.array<double>("MassMatrix");
1927 xt::pyarray<double>& dt_times_dH_minus_dL = args.array<double>("dt_times_dH_minus_dL");
1928 xt::pyarray<double>& min_u_bc = args.array<double>("min_u_bc");
1929 xt::pyarray<double>& max_u_bc = args.array<double>("max_u_bc");
1930 int LUMPED_MASS_MATRIX = args.scalar<int>("LUMPED_MASS_MATRIX");
1931 int STABILIZATION_TYPE = args.scalar<int>("STABILIZATION_TYPE");
1932 Rpos.resize(numDOFs,0.0);
1933 Rneg.resize(numDOFs,0.0);
1934 FluxCorrectionMatrix.resize(NNZ,0.0);
1936 // LOOP in DOFs //
1938 int ij=0;
1939 for (int i=0; i<numDOFs; i++)
1940 {
1941 //read some vectors
1942 double solHi = solH.data()[i];
1943 double solni = soln.data()[i];
1944 double mi = lumped_mass_matrix.data()[i];
1945 double uLowi = uLow.data()[i];
1946 double uDotLowi = (uLowi - solni)/dt;
1947 double mini=min_u_bc.data()[i], maxi=max_u_bc.data()[i]; // init min/max with value at BCs (NOTE: if no boundary then min=1E10, max=-1E10)
1948 if (GLOBAL_FCT==1)
1949 {
1950 mini = 0.;
1951 maxi = 1.;
1952 }
1953
1954 double Pposi=0, Pnegi=0;
1955 // LOOP OVER THE SPARSITY PATTERN (j-LOOP)//
1956 for (int offset=csrRowIndeces_DofLoops.data()[i]; offset<csrRowIndeces_DofLoops.data()[i+1]; offset++)
1957 {
1958 int j = csrColumnOffsets_DofLoops.data()[offset];
1959 double solnj = soln.data()[j];
1961 // COMPUTE THE BOUNDS //
1963 if (GLOBAL_FCT == 0)
1964 {
1965 mini = fmin(mini,solnj);
1966 maxi = fmax(maxi,solnj);
1967 }
1968 double uLowj = uLow.data()[j];
1969 double uDotLowj = (uLowj - solnj)/dt;
1970 // i-th row of flux correction matrix
1971 if (STABILIZATION_TYPE==4) // DK high-order, linearly stable anti-dif. flux
1972 {
1973 FluxCorrectionMatrix[ij] = dt*(MassMatrix.data()[ij]*(uDotLowi-uDotLowj)
1974 + dLow.data()[ij]*(uLowi-uLowj));
1975 }
1976 else
1977 {
1978 double ML_minus_MC =
1979 (LUMPED_MASS_MATRIX == 1 ? 0. : (i==j ? 1. : 0.)*mi - MassMatrix.data()[ij]);
1980 FluxCorrectionMatrix[ij] = ML_minus_MC * (solH.data()[j]-solnj - (solHi-solni))
1981 + dt_times_dH_minus_dL.data()[ij]*(solnj-solni);
1982 }
1983
1985 // COMPUTE P VECTORS //
1987 Pposi += FluxCorrectionMatrix[ij]*((FluxCorrectionMatrix[ij] > 0) ? 1. : 0.);
1988 Pnegi += FluxCorrectionMatrix[ij]*((FluxCorrectionMatrix[ij] < 0) ? 1. : 0.);
1989
1990 //update ij
1991 ij+=1;
1992 }
1994 // COMPUTE Q VECTORS //
1996 double Qposi = mi*(maxi-uLow.data()[i]);
1997 double Qnegi = mi*(mini-uLow.data()[i]);
1998
2000 // COMPUTE R VECTORS //
2002 Rpos[i] = ((Pposi==0) ? 1. : fmin(1.0,Qposi/Pposi));
2003 Rneg[i] = ((Pnegi==0) ? 1. : fmin(1.0,Qnegi/Pnegi));
2004 } // i DOFs
2005
2007 // COMPUTE LIMITERS //
2009 ij=0;
2010 for (int i=0; i<numDOFs; i++)
2011 {
2012 double ith_Limiter_times_FluxCorrectionMatrix = 0.;
2013 double Rposi = Rpos[i], Rnegi = Rneg[i];
2014 // LOOP OVER THE SPARSITY PATTERN (j-LOOP)//
2015 for (int offset=csrRowIndeces_DofLoops.data()[i]; offset<csrRowIndeces_DofLoops.data()[i+1]; offset++)
2016 {
2017 int j = csrColumnOffsets_DofLoops.data()[offset];
2018 double Lij = 1;
2019 Lij = ((FluxCorrectionMatrix[ij]>0) ? fmin(Rposi,Rneg[j]) : fmin(Rnegi,Rpos[j]));
2020 ith_Limiter_times_FluxCorrectionMatrix += Lij * FluxCorrectionMatrix[ij];
2021 //update ij
2022 ij+=1;
2023 }
2024 limited_solution.data()[i] = uLow.data()[i] + 1./lumped_mass_matrix.data()[i]*ith_Limiter_times_FluxCorrectionMatrix;
2025 }
2026 }
2027
2028 void calculateResidualEdgeBased(arguments_dict& args)
2029 {
2030 double dt = args.scalar<double>("dt");
2031 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
2032 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
2033 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
2034 xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
2035 double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
2036 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
2037 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
2038 xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
2039 xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
2040 xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
2041 xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
2042 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
2043 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
2044 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
2045 xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
2046 xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
2047 xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
2048 xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
2049 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
2050 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
2051 int nElements_global = args.scalar<int>("nElements_global");
2052 double useMetrics = args.scalar<double>("useMetrics");
2053 double alphaBDF = args.scalar<double>("alphaBDF");
2054 int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
2055 double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
2056 double sc_uref = args.scalar<double>("sc_uref");
2057 double sc_alpha = args.scalar<double>("sc_alpha");
2058 const xt::pyarray<double>& q_porosity = args.array<double>("q_porosity");
2059 const xt::pyarray<double>& porosity_dof = args.array<double>("porosity_dof");
2060 xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
2061 xt::pyarray<int>& r_l2g = args.array<int>("r_l2g");
2062 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
2063 double degree_polynomial = args.scalar<double>("degree_polynomial");
2064 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
2065 xt::pyarray<double>& u_dof_old = args.array<double>("u_dof_old");
2066 xt::pyarray<double>& velocity = args.array<double>("velocity");
2067 xt::pyarray<double>& q_m = args.array<double>("q_m");
2068 xt::pyarray<double>& q_u = args.array<double>("q_u");
2069 xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
2070 xt::pyarray<double>& q_dV = args.array<double>("q_dV");
2071 xt::pyarray<double>& q_dV_last = args.array<double>("q_dV_last");
2072 xt::pyarray<double>& cfl = args.array<double>("cfl");
2073 xt::pyarray<double>& edge_based_cfl = args.array<double>("edge_based_cfl");
2074 xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
2075 xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
2076 int offset_u = args.scalar<int>("offset_u");
2077 int stride_u = args.scalar<int>("stride_u");
2078 xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
2079 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
2080 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
2081 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
2082 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
2083 xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
2084 const xt::pyarray<double>& ebqe_porosity_ext = args.array<double>("ebqe_porosity_ext");
2085 xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
2086 xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
2087 xt::pyarray<int>& isFluxBoundary_u = args.array<int>("isFluxBoundary_u");
2088 xt::pyarray<double>& ebqe_bc_flux_u_ext = args.array<double>("ebqe_bc_flux_u_ext");
2089 xt::pyarray<double>& ebqe_phi = args.array<double>("ebqe_phi");
2090 double epsFact = args.scalar<double>("epsFact");
2091 xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
2092 xt::pyarray<double>& ebqe_flux = args.array<double>("ebqe_flux");
2093 int stage = args.scalar<int>("stage");
2094 xt::pyarray<double>& uTilde_dof = args.array<double>("uTilde_dof");
2095 double cE = args.scalar<double>("cE");
2096 double cMax = args.scalar<double>("cMax");
2097 double cK = args.scalar<double>("cK");
2098 double uL = args.scalar<double>("uL");
2099 double uR = args.scalar<double>("uR");
2100 int numDOFs = args.scalar<int>("numDOFs");
2101 int NNZ = args.scalar<int>("NNZ");
2102 xt::pyarray<int>& csrRowIndeces_DofLoops = args.array<int>("csrRowIndeces_DofLoops");
2103 xt::pyarray<int>& csrColumnOffsets_DofLoops = args.array<int>("csrColumnOffsets_DofLoops");
2104 xt::pyarray<int>& csrRowIndeces_CellLoops = args.array<int>("csrRowIndeces_CellLoops");
2105 xt::pyarray<int>& csrColumnOffsets_CellLoops = args.array<int>("csrColumnOffsets_CellLoops");
2106 xt::pyarray<int>& csrColumnOffsets_eb_CellLoops = args.array<int>("csrColumnOffsets_eb_CellLoops");
2107 xt::pyarray<double>& ML = args.array<double>("ML");
2108 int LUMPED_MASS_MATRIX = args.scalar<int>("LUMPED_MASS_MATRIX");
2109 int STABILIZATION_TYPE = args.scalar<int>("STABILIZATION_TYPE");
2110 int ENTROPY_TYPE = args.scalar<int>("ENTROPY_TYPE");
2111 xt::pyarray<double>& uLow = args.array<double>("uLow");
2112 xt::pyarray<double>& dLow = args.array<double>("dLow");
2113 xt::pyarray<double>& dt_times_dH_minus_dL = args.array<double>("dt_times_dH_minus_dL");
2114 xt::pyarray<double>& min_u_bc = args.array<double>("min_u_bc");
2115 xt::pyarray<double>& max_u_bc = args.array<double>("max_u_bc");
2116 xt::pyarray<double>& quantDOFs = args.array<double>("quantDOFs");
2117 // NOTE: This function follows a different (but equivalent) implementation of the smoothness based indicator than NCLS.h
2118 // Allocate space for the transport matrices
2119 // This is used for first order KUZMIN'S METHOD
2120 TransportMatrix.resize(NNZ,0.0);
2121 TransposeTransportMatrix.resize(NNZ,0.0);
2122 // compute entropy and init global_entropy_residual and boundary_integral
2123 psi.resize(numDOFs,0.0);
2124 eta.resize(numDOFs,0.0);
2125 global_entropy_residual.resize(numDOFs,0.0);
2126 boundary_integral.resize(numDOFs,0.0);
2127
2128 for (int i=0; i<numDOFs; i++)
2129 {
2130 // NODAL ENTROPY //
2131 if (STABILIZATION_TYPE==2) //EV stab
2132 {
2133 double porosity_times_solni = porosity_dof.data()[i]*u_dof_old.data()[i];
2134 eta[i] = ENTROPY_TYPE == 0 ? ENTROPY(porosity_times_solni,uL,uR) : ENTROPY_LOG(porosity_times_solni,uL,uR);
2135 global_entropy_residual[i]=0.;
2136 }
2137 boundary_integral[i]=0.;
2138 }
2139
2141 // ** LOOP IN CELLS FOR CELL BASED TERMS ** //
2143 // HERE WE COMPUTE:
2144 // * Time derivative term. porosity*u_t
2145 // * cell based CFL (for reference)
2146 // * Entropy residual
2147 // * Transport matrices
2148 for(int eN=0;eN<nElements_global;eN++)
2149 {
2150 //declare local storage for local contributions and initialize
2151 double
2152 elementResidual_u[nDOF_test_element],
2153 element_entropy_residual[nDOF_test_element];
2154 double elementTransport[nDOF_test_element][nDOF_trial_element];
2155 double elementTransposeTransport[nDOF_test_element][nDOF_trial_element];
2156 for (int i=0;i<nDOF_test_element;i++)
2157 {
2158 elementResidual_u[i]=0.0;
2159 element_entropy_residual[i]=0.0;
2160 for (int j=0;j<nDOF_trial_element;j++)
2161 {
2162 elementTransport[i][j]=0.0;
2163 elementTransposeTransport[i][j]=0.0;
2164 }
2165 }
2166 //loop over quadrature points and compute integrands
2167 for (int k=0;k<nQuadraturePoints_element;k++)
2168 {
2169 //compute indeces and declare local storage
2170 int eN_k = eN*nQuadraturePoints_element+k,
2171 eN_k_nSpace = eN_k*nSpace,
2172 eN_nDOF_trial_element = eN*nDOF_trial_element;
2173 double
2174 // for entropy residual
2175 aux_entropy_residual=0., DENTROPY_un, DENTROPY_uni,
2176 //for mass matrix contributions
2177 u=0.0, un=0.0, grad_un[nSpace], porosity_times_velocity[nSpace],
2178 u_test_dV[nDOF_trial_element],
2179 u_grad_trial[nDOF_trial_element*nSpace],
2180 u_grad_test_dV[nDOF_test_element*nSpace],
2181 //for general use
2182 jac[nSpace*nSpace], jacDet, jacInv[nSpace*nSpace],
2183 dV,x,y,z,xt,yt,zt,
2184 //VRANS
2185 porosity;
2186 //get the physical integration weight
2187 ck.calculateMapping_element(eN,
2188 k,
2189 mesh_dof.data(),
2190 mesh_l2g.data(),
2191 mesh_trial_ref.data(),
2192 mesh_grad_trial_ref.data(),
2193 jac,
2194 jacDet,
2195 jacInv,
2196 x,y,z);
2197 ck.calculateMappingVelocity_element(eN,
2198 k,
2199 mesh_velocity_dof.data(),
2200 mesh_l2g.data(),
2201 mesh_trial_ref.data(),
2202 xt,yt,zt);
2203 dV = fabs(jacDet)*dV_ref.data()[k];
2204 //get the solution (of Newton's solver). To compute time derivative term
2205 ck.valFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],u);
2206 //get the solution at quad point at tn and tnm1 for entropy viscosity
2207 ck.valFromDOF(u_dof_old.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],un);
2208 //get the solution gradients at tn for entropy viscosity
2209 ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
2210 ck.gradFromDOF(u_dof_old.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial,grad_un);
2211
2212 //precalculate test function products with integration weights for mass matrix terms
2213 for (int j=0;j<nDOF_trial_element;j++)
2214 {
2215 u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
2216 for (int I=0;I<nSpace;I++)
2217 u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
2218 }
2219
2220 //calculate time derivative at quadrature points
2221 if (q_dV_last.data()[eN_k] <= -100)
2222 q_dV_last.data()[eN_k] = dV;
2223 q_dV.data()[eN_k] = dV;
2224 //VRANS
2225 porosity = q_porosity.data()[eN_k];
2226 //
2227 //moving mesh
2228 //
2229 double mesh_velocity[3];
2230 mesh_velocity[0] = xt;
2231 mesh_velocity[1] = yt;
2232 mesh_velocity[2] = zt;
2233 //relative velocity at tn
2234 for (int I=0;I<nSpace;I++)
2235 porosity_times_velocity[I] = porosity*(velocity.data()[eN_k_nSpace+I]-MOVING_DOMAIN*mesh_velocity[I]);
2236
2238 // CALCULATE CELL BASED CFL //
2240 calculateCFL(elementDiameter.data()[eN]/degree_polynomial,porosity_times_velocity,cfl.data()[eN_k]);
2241
2243 // CALCULATE ENTROPY RESIDUAL AT QUAD POINT //
2245 if (STABILIZATION_TYPE==2) // EV stab
2246 {
2247 for (int I=0;I<nSpace;I++)
2248 aux_entropy_residual += porosity_times_velocity[I]*grad_un[I];
2249 DENTROPY_un = ENTROPY_TYPE==0 ? DENTROPY(porosity*un,uL,uR) : DENTROPY_LOG(porosity*un,uL,uR);
2250 }
2252 // ith-LOOP //
2254 for(int i=0;i<nDOF_test_element;i++)
2255 {
2256 // VECTOR OF ENTROPY RESIDUAL //
2257 int eN_i=eN*nDOF_test_element+i;
2258 if (STABILIZATION_TYPE==2) // EV stab
2259 {
2260 int gi = offset_u+stride_u*u_l2g.data()[eN_i]; //global i-th index
2261 double porosity_times_uni = porosity_dof.data()[gi]*u_dof_old.data()[gi];
2262 DENTROPY_uni = ENTROPY_TYPE == 0 ? DENTROPY(porosity_times_uni,uL,uR) : DENTROPY_LOG(porosity_times_uni,uL,uR);
2263 element_entropy_residual[i] += (DENTROPY_un - DENTROPY_uni)*aux_entropy_residual*u_test_dV[i];
2264 }
2265 elementResidual_u[i] += porosity*(u-un)*u_test_dV[i];
2267 // j-th LOOP // To construct transport matrices
2269 for(int j=0;j<nDOF_trial_element;j++)
2270 {
2271 int j_nSpace = j*nSpace;
2272 int i_nSpace = i*nSpace;
2273 elementTransport[i][j] += // -int[(vel.grad_wi)*wj*dx]
2274 ck.AdvectionJacobian_weak(porosity_times_velocity,
2275 u_trial_ref.data()[k*nDOF_trial_element+j],&u_grad_test_dV[i_nSpace]);
2276 elementTransposeTransport[i][j] += // -int[(vel.grad_wj)*wi*dx]
2277 ck.AdvectionJacobian_weak(porosity_times_velocity,
2278 u_trial_ref.data()[k*nDOF_trial_element+i],&u_grad_test_dV[j_nSpace]);
2279 }
2280 }//i
2281 //save solution for other models
2282 q_u.data()[eN_k] = u;
2283 q_m.data()[eN_k] = porosity*u;
2284 }
2286 // DISTRIBUTE // load cell based element into global residual
2288 for(int i=0;i<nDOF_test_element;i++)
2289 {
2290 int eN_i=eN*nDOF_test_element+i;
2291 int gi = offset_u+stride_u*u_l2g.data()[eN_i]; //global i-th index
2292
2293 // distribute global residual for (lumped) mass matrix
2294 globalResidual.data()[gi] += elementResidual_u[i];
2295 // distribute entropy_residual
2296 if (STABILIZATION_TYPE==2) // EV Stab
2297 global_entropy_residual[gi] += element_entropy_residual[i];
2298
2299 // distribute transport matrices
2300 for (int j=0;j<nDOF_trial_element;j++)
2301 {
2302 int eN_i_j = eN_i*nDOF_trial_element+j;
2303 TransportMatrix[csrRowIndeces_CellLoops.data()[eN_i] +
2304 csrColumnOffsets_CellLoops.data()[eN_i_j]] += elementTransport[i][j];
2305 TransposeTransportMatrix[csrRowIndeces_CellLoops.data()[eN_i] +
2306 csrColumnOffsets_CellLoops.data()[eN_i_j]]
2307 += elementTransposeTransport[i][j];
2308 }//j
2309 }//i
2310 }//elements
2311
2313 // ADD OUTFLOW BOUNDARY TERM TO TRANSPORT MATRICES AND COMPUTE INFLOW BOUNDARY INTEGRAL //
2315 // * Compute outflow boundary integral as a matrix; i.e., int_B[ (vel.normal)*wi*wj*dx]
2316 for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
2317 {
2318 double min_u_bc_local = 1E10, max_u_bc_local = -1E10;
2319 int ebN = exteriorElementBoundariesArray.data()[ebNE];
2320 int eN = elementBoundaryElementsArray.data()[ebN*2+0],
2321 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
2322 eN_nDOF_trial_element = eN*nDOF_trial_element;
2323 double elementResidual_u[nDOF_test_element];
2324 for (int i=0;i<nDOF_test_element;i++)
2325 elementResidual_u[i]=0.0;
2326 // loop on quad points
2327 for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
2328 {
2329 int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
2330 ebNE_kb_nSpace = ebNE_kb*nSpace,
2331 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
2332 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
2333 double
2334 u_ext=0.0, bc_u_ext=0.0,
2335 porosity_times_velocity[nSpace],
2336 flux_ext=0.0, dflux_ext=0.0,
2337 fluxTransport[nDOF_trial_element],
2338 jac_ext[nSpace*nSpace],
2339 jacDet_ext,
2340 jacInv_ext[nSpace*nSpace],
2341 boundaryJac[nSpace*(nSpace-1)],
2342 metricTensor[(nSpace-1)*(nSpace-1)],
2343 metricTensorDetSqrt,
2344 dS,
2345 u_test_dS[nDOF_test_element],
2346 normal[nSpace],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,porosity_ext;
2347 // calculate mappings
2348 ck.calculateMapping_elementBoundary(eN,
2349 ebN_local,
2350 kb,
2351 ebN_local_kb,
2352 mesh_dof.data(),
2353 mesh_l2g.data(),
2354 mesh_trial_trace_ref.data(),
2355 mesh_grad_trial_trace_ref.data(),
2356 boundaryJac_ref.data(),
2357 jac_ext,
2358 jacDet_ext,
2359 jacInv_ext,
2360 boundaryJac,
2361 metricTensor,
2362 metricTensorDetSqrt,
2363 normal_ref.data(),
2364 normal,
2365 x_ext,y_ext,z_ext);
2366 ck.calculateMappingVelocity_elementBoundary(eN,
2367 ebN_local,
2368 kb,
2369 ebN_local_kb,
2370 mesh_velocity_dof.data(),
2371 mesh_l2g.data(),
2372 mesh_trial_trace_ref.data(),
2373 xt_ext,yt_ext,zt_ext,
2374 normal,
2375 boundaryJac,
2376 metricTensor,
2377 integralScaling);
2378 dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref.data()[kb];
2379 //compute shape and solution information
2380 ck.valFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],u_ext);
2381 //precalculate test function products with integration weights
2382 for (int j=0;j<nDOF_trial_element;j++)
2383 u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
2384
2385 //VRANS
2386 porosity_ext = ebqe_porosity_ext.data()[ebNE_kb];
2387 //
2388 //moving mesh
2389 //
2390 double mesh_velocity[3];
2391 mesh_velocity[0] = xt_ext;
2392 mesh_velocity[1] = yt_ext;
2393 mesh_velocity[2] = zt_ext;
2394 //std::cout<<"mesh_velocity ext"<<std::endl;
2395 for (int I=0;I<nSpace;I++)
2396 porosity_times_velocity[I] = porosity_ext*(ebqe_velocity_ext.data()[ebNE_kb_nSpace+I] - MOVING_DOMAIN*mesh_velocity[I]);
2397 //
2398 //calculate the fluxes
2399 //
2400 double flow = 0.;
2401 for (int I=0; I < nSpace; I++)
2402 flow += normal[I]*porosity_times_velocity[I];
2403
2404 if (flow >= 0) //outflow. This is handled via the transport matrices. Then flux_ext=0 and dflux_ext!=0
2405 {
2406 dflux_ext = flow;
2407 flux_ext = 0;
2408 // save external u
2409 ebqe_u.data()[ebNE_kb] = u_ext;
2410 }
2411 else // inflow. This is handled via the boundary integral. Then flux_ext!=0 and dflux_ext=0
2412 {
2413 dflux_ext = 0;
2414 // save external u
2415 ebqe_u.data()[ebNE_kb] = isDOFBoundary_u.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+(1-isDOFBoundary_u.data()[ebNE_kb])*u_ext;
2416 if (isDOFBoundary_u.data()[ebNE_kb] == 1)
2417 flux_ext = ebqe_bc_u_ext.data()[ebNE_kb]*flow;
2418 else if (isFluxBoundary_u.data()[ebNE_kb] == 1)
2419 flux_ext = ebqe_bc_flux_u_ext.data()[ebNE_kb];
2420 else
2421 {
2422 std::cout<<"warning: VOF open boundary with no external trace, setting to zero for inflow"<<std::endl;
2423 flux_ext = 0.0;
2424 }
2425 }
2426
2427 for (int j=0;j<nDOF_trial_element;j++)
2428 {
2429 // elementResidual. This is to include the inflow boundary integral.
2430 // NOTE: here I assume that we use a Galerkin approach st nDOF_test_element = nDOF_trial_element
2431 elementResidual_u[j] += flux_ext*u_test_dS[j];
2432 int ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j;
2433 fluxTransport[j] = dflux_ext*u_trial_trace_ref.data()[ebN_local_kb_j];
2434 }//j
2436 // DISTRIBUTE OUTFLOW BOUNDARY TO TRANSPORT MATRICES //
2438 for (int i=0;i<nDOF_test_element;i++)
2439 {
2440 int eN_i = eN*nDOF_test_element+i;
2441 for (int j=0;j<nDOF_trial_element;j++)
2442 {
2443 int ebN_i_j = ebN*4*nDOF_test_X_trial_element + i*nDOF_trial_element + j;
2444 TransportMatrix[csrRowIndeces_CellLoops.data()[eN_i] + csrColumnOffsets_eb_CellLoops.data()[ebN_i_j]]
2445 += fluxTransport[j]*u_test_dS[i];
2446 TransposeTransportMatrix[csrRowIndeces_CellLoops.data()[eN_i] + csrColumnOffsets_eb_CellLoops.data()[ebN_i_j]]
2447 += fluxTransport[i]*u_test_dS[j];
2448 }//j
2449 }//i
2450 // local min/max at boundary
2451 min_u_bc_local = fmin(ebqe_u.data()[ebNE_kb], min_u_bc_local);
2452 max_u_bc_local = fmax(ebqe_u.data()[ebNE_kb], max_u_bc_local);
2453 }//kb
2454 // global min/max at boundary
2455 for (int i=0;i<nDOF_test_element;i++)
2456 {
2457 int eN_i = eN*nDOF_test_element+i;
2458 int gi = offset_u+stride_u*u_l2g.data()[eN_i]; //global i-th index
2459 globalResidual.data()[gi] += dt*elementResidual_u[i];
2460 boundary_integral[gi] += elementResidual_u[i];
2461 min_u_bc[gi] = fmin(min_u_bc_local,min_u_bc[gi]);
2462 max_u_bc[gi] = fmax(max_u_bc_local,max_u_bc[gi]);
2463 }
2464 }//ebNE
2465 // END OF ADDING BOUNDARY TERM TO TRANSPORT MATRICES and COMPUTING BOUNDARY INTEGRAL //
2466
2468 // COMPUTE SMOOTHNESS INDICATOR and NORMALIZE ENTROPY RESIDUAL //
2470 // NOTE: see NCLS.h for a different but equivalent implementation of this.
2471 int ij = 0;
2472 for (int i=0; i<numDOFs; i++)
2473 {
2474 double etaMaxi, etaMini;
2475 if (STABILIZATION_TYPE==2) //EV
2476 {
2477 // For eta min and max
2478 etaMaxi = fabs(eta[i]);
2479 etaMini = fabs(eta[i]);
2480 }
2481 double porosity_times_solni = porosity_dof.data()[i]*u_dof_old.data()[i];
2482 // for smoothness indicator //
2483 double alpha_numerator = 0., alpha_denominator = 0.;
2484 for (int offset=csrRowIndeces_DofLoops.data()[i]; offset<csrRowIndeces_DofLoops.data()[i+1]; offset++)
2485 { // First loop in j (sparsity pattern)
2486 int j = csrColumnOffsets_DofLoops.data()[offset];
2487 if (STABILIZATION_TYPE==2) //EV Stabilization
2488 {
2489 // COMPUTE ETA MIN AND ETA MAX //
2490 etaMaxi = fmax(etaMaxi,fabs(eta[j]));
2491 etaMini = fmin(etaMini,fabs(eta[j]));
2492 }
2493 double porosity_times_solnj = porosity_dof.data()[j]*u_dof_old.data()[j];
2494 alpha_numerator += porosity_times_solni - porosity_times_solnj;
2495 alpha_denominator += fabs(porosity_times_solni - porosity_times_solnj);
2496 //update ij
2497 ij+=1;
2498 }
2499 if (STABILIZATION_TYPE==2) //EV Stab
2500 {
2501 // Normalize entropy residual
2502 global_entropy_residual[i] *= etaMini == etaMaxi ? 0. : 2*cE/(etaMaxi-etaMini);
2503 quantDOFs[i] = fabs(global_entropy_residual[i]);
2504 }
2505
2506 double alphai = alpha_numerator/(alpha_denominator+1E-15);
2507 quantDOFs[i] = alphai;
2508
2509 if (POWER_SMOOTHNESS_INDICATOR==0)
2510 psi[i] = 1.0;
2511 else
2512 psi[i] = std::pow(alphai,POWER_SMOOTHNESS_INDICATOR); //NOTE: they use alpha^2 in the paper
2513 }
2515 // ** LOOP IN DOFs FOR EDGE BASED TERMS ** //
2517 ij=0;
2518 for (int i=0; i<numDOFs; i++)
2519 {
2520 // NOTE: Transport matrices already have the porosity considered. ---> Dissipation matrices as well.
2521 double solni = u_dof_old.data()[i]; // solution at time tn for the ith DOF
2522 double porosityi = porosity_dof.data()[i];
2523 double ith_dissipative_term = 0;
2524 double ith_low_order_dissipative_term = 0;
2525 double ith_flux_term = 0;
2526 double dLii = 0.;
2527
2528 // loop over the sparsity pattern of the i-th DOF
2529 for (int offset=csrRowIndeces_DofLoops.data()[i]; offset<csrRowIndeces_DofLoops.data()[i+1]; offset++)
2530 {
2531 int j = csrColumnOffsets_DofLoops.data()[offset];
2532 double solnj = u_dof_old.data()[j]; // solution at time tn for the jth DOF
2533 double porosityj = porosity_dof.data()[j];
2534 double dLowij, dLij, dEVij, dHij;
2535
2536 ith_flux_term += TransportMatrix[ij]*solnj;
2537 if (i != j) //NOTE: there is really no need to check for i!=j (see formula for ith_dissipative_term)
2538 {
2539 // artificial compression
2540 double solij = 0.5*(porosityi*solni+porosityj*solnj);
2541 double Compij = cK*fmax(solij*(1.0-solij),0.0)/(fabs(porosityi*solni-porosityj*solnj)+1E-14);
2542 // first-order dissipative operator
2543 dLowij = fmax(fabs(TransportMatrix[ij]),fabs(TransposeTransportMatrix[ij]));
2544 //dLij = fmax(0.,fmax(psi[i]*TransportMatrix[ij], // Approach by S. Badia
2545 // psi[j]*TransposeTransportMatrix[ij]));
2546 dLij = dLowij*fmax(psi[i],psi[j]); // Approach by JLG & BP
2547 if (STABILIZATION_TYPE==2) //EV Stab
2548 {
2549 // high-order (entropy viscosity) dissipative operator
2550 dEVij = fmax(fabs(global_entropy_residual[i]),fabs(global_entropy_residual[j]));
2551 dHij = fmin(dLowij,dEVij) * fmax(1.0-Compij,0.0); // artificial compression
2552 }
2553 else // smoothness based indicator
2554 {
2555 dHij = dLij * fmax(1.0-Compij,0.0); // artificial compression
2556 }
2557 //dissipative terms
2558 ith_dissipative_term += dHij*(solnj-solni);
2559 ith_low_order_dissipative_term += dLowij*(solnj-solni);
2560 //dHij - dLij. This matrix is needed during FCT step
2561 dt_times_dH_minus_dL[ij] = dt*(dHij - dLowij);
2562 dLii -= dLij;
2563 dLow[ij] = dLowij;
2564 }
2565 else //i==j
2566 {
2567 // NOTE: this is incorrect. Indeed, dLii = -sum_{j!=i}(dLij) and similarly for dCii.
2568 // However, it is irrelevant since during the FCT step we do (dL-dC)*(solnj-solni)
2569 dt_times_dH_minus_dL[ij]=0;
2570 dLow[ij]=0;
2571 }
2572 //update ij
2573 ij+=1;
2574 }
2575 double mi = ML.data()[i];
2576 // compute edge_based_cfl
2577 edge_based_cfl.data()[i] = 2.*fabs(dLii)/mi;
2578 uLow[i] = u_dof_old.data()[i] - dt/mi*(ith_flux_term
2579 + boundary_integral[i]
2580 - ith_low_order_dissipative_term);
2581
2582 // update residual
2583 if (LUMPED_MASS_MATRIX==1)
2584 globalResidual.data()[i] = u_dof_old.data()[i] - dt/mi*(ith_flux_term
2585 + boundary_integral[i]
2586 - ith_dissipative_term);
2587 else
2588 globalResidual.data()[i] += dt*(ith_flux_term - ith_dissipative_term);
2589 }
2590 }
2591 };//VOF
2592
2593 inline VOF_base* newVOF(int nSpaceIn,
2594 int nQuadraturePoints_elementIn,
2595 int nDOF_mesh_trial_elementIn,
2596 int nDOF_trial_elementIn,
2597 int nDOF_test_elementIn,
2598 int nQuadraturePoints_elementBoundaryIn,
2599 int CompKernelFlag)
2600 {
2601 if (nSpaceIn == 1)
2602 return proteus::chooseAndAllocateDiscretization1D<VOF_base,VOF,CompKernel>(nSpaceIn,
2603 nQuadraturePoints_elementIn,
2604 nDOF_mesh_trial_elementIn,
2605 nDOF_trial_elementIn,
2606 nDOF_test_elementIn,
2607 nQuadraturePoints_elementBoundaryIn,
2608 CompKernelFlag);
2609 else if (nSpaceIn == 2)
2610 return proteus::chooseAndAllocateDiscretization2D<VOF_base,VOF,CompKernel>(nSpaceIn,
2611 nQuadraturePoints_elementIn,
2612 nDOF_mesh_trial_elementIn,
2613 nDOF_trial_elementIn,
2614 nDOF_test_elementIn,
2615 nQuadraturePoints_elementBoundaryIn,
2616 CompKernelFlag);
2617 else
2618 return proteus::chooseAndAllocateDiscretization<VOF_base,VOF,CompKernel>(nSpaceIn,
2619 nQuadraturePoints_elementIn,
2620 nDOF_mesh_trial_elementIn,
2621 nDOF_trial_elementIn,
2622 nDOF_test_elementIn,
2623 nQuadraturePoints_elementBoundaryIn,
2624 CompKernelFlag);
2625 }
2626}//proteus
2627#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
equivalent_polynomials::GeneralizedFunctions_mix::useExact
bool useExact
Definition equivalent_polynomials.h:1289
equivalent_polynomials::GeneralizedFunctions_mix::set_quad
void set_quad(unsigned int q)
Definition equivalent_polynomials.h:1325
equivalent_polynomials::GeneralizedFunctions_mix::H
double H(double eps, double phi)
Definition equivalent_polynomials.h:1337
equivalent_polynomials::GeneralizedFunctions_mix::calculate
int calculate(const double *phi_dof, const double *phi_nodes, const double *xi_r, double ma, double mb, double jf, bool isBoundary, bool scale)
Definition equivalent_polynomials.h:1294
equivalent_polynomials::GeneralizedFunctions_mix::set_boundary_quad
void set_boundary_quad(unsigned int ebq)
Definition equivalent_polynomials.h:1331
proteus::VOF_base
Definition VOF.h:59
proteus::VOF_base::psi
std::valarray< double > psi
Definition VOF.h:65
proteus::VOF_base::TransportMatrix
std::valarray< double > TransportMatrix
Definition VOF.h:64
proteus::VOF_base::calculateResidualEdgeBased
virtual void calculateResidualEdgeBased(arguments_dict &args)=0
proteus::VOF_base::maxEntRes
std::valarray< double > maxEntRes
Definition VOF.h:66
proteus::VOF_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)=0
proteus::VOF_base::TransposeTransportMatrix
std::valarray< double > TransposeTransportMatrix
Definition VOF.h:64
proteus::VOF_base::FluxCorrectionMatrix
std::valarray< double > FluxCorrectionMatrix
Definition VOF.h:63
proteus::VOF_base::~VOF_base
virtual ~VOF_base()
Definition VOF.h:67
proteus::VOF_base::maxVel
std::valarray< double > maxVel
Definition VOF.h:66
proteus::VOF_base::global_entropy_residual
std::valarray< double > global_entropy_residual
Definition VOF.h:65
proteus::VOF_base::Rpos
std::valarray< double > Rpos
Definition VOF.h:62
proteus::VOF_base::boundary_integral
std::valarray< double > boundary_integral
Definition VOF.h:65
proteus::VOF_base::Rneg
std::valarray< double > Rneg
Definition VOF.h:62
proteus::VOF_base::eta
std::valarray< double > eta
Definition VOF.h:65
proteus::VOF_base::FCTStep
virtual void FCTStep(arguments_dict &args)=0
proteus::VOF_base::calculateResidualElementBased
virtual void calculateResidualElementBased(arguments_dict &args)=0
proteus::VOF::FCTStep
void FCTStep(arguments_dict &args)
Definition VOF.h:1913
proteus::VOF::exteriorNumericalAdvectiveFlux
void exteriorNumericalAdvectiveFlux(const int &isDOFBoundary_u, const int &isFluxBoundary_u, const double n[nSpace], const double &bc_u, const double &bc_flux_u, const double &u, const double velocity[nSpace], double &flux)
Definition VOF.h:181
proteus::VOF::ck
CompKernelType ck
Definition VOF.h:87
proteus::VOF::calculateCFL
void calculateCFL(const double &elementDiameter, const double df[nSpace], double &cfl)
Definition VOF.h:95
proteus::VOF::calculateJacobian
void calculateJacobian(arguments_dict &args)
Definition VOF.h:1246
proteus::VOF::calculateNumericalDiffusion
void calculateNumericalDiffusion(const double &shockCapturingDiffusion, const double &elementDiameter, const double &strong_residual, const double grad_u[nSpace], double &numDiff)
Definition VOF.h:161
proteus::VOF::calculateSubgridError_tau
void calculateSubgridError_tau(const double &Ct_sge, const double G[nSpace *nSpace], const double &A0, const double Ai[nSpace], double &tau_v, double &cfl)
Definition VOF.h:145
proteus::VOF::calculateResidualEdgeBased
void calculateResidualEdgeBased(arguments_dict &args)
Definition VOF.h:2028
proteus::VOF::calculateSubgridError_tau
void calculateSubgridError_tau(const double &elementDiameter, const double &dmt, const double dH[nSpace], double &cfl, double &tau)
Definition VOF.h:127
proteus::VOF::cutfem_local_boundaries
std::map< int, int > cutfem_local_boundaries
Definition VOF.h:85
proteus::VOF::evaluateCoefficients
void evaluateCoefficients(const double v[nSpace], const double &u, const double &porosity, double &m, double &dm, double f[nSpace], double df[nSpace])
Definition VOF.h:109
proteus::VOF::cutfem_boundaries
std::set< int > cutfem_boundaries
Definition VOF.h:84
proteus::VOF::gf_s
GeneralizedFunctions< nSpace, 2, nQuadraturePoints_element, nQuadraturePoints_elementBoundary > gf_s
Definition VOF.h:88
proteus::VOF::exteriorNumericalAdvectiveFluxDerivative
void exteriorNumericalAdvectiveFluxDerivative(const int &isDOFBoundary_u, const int &isFluxBoundary_u, const double n[nSpace], const double velocity[nSpace], double &dflux)
Definition VOF.h:225
proteus::VOF::nDOF_test_X_trial_element
const int nDOF_test_X_trial_element
Definition VOF.h:86
proteus::VOF::calculateResidualElementBased
void calculateResidualElementBased(arguments_dict &args)
Definition VOF.h:259
equivalent_polynomials.h
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::GeneralizedFunctions
equivalent_polynomials::GeneralizedFunctions_mix< nSpace, nP_ifem, nP, nQ, nEBQ, true > GeneralizedFunctions
Definition ADR.h:21
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::newVOF
VOF_base * newVOF(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition VOF.h:2593
proteus::DENTROPY
double DENTROPY(const double &phi, const double &dummyL, const double &dummyR)
Definition NCLS.h:25
proteus::chooseAndAllocateDiscretization1D
Model_Base * chooseAndAllocateDiscretization1D(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition ModelFactory.h:502
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