proteus 1.9.0
C/C++/Fortran libraries
Loading...
Searching...
No Matches
NCLS3P.h
Go to the documentation of this file.
1#ifndef NCLS3P_H
2#define NCLS3P_H
3#include <cmath>
4#include <iostream>
5#include "CompKernel.h"
6#include "ModelFactory.h"
7#include "ArgumentsDict.h"
8#include "xtensor-python/pyarray.hpp"
9
10#define cE 0.1
11#define cMax 0.1
12
13namespace proteus
14{
15 inline double ENTROPY(const double& u){
16 return u;
17 }
18 inline double DENTROPY(const double& u){
19 return 1.0;
20 }
21 inline double Sign(const double& z){
22 return (z >= 0.0 ? 1.0 : -1.0);
23 }
24}
25
26namespace proteus
27{
28 class cppNCLS3P_base
29 {
30 //The base class defining the interface
31 public:
32 virtual ~cppNCLS3P_base(){}
33 virtual void calculateResidual(arguments_dict& args)=0;
34 virtual void calculateJacobian(arguments_dict& args)=0;
35 virtual void calculateWaterline(arguments_dict& args)=0;
36 };
37
38 template<class CompKernelType,
39 int nSpace,
40 int nQuadraturePoints_element,
41 int nDOF_mesh_trial_element,
42 int nDOF_trial_element,
43 int nDOF_test_element,
44 int nQuadraturePoints_elementBoundary>
45 class cppNCLS3P : public cppNCLS3P_base
46 {
47 public:
48 const int nDOF_test_X_trial_element;
49 CompKernelType ck;
50 cppNCLS3P():
51 nDOF_test_X_trial_element(nDOF_test_element*nDOF_trial_element),
52 ck()
53 {}
54
55 inline
56 void calculateCFL(const double& elementDiameter,
57 const double df[nSpace],
58 double& cfl)
59 {
60 double h,nrm_v;
61 h = elementDiameter;
62 nrm_v=0.0;
63 for(int I=0;I<nSpace;I++)
64 nrm_v+=df[I]*df[I];
65 nrm_v = sqrt(nrm_v);
66 cfl = nrm_v/h;
67 }
68
69 inline void evaluateCoefficients(const double v[nSpace],
70 const double& u,
71 const double grad_u[nSpace],
72 double& m,
73 double& dm,
74 double& H,
75 double dH[nSpace])
76 {
77 m = u;
78 dm=1.0;
79 H = 0.0;
80 for (int I=0; I < nSpace; I++)
81 {
82 H += v[I]*grad_u[I];
83 dH[I] = v[I];
84 }
85 }
86
87 inline
88 void calculateSubgridError_tau(const double& elementDiameter,
89 const double& dmt,
90 const double dH[nSpace],
91 double& cfl,
92 double& tau)
93 {
94 double h,nrm_v,oneByAbsdt;
95 h = elementDiameter;
96 nrm_v=0.0;
97 for(int I=0;I<nSpace;I++)
98 nrm_v+=dH[I]*dH[I];
99 nrm_v = sqrt(nrm_v);
100 cfl = nrm_v/h;
101 oneByAbsdt = fabs(dmt);
102 tau = 1.0/(2.0*nrm_v/h + oneByAbsdt + 1.0e-8);
103 }
104
105
106 inline
107 void calculateSubgridError_tau( const double& Ct_sge,
108 const double G[nSpace*nSpace],
109 const double& A0,
110 const double Ai[nSpace],
111 double& tau_v,
112 double& cfl)
113 {
114 double v_d_Gv=0.0;
115 for(int I=0;I<nSpace;I++)
116 for (int J=0;J<nSpace;J++)
117 v_d_Gv += Ai[I]*G[I*nSpace+J]*Ai[J];
118
119 tau_v = 1.0/sqrt(Ct_sge*A0*A0 + v_d_Gv + 1.0e-8);
120 }
121
122 void exteriorNumericalFlux(const double n[nSpace],
123 const double& bc_u,
124 const double& u,
125 const double velocity[nSpace],
126 const double velocity_movingDomain[nSpace],
127 double& flux)
128 {
129 double flow_total=0.0,flow_fluid=0.0,flow_movingDomain=0.0;
130 for (int I=0; I < nSpace; I++)
131 {
132 flow_fluid += n[I]*velocity[I];
133 //flow_movingDomain -= n[I]*velocity_movingDomain[I];
134 }
135 flow_total = flow_fluid+flow_movingDomain;
136 if (flow_total > 0.0)
137 {
138 flux = u*flow_movingDomain;
139 }
140 else
141 {
142
143 flux = bc_u*flow_movingDomain - flow_fluid*(u-bc_u);
144 //std::cout<<"bc_u "<<bc_u<<" flow_fluid "<<flow_fluid<<" u "<<u<<std::endl;
145 }
146 }
147
148 inline
149 void exteriorNumericalFluxDerivative(const double n[nSpace],
150 const double velocity[nSpace],
151 const double velocity_movingDomain[nSpace],
152 double& dflux)
153 {
154 double flow_total=0.0,flow_fluid=0.0,flow_movingDomain=0.0;
155 for (int I=0; I < nSpace; I++)
156 {
157 flow_fluid += n[I]*velocity[I];
158 //flow_movingDomain -= n[I]*velocity_movingDomain[I];
159 }
160 flow_total=flow_fluid+flow_movingDomain;
161 if (flow_total > 0.0)
162 {
163 dflux = flow_movingDomain;
164 }
165 else
166 {
167 dflux = -flow_fluid;
168 }
169 }
170
171 void calculateResidual(arguments_dict& args)
172 {
173 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
174 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
175 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
176 xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
177 double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
178 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
179 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
180 xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
181 xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
182 xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
183 xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
184 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
185 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
186 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
187 xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
188 xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
189 xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
190 xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
191 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
192 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
193 int nElements_global = args.scalar<int>("nElements_global");
194 double useMetrics = args.scalar<double>("useMetrics");
195 double alphaBDF = args.scalar<double>("alphaBDF");
196 int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
197 double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
198 double sc_uref = args.scalar<double>("sc_uref");
199 double sc_alpha = args.scalar<double>("sc_alpha");
200 xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
201 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
202 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
203 xt::pyarray<double>& u_dof_old = args.array<double>("u_dof_old");
204 xt::pyarray<double>& velocity = args.array<double>("velocity");
205 xt::pyarray<double>& q_m = args.array<double>("q_m");
206 xt::pyarray<double>& q_u = args.array<double>("q_u");
207 xt::pyarray<double>& q_n = args.array<double>("q_n");
208 xt::pyarray<double>& q_dH = args.array<double>("q_dH");
209 xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
210 xt::pyarray<double>& q_dV = args.array<double>("q_dV");
211 xt::pyarray<double>& q_dV_last = args.array<double>("q_dV_last");
212 xt::pyarray<double>& cfl = args.array<double>("cfl");
213 xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
214 xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
215 int offset_u = args.scalar<int>("offset_u");
216 int stride_u = args.scalar<int>("stride_u");
217 xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
218 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
219 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
220 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
221 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
222 xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
223 xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
224 xt::pyarray<double>& ebqe_rd_u_ext = args.array<double>("ebqe_rd_u_ext");
225 xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
226 xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
227 xt::pyarray<double>& cell_interface_locator = args.array<double>("cell_interface_locator");
228 xt::pyarray<double>& interface_locator = args.array<double>("interface_locator");
229 int EXPLICIT_METHOD = args.scalar<int>("EXPLICIT_METHOD");
230 double degree_polynomial = args.scalar<double>("degree_polynomial");
231 int stage = args.scalar<int>("stage");
232 xt::pyarray<double>& uTilde_dof = args.array<double>("uTilde_dof");
233 double dt = args.scalar<double>("dt");
234 int PURE_BDF = args.scalar<int>("PURE_BDF");
235 double meanEntropy = 0., meanOmega = 0., maxEntropy = -1E10, minEntropy = 1E10;
236 double maxVel[nElements_global], maxEntRes[nElements_global];
237 //cek should this be read in?
238 double Ct_sge = 4.0;
239
240 //loop over elements to compute volume integrals and load them into element and global residual
241 //
242 //eN is the element index
243 //eN_k is the quadrature point index for a scalar
244 //eN_k_nSpace is the quadrature point index for a vector
245 //eN_i is the element test function index
246 //eN_j is the element trial function index
247 //eN_k_j is the quadrature point index for a trial function
248 //eN_k_i is the quadrature point index for a trial function
249 for(int eN=0;eN<nElements_global;eN++)
250 {
251 // init maxVel and maxEntRes
252 maxVel[eN] = 0.;
253 maxEntRes[eN] = 0.;
254 //declare local storage for element residual and initialize
255 double elementResidual_u[nDOF_test_element];
256 for (int i=0;i<nDOF_test_element;i++)
257 {
258 elementResidual_u[i]=0.0;
259 }//i
260 //loop over quadrature points and compute integrands
261 for (int k=0;k<nQuadraturePoints_element;k++)
262 {
263 //compute indeces and declare local storage
264 int eN_k = eN*nQuadraturePoints_element+k,
265 eN_k_nSpace = eN_k*nSpace,
266 eN_nDOF_trial_element = eN*nDOF_trial_element;
267 double u=0.0,grad_u[nSpace],grad_u_old[nSpace],
268 un=0.0, Hn=0.0, HTilde=0.0, grad_uTilde[nSpace],
269 m=0.0,dm=0.0,
270 H=0.0,dH[nSpace],
271 f[nSpace],df[nSpace],//for MOVING_DOMAIN
272 m_t=0.0,dm_t=0.0,
273 pdeResidual_u=0.0,
274 Lstar_u[nDOF_test_element],
275 subgridError_u=0.0,
276 tau=0.0,tau0=0.0,tau1=0.0,
277 numDiff0=0.0,numDiff1=0.0,
278 jac[nSpace*nSpace],
279 jacDet,
280 jacInv[nSpace*nSpace],
281 u_grad_trial[nDOF_trial_element*nSpace],
282 u_test_dV[nDOF_trial_element],
283 u_grad_test_dV[nDOF_test_element*nSpace],
284 dV,x,y,z,xt,yt,zt,
285 G[nSpace*nSpace],G_dd_G,tr_G;//,norm_Rv;
286 //
287 //compute solution and gradients at quadrature points
288 //
289 ck.calculateMapping_element(eN,
290 k,
291 mesh_dof.data(),
292 mesh_l2g.data(),
293 mesh_trial_ref.data(),
294 mesh_grad_trial_ref.data(),
295 jac,
296 jacDet,
297 jacInv,
298 x,y,z);
299 ck.calculateMappingVelocity_element(eN,
300 k,
301 mesh_velocity_dof.data(),
302 mesh_l2g.data(),
303 mesh_trial_ref.data(),
304 xt,yt,zt);
305 //get the physical integration weight
306 dV = fabs(jacDet)*dV_ref[k];
307 ck.calculateG(jacInv,G,G_dd_G,tr_G);
308 //get the trial function gradients
309 ck.gradTrialFromRef(&u_grad_trial_ref[k*nDOF_trial_element*nSpace],
310 jacInv,
311 u_grad_trial);
312 //get the solution
313 ck.valFromDOF(u_dof.data(),
314 &u_l2g[eN_nDOF_trial_element],
315 &u_trial_ref[k*nDOF_trial_element],
316 u);
317 ck.valFromDOF(u_dof_old.data(),
318 &u_l2g[eN_nDOF_trial_element],
319 &u_trial_ref[k*nDOF_trial_element],
320 un);
321 //get the solution gradients
322 ck.gradFromDOF(u_dof.data(),
323 &u_l2g[eN_nDOF_trial_element],
324 u_grad_trial,
325 grad_u);
326 ck.gradFromDOF(u_dof_old.data(),
327 &u_l2g[eN_nDOF_trial_element],
328 u_grad_trial,
329 grad_u_old);
330 ck.gradFromDOF(uTilde_dof.data(),
331 &u_l2g[eN_nDOF_trial_element],
332 u_grad_trial,
333 grad_uTilde);
334 //precalculate test function products with integration weights
335 for (int j=0;j<nDOF_trial_element;j++)
336 {
337 u_test_dV[j] = u_test_ref[k*nDOF_trial_element+j]*dV;
338 for (int I=0;I<nSpace;I++)
339 {
340 u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
341 }
342 }
343 //save solution at quadrature points for other models to use
344 q_u[eN_k]=u;
345 for (int I=0;I<nSpace;I++)
346 q_n[eN_k_nSpace+I] = grad_u[I];
347 //
348 //calculate pde coefficients at quadrature points
349 //
350 evaluateCoefficients(&velocity[eN_k_nSpace],
351 u,
352 grad_u,
353 m,
354 dm,
355 H,
356 dH);
357 //
358 //moving mesh
359 //
360 double mesh_velocity[3];
361 mesh_velocity[0] = xt;
362 mesh_velocity[1] = yt;
363 mesh_velocity[2] = zt;
364 for (int I=0;I<nSpace;I++)
365 {
366 f[I] = -MOVING_DOMAIN*m*mesh_velocity[I];
367 df[I] = -MOVING_DOMAIN*dm*mesh_velocity[I];
368 }
369 //
370 //calculate time derivative at quadrature points
371 //
372 if (q_dV_last[eN_k] <= -100)
373 q_dV_last[eN_k] = dV;
374 q_dV[eN_k] = dV;
375 ck.bdf(alphaBDF,
376 q_m_betaBDF[eN_k],
377 m,
378 dm,
379 m_t,
380 dm_t);
381 if (EXPLICIT_METHOD==1)
382 {
383 double normVel=0.;
384 double relVelocity[nSpace];
385 for (int I=0;I<nSpace;I++)
386 {
387 Hn += velocity[eN_k_nSpace+I]*grad_u_old[I];
388 HTilde += velocity[eN_k_nSpace+I]*grad_uTilde[I];
389 H += velocity[eN_k_nSpace+I]*grad_u[I];
390 relVelocity[I] = dH[I] - MOVING_DOMAIN*df[I];
391 normVel += relVelocity[I]*relVelocity[I];
392 }
393 normVel = std::sqrt(normVel);
394
395 // calculate CFL
396 calculateCFL(elementDiameter[eN]/degree_polynomial,relVelocity,cfl[eN_k]);
397
398 // compute max velocity at cell
399 maxVel[eN] = fmax(normVel,maxVel[eN]);
400
401 // entropy residual
402 double entRes = (ENTROPY(u)-ENTROPY(un))/dt + 0.5*(DENTROPY(u)*H +
403 DENTROPY(un)*Hn);
404 maxEntRes[eN] = fmax(maxEntRes[eN],fabs(entRes));
405
406 // Quantities for normalization factor //
407 meanEntropy += ENTROPY(u)*dV;
408 meanOmega += dV;
409 maxEntropy = fmax(maxEntropy,ENTROPY(u));
410 minEntropy = fmin(minEntropy,ENTROPY(u));
411 }
412 else
413 {
414 //
415 //calculate subgrid error (strong residual and adjoint)
416 //
417 //calculate strong residual
418 pdeResidual_u = ck.Mass_strong(m_t) +
419 ck.Hamiltonian_strong(dH,grad_u)+
420 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
421
422 //calculate adjoint
423 for (int i=0;i<nDOF_test_element;i++)
424 {
425 //int eN_k_i_nSpace = (eN_k*nDOF_trial_element+i)*nSpace;
426 //Lstar_u[i] = ck.Hamiltonian_adjoint(dH,&u_grad_test_dV[eN_k_i_nSpace]);
427 int i_nSpace = i*nSpace;
428 Lstar_u[i] = ck.Hamiltonian_adjoint(dH,&u_grad_test_dV[i_nSpace])
429 + MOVING_DOMAIN*ck.Advection_adjoint(df,&u_grad_test_dV[i_nSpace]);
430 }
431 //calculate tau and tau*Res
432 double subgridErrorVelocity[nSpace];
433 for (int I=0;I<nSpace;I++)
434 subgridErrorVelocity[I] = dH[I] - MOVING_DOMAIN*df[I];
435
436 calculateSubgridError_tau(elementDiameter[eN],
437 dm_t,
438 subgridErrorVelocity,//dH,
439 cfl[eN_k],
440 tau0);
441
442 calculateSubgridError_tau(Ct_sge,
443 G,
444 dm_t,
445 subgridErrorVelocity,//dH,
446 tau1,
447 cfl[eN_k]);
448
449 tau = useMetrics*tau1+(1.0-useMetrics)*tau0;
450
451 subgridError_u = -tau*pdeResidual_u;
452 //
453 //calculate shock capturing diffusion
454 //
455 ck.calculateNumericalDiffusion(shockCapturingDiffusion,
456 elementDiameter[eN],
457 pdeResidual_u,
458 grad_u,
459 numDiff0);
460 //ck.calculateNumericalDiffusion(shockCapturingDiffusion,G,pdeResidual_u,grad_u_old,numDiff1);
461 ck.calculateNumericalDiffusion(shockCapturingDiffusion,
462 sc_uref,
463 sc_alpha,
464 G,
465 G_dd_G,
466 pdeResidual_u,
467 grad_u,
468 numDiff1);
469 q_numDiff_u[eN_k] = useMetrics*numDiff1+(1.0-useMetrics)*numDiff0;
470 //std::cout<<tau<<" "<<q_numDiff_u[eN_k]<<std::endl;
471 }
472 //
473 //update element residual
474 //
475 double sign;
476 int same_sign=1;
477 for(int i=0;i<nDOF_test_element;i++)
478 {
479 int gi = offset_u+stride_u*u_l2g[eN*nDOF_test_element+i];
480 if (i==0)
481 sign = Sign(u_dof[gi]);
482 else if (same_sign==1)
483 {
484 same_sign = sign == Sign(u_dof[gi]) ? 1 : 0;
485 sign = Sign(u_dof[gi]);
486 }
487 //int eN_k_i=eN_k*nDOF_test_element+i,
488 // eN_k_i_nSpace = eN_k_i*nSpace;
489 int i_nSpace=i*nSpace;
490 if (EXPLICIT_METHOD==1)
491 {
492 if (stage == 1)
493 elementResidual_u[i] +=
494 //ck.Mass_weak(u-un,u_test_dV[i]) + // time derivative
495 ck.Mass_weak(dt*m_t,u_test_dV[i]) + // time derivative
496 1./3*dt*ck.Hamiltonian_weak(Hn,u_test_dV[i]) + // v*grad(phi)
497 1./9*dt*dt*ck.NumericalDiffusion(Hn,dH,&u_grad_test_dV[i_nSpace]) +
498 1./3*dt*ck.NumericalDiffusion(q_numDiff_u_last[eN_k],
499 grad_u_old,
500 &u_grad_test_dV[i_nSpace]);
501 // TODO: Add part about moving mesh
502 else
503 elementResidual_u[i] +=
504 //ck.Mass_weak(u-un,u_test_dV[i]) + // time derivative
505 ck.Mass_weak(dt*m_t,u_test_dV[i]) + // time derivative
506 dt*ck.Hamiltonian_weak(Hn,u_test_dV[i]) + // v*grad(phi)
507 0.5*dt*dt*ck.NumericalDiffusion(HTilde,dH,&u_grad_test_dV[i_nSpace]) +
508 dt*ck.NumericalDiffusion(q_numDiff_u_last[eN_k],
509 grad_u_old,
510 &u_grad_test_dV[i_nSpace]);
511 //elementResidual_u[i] += // semi-implicit Lax Wendroff
512 // ck.Mass_weak(u-un,u_test_dV[i]) +
513 // dt*ck.Hamiltonian_weak(Hn,u_test_dV[i]) +
514 // 0.5*dt*dt*ck.NumericalDiffusion(H,dH,&u_grad_test_dV[i_nSpace]);
515 }
516 else // Implicit SUPG with SHOCK CAPTURING
517 {
518 elementResidual_u[i] +=
519 ck.Mass_weak(m_t,u_test_dV[i]) +
520 ck.Hamiltonian_weak(H,u_test_dV[i]) +
521 MOVING_DOMAIN*ck.Advection_weak(f,&u_grad_test_dV[i_nSpace])+
522 (PURE_BDF == 1 ? 0. : 1.)*
523 (ck.SubgridError(subgridError_u,Lstar_u[i]) +
524 ck.NumericalDiffusion(q_numDiff_u_last[eN_k],
525 grad_u,
526 &u_grad_test_dV[i_nSpace]));
527 }
528 }//i
529 if (same_sign == 0) // This cell contains the interface
530 {
531 cell_interface_locator[eN] = 1.0;
532 for(int i=0;i<nDOF_test_element;i++)
533 {
534 int gi = offset_u+stride_u*u_l2g[eN*nDOF_test_element+i];
535 interface_locator[gi] = 1.0;
536 }
537 }
538 //
539 //cek/ido todo, get rid of m, since u=m
540 //save momentum for time history and velocity for subgrid error
541 //save solution for other models
542 //
543 q_m[eN_k] = m;
544 q_u[eN_k] = u;
545 for (int I=0;I<nSpace;I++)
546 {
547 int eN_k_I = eN_k*nSpace+I;
548 q_dH[eN_k_I] = dH[I];
549 }
550 }
551 //
552 //load element into global residual and save element residual
553 //
554 for(int i=0;i<nDOF_test_element;i++)
555 {
556 int eN_i=eN*nDOF_test_element+i;
557 globalResidual[offset_u+stride_u*u_l2g[eN_i]] += elementResidual_u[i];
558 }//i
559 }//elements
560 if (EXPLICIT_METHOD==1)
561 {
562 meanEntropy /= meanOmega;
563 double norm_factor = fmax(fabs(maxEntropy - meanEntropy), fabs(meanEntropy-minEntropy));
564 for(int eN=0;eN<nElements_global;eN++)
565 {
566 double hK=elementDiameter[eN]/degree_polynomial;
567 double linear_viscosity = cMax*hK*maxVel[eN];
568 double entropy_viscosity = cE*hK*hK*maxEntRes[eN]/norm_factor;
569 for (int k=0;k<nQuadraturePoints_element;k++)
570 {
571 int eN_k = eN*nQuadraturePoints_element+k;
572 q_numDiff_u[eN_k] = fmin(linear_viscosity,entropy_viscosity);
573 }
574 }
575 }
576 //
577 //loop over exterior element boundaries to calculate surface integrals and load into element and global residuals
578 //
579 //ebNE is the Exterior element boundary INdex
580 //ebN is the element boundary INdex
581 //eN is the element index
582 //std::cout <<nExteriorElementBoundaries_global<<std::endl;
583 for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
584 {
585 int ebN = exteriorElementBoundariesArray[ebNE],
586 eN = elementBoundaryElementsArray[ebN*2+0],
587 ebN_local = elementBoundaryLocalElementBoundariesArray[ebN*2+0],
588 eN_nDOF_trial_element = eN*nDOF_trial_element;
589 double elementResidual_u[nDOF_test_element];
590 for (int i=0;i<nDOF_test_element;i++)
591 {
592 elementResidual_u[i]=0.0;
593 }
594 for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
595 {
596 int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
597 ebNE_kb_nSpace = ebNE_kb*nSpace,
598 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
599 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
600 double u_ext=0.0,
601 grad_u_ext[nSpace],
602 m_ext=0.0,
603 dm_ext=0.0,
604 H_ext=0.0,
605 dH_ext[nSpace],
606 //f_ext[nSpace],//MOVING_DOMAIN
607 //df_ext[nSpace],//MOVING_DOMAIN
608 //flux_ext=0.0,
609 bc_u_ext=0.0,
610 bc_grad_u_ext[nSpace],
611 bc_m_ext=0.0,
612 bc_dm_ext=0.0,
613 bc_H_ext=0.0,
614 bc_dH_ext[nSpace],
615 jac_ext[nSpace*nSpace],
616 jacDet_ext,
617 jacInv_ext[nSpace*nSpace],
618 boundaryJac[nSpace*(nSpace-1)],
619 metricTensor[(nSpace-1)*(nSpace-1)],
620 metricTensorDetSqrt,
621 dS,
622 u_test_dS[nDOF_test_element],
623 u_grad_trial_trace[nDOF_trial_element*nSpace],
624 normal[nSpace],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
625 G[nSpace*nSpace],G_dd_G,tr_G,flux_ext;
626 //
627 //calculate the solution and gradients at quadrature points
628 //
629 //compute information about mapping from reference element to physical element
630 ck.calculateMapping_elementBoundary(eN,
631 ebN_local,
632 kb,
633 ebN_local_kb,
634 mesh_dof.data(),
635 mesh_l2g.data(),
636 mesh_trial_trace_ref.data(),
637 mesh_grad_trial_trace_ref.data(),
638 boundaryJac_ref.data(),
639 jac_ext,
640 jacDet_ext,
641 jacInv_ext,
642 boundaryJac,
643 metricTensor,
644 metricTensorDetSqrt,
645 normal_ref.data(),
646 normal,
647 x_ext,y_ext,z_ext);
648 ck.calculateMappingVelocity_elementBoundary(eN,
649 ebN_local,
650 kb,
651 ebN_local_kb,
652 mesh_velocity_dof.data(),
653 mesh_l2g.data(),
654 mesh_trial_trace_ref.data(),
655 xt_ext,yt_ext,zt_ext,
656 normal,
657 boundaryJac,
658 metricTensor,
659 integralScaling);
660 dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref[kb];
661 //get the metric tensor
662 //cek todo use symmetry
663 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
664 //compute shape and solution information
665 //shape
666 ck.gradTrialFromRef(&u_grad_trial_trace_ref[ebN_local_kb_nSpace*nDOF_trial_element],
667 jacInv_ext,
668 u_grad_trial_trace);
669 //solution and gradients
670 if (EXPLICIT_METHOD==1) // explicit
671 {
672 ck.valFromDOF(u_dof_old.data(),
673 &u_l2g[eN_nDOF_trial_element],
674 &u_trial_trace_ref[ebN_local_kb*nDOF_test_element],
675 u_ext);
676 ck.gradFromDOF(u_dof_old.data(),
677 &u_l2g[eN_nDOF_trial_element],
678 u_grad_trial_trace,
679 grad_u_ext);
680 }
681 else
682 {
683 ck.valFromDOF(u_dof.data(),
684 &u_l2g[eN_nDOF_trial_element],
685 &u_trial_trace_ref[ebN_local_kb*nDOF_test_element],
686 u_ext);
687 ck.gradFromDOF(u_dof.data(),
688 &u_l2g[eN_nDOF_trial_element],
689 u_grad_trial_trace,
690 grad_u_ext);
691 }
692 //precalculate test function products with integration weights
693 for (int j=0;j<nDOF_trial_element;j++)
694 {
695 u_test_dS[j] = u_test_trace_ref[ebN_local_kb*nDOF_test_element+j]*dS;
696 }
697 //
698 //load the boundary values
699 //
700 bc_u_ext = (isDOFBoundary_u[ebNE_kb]*ebqe_bc_u_ext[ebNE_kb]
701 +(1-isDOFBoundary_u[ebNE_kb])*ebqe_rd_u_ext[ebNE_kb]);
702 //
703 //calculate the pde coefficients using the solution and the boundary values for the solution
704 //
705 evaluateCoefficients(&ebqe_velocity_ext[ebNE_kb_nSpace],
706 u_ext,
707 grad_u_ext,
708 m_ext,
709 dm_ext,
710 H_ext,
711 dH_ext);
712 evaluateCoefficients(&ebqe_velocity_ext[ebNE_kb_nSpace],
713 bc_u_ext,
714 bc_grad_u_ext,
715 bc_m_ext,
716 bc_dm_ext,
717 bc_H_ext,
718 bc_dH_ext);
719 //
720 //moving mesh
721 //
722 double velocity_ext[nSpace];
723 double mesh_velocity[3];
724 mesh_velocity[0] = xt_ext;
725 mesh_velocity[1] = yt_ext;
726 mesh_velocity[2] = zt_ext;
727 for (int I=0;I<nSpace;I++)
728 velocity_ext[I] = - MOVING_DOMAIN*mesh_velocity[I];
729 //
730 //calculate the numerical fluxes
731 //
732 exteriorNumericalFlux(normal,
733 bc_u_ext,
734 u_ext,
735 dH_ext,
736 velocity_ext,
737 flux_ext);
738 ebqe_u[ebNE_kb] = u_ext;
739
740 if (EXPLICIT_METHOD==1)
741 if (stage==1)
742 flux_ext *= 1./3*dt;
743 else
744 flux_ext *= dt;
745 //std::cout<<u_ext<<ebqe_bc_u_ext
746 //
747 //update residuals
748 //
749 for (int i=0;i<nDOF_test_element;i++)
750 {
751 //int ebNE_kb_i = ebNE_kb*nDOF_test_element+i;
752 elementResidual_u[i] += ck.ExteriorElementBoundaryFlux(flux_ext,u_test_dS[i]);
753 }//i
754 }//kb
755 //
756 //update the element and global residual storage
757 //
758 for (int i=0;i<nDOF_test_element;i++)
759 {
760 int eN_i = eN*nDOF_test_element+i;
761 //globalResidual[offset_u+stride_u*u_l2g[eN_i]] += MOVING_DOMAIN*elementResidual_u[i];
762 globalResidual[offset_u+stride_u*u_l2g[eN_i]] += elementResidual_u[i];
763
764 }//i
765 }//ebNE
766 }
767
768 void calculateJacobian(arguments_dict& args)
769 {
770 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
771 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
772 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
773 xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
774 double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
775 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
776 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
777 xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
778 xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
779 xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
780 xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
781 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
782 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
783 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
784 xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
785 xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
786 xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
787 xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
788 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
789 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
790 int nElements_global = args.scalar<int>("nElements_global");
791 double useMetrics = args.scalar<double>("useMetrics");
792 double alphaBDF = args.scalar<double>("alphaBDF");
793 int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
794 double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
795 xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
796 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
797 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
798 xt::pyarray<double>& velocity = args.array<double>("velocity");
799 xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
800 xt::pyarray<double>& cfl = args.array<double>("cfl");
801 xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
802 xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
803 xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
804 xt::pyarray<double>& globalJacobian = args.array<double>("globalJacobian");
805 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
806 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
807 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
808 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
809 xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
810 xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
811 xt::pyarray<double>& ebqe_rd_u_ext = args.array<double>("ebqe_rd_u_ext");
812 xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
813 xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
814 int EXPLICIT_METHOD = args.scalar<int>("EXPLICIT_METHOD");
815 int PURE_BDF = args.scalar<int>("PURE_BDF");
816 double Ct_sge = 4.0;
817
818 //
819 //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
820 //
821 for(int eN=0;eN<nElements_global;eN++)
822 {
823 double elementJacobian_u_u[nDOF_test_element][nDOF_trial_element];
824 for (int i=0;i<nDOF_test_element;i++)
825 for (int j=0;j<nDOF_trial_element;j++)
826 {
827 elementJacobian_u_u[i][j]=0.0;
828 }
829 for (int k=0;k<nQuadraturePoints_element;k++)
830 {
831 int eN_k = eN*nQuadraturePoints_element+k, //index to a scalar at a quadrature point
832 eN_k_nSpace = eN_k*nSpace,
833 eN_nDOF_trial_element = eN*nDOF_trial_element; //index to a vector at a quadrature point
834
835 //declare local storage
836 double u=0.0,
837 grad_u[nSpace],
838 m=0.0,dm=0.0,
839 H=0.0,dH[nSpace],
840 f[nSpace],df[nSpace],//MOVING_MESH
841 m_t=0.0,dm_t=0.0,
842 dpdeResidual_u_u[nDOF_trial_element],
843 Lstar_u[nDOF_test_element],
844 dsubgridError_u_u[nDOF_trial_element],
845 tau=0.0,tau0=0.0,tau1=0.0,
846 jac[nSpace*nSpace],
847 jacDet,
848 jacInv[nSpace*nSpace],
849 u_grad_trial[nDOF_trial_element*nSpace],
850 dV,
851 u_test_dV[nDOF_test_element],
852 u_grad_test_dV[nDOF_test_element*nSpace],
853 x,y,z,xt,yt,zt,
854 G[nSpace*nSpace],G_dd_G,tr_G;
855 //
856 //calculate solution and gradients at quadrature points
857 //
858 //get jacobian, etc for mapping reference element
859 ck.calculateMapping_element(eN,
860 k,
861 mesh_dof.data(),
862 mesh_l2g.data(),
863 mesh_trial_ref.data(),
864 mesh_grad_trial_ref.data(),
865 jac,
866 jacDet,
867 jacInv,
868 x,y,z);
869 ck.calculateMappingVelocity_element(eN,
870 k,
871 mesh_velocity_dof.data(),
872 mesh_l2g.data(),
873 mesh_trial_ref.data(),
874 xt,yt,zt);
875 //get the physical integration weight
876 dV = fabs(jacDet)*dV_ref[k];
877 ck.calculateG(jacInv,G,G_dd_G,tr_G);
878 //get the trial function gradients
879 ck.gradTrialFromRef(&u_grad_trial_ref[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
880 //get the solution
881 ck.valFromDOF(u_dof.data(),&u_l2g[eN_nDOF_trial_element],&u_trial_ref[k*nDOF_trial_element],u);
882 //get the solution gradients
883 ck.gradFromDOF(u_dof.data(),&u_l2g[eN_nDOF_trial_element],u_grad_trial,grad_u);
884 //precalculate test function products with integration weights
885 for (int j=0;j<nDOF_trial_element;j++)
886 {
887 u_test_dV[j] = u_test_ref[k*nDOF_trial_element+j]*dV;
888 for (int I=0;I<nSpace;I++)
889 {
890 u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
891 }
892 }
893 //
894 //calculate pde coefficients and derivatives at quadrature points
895 //
896 evaluateCoefficients(&velocity[eN_k_nSpace],
897 u,
898 grad_u,
899 m,
900 dm,
901 H,
902 dH);
903 //
904 //moving mesh
905 //
906 double mesh_velocity[3];
907 mesh_velocity[0] = xt;
908 mesh_velocity[1] = yt;
909 mesh_velocity[2] = zt;
910 for (int I=0;I<nSpace;I++)
911 {
912 f[I] = -MOVING_DOMAIN*m*mesh_velocity[I];
913 df[I] = -MOVING_DOMAIN*dm*mesh_velocity[I];
914 }
915 //
916 //calculate time derivatives
917 //
918 ck.bdf(alphaBDF,
919 q_m_betaBDF[eN_k],
920 m,
921 dm,
922 m_t,
923 dm_t);
924 //
925 //calculate subgrid error contribution to the Jacobian (strong residual, adjoint, jacobian of strong residual)
926 //
927 //calculate the adjoint times the test functions
928 for (int i=0;i<nDOF_test_element;i++)
929 {
930 //int eN_k_i_nSpace = (eN_k*nDOF_trial_element+i)*nSpace;
931 //Lstar_u[i]=ck.Hamiltonian_adjoint(dH,&u_grad_test_dV[eN_k_i_nSpace]);
932 int i_nSpace = i*nSpace;
933 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]);
934
935 }
936 //calculate the Jacobian of strong residual
937 for (int j=0;j<nDOF_trial_element;j++)
938 {
939 //int eN_k_j=eN_k*nDOF_trial_element+j;
940 //int eN_k_j_nSpace = eN_k_j*nSpace;
941 int j_nSpace = j*nSpace;
942 dpdeResidual_u_u[j]=ck.MassJacobian_strong(dm_t,u_trial_ref[k*nDOF_trial_element+j]) +
943 ck.HamiltonianJacobian_strong(dH,&u_grad_trial[j_nSpace]) +
944 MOVING_DOMAIN*ck.AdvectionJacobian_strong(df,&u_grad_trial[j_nSpace]);
945
946 }
947 //tau and tau*Res
948 double subgridErrorVelocity[nSpace];
949 for (int I=0;I<nSpace;I++)
950 subgridErrorVelocity[I] = dH[I] - MOVING_DOMAIN*df[I];
951
952 calculateSubgridError_tau(elementDiameter[eN],
953 dm_t,
954 subgridErrorVelocity,//dH,
955 cfl[eN_k],
956 tau0);
957
958 calculateSubgridError_tau(Ct_sge,
959 G,
960 dm_t,
961 subgridErrorVelocity,//dH,
962 tau1,
963 cfl[eN_k]);
964 tau = useMetrics*tau1+(1.0-useMetrics)*tau0;
965
966 for(int j=0;j<nDOF_trial_element;j++)
967 dsubgridError_u_u[j] = -tau*dpdeResidual_u_u[j];
968 for(int i=0;i<nDOF_test_element;i++)
969 {
970 //int eN_k_i=eN_k*nDOF_test_element+i;
971 //int eN_k_i_nSpace=eN_k_i*nSpace;
972 for(int j=0;j<nDOF_trial_element;j++)
973 {
974 //int eN_k_j=eN_k*nDOF_trial_element+j;
975 //int eN_k_j_nSpace = eN_k_j*nSpace;
976 int j_nSpace = j*nSpace;
977 int i_nSpace = i*nSpace;
978 if (EXPLICIT_METHOD==1)
979 {
980 elementJacobian_u_u[i][j] +=
981 ck.MassJacobian_weak(1.0,
982 u_trial_ref[k*nDOF_trial_element+j],
983 u_test_dV[i]);
984 //ck.MassJacobian_weak(1.0, // semi-implicit Lax Wendroff
985 // u_trial_ref[k*nDOF_trial_element+j],
986 // u_test_dV[i]) +
987 //0.5*dt*dt*dV*
988 //ck.NumericalDiffusion(1.0,dH,&u_grad_trial[i_nSpace])*
989 //ck.NumericalDiffusion(1.0,dH,&u_grad_trial[j_nSpace]);
990 }
991 else // Implicit SUPG with SHOCK CAPTURING
992 elementJacobian_u_u[i][j] +=
993 ck.MassJacobian_weak(dm_t,
994 u_trial_ref[k*nDOF_trial_element+j],
995 u_test_dV[i]) +
996 ck.HamiltonianJacobian_weak(dH,&u_grad_trial[j_nSpace],u_test_dV[i]) +
997 MOVING_DOMAIN*ck.AdvectionJacobian_weak(df,
998 u_trial_ref[k*nDOF_trial_element+j],
999 &u_grad_test_dV[i_nSpace]) +
1000 (PURE_BDF == 1 ? 0. : 1.)*
1001 (ck.SubgridErrorJacobian(dsubgridError_u_u[j],Lstar_u[i]) +
1002 ck.NumericalDiffusionJacobian(q_numDiff_u_last[eN_k],
1003 &u_grad_trial[j_nSpace],
1004 &u_grad_test_dV[i_nSpace]));
1005 }//j
1006 }//i
1007 }//k
1008 //
1009 //load into element Jacobian into global Jacobian
1010 //
1011 for (int i=0;i<nDOF_test_element;i++)
1012 {
1013 int eN_i = eN*nDOF_test_element+i;
1014 for (int j=0;j<nDOF_trial_element;j++)
1015 {
1016 int eN_i_j = eN_i*nDOF_trial_element+j;
1017 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_u_u[eN_i_j]] += elementJacobian_u_u[i][j];
1018 }//j
1019 }//i
1020 }//elements
1021 //
1022 //loop over exterior element boundaries to compute the surface integrals and load them into the global Jacobian
1023 //
1024 if (EXPLICIT_METHOD==0) //(semi)-implicit
1025 for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
1026 {
1027 int ebN = exteriorElementBoundariesArray[ebNE];
1028 int eN = elementBoundaryElementsArray[ebN*2+0],
1029 ebN_local = elementBoundaryLocalElementBoundariesArray[ebN*2+0],
1030 eN_nDOF_trial_element = eN*nDOF_trial_element;
1031 for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
1032 {
1033 int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
1034 ebNE_kb_nSpace = ebNE_kb*nSpace,
1035 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
1036 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
1037
1038 double u_ext=0.0,
1039 grad_u_ext[nSpace],
1040 m_ext=0.0,
1041 dm_ext=0.0,
1042 H_ext=0.0,
1043 dH_ext[nSpace],
1044 bc_grad_u_ext[nSpace],
1045 bc_H_ext=0.0,
1046 bc_dH_ext[nSpace],
1047 //f_ext[nSpace],
1048 //df_ext[nSpace],
1049 dflux_u_u_ext=0.0,
1050 bc_u_ext=0.0,
1051 //bc_grad_u_ext[nSpace],
1052 bc_m_ext=0.0,
1053 bc_dm_ext=0.0,
1054 //bc_f_ext[nSpace],
1055 //bc_df_ext[nSpace],
1056 fluxJacobian_u_u[nDOF_trial_element],
1057 jac_ext[nSpace*nSpace],
1058 jacDet_ext,
1059 jacInv_ext[nSpace*nSpace],
1060 boundaryJac[nSpace*(nSpace-1)],
1061 metricTensor[(nSpace-1)*(nSpace-1)],
1062 metricTensorDetSqrt,
1063 dS,
1064 u_test_dS[nDOF_test_element],
1065 u_grad_trial_trace[nDOF_trial_element*nSpace],
1066 normal[nSpace],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
1067 G[nSpace*nSpace],G_dd_G,tr_G;
1068 //
1069 //calculate the solution and gradients at quadrature points
1070 //
1071 // u_ext=0.0;
1072 // for (int I=0;I<nSpace;I++)
1073 // {
1074 // grad_u_ext[I] = 0.0;
1075 // bc_grad_u_ext[I] = 0.0;
1076 // }
1077 // for (int j=0;j<nDOF_trial_element;j++)
1078 // {
1079 // int eN_j = eN*nDOF_trial_element+j,
1080 // ebNE_kb_j = ebNE_kb*nDOF_trial_element+j,
1081 // ebNE_kb_j_nSpace= ebNE_kb_j*nSpace;
1082 // u_ext += valFromDOF_c(u_dof[u_l2g[eN_j]],u_trial_ext[ebNE_kb_j]);
1083
1084 // for (int I=0;I<nSpace;I++)
1085 // {
1086 // grad_u_ext[I] += gradFromDOF_c(u_dof[u_l2g[eN_j]],u_grad_trial_ext[ebNE_kb_j_nSpace+I]);
1087 // }
1088 // }
1089 ck.calculateMapping_elementBoundary(eN,
1090 ebN_local,
1091 kb,
1092 ebN_local_kb,
1093 mesh_dof.data(),
1094 mesh_l2g.data(),
1095 mesh_trial_trace_ref.data(),
1096 mesh_grad_trial_trace_ref.data(),
1097 boundaryJac_ref.data(),
1098 jac_ext,
1099 jacDet_ext,
1100 jacInv_ext,
1101 boundaryJac,
1102 metricTensor,
1103 metricTensorDetSqrt,
1104 normal_ref.data(),
1105 normal,
1106 x_ext,y_ext,z_ext);
1107 ck.calculateMappingVelocity_elementBoundary(eN,
1108 ebN_local,
1109 kb,
1110 ebN_local_kb,
1111 mesh_velocity_dof.data(),
1112 mesh_l2g.data(),
1113 mesh_trial_trace_ref.data(),
1114 xt_ext,yt_ext,zt_ext,
1115 normal,
1116 boundaryJac,
1117 metricTensor,
1118 integralScaling);
1119 dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref[kb];
1120 //dS = metricTensorDetSqrt*dS_ref[kb];
1121 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
1122 //compute shape and solution information
1123 //shape
1124 ck.gradTrialFromRef(&u_grad_trial_trace_ref[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,u_grad_trial_trace);
1125 //solution and gradients
1126 ck.valFromDOF(u_dof.data(),&u_l2g[eN_nDOF_trial_element],&u_trial_trace_ref[ebN_local_kb*nDOF_test_element],u_ext);
1127 ck.gradFromDOF(u_dof.data(),&u_l2g[eN_nDOF_trial_element],u_grad_trial_trace,grad_u_ext);
1128 //precalculate test function products with integration weights
1129 for (int j=0;j<nDOF_trial_element;j++)
1130 {
1131 u_test_dS[j] = u_test_trace_ref[ebN_local_kb*nDOF_test_element+j]*dS;
1132 }
1133 //
1134 //load the boundary values
1135 //
1136 bc_u_ext = isDOFBoundary_u[ebNE_kb]*ebqe_bc_u_ext[ebNE_kb]+(1-isDOFBoundary_u[ebNE_kb])*ebqe_rd_u_ext[ebNE_kb];
1137 //
1138 //calculate the internal and external trace of the pde coefficients
1139 //
1140 evaluateCoefficients(&ebqe_velocity_ext[ebNE_kb_nSpace],
1141 u_ext,
1142 grad_u_ext,
1143 m_ext,
1144 dm_ext,
1145 H_ext,
1146 dH_ext);
1147 evaluateCoefficients(&ebqe_velocity_ext[ebNE_kb_nSpace],
1148 bc_u_ext,
1149 bc_grad_u_ext,
1150 bc_m_ext,
1151 bc_dm_ext,
1152 bc_H_ext,
1153 bc_dH_ext);
1154 //
1155 //moving domain
1156 //
1157 double velocity_ext[nSpace];
1158 double mesh_velocity[3];
1159 mesh_velocity[0] = xt_ext;
1160 mesh_velocity[1] = yt_ext;
1161 mesh_velocity[2] = zt_ext;
1162 for (int I=0;I<nSpace;I++)
1163 {
1164 velocity_ext[I] = - MOVING_DOMAIN*mesh_velocity[I];
1165 }
1166 //
1167 //calculate the numerical fluxes
1168 //
1169 exteriorNumericalFluxDerivative(normal,
1170 dH_ext,
1171 velocity_ext,
1172 dflux_u_u_ext);
1173 //
1174 //calculate the flux jacobian
1175 //
1176 for (int j=0;j<nDOF_trial_element;j++)
1177 {
1178 //int ebNE_kb_j = ebNE_kb*nDOF_trial_element+j;
1179 int ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j;
1180
1181 fluxJacobian_u_u[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_u_u_ext,u_trial_trace_ref[ebN_local_kb_j]);
1182 }//j
1183 //
1184 //update the global Jacobian from the flux Jacobian
1185 //
1186 for (int i=0;i<nDOF_test_element;i++)
1187 {
1188 int eN_i = eN*nDOF_test_element+i;
1189 //int ebNE_kb_i = ebNE_kb*nDOF_test_element+i;
1190 for (int j=0;j<nDOF_trial_element;j++)
1191 {
1192 int ebN_i_j = ebN*4*nDOF_test_X_trial_element + i*nDOF_trial_element + j;
1193 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] += fluxJacobian_u_u[j]*u_test_dS[i];
1194 //globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] += MOVING_DOMAIN*fluxJacobian_u_u[j]*u_test_dS[i];
1195 }//j
1196 }//i
1197 }//kb
1198 }//ebNE
1199 }//computeJacobian
1200
1201 void calculateWaterline(arguments_dict& args)
1202 {
1203 xt::pyarray<int>& wlc = args.array<int>("wlc");
1204 xt::pyarray<double>& waterline = args.array<double>("waterline");
1205 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1206 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1207 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1208 xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
1209 double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
1210 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1211 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1212 xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1213 xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1214 xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1215 xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
1216 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
1217 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
1218 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
1219 xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
1220 xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
1221 xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
1222 xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
1223 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
1224 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
1225 int nElements_global = args.scalar<int>("nElements_global");
1226 double useMetrics = args.scalar<double>("useMetrics");
1227 double alphaBDF = args.scalar<double>("alphaBDF");
1228 int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
1229 double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
1230 double sc_uref = args.scalar<double>("sc_uref");
1231 double sc_alpha = args.scalar<double>("sc_alpha");
1232 xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1233 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1234 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1235 xt::pyarray<double>& u_dof_old = args.array<double>("u_dof_old");
1236 xt::pyarray<double>& velocity = args.array<double>("velocity");
1237 xt::pyarray<double>& q_m = args.array<double>("q_m");
1238 xt::pyarray<double>& q_u = args.array<double>("q_u");
1239 xt::pyarray<double>& q_n = args.array<double>("q_n");
1240 xt::pyarray<double>& q_dH = args.array<double>("q_dH");
1241 xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
1242 xt::pyarray<double>& cfl = args.array<double>("cfl");
1243 xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
1244 xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
1245 int offset_u = args.scalar<int>("offset_u");
1246 int stride_u = args.scalar<int>("stride_u");
1247 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
1248 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
1249 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
1250 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
1251 xt::pyarray<int>& elementBoundaryMaterialTypes = args.array<int>("elementBoundaryMaterialTypes");
1252 xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
1253 xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
1254 xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
1255 xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
1256 // Tetrehedral elements specific extraction routine for waterline extraction
1257 // Loops over boundaries and checks if boundary is infact a hull (hardwired check if mattype > 6)
1258 // Extracts the nodal values of boundary triangle (4th point is dropped = hardwired assumption we are dealing with linear tet)
1259 // Then computes an average value and position for both negative and positive values
1260 // If both positive and negative values re found, and we are actually in a triangle containing the interface
1261 // a linear interpolation of negative and positive average is reported as interface location (exact in case of linear tets)
1262
1263 for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
1264 {
1265 int ebN = exteriorElementBoundariesArray[ebNE];
1266 int eN = elementBoundaryElementsArray[ebN*2+0];
1267 int bN = elementBoundaryLocalElementBoundariesArray[ebN*2+0];
1268
1269 if (elementBoundaryMaterialTypes[ebN] >6)
1270 {
1271
1272 double val,x,y,z;
1273 int pos=0, neg=0;
1274 double xn=0.0, yn=0.0, zn=0.0, vn=0.0;
1275 double xp=0.0, yp=0.0, zp=0.0, vp=0.0;
1276
1277 for (int j=0;j<nDOF_trial_element;j++)
1278 {
1279 if (j != bN) {
1280 int eN_nDOF_trial_element_j = eN*nDOF_trial_element + j;
1281 int eN_nDOF_mesh_trial_element_j = eN*nDOF_mesh_trial_element + j;
1282 val = u_dof[u_l2g[eN_nDOF_trial_element_j]];
1283 x = mesh_dof[mesh_l2g[eN_nDOF_mesh_trial_element_j]*3+0];
1284 y = mesh_dof[mesh_l2g[eN_nDOF_mesh_trial_element_j]*3+1];
1285 z = mesh_dof[mesh_l2g[eN_nDOF_mesh_trial_element_j]*3+2];
1286
1287 if (val < 0.0)
1288 {
1289 neg++;
1290 vn+=val;
1291 xn+=x;
1292 yn+=y;
1293 zn+=z;
1294 }
1295 else
1296 {
1297 pos++;
1298 vp+=val;
1299 xp+=x;
1300 yp+=y;
1301 zp+=z;
1302 }
1303 }
1304 } // trail for
1305
1306
1307 if ((pos > 0) && (neg > 0) )
1308 {
1309 vp /= pos;
1310 vn /= neg;
1311
1312 double alpha = vp/(vp -vn);
1313
1314 waterline[wlc[0]*3 + 0] = alpha*(xn/neg) + (1.0-alpha)*(xp/pos);
1315 waterline[wlc[0]*3 + 1] = alpha*(yn/neg) + (1.0-alpha)*(yp/pos);
1316 waterline[wlc[0]*3 + 2] = alpha*(zn/neg) + (1.0-alpha)*(zp/pos);
1317 wlc[0]++;
1318
1319 } // end value if
1320
1321 } // end bnd mat check
1322
1323 }//ebNE
1324
1325 //std::cout<<"CPP WLC "<<wlc[0]<<std::endl;
1326 } // calcWaterline
1327
1328 };//cppNCLS3P
1329
1330 inline cppNCLS3P_base* newNCLS3P(int nSpaceIn,
1331 int nQuadraturePoints_elementIn,
1332 int nDOF_mesh_trial_elementIn,
1333 int nDOF_trial_elementIn,
1334 int nDOF_test_elementIn,
1335 int nQuadraturePoints_elementBoundaryIn,
1336 int CompKernelFlag)
1337 {
1338 if (nSpaceIn == 2)
1339 return proteus::chooseAndAllocateDiscretization2D<cppNCLS3P_base,cppNCLS3P,CompKernel>(nSpaceIn,
1340 nQuadraturePoints_elementIn,
1341 nDOF_mesh_trial_elementIn,
1342 nDOF_trial_elementIn,
1343 nDOF_test_elementIn,
1344 nQuadraturePoints_elementBoundaryIn,
1345 CompKernelFlag);
1346 else
1347 return proteus::chooseAndAllocateDiscretization<cppNCLS3P_base,cppNCLS3P,CompKernel>(nSpaceIn,
1348 nQuadraturePoints_elementIn,
1349 nDOF_mesh_trial_elementIn,
1350 nDOF_trial_elementIn,
1351 nDOF_test_elementIn,
1352 nQuadraturePoints_elementBoundaryIn,
1353 CompKernelFlag);
1354 /* return proteus::chooseAndAllocateDiscretization<cppNCLS3P_base,cppNCLS3P>(nSpaceIn, */
1355 /* nQuadraturePoints_elementIn, */
1356 /* nDOF_mesh_trial_elementIn, */
1357 /* nDOF_trial_elementIn, */
1358 /* nDOF_test_elementIn, */
1359 /* nQuadraturePoints_elementBoundaryIn, */
1360 /* CompKernelFlag); */
1361 }
1362}//proteus
1363#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
cMax
#define cMax
Definition NCLS3P.h:11
cE
#define cE
Definition NCLS3P.h:10
proteus::cppNCLS3P_base
Definition NCLS3P.h:29
proteus::cppNCLS3P_base::~cppNCLS3P_base
virtual ~cppNCLS3P_base()
Definition NCLS3P.h:32
proteus::cppNCLS3P_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)=0
proteus::cppNCLS3P_base::calculateWaterline
virtual void calculateWaterline(arguments_dict &args)=0
proteus::cppNCLS3P_base::calculateResidual
virtual void calculateResidual(arguments_dict &args)=0
proteus::cppNCLS3P::calculateCFL
void calculateCFL(const double &elementDiameter, const double df[nSpace], double &cfl)
Definition NCLS3P.h:56
proteus::cppNCLS3P::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 NCLS3P.h:69
proteus::cppNCLS3P::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 NCLS3P.h:107
proteus::cppNCLS3P::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 NCLS3P.h:122
proteus::cppNCLS3P::calculateJacobian
void calculateJacobian(arguments_dict &args)
Definition NCLS3P.h:768
proteus::cppNCLS3P::exteriorNumericalFluxDerivative
void exteriorNumericalFluxDerivative(const double n[nSpace], const double velocity[nSpace], const double velocity_movingDomain[nSpace], double &dflux)
Definition NCLS3P.h:149
proteus::cppNCLS3P::calculateWaterline
void calculateWaterline(arguments_dict &args)
Definition NCLS3P.h:1201
proteus::cppNCLS3P::nDOF_test_X_trial_element
const int nDOF_test_X_trial_element
Definition NCLS3P.h:48
proteus::cppNCLS3P::calculateResidual
void calculateResidual(arguments_dict &args)
Definition NCLS3P.h:171
proteus::cppNCLS3P::ck
CompKernelType ck
Definition NCLS3P.h:49
proteus::cppNCLS3P::calculateSubgridError_tau
void calculateSubgridError_tau(const double &elementDiameter, const double &dmt, const double dH[nSpace], double &cfl, double &tau)
Definition NCLS3P.h:88
df
double df(double C, double b, double a, int q, int r)
Definition analyticalSolutions.c:2209
sign
#define sign(x, y)
Definition jf.h:44
Sign
#define Sign(x)
Definition jf.h:45
proteus
Definition ADR.h:19
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::newNCLS3P
cppNCLS3P_base * newNCLS3P(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition NCLS3P.h:1330
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