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