proteus 1.9.0
C/C++/Fortran libraries
Loading...
Searching...
No Matches
TADR.h
Go to the documentation of this file.
1#ifndef TADR_H
2#define TADR_H
3#include <cmath>
4#include <iomanip>
5#include <iostream>
6#include <valarray>
7#include "CompKernel.h"
8#include "ModelFactory.h"
9#include "ArgumentsDict.h"
10#include "xtensor-python/pyarray.hpp"
11#define nnz nSpace
12
13namespace py = pybind11;
14
15#define POWER_SMOOTHNESS_INDICATOR 2
16#define IS_BETAij_ONE 0
17
18// Cell based methods:
19// * Galerkin (unstabilized)
20// * VMS(SUPG) with BDF1 or BDF2 time integration
21// * Explicit Taylor Galerkin with EV stabilization
22// Edge based methods.
23// Low order via D. Kuzmin's
24// High order methods: Smoothness indicator with MC, EV commutator with MC, D.K with ML
25// Zalesak's FCT
26
27namespace proteus
28{
29 // ImplicitEV=5: backward-Euler implicit edge-based scheme (low-order graph
30 // dissipation, no FCT) modeled on Richards STABILIZATION_TYPE==2. Unlike the
31 // explicit edge-based types (2,3,4) the advection/diffusion are evaluated at
32 // the current Newton iterate and contribute to the Jacobian, so it is NOT
33 // CFL-limited. The Jacobian is assembled in calculateJacobian (a separate
34 // function), NOT inside calculateResidual as Richards does.
35 enum class STABILIZATION : int { Galerkin=-1, VMS=0, TaylorGalerkinEV=1, EntropyViscosity=2, SmoothnessIndicator=3, Kuzmin=4, ImplicitEV=5};
36 enum class ENTROPY : int { POWER=0, LOG=1};
37 enum class DISPERSION : int { Constant=0, PowerLawSaturation=1, VelocityBased=2};
38 // Power entropy //
39 inline double EPOWER(const double& phi, const double& phiL, const double& phiR)
40 {
41 return 1./2.*std::pow(fabs(phi),2.);
42 }
43 inline double DEPOWER(const double& phi, const double& phiL, const double& phiR)
44 {
45 return fabs(phi)*(phi>=0 ? 1 : -1);
46 }
47 // Log entropy // for level set from 0 to 1
48 inline double ELOG(const double& phi, const double& phiL, const double& phiR)
49 {
50 return std::log(fabs((phi-phiL)*(phiR-phi))+1E-14);
51 }
52 inline double DELOG(const double& phi, const double& phiL, const double& phiR)
53 {
54 return (phiL+phiR-2*phi)*((phi-phiL)*(phiR-phi)>=0 ? 1 : -1)/(fabs((phi-phiL)*(phiR-phi))+1E-14);
55 }
56}
57
58namespace proteus
59{
60 class TADR_base
61 {
62 //The base class defining the interface
63 public:
64 std::valarray<double> Rpos, Rneg;
65 std::valarray<double> FluxCorrectionMatrix;
66 std::valarray<double> TransportMatrix, DiffusionMatrix, TransposeTransportMatrix;
67 std::valarray<double> psi, eta, global_entropy_residual, boundary_integral;
68 std::valarray<double> m_dof, theta_dof_proj, rho_dof_proj, ML_mass_proj;
69 std::valarray<double> maxVel,maxEntRes;
70 virtual ~TADR_base(){}
71 virtual void calculateResidual(arguments_dict& args)=0;
72 virtual void calculateJacobian(arguments_dict& args)=0;
73 virtual void invert(arguments_dict& args)=0;
74 virtual void FCTStep(arguments_dict& args)=0;
75 };
76
77 template<class CompKernelType,
78 int nSpace,
79 int nQuadraturePoints_element,
80 int nDOF_mesh_trial_element,
81 int nDOF_trial_element,
82 int nDOF_test_element,
83 int nQuadraturePoints_elementBoundary>
84 class TADR : public TADR_base
85 {
86 public:
87 const int nDOF_test_X_trial_element;
88 CompKernelType ck;
89 TADR():
90 nDOF_test_X_trial_element(nDOF_test_element*nDOF_trial_element),
91 ck()
92 {}
93
94
95
96 inline
97 void calculateCFL(const double& elementDiameter,
98 const double df[nSpace],
99 double& cfl)
100 {
101 double h,nrm_v;
102 h = elementDiameter;
103 nrm_v=0.0;
104 for(int I=0;I<nSpace;I++)
105 nrm_v+=df[I]*df[I];
106 nrm_v = sqrt(nrm_v);
107 cfl = nrm_v/h;
108 }
109
110
111
112 inline
113void evaluateCoefficients(const int rowptr[nSpace],
114 const int colind[nnz],
115 const double v[nSpace],
116 const double alpha_L,
117 const double alpha_T,
118 const double Dm,
119 const double thetaW,
120 const double rho_f,
121 const double rho_s,
122 const double& u,
123 double& rho_out,
124 double& m,
125 double& dm,
126 double f[nSpace],
127 double df[nSpace],
128 double a[nnz],
129 double da[nnz])
130 {
131 const double epsilon = (rho_s - rho_f)/rho_f;
132 rho_out = rho_f*(1.0 + epsilon*u);
133 const double drho_du = rho_f*epsilon;
134 const double rho_transport = rho_out;
135 m = thetaW*rho_transport*u;
136 dm = thetaW*(rho_transport + u*drho_du);
137 double v_mag = 0.0;
138 for (int I=0; I<nSpace; I++)
139 v_mag += v[I]*v[I];
140 v_mag = std::sqrt(v_mag);
141 const double v_pore_mag = (thetaW > 1.0e-8) ? v_mag / thetaW : 0.0;
142
143 double alpha_L_eff = alpha_L;
144 double alpha_T_eff = alpha_T;
145 double v_unit[nSpace] = {0.0};
146 if (v_mag > 1.0e-10)
147 for (int I=0; I<nSpace; I++)
148 v_unit[I] = v[I]/v_mag;
149
150 for (int I=0; I < nSpace; I++)
151 {
152 f[I] = rho_transport*v[I]*u;
153 df[I] = (rho_transport + u*drho_du)*v[I];
154 for (int ii = rowptr[I]; ii < rowptr[I + 1]; ii++)
155 {
156 const int J = colind[ii];
157 const double deltaIJ = (I == J) ? 1.0 : 0.0;
158 const double dispersion_tensor =
159 Dm*deltaIJ +
160 alpha_L_eff*v_unit[I]*v_unit[J]*v_pore_mag +
161 alpha_T_eff*v_pore_mag*(deltaIJ - v_unit[I]*v_unit[J]);
162 a[ii] = thetaW*rho_transport*dispersion_tensor;
163 da[ii] = thetaW*drho_du*dispersion_tensor;
164 }
165
166 }
167 }
168inline
169double inversevaluateCoefficients(const double storage,
170 const double porosity,
171 const double rho_f,
172 const double rho_s)
173 {
174 const double mass_scale = std::max(porosity*rho_f, 1.0e-14);
175 const double epsilon = (rho_s - rho_f)/rho_f;
176 const double rhs = storage/mass_scale;
177 if (std::fabs(epsilon) < 1.0e-14)
178 return rhs;
179 const double discriminant = std::max(1.0 + 4.0*epsilon*rhs, 0.0);
180 const double sqrt_discriminant = std::sqrt(discriminant);
181 if (epsilon > 0.0)
182 return 2.0*rhs/(1.0 + sqrt_discriminant);
183 return (-1.0 + sqrt_discriminant)/(2.0*epsilon);
184 }
185
186inline
187 void exteriorNumericalDiffusiveFlux(int* rowptr,
188 int* colind,
189 const int& isDOFBoundary,
190 const int& isDiffusiveFluxBoundary,
191 const double n[nSpace],
192 double* bc_a,
193 const double& bc_u,
194 const double& bc_flux,
195 double* a,
196 const double grad_potential[nSpace],
197 const double& u,
198 const double& penalty,
199 double& flux)
200 {
201 double diffusiveVelocityComponent_I;
202 double penaltyFlux;
203 double max_a;
204 if (isDiffusiveFluxBoundary == 1)
205 {
206 flux = bc_flux;
207 }
208 else if (isDOFBoundary == 1)
209 {
210 flux = 0.0;
211 max_a = 0.0;
212 for (int I = 0; I < nSpace; I++)
213 {
214 diffusiveVelocityComponent_I = 0.0;
215 for (int m = rowptr[I]; m < rowptr[I+1]; m++)
216 {
217 diffusiveVelocityComponent_I -= a[m] * grad_potential[colind[m]];
218 max_a = fmax(max_a, a[m]);
219 }
220 flux += diffusiveVelocityComponent_I * n[I];
221 }
222 penaltyFlux = max_a * penalty * (u - bc_u);
223 flux += penaltyFlux;
224 }
225 else
226 {
227 //std::cerr << "warning, diffusion term with no boundary condition set, setting diffusive flux to 0.0" << std::endl;
228 flux = 0.0;
229 }
230 }
231
232 inline
233 double ExteriorNumericalDiffusiveFluxJacobian(int* rowptr,
234 int* colind,
235 const int& isDOFBoundary,
236 const int& isDiffusiveFluxBoundary,
237 const double n[nSpace],
238 double* a,
239 const double& v,
240 const double grad_v[nSpace],
241 const double& penalty)
242 {
243 double dvel_I;
244 double tmp = 0.0;
245 double max_a = 0.0;
246 if ((isDiffusiveFluxBoundary == 0) && (isDOFBoundary == 1))
247 {
248 for (int I = 0; I < nSpace; I++)
249 {
250 dvel_I = 0.0;
251 for (int m = rowptr[I]; m < rowptr[I + 1]; m++)
252 {
253 dvel_I -= a[m] * grad_v[colind[m]];
254 max_a = fmax(max_a, a[m]);
255 }
256 tmp += dvel_I * n[I];
257 }
258 tmp += max_a * penalty * v;
259 }
260 return tmp;
261 }
262
263 inline
264 void calculateSubgridError_tau(const double& elementDiameter,
265 const double& dmt,
266 const double df[nSpace],
267 double& cfl,
268 double& tau)
269 {
270 //regular elements
271 double h,nrm_v,oneByAbsdt;
272 h = elementDiameter;
273 nrm_v=0.0;
274 for(int I=0;I<nSpace;I++)
275 nrm_v+=df[I]*df[I];
276 nrm_v = sqrt(nrm_v);
277 cfl = nrm_v/h;
278 oneByAbsdt = fabs(dmt);
279 tau = 1.0/(2.0*nrm_v/h + oneByAbsdt + 1.0e-8);
280 }
281
282 inline
283 void calculateSubgridError_tau(const double& Ct_sge,
284 const double G[nSpace*nSpace],
285 const double& A0,
286 const double Ai[nSpace],
287 double& tau_v,
288 double& cfl)
289 {
290 //metric-based tau for arbitrarily shaped elements
291 double v_d_Gv=0.0;
292 for(int I=0;I<nSpace;I++)
293 {for (int J=0;J<nSpace;J++)
294 v_d_Gv += Ai[I]*G[I*nSpace+J]*Ai[J];
295 v_d_Gv += Ai[I]*G[I*nSpace+I];
296 for(int J=0;J<nSpace;J++)
297 {
298 if(J!=I)
299 v_d_Gv += 2.0*Ai[I]*G[I*nSpace+J];
300 }
301 }
302 tau_v = 1.0/sqrt(Ct_sge*A0*A0 + v_d_Gv + 1.0e-8);
303 }
304
305 inline
306 void calculateNumericalDiffusion(const double& shockCapturingDiffusion,
307 const double& elementDiameter,
308 const double& strong_residual,
309 const double grad_u[nSpace],
310 double& numDiff)
311 {
312 double h,
313 num,
314 den,
315 n_grad_u;
316 h = elementDiameter;
317 n_grad_u = 0.0;
318 for (int I=0;I<nSpace;I++)
319 n_grad_u += grad_u[I]*grad_u[I];
320 num = shockCapturingDiffusion*0.5*h*fabs(strong_residual);
321 den = sqrt(n_grad_u) + 1.0e-8;
322 numDiff = num/den;
323 }
324
325 inline
326 void exteriorNumericalAdvectiveFlux(const int& isDOFBoundary_u,
327 const int& isFluxBoundary_u,
328 const int& forceStrongConditions,
329 const double n[nSpace],
330 const double& bc_flux_u,
331 const double f[nSpace],
332 const double bc_f[nSpace],
333 const double velocity[nSpace],
334 double& flux)
335 {
336
337 double flow=0.0;
338 for (int I=0; I < nSpace; I++)
339 flow += n[I]*velocity[I];
340
341 if (isDOFBoundary_u == 1)
342 {
343 flux = 0.0;
344 if (forceStrongConditions == 1)
345 for (int I=0; I < nSpace; I++) flux += n[I]*bc_f[I];
346 else if (flow >= 0.0)
347 for (int I=0; I < nSpace; I++) flux += n[I]*f[I];
348 else
349 for (int I=0; I < nSpace; I++) flux += n[I]*bc_f[I];
350 }
351 else if (isFluxBoundary_u == 1)
352 {
353 flux = bc_flux_u;
354 }
355 else
356 {
357 flux = 0.0;
358 if (flow >= 0.0)
359 for (int I=0; I < nSpace; I++) flux += n[I]*f[I];
360 // else: open boundary with inflow, no external trace — flux = 0
361 }
362 }
363
364 inline
365 void exteriorNumericalAdvectiveFluxDerivative(const int& isDOFBoundary_u,
366 const int& isFluxBoundary_u,
367 const int& forceStrongConditions,
368 const double n[nSpace],
369 const double velocity[nSpace],
370 double& dflux)
371 {
372 double flow=0.0;
373 for (int I=0; I < nSpace; I++)
374 flow += n[I]*velocity[I];
375
376 dflux=0.0;//default to no flux
377 if (isDOFBoundary_u == 1)
378 {
379 if (forceStrongConditions == 1)
380 {
381 dflux = 0.0;
382 }
383 else if (flow >= 0.0)
384 {
385 dflux = flow;
386 }
387 else
388 {
389 dflux = 0.0;
390 }
391 }
392 else if (isFluxBoundary_u == 1)
393 {
394 dflux = 0.0;
395 }
396 else
397 {
398 if (flow >= 0.0)
399 {
400 dflux = flow;
401 }
402 }
403 }
404inline
405
406 void exteriorNumericalDiffusiveFluxDerivative(const int& isDOFBoundary,
407 const int& isDiffusiveFluxBoundary,
408 const int rowptr[nSpace],
409 const int colind[nnz],
410 const double n[nSpace],
411 const double a[nnz],
412 const double da[nnz],
413 const double grad_psi[nSpace],
414 const double grad_v[nSpace],
415 const double& v,
416 const double penalty,
417 double& fluxJacobian)
418{
419 if (isDiffusiveFluxBoundary == 0 && isDOFBoundary == 1)
420 {
421 fluxJacobian = 0.0;
422 double max_a = 0.0;
423 for (int I = 0; I < nSpace; I++) {
424 for(int m=rowptr[I]; m<rowptr[I+1]; m++)
425 {
426 max_a = fmax(max_a, a[m]);
427 fluxJacobian -= (a[m] * grad_v[colind[m]] + da[m] * v * grad_psi[colind[m]]) * n[I];
428 }
429 fluxJacobian += max_a * penalty * v;
430 }
431 }
432 else
433 {
434 fluxJacobian = 0.0;
435 }
436}
437
438
439 void calculateResidual(arguments_dict& args)
440 {
441 double dt = args.scalar<double>("dt");
442 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
443 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
444 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
445 xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
446 double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
447 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
448 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
449 xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
450 xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
451 xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
452 xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
453 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
454 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
455 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
456 xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
457 xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
458 xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
459 xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
460 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
461 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
462 int nElements_global = args.scalar<int>("nElements_global");
463 double useMetrics = args.scalar<double>("useMetrics");
464 double alphaBDF = args.scalar<double>("alphaBDF");
465 int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
466 double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
467 double sc_uref = args.scalar<double>("sc_uref");
468 double sc_alpha = args.scalar<double>("sc_alpha");
469 xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
470 xt::pyarray<int>& r_l2g = args.array<int>("r_l2g");
471 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
472 double degree_polynomial = args.scalar<double>("degree_polynomial");
473 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
474 xt::pyarray<double>& u_dof_old = args.array<double>("u_dof_old");
475 xt::pyarray<double>& velocity = args.array<double>("velocity");
476 xt::pyarray<double>& velocity_old = args.array<double>("velocity_old");
477 xt::pyarray<double>& q_m = args.array<double>("q_m");
478 xt::pyarray<double>& q_u = args.array<double>("q_u");
479 xt::pyarray<double>& q_porosity = args.array<double>("q_porosity");
480 xt::pyarray<double>& q_porosity_old = args.array<double>("q_porosity_old");
481 xt::pyarray<double>& q_rho = args.array<double>("q_rho");
482 xt::pyarray<double>& q_rho_old = args.array<double>("q_rho_old");
483
484 xt::pyarray<double>& q_r = args.array<double>("q_r");
485 const double alpha_L = args.scalar<double>("alpha_L");
486 const double alpha_T = args.scalar<double>("alpha_T");
487 const double Dm = args.scalar<double>("Dm");
488 // const int dispersion_type_int = args.scalar<int>("dispersion_type");
489 // const double theta_s = args.scalar<double>("theta_s");
490 // const double theta_r = args.scalar<double>("theta_r");
491 // const double power_law_exponent = args.scalar<double>("power_law_exponent");
492 // const double velocity_exponent = args.scalar<double>("velocity_exponent");
493 const double rho_f = args.scalar<double>("rho_f");
494 const double rho_s = args.scalar<double>("rho_s");
495 int forceStrongConditions = args.scalar<int>("forceStrongConditions");
496 // DISPERSION DISPERSION_TYPE = static_cast<DISPERSION>(dispersion_type_int);
497 xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
498 xt::pyarray<double>& q_dV = args.array<double>("q_dV");
499 xt::pyarray<double>& q_dV_last = args.array<double>("q_dV_last");
500 xt::pyarray<double>& cfl = args.array<double>("cfl");
501 xt::pyarray<double>& edge_based_cfl = args.array<double>("edge_based_cfl");
502 xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
503 xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
504 int offset_u = args.scalar<int>("offset_u");
505 int stride_u = args.scalar<int>("stride_u");
506 xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
507 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
508 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
509 xt::pyarray<int>& elementBoundaryMaterialTypes = args.array<int>("elementBoundaryMaterialTypes");
510 xt::pyarray<int>& isExteriorBoundaryPhysical = args.array<int>("isExteriorBoundaryPhysical");
511 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
512 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
513 xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
514 xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
515 xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
516 xt::pyarray<int>& isFluxBoundary_u = args.array<int>("isFluxBoundary_u");
517 xt::pyarray<double>& ebqe_bc_flux_u_ext = args.array<double>("ebqe_bc_flux_u_ext");
518 xt::pyarray<double>& ebqe_bc_diffusiveFlux_u_ext = args.array<double>("ebqe_bc_diffusiveFlux_u_ext");
519 xt::pyarray<double>& ebqe_porosity = args.array<double>("ebqe_porosity");
520 xt::pyarray<double>& ebqe_rho = args.array<double>("ebqe_rho");
521
522 double epsFact = args.scalar<double>("epsFact");
523 xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
524 xt::pyarray<double>& ebqe_flux = args.array<double>("ebqe_flux");
525 int stage = args.scalar<int>("stage");
526 xt::pyarray<double>& uTilde_dof = args.array<double>("uTilde_dof");
527 double cE = args.scalar<double>("cE");
528 double cMax = args.scalar<double>("cMax");
529 double cK = args.scalar<double>("cK");
530 double uL = args.scalar<double>("uL");
531 double uR = args.scalar<double>("uR");
532 int numDOFs = args.scalar<int>("numDOFs");
533 int NNZ = args.scalar<int>("NNZ");
534 xt::pyarray<int>& csrRowIndeces_DofLoops = args.array<int>("csrRowIndeces_DofLoops");
535 xt::pyarray<int>& csrColumnOffsets_DofLoops = args.array<int>("csrColumnOffsets_DofLoops");
536 xt::pyarray<int>& csrRowIndeces_CellLoops = args.array<int>("csrRowIndeces_CellLoops");
537 xt::pyarray<int>& csrColumnOffsets_CellLoops = args.array<int>("csrColumnOffsets_CellLoops");
538 xt::pyarray<int>& csrColumnOffsets_eb_CellLoops = args.array<int>("csrColumnOffsets_eb_CellLoops");
539 xt::pyarray<double>& ML = args.array<double>("ML");
540 int LUMPED_MASS_MATRIX = args.scalar<int>("LUMPED_MASS_MATRIX");
541 STABILIZATION STABILIZATION_TYPE = static_cast<STABILIZATION>(args.scalar<int>("STABILIZATION_TYPE"));
542 ENTROPY ENTROPY_TYPE = static_cast<ENTROPY>(args.scalar<int>("ENTROPY_TYPE"));
543 //STABILIZATION STABILIZATION_TYPE{args.scalar<int>("STABILIZATION_TYPE")};
544 //ENTROPY ENTROPY_TYPE{args.scalar<int>("ENTROPY_TYPE")};
545 xt::pyarray<double>& uLow = args.array<double>("uLow");
546 xt::pyarray<double>& dLow = args.array<double>("dLow");
547 xt::pyarray<double>& dt_times_dH_minus_dL = args.array<double>("dt_times_dH_minus_dL");
548 xt::pyarray<double>& min_u_bc = args.array<double>("min_u_bc");
549 xt::pyarray<double>& max_u_bc = args.array<double>("max_u_bc");
550 xt::pyarray<double>& quantDOFs = args.array<double>("quantDOFs");
551 // Stage 3 (kinetic dissolution). Adds R_diss = k_d * S_n * S_w *
552 // (c_sat - c) per DOF to the mass update, scaled by theta_w * rho_w
553 // (so it has mass-rate units). When the flow model is single-phase
554 // (Richards), Sn_dof is zeros and R_diss vanishes.
555 xt::pyarray<double>& Sn_dof = args.array<double>("Sn_dof");
556 const double k_d = args.scalar<double>("k_d");
557 const double c_sat = args.scalar<double>("c_sat");
559 xt::pyarray<int>& a_rowptr = args.array<int>("a_rowptr");
560 xt::pyarray<int>& a_colind = args.array<int>("a_colind");
561 //xt::pyarray<double>& D = args.array<double>("D");
562 //initializeDToZero(D);
564 xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.array<int>("isDiffusiveFluxBoundary_u");
565 xt::pyarray<int>& isAdvectiveFluxBoundary_u = args.array<int>("isAdvectiveFluxBoundary_u");
566 xt::pyarray<double>& ebqe_bc_advectiveFlux_u_ext = args.array<double>("ebqe_bc_advectiveFlux_u_ext");
567 xt::pyarray<double>& ebqe_penalty_ext = args.array<double>("ebqe_penalty_ext");
569 double physicalDiffusion = args.scalar<double>("physicalDiffusion");
570 const double eb_adjoint_sigma = args.scalar<double>("eb_adjoint_sigma");
571
572 double meanEntropy = 0., meanOmega = 0., maxEntropy = -1E10, minEntropy = 1E10;
573 const double eps_rho = (rho_s - rho_f)/rho_f;
574 const double zL_mass = uL + eps_rho*uL*uL;
575 const double zR_mass = uR + eps_rho*uR*uR;
576 maxVel.resize(nElements_global, 0.0);
577 maxEntRes.resize(nElements_global, 0.0);
578 double Ct_sge = 4.0;
579 if (STABILIZATION_TYPE==STABILIZATION::EntropyViscosity or
580 STABILIZATION_TYPE==STABILIZATION::SmoothnessIndicator or
581 STABILIZATION_TYPE==STABILIZATION::Kuzmin or
582 STABILIZATION_TYPE==STABILIZATION::ImplicitEV)
583 {
584 TransportMatrix.resize(NNZ,0.0);
585 DiffusionMatrix.resize(NNZ,0.0);
586 TransposeTransportMatrix.resize(NNZ,0.0);
587 m_dof.resize(numDOFs,0.0);
588 theta_dof_proj.resize(numDOFs,0.0);
589 rho_dof_proj.resize(numDOFs,0.0);
590 ML_mass_proj.resize(numDOFs,0.0);
591 // ONLY the porosity gets a quadrature->DOF projection: theta lives at
592 // quadrature points and has no nodal representation. rho and the
593 // conservative variable m are evaluated DIRECTLY at the DOF from
594 // u_dof_old. Projecting them would apply (M*u)_i/ML_i -- a smoothing
595 // filter, because ML_i is exactly the row sum of the consistent mass
596 // matrix -- so m_dof would NOT equal u_dof_old even at rho=theta=1 and
597 // the scheme would not reduce to the constant-density one. Evaluating
598 // at the DOF also makes inversevaluateCoefficients(m_dof[i],...) return
599 // u_dof_old[i] exactly, and makes this m^n consistent with the nodal
600 // m^{n+1} = theta_i*rho(c_i)*c_i that the ImplicitEV branch builds
601 // (otherwise ML_i*(m^{n+1}-m^n)/dt is nonzero at steady state).
602 for (int eN=0; eN<nElements_global; eN++)
603 for (int k=0; k<nQuadraturePoints_element; k++)
604 {
605 int eN_k = eN*nQuadraturePoints_element + k;
606 double jac[nSpace*nSpace], jacDet, jacInv[nSpace*nSpace], x, y, z;
607 ck.calculateMapping_element(eN,
608 k,
609 mesh_dof.data(),
610 mesh_l2g.data(),
611 mesh_trial_ref.data(),
612 mesh_grad_trial_ref.data(),
613 jac,
614 jacDet,
615 jacInv,
616 x,y,z);
617 const double dV = fabs(jacDet)*dV_ref.data()[k];
618 const double theta_k = q_porosity_old.data()[eN_k];
619 for (int i=0; i<nDOF_test_element; i++)
620 {
621 int eN_i = eN*nDOF_test_element+i;
622 const int gi = u_l2g.data()[eN_i];
623 const double w = u_test_ref.data()[k*nDOF_trial_element+i]*dV;
624 theta_dof_proj[gi] += theta_k*w;
625 ML_mass_proj[gi] += w;
626 }
627 }
628 for (int i=0; i<numDOFs; i++)
629 {
630 if (ML_mass_proj[i] > 1.0e-14)
631 theta_dof_proj[i] /= ML_mass_proj[i];
632 else
633 theta_dof_proj[i] = 1.0;
634 // nodal density and conservative variable at t^n
635 const double un_i = u_dof_old.data()[i];
636 rho_dof_proj[i] = rho_f*(1.0 + eps_rho*un_i);
637 m_dof[i] = theta_dof_proj[i]*rho_dof_proj[i]*un_i;
638 }
639 // compute entropy and init global_entropy_residual and boundary_integral
640 psi.resize(numDOFs,0.0);
641 eta.resize(numDOFs,0.0);
642 global_entropy_residual.resize(numDOFs,0.0);
643 boundary_integral.resize(numDOFs,0.0);
644 for (int i=0; i<numDOFs; i++)
645 {
646 // NODAL ENTROPY //
647 if (STABILIZATION_TYPE==STABILIZATION::EntropyViscosity) //EV stab
648 {
649 if (ENTROPY_TYPE == ENTROPY::POWER)
650 eta[i] = EPOWER(m_dof[i],uL,uR);
651 else
652 {
653 const double mass_scale_i = fmax(theta_dof_proj[i]*rho_f, 1.0e-14);
654 const double z_i = m_dof[i]/mass_scale_i;
655 eta[i] = ELOG(z_i,zL_mass,zR_mass);
656 }
657 global_entropy_residual[i]=0.;
658 }
659 boundary_integral[i]=0.;
660 }
661 }
662 //
663 //loop over elements to compute volume integrals and load them into element and global residual
664 //
665 //eN is the element index
666 //eN_k is the quadrature point index for a scalar
667 //eN_k_nSpace is the quadrature point index for a vector
668 //eN_i is the element test function index
669 //eN_j is the element trial function index
670 //eN_k_j is the quadrature point index for a trial function
671 //eN_k_i is the quadrature point index for a trial function
672 for(int eN=0;eN<nElements_global;eN++)
673 {
674 //declare local storage for element residual and initialize
675 double
676 elementResidual_u[nDOF_test_element],
677 element_entropy_residual[nDOF_test_element];
678 double elementTransport[nDOF_test_element][nDOF_trial_element];
679 double elementDiffusion[nDOF_test_element][nDOF_trial_element];
680 double elementTransposeTransport[nDOF_test_element][nDOF_trial_element];
681 for (int i=0;i<nDOF_test_element;i++)
682 {
683 elementResidual_u[i]=0.0;
684 }//i
685 if (STABILIZATION_TYPE==STABILIZATION::EntropyViscosity or
686 STABILIZATION_TYPE==STABILIZATION::SmoothnessIndicator or
687 STABILIZATION_TYPE==STABILIZATION::Kuzmin or
688 STABILIZATION_TYPE==STABILIZATION::ImplicitEV)
689 {
690 for (int i=0;i<nDOF_test_element;i++)
691 {
692 element_entropy_residual[i]=0.0;
693 for (int j=0;j<nDOF_trial_element;j++)
694 {
695 elementTransport[i][j]=0.0;
696 elementDiffusion[i][j]=0.0;
697 elementTransposeTransport[i][j]=0.0;
698 }
699 }
700 }
701 //loop over quadrature points and compute integrands
702 for (int k=0;k<nQuadraturePoints_element;k++)
703 {
704 //compute indeces and declare local storage
705 int eN_k = eN*nQuadraturePoints_element+k,
706 eN_k_nSpace = eN_k*nSpace,
707 eN_nDOF_trial_element = eN*nDOF_trial_element;
708 //int index_D = eN_k * a_rowptr.data()[nSpace];
709 double
710 entVisc_minus_artComp,
711 u=0.0,un=0.0,
712 grad_u[nSpace],grad_u_old[nSpace],grad_uTilde[nSpace],
713 rho_out=0.0,rho_out_old=0.0,
714 m=0.0,dm=0.0,mn=0.0,dmn=0.0,
715 H=0.0,Hn=0.0,HTilde=0.0,
716 f[nSpace],fn[nSpace],df[nSpace],dfn[nSpace],
718 //a[nSpace], da[nSpace], an[nSpace], dan[nSpace],
719 a[nnz], da[nnz], an[nnz], dan[nnz],
720 m_t=0.0,dm_t=0.0,
721 pdeResidual_u=0.0,
722 Lstar_u[nDOF_test_element],
723 subgridError_u=0.0,
724 tau=0.0,tau0=0.0,tau1=0.0,
725 numDiff0=0.0,numDiff1=0.0,
726 jac[nSpace*nSpace],
727 jacDet,
728 jacInv[nSpace*nSpace],
729 u_grad_trial[nDOF_trial_element*nSpace],
730 u_test_dV[nDOF_trial_element],
731 u_grad_test_dV[nDOF_test_element*nSpace],
732 dV,x,y,z,xt,yt,zt,
733 G[nSpace*nSpace],G_dd_G,tr_G,
734 // for entropy residual
735 aux_entropy_residual=0.0, DENTROPY_un, DENTROPY_uni;//norm_Rv;
736
737 ck.calculateMapping_element(eN,
738 k,
739 mesh_dof.data(),
740 mesh_l2g.data(),
741 mesh_trial_ref.data(),
742 mesh_grad_trial_ref.data(),
743 jac,
744 jacDet,
745 jacInv,
746 x,y,z);
747 ck.calculateMappingVelocity_element(eN,
748 k,
749 mesh_velocity_dof.data(),
750 mesh_l2g.data(),
751 mesh_trial_ref.data(),
752 xt,yt,zt);
753 //get the physical integration weight
754 dV = fabs(jacDet)*dV_ref.data()[k];
755 ck.calculateG(jacInv,G,G_dd_G,tr_G);
756 //get the trial function gradients
757 ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],
758 jacInv,
759 u_grad_trial);
760 //get the solution
761 ck.valFromDOF(u_dof.data(),
762 &u_l2g.data()[eN_nDOF_trial_element],
763 &u_trial_ref.data()[k*nDOF_trial_element],
764 u);
765 ck.valFromDOF(u_dof_old.data(),
766 &u_l2g.data()[eN_nDOF_trial_element],
767 &u_trial_ref.data()[k*nDOF_trial_element],
768 un);
769 //get the solution gradients
770 ck.gradFromDOF(u_dof.data(),
771 &u_l2g.data()[eN_nDOF_trial_element],
772 u_grad_trial,
773 grad_u);
774 ck.gradFromDOF(u_dof_old.data(),
775 &u_l2g.data()[eN_nDOF_trial_element],
776 u_grad_trial,
777 grad_u_old);
778 ck.gradFromDOF(uTilde_dof.data(),
779 &u_l2g.data()[eN_nDOF_trial_element],
780 u_grad_trial,
781 grad_uTilde);
782 //precalculate test function products with integration weights
783 for (int j=0;j<nDOF_trial_element;j++)
784 {
785 u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
786 for (int I=0;I<nSpace;I++)
787 {
788 u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
789 }
790 }
791
792 //
793 //
794 //calculate pde coefficients at quadrature points
795
796 evaluateCoefficients(a_rowptr.data(),
797 a_colind.data(),
798 &velocity.data()[eN_k_nSpace],
799 alpha_L,
800 alpha_T,
801 Dm,
802 q_porosity.data()[eN*nQuadraturePoints_element+k],
803 rho_f,
804 rho_s,
805 u,
806 rho_out,
807 m,
808 dm,
809 f,
810 df,
811 a,
812 da);
813 q_rho.data()[eN_k] = rho_out;
814
815 evaluateCoefficients(a_rowptr.data(),
816 a_colind.data(),
817 &velocity_old.data()[eN_k_nSpace],
818 alpha_L,
819 alpha_T,
820 Dm,
821 q_porosity_old.data()[eN*nQuadraturePoints_element+k],
822 rho_f,
823 rho_s,
824 un,
825 rho_out_old,
826 mn,
827 dmn,
828 fn,
829 dfn,
830 an,
831 dan);
832 //an= &q_a.data()[eN_k * sd_rowptr.data()[nSpace]];
833
834 //
835 //moving mesh
836 //
837 double mesh_velocity[3];
838 mesh_velocity[0] = xt;
839 mesh_velocity[1] = yt;
840 mesh_velocity[2] = zt;
841
842 for (int I=0;I<nSpace;I++)
843 {
844 f[I] -= MOVING_DOMAIN*m*mesh_velocity[I];
845 df[I] -= MOVING_DOMAIN*dm*mesh_velocity[I];
846 fn[I] -= MOVING_DOMAIN*mn*mesh_velocity[I];
847 dfn[I] -= MOVING_DOMAIN*dmn*mesh_velocity[I];
848 }
849 //
850 //calculate time derivative at quadrature points
851 //
852 if (q_dV_last.data()[eN_k] <= -100)
853 q_dV_last.data()[eN_k] = dV;
854 q_dV.data()[eN_k] = dV;
855 ck.bdf(alphaBDF,
856 q_m_betaBDF.data()[eN_k]*q_dV_last.data()[eN_k]/dV,//ensure prior mass integral is correct for m_t with BDF1
857 m,
858 dm,
859 m_t,
860 dm_t);
861
862 const double thetaW_k = std::max(q_porosity_old.data()[eN_k], 1.0e-8);
863 double dfn_pore[nSpace];
864 for (int I=0; I<nSpace; I++) dfn_pore[I] = dfn[I] / thetaW_k;
865
866 if (STABILIZATION_TYPE==STABILIZATION::TaylorGalerkinEV)
867 {
868 double normVel=0., norm_grad_un=0.;
869 for (int I=0;I<nSpace;I++)
870 {
871 Hn += dfn[I]*grad_u_old[I];
872 HTilde += dfn[I]*grad_uTilde[I];
873 fn[I] = dfn[I]*un-MOVING_DOMAIN*m*mesh_velocity[I];//cek check this for moving domain
874 H += dfn[I]*grad_u[I];
875 normVel += dfn[I]*df[I];
876 norm_grad_un += grad_u_old[I]*grad_u_old[I];
877 }
878 normVel = std::sqrt(normVel);
879 norm_grad_un = std::sqrt(norm_grad_un)+1E-10;
880
881 // calculate CFL
882 calculateCFL(elementDiameter.data()[eN]/degree_polynomial,dfn_pore,cfl.data()[eN_k]);
883
884
885 // compute max velocity at cell
886 maxVel[eN] = fmax(normVel,maxVel[eN]);
887
888 // Strong entropy residual
889 double entRes = (EPOWER(u,0,1)-EPOWER(un,0,1))/dt + 0.5*(DEPOWER(u,0,1)*H + DEPOWER(un,0,1)*Hn);
890 maxEntRes[eN] = fmax(maxEntRes[eN],fabs(entRes));
891
892 // Quantities for normalization factor //
893 meanEntropy += EPOWER(u,0,1)*dV;
894 meanOmega += dV;
895 maxEntropy = fmax(maxEntropy,EPOWER(u,0,1));
896 minEntropy = fmin(minEntropy,EPOWER(u,0,1));
897
898 // artificial compression
899 double hK=elementDiameter.data()[eN]/degree_polynomial;
900 entVisc_minus_artComp = fmax(1-cK*fmax(un*(1-un),0)/hK/norm_grad_un,0);
901 }
902 else if (STABILIZATION_TYPE==STABILIZATION::VMS)
903 {
904 //
905 //calculate subgrid error (strong residual and adjoint)
906 //
907 //calculate strong residual
908 pdeResidual_u = ck.Mass_strong(m_t) + ck.Advection_strong(df,grad_u);
909 //calculate adjoint
910 for (int i=0;i<nDOF_test_element;i++)
911 {
912 int i_nSpace = i*nSpace;
913 Lstar_u[i] = ck.Advection_adjoint(df,&u_grad_test_dV[i_nSpace]);
914 }
915 //calculate tau and tau*Res
916 calculateSubgridError_tau(elementDiameter.data()[eN],dm_t,df,cfl.data()[eN_k],tau0);
917 calculateSubgridError_tau(Ct_sge,
918 G,
919 dm_t,
920 df,
921 tau1,
922 cfl.data()[eN_k]);
923 tau = useMetrics*tau1+(1.0-useMetrics)*tau0;
924
925 subgridError_u = -tau*pdeResidual_u;
926 //
927 //calculate shock capturing diffusion
928 //
929 ck.calculateNumericalDiffusion(shockCapturingDiffusion,
930 elementDiameter.data()[eN],
931 pdeResidual_u,
932 grad_u,
933 numDiff0);
934 ck.calculateNumericalDiffusion(shockCapturingDiffusion,
935 sc_uref,
936 sc_alpha,
937 G,
938 G_dd_G,
939 pdeResidual_u,
940 grad_u,
941 numDiff1);
942 q_numDiff_u.data()[eN_k] = useMetrics*numDiff1+(1.0-useMetrics)*numDiff0;
943 }
944 else if (STABILIZATION_TYPE==STABILIZATION::EntropyViscosity)
945 {
946 aux_entropy_residual = m_t;
947 for (int I=0;I<nSpace;I++)
948 aux_entropy_residual += dfn[I]*grad_u_old[I];
949 if (ENTROPY_TYPE==ENTROPY::POWER)
950 DENTROPY_un = DEPOWER(mn,uL,uR);
951 else
952 {
953 const double mass_scale_k = fmax(q_porosity_old.data()[eN_k]*rho_f, 1.0e-14);
954 const double z_n = mn/mass_scale_k;
955 DENTROPY_un = DELOG(z_n,zL_mass,zR_mass)/mass_scale_k;
956 }
957 calculateCFL(elementDiameter.data()[eN]/degree_polynomial,dfn_pore,cfl.data()[eN_k]);
958 }
959 else
960 calculateCFL(elementDiameter.data()[eN]/degree_polynomial,dfn_pore,cfl.data()[eN_k]);
961
962 for(int i=0;i<nDOF_test_element;i++)
963 {
964 int i_nSpace=i*nSpace;
965 if (STABILIZATION_TYPE==STABILIZATION::TaylorGalerkinEV)
966 {
967 if (stage == 1)
968 elementResidual_u[i] +=
969 ck.Mass_weak(dt*m_t,u_test_dV[i]) + // time derivative
970 1./3*dt*(ck.Advection_weak(fn,&u_grad_test_dV[i_nSpace]) +
971 ck.Diffusion_weak(a_rowptr.data(),a_colind.data(),a,grad_u,&u_grad_test_dV[i_nSpace])+
972 ck.NumericalDiffusion(physicalDiffusion, grad_u_old, &u_grad_test_dV[i_nSpace])) +
973 1./9*dt*dt*ck.NumericalDiffusion(Hn,dfn,&u_grad_test_dV[i_nSpace]) +
974 1./3*dt*entVisc_minus_artComp*ck.NumericalDiffusion(q_numDiff_u_last.data()[eN_k]+physicalDiffusion,
975 grad_u_old,
976 &u_grad_test_dV[i_nSpace]);
977 // TODO: Add part about moving mesh
978 else //stage 2
979 elementResidual_u[i] +=
980 ck.Mass_weak(dt*m_t,u_test_dV[i]) + // time derivative
981 dt*(ck.Advection_weak(fn,&u_grad_test_dV[i_nSpace]) +
982 ck.Diffusion_weak(a_rowptr.data(),a_colind.data(),an,grad_u,&u_grad_test_dV[i_nSpace])+
983 ck.NumericalDiffusion(physicalDiffusion, grad_u_old, &u_grad_test_dV[i_nSpace])) +
984 0.5*dt*dt*ck.NumericalDiffusion(HTilde,dfn,&u_grad_test_dV[i_nSpace]) +
985 dt*entVisc_minus_artComp*ck.NumericalDiffusion(q_numDiff_u_last.data()[eN_k]+physicalDiffusion,
986 grad_u_old,
987 &u_grad_test_dV[i_nSpace]);
988 }
989 else if (STABILIZATION_TYPE==STABILIZATION::VMS)
990 {
991 elementResidual_u[i] +=
992 ck.Mass_weak(m_t,u_test_dV[i]) +
993 ck.Advection_weak(f,&u_grad_test_dV[i_nSpace]) +
994 ck.Diffusion_weak(a_rowptr.data(),a_colind.data(),a,grad_u,&u_grad_test_dV[i_nSpace]) +
995 ck.SubgridError(subgridError_u,Lstar_u[i]) +
996 ck.NumericalDiffusion(q_numDiff_u_last.data()[eN_k] + physicalDiffusion,//todo add full sparse diffusion terms
997 grad_u,
998 &u_grad_test_dV[i_nSpace]);
999 }
1000 else if(STABILIZATION_TYPE==STABILIZATION::EntropyViscosity or
1001 STABILIZATION_TYPE==STABILIZATION::SmoothnessIndicator or
1002 STABILIZATION_TYPE==STABILIZATION::Kuzmin or
1003 STABILIZATION_TYPE==STABILIZATION::ImplicitEV)
1004 {
1005 int eN_i=eN*nDOF_test_element+i;
1006 if (STABILIZATION_TYPE==STABILIZATION::EntropyViscosity) // EV stab
1007 {
1008 element_entropy_residual[i] += DENTROPY_un*aux_entropy_residual*u_test_dV[i];
1009 }
1010 elementResidual_u[i] += (u-un)*u_test_dV[i];
1011 double* adv_df = (STABILIZATION_TYPE==STABILIZATION::ImplicitEV) ? df : dfn;
1012 for(int j=0;j<nDOF_trial_element;j++)
1013 {
1014 int j_nSpace = j*nSpace;
1015 int i_nSpace = i*nSpace;
1016 elementTransport[i][j] +=
1017 ck.AdvectionJacobian_weak(adv_df,
1018 u_trial_ref.data()[k*nDOF_trial_element+j],
1019 &u_grad_test_dV[i_nSpace])
1020 +
1021 ck.SimpleDiffusionJacobian_weak(a_rowptr.data(),
1022 a_colind.data(),
1023 a,
1024 &u_grad_trial[j_nSpace],
1025 &u_grad_test_dV[i_nSpace]);
1026
1027
1028
1029
1030 elementDiffusion[i][j] += ck.NumericalDiffusionJacobian(physicalDiffusion,
1031 &u_grad_trial[j_nSpace],
1032 &u_grad_test_dV[i_nSpace]);
1033 elementTransposeTransport[i][j] += ck.AdvectionJacobian_weak(adv_df,
1034 u_trial_ref.data()[k*nDOF_trial_element+i],
1035 &u_grad_test_dV[j_nSpace])+
1036 ck.SimpleDiffusionJacobian_weak(a_rowptr.data(),
1037 a_colind.data(),
1038 a,
1039 &u_grad_trial[j_nSpace],
1040 &u_grad_test_dV[i_nSpace]);
1041
1042
1043 }
1044 }
1045 else
1046 {
1047 elementResidual_u[i] +=
1048 ck.Mass_weak(m_t,u_test_dV[i]) +
1049 ck.Advection_weak(f,&u_grad_test_dV[i_nSpace])+
1050 ck.Diffusion_weak(a_rowptr.data(),a_colind.data(),a,grad_u,&u_grad_test_dV[i_nSpace]);
1051 // +
1052 // ck.NumericalDiffusion(physicalDiffusion,//todo add full sparse diffusion terms
1053 // grad_u,
1054 // &u_grad_test_dV[i_nSpace]);
1055 }
1056 }//i
1057 //
1058 //save solution for other models
1059 //
1060 q_u.data()[eN_k] = u;
1061 q_m.data()[eN_k] = m;
1062 //logInteriorState("storeQuadrature", eN, k, x, y, z, u, m, f, a);
1063
1064 }//k
1065 //
1066 //load element into global residual and save element residual
1067 //
1068 for(int i=0;i<nDOF_test_element;i++)
1069 {
1070 int eN_i=eN*nDOF_test_element+i;
1071 int gi = offset_u+stride_u*u_l2g.data()[eN_i]; //global i-th index
1072 globalResidual.data()[gi] += elementResidual_u[i];
1073 if (STABILIZATION_TYPE==STABILIZATION::EntropyViscosity or
1074 STABILIZATION_TYPE==STABILIZATION::SmoothnessIndicator or
1075 STABILIZATION_TYPE==STABILIZATION::Kuzmin or
1076 STABILIZATION_TYPE==STABILIZATION::ImplicitEV)
1077 {
1078
1079 // distribute entropy_residual
1080 if (STABILIZATION_TYPE==STABILIZATION::EntropyViscosity) // EV Stab
1081 global_entropy_residual[gi] += element_entropy_residual[i];
1082 // distribute transport matrices
1083 for (int j=0;j<nDOF_trial_element;j++)
1084 {
1085 int eN_i_j = eN_i*nDOF_trial_element+j;
1086 TransportMatrix[csrRowIndeces_CellLoops.data()[eN_i] +
1087 csrColumnOffsets_CellLoops.data()[eN_i_j]] += elementTransport[i][j];
1088 DiffusionMatrix[csrRowIndeces_CellLoops.data()[eN_i] +
1089 csrColumnOffsets_CellLoops.data()[eN_i_j]] += elementDiffusion[i][j];
1090 TransposeTransportMatrix[csrRowIndeces_CellLoops.data()[eN_i] +
1091 csrColumnOffsets_CellLoops.data()[eN_i_j]]
1092 += elementTransposeTransport[i][j];
1093 }//j
1094 }//edge-based
1095 }//i
1096 }//eN
1097 //
1098 //loop over exterior element boundaries to calculate surface integrals and load into element and global residuals
1099 //
1100 //ebNE is the Exterior element boundary INdex
1101 //ebN is the element boundary INdex
1102 //eN is the element index
1103 for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
1104 {
1105 double min_u_bc_local = 1E10, max_u_bc_local = -1E10;
1106 int ebN = exteriorElementBoundariesArray.data()[ebNE],
1107 eN = elementBoundaryElementsArray.data()[ebN*2+0],
1108 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
1109 eN_nDOF_trial_element = eN*nDOF_trial_element;
1110 const int eN_out = elementBoundaryElementsArray.data()[ebN*2+1];
1111 const int ebFlag = elementBoundaryMaterialTypes.data()[ebN];
1112 // Only integrate true physical exterior boundaries; skip partition/non-physical faces.
1113 if (ebFlag <= 0 || isExteriorBoundaryPhysical.data()[ebNE] == 0 || eN_out >= 0)
1114 {
1115 continue;
1116 }
1117 double elementResidual_u[nDOF_test_element],
1118 fluxTransport[nDOF_test_element][nDOF_trial_element];
1119 for (int i=0;i<nDOF_test_element;i++)
1120 {
1121 elementResidual_u[i]=0.0;
1122 for (int j=0;j<nDOF_trial_element;j++)
1123 fluxTransport[i][j] = 0.0;
1124 }
1125 for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
1126 {
1127 int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
1128 ebNE_kb_nSpace = ebNE_kb*nSpace,
1129 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
1130 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
1131 double u_ext=0.0,
1132 grad_u_ext[nSpace],
1133 m_ext=0.0,
1134 dm_ext=0.0,
1135 f_ext[nSpace],
1136 df_ext[nSpace],
1138
1139 a_ext[nnz],
1140 da_ext[nnz],
1141
1142 bc_a_ext[nnz],
1143 bc_da_ext[nnz],
1144
1146 flux_ext=0.0,
1147 dflux_u_u_ext=0.0,
1148 bc_u_ext=0.0,
1149 bc_m_ext=0.0,
1150 bc_dm_ext=0.0,
1151
1152 flux_adv_ext=0.0,
1153 flux_diff_ext=0.0,
1154 difffluxjacobian_ext=0.0,
1155 bc_f_ext[nSpace],
1156 bc_df_ext[nSpace],
1157 jac_ext[nSpace*nSpace],
1158 jacDet_ext,
1159 jacInv_ext[nSpace*nSpace],
1160 boundaryJac[nSpace*(nSpace-1)],
1161 metricTensor[(nSpace-1)*(nSpace-1)],
1162 metricTensorDetSqrt,
1163 dS,
1164 u_test_dS[nDOF_test_element],
1165 u_grad_trial_trace[nDOF_trial_element*nSpace],
1166 u_grad_test_dS[nDOF_trial_element*nSpace],
1167 normal[nSpace],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
1168
1169 G[nSpace*nSpace],G_dd_G,tr_G;
1170
1171 //
1172 //calculate the solution and gradients at quadrature points
1173 //
1174 //compute information about mapping from reference element to physical element
1175 ck.calculateMapping_elementBoundary(eN,
1176 ebN_local,
1177 kb,
1178 ebN_local_kb,
1179 mesh_dof.data(),
1180 mesh_l2g.data(),
1181 mesh_trial_trace_ref.data(),
1182 mesh_grad_trial_trace_ref.data(),
1183 boundaryJac_ref.data(),
1184 jac_ext,
1185 jacDet_ext,
1186 jacInv_ext,
1187 boundaryJac,
1188 metricTensor,
1189 metricTensorDetSqrt,
1190 normal_ref.data(),
1191 normal,
1192 x_ext,y_ext,z_ext);
1193 ck.calculateMappingVelocity_elementBoundary(eN,
1194 ebN_local,
1195 kb,
1196 ebN_local_kb,
1197 mesh_velocity_dof.data(),
1198 mesh_l2g.data(),
1199 mesh_trial_trace_ref.data(),
1200 xt_ext,yt_ext,zt_ext,
1201 normal,
1202 boundaryJac,
1203 metricTensor,
1204 integralScaling);
1205 dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref.data()[kb];
1206 //get the metric tensor
1207 //cek todo use symmetry
1208 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
1209 //compute shape and solution information
1210 //shape
1211 ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],
1212 jacInv_ext,
1213 u_grad_trial_trace);
1214 //solution and gradients
1215 if (STABILIZATION_TYPE==STABILIZATION::TaylorGalerkinEV) //explicit
1216 {
1217 ck.valFromDOF(u_dof_old.data(),
1218 &u_l2g.data()[eN_nDOF_trial_element],
1219 &u_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],
1220 u_ext);
1221 ck.gradFromDOF(u_dof_old.data(),
1222 &u_l2g.data()[eN_nDOF_trial_element],
1223 u_grad_trial_trace,
1224 grad_u_ext);
1225 }
1226 else
1227 {
1228 ck.valFromDOF(u_dof.data(),
1229 &u_l2g.data()[eN_nDOF_trial_element],
1230 &u_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],
1231 u_ext);
1232 ck.gradFromDOF(u_dof.data(),
1233 &u_l2g.data()[eN_nDOF_trial_element],
1234 u_grad_trial_trace,
1235 grad_u_ext);
1236 }
1237 //precalculate test function products with integration weights
1238 for (int j=0;j<nDOF_trial_element;j++)
1239 {
1240 u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
1241 for (int I=0;I<nSpace;I++)
1242 u_grad_test_dS[j*nSpace+I] = u_grad_trial_trace[j*nSpace+I]*dS;//cek hack, using trial
1243
1244 }
1245 //
1246 //load the boundary values
1247 //
1248 bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+
1249 (1-isDOFBoundary_u.data()[ebNE_kb])*u_ext;
1250
1251
1252 //
1253 //
1254 //calculate the pde coefficients using the solution and the boundary values for the solution
1255 //
1256 double rho_out_ext=0.0,rho_out_bc=0.0;
1257 evaluateCoefficients(a_rowptr.data(),
1258 a_colind.data(),
1259 &ebqe_velocity_ext.data()[ebNE_kb_nSpace],
1260 alpha_L,
1261 alpha_T,
1262 Dm,
1263 ebqe_porosity.data()[ebNE_kb],
1264 rho_f,
1265 rho_s,
1266 u_ext,
1267 rho_out_ext,
1268 m_ext,
1269 dm_ext,
1270 f_ext,
1271 df_ext,
1272 a_ext,
1273 da_ext);
1274 ebqe_rho.data()[ebNE_kb] = rho_out_ext;
1275
1276 evaluateCoefficients(a_rowptr.data(),
1277 a_colind.data(),
1278 &ebqe_velocity_ext.data()[ebNE_kb_nSpace],
1279 alpha_L,
1280 alpha_T,
1281 Dm,
1282 ebqe_porosity.data()[ebNE_kb],
1283 rho_f,
1284 rho_s,
1285 bc_u_ext,
1286 rho_out_bc,
1287 bc_m_ext,
1288 bc_dm_ext,
1289 bc_f_ext,
1290 bc_df_ext,
1291 bc_a_ext,
1292 bc_da_ext);
1293 //moving mesh
1294 //
1295 double mesh_velocity[3];
1296 mesh_velocity[0] = xt_ext;
1297 mesh_velocity[1] = yt_ext;
1298 mesh_velocity[2] = zt_ext;
1299
1300 for (int I=0;I<nSpace;I++)
1301 {
1302 f_ext[I] -= MOVING_DOMAIN*m_ext*mesh_velocity[I];
1303 df_ext[I] -= MOVING_DOMAIN*dm_ext*mesh_velocity[I];
1304 bc_f_ext[I] -= MOVING_DOMAIN*bc_m_ext*mesh_velocity[I];
1305 bc_df_ext[I] -= MOVING_DOMAIN*bc_dm_ext*mesh_velocity[I];
1306 }
1307 //
1308 //calculate the numerical fluxes
1309 //
1310 exteriorNumericalAdvectiveFlux(isDOFBoundary_u.data()[ebNE_kb],
1311 isFluxBoundary_u.data()[ebNE_kb],
1312 forceStrongConditions,
1313 normal,
1314 ebqe_bc_flux_u_ext.data()[ebNE_kb],
1315 f_ext,
1316 bc_f_ext,
1317 df_ext,
1318 flux_adv_ext);
1319 exteriorNumericalDiffusiveFlux(a_rowptr.data(),
1320 a_colind.data(),
1321 isDOFBoundary_u.data()[ebNE_kb],
1322 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
1323 normal,
1324 a_ext,
1325 bc_u_ext,
1326 ebqe_bc_diffusiveFlux_u_ext.data()[ebNE_kb],
1327 a_ext,
1328 grad_u_ext,
1329 u_ext,
1330 ebqe_penalty_ext.data()[ebNE_kb],
1331 flux_diff_ext);
1332 flux_ext = flux_adv_ext + flux_diff_ext;
1333 double boundary_flow = 0.0;
1334 for (int I=0; I<nSpace; I++)
1335 boundary_flow += normal[I]*df_ext[I];
1336
1337 ebqe_flux.data()[ebNE_kb] = flux_ext;
1338 if (isDOFBoundary_u.data()[ebNE_kb] == 1)
1339 ebqe_u.data()[ebNE_kb] = bc_u_ext;
1340 else if (boundary_flow >= 0.0)
1341 ebqe_u.data()[ebNE_kb] = u_ext;
1342 else
1343 ebqe_u.data()[ebNE_kb] = bc_u_ext;
1344 if (STABILIZATION_TYPE==STABILIZATION::TaylorGalerkinEV)
1345 {
1346 if (stage == 1)
1347 flux_ext *= 1./3*dt;
1348 else
1349 flux_ext *= dt;
1350 }
1351
1352 //
1353 //update residuals
1354 //
1355 for (int i=0;i<nDOF_test_element;i++)
1356 {
1357 if (STABILIZATION_TYPE == STABILIZATION::Galerkin or STABILIZATION_TYPE == STABILIZATION::VMS or
1358 STABILIZATION_TYPE == STABILIZATION::TaylorGalerkinEV or STABILIZATION_TYPE == STABILIZATION::ImplicitEV)
1359 {
1360 elementResidual_u[i] += ck.ExteriorElementBoundaryFlux(flux_ext,u_test_dS[i])+
1361 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_u.data()[ebNE_kb],
1362 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
1363 eb_adjoint_sigma,
1364 u_ext,
1365 bc_u_ext,
1366 normal,
1367 a_rowptr.data(),
1368 a_colind.data(),
1369 a_ext,
1370 &u_grad_test_dS[i*nSpace]);
1371 }
1372 else if (STABILIZATION_TYPE == STABILIZATION::EntropyViscosity or
1373 STABILIZATION_TYPE == STABILIZATION::SmoothnessIndicator or
1374 STABILIZATION_TYPE == STABILIZATION::Kuzmin)
1375 {
1376 const double boundaryAdvectiveContribution =
1377 ck.ExteriorElementBoundaryFlux(flux_ext,u_test_dS[i]);
1378 const double boundaryDiffusiveContribution =
1379 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_u.data()[ebNE_kb],
1380 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
1381 eb_adjoint_sigma,
1382 u_ext,
1383 bc_u_ext,
1384 normal,
1385 a_rowptr.data(),
1386 a_colind.data(),
1387 a_ext,
1388 &u_grad_test_dS[i*nSpace]);
1389 const double boundaryResidualContribution =
1390 boundaryAdvectiveContribution + boundaryDiffusiveContribution;
1391 exteriorNumericalAdvectiveFluxDerivative(isDOFBoundary_u.data()[ebNE_kb],
1392 isFluxBoundary_u.data()[ebNE_kb],
1393 forceStrongConditions,
1394 normal,
1395 df_ext,
1396 dflux_u_u_ext);
1397
1398 if (dflux_u_u_ext> 0.0)
1399 {
1400 double boundaryTransportContribution = 0.0;
1401 for (int j=0;j<nDOF_trial_element;j++)
1402 {
1403 int ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j;
1404 double advJacobian_ext = 0.0, diffJacobian_ext = 0.0;
1405 exteriorNumericalAdvectiveFluxDerivative(isDOFBoundary_u.data()[ebNE_kb],
1406 isFluxBoundary_u.data()[ebNE_kb],
1407 forceStrongConditions,
1408 normal,
1409 df_ext,
1410 advJacobian_ext);
1411 exteriorNumericalDiffusiveFluxDerivative(isDOFBoundary_u.data()[ebNE_kb],
1412 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
1413 a_rowptr.data(),
1414 a_colind.data(),
1415 normal,
1416 a_ext,
1417 da_ext,
1418 grad_u_ext,
1419 &u_grad_trial_trace[j*nSpace],
1420 u_trial_trace_ref.data()[ebN_local_kb_j],
1421 ebqe_penalty_ext.data()[ebNE_kb],
1422 diffJacobian_ext);
1423 difffluxjacobian_ext = advJacobian_ext*u_trial_trace_ref.data()[ebN_local_kb_j]
1424 + diffJacobian_ext;
1425 const double localFluxTransportContribution =
1426 difffluxjacobian_ext*u_test_dS[i];
1427 fluxTransport[i][j] += localFluxTransportContribution;
1428 boundaryTransportContribution +=
1429 localFluxTransportContribution*u_dof_old.data()[u_l2g.data()[eN_nDOF_trial_element+j]];
1430 }
1431 elementResidual_u[i] += boundaryResidualContribution - boundaryTransportContribution;
1432 }
1433 else
1434 {
1435 elementResidual_u[i] += boundaryResidualContribution;
1436 // Upwind Nitsche penalty for advection-dominated Dirichlet
1437 // inflow. At inflow (boundary_flow = v.n < 0) the upwind
1438 // advective flux uses bc_u_ext and has zero derivative wrt
1439 // the interior u_ext -- so the existing IIPG penalty is the
1440 // only thing pulling u_ext toward bc_u_ext, and it scales as
1441 // max_a*penalty/h which becomes vanishingly small when D_m
1442 // is small. Add a Nitsche term that scales with |v.n| so
1443 // BC enforcement is independent of the diffusion coefficient.
1444 // Sign: at inflow with u_ext > bc_u_ext (overshoot), this
1445 // contributes positively to elementResidual_u[i] (= positive
1446 // boundary flux out of node i), which reduces mLow at the BC
1447 // DOF and pulls c back to bc_u_ext. Mass-conservative: the
1448 // term is integrated weakly with u_test_dS like the rest of
1449 // the boundary residual; sum over all faces telescopes the
1450 // weak Dirichlet to a consistent transport balance.
1451 if (isDOFBoundary_u.data()[ebNE_kb] == 1 && boundary_flow < 0.0)
1452 {
1453 const double upwind_penalty_rate = -boundary_flow; // |v.n|
1454 elementResidual_u[i] += upwind_penalty_rate
1455 * (u_ext - bc_u_ext)
1456 * u_test_dS[i];
1457 }
1458 }
1459 }
1460 }//i
1461 // local min/max at boundary.
1462 // At Dirichlet faces use the BC value (ebqe_bc_u_ext), not the
1463 // current solution trace (ebqe_u): the trace can drift off c_sat
1464 // under weak Nitsche enforcement, and feeding that drifted value
1465 // into min/max_u_bc pollutes the FCT bound at every interior
1466 // neighbor. Using the BC value keeps the bound tight at c_sat
1467 // (combined with bc_mask in the FCT step, this gives the bounded
1468 // + mass-conservative recipe that mphase_co2 uses).
1469 const double u_for_bound =
1470 isDOFBoundary_u.data()[ebNE_kb]
1471 ? ebqe_bc_u_ext.data()[ebNE_kb]
1472 : ebqe_u.data()[ebNE_kb];
1473 min_u_bc_local = fmin(u_for_bound, min_u_bc_local);
1474 max_u_bc_local = fmax(u_for_bound, max_u_bc_local);
1475 }//kb
1476 //
1477 //update the element and global residual storage
1478 //
1479 for (int i=0;i<nDOF_test_element;i++)
1480 {
1481 int eN_i = eN*nDOF_test_element+i;
1482 int gi = offset_u+stride_u*u_l2g.data()[eN_i]; //global i-th index
1483 if (STABILIZATION_TYPE==STABILIZATION::EntropyViscosity or STABILIZATION_TYPE==STABILIZATION::SmoothnessIndicator or STABILIZATION_TYPE==STABILIZATION::Kuzmin)
1484 {
1485 globalResidual.data()[gi] += dt*elementResidual_u[i];
1486 boundary_integral[gi] += elementResidual_u[i];
1487 min_u_bc[gi] = fmin(min_u_bc_local,min_u_bc[gi]);
1488 max_u_bc[gi] = fmax(max_u_bc_local,max_u_bc[gi]);
1489 for (int j=0;j<nDOF_trial_element;j++)
1490 {
1491 int ebN_i_j = ebN*4*nDOF_test_X_trial_element + i*nDOF_trial_element + j;
1492 TransportMatrix[csrRowIndeces_CellLoops.data()[eN_i] + csrColumnOffsets_eb_CellLoops.data()[ebN_i_j]]
1493 += fluxTransport[i][j];
1494 TransposeTransportMatrix[csrRowIndeces_CellLoops.data()[eN_i] + csrColumnOffsets_eb_CellLoops.data()[ebN_i_j]]
1495 += fluxTransport[j][i];
1496 }//j
1497 }
1498 else if (STABILIZATION_TYPE==STABILIZATION::ImplicitEV)
1499 {
1500 // Implicit boundary: stash the consistent boundary residual in
1501 // boundary_integral (the edge loop adds + boundary_integral[i]
1502 // and OVERWRITES globalResidual[i], so we must NOT add here).
1503 // The boundary flux Jacobian is assembled in calculateJacobian.
1504 boundary_integral[gi] += elementResidual_u[i];
1505 // FCT bounds at the boundary (used by the explicit FCT
1506 // post-step when FCT=True; harmless otherwise).
1507 min_u_bc[gi] = fmin(min_u_bc_local,min_u_bc[gi]);
1508 max_u_bc[gi] = fmax(max_u_bc_local,max_u_bc[gi]);
1509 }
1510 else
1511 {
1512 globalResidual.data()[offset_u+stride_u*r_l2g.data()[eN_i]] += elementResidual_u[i];
1513 }
1514 }//i
1515 }//ebNE
1516 if (STABILIZATION_TYPE==STABILIZATION::TaylorGalerkinEV)
1517 {
1518 meanEntropy /= meanOmega;
1519 double norm_factor = fmax(fabs(maxEntropy - meanEntropy), fabs(meanEntropy-minEntropy));
1520 for(int eN=0;eN<nElements_global;eN++)
1521 {
1522 double hK=elementDiameter.data()[eN]/degree_polynomial;
1523 double linear_viscosity = cMax*hK*maxVel[eN];
1524 double entropy_viscosity = cE*hK*hK*maxEntRes[eN]/norm_factor;
1525 for (int k=0;k<nQuadraturePoints_element;k++)
1526 {
1527 int eN_k = eN*nQuadraturePoints_element+k;
1528 q_numDiff_u.data()[eN_k] = fmin(linear_viscosity,entropy_viscosity);
1529 }
1530 }
1531 }
1532 //edge based stabilization
1533 else if (STABILIZATION_TYPE==STABILIZATION::EntropyViscosity or
1534 STABILIZATION_TYPE==STABILIZATION::SmoothnessIndicator or
1535 STABILIZATION_TYPE==STABILIZATION::Kuzmin)
1536 {
1538 // COMPUTE SMOOTHNESS INDICATOR and NORMALIZE ENTROPY RESIDUAL //
1540 // NOTE: see NCLS.h for a different but equivalent implementation of this.
1541 //cek todo: can these loops over numDOFs be collapsed?
1542 int ij = 0;
1543 for (int i=0; i<numDOFs; i++)
1544 {
1545 double etaMaxi, etaMini;
1546 if (STABILIZATION_TYPE==STABILIZATION::EntropyViscosity) //EV
1547 {
1548 // For eta min and max
1549 etaMaxi = fabs(eta[i]);
1550 etaMini = fabs(eta[i]);
1551 }
1552 // for smoothness indicator //
1553 double alpha_numerator = 0., alpha_denominator = 0.;
1554 for (int offset=csrRowIndeces_DofLoops.data()[i]; offset<csrRowIndeces_DofLoops.data()[i+1]; offset++)
1555 { // First loop in j (sparsity pattern)
1556 int j = csrColumnOffsets_DofLoops.data()[offset];
1557 if (STABILIZATION_TYPE==STABILIZATION::EntropyViscosity) //EV Stabilization
1558 {
1559 // COMPUTE ETA MIN AND ETA MAX //
1560 etaMaxi = fmax(etaMaxi,fabs(eta[j]));
1561 etaMini = fmin(etaMini,fabs(eta[j]));
1562 }
1563 // Sense smoothness on the projected nodal conservative variable
1564 // m = theta * rho(u) * u so the edge indicator is aligned with
1565 // the variable-density storage used by the PDE residual.
1566 const double mi = m_dof[i];
1567 const double mj = m_dof[j];
1568 alpha_numerator += mi - mj;
1569 alpha_denominator += fabs(mi - mj);
1570 //update ij
1571 ij+=1;
1572 }
1573 if (STABILIZATION_TYPE==STABILIZATION::EntropyViscosity) //EV Stab
1574 {
1575 // Normalize entropy residual
1576 global_entropy_residual[i] *= etaMini == etaMaxi ? 0. : 2*cE/(etaMaxi-etaMini);
1577 quantDOFs[i] = fabs(global_entropy_residual[i]);
1578 }
1579
1580 double alphai = alpha_numerator/(alpha_denominator+1E-15);
1581 quantDOFs[i] = alphai;
1582
1583
1584 if (POWER_SMOOTHNESS_INDICATOR==0)
1585 psi[i] = 1.0;
1586 else
1587 psi[i] = std::pow(alphai,POWER_SMOOTHNESS_INDICATOR); //NOTE: they use alpha^2 in the paper
1588 }
1590 // ** LOOP IN DOFs FOR EDGE BASED TERMS ** //
1592 ij=0;
1593 for (int i=0; i<numDOFs; i++)
1594 {
1595 const double mi_mass = m_dof[i];
1596 const double theta_i = fmax(theta_dof_proj[i], 1.0e-14);
1597 const double ui_mass = inversevaluateCoefficients(mi_mass, theta_i, rho_f, rho_s);
1598 const double rho_i = fmax(rho_dof_proj[i], rho_f);
1599 const double drho_du = rho_s - rho_f;
1600 const double dmdu_i = theta_i * (rho_i + ui_mass*drho_du);
1601 double ith_dissipative_term_mass = 0;
1602 double ith_low_order_dissipative_term_mass = 0;
1603 double ith_flux_term_mass = 0;
1604 // Row sum of the high-order graph viscosity, dLii = -sum_{j!=i} dLij.
1605 // Feeds edge_based_cfl below.
1606 double dLii = 0.;
1607
1608 // loop over the sparsity pattern of the i-th DOF
1609 for (int offset=csrRowIndeces_DofLoops.data()[i]; offset<csrRowIndeces_DofLoops.data()[i+1]; offset++)
1610 {
1611 int j = csrColumnOffsets_DofLoops.data()[offset];
1612 const double mj_mass = m_dof[j];
1613 const double theta_j = fmax(theta_dof_proj[j], 1.0e-14);
1614 const double uj_mass = inversevaluateCoefficients(mj_mass, theta_j, rho_f, rho_s);
1615 double dLowij, dLij, dEVij, dHij;
1616
1617 ith_flux_term_mass += (TransportMatrix[ij] + DiffusionMatrix[ij])*uj_mass;
1618
1619 if (i != j)
1620 {
1621 double solij = 0.5*(ui_mass+uj_mass);
1622 double Compij = cK*fmax(solij*(1.0-solij),0.0)/(fabs(ui_mass-uj_mass)+1E-14);
1623 dLowij = fmax(fabs(TransportMatrix[ij]),fabs(TransposeTransportMatrix[ij]));
1624
1625 dLij = dLowij*fmax(psi[i],psi[j]); // Approach by JLG & BP
1626
1627 if (STABILIZATION_TYPE==STABILIZATION::EntropyViscosity) //EV Stab
1628 {
1629 // high-order (entropy viscosity) dissipative operator
1630 dEVij = fmax(fabs(global_entropy_residual[i]),fabs(global_entropy_residual[j]));
1631 dHij = fmin(dLowij,dEVij) * fmax(1.0-Compij,0.0); // artificial compression
1632 }
1633 else // smoothness based indicator
1634 {
1635 dHij = dLij * fmax(1.0-Compij,0.0); // artificial compression
1636 }
1637 // Dissipative terms. The DIFFERENCES are taken in c, not in
1638 // m: dLow comes from the transport matrix, which already
1639 // carries rho (df/du = (rho + u*drho/du)*v), so dLow*(c_j-c_i)
1640 // is a mass flux, dimensionally consistent with
1641 // ith_flux_term_mass. Differencing m instead multiplies the
1642 // dissipation by a spurious extra theta*rho (~1e3 in SI).
1643 // Note this pair IS already the monotone upwind operator:
1644 // sum_j (T+D)_ij c_j - sum_{j!=i} dLow_ij (c_j - c_i)
1645 // = -sum_{j!=i} (dLow_ij - T_ij - D_ij)(c_j - c_i)
1646 // with dLow_ij - T_ij - D_ij >= 0 because dLow_ij >= |T_ij|,
1647 // i.e. an M-matrix. It must NOT be replaced by
1648 // -sum max(0,-T_ij)(c_j - c_i): for the skew-symmetric
1649 // advection matrix that is exactly HALF the upwind flux.
1650 ith_dissipative_term_mass += dHij*(uj_mass-ui_mass);
1651 ith_low_order_dissipative_term_mass += dLowij*(uj_mass-ui_mass);
1652 //dHij - dLij. This matrix is needed during FCT step
1653 dt_times_dH_minus_dL[ij] = dt*(dHij - dLowij);
1654
1655 dLii -= dLij;
1656 dLow[ij] = dLowij;
1657
1658 }
1659 else //i==j
1660 {
1661 // NOTE: this is incorrect. Indeed, dLii = -sum_{j!=i}(dLij) and similarly for dCii.
1662 // However, it is irrelevant since during the FCT step we do (dL-dC)*(solnj-solni)
1663 dt_times_dH_minus_dL[ij]=0;
1664 dLow[ij]=0;
1665 }
1666 //update ij
1667 ij+=1;
1668 }
1669 double mi = ML.data()[i];
1670 // boundary_integral was assembled from m-space fluxes (rho*v*u and
1671 // a=theta*rho*Disp); no dmdu_i lifting is needed.
1672 const double boundary_integral_mass = boundary_integral[i];
1673 // compute edge_based_cfl
1674 // 2|dLii|/mi is the constant-density edge CFL. The low-order step
1675 // advances m = theta*rho(c)*c, so recovering c divides the mass
1676 // change by the storage Jacobian dm/dc = dmdu_i; without that
1677 // factor the stable dt is over-predicted in low-water-content (high
1678 // gas saturation) zones where dmdu_i -> 0, and uLow overshoots out
1679 // of [0, c_sat]. dmdu_i = 1 at rho=theta=1, so this reduces to the
1680 // original 2|dLii|/mi.
1681 edge_based_cfl.data()[i] = 2.*fabs(dLii)/(mi * fmax(dmdu_i, 1.0e-14));
1682
1683 // Stage 3 kinetic dissolution source at node i (mass-rate form):
1684 // R_diss_i = theta_w_i * rho_w(u_i) * k_d * S_n_i * (c_sat - u_i)
1685 // Added directly to mLow_i / mHigh_i since the lumped-mass time
1686 // discretization gives dm/dt = R_diss with no further scaling.
1687 // Sign of (c_sat - u_i) ensures R_diss > 0 when undersaturated
1688 // (dissolution) and < 0 when supersaturated (exsolution).
1689 // The S_w (interfacial-area) factor a_gw ~ S_n*S_w was dropped:
1690 // it makes R_diss vanish at S_n = 1, so a gas pool stops
1691 // dissolving once it consolidates. theta_w is kept so the CO2
1692 // still goes into the brine that exists (theta_w > 0 down to
1693 // residual S_wr). MUST stay consistent with mphase_co2.h.
1694 const double S_n_i = Sn_dof.data()[i];
1695 const double rho_w_i = rho_f * (1.0 + ((rho_s - rho_f)/rho_f) * ui_mass);
1696 const double R_diss_i = theta_i * rho_w_i * k_d * S_n_i
1697 * (c_sat - ui_mass);
1698
1699 const double mLow_i = mi_mass - dt/mi*(ith_flux_term_mass
1700 + boundary_integral_mass
1701 - ith_low_order_dissipative_term_mass)
1702 + dt * R_diss_i;
1703 uLow[i] = inversevaluateCoefficients(mLow_i, theta_i, rho_f, rho_s);
1704
1705 // update residual
1706 if (LUMPED_MASS_MATRIX==1)
1707 {
1708 const double mHigh_i = mi_mass - dt/mi*(ith_flux_term_mass
1709 + boundary_integral_mass
1710 - ith_dissipative_term_mass)
1711 + dt * R_diss_i;
1712 globalResidual.data()[i] = inversevaluateCoefficients(mHigh_i, theta_i, rho_f, rho_s);
1713 }
1714 else
1715 globalResidual.data()[i] += dt*(ith_flux_term_mass - ith_dissipative_term_mass - R_diss_i);//cek todo: shouldn't this have boundaryIntegral?
1716 }//i
1717 }//edge-based
1718 else if (STABILIZATION_TYPE==STABILIZATION::ImplicitEV)
1719 {
1721 // IMPLICIT edge-based scheme (backward Euler, low-order graph //
1722 // dissipation, no FCT). Mirrors Richards STABILIZATION_TYPE==2 //
1723 // but evaluates the advection/diffusion at the CURRENT Newton //
1724 // iterate (u_dof) so they contribute to the Jacobian (assembled //
1725 // separately in calculateJacobian) and the scheme is NOT //
1726 // CFL-limited. Per-DOF residual: //
1727 // //
1728 // R_i = ML_i*(m_i^{n+1} - m_i^n)/dt //
1729 // + sum_j (T_ij + D_ij)*c_j (consistent flux) //
1730 // - sum_{j!=i} dLow_ij*(m_j^{n+1}-m_i^{n+1}) (graph diss) //
1731 // + boundary_integral_i (0 for closed BCs) //
1732 // - ML_i*R_diss_i (dissolution) //
1733 // //
1734 // m^{n+1}_i = theta_i*rho(c_i)*c_i (nodal, porosity lagged to //
1735 // theta_dof_proj); m^n_i = m_dof[i] (old projected mass). //
1736 // The advection coefficient in T uses the CURRENT-time velocity //
1737 // (q_v, refreshed from mphase_co2 before TADR solves) and the //
1738 // current iterate u -- fully implicit backward Euler. //
1739 // //
1740 // FLUX IS UPWINDED: instead of the central sum_j(T_ij+D_ij)c_j //
1741 // plus symmetric dissipation, the advective/diffusive flux is the//
1742 // first-order upwind low-order operator (the same one TADR's //
1743 // explicit branch builds for uLow), evaluated at the current c: //
1744 // //
1745 // F_i = sum_{j!=i} -a_ij*(c_j - c_i), //
1746 // a_ij = max(0,-T_ij)*(rho_up/rho_f) + max(0,-D_ij) >= 0, //
1747 // rho_up = rho(c_i) if j->i is inflow (-T_ij*(c_j-c_i)<=0) //
1748 // else rho(c_j). //
1749 // //
1750 // a_ij>=0 makes the spatial operator an M-matrix (monotone, DMP);//
1751 // diffusion (D symmetric, row-sum 0) is reproduced exactly since //
1752 // -max(0,-D_ij)(c_j-c_i)=D_ij(c_j-c_i). Conservative by global //
1753 // telescoping of the edge fluxes. a_ij (frozen) is stashed in //
1754 // dLow[ij] for calculateJacobian. m^{n+1}_i=theta_i*rho(c_i)*c_i //
1755 // (nodal, porosity lagged); m^n_i=m_dof[i]. boundary_integral //
1756 // stays 0 (closed system); implicit boundary flux is a TODO. //
1758 const double drho_du = rho_s - rho_f;
1759 // nodal current conservative variable m^{n+1}_i = theta_i*rho(c_i)*c_i
1760 std::valarray<double> m_new(numDOFs);
1761 for (int i=0; i<numDOFs; i++)
1762 {
1763 const double ci = u_dof.data()[i];
1764 const double theta_i = fmax(theta_dof_proj[i], 1.0e-14);
1765 const double rho_ci = rho_f*(1.0 + (drho_du/rho_f)*ci);
1766 m_new[i] = theta_i*rho_ci*ci;
1767 }
1768 int ij=0;
1769 for (int i=0; i<numDOFs; i++)
1770 {
1771 const double ci = u_dof.data()[i];
1772 const double theta_i = fmax(theta_dof_proj[i], 1.0e-14);
1773 const double rho_ci = rho_f*(1.0 + (drho_du/rho_f)*ci);
1774 const double mi_new = m_new[i];
1775 const double mn_i = m_dof[i];
1776 const double MLi = ML.data()[i];
1777
1778 double ith_upwind_flux_term_mass = 0.0;
1779 for (int offset=csrRowIndeces_DofLoops.data()[i]; offset<csrRowIndeces_DofLoops.data()[i+1]; offset++)
1780 {
1781 int j = csrColumnOffsets_DofLoops.data()[offset];
1782 if (i != j)
1783 {
1784 const double cj = u_dof.data()[j];
1785 const double rho_cj = rho_f*(1.0 + (drho_du/rho_f)*cj);
1786 const double T_ij = TransportMatrix[ij];
1787 const double D_ij = DiffusionMatrix[ij];
1788 const double delta_c = cj - ci;
1789 const double T_neg = fmax(0.0, -T_ij);
1790 const double D_neg = fmax(0.0, -D_ij);
1791 const double rho_up = (-T_ij*delta_c <= 0.0) ? rho_ci : rho_cj;
1792 // frozen edge coefficient a_ij >= 0 (advection upwind + diffusion)
1793 const double a_ij = T_neg*(rho_up/rho_f) + D_neg;
1794 ith_upwind_flux_term_mass += -a_ij*delta_c;
1795 // dLow stores the SYMMETRIC graph viscosity used by the
1796 // explicit FCT post-step (Kuzmin antidiffusive flux, which
1797 // must be antisymmetric per edge -> needs symmetric dLow,
1798 // NOT the directional a_ij). calculateJacobian recomputes
1799 // a_ij directly from T,D, so it no longer reads dLow.
1800 dLow.data()[ij] = fmax(fabs(T_ij), fabs(TransposeTransportMatrix[ij]));
1801 // High-order (smoothness-compressed) graph dissipation for the
1802 // FCT post-step, mirroring the explicit EntropyViscosity branch
1803 // (dHij = dLowij*(1-Compij)). Store dt*(dH - dLow) <= 0 so the
1804 // FCTStep removes ONLY the EXCESS low-order dissipation
1805 // (dLow - dH) and KEEPS dH. Removing the full dLow (as the
1806 // explicit Kuzmin flux does) makes the antidiffusion target the
1807 // dissipation-free Galerkin solution, which over-sharpens the
1808 // Dirichlet source front back to the initial condition -> the
1809 // observed FCT=True freeze. globalResidual is NOT touched here,
1810 // so calculateJacobian needs no change.
1811 {
1812 // Standard Kuzmin implicit FEM-FCT high-order target: ZERO
1813 // artificial dissipation (d^H = 0). The full low-order graph
1814 // viscosity dLow is antidiffused; the Zalesak limiter
1815 // (min/max_u_bc bounds) provides boundedness. The earlier
1816 // dH = dLow*(1-Comp) "smoothness compression" throttled the
1817 // antidiffusion so the consistent-mass term dominated and the
1818 // front stayed diffuse (verified in the 1D FCT numpy replica);
1819 // removing it recovers the textbook antidiffusive flux
1820 // F_ij = M~_ij[(m_j^H-m_j^n)-(m_i^H-m_i^n)]
1821 // + dt (d^H-d^L)_ij (m_j^n-m_i^n), d^H=0.
1822 const double dHij = 0.0;
1823 dt_times_dH_minus_dL.data()[ij] = dt*(dHij - dLow.data()[ij]);
1824 }
1825 }
1826 else
1827 {
1828 dLow.data()[ij] = 0.0;
1829 dt_times_dH_minus_dL.data()[ij] = 0.0;
1830 }
1831 ij += 1;
1832 }
1833 // Stage-3 kinetic dissolution source (mass-rate form); MUST match
1834 // mphase_co2.h and the explicit edge-based branch above.
1835 const double S_n_i = Sn_dof.data()[i];
1836 const double rho_w_i = rho_ci; // rho_w(c_i) = rho_f*(1 + eps*c_i)
1837 const double R_diss_i = theta_i * rho_w_i * k_d * S_n_i * (c_sat - ci);
1838
1839 globalResidual.data()[i] = MLi*(mi_new - mn_i)/dt
1840 + ith_upwind_flux_term_mass
1841 + boundary_integral[i]
1842 - MLi*R_diss_i;
1843 }//i
1844 }//implicit edge-based
1845 }
1846
1847 void invert(arguments_dict& args)
1848 {
1849 int numDOFs = args.scalar<int>("numDOFs");
1850 xt::pyarray<double>& mIn = args.array<double>("mIn");
1851 xt::pyarray<double>& uOut = args.array<double>("uOut");
1852 xt::pyarray<double>& nodal_porosity = args.array<double>("nodal_porosity");
1853 const double rho_f = args.scalar<double>("rho_f");
1854 const double rho_s = args.scalar<double>("rho_s");
1855 for (int i=0; i<numDOFs; i++)
1856 uOut.data()[i] = inversevaluateCoefficients(mIn.data()[i], nodal_porosity.data()[i], rho_f, rho_s);
1857 }
1858
1859
1860 void calculateJacobian(arguments_dict& args)
1861 {
1862 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1863 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1864 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1865 xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
1866 double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
1867 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1868 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1869 xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1870 xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1871 xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1872 xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
1873 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
1874 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
1875 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
1876 xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
1877 xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
1878 xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
1879 xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
1880 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
1881 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
1882 int nElements_global = args.scalar<int>("nElements_global");
1883 double useMetrics = args.scalar<double>("useMetrics");
1884 double alphaBDF = args.scalar<double>("alphaBDF");
1885 int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
1886 double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
1887 xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1888 xt::pyarray<int>& r_l2g = args.array<int>("r_l2g");
1889 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1890 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1891 xt::pyarray<double>& velocity = args.array<double>("velocity");
1892 xt::pyarray<double>& q_porosity = args.array<double>("q_porosity");
1893 xt::pyarray<double>& q_rho = args.array<double>("q_rho");
1894 xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
1895 xt::pyarray<double>& cfl = args.array<double>("cfl");
1896 xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
1897 xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
1898 xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
1899 xt::pyarray<double>& globalJacobian = args.array<double>("globalJacobian");
1900 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
1901 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
1902 xt::pyarray<int>& elementBoundaryMaterialTypes = args.array<int>("elementBoundaryMaterialTypes");
1903 xt::pyarray<int>& isExteriorBoundaryPhysical = args.array<int>("isExteriorBoundaryPhysical");
1904 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
1905 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
1906 xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
1907 xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
1908 xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
1909 xt::pyarray<int>& isFluxBoundary_u = args.array<int>("isFluxBoundary_u");
1910 xt::pyarray<double>& ebqe_bc_flux_u_ext = args.array<double>("ebqe_bc_flux_u_ext");
1911 xt::pyarray<double>& ebqe_porosity = args.array<double>("ebqe_porosity");
1912 xt::pyarray<double>& ebqe_rho = args.array<double>("ebqe_rho");
1913 xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
1914 STABILIZATION STABILIZATION_TYPE = static_cast<STABILIZATION>(args.scalar<int>("STABILIZATION_TYPE"));
1915// ENTROPY ENTROPY_TYPE = static_cast<ENTROPY>(args.scalar<int>("ENTROPY_TYPE"));
1916// STABILIZATION STABILIZATION_TYPE{args.scalar<int>("STABILIZATION_TYPE")};
1917 double physicalDiffusion = args.scalar<double>("physicalDiffusion");
1918 const double alpha_L = args.scalar<double>("alpha_L");
1919 const double alpha_T = args.scalar<double>("alpha_T");
1920 const double Dm = args.scalar<double>("Dm");
1921 int forceStrongConditions = args.scalar<int>("forceStrongConditions");
1922 // const int dispersion_type_int = args.scalar<int>("dispersion_type");
1923 // const double theta_s = args.scalar<double>("theta_s");
1924 // const double theta_r = args.scalar<double>("theta_r");
1925 // const double power_law_exponent = args.scalar<double>("power_law_exponent");
1926 // const double velocity_exponent = args.scalar<double>("velocity_exponent");
1927 const double rho_f = args.scalar<double>("rho_f");
1928 const double rho_s = args.scalar<double>("rho_s");
1929
1930 // DISPERSION DISPERSION_TYPE = static_cast<DISPERSION>(dispersion_type_int);
1931
1932
1933 double Ct_sge = 4.0;
1934
1935
1936
1938 xt::pyarray<int>& a_rowptr = args.array<int>("a_rowptr");
1939 xt::pyarray<int>& a_colind = args.array<int>("a_colind");
1940 //xt::pyarray<double>& D = args.array<double>("D");
1942 xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.array<int>("isDiffusiveFluxBoundary_u");
1943 xt::pyarray<double>& ebqe_penalty_ext = args.array<double>("ebqe_penalty_ext");
1944
1946 // ImplicitEV (STAB=5): assemble the implicit edge-based Jacobian in //
1947 // THIS function (unlike Richards, which builds it inside the //
1948 // residual). It reads the transport/diffusion matrices and nodal //
1949 // projections (theta_dof_proj) that the immediately preceding //
1950 // calculateResidual call filled -- Proteus Newton evaluates F(u) //
1951 // then J(u) at the same iterate -- plus the frozen per-edge upwind //
1952 // coefficient a_ij the residual stashed in dLow[ij]. Single-comp P1 //
1953 // => the DOF-loop CSR coincides with the matrix CSR, so globalJacobian//
1954 // is indexed by the same running offset 'ij' as the residual //
1955 // (identical to the Richards STAB==2 convention). //
1956 // //
1957 // Differentiating the upwind flux F_i = sum_{j!=i} -a_ij*(c_j-c_i) //
1958 // with a_ij frozen gives a clean M-matrix: //
1959 // J_ij = -a_ij (j != i, <= 0) //
1960 // J_ii = ML_i*dmdu_i/dt + sum_{j!=i} a_ij - ML_i*dRdiss_i/dc_i //
1961 // dmdu_i = theta_i*(rho(c_i)+c_i*drho_du), rho(c)=rho_f(1+eps*c). //
1962 // a_ij carries the upwind density/diffusion weighting; freezing it //
1963 // (and the O(eps) rho_up dependence) is the only inexactness //
1964 // (standard inexact-Newton, as Richards freezes Kr/psi). //
1966 if (STABILIZATION_TYPE==STABILIZATION::ImplicitEV)
1967 {
1968 const double dt = args.scalar<double>("dt");
1969 const int numDOFs = args.scalar<int>("numDOFs");
1970 xt::pyarray<double>& ML = args.array<double>("ML");
1971 xt::pyarray<double>& dLow = args.array<double>("dLow");
1972 xt::pyarray<double>& Sn_dof = args.array<double>("Sn_dof");
1973 const double k_d = args.scalar<double>("k_d");
1974 const double c_sat = args.scalar<double>("c_sat");
1975 xt::pyarray<int>& csrRowIndeces_DofLoops = args.array<int>("csrRowIndeces_DofLoops");
1976 xt::pyarray<int>& csrColumnOffsets_DofLoops = args.array<int>("csrColumnOffsets_DofLoops");
1977 const double drho_du = rho_s - rho_f;
1978
1979 // nodal storage Jacobian dmdu_i = theta_i*(rho(c_i) + c_i*drho_du)
1980 std::valarray<double> dmdu(numDOFs);
1981 for (int i=0; i<numDOFs; i++)
1982 {
1983 const double ci = u_dof.data()[i];
1984 const double theta_i = fmax(theta_dof_proj[i], 1.0e-14);
1985 const double rho_ci = rho_f*(1.0 + (drho_du/rho_f)*ci);
1986 dmdu[i] = theta_i*(rho_ci + ci*drho_du);
1987 }
1988 int ij=0;
1989 for (int i=0; i<numDOFs; i++)
1990 {
1991 const double ci = u_dof.data()[i];
1992 const double theta_i = fmax(theta_dof_proj[i], 1.0e-14);
1993 const double rho_ci = rho_f*(1.0 + (drho_du/rho_f)*ci);
1994 const double MLi = ML.data()[i];
1995 const double S_n_i = Sn_dof.data()[i];
1996 // dR_diss_i/dc_i = theta_i*k_d*S_n_i*[ drho_du*(c_sat-c_i) - rho(c_i) ]
1997 const double dRdiss_dci = theta_i*k_d*S_n_i*( drho_du*(c_sat - ci) - rho_ci );
1998
1999 int diag_ij = -1;
2000 double sum_a = 0.0;
2001 for (int offset=csrRowIndeces_DofLoops.data()[i]; offset<csrRowIndeces_DofLoops.data()[i+1]; offset++)
2002 {
2003 int j = csrColumnOffsets_DofLoops.data()[offset];
2004 if (i != j)
2005 {
2006 // off-diagonal upwind flux: J_ij = -a_ij. Recompute a_ij
2007 // from the (frozen) transport/diffusion matrices and the
2008 // current iterate -- dLow now stores the symmetric FCT
2009 // viscosity, not a_ij. Matches the residual exactly.
2010 const double T_ij = TransportMatrix[ij];
2011 const double D_ij = DiffusionMatrix[ij];
2012 const double cj = u_dof.data()[j];
2013 const double rho_cj= rho_f*(1.0 + (drho_du/rho_f)*cj);
2014 const double T_neg = fmax(0.0, -T_ij);
2015 const double D_neg = fmax(0.0, -D_ij);
2016 const double rho_up= (-T_ij*(cj - ci) <= 0.0) ? rho_ci : rho_cj;
2017 const double a_ij = T_neg*(rho_up/rho_f) + D_neg;
2018 globalJacobian.data()[ij] += -a_ij;
2019 sum_a += a_ij;
2020 }
2021 else
2022 diag_ij = ij;
2023 ij += 1;
2024 }
2025 // diagonal: storage + upwind self term (sum a_ij) - dissolution
2026 globalJacobian.data()[diag_ij] += MLi*dmdu[i]/dt
2027 + sum_a
2028 - MLi*dRdiss_dci;
2029 }//i
2030
2031 //
2032 // Implicit boundary flux Jacobian: d(boundary flux)/du -> globalJacobian
2033 // so the boundary is implicit (like Richards), consistent with the
2034 // boundary residual the calculateResidual ImplicitEV path stashes in
2035 // boundary_integral. fluxJacobian accumulates over the face
2036 // quadrature and is loaded ONCE per face (the residual side likewise
2037 // distributes once). For closed BCs (zero adv/diff flux, no
2038 // Dirichlet) every term here is identically 0.
2039 //
2040 for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
2041 {
2042 int ebN = exteriorElementBoundariesArray.data()[ebNE];
2043 const int eN_out = elementBoundaryElementsArray.data()[ebN*2+1];
2044 const int ebFlag = elementBoundaryMaterialTypes.data()[ebN];
2045 if (ebFlag <= 0 || isExteriorBoundaryPhysical.data()[ebNE] == 0 || eN_out >= 0)
2046 continue;
2047 int eN = elementBoundaryElementsArray.data()[ebN*2+0],
2048 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
2049 eN_nDOF_trial_element = eN*nDOF_trial_element;
2050 double fluxJacobian_u_u[nDOF_test_element][nDOF_trial_element];
2051 for (int i=0;i<nDOF_test_element;i++)
2052 for (int j=0;j<nDOF_trial_element;j++)
2053 fluxJacobian_u_u[i][j]=0.0;
2054 for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
2055 {
2056 int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
2057 ebNE_kb_nSpace = ebNE_kb*nSpace,
2058 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
2059 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
2060 double u_ext=0.0, grad_u_ext[nSpace], m_ext=0.0, dm_ext=0.0,
2061 f_ext[nSpace], df_ext[nSpace], a_ext[nnz], da_ext[nnz],
2062 bc_a_ext[nnz], bc_da_ext[nnz], difffluxjacobian_ext=0.0,
2063 bc_u_ext=0.0, bc_m_ext=0.0, bc_dm_ext=0.0,
2064 bc_f_ext[nSpace], bc_df_ext[nSpace],
2065 jac_ext[nSpace*nSpace], jacDet_ext, jacInv_ext[nSpace*nSpace],
2066 boundaryJac[nSpace*(nSpace-1)], metricTensor[(nSpace-1)*(nSpace-1)],
2067 metricTensorDetSqrt, dS, u_test_dS[nDOF_test_element],
2068 u_grad_trial_trace[nDOF_trial_element*nSpace],
2069 u_grad_test_dS[nDOF_trial_element*nSpace], normal[nSpace],
2070 x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
2071 G[nSpace*nSpace],G_dd_G,tr_G;
2072 ck.calculateMapping_elementBoundary(eN,ebN_local,kb,ebN_local_kb,mesh_dof.data(),mesh_l2g.data(),mesh_trial_trace_ref.data(),mesh_grad_trial_trace_ref.data(),boundaryJac_ref.data(),jac_ext,jacDet_ext,jacInv_ext,boundaryJac,metricTensor,metricTensorDetSqrt,normal_ref.data(),normal,x_ext,y_ext,z_ext);
2073 ck.calculateMappingVelocity_elementBoundary(eN,ebN_local,kb,ebN_local_kb,mesh_velocity_dof.data(),mesh_l2g.data(),mesh_trial_trace_ref.data(),xt_ext,yt_ext,zt_ext,normal,boundaryJac,metricTensor,integralScaling);
2074 dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref.data()[kb];
2075 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
2076 ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,u_grad_trial_trace);
2077 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);
2078 ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial_trace,grad_u_ext);
2079 for (int j=0;j<nDOF_trial_element;j++)
2080 {
2081 u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
2082 for (int I=0;I<nSpace;I++)
2083 u_grad_test_dS[j*nSpace+I]= u_grad_trial_trace[j*nSpace+I]*dS;
2084 }
2085 bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+(1-isDOFBoundary_u.data()[ebNE_kb])*u_ext;
2086 double rho_out_ext=0.0,rho_out_bc=0.0;
2087 evaluateCoefficients(a_rowptr.data(),a_colind.data(),&ebqe_velocity_ext.data()[ebNE_kb_nSpace],alpha_L,alpha_T,Dm,ebqe_porosity.data()[ebNE_kb],rho_f,rho_s,u_ext,rho_out_ext,m_ext,dm_ext,f_ext,df_ext,a_ext,da_ext);
2088 evaluateCoefficients(a_rowptr.data(),a_colind.data(),&ebqe_velocity_ext.data()[ebNE_kb_nSpace],alpha_L,alpha_T,Dm,ebqe_porosity.data()[ebNE_kb],rho_f,rho_s,bc_u_ext,rho_out_bc,bc_m_ext,bc_dm_ext,bc_f_ext,bc_df_ext,bc_a_ext,bc_da_ext);
2089 double mesh_velocity[3]; mesh_velocity[0]=xt_ext; mesh_velocity[1]=yt_ext; mesh_velocity[2]=zt_ext;
2090 for (int I=0;I<nSpace;I++)
2091 {
2092 f_ext[I] -= MOVING_DOMAIN*m_ext*mesh_velocity[I];
2093 df_ext[I] -= MOVING_DOMAIN*dm_ext*mesh_velocity[I];
2094 bc_f_ext[I] -= MOVING_DOMAIN*bc_m_ext*mesh_velocity[I];
2095 bc_df_ext[I] -= MOVING_DOMAIN*bc_dm_ext*mesh_velocity[I];
2096 }
2097 for (int i=0;i<nDOF_test_element;i++)
2098 for (int j=0;j<nDOF_trial_element;j++)
2099 {
2100 int ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j;
2101 double advJacobian_ext = 0.0, diffJacobian_ext = 0.0;
2102 exteriorNumericalAdvectiveFluxDerivative(isDOFBoundary_u.data()[ebNE_kb],isFluxBoundary_u.data()[ebNE_kb],forceStrongConditions,normal,df_ext,advJacobian_ext);
2103 exteriorNumericalDiffusiveFluxDerivative(isDOFBoundary_u.data()[ebNE_kb],isDiffusiveFluxBoundary_u.data()[ebNE_kb],a_rowptr.data(),a_colind.data(),normal,a_ext,da_ext,grad_u_ext,&u_grad_trial_trace[j*nSpace],u_trial_trace_ref.data()[ebN_local_kb_j],ebqe_penalty_ext.data()[ebNE_kb],diffJacobian_ext);
2104 difffluxjacobian_ext = advJacobian_ext*u_trial_trace_ref.data()[ebN_local_kb_j] + diffJacobian_ext;
2105 fluxJacobian_u_u[i][j] += difffluxjacobian_ext*u_test_dS[i];
2106 }//j
2107 }//kb
2108 for (int i=0;i<nDOF_test_element;i++)
2109 {
2110 int eN_i = eN*nDOF_test_element+i;
2111 for (int j=0;j<nDOF_trial_element;j++)
2112 {
2113 int ebN_i_j = ebN*4*nDOF_test_X_trial_element + i*nDOF_trial_element + j;
2114 globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j]] += fluxJacobian_u_u[i][j];
2115 }//j
2116 }//i
2117 }//ebNE
2118 return; // skip the element/boundary Jacobian loops below
2119 }
2120
2121 //
2122 //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
2123 //
2124 for(int eN=0;eN<nElements_global;eN++)
2125 {
2126 double elementJacobian_u_u[nDOF_test_element][nDOF_trial_element];
2127 for (int i=0;i<nDOF_test_element;i++)
2128 for (int j=0;j<nDOF_trial_element;j++)
2129 {
2130 elementJacobian_u_u[i][j]=0.0;
2131 }
2132 for (int k=0;k<nQuadraturePoints_element;k++)
2133 {
2134 int eN_k = eN*nQuadraturePoints_element+k, //index to a scalar at a quadrature point
2135 eN_k_nSpace = eN_k*nSpace,
2136 eN_nDOF_trial_element = eN*nDOF_trial_element; //index to a vector at a quadrature point
2137
2138 //declare local storage
2139 double u=0.0,
2140 grad_u[nSpace],
2141 m=0.0,dm=0.0,
2142 f[nSpace],df[nSpace],
2143 a[nnz],da[nnz],
2144
2145 m_t=0.0,dm_t=0.0,
2146 dpdeResidual_u_u[nDOF_trial_element],
2147 Lstar_u[nDOF_test_element],
2148 dsubgridError_u_u[nDOF_trial_element],
2149 tau=0.0,tau0=0.0,tau1=0.0,
2150 jac[nSpace*nSpace],
2151 jacDet,
2152 jacInv[nSpace*nSpace],
2153 u_grad_trial[nDOF_trial_element*nSpace],
2154 dV,
2155 u_test_dV[nDOF_test_element],
2156 u_grad_test_dV[nDOF_test_element*nSpace],
2157 x,y,z,xt,yt,zt,
2158 G[nSpace*nSpace],G_dd_G,tr_G;
2159 //
2160 //calculate solution and gradients at quadrature points
2161 //
2162 //get jacobian, etc for mapping reference element
2163 ck.calculateMapping_element(eN,
2164 k,
2165 mesh_dof.data(),
2166 mesh_l2g.data(),
2167 mesh_trial_ref.data(),
2168 mesh_grad_trial_ref.data(),
2169 jac,
2170 jacDet,
2171 jacInv,
2172 x,y,z);
2173 ck.calculateMappingVelocity_element(eN,
2174 k,
2175 mesh_velocity_dof.data(),
2176 mesh_l2g.data(),
2177 mesh_trial_ref.data(),
2178 xt,yt,zt);
2179 //get the physical integration weight
2180 dV = fabs(jacDet)*dV_ref.data()[k];
2181 ck.calculateG(jacInv,G,G_dd_G,tr_G);
2182 //get the trial function gradients
2183 ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
2184 //get the solution
2185 ck.valFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],u);
2186 //get the solution gradients
2187 ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial,grad_u);
2188 //precalculate test function products with integration weights
2189 for (int j=0;j<nDOF_trial_element;j++)
2190 {
2191 u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
2192 for (int I=0;I<nSpace;I++)
2193 {
2194 u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
2195 }
2196 }
2197 //
2198 //calculate pde coefficients and derivatives at quadrature points
2199 //
2200
2201 double rho_out=0.0;
2202 evaluateCoefficients(a_rowptr.data(),
2203 a_colind.data(),
2204 &velocity.data()[eN_k_nSpace],
2205 alpha_L,
2206 alpha_T,
2207 Dm,
2208 q_porosity.data()[eN*nQuadraturePoints_element+k],
2209 rho_f,
2210 rho_s,
2211 u,
2212 rho_out,
2213 m,
2214 dm,
2215 f,
2216 df,
2217 a,
2218 da);
2219 //
2220 //moving mesh
2221 //
2222 double mesh_velocity[3];
2223 mesh_velocity[0] = xt;
2224 mesh_velocity[1] = yt;
2225 mesh_velocity[2] = zt;
2226
2227 for(int I=0;I<nSpace;I++)
2228 {
2229 f[I] -= MOVING_DOMAIN*m*mesh_velocity[I];
2230 df[I] -= MOVING_DOMAIN*dm*mesh_velocity[I];
2231 }
2232 //
2233 //calculate time derivatives
2234 //
2235 ck.bdf(alphaBDF,
2236 q_m_betaBDF.data()[eN_k],//since m_t isn't used, we don't have to correct mass
2237 m,
2238 dm,
2239 m_t,
2240 dm_t);
2241 if (STABILIZATION_TYPE == STABILIZATION::VMS)
2242 {
2243 //
2244 //calculate subgrid error contribution to the Jacobian (strong residual, adjoint, jacobian of strong residual)
2245 //
2246 //calculate the adjoint times the test functions
2247 for (int i=0;i<nDOF_test_element;i++)
2248 {
2249 int i_nSpace = i*nSpace;
2250 Lstar_u[i]=ck.Advection_adjoint(df,&u_grad_test_dV[i_nSpace]);
2251 }
2252 //calculate the Jacobian of strong residual
2253 for (int j=0;j<nDOF_trial_element;j++)
2254 {
2255 int j_nSpace = j*nSpace;
2256 dpdeResidual_u_u[j]= ck.MassJacobian_strong(dm_t,u_trial_ref.data()[k*nDOF_trial_element+j]) +
2257 ck.AdvectionJacobian_strong(df,&u_grad_trial[j_nSpace]);
2258 }
2259 //tau and tau*Res
2260 calculateSubgridError_tau(elementDiameter.data()[eN],
2261 dm_t,
2262 df,
2263 cfl.data()[eN_k],
2264 tau0);
2265
2266 calculateSubgridError_tau(Ct_sge,
2267 G,
2268 dm_t,
2269 df,
2270 tau1,
2271 cfl.data()[eN_k]);
2272 tau = useMetrics*tau1+(1.0-useMetrics)*tau0;
2273
2274 for(int j=0;j<nDOF_trial_element;j++)
2275 dsubgridError_u_u[j] = -tau*dpdeResidual_u_u[j];
2276 }
2277 for(int i=0;i<nDOF_test_element;i++)
2278 {
2279 for(int j=0;j<nDOF_trial_element;j++)
2280 {
2281 int j_nSpace = j*nSpace;
2282 int i_nSpace = i*nSpace;
2283 if (STABILIZATION_TYPE==STABILIZATION::Galerkin)
2284 {
2285 elementJacobian_u_u[i][j] +=
2286 ck.MassJacobian_weak(dm_t,
2287 u_trial_ref.data()[k*nDOF_trial_element+j],
2288 u_test_dV[i]) +
2289 ck.AdvectionJacobian_weak(df,
2290 u_trial_ref.data()[k*nDOF_trial_element+j],
2291 &u_grad_test_dV[i_nSpace]) +
2292 ck.DiffusionJacobian_weak(a_rowptr.data(),a_colind.data(),a,da,
2293 grad_u,&u_grad_test_dV[i_nSpace],1.0,
2294 u_trial_ref.data()[k*nDOF_trial_element+j],&u_grad_trial[j_nSpace])
2295 +
2296 ck.NumericalDiffusionJacobian(physicalDiffusion,
2297 &u_grad_trial[j_nSpace],
2298 &u_grad_test_dV[i_nSpace]); //implicit
2299 }
2300 else if (STABILIZATION_TYPE==STABILIZATION::VMS)
2301 {
2302 elementJacobian_u_u[i][j] +=
2303 ck.MassJacobian_weak(dm_t,
2304 u_trial_ref.data()[k*nDOF_trial_element+j],
2305 u_test_dV[i]) +
2306 ck.AdvectionJacobian_weak(df,
2307 u_trial_ref.data()[k*nDOF_trial_element+j],
2308 &u_grad_test_dV[i_nSpace]) +
2309 ck.DiffusionJacobian_weak(a_rowptr.data(),a_colind.data(),a,da,
2310 grad_u,&u_grad_test_dV[i_nSpace],1.0,
2311 u_trial_ref.data()[k*nDOF_trial_element+j],&u_grad_trial[j_nSpace])+
2312
2313 ck.SubgridErrorJacobian(dsubgridError_u_u[j],Lstar_u[i]) +
2314 ck.NumericalDiffusionJacobian(q_numDiff_u_last.data()[eN_k] + physicalDiffusion,
2315 &u_grad_trial[j_nSpace],
2316 &u_grad_test_dV[i_nSpace]); //implicit
2317 }
2318 else if (STABILIZATION_TYPE==STABILIZATION::TaylorGalerkinEV or
2319 STABILIZATION_TYPE==STABILIZATION::EntropyViscosity or
2320 STABILIZATION_TYPE==STABILIZATION::SmoothnessIndicator or
2321 STABILIZATION_TYPE==STABILIZATION::Kuzmin)
2322 {
2323 elementJacobian_u_u[i][j] +=
2324 ck.MassJacobian_weak(1.0,
2325 u_trial_ref.data()[k*nDOF_trial_element+j],
2326 u_test_dV[i]);
2327 }
2328 }//j
2329 }//i
2330 }//k
2331 //
2332 //load into element Jacobian into global Jacobian
2333 //
2334 for (int i=0;i<nDOF_test_element;i++)
2335 {
2336 int eN_i = eN*nDOF_test_element+i;
2337 for (int j=0;j<nDOF_trial_element;j++)
2338 {
2339 int eN_i_j = eN_i*nDOF_trial_element+j;
2340 globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_u_u.data()[eN_i_j]] += elementJacobian_u_u[i][j];
2341 }//j
2342 }//i
2343 }//elements
2344 //
2345 //loop over exterior element boundaries to compute the surface integrals and load them into the global Jacobian
2346 //
2347 if (STABILIZATION_TYPE==STABILIZATION::VMS or STABILIZATION_TYPE==STABILIZATION::Galerkin)
2348 {
2349 for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
2350 {
2351 int ebN = exteriorElementBoundariesArray.data()[ebNE];
2352 const int eN_out = elementBoundaryElementsArray.data()[ebN*2+1];
2353 const int ebFlag = elementBoundaryMaterialTypes.data()[ebN];
2354 if (ebFlag <= 0 || isExteriorBoundaryPhysical.data()[ebNE] == 0 || eN_out >= 0)
2355 continue;
2356 int eN = elementBoundaryElementsArray.data()[ebN*2+0],
2357 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
2358 eN_nDOF_trial_element = eN*nDOF_trial_element;
2359 double fluxJacobian_u_u[nDOF_test_element][nDOF_trial_element];
2360 for (int i=0;i<nDOF_test_element;i++)
2361 for (int j=0;j<nDOF_trial_element;j++)
2362 {
2363 fluxJacobian_u_u[i][j]=0.0;
2364 }
2365 for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
2366 {
2367 int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
2368 ebNE_kb_nSpace = ebNE_kb*nSpace,
2369 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
2370 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
2371 double u_ext=0.0,
2372 grad_u_ext[nSpace],
2373 m_ext=0.0,
2374 dm_ext=0.0,
2375 f_ext[nSpace],
2376 df_ext[nSpace],
2377
2378 a_ext[nnz],
2379 da_ext[nnz],
2380 bc_a_ext[nnz],
2381 bc_da_ext[nnz],
2382
2383 dflux_u_u_ext=0.0,
2384 difffluxjacobian_ext=0.0,
2385 bc_u_ext=0.0,
2386 //bc_grad_u_ext[nSpace],
2387 bc_m_ext=0.0,
2388 bc_dm_ext=0.0,
2389 bc_f_ext[nSpace],
2390 bc_df_ext[nSpace],
2392 diffusiveFluxJacobian_u_u[nDOF_trial_element],
2394 jac_ext[nSpace*nSpace],
2395 jacDet_ext,
2396 jacInv_ext[nSpace*nSpace],
2397 boundaryJac[nSpace*(nSpace-1)],
2398 metricTensor[(nSpace-1)*(nSpace-1)],
2399 metricTensorDetSqrt,
2400 dS,
2401 u_test_dS[nDOF_test_element],
2402 u_grad_trial_trace[nDOF_trial_element*nSpace],
2403 u_grad_test_dS[nDOF_trial_element*nSpace],
2404 normal[nSpace],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
2405 //
2406 G[nSpace*nSpace],G_dd_G,tr_G;
2407
2408 //
2409 //calculate the solution and gradients at quadrature points
2410 //
2411 ck.calculateMapping_elementBoundary(eN,
2412 ebN_local,
2413 kb,
2414 ebN_local_kb,
2415 mesh_dof.data(),
2416 mesh_l2g.data(),
2417 mesh_trial_trace_ref.data(),
2418 mesh_grad_trial_trace_ref.data(),
2419 boundaryJac_ref.data(),
2420 jac_ext,
2421 jacDet_ext,
2422 jacInv_ext,
2423 boundaryJac,
2424 metricTensor,
2425 metricTensorDetSqrt,
2426 normal_ref.data(),
2427 normal,
2428 x_ext,y_ext,z_ext);
2429 ck.calculateMappingVelocity_elementBoundary(eN,
2430 ebN_local,
2431 kb,
2432 ebN_local_kb,
2433 mesh_velocity_dof.data(),
2434 mesh_l2g.data(),
2435 mesh_trial_trace_ref.data(),
2436 xt_ext,yt_ext,zt_ext,
2437 normal,
2438 boundaryJac,
2439 metricTensor,
2440 integralScaling);
2441 dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref.data()[kb];
2442 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
2443 //compute shape and solution information
2444 //shape
2445 ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,u_grad_trial_trace);
2446 //solution and gradients
2447 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);
2448 ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial_trace,grad_u_ext);
2449 //precalculate test function products with integration weights
2450 for (int j=0;j<nDOF_trial_element;j++)
2451 {
2452 u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
2453 for (int I=0;I<nSpace;I++)
2454 {
2455 u_grad_test_dS[j*nSpace+I]= u_grad_trial_trace[j*nSpace+I]*dS;
2456 }
2457 }
2458 //
2459 //load the boundary values
2460 //
2461 bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+(1-isDOFBoundary_u.data()[ebNE_kb])*u_ext;
2462
2463 //
2464 //
2465 //calculate the internal and external trace of the pde coefficients
2466 //
2467
2468 double rho_out_ext=0.0,rho_out_bc=0.0;
2469 evaluateCoefficients(a_rowptr.data(),
2470 a_colind.data(),
2471 &ebqe_velocity_ext.data()[ebNE_kb_nSpace],
2472 alpha_L,
2473 alpha_T,
2474 Dm,
2475 ebqe_porosity.data()[ebNE_kb],
2476 rho_f,
2477 rho_s,
2478 u_ext,
2479 rho_out_ext,
2480 m_ext,
2481 dm_ext,
2482 f_ext,
2483 df_ext,
2484 a_ext,
2485 da_ext
2486 );
2487
2488 evaluateCoefficients(a_rowptr.data(),
2489 a_colind.data(),
2490 &ebqe_velocity_ext.data()[ebNE_kb_nSpace],
2491 alpha_L,
2492 alpha_T,
2493 Dm,
2494 ebqe_porosity.data()[ebNE_kb],
2495 rho_f,
2496 rho_s,
2497 bc_u_ext,
2498 rho_out_bc,
2499 bc_m_ext,
2500 bc_dm_ext,
2501 bc_f_ext,
2502 bc_df_ext,
2503 bc_a_ext,
2504 bc_da_ext);
2505 //
2506 //moving domain
2507 //
2508 double mesh_velocity[3];
2509 mesh_velocity[0] = xt_ext;
2510 mesh_velocity[1] = yt_ext;
2511 mesh_velocity[2] = zt_ext;
2512 for (int I=0;I<nSpace;I++)
2513 {
2514 f_ext[I] -= MOVING_DOMAIN*m_ext*mesh_velocity[I];
2515 df_ext[I] -= MOVING_DOMAIN*dm_ext*mesh_velocity[I];
2516 bc_f_ext[I] -= MOVING_DOMAIN*bc_m_ext*mesh_velocity[I];
2517 bc_df_ext[I] -= MOVING_DOMAIN*bc_dm_ext*mesh_velocity[I];
2518 }
2519 //
2520 //calculate the numerical fluxes
2521 //
2522 //
2523 //calculate the flux jacobian
2524 //
2525 for (int i=0;i<nDOF_test_element;i++)
2526 for (int j=0;j<nDOF_trial_element;j++)
2527 {
2528 int ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j;
2529
2530 double advJacobian_ext = 0.0, diffJacobian_ext = 0.0;
2531 exteriorNumericalAdvectiveFluxDerivative(isDOFBoundary_u.data()[ebNE_kb],
2532 isFluxBoundary_u.data()[ebNE_kb],
2533 forceStrongConditions,
2534 normal,
2535 df_ext,
2536 advJacobian_ext);
2537 exteriorNumericalDiffusiveFluxDerivative(isDOFBoundary_u.data()[ebNE_kb],
2538 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
2539 a_rowptr.data(),
2540 a_colind.data(),
2541 normal,
2542 a_ext,
2543 da_ext,
2544 grad_u_ext,
2545 &u_grad_trial_trace[j*nSpace],
2546 u_trial_trace_ref.data()[ebN_local_kb_j],
2547 ebqe_penalty_ext.data()[ebNE_kb],
2548 diffJacobian_ext);
2549 difffluxjacobian_ext = advJacobian_ext*u_trial_trace_ref.data()[ebN_local_kb_j]
2550 + diffJacobian_ext;
2551 fluxJacobian_u_u[i][j] += difffluxjacobian_ext*u_test_dS[i];
2552 }//j
2553
2554 //
2555 //update the global Jacobian from the flux Jacobian
2556 //
2557 for (int i=0;i<nDOF_test_element;i++)
2558 {
2559 int eN_i = eN*nDOF_test_element+i;
2560 for (int j=0;j<nDOF_trial_element;j++)
2561 {
2562 int ebN_i_j = ebN*4*nDOF_test_X_trial_element + i*nDOF_trial_element + j;
2563 globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j]] += fluxJacobian_u_u[i][j];
2564//
2565 }//j
2566 }//i
2567 }//kb
2568 }//ebNE
2569 }//VMS and Galerkin
2570 }//computeJacobian
2571
2572 void FCTStep(arguments_dict& args)
2573 {
2574 double dt = args.scalar<double>("dt");
2575 int NNZ = args.scalar<int>("NNZ");
2576 int numDOFs = args.scalar<int>("numDOFs");
2577 xt::pyarray<double>& lumped_mass_matrix = args.array<double>("lumped_mass_matrix");
2578 xt::pyarray<double>& soln = args.array<double>("soln");
2579 xt::pyarray<double>& solH = args.array<double>("solH");
2580 xt::pyarray<double>& uLow = args.array<double>("uLow");
2581 xt::pyarray<double>& dLow = args.array<double>("dLow");
2582 xt::pyarray<double>& limited_solution = args.array<double>("limited_solution");
2583 // bc_mask: 0.0 at Dirichlet DOFs, 1.0 elsewhere. Mirrors mphase_co2's
2584 // pattern. Multiplied into the antidiffusive flux correction below so
2585 // Dirichlet DOFs stay at their low-order value (which already carries the
2586 // Nitsche BC contribution from the boundary residual) instead of being
2587 // antidiffused past the BC value.
2588 xt::pyarray<double>& bc_mask = args.array<double>("bc_mask");
2589 xt::pyarray<int>& csrRowIndeces_DofLoops = args.array<int>("csrRowIndeces_DofLoops");
2590 xt::pyarray<int>& csrColumnOffsets_DofLoops = args.array<int>("csrColumnOffsets_DofLoops");
2591 xt::pyarray<double>& MassMatrix = args.array<double>("MassMatrix");
2592 xt::pyarray<double>& dt_times_dH_minus_dL = args.array<double>("dt_times_dH_minus_dL");
2593 xt::pyarray<double>& min_u_bc = args.array<double>("min_u_bc");
2594 xt::pyarray<double>& max_u_bc = args.array<double>("max_u_bc");
2595 int LUMPED_MASS_MATRIX = args.scalar<int>("LUMPED_MASS_MATRIX");
2596 // projection arrays replacing nodal_porosity
2597 xt::pyarray<double>& q_porosity_old_fct = args.array<double>("q_porosity_old_fct");
2598 xt::pyarray<double>& q_rho_fct = args.array<double>("q_rho_fct");
2599 xt::pyarray<double>& q_dV_fct = args.array<double>("q_dV_fct");
2600 xt::pyarray<int>& u_l2g_fct = args.array<int>("u_l2g_fct");
2601 xt::pyarray<double>& u_test_ref_fct = args.array<double>("u_test_ref_fct");
2602 xt::pyarray<double>& theta_dof_out = args.array<double>("theta_dof_out");
2603 int nElements_global_fct = args.scalar<int>("nElements_global_fct");
2604 int nQuadraturePoints_element_fct = args.scalar<int>("nQuadraturePoints_element_fct");
2605 int nDOF_trial_element_fct = args.scalar<int>("nDOF_trial_element_fct");
2606 const double rho_f = args.scalar<double>("rho_f");
2607 const double rho_s = args.scalar<double>("rho_s");
2608// STABILIZATION STABILIZATION_TYPE{args.scalar<int>("STABILIZATION_TYPE")};
2609 STABILIZATION STABILIZATION_TYPE = static_cast<STABILIZATION>(args.scalar<int>("STABILIZATION_TYPE"));
2610// ENTROPY ENTROPY_TYPE = static_cast<ENTROPY>(args.scalar<int>("ENTROPY_TYPE"));
2611 // --- L2 projection: quadrature → DOF ---
2612 std::vector<double> theta_dof(numDOFs, 0.0);
2613 std::vector<double> rho_dof(numDOFs, 0.0);
2614 std::vector<double> ML_proj(numDOFs, 0.0);
2615 for (int eN = 0; eN < nElements_global_fct; eN++) {
2616 for (int k = 0; k < nQuadraturePoints_element_fct; k++) {
2617 const int eN_k = eN * nQuadraturePoints_element_fct + k;
2618 const double dV_k = q_dV_fct.data()[eN_k];
2619 const double theta_k = q_porosity_old_fct.data()[eN_k];
2620 const double rho_k = q_rho_fct.data()[eN_k];
2621 for (int i = 0; i < nDOF_trial_element_fct; i++) {
2622 const int gi = u_l2g_fct.data()[eN * nDOF_trial_element_fct + i];
2623 const double w = u_test_ref_fct.data()[k * nDOF_trial_element_fct + i] * dV_k;
2624 theta_dof[gi] += theta_k * w;
2625 rho_dof[gi] += rho_k * w;
2626 ML_proj[gi] += w;
2627 }
2628 }
2629 }
2630 for (int i = 0; i < numDOFs; i++) {
2631 if (ML_proj[i] > 1.0e-14) { theta_dof[i] /= ML_proj[i]; rho_dof[i] /= ML_proj[i]; }
2632 else { theta_dof[i] = 1.0; rho_dof[i] = rho_f; }
2633 }
2634 Rpos.resize(numDOFs,0.0);
2635 Rneg.resize(numDOFs,0.0);
2636 FluxCorrectionMatrix.resize(NNZ,0.0);
2637 int ij=0;
2638 //loop over nodes (i)
2639 for (int i=0; i<numDOFs; i++)
2640 {
2641 //read some vectors
2642 double solHi = solH.data()[i];
2643 double solni = soln.data()[i];
2644 const double lumped_volume = lumped_mass_matrix.data()[i];
2645 double uLowi = uLow.data()[i];
2646 const double theta_i = theta_dof[i];
2647 double mLowi = theta_i*rho_f*(1.0 + ((rho_s-rho_f)/rho_f)*uLowi)*uLowi;
2648 double solHmi = theta_i*rho_f*(1.0 + ((rho_s-rho_f)/rho_f)*solHi)*solHi;
2649 double solnmi = theta_i*rho_f*(1.0 + ((rho_s-rho_f)/rho_f)*solni)*solni;
2650 double uDotLowi = (mLowi - solnmi)/dt;
2651 double mini=min_u_bc.data()[i], maxi=max_u_bc.data()[i]; // init min/max with value at BCs (NOTE: if no boundary then min=1E10, max=-1E10)
2652 double Pposi=0, Pnegi=0;
2653 // Loop over neighbors (j)
2654 for (int offset=csrRowIndeces_DofLoops.data()[i]; offset<csrRowIndeces_DofLoops.data()[i+1]; offset++)
2655 {
2656 assert(offset == ij); // (CSR matrix consistency)
2657 int j = csrColumnOffsets_DofLoops.data()[offset];
2658 double solnj = soln.data()[j];
2659 double uLowj = uLow.data()[j];
2661 // COMPUTE THE BOUNDS //
2663
2664 // Explicit paths bound against soln (as VOF.h/VOF3P.h do);
2665 // ImplicitEV, whose limiter works on mass, bounds against uLow.
2666 const double bound_j =
2667 (STABILIZATION_TYPE == STABILIZATION::ImplicitEV) ? uLowj : solnj;
2668 mini = fmin(mini,bound_j);
2669 maxi = fmax(maxi,bound_j);
2670 const double theta_j = theta_dof[j];
2671 double mLowj = theta_j*rho_f*(1.0 + ((rho_s-rho_f)/rho_f)*uLowj)*uLowj;
2672 double solHmj = theta_j*rho_f*(1.0 + ((rho_s-rho_f)/rho_f)*solH.data()[j])*solH.data()[j];
2673 double solnmj = theta_j*rho_f*(1.0 + ((rho_s-rho_f)/rho_f)*solnj)*solnj;
2674 double uDotLowj = (mLowj - solnmj)/dt;
2675 // i-th row of flux correction matrix
2676 if (STABILIZATION_TYPE == STABILIZATION::Kuzmin)
2677 {
2678 // Explicit Kuzmin antidiffusive flux: consistent-mass correction +
2679 // full removal of the symmetric graph viscosity dLow*(mLow_i-mLow_j).
2680 FluxCorrectionMatrix[ij] = dt*(MassMatrix.data()[ij]*(uDotLowi-uDotLowj)
2681 + dLow.data()[ij]*(mLowi-mLowj));
2682 }
2683 else if (STABILIZATION_TYPE == STABILIZATION::ImplicitEV)
2684 {
2685 // Standard Kuzmin implicit FEM-FCT antidiffusive flux:
2686 // F_ij = M~_ij[(m_j^H - m_j^n) - (m_i^H - m_i^n)]
2687 // + dt (d^H - d^L)_ij (m_j^n - m_i^n)
2688 // For ImplicitEV the Newton solve gives the low-order implicit
2689 // solution, so the high-order state is taken as m^H := m^L (uLow):
2690 // first term -M~_ij[(m_j^L-m_j^n)-(m_i^L-m_i^n)]
2691 // = dt*MC_ij*(uDotLow_i - uDotLow_j) (consistent mass)
2692 // second term dt_times_dH_minus_dL_ij*(m_j^n - m_i^n)
2693 // with dt_times_dH_minus_dL = dt*(d^H-d^L) = -dt*dLow.
2694 // The dissipation antidiffusion now acts on the OLD-time mass m^n
2695 // (solnm), matching the textbook Kuzmin form and the STAB=2
2696 // (EntropyViscosity) `else` branch below. The previous code used the
2697 // diffused low-order m^L here and threw away most of dLow via the Comp
2698 // throttle, so the front stayed under-sharpened. solnm carries the
2699 // density nonlinearity m^s(c)=theta*rho_f*(c + eps c^2); the inverse
2700 // c = 2r/(1+sqrt(1+4 eps r)) is applied in invert().
2701 FluxCorrectionMatrix[ij] = dt*MassMatrix.data()[ij]*(uDotLowi-uDotLowj)
2702 + dt_times_dH_minus_dL.data()[ij]*(solnmj-solnmi);
2703 }
2704 else
2705 {
2706 double ML_minus_MC =
2707 (LUMPED_MASS_MATRIX == 1 ? 0. : (i==j ? 1. : 0.)*lumped_volume - MassMatrix.data()[ij]);
2708 FluxCorrectionMatrix[ij] = ML_minus_MC * (solHmj-solnmj - (solHmi-solnmi))
2709 + dt_times_dH_minus_dL.data()[ij]*(solnmj-solnmi);
2710 }
2711 Pposi += FluxCorrectionMatrix[ij]*((FluxCorrectionMatrix[ij] > 0) ? 1. : 0.);
2712 Pnegi += FluxCorrectionMatrix[ij]*((FluxCorrectionMatrix[ij] < 0) ? 1. : 0.);
2713 ij+=1;
2714 }//j
2715 const double Qposi = lumped_volume*(theta_i*rho_f*(1.0 + ((rho_s-rho_f)/rho_f)*maxi)*maxi - mLowi);
2716 const double Qnegi = lumped_volume*(theta_i*rho_f*(1.0 + ((rho_s-rho_f)/rho_f)*mini)*mini - mLowi);
2717 Rpos[i] = ((Pposi==0) ? 1. : fmax(0.0, fmin(1.0,Qposi/Pposi)));
2718 Rneg[i] = ((Pnegi==0) ? 1. : fmax(0.0, fmin(1.0,Qnegi/Pnegi)));
2719 }//i
2720 ij=0;
2721 for (int i=0; i<numDOFs; i++)
2722 {
2723 double ith_Limiter_times_FluxCorrectionMatrix = 0.;
2724 double Rposi = Rpos[i], Rnegi = Rneg[i];
2725 const double lumped_volume = lumped_mass_matrix.data()[i];
2726 for (int offset=csrRowIndeces_DofLoops.data()[i]; offset<csrRowIndeces_DofLoops.data()[i+1]; offset++)
2727 {
2728 assert(offset == ij); // (CSR matrix consistency
2729 int j = csrColumnOffsets_DofLoops.data()[offset];
2730 double Lij = 1;
2731 Lij = ((FluxCorrectionMatrix[ij]>0) ? fmin(Rposi,Rneg[j]) : fmin(Rnegi,Rpos[j]));
2732 ith_Limiter_times_FluxCorrectionMatrix += Lij * FluxCorrectionMatrix[ij];
2733 ij+=1;
2734 }
2735 const double uLowi = uLow.data()[i];
2736 const double theta_i = theta_dof[i];
2737 const double mLowi = theta_i*rho_f*(1.0 + ((rho_s-rho_f)/rho_f)*uLowi)*uLowi;
2738 // bc_mask.data()[i] is 0 at Dirichlet DOFs (freeze at mLow) and 1 elsewhere.
2739 limited_solution.data()[i] = mLowi + bc_mask.data()[i] * (1./lumped_volume * ith_Limiter_times_FluxCorrectionMatrix);
2740 theta_dof_out.data()[i] = theta_dof[i];
2741 }
2742 }//FCTStep
2743 };//TADR
2744
2745inline TADR_base* newTADR(int nSpaceIn,
2746 int nQuadraturePoints_elementIn,
2747 int nDOF_mesh_trial_elementIn,
2748 int nDOF_trial_elementIn,
2749 int nDOF_test_elementIn,
2750 int nQuadraturePoints_elementBoundaryIn,
2751 int CompKernelFlag)
2752{
2753 if (nSpaceIn == 1)
2754 return proteus::chooseAndAllocateDiscretization1D<TADR_base,TADR,CompKernel>(nSpaceIn,
2755 nQuadraturePoints_elementIn,
2756 nDOF_mesh_trial_elementIn,
2757 nDOF_trial_elementIn,
2758 nDOF_test_elementIn,
2759 nQuadraturePoints_elementBoundaryIn,
2760 CompKernelFlag);
2761 else if (nSpaceIn == 2)
2762 return proteus::chooseAndAllocateDiscretization2D<TADR_base,TADR,CompKernel>(nSpaceIn,
2763 nQuadraturePoints_elementIn,
2764 nDOF_mesh_trial_elementIn,
2765 nDOF_trial_elementIn,
2766 nDOF_test_elementIn,
2767 nQuadraturePoints_elementBoundaryIn,
2768 CompKernelFlag);
2769 else
2770 return proteus::chooseAndAllocateDiscretization<TADR_base,TADR,CompKernel>(nSpaceIn,
2771 nQuadraturePoints_elementIn,
2772 nDOF_mesh_trial_elementIn,
2773 nDOF_trial_elementIn,
2774 nDOF_test_elementIn,
2775 nQuadraturePoints_elementBoundaryIn,
2776 CompKernelFlag);
2777}
2778}//proteus
2779#endif
n
Int n
Definition Headers.h:28
H
Double H
Definition Headers.h:65
u
Double u
Definition Headers.h:89
num
Int num
Definition Headers.h:32
z
Double * z
Definition Headers.h:49
v
Double v
Definition Headers.h:95
cMax
#define cMax
Definition NCLS3P.h:11
cE
#define cE
Definition NCLS3P.h:10
proteus::TADR_base
Definition TADR.h:61
proteus::TADR_base::theta_dof_proj
std::valarray< double > theta_dof_proj
Definition TADR.h:68
proteus::TADR_base::m_dof
std::valarray< double > m_dof
Definition TADR.h:68
proteus::TADR_base::boundary_integral
std::valarray< double > boundary_integral
Definition TADR.h:67
proteus::TADR_base::eta
std::valarray< double > eta
Definition TADR.h:67
proteus::TADR_base::DiffusionMatrix
std::valarray< double > DiffusionMatrix
Definition TADR.h:66
proteus::TADR_base::Rneg
std::valarray< double > Rneg
Definition TADR.h:64
proteus::TADR_base::maxEntRes
std::valarray< double > maxEntRes
Definition TADR.h:69
proteus::TADR_base::calculateResidual
virtual void calculateResidual(arguments_dict &args)=0
proteus::TADR_base::ML_mass_proj
std::valarray< double > ML_mass_proj
Definition TADR.h:68
proteus::TADR_base::TransposeTransportMatrix
std::valarray< double > TransposeTransportMatrix
Definition TADR.h:66
proteus::TADR_base::FCTStep
virtual void FCTStep(arguments_dict &args)=0
proteus::TADR_base::maxVel
std::valarray< double > maxVel
Definition TADR.h:69
proteus::TADR_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)=0
proteus::TADR_base::invert
virtual void invert(arguments_dict &args)=0
proteus::TADR_base::~TADR_base
virtual ~TADR_base()
Definition TADR.h:70
proteus::TADR_base::Rpos
std::valarray< double > Rpos
Definition TADR.h:64
proteus::TADR_base::TransportMatrix
std::valarray< double > TransportMatrix
Definition TADR.h:66
proteus::TADR_base::FluxCorrectionMatrix
std::valarray< double > FluxCorrectionMatrix
Definition TADR.h:65
proteus::TADR_base::rho_dof_proj
std::valarray< double > rho_dof_proj
Definition TADR.h:68
proteus::TADR_base::global_entropy_residual
std::valarray< double > global_entropy_residual
Definition TADR.h:67
proteus::TADR_base::psi
std::valarray< double > psi
Definition TADR.h:67
proteus::TADR::inversevaluateCoefficients
double inversevaluateCoefficients(const double storage, const double porosity, const double rho_f, const double rho_s)
Definition TADR.h:169
proteus::TADR::calculateResidual
void calculateResidual(arguments_dict &args)
Definition TADR.h:439
proteus::TADR::ExteriorNumericalDiffusiveFluxJacobian
double ExteriorNumericalDiffusiveFluxJacobian(int *rowptr, int *colind, const int &isDOFBoundary, const int &isDiffusiveFluxBoundary, const double n[nSpace], double *a, const double &v, const double grad_v[nSpace], const double &penalty)
Definition TADR.h:233
proteus::TADR::exteriorNumericalAdvectiveFluxDerivative
void exteriorNumericalAdvectiveFluxDerivative(const int &isDOFBoundary_u, const int &isFluxBoundary_u, const int &forceStrongConditions, const double n[nSpace], const double velocity[nSpace], double &dflux)
Definition TADR.h:365
proteus::TADR::calculateJacobian
void calculateJacobian(arguments_dict &args)
Definition TADR.h:1860
proteus::TADR::exteriorNumericalDiffusiveFluxDerivative
void exteriorNumericalDiffusiveFluxDerivative(const int &isDOFBoundary, const int &isDiffusiveFluxBoundary, const int rowptr[nSpace], const int colind[nnz], const double n[nSpace], const double a[nnz], const double da[nnz], const double grad_psi[nSpace], const double grad_v[nSpace], const double &v, const double penalty, double &fluxJacobian)
Definition TADR.h:406
proteus::TADR::nDOF_test_X_trial_element
const int nDOF_test_X_trial_element
Definition TADR.h:87
proteus::TADR::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 TADR.h:283
proteus::TADR::exteriorNumericalAdvectiveFlux
void exteriorNumericalAdvectiveFlux(const int &isDOFBoundary_u, const int &isFluxBoundary_u, const int &forceStrongConditions, const double n[nSpace], const double &bc_flux_u, const double f[nSpace], const double bc_f[nSpace], const double velocity[nSpace], double &flux)
Definition TADR.h:326
proteus::TADR::FCTStep
void FCTStep(arguments_dict &args)
Definition TADR.h:2572
proteus::TADR::exteriorNumericalDiffusiveFlux
void exteriorNumericalDiffusiveFlux(int *rowptr, int *colind, const int &isDOFBoundary, const int &isDiffusiveFluxBoundary, const double n[nSpace], double *bc_a, const double &bc_u, const double &bc_flux, double *a, const double grad_potential[nSpace], const double &u, const double &penalty, double &flux)
Definition TADR.h:187
proteus::TADR::invert
void invert(arguments_dict &args)
Definition TADR.h:1847
proteus::TADR::calculateSubgridError_tau
void calculateSubgridError_tau(const double &elementDiameter, const double &dmt, const double df[nSpace], double &cfl, double &tau)
Definition TADR.h:264
proteus::TADR::calculateNumericalDiffusion
void calculateNumericalDiffusion(const double &shockCapturingDiffusion, const double &elementDiameter, const double &strong_residual, const double grad_u[nSpace], double &numDiff)
Definition TADR.h:306
proteus::TADR::ck
CompKernelType ck
Definition TADR.h:88
proteus::TADR::evaluateCoefficients
void evaluateCoefficients(const int rowptr[nSpace], const int colind[nnz], const double v[nSpace], const double alpha_L, const double alpha_T, const double Dm, const double thetaW, const double rho_f, const double rho_s, const double &u, double &rho_out, double &m, double &dm, double f[nSpace], double df[nSpace], double a[nnz], double da[nnz])
Definition TADR.h:113
proteus::TADR::calculateCFL
void calculateCFL(const double &elementDiameter, const double df[nSpace], double &cfl)
Definition TADR.h:97
df
double df(double C, double b, double a, int q, int r)
Definition analyticalSolutions.c:2209
w
#define w(x)
Definition jf.h:22
POWER_SMOOTHNESS_INDICATOR
#define POWER_SMOOTHNESS_INDICATOR
Definition m_comp_co2.h:22
nnz
#define nnz
Definition m_comp_co2.h:19
proteus
Definition ADR.h:19
proteus::newTADR
TADR_base * newTADR(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition TADR.h:2745
proteus::DEPOWER
double DEPOWER(const double &phi, const double &phiL, const double &phiR)
Definition TADR.h:43
proteus::phi
double phi(const double &g, const double &h, const double &hL, const double &hR, const double &uL, const double &uR)
Definition SW2DCV.h:62
proteus::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::DELOG
double DELOG(const double &phi, const double &phiL, const double &phiR)
Definition TADR.h:52
proteus::ELOG
double ELOG(const double &phi, const double &phiL, const double &phiR)
Definition TADR.h:48
proteus::EPOWER
double EPOWER(const double &phi, const double &phiL, const double &phiR)
Definition TADR.h:39
proteus::chooseAndAllocateDiscretization1D
Model_Base * chooseAndAllocateDiscretization1D(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition ModelFactory.h:502
proteus::chooseAndAllocateDiscretization2D
Model_Base * chooseAndAllocateDiscretization2D(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nDOF_v_trial_elementIn, int nDOF_v_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition MixedModelFactory.h:309
proteus::DISPERSION
DISPERSION
Definition TADR.h:37
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