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