proteus 1.9.0
C/C++/Fortran libraries
Loading...
Searching...
No Matches
RDLS.h
Go to the documentation of this file.
1#ifndef RDLS_H
2#define RDLS_H
3#include <cmath>
4#include <iostream>
5#include <valarray>
6#include "CompKernel.h"
7#include "ModelFactory.h"
8#include "equivalent_polynomials.h"
9#include "ArgumentsDict.h"
10#include "xtensor-python/pyarray.hpp"
11
12namespace py = pybind11;
13
14#define SINGLE_POTENTIAL 1
15
16namespace proteus
17{
18 inline double heaviside(const double &z)
19 {
20 return (z>0 ? 1. : (z<0 ? 0. : 0.5));
21 }
22
23 template<int nSpace, int nP, int nQ, int nEBQ>
24 // The trailing flag is the IFEM gate. It has always been false here -- by
25 // default rather than by statement -- and must stay false: with it set, an
26 // element whose interface passes through an edge or corner node takes the
27 // IFEM branch in Simplex::set_quad, which forces D to 0 and H/ImH to a hard
28 // 0/1 instead of the moment fit. That deletes the interface measure this
29 // model integrates over. Written out so the choice is visible at the call
30 // site and cannot change underneath us if the template default changes.
31 using GeneralizedFunctions = equivalent_polynomials::GeneralizedFunctions_mix<nSpace, nP, nP, nQ, nEBQ, false>;
32
33
34 class RDLS_base
35 {
36 public:
37 std::valarray<double> weighted_lumped_mass_matrix;
38 virtual ~RDLS_base(){}
39 virtual void calculateResidual(arguments_dict& args, bool useExact)=0;
40 virtual void calculateJacobian(arguments_dict& args, bool useExact)=0;
41 virtual void calculateResidual_ellipticRedist(arguments_dict& args, bool useExact)=0;
42 virtual void calculateJacobian_ellipticRedist(arguments_dict& args, bool useExact)=0;
43 virtual void normalReconstruction(arguments_dict& args)=0;
44 virtual std::tuple<double, double, double> calculateMetricsAtEOS(arguments_dict& args)=0;
45 };
46
47 template<class CompKernelType,
48 int nSpace,
49 int nQuadraturePoints_element,
50 int nDOF_mesh_trial_element,
51 int nDOF_trial_element,
52 int nDOF_test_element,
53 int nQuadraturePoints_elementBoundary>
54 class RDLS : public RDLS_base
55 {
56 public:
57 const int nDOF_test_X_trial_element;
58 CompKernelType ck;
59 GeneralizedFunctions<nSpace,2,nQuadraturePoints_element,nQuadraturePoints_elementBoundary> gf,gfu;
60 RDLS():
61 nDOF_test_X_trial_element(nDOF_test_element*nDOF_trial_element),
62 ck()
63 {}
64
65 inline
66 void evaluateCoefficients(const double& eps,
67 const double& u_levelSet,
68 const double& u,
69 const double grad_u[nSpace],
70 double& m,
71 double& dm,
72 double& H,
73 double dH[nSpace],
74 double& r)
75 {
76 int I;
77 double normGradU=0.0,Si=0.0;
78 m = u;
79 dm=1.0;
80 H = 0.0;
81 Si= gf.H(eps,u_levelSet) - gf.ImH(eps,u_levelSet);
82 /* if (u_levelSet > 0.0) */
83 /* Si=1.0; */
84 /* else if (u_levelSet < 0.0) */
85 /* Si = -1.0; */
86 /* else */
87 /* Si=0.0; */
88 r = -Si;
89 for (I=0; I < nSpace; I++)
90 {
91 normGradU += grad_u[I]*grad_u[I];
92 }
93 normGradU = sqrt(normGradU);
94 H = Si*normGradU;
95 for (I=0; I < nSpace; I++)
96 {
97 dH[I] = Si*grad_u[I]/(normGradU+1.0e-12);
98 }
99 }
100
101 inline
102 void calculateSubgridError_tau(const double& elementDiameter,
103 const double& dmt,
104 const double dH[nSpace],
105 double& cfl,
106 double& tau)
107 {
108 double h,nrm_v,oneByAbsdt;
109 h = elementDiameter;
110 nrm_v=0.0;
111 for(int I=0;I<nSpace;I++)
112 nrm_v+=dH[I]*dH[I];
113 nrm_v = sqrt(nrm_v);
114 cfl = nrm_v/h;
115 oneByAbsdt = dmt;
116 //'1'
117 //tau = 1.0/(2.0*nrm_v/h + oneByAbsdt + 1.0e-8);
118 //'2'
119 tau = 1.0/sqrt(4.0*nrm_v*nrm_v/(h*h) + oneByAbsdt*oneByAbsdt + 1.0e-8);
120 }
121
122 inline
123 void calculateSubgridError_tau( const double G[nSpace*nSpace],
124 const double Ai[nSpace],
125 double& tau_v,
126 double& q_cfl)
127 {
128 double v_d_Gv=0.0;
129 for(int I=0;I<nSpace;I++)
130 for (int J=0;J<nSpace;J++)
131 v_d_Gv += Ai[I]*G[I*nSpace+J]*Ai[J];
132
133 tau_v = 1.0/(sqrt(v_d_Gv) + 1.0e-8);
134 }
135
136#undef CKDEBUG
137 void calculateResidual(arguments_dict& args,
138 bool useExact)
139 {
140 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
141 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
142 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
143 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
144 xt::pyarray<double>& x_ref = args.array<double>("x_ref");
145 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
146 xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
147 xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
148 xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
149 xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
150 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
151 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
152 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
153 xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
154 xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
155 xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
156 xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
157 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
158 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
159 int nElements_global = args.scalar<int>("nElements_global");
160 double useMetrics = args.scalar<double>("useMetrics");
161 double alphaBDF = args.scalar<double>("alphaBDF");
162 double epsFact_redist = args.scalar<double>("epsFact_redist");
163 double backgroundDiffusionFactor = args.scalar<double>("backgroundDiffusionFactor");
164 double weakDirichletFactor = args.scalar<double>("weakDirichletFactor");
165 int freezeLevelSet = args.scalar<int>("freezeLevelSet");
166 int useTimeIntegration = args.scalar<int>("useTimeIntegration");
167 int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
168 int lag_subgridError = args.scalar<int>("lag_subgridError");
169 double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
170 xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
171 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
172 xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
173 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
174 xt::pyarray<double>& phi_dof = args.array<double>("phi_dof");
175 xt::pyarray<double>& phi_ls = args.array<double>("phi_ls");
176 xt::pyarray<double>& q_m = args.array<double>("q_m");
177 xt::pyarray<double>& q_u = args.array<double>("q_u");
178 xt::pyarray<double>& q_n = args.array<double>("q_n");
179 xt::pyarray<double>& q_dH = args.array<double>("q_dH");
180 xt::pyarray<double>& u_weak_internal_bc_dofs = args.array<double>("u_weak_internal_bc_dofs");
181 xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
182 xt::pyarray<double>& q_dH_last = args.array<double>("q_dH_last");
183 xt::pyarray<double>& q_cfl = args.array<double>("q_cfl");
184 xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
185 xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
186 xt::pyarray<int>& weakDirichletConditionFlags = args.array<int>("weakDirichletConditionFlags");
187 int offset_u = args.scalar<int>("offset_u");
188 int stride_u = args.scalar<int>("stride_u");
189 xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
190 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
191 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
192 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
193 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
194 xt::pyarray<double>& ebqe_phi_ls_ext = args.array<double>("ebqe_phi_ls_ext");
195 xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
196 xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
197 xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
198 xt::pyarray<double>& ebqe_n = args.array<double>("ebqe_n");
199 int ELLIPTIC_REDISTANCING = args.scalar<int>("ELLIPTIC_REDISTANCING");
200 double backgroundDissipationEllipticRedist = args.scalar<double>("backgroundDissipationEllipticRedist");
201 xt::pyarray<double>& lumped_qx = args.array<double>("lumped_qx");
202 xt::pyarray<double>& lumped_qy = args.array<double>("lumped_qy");
203 xt::pyarray<double>& lumped_qz = args.array<double>("lumped_qz");
204 double alpha = args.scalar<double>("alpha");
205 double circbc=0.0,circ=0.0;
206 gf.useExact=useExact;
207 gfu.useExact=useExact;
208#ifdef CKDEBUG
209 std::cout<<"stuff"<<"\t"
210 <<alphaBDF<<"\t"
211 <<epsFact_redist<<"\t"
212 <<freezeLevelSet<<"\t"
213 <<useTimeIntegration<<"\t"
214 <<lag_shockCapturing<< "\t"
215 <<lag_subgridError<<std::endl;
216#endif
217 //
218 //loop over elements to compute volume integrals and load them into element and global residual
219 //
220 //eN is the element index
221 //eN_k is the quadrature point index for a scalar
222 //eN_k_nSpace is the quadrature point index for a vector
223 //eN_i is the element test function index
224 //eN_j is the element trial function index
225 //eN_k_j is the quadrature point index for a trial function
226 //eN_k_i is the quadrature point index for a trial function
227 double timeIntegrationScale = 1.0;
228 if (useTimeIntegration == 0)
229 timeIntegrationScale = 0.0;
230 double lag_shockCapturingScale = 1.0;
231 if (lag_shockCapturing == 0)
232 lag_shockCapturingScale = 0.0;
233 for(int eN=0;eN<nElements_global;eN++)
234 {
235 //declare local storage for element residual and initialize
236 int dummy_l2g[nDOF_mesh_trial_element];
237 double elementResidual_u[nDOF_test_element],element_phi[nDOF_trial_element];
238 double epsilon_redist,h_phi, dir[nSpace], norm;
239 for (int i=0;i<nDOF_test_element;i++)
240 {
241 int eN_i=eN*nDOF_trial_element+i;
242 elementResidual_u[i]=0.0;
243 element_phi[i] = phi_dof.data()[u_l2g.data()[eN_i]];
244 dummy_l2g[i] = i;
245 }//i
246 double element_nodes[nDOF_mesh_trial_element*3];
247 for (int i=0;i<nDOF_mesh_trial_element;i++)
248 {
249 int eN_i=eN*nDOF_mesh_trial_element+i;
250 for(int I=0;I<3;I++)
251 element_nodes[i*3 + I] = mesh_dof.data()[mesh_l2g.data()[eN_i]*3 + I];
252 }//i
253 gf.calculate(element_phi, element_nodes, x_ref.data(),false);
254 /* for (int i=0;i<nDOF_test_element;i++) */
255 /* { */
256 /* int eN_i=eN*nDOF_trial_element+i; */
257 /* double eps=1.0e-4; */
258 /* if((fabs(gf.exact.phi_dof_corrected[i]) < eps) || */
259 /* (fabs(phi_dof.data()[u_l2g.data()[eN_i]]) < eps)) */
260 /* std::cout<<"Warning "<<gf.exact.phi_dof_corrected[i]<<'\t' */
261 /* <<phi_dof.data()[u_l2g.data()[eN_i]]<<'\t' */
262 /* <<u_l2g.data()[eN_i]<<std::endl; */
263 /* } */
264 //loop over quadrature points and compute integrands
265 for (int k=0;k<nQuadraturePoints_element;k++)
266 {
267 gf.set_quad(k);
268 //compute indeces and declare local storage
269 int eN_k = eN*nQuadraturePoints_element+k,
270 eN_k_nSpace = eN_k*nSpace,
271 eN_nDOF_trial_element = eN*nDOF_trial_element;
272 double u=0.0,grad_u[nSpace],u0=0.0, grad_phi[nSpace],
273 m=0.0,dm=0.0,
274 H=0.0,dH[nSpace],
275 m_t=0.0,dm_t=0.0,
276 r=0.0,
277 dH_tau[nSpace],//dH if not lagging or q_dH_last if lagging tau
278 dH_strong[nSpace],//dH if not lagging or q_dH_last if lagging strong residual and adjoint
279 pdeResidual_u=0.0,
280 Lstar_u[nDOF_test_element],
281 subgridError_u=0.0,
282 tau=0.0,tau0=0.0,tau1=0.0,
283 numDiff0=0.0,numDiff1=0.0,
284 nu_sc=0.0,
285 jac[nSpace*nSpace],
286 jacDet,
287 jacInv[nSpace*nSpace],
288 u_grad_trial[nDOF_trial_element*nSpace],
289 u_test_dV[nDOF_trial_element],
290 u_grad_test_dV[nDOF_test_element*nSpace],
291 dV,x,y,z,
292 G[nSpace*nSpace],G_dd_G,tr_G;
293 ck.calculateMapping_element(eN,
294 k,
295 mesh_dof.data(),
296 mesh_l2g.data(),
297 mesh_trial_ref.data(),
298 mesh_grad_trial_ref.data(),
299 jac,
300 jacDet,
301 jacInv,
302 x,y,z);
303 ck.calculateH_element(eN,
304 k,
305 nodeDiametersArray.data(),
306 mesh_l2g.data(),
307 mesh_trial_ref.data(),
308 h_phi);
309 //get the physical integration weight
310 dV = fabs(jacDet)*dV_ref.data()[k];
311 ck.calculateG(jacInv,G,G_dd_G,tr_G);
312
313 //get the trial function gradients
314 ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
315 //get the solution
316 ck.valFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],u);
317 ck.valFromDOF(gf.exact.phi_dof_corrected,dummy_l2g,&u_trial_ref.data()[k*nDOF_trial_element],u0);
318 if (freezeLevelSet)
319 {
320 u0 = phi_ls.data()[eN_k];
321 }
322 /* double DX=(x-0.5),DY=(y-0.75); */
323 /* double radius = std::sqrt(DX*DX+DY*DY); */
324 /* double theta = std::atan2(DY,DX); */
325 /* double kp=10.0, scaling=1.0, rp=1; */
326 /* u0 = scaling*std::pow((0.15+(0.015/2.)*std::cos(kp*theta) - radius),rp); */
327
328 //u0 = 0.15 - std::sqrt(DX*DX + DY*DY);
329 //u0 = phi_ls.data()[eN_k];//cek debug--set to exact input
330 /* //cek debug */
331 /* double u0_test; */
332 /* ck.valFromDOF(phi_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],u0_test); */
333 /* if (u0 != u0_test) */
334 /* std::cout<<"eN "<<eN<<" k "<<k<<" u0 "<<u0<<" u0_test "<<u0_test<<std::endl; */
335 //get the solution gradients
336 ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial,grad_u);
337 //precalculate test function products with integration weights
338 for (int j=0;j<nDOF_trial_element;j++)
339 {
340 u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
341 for (int I=0;I<nSpace;I++)
342 {
343 u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
344 }
345 }
346 //
347 //calculate pde coefficients at quadrature points
348 //
349 /* norm = 1.0e-8; */
350 /* for (int I=0;I<nSpace;I++) */
351 /* norm += grad_u[I]*grad_u[I]; */
352 /* norm = sqrt(norm); */
353
354 /* for (int I=0;I<nSpace;I++) */
355 /* dir[I] = grad_u[I]/norm; */
356
357 /* ck.calculateGScale(G,dir,h_phi); */
358
359 epsilon_redist = epsFact_redist*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
360
361 evaluateCoefficients(epsilon_redist,
362 u0,
363 u,
364 grad_u,
365 m,
366 dm,
367 H,
368 dH,
369 r);
370 //TODO allow not lagging of subgrid error etc,
371 //remove conditional?
372 //default no lagging
373 for (int I=0; I < nSpace; I++)
374 {
375 dH_tau[I] = dH[I];
376 dH_strong[I] = dH[I];
377 }
378 if (lag_subgridError > 0)
379 {
380 for (int I=0; I < nSpace; I++)
381 {
382 dH_tau[I] = q_dH_last.data()[eN_k_nSpace+I];
383 }
384 }
385 if (lag_subgridError > 1)
386 {
387 for (int I=0; I < nSpace; I++)
388 {
389 dH_strong[I] = q_dH_last.data()[eN_k_nSpace+I];
390 }
391 }
392 //save mass for time history and dH for subgrid error
393 //save solution for other models
394 //
395 q_m.data()[eN_k] = m;
396 q_u.data()[eN_k] = u;
397 for (int I=0;I<nSpace;I++)
398 q_n.data()[eN_k_nSpace+I] = dir[I];
399
400 for (int I=0;I<nSpace;I++)
401 {
402 int eN_k_nSpace_I = eN_k_nSpace+I;
403 q_dH.data()[eN_k_nSpace_I] = dH[I];
404 }
405
406 //
407 //moving mesh
408 //
409 //omit for now
410 //
411 //calculate time derivative at quadrature points
412 //
413 ck.bdf(alphaBDF,
414 q_m_betaBDF.data()[eN_k],
415 m,
416 dm,
417 m_t,
418 dm_t);
419
420 //TODO add option to skip if not doing time integration (Newton stead-state solve)
421 m *= timeIntegrationScale; dm *= timeIntegrationScale; m_t *= timeIntegrationScale;
422 dm_t *= timeIntegrationScale;
423#ifdef CKDEBUG
424 std::cout<<"alpha "<<alphaBDF<<"\t"<<q_m_betaBDF.data()[eN_k]<<"\t"<<m_t<<'\t'<<m<<'\t'<<alphaBDF*m<<std::endl;
425#endif
426 //
427 //calculate subgrid error (strong residual and adjoint)
428 //
429 //calculate strong residual
430 pdeResidual_u = ck.Mass_strong(m_t) +
431 ck.Hamiltonian_strong(dH_strong,grad_u) + //would need dH if not lagging
432 ck.Reaction_strong(r);
433#ifdef CKDEBUG
434 std::cout<<"dH_strong "<<dH_strong[0]<<'\t'<<dH_strong[1]<<'\t'<<dH_strong[2]<<std::endl;
435#endif
436 //calculate adjoint
437 for (int i=0;i<nDOF_test_element;i++)
438 {
439 //int eN_k_i_nSpace = (eN_k*nDOF_trial_element+i)*nSpace;
440 int i_nSpace=i*nSpace;
441 Lstar_u[i] = ck.Hamiltonian_adjoint(dH_strong,&u_grad_test_dV[i_nSpace]);
442 //reaction is constant
443 }
444 //calculate tau and tau*Res
445 calculateSubgridError_tau(elementDiameter.data()[eN],
446 dm_t,dH_tau,
447 q_cfl.data()[eN_k],
448 tau0);
449 calculateSubgridError_tau(G,
450 dH_tau,
451 tau1,
452 q_cfl.data()[eN_k]);
453
454 tau = useMetrics*tau1+(1.0-useMetrics)*tau0;
455
456 //std::cout<<tau<<std::endl;
457
458
459 subgridError_u = -tau*pdeResidual_u;
460 //
461 //calculate shock capturing diffusion
462 //
463 ck.calculateNumericalDiffusion(shockCapturingDiffusion,elementDiameter.data()[eN],pdeResidual_u,grad_u,numDiff0);
464 ck.calculateNumericalDiffusion(shockCapturingDiffusion,G,pdeResidual_u,grad_u,numDiff1);
465
466 q_numDiff_u.data()[eN_k] = useMetrics*numDiff1+(1.0-useMetrics)*numDiff0;
467
468#ifdef CKDEBUG
469 std::cout<<"q_numDiff_u[eN_k] "<<q_numDiff_u.data()[eN_k]<<" q_numDiff_u_last[eN_k] "<<q_numDiff_u_last.data()[eN_k]<<" lag "<<lag_shockCapturingScale<<std::endl;
470#endif
471 nu_sc = q_numDiff_u.data()[eN_k]*(1.0-lag_shockCapturingScale) + q_numDiff_u_last.data()[eN_k]*lag_shockCapturingScale;
472 /* double epsilon_background_diffusion = 3.0*h_phi;//2.0*epsFact_redist*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]); */
473 /* if (fabs(u0) > epsilon_background_diffusion) */
474 /* nu_sc += backgroundDiffusionFactor*h_phi; */
475
476 double epsilon_background_diffusion = 2.0*epsFact_redist*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
477 if (fabs(phi_ls.data()[eN_k]) > epsilon_background_diffusion)
478 nu_sc += backgroundDiffusionFactor*elementDiameter.data()[eN]; //
479 //update element residual
480 //
481 for(int i=0;i<nDOF_test_element;i++)
482 {
483 int i_nSpace = i*nSpace;
484 double FREEZE=double(freezeLevelSet);
485 //assert(FREEZE==0.0);
486 //int eN_k_i=eN_k*nDOF_test_element+i;
487 //int eN_k_i_nSpace = eN_k_i*nSpace;
488
489#ifdef CKDEBUG
490 std::cout<<"shock capturing input nu_sc "<<nu_sc<<'\t'<<grad_u[0]<<'\t'<<grad_u[1]<<'\t'<<grad_u[1]<<'\t'<<u_grad_test_dV[i_nSpace]<<std::endl;
491#endif
492 //std::cout<<element_phi[i]<<'\t'<<element_nodes[i*3+0]<<'\t'<<element_nodes[i*3+1]<<std::endl;
493 //std::cout<<"eN_k "<<eN_k<<" D "<<gf.D(epsilon_redist,phi_ls.data()[eN_k])<<std::endl;
494 circbc += ck.Reaction_weak(gf.D(epsilon_redist,u0)*(u-u0), u_test_dV[i]);
495 elementResidual_u[i] += ck.Mass_weak(m_t,u_test_dV[i]) +
496 ck.Hamiltonian_weak(H,u_test_dV[i]) +
497 ck.Reaction_weak(r,u_test_dV[i]) +
498 (1.0-FREEZE)*(weakDirichletFactor/h_phi)*ck.Reaction_weak(gf.D(epsilon_redist,u0)*(u0-u),
499 u_test_dV[i]) +
500 ck.SubgridError(subgridError_u,Lstar_u[i]) +
501 ck.NumericalDiffusion(nu_sc,grad_u,&u_grad_test_dV[i_nSpace]);
502#ifdef CKDEBUG
503 std::cout<<ck.Mass_weak(m_t,u_test_dV[i])<<'\t'
504 <<ck.Hamiltonian_weak(H,u_test_dV[i]) <<'\t'
505 <<ck.Reaction_weak(r,u_test_dV[i])<<'\t'
506 <<ck.SubgridError(subgridError_u,Lstar_u[i])<<'\t'
507 <<ck.NumericalDiffusion(nu_sc,grad_u,&u_grad_test_dV[i_nSpace])<<std::endl;
508#endif
509 }//i
510 //
511 }//k
512 //
513 //apply weak constraints for unknowns near zero level set
514 //
515 //
516 if (freezeLevelSet)
517 {
518 for (int j = 0; j < nDOF_trial_element; j++)
519 {
520 const int eN_j = eN*nDOF_trial_element+j;
521 const int J = u_l2g.data()[eN_j];
522 //if (weakDirichletConditionFlags.data()[J] == 1)
523 if (fabs(u_weak_internal_bc_dofs.data()[J]) < epsilon_redist)
524 {
525 elementResidual_u[j] = (u_dof.data()[J]-u_weak_internal_bc_dofs.data()[J])*weakDirichletFactor*elementDiameter.data()[eN];
526 }
527 }//j
528 }//freeze
529
530 //
531 //load element into global residual and save element residual
532 //
533 for(int i=0;i<nDOF_test_element;i++)
534 {
535 int eN_i=eN*nDOF_test_element+i;
536
537#ifdef CKDEBUG
538 std::cout<<"element residual i = "<<i<<"\t"<<elementResidual_u[i]<<std::endl;
539#endif
540 globalResidual.data()[offset_u+stride_u*u_l2g.data()[eN_i]]+=elementResidual_u[i];
541 }//i
542 }//elements
543 //
544 //loop over exterior element boundaries
545 //
546 //ebNE is the Exterior element boundary INdex
547 //ebN is the element boundary INdex
548 //eN is the element index
549 for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
550 {
551 int ebN = exteriorElementBoundariesArray.data()[ebNE],
552 eN = elementBoundaryElementsArray.data()[ebN*2+0],
553 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
554 eN_nDOF_trial_element = eN*nDOF_trial_element;
555 double epsilon_redist, h_phi;
556 for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
557 {
558 int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
559 ebNE_kb_nSpace = ebNE_kb*nSpace,
560 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
561 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
562 double u_ext=0.0,
563 grad_u_ext[nSpace],
564 jac_ext[nSpace*nSpace],
565 jacDet_ext,
566 jacInv_ext[nSpace*nSpace],
567 boundaryJac[nSpace*(nSpace-1)],
568 metricTensor[(nSpace-1)*(nSpace-1)],
569 metricTensorDetSqrt,
570 u_test_dS[nDOF_test_element],
571 u_grad_trial_trace[nDOF_trial_element*nSpace],
572 normal[nSpace],x_ext,y_ext,z_ext,
573 dir[nSpace],norm;
574 ck.calculateMapping_elementBoundary(eN,
575 ebN_local,
576 kb,
577 ebN_local_kb,
578 mesh_dof.data(),
579 mesh_l2g.data(),
580 mesh_trial_trace_ref.data(),
581 mesh_grad_trial_trace_ref.data(),
582 boundaryJac_ref.data(),
583 jac_ext,
584 jacDet_ext,
585 jacInv_ext,
586 boundaryJac,
587 metricTensor,
588 metricTensorDetSqrt,
589 normal_ref.data(),
590 normal,
591 x_ext,y_ext,z_ext);
592 //compute shape and solution information
593 //shape
594 ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,u_grad_trial_trace);
595 //solution and gradients
596 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);
597 ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial_trace,grad_u_ext);
598 norm = 1.0e-8;
599 for (int I=0;I<nSpace;I++)
600 norm += grad_u_ext[I]*grad_u_ext[I];
601 norm = sqrt(norm);
602 for (int I=0;I<nSpace;I++)
603 dir[I] = grad_u_ext[I]/norm;
604 //save for other models
605 ebqe_u.data()[ebNE_kb] = u_ext;
606 for (int I=0;I<nSpace;I++)
607 ebqe_n.data()[ebNE_kb_nSpace+I] = dir[I];
608 }//kb
609 }//ebNE
610 //std::cout<<"Circ Res"<<circbc<<'\t'<<circ<<std::endl;
611 }
612
613 //for now assumes that using time integration
614 //and so lags stabilization and subgrid error
615 //extern "C" void calculateJacobian_RDLSV2(int nElements_global,
616 void calculateJacobian(arguments_dict& args,
617 bool useExact)
618 {
619 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
620 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
621 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
622 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
623 xt::pyarray<double>& x_ref = args.array<double>("x_ref");
624 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
625 xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
626 xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
627 xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
628 xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
629 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
630 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
631 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
632 xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
633 xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
634 xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
635 xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
636 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
637 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
638 int nElements_global = args.scalar<int>("nElements_global");
639 double useMetrics = args.scalar<double>("useMetrics");
640 double alphaBDF = args.scalar<double>("alphaBDF");
641 double epsFact_redist = args.scalar<double>("epsFact_redist");
642 double backgroundDiffusionFactor = args.scalar<double>("backgroundDiffusionFactor");
643 double weakDirichletFactor = args.scalar<double>("weakDirichletFactor");
644 int freezeLevelSet = args.scalar<int>("freezeLevelSet");
645 int useTimeIntegration = args.scalar<int>("useTimeIntegration");
646 int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
647 int lag_subgridError = args.scalar<int>("lag_subgridError");
648 double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
649 xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
650 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
651 xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
652 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
653 xt::pyarray<double>& phi_dof = args.array<double>("phi_dof");
654 xt::pyarray<double>& u_weak_internal_bc_dofs = args.array<double>("u_weak_internal_bc_dofs");
655 xt::pyarray<double>& phi_ls = args.array<double>("phi_ls");
656 xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
657 xt::pyarray<double>& q_dH_last = args.array<double>("q_dH_last");
658 xt::pyarray<double>& q_cfl = args.array<double>("q_cfl");
659 xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
660 xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
661 xt::pyarray<int>& weakDirichletConditionFlags = args.array<int>("weakDirichletConditionFlags");
662 xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
663 xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
664 xt::pyarray<double>& globalJacobian = args.array<double>("globalJacobian");
665 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
666 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
667 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
668 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
669 xt::pyarray<double>& ebqe_phi_ls_ext = args.array<double>("ebqe_phi_ls_ext");
670 xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
671 xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
672 xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
673 int ELLIPTIC_REDISTANCING = args.scalar<int>("ELLIPTIC_REDISTANCING");
674 double backgroundDissipationEllipticRedist = args.scalar<double>("backgroundDissipationEllipticRedist");
675 double alpha = args.scalar<double>("alpha");
676 double circ=0.0;
677 gf.useExact=useExact;
678 //
679 //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
680 //
681 double timeIntegrationScale = 1.0;
682 if (useTimeIntegration == 0)
683 timeIntegrationScale = 0.0;
684 double lag_shockCapturingScale = 1.0;
685 if (lag_shockCapturing == 0)
686 lag_shockCapturingScale = 0.0;
687 for(int eN=0;eN<nElements_global;eN++)
688 {
689 int dummy_l2g[nDOF_mesh_trial_element];
690 double elementJacobian_u_u[nDOF_test_element][nDOF_trial_element],element_phi[nDOF_trial_element];
691 double epsilon_redist,h_phi, dir[nSpace], norm;
692 for (int i=0;i<nDOF_test_element;i++)
693 {
694 int eN_i=eN*nDOF_trial_element+i;
695 element_phi[i] = phi_dof.data()[u_l2g.data()[eN_i]];
696 dummy_l2g[i] = i;
697 for (int j=0;j<nDOF_trial_element;j++)
698 {
699 elementJacobian_u_u[i][j]=0.0;
700 }
701 }
702 double element_nodes[nDOF_mesh_trial_element*3];
703 for (int i=0;i<nDOF_mesh_trial_element;i++)
704 {
705 int eN_i=eN*nDOF_mesh_trial_element+i;
706 for(int I=0;I<3;I++)
707 element_nodes[i*3 + I] = mesh_dof.data()[mesh_l2g.data()[eN_i]*3 + I];
708 }//i
709 gf.calculate(element_phi, element_nodes, x_ref.data(),false);
710 for (int k=0;k<nQuadraturePoints_element;k++)
711 {
712 gf.set_quad(k);
713 int eN_k = eN*nQuadraturePoints_element+k, //index to a scalar at a quadrature point
714 eN_k_nSpace = eN_k*nSpace,
715 eN_nDOF_trial_element = eN*nDOF_trial_element; //index to a vector at a quadrature point
716
717 //declare local storage
718 double u=0.0,u0=0.0,
719 grad_u[nSpace],
720 m=0.0,dm=0.0,
721 H=0.0,dH[nSpace],
722 m_t=0.0,dm_t=0.0,r=0.0,
723 dH_tau[nSpace],//dH or dH_last if lagging for tau formula
724 dH_strong[nSpace],//dH or dH_last if lagging for strong residual and adjoint
725 dpdeResidual_u_u[nDOF_trial_element],
726 Lstar_u[nDOF_test_element],
727 dsubgridError_u_u[nDOF_trial_element],
728 tau=0.0,tau0=0.0,tau1=0.0,
729 nu_sc=0.0,
730 jac[nSpace*nSpace],
731 jacDet,
732 jacInv[nSpace*nSpace],
733 u_grad_trial[nDOF_trial_element*nSpace],
734 dV,
735 u_test_dV[nDOF_test_element],
736 u_grad_test_dV[nDOF_test_element*nSpace],
737 x,y,z,
738 G[nSpace*nSpace],G_dd_G,tr_G;
739 //
740 //calculate solution and gradients at quadrature points
741 //
742 ck.calculateMapping_element(eN,
743 k,
744 mesh_dof.data(),
745 mesh_l2g.data(),
746 mesh_trial_ref.data(),
747 mesh_grad_trial_ref.data(),
748 jac,
749 jacDet,
750 jacInv,
751 x,y,z);
752 ck.calculateH_element(eN,
753 k,
754 nodeDiametersArray.data(),
755 mesh_l2g.data(),
756 mesh_trial_ref.data(),
757 h_phi);
758 //get the physical integration weight
759 dV = fabs(jacDet)*dV_ref.data()[k];
760 ck.calculateG(jacInv,G,G_dd_G,tr_G);
761 //get the trial function gradients
762 ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
763 //get the solution
764 ck.valFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],u);
765 ck.valFromDOF(gf.exact.phi_dof_corrected,dummy_l2g,&u_trial_ref.data()[k*nDOF_trial_element],u0);
766 if (freezeLevelSet)
767 u0 = phi_ls.data()[eN_k];
768 //u0 = phi_ls.data()[eN_k];//cek debug
769 /* double DX=(x-0.5),DY=(y-0.75); */
770 /* double radius = std::sqrt(DX*DX+DY*DY); */
771 /* double theta = std::atan2(DY,DX); */
772 /* double kp=10.0, scaling=1.0, rp=1; */
773 /* u0 = scaling*std::pow((0.15+(0.015/2.)*std::cos(kp*theta) - radius),rp); */
774 //u0 = 0.15 - std::sqrt(DX*DX + DY*DY);
775 //u0 = phi_ls.data()[eN_k];//cek debug--set to exact input
776 /* double u0_test=0; */
777 /* ck.valFromDOF(phi_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],u0_test); */
778 /* if (u0 != u0_test) */
779 /* std::cout<<"JAC eN "<<eN<<" k "<<k<<" u0 "<<u0<<" u0_test "<<u0_test<<std::endl; */
780
781 //get the solution gradients
782 ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial,grad_u);
783 //precalculate test function products with integration weights
784 for (int j=0;j<nDOF_trial_element;j++)
785 {
786 u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
787 for (int I=0;I<nSpace;I++)
788 {
789 u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
790 }
791 }
792 //
793 //calculate pde coefficients and derivatives at quadrature points
794 //
795 /* norm = 1.0e-8; */
796 /* for (int I=0;I<nSpace;I++) */
797 /* norm += grad_u[I]*grad_u[I]; */
798 /* norm = sqrt(norm); */
799 /* for (int I=0;I<nSpace;I++) */
800 /* dir[I] = grad_u[I]/norm; */
801
802 /* ck.calculateGScale(G,dir,h_phi); */
803
804 epsilon_redist = epsFact_redist*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
805
806 evaluateCoefficients(epsilon_redist,
807 u0,
808 u,
809 grad_u,
810 m,
811 dm,
812 H,
813 dH,
814 r);
815 //TODO allow not lagging of subgrid error etc
816 //remove conditional?
817 //default no lagging
818 for (int I=0; I < nSpace; I++)
819 {
820 dH_tau[I] = dH[I];
821 dH_strong[I] = dH[I];
822 }
823 if (lag_subgridError > 0)
824 {
825 for (int I=0; I < nSpace; I++)
826 {
827 dH_tau[I] = q_dH_last.data()[eN_k_nSpace+I];
828 }
829 }
830 if (lag_subgridError > 1)
831 {
832 for (int I=0; I < nSpace; I++)
833 {
834 dH_strong[I] = q_dH_last.data()[eN_k_nSpace+I];
835 }
836 }
837 //
838 //moving mesh
839 //
840 //omit for now
841 //
842 //calculate time derivatives
843 //
844 ck.bdf(alphaBDF,
845 q_m_betaBDF.data()[eN_k],
846 m,
847 dm,
848 m_t,
849 dm_t);
850 //TODO add option to skip if not doing time integration (Newton stead-state solve)
851 m *= timeIntegrationScale; dm *= timeIntegrationScale; m_t *= timeIntegrationScale;
852 dm_t *= timeIntegrationScale;
853
854 //
855 //calculate subgrid error contribution to the Jacobian (strong residual, adjoint, jacobian of strong residual)
856 //
857 //calculate the adjoint times the test functions
858 for (int i=0;i<nDOF_test_element;i++)
859 {
860 //int eN_k_i_nSpace = (eN_k*nDOF_trial_element+i)*nSpace;
861 int i_nSpace=i*nSpace;
862
863 Lstar_u[i]=ck.Hamiltonian_adjoint(dH_strong,&u_grad_test_dV[i_nSpace]);
864
865 }
866 //calculate the Jacobian of strong residual
867 for (int j=0;j<nDOF_trial_element;j++)
868 {
869 //int eN_k_j=eN_k*nDOF_trial_element+j;
870 //int eN_k_j_nSpace = eN_k_j*nSpace;
871 int j_nSpace = j*nSpace;
872 dpdeResidual_u_u[j]=ck.MassJacobian_strong(dm_t,u_trial_ref.data()[k*nDOF_trial_element+j]) +
873 ck.HamiltonianJacobian_strong(dH_strong,&u_grad_trial[j_nSpace]);
874 }
875 //tau and tau*Res
876 calculateSubgridError_tau(elementDiameter.data()[eN],
877 dm_t,
878 dH_tau,
879 q_cfl.data()[eN_k],
880 tau0);
881 calculateSubgridError_tau(G,
882 dH_tau,
883 tau1,
884 q_cfl.data()[eN_k]);
885
886 tau = useMetrics*tau1+(1.0-useMetrics)*tau0;
887
888 for (int j=0;j<nDOF_trial_element;j++)
889 dsubgridError_u_u[j] = -tau*dpdeResidual_u_u[j];
890
891 nu_sc = q_numDiff_u.data()[eN_k]*(1.0-lag_shockCapturingScale) + q_numDiff_u_last.data()[eN_k]*lag_shockCapturingScale;
892 /* double epsilon_background_diffusion = 3.0*epsFact_redist*h_phi; */
893 /* if (fabs(u0) > epsilon_background_diffusion) */
894 /* nu_sc += backgroundDiffusionFactor*h_phi; */
895 double epsilon_background_diffusion = 2.0*epsFact_redist*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
896 if (fabs(phi_ls.data()[eN_k]) > epsilon_background_diffusion)
897 nu_sc += backgroundDiffusionFactor*elementDiameter.data()[eN]; //
898 for(int i=0;i<nDOF_test_element;i++)
899 {
900 //int eN_k_i=eN_k*nDOF_test_element+i;
901 //int eN_k_i_nSpace=eN_k_i*nSpace;
902 circ += ck.Reaction_weak(gf.D(epsilon_redist, phi_ls.data()[eN_k]), u_test_dV[i]);
903 for(int j=0;j<nDOF_trial_element;j++)
904 {
905 //int eN_k_j=eN_k*nDOF_trial_element+j;
906 //int eN_k_j_nSpace = eN_k_j*nSpace;
907 int j_nSpace = j*nSpace;
908 int i_nSpace = i*nSpace;
909 double FREEZE=double(freezeLevelSet);
910 //std::cout<<element_phi[i]<<'\t'<<element_nodes[i*3+0]<<'\t'<<element_nodes[i*3+1]<<std::endl;
911 //std::cout<<"eN_k "<<eN_k<<" D-J "<<gf.D(epsilon_redist,phi_ls.data()[eN_k])<<std::endl;
912 elementJacobian_u_u[i][j] += ck.MassJacobian_weak(dm_t,u_trial_ref.data()[k*nDOF_trial_element+j],u_test_dV[i]) +
913 ck.HamiltonianJacobian_weak(dH,&u_grad_trial[j_nSpace],u_test_dV[i]) +
914 (1.0-FREEZE)*(weakDirichletFactor/h_phi)*ck.ReactionJacobian_weak(-gf.D(epsilon_redist,u0),
915 u_trial_ref.data()[k*nDOF_trial_element+j],
916 u_test_dV[i]) +
917 ck.SubgridErrorJacobian(dsubgridError_u_u[j],Lstar_u[i]) +
918 ck.NumericalDiffusionJacobian(nu_sc,&u_grad_trial[j_nSpace],&u_grad_test_dV[i_nSpace]);
919 }//j
920 }//i
921 }//k
922 //
923 //load into element Jacobian into global Jacobian
924 //
925
926 //now try to account for weak dirichlet conditions in interior (frozen level set values)
927 if (freezeLevelSet)
928 {
929 //assume correspondence between dof and equations
930 for (int j = 0; j < nDOF_trial_element; j++)
931 {
932 const int J = u_l2g.data()[eN*nDOF_trial_element+j];
933 if (fabs(u_weak_internal_bc_dofs.data()[J]) < epsilon_redist)
934 //if (weakDirichletConditionFlags[J] == 1)
935 //if (fabs(gf.exact.phi_dof_corrected[j]) < epsilon_redist)
936 {
937 for (int jj=0; jj < nDOF_trial_element; jj++)
938 elementJacobian_u_u[j][jj] = 0.0;
939 elementJacobian_u_u[j][j] = weakDirichletFactor*elementDiameter.data()[eN];
940 }
941 }
942 }
943 for (int i=0;i<nDOF_test_element;i++)
944 {
945 int eN_i = eN*nDOF_test_element+i;
946 for (int j=0;j<nDOF_trial_element;j++)
947 {
948 int eN_i_j = eN_i*nDOF_trial_element+j;
949#ifdef CKDEBUG
950 std::cout<<"element jacobian i = "<<i<<"\t"<<"j = "<<j<<"\t"<<elementJacobian_u_u[i][j]<<std::endl;
951#endif
952 globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_u_u.data()[eN_i_j]] += elementJacobian_u_u[i][j];
953 }//j
954 }//i
955 }//elements
956 /* //cek todo should get rid of this, see res */
957 /* // */
958 /* //loop over exterior element boundaries to compute the surface integrals and load them into the global Jacobian */
959 /* // */
960 /* gf.useExact=false;//exact Heaviside integration not implemented for boundaries yet */
961 /* for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++) */
962 /* { */
963 /* int ebN = exteriorElementBoundariesArray.data()[ebNE]; */
964 /* int eN = elementBoundaryElementsArray.data()[ebN*2+0], */
965 /* ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0], */
966 /* eN_nDOF_trial_element = eN*nDOF_trial_element; */
967 /* double epsilon_redist,h_phi; */
968 /* for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++) */
969 /* { */
970 /* int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb; */
971 /* //int ebNE_kb_nSpace = ebNE_kb*nSpace; */
972 /* int ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb; */
973 /* int ebN_local_kb_nSpace = ebN_local_kb*nSpace; */
974
975 /* double u_ext=0.0, */
976 /* grad_u_ext[nSpace], */
977 /* m_ext=0.0, */
978 /* dm_ext=0.0, */
979 /* H_ext=0.0, */
980 /* dH_ext[nSpace], */
981 /* r_ext=0.0, */
982 /* // flux_ext=0.0, */
983 /* dflux_u_u_ext=0.0, */
984 /* bc_u_ext=0.0, */
985 /* bc_grad_u_ext[nSpace], */
986 /* bc_m_ext=0.0, */
987 /* bc_dm_ext=0.0, */
988 /* bc_H_ext=0.0, */
989 /* bc_dH_ext[nSpace], */
990 /* bc_r_ext=0.0, */
991 /* fluxJacobian_u_u[nDOF_trial_element], */
992 /* jac_ext[nSpace*nSpace], */
993 /* jacDet_ext, */
994 /* jacInv_ext[nSpace*nSpace], */
995 /* boundaryJac[nSpace*(nSpace-1)], */
996 /* metricTensor[(nSpace-1)*(nSpace-1)], */
997 /* metricTensorDetSqrt, */
998 /* dS, */
999 /* u_test_dS[nDOF_test_element], */
1000 /* u_grad_trial_trace[nDOF_trial_element*nSpace], */
1001 /* normal[nSpace],x_ext,y_ext,z_ext, */
1002 /* G[nSpace*nSpace],G_dd_G,tr_G, dir[nSpace],norm; */
1003 /* // */
1004 /* //calculate the solution and gradients at quadrature points */
1005 /* // */
1006 /* ck.calculateMapping_elementBoundary(eN, */
1007 /* ebN_local, */
1008 /* kb, */
1009 /* ebN_local_kb, */
1010 /* mesh_dof.data(), */
1011 /* mesh_l2g.data(), */
1012 /* mesh_trial_trace_ref.data(), */
1013 /* mesh_grad_trial_trace_ref.data(), */
1014 /* boundaryJac_ref.data(), */
1015 /* jac_ext, */
1016 /* jacDet_ext, */
1017 /* jacInv_ext, */
1018 /* boundaryJac, */
1019 /* metricTensor, */
1020 /* metricTensorDetSqrt, */
1021 /* normal_ref.data(), */
1022 /* normal, */
1023 /* x_ext,y_ext,z_ext); */
1024 /* dS = metricTensorDetSqrt*dS_ref.data()[kb]; */
1025 /* ck.calculateG(jacInv_ext,G,G_dd_G,tr_G); */
1026 /* //compute shape and solution information */
1027 /* //shape */
1028 /* ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,u_grad_trial_trace); */
1029 /* //solution and gradients */
1030 /* 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); */
1031 /* ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial_trace,grad_u_ext); */
1032 /* //precalculate test function products with integration weights */
1033 /* for (int j=0;j<nDOF_trial_element;j++) */
1034 /* { */
1035 /* u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS; */
1036 /* } */
1037
1038 /* norm = 1.0e-8; */
1039 /* for (int I=0;I<nSpace;I++) */
1040 /* norm += grad_u_ext[I]*grad_u_ext[I]; */
1041 /* norm = sqrt(norm); */
1042 /* for(int I=0;I<nSpace;I++) */
1043 /* dir[I] = grad_u_ext[I]/norm; */
1044
1045 /* ck.calculateGScale(G,dir,h_phi); */
1046 /* epsilon_redist = epsFact_redist*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]); */
1047
1048 /* // */
1049 /* //load the boundary values */
1050 /* // */
1051 /* bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+(1-isDOFBoundary_u.data()[ebNE_kb])*u_ext; */
1052 /* // */
1053 /* //calculate the internal and external trace of the pde coefficients */
1054 /* // */
1055 /* evaluateCoefficients(epsilon_redist, */
1056 /* ebqe_phi_ls_ext.data()[ebNE_kb], */
1057 /* u_ext, */
1058 /* grad_u_ext, */
1059 /* m_ext, */
1060 /* dm_ext, */
1061 /* H_ext, */
1062 /* dH_ext, */
1063 /* r_ext); */
1064 /* evaluateCoefficients(epsilon_redist, */
1065 /* ebqe_phi_ls_ext.data()[ebNE_kb], */
1066 /* bc_u_ext, */
1067 /* bc_grad_u_ext, */
1068 /* bc_m_ext, */
1069 /* bc_dm_ext, */
1070 /* bc_H_ext, */
1071 /* bc_dH_ext, */
1072 /* bc_r_ext); */
1073 /* // */
1074 /* //calculate the numerical fluxes */
1075 /* // */
1076 /* //DoNothing for now */
1077 /* // */
1078 /* //calculate the flux jacobian */
1079 /* // */
1080 /* for (int j=0;j<nDOF_trial_element;j++) */
1081 /* { */
1082 /* //int ebNE_kb_j = ebNE_kb*nDOF_trial_element+j; */
1083 /* //int ebNE_kb_j_nSpace = ebNE_kb_j*nSpace; */
1084 /* //int j_nSpace = j*nSpace; */
1085 /* int ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j; */
1086
1087 /* fluxJacobian_u_u[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_u_u_ext,u_trial_trace_ref.data()[ebN_local_kb_j]); */
1088 /* }//j */
1089 /* // */
1090 /* //update the global Jacobian from the flux Jacobian */
1091 /* // */
1092 /* for (int i=0;i<nDOF_test_element;i++) */
1093 /* { */
1094 /* int eN_i = eN*nDOF_test_element+i; */
1095 /* //int ebNE_kb_i = ebNE_kb*nDOF_test_element+i; */
1096 /* for (int j=0;j<nDOF_trial_element;j++) */
1097 /* { */
1098 /* int ebN_i_j = ebN*4*nDOF_test_X_trial_element + i*nDOF_trial_element + j; */
1099 /* //mwf debug */
1100 /* assert(fluxJacobian_u_u[j] == 0.0); */
1101 /* globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j]] += fluxJacobian_u_u[j]*u_test_dS[i]; */
1102 /* }//j */
1103 /* }//i */
1104 /* }//kb */
1105 /* }//ebNE */
1106 /* gf.useExact = useExact;//just to be safe */
1107 }//computeJacobian
1108
1109 void calculateResidual_ellipticRedist(arguments_dict& args,
1110 bool useExact)
1111 {
1112 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1113 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1114 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1115 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1116 xt::pyarray<double>& x_ref = args.array<double>("x_ref");
1117 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1118 xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1119 xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1120 xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1121 xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
1122 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
1123 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
1124 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
1125 xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
1126 xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
1127 xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
1128 xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
1129 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
1130 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
1131 int nElements_global = args.scalar<int>("nElements_global");
1132 double useMetrics = args.scalar<double>("useMetrics");
1133 double alphaBDF = args.scalar<double>("alphaBDF");
1134 double epsFact_redist = args.scalar<double>("epsFact_redist");
1135 double backgroundDiffusionFactor = args.scalar<double>("backgroundDiffusionFactor");
1136 double weakDirichletFactor = args.scalar<double>("weakDirichletFactor");
1137 int freezeLevelSet = args.scalar<int>("freezeLevelSet");
1138 int useTimeIntegration = args.scalar<int>("useTimeIntegration");
1139 int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
1140 int lag_subgridError = args.scalar<int>("lag_subgridError");
1141 double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
1142 xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1143 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1144 xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
1145 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1146 xt::pyarray<double>& phi_dof = args.array<double>("phi_dof");
1147 xt::pyarray<double>& phi_ls = args.array<double>("phi_ls");
1148 xt::pyarray<double>& q_m = args.array<double>("q_m");
1149 xt::pyarray<double>& q_u = args.array<double>("q_u");
1150 xt::pyarray<double>& q_n = args.array<double>("q_n");
1151 xt::pyarray<double>& q_dH = args.array<double>("q_dH");
1152 xt::pyarray<double>& u_weak_internal_bc_dofs = args.array<double>("u_weak_internal_bc_dofs");
1153 xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
1154 xt::pyarray<double>& q_dH_last = args.array<double>("q_dH_last");
1155 xt::pyarray<double>& q_cfl = args.array<double>("q_cfl");
1156 xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
1157 xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
1158 xt::pyarray<int>& weakDirichletConditionFlags = args.array<int>("weakDirichletConditionFlags");
1159 int offset_u = args.scalar<int>("offset_u");
1160 int stride_u = args.scalar<int>("stride_u");
1161 xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
1162 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
1163 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
1164 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
1165 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
1166 xt::pyarray<double>& ebqe_phi_ls_ext = args.array<double>("ebqe_phi_ls_ext");
1167 xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
1168 xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
1169 xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
1170 xt::pyarray<double>& ebqe_n = args.array<double>("ebqe_n");
1171 int ELLIPTIC_REDISTANCING = args.scalar<int>("ELLIPTIC_REDISTANCING");
1172 double backgroundDissipationEllipticRedist = args.scalar<double>("backgroundDissipationEllipticRedist");
1173 xt::pyarray<double>& lumped_qx = args.array<double>("lumped_qx");
1174 xt::pyarray<double>& lumped_qy = args.array<double>("lumped_qy");
1175 xt::pyarray<double>& lumped_qz = args.array<double>("lumped_qz");
1176 double alpha = args.scalar<double>("alpha");
1177 gf.useExact=useExact;
1178 //
1179 //loop over elements to compute volume integrals and load them into element and global residual
1180 //
1181 //eN is the element index
1182 //eN_k is the quadrature point index for a scalar
1183 //eN_k_nSpace is the quadrature point index for a vector
1184 //eN_i is the element test function index
1185 //eN_j is the element trial function index
1186 //eN_k_j is the quadrature point index for a trial function
1187 //eN_k_i is the quadrature point index for a trial function
1188 for(int eN=0;eN<nElements_global;eN++)
1189 {
1190 //declare local storage for element residual and initialize
1191 double elementResidual_u[nDOF_test_element],element_phi[nDOF_trial_element];
1192 double epsilon_redist,h_phi, norm;
1193 for (int i=0;i<nDOF_test_element;i++)
1194 {
1195 int eN_i=eN*nDOF_trial_element+i;
1196 elementResidual_u[i]=0.0;
1197 element_phi[i] = phi_dof.data()[u_l2g.data()[eN_i]];
1198 }//i
1199 double element_nodes[nDOF_mesh_trial_element*3];
1200 for (int i=0;i<nDOF_mesh_trial_element;i++)
1201 {
1202 int eN_i=eN*nDOF_mesh_trial_element+i;
1203 for(int I=0;I<3;I++)
1204 element_nodes[i*3 + I] = mesh_dof.data()[mesh_l2g.data()[eN_i]*3 + I];
1205 }//i
1206 gf.calculate(element_phi, element_nodes, x_ref.data(),false);
1207 //loop over quadrature points and compute integrands
1208 for (int k=0;k<nQuadraturePoints_element;k++)
1209 {
1210 gf.set_quad(k);
1211 //compute indeces and declare local storage
1212 int eN_k = eN*nQuadraturePoints_element+k,
1213 eN_k_nSpace = eN_k*nSpace,
1214 eN_nDOF_trial_element = eN*nDOF_trial_element;
1215 double
1216 coeff, delta,
1217 qx, qy, qz, normalReconstruction[nSpace],
1218 u=0,grad_u[nSpace],
1219 m=0.0,
1220 jac[nSpace*nSpace], jacDet, jacInv[nSpace*nSpace],
1221 u_grad_trial[nDOF_trial_element*nSpace],
1222 u_test_dV[nDOF_trial_element], u_grad_test_dV[nDOF_test_element*nSpace],
1223 dV,x,y,z,G[nSpace*nSpace],G_dd_G,tr_G;
1224 ck.calculateMapping_element(eN,
1225 k,
1226 mesh_dof.data(),
1227 mesh_l2g.data(),
1228 mesh_trial_ref.data(),
1229 mesh_grad_trial_ref.data(),
1230 jac,
1231 jacDet,
1232 jacInv,
1233 x,y,z);
1234 ck.calculateH_element(eN,
1235 k,
1236 nodeDiametersArray.data(),
1237 mesh_l2g.data(),
1238 mesh_trial_ref.data(),
1239 h_phi);
1240 //get the physical integration weight
1241 dV = fabs(jacDet)*dV_ref.data()[k];
1242 ck.calculateG(jacInv,G,G_dd_G,tr_G);
1243 //get the trial function gradients
1244 ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],
1245 jacInv,
1246 u_grad_trial);
1247 //get the solution
1248 ck.valFromDOF(u_dof.data(),
1249 &u_l2g.data()[eN_nDOF_trial_element],
1250 &u_trial_ref.data()[k*nDOF_trial_element],
1251 u);
1252 //get the solution gradients
1253 ck.gradFromDOF(u_dof.data(),
1254 &u_l2g.data()[eN_nDOF_trial_element],
1255 u_grad_trial,
1256 grad_u);
1257 if (ELLIPTIC_REDISTANCING > 1)
1258 { // use linear elliptic re-distancing via C0 normal reconstruction
1259 ck.valFromDOF(lumped_qx.data(),
1260 &u_l2g.data()[eN_nDOF_trial_element],
1261 &u_trial_ref.data()[k*nDOF_trial_element],
1262 qx);
1263 ck.valFromDOF(lumped_qy.data(),
1264 &u_l2g.data()[eN_nDOF_trial_element],
1265 &u_trial_ref.data()[k*nDOF_trial_element],
1266 qy);
1267 ck.valFromDOF(lumped_qz.data(),
1268 &u_l2g.data()[eN_nDOF_trial_element],
1269 &u_trial_ref.data()[k*nDOF_trial_element],
1270 qz);
1271 normalReconstruction[0] = qx;
1272 normalReconstruction[1] = qy;
1273 if (nSpace == 3)
1274 normalReconstruction[2] = qz;
1275 }
1276 //precalculate test function products with integration weights
1277 for (int j=0;j<nDOF_trial_element;j++)
1278 {
1279 u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
1280 for (int I=0;I<nSpace;I++)
1281 u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;
1282 }
1283 // MOVING MESH. Omit for now //
1284 // COMPUTE NORM OF GRAD(u) //
1285 double norm_grad_u = 0.;
1286 for (int I=0;I<nSpace;I++)
1287 norm_grad_u += grad_u[I]*grad_u[I];
1288 norm_grad_u = std::sqrt(norm_grad_u) + 1.0E-10;
1289
1290 // SAVE MASS AND SOLUTION FOR OTHER MODELS //
1291 q_m.data()[eN_k] = u; //m=u
1292 q_u.data()[eN_k] = u;
1293 for (int I=0;I<nSpace;I++)
1294 q_n.data()[eN_k_nSpace+I] = grad_u[I]/norm_grad_u;
1295
1296 // COMPUTE COEFFICIENTS //
1297 if (SINGLE_POTENTIAL == 1)
1298 coeff = 1.0-1.0/norm_grad_u; //single potential
1299 else // double potential
1300 coeff = 1.0+2*std::pow(norm_grad_u,2)-3*norm_grad_u;
1301
1302 // COMPUTE DELTA FUNCTION //
1303 epsilon_redist = epsFact_redist*(useMetrics*h_phi
1304 +(1.0-useMetrics)*elementDiameter.data()[eN]);
1305 delta = gf.D(epsilon_redist,phi_ls.data()[eN_k]);
1306
1307 // COMPUTE STRONG RESIDUAL //
1308 double Si = -1.0+2.0*gf.H(epsilon_redist,phi_ls.data()[eN_k]);
1309 double residualEikonal = Si*(norm_grad_u-1.0);
1310 double backgroundDissipation = backgroundDissipationEllipticRedist*elementDiameter.data()[eN];
1311
1312 // UPDATE ELEMENT RESIDUAL //
1313 for(int i=0;i<nDOF_test_element;i++)
1314 {
1315 int i_nSpace = i*nSpace;
1316 // global i-th index
1317 int gi = offset_u+stride_u*u_l2g.data()[eN*nDOF_test_element+i];
1318
1319 if (ELLIPTIC_REDISTANCING > 1) // (NON)LINEAR VIA C0 NORMAL RECONSTRUCTION
1320 {
1321 elementResidual_u[i] +=
1322 residualEikonal*u_test_dV[i]
1323 +ck.NumericalDiffusion(1.0+backgroundDissipation,
1324 grad_u,
1325 &u_grad_test_dV[i_nSpace])
1326 -ck.NumericalDiffusion(1.0,
1327 normalReconstruction,
1328 &u_grad_test_dV[i_nSpace])
1329 +alpha*(u_dof.data()[gi]-phi_dof.data()[gi])*delta*u_test_dV[i]; // BCs
1330 }
1331 else // =1. Nonlinear via single or double pot.
1332 {
1333 elementResidual_u[i] +=
1334 residualEikonal*u_test_dV[i]
1335 +ck.NumericalDiffusion(coeff+backgroundDissipation,
1336 grad_u,
1337 &u_grad_test_dV[i_nSpace])
1338 + alpha*(u_dof.data()[gi]-phi_dof.data()[gi])*delta*u_test_dV[i]; // BCs
1339 }
1340 }//i
1341 }//k
1342 //
1343 //load element into global residual and save element residual
1344 //
1345 for(int i=0;i<nDOF_test_element;i++)
1346 {
1347 int eN_i=eN*nDOF_test_element+i;
1348 globalResidual.data()[offset_u+stride_u*u_l2g.data()[eN_i]]+=elementResidual_u[i];
1349 }//i
1350 }//elements
1351 //
1352 //loop over exterior element boundaries to save soln at quad points
1353 //
1354 for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
1355 {
1356 int ebN = exteriorElementBoundariesArray.data()[ebNE],
1357 eN = elementBoundaryElementsArray.data()[ebN*2+0],
1358 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
1359 eN_nDOF_trial_element = eN*nDOF_trial_element;
1360 for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
1361 {
1362 int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
1363 ebNE_kb_nSpace = ebNE_kb*nSpace,
1364 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
1365 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
1366 double
1367 u_ext=0.0,
1368 grad_u_ext[nSpace],
1369 jac_ext[nSpace*nSpace],jacDet_ext,jacInv_ext[nSpace*nSpace],
1370 boundaryJac[nSpace*(nSpace-1)],
1371 metricTensor[(nSpace-1)*(nSpace-1)],metricTensorDetSqrt,
1372 u_grad_trial_trace[nDOF_trial_element*nSpace],
1373 normal[nSpace],x_ext,y_ext,z_ext,
1374 dir[nSpace],norm;
1375 ck.calculateMapping_elementBoundary(eN,
1376 ebN_local,
1377 kb,
1378 ebN_local_kb,
1379 mesh_dof.data(),
1380 mesh_l2g.data(),
1381 mesh_trial_trace_ref.data(),
1382 mesh_grad_trial_trace_ref.data(),
1383 boundaryJac_ref.data(),
1384 jac_ext,
1385 jacDet_ext,
1386 jacInv_ext,
1387 boundaryJac,
1388 metricTensor,
1389 metricTensorDetSqrt,
1390 normal_ref.data(),
1391 normal,
1392 x_ext,y_ext,z_ext);
1393 //compute shape and solution information
1394 //shape
1395 ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],
1396 jacInv_ext,
1397 u_grad_trial_trace);
1398 //solution and gradients
1399 ck.valFromDOF(u_dof.data(),
1400 &u_l2g.data()[eN_nDOF_trial_element],
1401 &u_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],
1402 u_ext);
1403 ck.gradFromDOF(u_dof.data(),
1404 &u_l2g.data()[eN_nDOF_trial_element],
1405 u_grad_trial_trace,
1406 grad_u_ext);
1407 norm = 0;
1408 for (int I=0;I<nSpace;I++)
1409 norm += grad_u_ext[I]*grad_u_ext[I];
1410 norm = sqrt(norm) + 1.0E-10;
1411 for (int I=0;I<nSpace;I++)
1412 dir[I] = grad_u_ext[I]/norm;
1413
1414 //save for other models
1415 ebqe_u.data()[ebNE_kb] = u_ext;
1416
1417 for (int I=0;I<nSpace;I++)
1418 ebqe_n.data()[ebNE_kb_nSpace+I] = dir[I];
1419 }//kb
1420 }//ebNE
1421 }
1422
1423 void calculateJacobian_ellipticRedist(arguments_dict& args,
1424 bool useExact)
1425 {
1426 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1427 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1428 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1429 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1430 xt::pyarray<double>& x_ref = args.array<double>("x_ref");
1431 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1432 xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1433 xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1434 xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1435 xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
1436 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
1437 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
1438 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
1439 xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
1440 xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
1441 xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
1442 xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
1443 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
1444 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
1445 int nElements_global = args.scalar<int>("nElements_global");
1446 double useMetrics = args.scalar<double>("useMetrics");
1447 double alphaBDF = args.scalar<double>("alphaBDF");
1448 double epsFact_redist = args.scalar<double>("epsFact_redist");
1449 double backgroundDiffusionFactor = args.scalar<double>("backgroundDiffusionFactor");
1450 double weakDirichletFactor = args.scalar<double>("weakDirichletFactor");
1451 int freezeLevelSet = args.scalar<int>("freezeLevelSet");
1452 int useTimeIntegration = args.scalar<int>("useTimeIntegration");
1453 int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
1454 int lag_subgridError = args.scalar<int>("lag_subgridError");
1455 double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
1456 xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1457 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1458 xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
1459 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1460 xt::pyarray<double>& phi_dof = args.array<double>("phi_dof");
1461 xt::pyarray<double>& u_weak_internal_bc_dofs = args.array<double>("u_weak_internal_bc_dofs");
1462 xt::pyarray<double>& phi_ls = args.array<double>("phi_ls");
1463 xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
1464 xt::pyarray<double>& q_dH_last = args.array<double>("q_dH_last");
1465 xt::pyarray<double>& q_cfl = args.array<double>("q_cfl");
1466 xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
1467 xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
1468 xt::pyarray<int>& weakDirichletConditionFlags = args.array<int>("weakDirichletConditionFlags");
1469 xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
1470 xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
1471 xt::pyarray<double>& globalJacobian = args.array<double>("globalJacobian");
1472 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
1473 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
1474 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
1475 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
1476 xt::pyarray<double>& ebqe_phi_ls_ext = args.array<double>("ebqe_phi_ls_ext");
1477 xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
1478 xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
1479 xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
1480 int ELLIPTIC_REDISTANCING = args.scalar<int>("ELLIPTIC_REDISTANCING");
1481 double backgroundDissipationEllipticRedist = args.scalar<double>("backgroundDissipationEllipticRedist");
1482 double alpha = args.scalar<double>("alpha");
1483 gf.useExact=useExact;
1484 //
1485 //loop over elements
1486 //
1487 for(int eN=0;eN<nElements_global;eN++)
1488 {
1489 double elementJacobian_u_u[nDOF_test_element][nDOF_trial_element],element_phi[nDOF_trial_element];
1490 double epsilon_redist,h_phi, norm;
1491 for (int i=0;i<nDOF_test_element;i++)
1492 {
1493 int eN_i=eN*nDOF_trial_element+i;
1494 element_phi[i] = phi_dof.data()[u_l2g.data()[eN_i]];
1495 for (int j=0;j<nDOF_trial_element;j++)
1496 {
1497 elementJacobian_u_u[i][j]=0.0;
1498 }
1499 }
1500 double element_nodes[nDOF_mesh_trial_element*3];
1501 for (int i=0;i<nDOF_mesh_trial_element;i++)
1502 {
1503 int eN_i=eN*nDOF_mesh_trial_element+i;
1504 for(int I=0;I<3;I++)
1505 element_nodes[i*3 + I] = mesh_dof.data()[mesh_l2g.data()[eN_i]*3 + I];
1506 }//i
1507 gf.calculate(element_phi, element_nodes, x_ref.data(),false);
1508 for (int k=0;k<nQuadraturePoints_element;k++)
1509 {
1510 gf.set_quad(k);
1511 int eN_k = eN*nQuadraturePoints_element+k, //index to a scalar at a quadrature point
1512 eN_k_nSpace = eN_k*nSpace,
1513 eN_nDOF_trial_element = eN*nDOF_trial_element; //index to a vector at a quadrature point
1514
1515 //declare local storage
1516 double
1517 coeff1, coeff2, delta,
1518 grad_u[nSpace],
1519 jac[nSpace*nSpace], jacDet, jacInv[nSpace*nSpace],
1520 u_grad_trial[nDOF_trial_element*nSpace],
1521 dV, u_test_dV[nDOF_test_element], u_grad_test_dV[nDOF_test_element*nSpace],
1522 x,y,z,G[nSpace*nSpace],G_dd_G,tr_G;
1523 //
1524 //calculate solution and gradients at quadrature points
1525 //
1526 ck.calculateMapping_element(eN,
1527 k,
1528 mesh_dof.data(),
1529 mesh_l2g.data(),
1530 mesh_trial_ref.data(),
1531 mesh_grad_trial_ref.data(),
1532 jac,
1533 jacDet,
1534 jacInv,
1535 x,y,z);
1536 ck.calculateH_element(eN,
1537 k,
1538 nodeDiametersArray.data(),
1539 mesh_l2g.data(),
1540 mesh_trial_ref.data(),
1541 h_phi);
1542 //get the physical integration weight
1543 dV = fabs(jacDet)*dV_ref.data()[k];
1544 ck.calculateG(jacInv,G,G_dd_G,tr_G);
1545 //get the trial function gradients
1546 ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],
1547 jacInv,
1548 u_grad_trial);
1549 //get the solution gradients
1550 ck.gradFromDOF(u_dof.data(),
1551 &u_l2g.data()[eN_nDOF_trial_element],
1552 u_grad_trial,
1553 grad_u);
1554 //precalculate test function products with integration weights
1555 for (int j=0;j<nDOF_trial_element;j++)
1556 {
1557 u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
1558 for (int I=0;I<nSpace;I++)
1559 u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;
1560 }
1561 // MOVING MESH. Omit for now //
1562 // COMPUTE NORM OF GRAD(u) //
1563 double norm_grad_u = 0;
1564 for(int I=0;I<nSpace;I++)
1565 norm_grad_u += grad_u[I]*grad_u[I];
1566 norm_grad_u = std::sqrt(norm_grad_u) + 1.0E-10;
1567
1568 // COMPUTE COEFFICIENTS //
1569 if (SINGLE_POTENTIAL==1)
1570 {
1571 coeff1 = 0.; //-1./norm_grad_u;
1572 coeff2 = 1./std::pow(norm_grad_u,3);
1573 }
1574 else
1575 {
1576 coeff1 = fmax(1.0E-10, 2*std::pow(norm_grad_u,2)-3*norm_grad_u);
1577 coeff2 = fmax(1.0E-10, 4.-3./norm_grad_u);
1578 }
1579
1580 // COMPUTE DELTA FUNCTION //
1581 epsilon_redist = epsFact_redist*(useMetrics*h_phi
1582 +(1.0-useMetrics)*elementDiameter.data()[eN]);
1583 delta = gf.D(epsilon_redist,phi_ls.data()[eN_k]);
1584
1585 // COMPUTE STRONG Jacobian //
1586 double Si = -1.0+2.0*gf.H(epsilon_redist,phi_ls.data()[eN_k]);
1587 double dH[nSpace];
1588 for (int I=0; I<nSpace;I++)
1589 dH[I] = Si*grad_u[I]/norm_grad_u;
1590 double backgroundDissipation = backgroundDissipationEllipticRedist*elementDiameter.data()[eN];
1591
1592 // LOOP IN I-DOFs //
1593 for(int i=0;i<nDOF_test_element;i++)
1594 {
1595 int i_nSpace = i*nSpace;
1596 for(int j=0;j<nDOF_trial_element;j++)
1597 {
1598 int j_nSpace = j*nSpace;
1599 elementJacobian_u_u[i][j] +=
1600 ck.HamiltonianJacobian_weak(dH,&u_grad_trial[j_nSpace],u_test_dV[i])
1601 +ck.NumericalDiffusionJacobian(1.0+backgroundDissipation,
1602 &u_grad_trial[j_nSpace],
1603 &u_grad_test_dV[i_nSpace])
1604 + (ELLIPTIC_REDISTANCING == 1 ? 1. : 0.)*
1605 ( ck.NumericalDiffusionJacobian(coeff1,
1606 &u_grad_trial[j_nSpace],
1607 &u_grad_test_dV[i_nSpace])
1608 + coeff2*dV*
1609 ck.NumericalDiffusion(1.0,grad_u,&u_grad_trial[i_nSpace])*
1610 ck.NumericalDiffusion(1.0,grad_u,&u_grad_trial[j_nSpace]) )
1611 + (i == j ? alpha*delta*u_test_dV[i] : 0.); //lumped
1612 }//j
1613 }//i
1614 }//k
1615 //
1616 //load into element Jacobian into global Jacobian
1617 //
1618 for (int i=0;i<nDOF_test_element;i++)
1619 {
1620 int eN_i = eN*nDOF_test_element+i;
1621 for (int j=0;j<nDOF_trial_element;j++)
1622 {
1623 int eN_i_j = eN_i*nDOF_trial_element+j;
1624 globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i]
1625 + csrColumnOffsets_u_u.data()[eN_i_j]] += elementJacobian_u_u[i][j];
1626 }//j
1627 }//i
1628 }//elements
1629 }//computeJacobian
1630
1631 void normalReconstruction(arguments_dict& args)
1632 {
1633 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1634 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1635 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1636 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1637 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1638 xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1639 xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1640 xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1641 int nElements_global = args.scalar<int>("nElements_global");
1642 xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1643 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1644 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1645 int offset_u = args.scalar<int>("offset_u");
1646 int stride_u = args.scalar<int>("stride_u");
1647 int numDOFs = args.scalar<int>("numDOFs");
1648 xt::pyarray<double>& lumped_qx = args.array<double>("lumped_qx");
1649 xt::pyarray<double>& lumped_qy = args.array<double>("lumped_qy");
1650 xt::pyarray<double>& lumped_qz = args.array<double>("lumped_qz");
1651 weighted_lumped_mass_matrix.resize(numDOFs,0.0);
1652 for (int i=0; i<numDOFs; i++)
1653 {
1654 // output vectors
1655 lumped_qx.data()[i]=0.;
1656 lumped_qy.data()[i]=0.;
1657 lumped_qz.data()[i]=0.;
1658 // auxiliary vectors
1659 weighted_lumped_mass_matrix[i]=0.;
1660 }
1661 for(int eN=0;eN<nElements_global;eN++)
1662 {
1663 //declare local storage for local contributions and initialize
1664 double
1665 element_weighted_lumped_mass_matrix[nDOF_test_element],
1666 element_rhsx_normal_reconstruction[nDOF_test_element],
1667 element_rhsy_normal_reconstruction[nDOF_test_element],
1668 element_rhsz_normal_reconstruction[nDOF_test_element];
1669 for (int i=0;i<nDOF_test_element;i++)
1670 {
1671 element_weighted_lumped_mass_matrix[i]=0.0;
1672 element_rhsx_normal_reconstruction[i]=0.0;
1673 element_rhsy_normal_reconstruction[i]=0.0;
1674 element_rhsz_normal_reconstruction[i]=0.0;
1675 }
1676 //loop over quadrature points and compute integrands
1677 for (int k=0;k<nQuadraturePoints_element;k++)
1678 {
1679 //compute indeces and declare local storage
1680 int eN_k = eN*nQuadraturePoints_element+k,
1681 eN_k_nSpace = eN_k*nSpace,
1682 eN_nDOF_trial_element = eN*nDOF_trial_element;
1683 double
1684 //for mass matrix contributions
1685 grad_u[nSpace],
1686 u_grad_trial[nDOF_trial_element*nSpace],
1687 u_test_dV[nDOF_trial_element],
1688 //for general use
1689 jac[nSpace*nSpace], jacDet, jacInv[nSpace*nSpace],
1690 dV,x,y,z;
1691 //get the physical integration weight
1692 ck.calculateMapping_element(eN,
1693 k,
1694 mesh_dof.data(),
1695 mesh_l2g.data(),
1696 mesh_trial_ref.data(),
1697 mesh_grad_trial_ref.data(),
1698 jac,
1699 jacDet,
1700 jacInv,
1701 x,y,z);
1702 dV = fabs(jacDet)*dV_ref.data()[k];
1703 ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],
1704 jacInv,
1705 u_grad_trial);
1706 ck.gradFromDOF(u_dof.data(),
1707 &u_l2g.data()[eN_nDOF_trial_element],u_grad_trial,
1708 grad_u);
1709 //precalculate test function products with integration weights for mass matrix terms
1710 for (int j=0;j<nDOF_trial_element;j++)
1711 u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
1712
1713 double rhsx = grad_u[0];
1714 double rhsy = grad_u[1];
1715 double rhsz = 0;
1716 if (nSpace==3)
1717 rhsz = grad_u[2];
1718
1719 double norm_grad_u = 0;
1720 for (int I=0;I<nSpace; I++)
1721 norm_grad_u += grad_u[I]*grad_u[I];
1722 norm_grad_u = std::sqrt(norm_grad_u) + 1.0E-10;
1723
1724 for(int i=0;i<nDOF_test_element;i++)
1725 {
1726 element_weighted_lumped_mass_matrix[i] += norm_grad_u*u_test_dV[i];
1727 element_rhsx_normal_reconstruction[i] += rhsx*u_test_dV[i];
1728 element_rhsy_normal_reconstruction[i] += rhsy*u_test_dV[i];
1729 element_rhsz_normal_reconstruction[i] += rhsz*u_test_dV[i];
1730 }
1731 } //k
1732 // DISTRIBUTE //
1733 for(int i=0;i<nDOF_test_element;i++)
1734 {
1735 int eN_i=eN*nDOF_test_element+i;
1736 int gi = offset_u+stride_u*u_l2g.data()[eN_i]; //global i-th index
1737
1738 weighted_lumped_mass_matrix[gi] += element_weighted_lumped_mass_matrix[i];
1739 lumped_qx.data()[gi] += element_rhsx_normal_reconstruction[i];
1740 lumped_qy.data()[gi] += element_rhsy_normal_reconstruction[i];
1741 lumped_qz.data()[gi] += element_rhsz_normal_reconstruction[i];
1742 }//i
1743 }//elements
1744 // COMPUTE LUMPED L2 PROJECTION
1745 for (int i=0; i<numDOFs; i++)
1746 {
1747 // normal reconstruction
1748 double weighted_mi = weighted_lumped_mass_matrix[i];
1749 lumped_qx.data()[i] /= weighted_mi;
1750 lumped_qy.data()[i] /= weighted_mi;
1751 lumped_qz.data()[i] /= weighted_mi;
1752 }
1753 }
1754
1755 std::tuple<double, double, double> calculateMetricsAtEOS(arguments_dict& args)
1756 {
1757 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1758 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1759 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1760 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1761 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1762 xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1763 xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1764 xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1765 int nElements_global = args.scalar<int>("nElements_global");
1766 xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1767 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1768 double degree_polynomial = args.scalar<double>("degree_polynomial");
1769 double epsFact_redist = args.scalar<double>("epsFact_redist");
1770 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1771 xt::pyarray<double>& u_exact = args.array<double>("u_exact");
1772 int offset_u = args.scalar<int>("offset_u");
1773 int stride_u = args.scalar<int>("stride_u)");
1774 double global_V = 0.;
1775 double global_V0 = 0.;
1776 double global_I_err = 0.0;
1777 double global_V_err = 0.0;
1778 double global_D_err = 0.0;
1780 // ** LOOP IN CELLS //
1782 for(int eN=0;eN<nElements_global;eN++)
1783 {
1784 //declare local storage for local contributions and initialize
1785 double
1786 elementResidual_u[nDOF_test_element];
1787 double
1788 cell_mass_error = 0., cell_mass_exact = 0.,
1789 cell_I_err = 0.,
1790 cell_V = 0., cell_V0 = 0.,
1791 cell_D_err = 0.;
1792
1793 //loop over quadrature points and compute integrands
1794 for (int k=0;k<nQuadraturePoints_element;k++)
1795 {
1796 //compute indeces and declare local storage
1797 int eN_k = eN*nQuadraturePoints_element+k,
1798 eN_k_nSpace = eN_k*nSpace,
1799 eN_nDOF_trial_element = eN*nDOF_trial_element;
1800 double
1801 u, uh,
1802 u_grad_trial[nDOF_trial_element*nSpace],
1803 grad_uh[nSpace],
1804 //for general use
1805 jac[nSpace*nSpace], jacDet, jacInv[nSpace*nSpace],
1806 dV,x,y,z;
1807 //get the physical integration weight
1808 ck.calculateMapping_element(eN,
1809 k,
1810 mesh_dof.data(),
1811 mesh_l2g.data(),
1812 mesh_trial_ref.data(),
1813 mesh_grad_trial_ref.data(),
1814 jac,
1815 jacDet,
1816 jacInv,
1817 x,y,z);
1818 dV = fabs(jacDet)*dV_ref.data()[k];
1819 // get functions at quad points
1820 ck.valFromDOF(u_dof.data(),
1821 &u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],
1822 uh);
1823 u = u_exact.data()[eN_k];
1824 // get gradients
1825 ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],
1826 jacInv,
1827 u_grad_trial);
1828 ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial,grad_uh);
1829
1830 double epsHeaviside = epsFact_redist*elementDiameter.data()[eN]/degree_polynomial;
1831 // compute (smoothed) heaviside functions //
1832 double Hu = heaviside(u);
1833 double Huh = heaviside(uh);
1834 // compute cell metrics //
1835 cell_I_err += fabs(Hu - Huh)*dV;
1836 cell_V += Huh*dV;
1837 cell_V0 += Hu*dV;
1838
1839 double norm2_grad_uh = 0.;
1840 for (int I=0; I<nSpace; I++)
1841 norm2_grad_uh += grad_uh[I]*grad_uh[I];
1842 cell_D_err += std::pow(std::sqrt(norm2_grad_uh) - 1, 2.)*dV;
1843 }
1844 global_V += cell_V;
1845 global_V0 += cell_V0;
1846 // metrics //
1847 global_I_err += cell_I_err;
1848 global_D_err += cell_D_err;
1849 }//elements
1850 global_V_err = fabs(global_V0 - global_V)/global_V0;
1851 global_D_err *= 0.5;
1852 return std::tuple<double, double, double>(global_I_err, global_V_err, global_D_err);
1853 }
1854
1855 };//RDLS
1856 inline RDLS_base* newRDLS(int nSpaceIn,
1857 int nQuadraturePoints_elementIn,
1858 int nDOF_mesh_trial_elementIn,
1859 int nDOF_trial_elementIn,
1860 int nDOF_test_elementIn,
1861 int nQuadraturePoints_elementBoundaryIn,
1862 int CompKernelFlag)
1863 {
1864 if (nSpaceIn == 2)
1865 return proteus::chooseAndAllocateDiscretization2D<RDLS_base,RDLS,CompKernel>(nSpaceIn,
1866 nQuadraturePoints_elementIn,
1867 nDOF_mesh_trial_elementIn,
1868 nDOF_trial_elementIn,
1869 nDOF_test_elementIn,
1870 nQuadraturePoints_elementBoundaryIn,
1871 CompKernelFlag);
1872 else
1873 return proteus::chooseAndAllocateDiscretization<RDLS_base,RDLS,CompKernel>(nSpaceIn,
1874 nQuadraturePoints_elementIn,
1875 nDOF_mesh_trial_elementIn,
1876 nDOF_trial_elementIn,
1877 nDOF_test_elementIn,
1878 nQuadraturePoints_elementBoundaryIn,
1879 CompKernelFlag);
1880 }
1881
1882}//proteus
1883
1884#endif
r
Double r
Definition Headers.h:83
H
Double H
Definition Headers.h:65
u
Double u
Definition Headers.h:89
z
Double * z
Definition Headers.h:49
coeff
Double * coeff
Definition Headers.h:42
SINGLE_POTENTIAL
#define SINGLE_POTENTIAL
Definition RDLS.h:14
equivalent_polynomials::GeneralizedFunctions_mix
Definition equivalent_polynomials.h:1285
equivalent_polynomials::GeneralizedFunctions_mix::useExact
bool useExact
Definition equivalent_polynomials.h:1289
equivalent_polynomials::GeneralizedFunctions_mix::set_quad
void set_quad(unsigned int q)
Definition equivalent_polynomials.h:1325
equivalent_polynomials::GeneralizedFunctions_mix::H
double H(double eps, double phi)
Definition equivalent_polynomials.h:1337
equivalent_polynomials::GeneralizedFunctions_mix::exact
Simplex< nSpace, nP_ifem, nP, nQ, nEBQ, useIfemBasis > exact
Definition equivalent_polynomials.h:1288
equivalent_polynomials::GeneralizedFunctions_mix::ImH
double ImH(double eps, double phi)
Definition equivalent_polynomials.h:1345
equivalent_polynomials::GeneralizedFunctions_mix::D
double D(double eps, double phi)
Definition equivalent_polynomials.h:1353
equivalent_polynomials::GeneralizedFunctions_mix::calculate
int calculate(const double *phi_dof, const double *phi_nodes, const double *xi_r, double ma, double mb, double jf, bool isBoundary, bool scale)
Definition equivalent_polynomials.h:1294
proteus::RDLS_base
Definition RDLS.h:35
proteus::RDLS_base::normalReconstruction
virtual void normalReconstruction(arguments_dict &args)=0
proteus::RDLS_base::calculateJacobian_ellipticRedist
virtual void calculateJacobian_ellipticRedist(arguments_dict &args, bool useExact)=0
proteus::RDLS_base::weighted_lumped_mass_matrix
std::valarray< double > weighted_lumped_mass_matrix
Definition RDLS.h:37
proteus::RDLS_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args, bool useExact)=0
proteus::RDLS_base::calculateMetricsAtEOS
virtual std::tuple< double, double, double > calculateMetricsAtEOS(arguments_dict &args)=0
proteus::RDLS_base::calculateResidual_ellipticRedist
virtual void calculateResidual_ellipticRedist(arguments_dict &args, bool useExact)=0
proteus::RDLS_base::calculateResidual
virtual void calculateResidual(arguments_dict &args, bool useExact)=0
proteus::RDLS_base::~RDLS_base
virtual ~RDLS_base()
Definition RDLS.h:38
proteus::RDLS::evaluateCoefficients
void evaluateCoefficients(const double &eps, const double &u_levelSet, const double &u, const double grad_u[nSpace], double &m, double &dm, double &H, double dH[nSpace], double &r)
Definition RDLS.h:66
proteus::RDLS::calculateMetricsAtEOS
std::tuple< double, double, double > calculateMetricsAtEOS(arguments_dict &args)
Definition RDLS.h:1755
proteus::RDLS::calculateResidual_ellipticRedist
void calculateResidual_ellipticRedist(arguments_dict &args, bool useExact)
Definition RDLS.h:1109
proteus::RDLS::normalReconstruction
void normalReconstruction(arguments_dict &args)
Definition RDLS.h:1631
proteus::RDLS::calculateSubgridError_tau
void calculateSubgridError_tau(const double G[nSpace *nSpace], const double Ai[nSpace], double &tau_v, double &q_cfl)
Definition RDLS.h:123
proteus::RDLS::ck
CompKernelType ck
Definition RDLS.h:58
proteus::RDLS::gfu
GeneralizedFunctions< nSpace, 2, nQuadraturePoints_element, nQuadraturePoints_elementBoundary > gfu
Definition RDLS.h:59
proteus::RDLS::calculateJacobian_ellipticRedist
void calculateJacobian_ellipticRedist(arguments_dict &args, bool useExact)
Definition RDLS.h:1423
proteus::RDLS::calculateResidual
void calculateResidual(arguments_dict &args, bool useExact)
Definition RDLS.h:137
proteus::RDLS::nDOF_test_X_trial_element
const int nDOF_test_X_trial_element
Definition RDLS.h:57
proteus::RDLS::calculateSubgridError_tau
void calculateSubgridError_tau(const double &elementDiameter, const double &dmt, const double dH[nSpace], double &cfl, double &tau)
Definition RDLS.h:102
proteus::RDLS::gf
GeneralizedFunctions< nSpace, 2, nQuadraturePoints_element, nQuadraturePoints_elementBoundary > gf
Definition RDLS.h:59
proteus::RDLS::calculateJacobian
void calculateJacobian(arguments_dict &args, bool useExact)
Definition RDLS.h:616
equivalent_polynomials.h
proteus
Definition ADR.h:19
proteus::GeneralizedFunctions
equivalent_polynomials::GeneralizedFunctions_mix< nSpace, nP_ifem, nP, nQ, nEBQ, true > GeneralizedFunctions
Definition ADR.h:21
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::heaviside
double heaviside(const double &z)
Definition CLSVOF.h:20
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::newRDLS
RDLS_base * newRDLS(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition RDLS.h:1856
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