proteus 1.9.0
C/C++/Fortran libraries
Loading...
Searching...
No Matches
Richards.h
Go to the documentation of this file.
1#ifndef Richards_H
2#define Richards_H
3#include <cmath>
4#include <iostream>
5#include <valarray>
6#include "CompKernel.h"
7#include "ModelFactory.h"
8#include "../pskRelations.h"
9#include "../mprans/ArgumentsDict.h"
10#include "xtensor-python/pyarray.hpp"
11#define nnz nSpace
12
13namespace py = pybind11;
14#define POWER_SMOOTHNESS_INDICATOR 2
15#define IS_BETAij_ONE 0
16#define GLOBAL_FCT 0
17namespace proteus
18{
19enum class STABILIZATION : int {
20 Galerkin = 0,
21 EV_Stab = 1,
22 EntropyViscosity = 2,
23 Implicit_FCT = 3
24};
25
26namespace richards
27{
28//cek todo: revisit entry for mass transport form
29// Power entropy //
30inline double ENTROPY(const double &phi, const double &phiL, const double &phiR)
31{
32 return 1. / 2. * std::pow(fabs(phi), 2.);
33}
34inline double DENTROPY(const double &phi, const double &phiL, const double &phiR)
35{
36 return fabs(phi) * (phi >= 0 ? 1 : -1);
37}
38// Log entropy // for level set from 0 to 1
39inline double ENTROPY_LOG(const double &phi, const double &phiL, const double &phiR)
40{
41 return std::log(fabs((phi - phiL) * (phiR - phi)) + 1E-14);
42}
43inline double DENTROPY_LOG(const double &phi, const double &phiL, const double &phiR)
44{
45 return (phiL + phiR - 2 * phi) * ((phi - phiL) * (phiR - phi) >= 0 ? 1 : -1) / (fabs((phi - phiL) * (phiR - phi)) + 1E-14);
46}
47} // namespace richards
48} // namespace proteus
49namespace proteus
50{
51namespace richards
52{
53class Richards_base {
54 //The base class defining the interface
55public:
56 virtual ~Richards_base() { double anb_seepage_flux = 1e-16; }
57 virtual void calculateResidual(arguments_dict &args) = 0;
58 virtual void calculateJacobian(arguments_dict &args) = 0;
59 virtual void invert(arguments_dict &args) = 0;
60 virtual void FCTStep(arguments_dict &args) = 0;
61 virtual void kth_FCT_step(arguments_dict &args) = 0;
62 virtual void calculateResidual_entropy_viscosity(arguments_dict &args) = 0;
63 virtual void calculateMassMatrix(arguments_dict &args) = 0;
64};
65
66template <class CompKernelType, int nSpace, int nQuadraturePoints_element, int nDOF_mesh_trial_element, int nDOF_trial_element, int nDOF_test_element, int nQuadraturePoints_elementBoundary>
67class Richards : public Richards_base {
68public:
69 const int nDOF_test_X_trial_element;
70 CompKernelType ck;
71 // Per-DOF density projected from q_rho in calculateResidual_entropy_viscosity.
72 // Reused by invert() so the m -> u inversion uses the same variable density
73 // that built the forward mass.
74 std::vector<double> rho_dof_member;
75 // Pore size distribution / relative permeability model selected from Python
76 // (Coefficients.PSK_type). 0 = van Genuchten-Mualem (default), 1 = Brooks-
77 // Corey-Burdine, 2 = Brooks-Corey-Mualem, 3 = Gardner. Refreshed from args at
78 // the top of every kernel that evaluates the closure.
79 int PSK_TYPE_member = 0;
80 Richards() : nDOF_test_X_trial_element(nDOF_test_element * nDOF_trial_element), ck() { }
81 inline void evaluateCoefficients(const int rowptr[nSpace], const int colind[nnz], const double rho0, const double rho_transport, const double beta, const double gravity[nSpace], const double alpha, const double n_vg, const double thetaR, const double thetaSR, const double KWs[nnz], const double &u, double &m, double &dm, double f[nSpace], double df[nSpace], double a[nnz], double da[nnz], double as[nnz], double &kr, double &dkr, double &thetaW_out)
82 {
83 const double psiC = -u;
84 double thetaW, DthetaW_DpsiC, KWr, DKWr_DpsiC;
85 if (PSK_TYPE_member == 1 || PSK_TYPE_member == 2) {
86 // Same Brooks-Corey retention curve either way; PSK_TYPE only picks which
87 // k_rw closure supplies the exponent (see pskRelations.h).
88 proteus::richards::psk::bc_wetting(
89 psiC, alpha, n_vg, thetaR, thetaSR,
90 thetaW, DthetaW_DpsiC, KWr, DKWr_DpsiC,
91 PSK_TYPE_member == 2 ? proteus::richards::psk::bc_kr::mualem
92 : proteus::richards::psk::bc_kr::burdine);
93 } else if (PSK_TYPE_member == 3) {
94 proteus::richards::psk::gardner_wetting(
95 psiC, alpha, n_vg, thetaR, thetaSR,
96 thetaW, DthetaW_DpsiC, KWr, DKWr_DpsiC);
97 } else {
98 proteus::richards::psk::vgm_wetting(
99 psiC, alpha, n_vg, thetaR, thetaSR,
100 thetaW, DthetaW_DpsiC, KWr, DKWr_DpsiC);
101 }
102 thetaW_out = thetaW;
103 // Density uses transported salinity scaled by the compressibility factor.
104 const double rhom = rho_transport * exp(beta * u);
105 const double drhom = beta * rhom;
106 m = rhom * thetaW;
107 dm = -rhom * DthetaW_DpsiC + drhom * thetaW;
108 const double rho_ratio = rhom / rho0;
109 for (int I = 0; I < nSpace; I++) {
110 f[I] = 0.0;
111 df[I] = 0.0;
112 for (int ii = rowptr[I]; ii < rowptr[I + 1]; ii++) {
113 a[ii] = rhom * KWr * KWs[ii];
114 da[ii] = (drhom * KWr - rhom * DKWr_DpsiC) * KWs[ii];
115 f[I] += rhom * rho_ratio * KWr * KWs[ii] * gravity[colind[ii]];
116 df[I] += (drhom * rho_ratio * KWr +
117 rhom * (drhom / rho0) * KWr -
118 rhom * rho_ratio * DKWr_DpsiC) * KWs[ii] * gravity[colind[ii]];
119 as[ii] = rhom * KWs[ii];
120 kr = KWr;
121 dkr = -DKWr_DpsiC;
122 }
123 }
124 }
125
126 inline void evaluateInverseCoefficients(const int rowptr[nSpace], const int colind[nnz], const double rho, const double beta, const double gravity[nSpace], const double alpha, const double n_vg, const double thetaR, const double thetaSR, const double KWs[nnz], double &u, const double &m, const double &dm, const double f[nSpace], const double df[nSpace], const double a[nnz], const double da[nnz])
127 {
128 (void)rowptr; (void)colind; (void)beta; (void)gravity; (void)KWs;
129 (void)dm; (void)f; (void)df; (void)a; (void)da;
130
131 // Both Brooks-Corey codes share one retention curve, so both invert with
132 // bc_*; PSK_TYPE 1 vs 2 differs only in k_rw, which plays no part here.
133 if (PSK_TYPE_member == 1 || PSK_TYPE_member == 2) {
134 proteus::richards::psk::bc_invert_analytic(m, rho, alpha, n_vg, thetaR, thetaSR, u);
135 } else if (PSK_TYPE_member == 3) {
136 proteus::richards::psk::gardner_invert_analytic(m, rho, alpha, n_vg, thetaR, thetaSR, u);
137 } else {
138 proteus::richards::psk::vgm_invert_analytic(m, rho, alpha, n_vg, thetaR, thetaSR, u);
139 }
140 }
141
142inline void evaluateInverseCoefficients_Newton(const int rowptr[nSpace],
143 const int colind[nnz],
144 const double rho,
145 const double beta,
146 const double gravity[nSpace],
147 const double alpha,
148 const double n_vg,
149 const double thetaR,
150 const double thetaSR,
151 const double KWs[nnz],
152 double &u, // in/out
153 const double &m,
154 const double &dm,
155 const double f[nSpace],
156 const double df[nSpace],
157 const double a[nnz],
158 const double da[nnz])
159{
160 (void)rowptr; (void)colind; (void)gravity; (void)KWs;
161 (void)dm; (void)f; (void)df; (void)a; (void)da;
162
163 // Newton inversion of the full forward mass m = rho(u)*theta_w(u), so the
164 // exp(beta*u) factor the analytic inverse drops is included. The retention
165 // curve it solves against is the selected PSK model.
166 if (PSK_TYPE_member == 1 || PSK_TYPE_member == 2) {
167 proteus::richards::psk::bc_invert_newton(m, rho, beta, alpha, n_vg, thetaR, thetaSR, u);
168 } else if (PSK_TYPE_member == 3) {
169 proteus::richards::psk::gardner_invert_newton(m, rho, beta, alpha, n_vg, thetaR, thetaSR, u);
170 } else {
171 proteus::richards::psk::vgm_invert_newton(m, rho, beta, alpha, n_vg, thetaR, thetaSR, u);
172 }
173}
174
175 inline void calculateCFL(const double &elementDiameter, const double df[nSpace], double &cfl)
176 {
177 double h, nrm_v;
178 h = elementDiameter;
179 nrm_v = 0.0;
180 for (int I = 0; I < nSpace; I++) nrm_v += df[I] * df[I];
181 nrm_v = sqrt(nrm_v);
182 cfl = nrm_v / h;
183 }
184
185 inline void calculateSubgridError_tau(const double &elementDiameter, const double &dmt, const double dH[nSpace], double &cfl, double &tau)
186 {
187 double h, nrm_v, oneByAbsdt;
188 h = elementDiameter;
189 nrm_v = 0.0;
190 for (int I = 0; I < nSpace; I++) nrm_v += dH[I] * dH[I];
191 nrm_v = sqrt(nrm_v);
192 cfl = nrm_v / h;
193 oneByAbsdt = fabs(dmt);
194 tau = 1.0 / (2.0 * nrm_v / h + oneByAbsdt + 1.0e-8);
195 }
196
197 inline void calculateSubgridError_tau(const double &Ct_sge, const double G[nSpace * nSpace], const double &A0, const double Ai[nSpace], double &tau_v, double &cfl)
198 {
199 double v_d_Gv = 0.0;
200 for (int I = 0; I < nSpace; I++)
201 for (int J = 0; J < nSpace; J++) v_d_Gv += Ai[I] * G[I * nSpace + J] * Ai[J];
202 tau_v = 1.0 / sqrt(Ct_sge * A0 * A0 + v_d_Gv);
203 }
204
205 inline void calculateNumericalDiffusion(const double &shockCapturingDiffusion, const double &elementDiameter, const double &strong_residual, const double grad_u[nSpace], double &numDiff)
206 {
207 double h, num, den, n_grad_u;
208 h = elementDiameter;
209 n_grad_u = 0.0;
210 for (int I = 0; I < nSpace; I++) n_grad_u += grad_u[I] * grad_u[I];
211 num = shockCapturingDiffusion * 0.5 * h * fabs(strong_residual);
212 den = sqrt(n_grad_u) + 1.0e-8;
213 numDiff = num / den;
214 }
215
216 inline void exteriorNumericalFlux(const double &bc_flux, int rowptr[nSpace], int colind[nnz], int isSeepageFace, int &isDOFBoundary, double n[nSpace], double bc_u, double K[nnz], double grad_psi[nSpace], double u, double K_rho_g[nSpace], double penalty, double &flux)
217 {
218 double v_I, bc_u_seepage = 0.0;
219 if (isSeepageFace || isDOFBoundary) {
220 flux = 0.0;
221 for (int I = 0; I < nSpace; I++) {
222 //gravity
223 v_I = K_rho_g[I];
224 //pressure head
225 for (int m = rowptr[I]; m < rowptr[I + 1]; m++) { v_I -= K[m] * grad_psi[colind[m]]; }
226 flux += v_I * n[I];
227 }
228 if (isSeepageFace) bc_u = bc_u_seepage;
229 flux += penalty * (u - bc_u);
230 //flux -= penalty * bc_u;
231 if (isSeepageFace) {
232 if (flux > 0.0) {
233 isDOFBoundary = 1;
234 bc_u = bc_u_seepage;
235 } else {
236 isDOFBoundary = 0;
237 flux = 0.0;
238 }
239 }
240 } else flux = bc_flux;
241 }
242
243 void exteriorNumericalFluxJacobian(const int rowptr[nSpace], const int colind[nnz], const int isDOFBoundary, const double n[nSpace], const double K[nnz], const double dK[nnz], const double grad_psi[nSpace], const double grad_v[nSpace], const double dK_rho_g[nSpace], const double v, const double penalty, double &fluxJacobian)
244 {
245 if (isDOFBoundary) {
246 fluxJacobian = 0.0;
247 for (int I = 0; I < nSpace; I++) {
248 //gravity
249 fluxJacobian += dK_rho_g[I] * v * n[I];
250 //pressure head
251 for (int m = rowptr[I]; m < rowptr[I + 1]; m++) { fluxJacobian -= (K[m] * grad_v[colind[m]] + dK[m] * v * grad_psi[colind[m]]) * n[I]; }
252 }
253 //Dirichlet penalty
254 fluxJacobian += penalty * v;
255 } else fluxJacobian = 0.0;
256 }
257
258inline void exteriorNumericalFlux2(const double &bc_flux, int rowptr[nSpace], int colind[nnz], int isSeepageFace, int &isDOFBoundary, double n[nSpace], double bc_u, double K[nnz], double grad_psi[nSpace], double u, double K_rho_g[nSpace], double penalty, double &flux, double &bflux)
259 {
260 double v_I, bc_u_seepage = 0.0;
261 if (isSeepageFace || isDOFBoundary) {
262 flux = 0.0;
263 bflux = 0.0;
264 for (int I = 0; I < nSpace; I++) {
265 //gravity
266 v_I = K_rho_g[I];
267 //pressure head
268 for (int m = rowptr[I]; m < rowptr[I + 1]; m++) { v_I -= K[m] * grad_psi[colind[m]]; }
269 flux += v_I * n[I];
270 }
271 if (isSeepageFace) bc_u = bc_u_seepage;
272 flux += penalty * (u - bc_u);
273 bflux += penalty * (u - bc_u);
274 if (isSeepageFace) {
275 if (flux > 0.0) {
276 isDOFBoundary = 1;
277 } else {
278 isDOFBoundary = 0;
279 flux = 0.0;
280 bflux = 0.0;
281 }
282 }
283 } else {
284 flux = bc_flux;
285 bflux = bc_flux;
286 }
287 }
288
289 void exteriorNumericalFluxJacobian2(const int rowptr[nSpace], const int colind[nnz], const int isDOFBoundary, const double n[nSpace], const double Ks[nnz], const double K[nnz], const double dK[nnz], const double grad_psi[nSpace], const double grad_v[nSpace], const double dK_rho_g[nSpace], const double v, const double penalty, double &fluxJacobian, double &bfluxJacobian)
290 {
291 if (isDOFBoundary) {
292 fluxJacobian = 0.0;
293 bfluxJacobian = 0.0;
294 for (int I = 0; I < nSpace; I++) {
295 for (int m = rowptr[I]; m < rowptr[I + 1]; m++) {
296 fluxJacobian -= Ks[m] * grad_v[colind[m]] * n[I];
297 }
298 }
299 //Dirichlet penalty
300 bfluxJacobian = penalty * v;
301 } else {
302 fluxJacobian = 0.0;
303 bfluxJacobian = 0.0;
304 }
305 }
306
307 double seepagefluxcalculator(double anb_seepage_flux, int isSeepageFace, double dS, double flux_ext)
308 {
309 if (isSeepageFace) { anb_seepage_flux += flux_ext * dS; }
310 return anb_seepage_flux;
311 }
312
313 void calculateResidual(arguments_dict &args)
314 {
315 xt::pyarray<double> &mesh_trial_ref = args.array<double>("mesh_trial_ref");
316 xt::pyarray<double> &mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
317 xt::pyarray<double> &mesh_dof = args.array<double>("mesh_dof");
318 xt::pyarray<double> &mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
319 double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
320 xt::pyarray<int> &mesh_l2g = args.array<int>("mesh_l2g");
321 xt::pyarray<double> &dV_ref = args.array<double>("dV_ref");
322 xt::pyarray<double> &u_trial_ref = args.array<double>("u_trial_ref");
323 xt::pyarray<double> &u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
324 xt::pyarray<double> &u_test_ref = args.array<double>("u_test_ref");
325 xt::pyarray<double> &u_grad_test_ref = args.array<double>("u_grad_test_ref");
326 xt::pyarray<double> &mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
327 xt::pyarray<double> &mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
328 xt::pyarray<double> &dS_ref = args.array<double>("dS_ref");
329 xt::pyarray<double> &u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
330 xt::pyarray<double> &u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
331 xt::pyarray<double> &u_test_trace_ref = args.array<double>("u_test_trace_ref");
332 xt::pyarray<double> &u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
333 xt::pyarray<double> &normal_ref = args.array<double>("normal_ref");
334 xt::pyarray<double> &boundaryJac_ref = args.array<double>("boundaryJac_ref");
335 int nElements_global = args.scalar<int>("nElements_global");
336 xt::pyarray<double> &ebqe_penalty_ext = args.array<double>("ebqe_penalty_ext");
337 xt::pyarray<int> &elementMaterialTypes = args.array<int>("elementMaterialTypes");
338 xt::pyarray<int> &isSeepageFace = args.array<int>("isSeepageFace");
339 xt::pyarray<int> &a_rowptr = args.array<int>("a_rowptr");
340 xt::pyarray<int> &a_colind = args.array<int>("a_colind");
341 double rho = args.scalar<double>("rho");
342 double beta = args.scalar<double>("beta");
343
345 xt::pyarray<double> &q_rho = args.array<double>("q_rho");
346 xt::pyarray<double> &ebqe_rho = args.array<double>("ebqe_rho");
347
348 xt::pyarray<double> &gravity = args.array<double>("gravity");
349 xt::pyarray<double> &alpha = args.array<double>("alpha");
350 xt::pyarray<double> &n = args.array<double>("n");
351 xt::pyarray<double> &thetaR = args.array<double>("thetaR");
352 xt::pyarray<double> &thetaSR = args.array<double>("thetaSR");
353 PSK_TYPE_member = args.scalar<int>("PSK_TYPE");
354 xt::pyarray<double> &KWs = args.array<double>("KWs");
355 double useMetrics = args.scalar<double>("useMetrics");
356 double alphaBDF = args.scalar<double>("alphaBDF");
357 int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
358 double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
359 double sc_uref = args.scalar<double>("sc_uref");
360 double sc_alpha = args.scalar<double>("sc_alpha");
361 xt::pyarray<int> &u_l2g = args.array<int>("u_l2g");
362 xt::pyarray<double> &elementDiameter = args.array<double>("elementDiameter");
363 xt::pyarray<double> &u_dof = args.array<double>("u_dof");
364 xt::pyarray<double> &u_dof_old = args.array<double>("u_dof_old");
365 xt::pyarray<double> &velocity = args.array<double>("velocity");
366 xt::pyarray<double> &q_m = args.array<double>("q_m");
367 xt::pyarray<double> &q_theta = args.array<double>("q_theta");
368 xt::pyarray<double> &q_u = args.array<double>("q_u");
369 xt::pyarray<double> &q_dV = args.array<double>("q_dV");
370 xt::pyarray<double> &q_m_betaBDF = args.array<double>("q_m_betaBDF");
371 xt::pyarray<double> &cfl = args.array<double>("cfl");
372 xt::pyarray<double> &q_numDiff_u = args.array<double>("q_numDiff_u");
373 xt::pyarray<double> &q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
374 int offset_u = args.scalar<int>("offset_u");
375 int stride_u = args.scalar<int>("stride_u");
376 xt::pyarray<double> &globalResidual = args.array<double>("globalResidual");
377 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
378 xt::pyarray<int> &exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
379 xt::pyarray<int> &elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
380 xt::pyarray<int> &elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
381 xt::pyarray<double> &ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
382 xt::pyarray<int> &isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
383 xt::pyarray<double> &ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
384 xt::pyarray<int> &isFluxBoundary_u = args.array<int>("isFluxBoundary_u");
385 xt::pyarray<double> &ebqe_bc_flux_ext = args.array<double>("ebqe_bc_flux_ext");
386 xt::pyarray<double> &ebqe_phi = args.array<double>("ebqe_phi");
387 double epsFact = args.scalar<double>("epsFact");
388 xt::pyarray<double> &ebqe_u = args.array<double>("ebqe_u");
389 xt::pyarray<double> &ebqe_theta = args.array<double>("ebqe_theta");
390 xt::pyarray<double> &ebqe_flux = args.array<double>("ebqe_flux");
391 // VMS
392 double VMS = args.scalar<double>("VMS");
393 // PARAMETERS FOR EDGE BASED STABILIZATION
394 double cE = args.scalar<double>("cE");
395 double cK = args.scalar<double>("cK");
396 // PARAMETERS FOR LOG BASED ENTROPY FUNCTION
397 double uL = args.scalar<double>("uL");
398 double uR = args.scalar<double>("uR");
399 // PARAMETERS FOR EDGE VISCOSITY
400 int numDOFs = args.scalar<int>("numDOFs");
401 int NNZ = args.scalar<int>("NNZ");
402 xt::pyarray<int> &csrRowIndeces_DofLoops = args.array<int>("csrRowIndeces_DofLoops");
403 xt::pyarray<int> &csrColumnOffsets_DofLoops = args.array<int>("csrColumnOffsets_DofLoops");
404 xt::pyarray<int> &csrRowIndeces_CellLoops = args.array<int>("csrRowIndeces_CellLoops");
405 xt::pyarray<int> &csrColumnOffsets_CellLoops = args.array<int>("csrColumnOffsets_CellLoops");
406 xt::pyarray<int> &csrColumnOffsets_eb_CellLoops = args.array<int>("csrColumnOffsets_eb_CellLoops");
407 // C matrices
408 xt::pyarray<double> &Cx = args.array<double>("Cx");
409 xt::pyarray<double> &Cy = args.array<double>("Cy");
410 xt::pyarray<double> &Cz = args.array<double>("Cz");
411 xt::pyarray<double> &CTx = args.array<double>("CTx");
412 xt::pyarray<double> &CTy = args.array<double>("CTy");
413 xt::pyarray<double> &CTz = args.array<double>("CTz");
414 xt::pyarray<double> &ML = args.array<double>("ML");
415 xt::pyarray<double> &delta_x_ij = args.array<double>("delta_x_ij");
416 // PARAMETERS FOR 1st or 2nd ORDER MPP METHOD
417 int LUMPED_MASS_MATRIX = args.scalar<int>("LUMPED_MASS_MATRIX");
418 STABILIZATION STABILIZATION_TYPE{static_cast<STABILIZATION>(args.scalar<int>("STABILIZATION_TYPE"))};
419 int ENTROPY_TYPE = args.scalar<int>("ENTROPY_TYPE");
420 // FOR FCT
421 xt::pyarray<double> &dLow = args.array<double>("dLow");
422 xt::pyarray<double> &fluxMatrix = args.array<double>("fluxMatrix");
423 // AUX QUANTITIES OF INTEREST
424 xt::pyarray<double> &quantDOFs = args.array<double>("quantDOFs");
425
426 assert(a_rowptr.data()[nSpace] == nnz);
427 assert(a_rowptr.data()[nSpace] == nSpace);
428 //cek should this be read in?
429 double Ct_sge = 4.0;
430
431 xt::pyarray<double> &anb_seepage_flux_n = args.array<double>("anb_seepage_flux_n");
432
433 xt::pyarray<double> &velocity_couple = args.array<double>("velocity_couple");
434 xt::pyarray<double> &ebqe_velocity_ext_couple = args.array<double>("ebqe_velocity_ext_couple");
435
436 // xt::pyarray<double> &q_x = args.array<double>("q_x");
437 // xt::pyarray<double> &ebqe_x = args.array<double>("ebqe_x");
438
439 //double anb_seepage_flux=0.0;
440 double &anb_seepage_flux(args.scalar<double>("anb_seepage_flux"));
441 xt::pyarray<double> &q_velocity = args.array<double>("q_velocity");
442 anb_seepage_flux = 0.0;
443
444 //loop over elements to compute volume integrals and load them into element and global residual
445 //
446 //eN is the element index
447 //eN_k is the quadrature point index for a scalar
448 //eN_k_nSpace is the quadrature point index for a vector
449 //eN_i is the element test function index
450 //eN_j is the element trial function index
451 //eN_k_j is the quadrature point index for a trial function
452 //eN_k_i is the quadrature point index for a trial function
453 for (int eN = 0; eN < nElements_global; eN++) {
454 //declare local storage for element residual and initialize
455 double elementResidual_u[nDOF_test_element];
456 for (int i = 0; i < nDOF_test_element; i++) { elementResidual_u[i] = 0.0; } //i
457 //loop over quadrature points and compute integrands
458 for (int k = 0; k < nQuadraturePoints_element; k++) {
459 //compute indeces and declare local storage
460 int eN_k = eN * nQuadraturePoints_element + k, eN_k_nSpace = eN_k * nSpace, eN_nDOF_trial_element = eN * nDOF_trial_element;
461 double u = 0.0, grad_u[nSpace], grad_u_old[nSpace], m = 0.0, dm = 0.0, f[nSpace], df[nSpace], a[nnz], da[nnz], as[nnz], m_t = 0.0, dm_t = 0.0, pdeResidual_u = 0.0, Lstar_u[nDOF_test_element], subgridError_u = 0.0, tau = 0.0, tau0 = 0.0, tau1 = 0.0, numDiff0 = 0.0, numDiff1 = 0.0, jac[nSpace * nSpace], jacDet, jacInv[nSpace * nSpace], u_grad_trial[nDOF_trial_element * nSpace], u_test_dV[nDOF_trial_element], u_grad_test_dV[nDOF_test_element * nSpace], dV, x, y, z, xt, yt, zt, G[nSpace * nSpace], G_dd_G, tr_G, norm_Rv;
462 //
463 //compute solution and gradients at quadrature points
464 //
465 ck.calculateMapping_element(eN, k, mesh_dof.data(), mesh_l2g.data(), mesh_trial_ref.data(), mesh_grad_trial_ref.data(), jac, jacDet, jacInv, x, y, z);
466 ck.calculateMappingVelocity_element(eN, k, mesh_velocity_dof.data(), mesh_l2g.data(), mesh_trial_ref.data(), xt, yt, zt);
467 //get the physical integration weight
468 dV = fabs(jacDet) * dV_ref.data()[k];
469 q_dV.data()[eN_k] = dV;
470 ck.calculateG(jacInv, G, G_dd_G, tr_G);
471 //get the trial function gradients
472 ck.gradTrialFromRef(&u_grad_trial_ref.data()[k * nDOF_trial_element * nSpace], jacInv, u_grad_trial);
473 //get the solution
474 ck.valFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element], &u_trial_ref.data()[k * nDOF_trial_element], u);
475 //get the solution gradients
476 ck.gradFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element], u_grad_trial, grad_u);
477
478 // //populate q_x
479 // const int eN_k_3d = eN_k * 3;
480 // q_x.data()[eN_k_3d + 0] = x;
481 // q_x.data()[eN_k_3d + 1] = y;
482 // q_x.data()[eN_k_3d + 2] = z;
483
484 //precalculate test function products with integration weights
485 for (int j = 0; j < nDOF_trial_element; j++) {
486 u_test_dV[j] = u_test_ref.data()[k * nDOF_trial_element + j] * dV;
487 for (int I = 0; I < nSpace; I++) {
488 u_grad_test_dV[j * nSpace + I] = u_grad_trial[j * nSpace + I] * dV; //cek warning won't work for Petrov-Galerkin
489 }
490 }
491 //
492 //calculate pde coefficients at quadrature points
493 //
494 double Kr, dKr, thetaW;
495 const double rho_local = q_rho.data()[eN_k];
496 const double rho_velocity = std::fabs(rho_local) > 1.0e-12 ? rho_local : rho;
497 evaluateCoefficients(a_rowptr.data(), a_colind.data(), rho, rho_local, beta, gravity.data(), alpha.data()[elementMaterialTypes.data()[eN]], n.data()[elementMaterialTypes.data()[eN]], thetaR.data()[elementMaterialTypes.data()[eN]],
498 thetaSR.data()[elementMaterialTypes.data()[eN]], &KWs.data()[elementMaterialTypes.data()[eN] * nnz], u, m, dm, f, df, a, da, as, Kr, dKr, thetaW);
499 q_theta.data()[eN_k] = thetaW;
500
501
502 for (int I = 0; I < nSpace; ++I) {
503 q_velocity.data()[eN_k_nSpace + I] = grad_u[I];
504 }
505 // Darcy Velocity
506 double pressure_gradient[nSpace];
507 const double rho_ratio = rho_velocity / rho;
508 for (int J=0; J<nSpace; ++J)
509 pressure_gradient[J] = grad_u[J] - rho_ratio * gravity.data()[J];
510 // for each row I, acc = sum_j (a_{Ij}/rho) * gp[j]
511 for (int I=0; I<nSpace; ++I) {
512 double acc = 0.0;
513 for (int ii = a_rowptr.data()[I]; ii < a_rowptr.data()[I+1]; ++ii) {
514 const int J = a_colind.data()[ii];
515 acc += (a[ii] / rho_velocity) * pressure_gradient[J];
516 }
517 velocity.data()[eN_k_nSpace + I] = -acc;
518 velocity_couple.data()[eN_k_nSpace + I] = -acc ;
519 }
520 //
521 //calculate time derivative at quadrature points
522 //
523 ck.bdf(alphaBDF, q_m_betaBDF.data()[eN_k], m, dm, m_t, dm_t);
524 //
525 //calculate subgrid error (strong residual and adjoint)
526 //
527 //calculate strong residual
528 pdeResidual_u = ck.Mass_strong(m_t) + ck.Advection_strong(df, grad_u);
529 //calculate adjoint
530 for (int i = 0; i < nDOF_test_element; i++) {
531 int i_nSpace = i * nSpace;
532 Lstar_u[i] = ck.Advection_adjoint(df, &u_grad_test_dV[i_nSpace]);
533 }
534 //calculate tau and tau*Res
535 calculateSubgridError_tau(elementDiameter[eN], dm_t, df, cfl[eN_k], tau0);
536 calculateSubgridError_tau(Ct_sge, G, dm_t, df, tau1, cfl[eN_k]);
537
538 tau = useMetrics * tau1 + (1.0 - useMetrics) * tau0;
539
540 subgridError_u = -tau * pdeResidual_u;
541 // //
542 // //calculate shock capturing diffusion
543 // //
544 // ck.calculateNumericalDiffusion(shockCapturingDiffusion, elementDiameter[eN], pdeResidual_u, grad_u, numDiff0);
545 // ck.calculateNumericalDiffusion(shockCapturingDiffusion, sc_uref, sc_alpha, G, G_dd_G, pdeResidual_u, grad_u, numDiff1);
546 // q_numDiff_u[eN_k] = useMetrics * numDiff1 + (1.0 - useMetrics) * numDiff0;
547 //
548 //update element residual
549 //
550 for (int i = 0; i < nDOF_test_element; i++) {
551 int eN_k_i = eN_k * nDOF_test_element + i, eN_k_i_nSpace = eN_k_i * nSpace, i_nSpace = i * nSpace;
552
553 elementResidual_u[i] += ck.Mass_weak(m_t, u_test_dV[i]) + ck.Advection_weak(f, &u_grad_test_dV[i_nSpace]) + ck.Diffusion_weak(a_rowptr.data(), a_colind.data(), a, grad_u, &u_grad_test_dV[i_nSpace]) + VMS * ck.SubgridError(subgridError_u, Lstar_u[i]) + VMS * ck.NumericalDiffusion(q_numDiff_u_last[eN_k], grad_u, &u_grad_test_dV[i_nSpace]);
554 } //i
555 //
556 q_m.data()[eN_k] = m;
557 q_u.data()[eN_k] = u;
558 }
559 //
560 //load element into global residual and save element residual
561 //
562 for (int i = 0; i < nDOF_test_element; i++) {
563 int eN_i = eN * nDOF_test_element + i;
564
565 globalResidual.data()[offset_u + stride_u * u_l2g.data()[eN_i]] += elementResidual_u[i];
566 } //i
567 } //elements
568 //
569 //loop over exterior element boundaries to calculate surface integrals and load into element and global residuals
570 //
571 //ebNE is the Exterior element boundary INdex
572 //ebN is the element boundary INdex
573 //eN is the element index
574 for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++) {
575 int ebN = exteriorElementBoundariesArray.data()[ebNE], eN = elementBoundaryElementsArray.data()[ebN * 2 + 0], ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN * 2 + 0], eN_nDOF_trial_element = eN * nDOF_trial_element;
576 double elementResidual_u[nDOF_test_element];
577 for (int i = 0; i < nDOF_test_element; i++) { elementResidual_u[i] = 0.0; }
578 for (int kb = 0; kb < nQuadraturePoints_elementBoundary; kb++) {
579 int ebNE_kb = ebNE * nQuadraturePoints_elementBoundary + kb, ebNE_kb_nSpace = ebNE_kb * nSpace, ebN_local_kb = ebN_local * nQuadraturePoints_elementBoundary + kb, ebN_local_kb_nSpace = ebN_local_kb * nSpace;
580 double u_ext = 0.0, grad_u_ext[nSpace], m_ext = 0.0, dm_ext = 0.0, f_ext[nSpace], df_ext[nSpace], a_ext[nnz], da_ext[nnz], as_ext[nnz], flux_ext = 0.0,
581 //anb_seepage_flux=0.0, // for flux calculation
582 bc_u_ext = 0.0, bc_grad_u_ext[nSpace], bc_m_ext = 0.0, bc_dm_ext = 0.0, bc_f_ext[nSpace], bc_df_ext[nSpace], bc_a_ext[nnz], bc_da_ext[nnz], bc_as_ext[nnz], jac_ext[nSpace * nSpace], jacDet_ext, jacInv_ext[nSpace * nSpace], boundaryJac[nSpace * (nSpace - 1)], metricTensor[(nSpace - 1) * (nSpace - 1)], metricTensorDetSqrt, dS, u_test_dS[nDOF_test_element], u_grad_trial_trace[nDOF_trial_element * nSpace], normal[3], x_ext, y_ext, z_ext, xt_ext, yt_ext, zt_ext, integralScaling, G[nSpace * nSpace], G_dd_G, tr_G;
583 //
584 //calculate the solution and gradients at quadrature points
585 //
586 //compute information about mapping from reference element to physical element
587 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,
588 normal_ref.data(), normal, x_ext, y_ext, z_ext);
589 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);
590 dS = ((1.0 - MOVING_DOMAIN) * metricTensorDetSqrt + MOVING_DOMAIN * integralScaling) * dS_ref.data()[kb];
591 //get the metric tensor
592 //cek todo use symmetry
593 ck.calculateG(jacInv_ext, G, G_dd_G, tr_G);
594 //compute shape and solution information
595 //shape
596 ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace * nDOF_trial_element], jacInv_ext, u_grad_trial_trace);
597 //solution and gradient
598 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);
599 ck.gradFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element], u_grad_trial_trace, grad_u_ext);
600
601 //populate ebqe_x
602 // const int ebNE_kb_3d = ebNE_kb * 3;
603 // ebqe_x.data()[ebNE_kb_3d + 0] = x_ext;
604 // ebqe_x.data()[ebNE_kb_3d + 1] = y_ext;
605 // ebqe_x.data()[ebNE_kb_3d + 2] = z_ext;
606
607 //precalculate test function products with integration weights
608 for (int j = 0; j < nDOF_trial_element; j++) { u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb * nDOF_test_element + j] * dS; }
609 //
610 //load the boundary values
611 //
612 bc_u_ext = isDOFBoundary_u.data()[ebNE_kb] * ebqe_bc_u_ext.data()[ebNE_kb] + (1 - isDOFBoundary_u.data()[ebNE_kb]) * u_ext;
613 //
614 //calculate the pde coefficients using the solution and the boundary values for the solution
615 //
616 const double rho_ext = ebqe_rho.data()[ebNE_kb];
617 const double rho_velocity_ext = std::fabs(rho_ext) > 1.0e-12 ? rho_ext : rho;
618 double Kr, dKr, thetaW_ext, thetaW_bc;
619 evaluateCoefficients(a_rowptr.data(), a_colind.data(), rho, rho_ext, beta, gravity.data(), alpha.data()[elementMaterialTypes.data()[eN]], n.data()[elementMaterialTypes.data()[eN]], thetaR.data()[elementMaterialTypes.data()[eN]],
620 thetaSR.data()[elementMaterialTypes.data()[eN]], &KWs.data()[elementMaterialTypes.data()[eN] * nnz], u_ext, m_ext, dm_ext, f_ext, df_ext, a_ext, da_ext, as_ext, Kr, dKr, thetaW_ext);
621 evaluateCoefficients(a_rowptr.data(), a_colind.data(), rho, rho_ext, beta, gravity.data(), alpha.data()[elementMaterialTypes.data()[eN]], n.data()[elementMaterialTypes.data()[eN]], thetaR.data()[elementMaterialTypes.data()[eN]],
622 thetaSR.data()[elementMaterialTypes.data()[eN]], &KWs.data()[elementMaterialTypes.data()[eN] * nnz], bc_u_ext, bc_m_ext, bc_dm_ext, bc_f_ext, bc_df_ext, bc_a_ext, bc_da_ext, bc_as_ext, Kr, dKr, thetaW_bc);
623 ebqe_theta.data()[ebNE_kb] = thetaW_ext;
624
625 //
626 //Calculate Darcy velocity on exterior face : v_ext = -(a_ext/rho) * (grad_u_ext + gravity) ---
627 //
628 double ext_pressure_gradient[nSpace];
629 const double rho_ratio_ext = rho_velocity_ext / rho;
630 for (int J=0; J<nSpace; ++J)
631 ext_pressure_gradient[J] = grad_u_ext[J] - rho_ratio_ext * gravity.data()[J];
632
633 for (int I=0; I<nSpace; ++I) {
634 double acc = 0.0;
635 for (int ii = a_rowptr.data()[I]; ii < a_rowptr.data()[I+1]; ++ii) {
636 const int J = a_colind.data()[ii];
637 acc += (a_ext[ii] / rho_velocity_ext) * ext_pressure_gradient[J];
638 }
639 ebqe_velocity_ext.data()[ebNE_kb_nSpace + I] = -acc;
640 ebqe_velocity_ext_couple.data()[ebNE_kb_nSpace + I] = -acc ; // store vector at this boundary qp
641 }
642
643
644 //
645 //calculate the numerical fluxes
646 //
647 exteriorNumericalFlux(ebqe_bc_flux_ext[ebNE_kb], a_rowptr.data(), a_colind.data(),
648 isSeepageFace.data()[ebNE], //tricky, this is a face flag not face quad
649 isDOFBoundary_u.data()[ebNE_kb], normal, bc_u_ext, a_ext, grad_u_ext, u_ext, f_ext,
650 ebqe_penalty_ext.data()[ebNE_kb], // penalty,
651 flux_ext);
652 ebqe_flux.data()[ebNE_kb] = flux_ext;
653
654 anb_seepage_flux = seepagefluxcalculator(anb_seepage_flux, isSeepageFace.data()[ebNE], dS, flux_ext);
655 anb_seepage_flux_n.data()[0] = anb_seepage_flux;
656 ebqe_u.data()[ebNE_kb] = u_ext;
657 //
658 //update residuals
659 //
660 for (int i = 0; i < nDOF_test_element; i++) {
661 elementResidual_u[i] += ck.ExteriorElementBoundaryFlux(flux_ext, u_test_dS[i]);
662 } //i
663 } //kb
664
665 //
666 //update the element and global residual storage
667 //
668 for (int i = 0; i < nDOF_test_element; i++) {
669 int eN_i = eN * nDOF_test_element + i;
670 globalResidual.data()[offset_u + stride_u * u_l2g.data()[eN_i]] += elementResidual_u[i];
671 } //i
672 } //ebNE
673 }
674
675 void calculateJacobian(arguments_dict &args)
676 {
677 xt::pyarray<double> &mesh_trial_ref = args.array<double>("mesh_trial_ref");
678 xt::pyarray<double> &mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
679 xt::pyarray<double> &mesh_dof = args.array<double>("mesh_dof");
680 xt::pyarray<double> &mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
681 double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
682 xt::pyarray<int> &mesh_l2g = args.array<int>("mesh_l2g");
683 xt::pyarray<double> &dV_ref = args.array<double>("dV_ref");
684 xt::pyarray<double> &u_trial_ref = args.array<double>("u_trial_ref");
685 xt::pyarray<double> &u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
686 xt::pyarray<double> &u_test_ref = args.array<double>("u_test_ref");
687 xt::pyarray<double> &u_grad_test_ref = args.array<double>("u_grad_test_ref");
688 xt::pyarray<double> &mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
689 xt::pyarray<double> &mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
690 xt::pyarray<double> &dS_ref = args.array<double>("dS_ref");
691 xt::pyarray<double> &u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
692 xt::pyarray<double> &u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
693 xt::pyarray<double> &u_test_trace_ref = args.array<double>("u_test_trace_ref");
694 xt::pyarray<double> &u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
695 xt::pyarray<double> &normal_ref = args.array<double>("normal_ref");
696 xt::pyarray<double> &boundaryJac_ref = args.array<double>("boundaryJac_ref");
697 int nElements_global = args.scalar<int>("nElements_global");
698 xt::pyarray<double> &ebqe_penalty_ext = args.array<double>("ebqe_penalty_ext");
699 xt::pyarray<int> &elementMaterialTypes = args.array<int>("elementMaterialTypes");
700 xt::pyarray<int> &isSeepageFace = args.array<int>("isSeepageFace");
701 xt::pyarray<int> &a_rowptr = args.array<int>("a_rowptr");
702 xt::pyarray<int> &a_colind = args.array<int>("a_colind");
703 double rho = args.scalar<double>("rho");
704 double beta = args.scalar<double>("beta");
705
707 xt::pyarray<double> &q_rho = args.array<double>("q_rho");
708 xt::pyarray<double> &ebqe_rho = args.array<double>("ebqe_rho");
710
711 xt::pyarray<double> &gravity = args.array<double>("gravity");
712 xt::pyarray<double> &alpha = args.array<double>("alpha");
713 xt::pyarray<double> &n = args.array<double>("n");
714 xt::pyarray<double> &thetaR = args.array<double>("thetaR");
715 xt::pyarray<double> &thetaSR = args.array<double>("thetaSR");
716 PSK_TYPE_member = args.scalar<int>("PSK_TYPE");
717 xt::pyarray<double> &KWs = args.array<double>("KWs");
718 double useMetrics = args.scalar<double>("useMetrics");
719 double alphaBDF = args.scalar<double>("alphaBDF");
720 int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
721 double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
722 // VMS
723 double VMS = args.scalar<double>("VMS");
724 xt::pyarray<int> &u_l2g = args.array<int>("u_l2g");
725 xt::pyarray<double> &elementDiameter = args.array<double>("elementDiameter");
726 xt::pyarray<double> &u_dof = args.array<double>("u_dof");
727 xt::pyarray<double> &velocity = args.array<double>("velocity");
728 xt::pyarray<double> &q_m_betaBDF = args.array<double>("q_m_betaBDF");
729 xt::pyarray<double> &cfl = args.array<double>("cfl");
730 xt::pyarray<double> &q_numDiff_u = args.array<double>("q_numDiff_u");
731 xt::pyarray<double> &q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
732 xt::pyarray<int> &csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
733 xt::pyarray<int> &csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
734 xt::pyarray<double> &globalJacobian = args.array<double>("globalJacobian");
735 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
736 xt::pyarray<int> &exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
737 xt::pyarray<int> &elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
738 xt::pyarray<int> &elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
739 xt::pyarray<double> &ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
740 xt::pyarray<int> &isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
741 xt::pyarray<double> &ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
742 xt::pyarray<int> &isFluxBoundary_u = args.array<int>("isFluxBoundary_u");
743 xt::pyarray<double> &ebqe_bc_flux_ext = args.array<double>("ebqe_bc_flux_ext");
744 xt::pyarray<int> &csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
745 int LUMPED_MASS_MATRIX = args.scalar<int>("LUMPED_MASS_MATRIX");
746 assert(a_rowptr.data()[nSpace] == nnz);
747 assert(a_rowptr.data()[nSpace] == nSpace);
748 double Ct_sge = 4.0;
749
750 //
751 //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
752 //
753 for (int eN = 0; eN < nElements_global; eN++) {
754 double elementJacobian_u_u[nDOF_test_element][nDOF_trial_element];
755 for (int i = 0; i < nDOF_test_element; i++) {
756 for (int j = 0; j < nDOF_trial_element; j++) { elementJacobian_u_u[i][j] = 0.0; }
757 }
758 for (int k = 0; k < nQuadraturePoints_element; k++) {
759 int eN_k = eN * nQuadraturePoints_element + k, //index to a scalar at a quadrature point
760 eN_k_nSpace = eN_k * nSpace,
761 eN_nDOF_trial_element = eN * nDOF_trial_element; //index to a vector at a quadrature point
762
763 //declare local storage
764 double u = 0.0, grad_u[nSpace], m = 0.0, dm = 0.0, f[nSpace], df[nSpace], a[nnz], da[nnz], as[nnz], m_t = 0.0, dm_t = 0.0, dpdeResidual_u_u[nDOF_trial_element], Lstar_u[nDOF_test_element], dsubgridError_u_u[nDOF_trial_element], tau = 0.0, tau0 = 0.0, tau1 = 0.0, jac[nSpace * nSpace], jacDet, jacInv[nSpace * nSpace], u_grad_trial[nDOF_trial_element * nSpace], dV, u_test_dV[nDOF_test_element], u_grad_test_dV[nDOF_test_element * nSpace], x, y, z, xt, yt, zt, G[nSpace * nSpace], G_dd_G, tr_G;
765 //
766 //calculate solution and gradients at quadrature points
767 //
768 //get jacobian, etc for mapping reference element
769 ck.calculateMapping_element(eN, k, mesh_dof.data(), mesh_l2g.data(), mesh_trial_ref.data(), mesh_grad_trial_ref.data(), jac, jacDet, jacInv, x, y, z);
770 ck.calculateMappingVelocity_element(eN, k, mesh_velocity_dof.data(), mesh_l2g.data(), mesh_trial_ref.data(), xt, yt, zt);
771 //get the physical integration weight
772 dV = fabs(jacDet) * dV_ref.data()[k];
773 ck.calculateG(jacInv, G, G_dd_G, tr_G);
774 //get the trial function gradients
775 ck.gradTrialFromRef(&u_grad_trial_ref.data()[k * nDOF_trial_element * nSpace], jacInv, u_grad_trial);
776 //get the solution
777 ck.valFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element], &u_trial_ref.data()[k * nDOF_trial_element], u);
778 //get the solution gradients
779 ck.gradFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element], u_grad_trial, grad_u);
780 //precalculate test function products with integration weights
781 for (int j = 0; j < nDOF_trial_element; j++) {
782 u_test_dV[j] = u_test_ref.data()[k * nDOF_trial_element + j] * dV;
783 for (int I = 0; I < nSpace; I++) {
784 u_grad_test_dV[j * nSpace + I] = u_grad_trial[j * nSpace + I] * dV; //cek warning won't work for Petrov-Galerkin
785 }
786 }
787 //
788 //calculate pde coefficients and derivatives at quadrature points
789 //
790 double Kr, dKr, thetaW;
791 //const double rho_local = q_rho.data()[eN_k];
792
793 evaluateCoefficients(a_rowptr.data(), a_colind.data(), rho, q_rho.data()[eN_k], beta, gravity.data(), alpha.data()[elementMaterialTypes.data()[eN]], n.data()[elementMaterialTypes.data()[eN]], thetaR.data()[elementMaterialTypes.data()[eN]],
794 thetaSR.data()[elementMaterialTypes.data()[eN]], &KWs.data()[elementMaterialTypes.data()[eN] * nnz], u, m, dm, f, df, a, da, as, Kr, dKr, thetaW);
795 //
796 //calculate time derivatives
797 //
798 ck.bdf(alphaBDF, q_m_betaBDF.data()[eN_k], m, dm, m_t, dm_t);
799 //
800 //calculate subgrid error contribution to the Jacobian (strong residual, adjoint, jacobian of strong residual)
801 //
802 //calculate the adjoint times the test functions
803 for (int i = 0; i < nDOF_test_element; i++) {
804 int i_nSpace = i * nSpace;
805 Lstar_u[i] = ck.Advection_adjoint(df, &u_grad_test_dV[i_nSpace]);
806 }
807 //calculate the Jacobian of strong residual
808 for (int j = 0; j < nDOF_trial_element; j++) {
809 int j_nSpace = j * nSpace;
810 dpdeResidual_u_u[j] = ck.MassJacobian_strong(dm_t, u_trial_ref[k * nDOF_trial_element + j]) + ck.AdvectionJacobian_strong(df, &u_grad_trial[j_nSpace]);
811 }
812 //tau and tau*Res
813 calculateSubgridError_tau(elementDiameter[eN], dm_t, df, cfl[eN_k], tau0);
814 calculateSubgridError_tau(Ct_sge, G, dm_t, df, tau1, cfl[eN_k]);
815 tau = useMetrics * tau1 + (1.0 - useMetrics) * tau0;
816 for (int j = 0; j < nDOF_trial_element; j++) dsubgridError_u_u[j] = -tau * dpdeResidual_u_u[j];
817 for (int i = 0; i < nDOF_test_element; i++) {
818 for (int j = 0; j < nDOF_trial_element; j++) {
819 int j_nSpace = j * nSpace;
820 int i_nSpace = i * nSpace;
821 elementJacobian_u_u[i][j] += ck.MassJacobian_weak(dm_t, u_trial_ref.data()[k * nDOF_trial_element + j], u_test_dV[i]) + ck.AdvectionJacobian_weak(df, u_trial_ref.data()[k * nDOF_trial_element + j], &u_grad_test_dV[i_nSpace]) +
822 ck.DiffusionJacobian_weak(a_rowptr.data(), a_colind.data(), a, da, grad_u, &u_grad_test_dV[i_nSpace], 1.0, u_trial_ref.data()[k * nDOF_trial_element + j], &u_grad_trial[j_nSpace]) + VMS * ck.SubgridErrorJacobian(dsubgridError_u_u[j], Lstar_u[i]) + VMS * ck.NumericalDiffusionJacobian(q_numDiff_u_last[eN_k], &u_grad_trial[j_nSpace], &u_grad_test_dV[i_nSpace]);
823 } //j
824 } //i
825 } //k
826 //
827 //load into element Jacobian into global Jacobian
828 //
829 for (int i = 0; i < nDOF_test_element; i++) {
830 int eN_i = eN * nDOF_test_element + i;
831 for (int j = 0; j < nDOF_trial_element; j++) {
832 int eN_i_j = eN_i * nDOF_trial_element + j;
833 globalJacobian.data()[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_u_u[eN_i_j]] += elementJacobian_u_u[i][j];
834 } //j
835 } //i
836 } //elements
837 //
838 //loop over exterior element boundaries to compute the surface integrals and load them into the global Jacobian
839 //
840 for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++) {
841 int ebN = exteriorElementBoundariesArray.data()[ebNE];
842 int eN = elementBoundaryElementsArray.data()[ebN * 2 + 0], ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN * 2 + 0], eN_nDOF_trial_element = eN * nDOF_trial_element;
843 for (int kb = 0; kb < nQuadraturePoints_elementBoundary; kb++) {
844 int ebNE_kb = ebNE * nQuadraturePoints_elementBoundary + kb, ebNE_kb_nSpace = ebNE_kb * nSpace, ebN_local_kb = ebN_local * nQuadraturePoints_elementBoundary + kb, ebN_local_kb_nSpace = ebN_local_kb * nSpace;
845
846 double u_ext = 0.0, grad_u_ext[nSpace], m_ext = 0.0, dm_ext = 0.0, f_ext[nSpace], df_ext[nSpace], a_ext[nnz], da_ext[nnz], as_ext[nnz], dflux_u_u_ext = 0.0, bc_u_ext = 0.0,
847 //bc_grad_u_ext[nSpace],
848 bc_m_ext = 0.0, bc_dm_ext = 0.0, bc_f_ext[nSpace], bc_df_ext[nSpace], bc_a_ext[nnz], bc_da_ext[nnz], bc_as_ext[nnz], fluxJacobian_u_u[nDOF_trial_element], jac_ext[nSpace * nSpace], jacDet_ext, jacInv_ext[nSpace * nSpace], boundaryJac[nSpace * (nSpace - 1)], metricTensor[(nSpace - 1) * (nSpace - 1)], metricTensorDetSqrt, dS, u_test_dS[nDOF_test_element], u_grad_trial_trace[nDOF_trial_element * nSpace], normal[3], x_ext, y_ext, z_ext, xt_ext, yt_ext, zt_ext, integralScaling, G[nSpace * nSpace], G_dd_G, tr_G;
849 //
850 //calculate the solution and gradients at quadrature points
851 //
852 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,
853 normal_ref.data(), normal, x_ext, y_ext, z_ext);
854 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);
855 dS = ((1.0 - MOVING_DOMAIN) * metricTensorDetSqrt + MOVING_DOMAIN * integralScaling) * dS_ref.data()[kb];
856 ck.calculateG(jacInv_ext, G, G_dd_G, tr_G);
857 //compute shape and solution information
858 //shape
859 ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace * nDOF_trial_element], jacInv_ext, u_grad_trial_trace);
860 //solution and gradients
861 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);
862 ck.gradFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element], u_grad_trial_trace, grad_u_ext);
863 //precalculate test function products with integration weights
864 for (int j = 0; j < nDOF_trial_element; j++) { u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb * nDOF_test_element + j] * dS; }
865 //
866 //load the boundary values
867 //
868 bc_u_ext = isDOFBoundary_u.data()[ebNE_kb] * ebqe_bc_u_ext.data()[ebNE_kb] + (1 - isDOFBoundary_u.data()[ebNE_kb]) * u_ext;
869 //
870 //calculate the internal and external trace of the pde coefficients
871 //
872 double Kr, dKr, thetaW, thetaW_bc;
873 const double rho_ext = ebqe_rho.data()[ebNE_kb];
874
875 evaluateCoefficients(a_rowptr.data(), a_colind.data(), rho, rho_ext, beta, gravity.data(), alpha.data()[elementMaterialTypes.data()[eN]], n.data()[elementMaterialTypes.data()[eN]], thetaR.data()[elementMaterialTypes.data()[eN]],
876 thetaSR.data()[elementMaterialTypes.data()[eN]], &KWs.data()[elementMaterialTypes.data()[eN] * nnz], u_ext, m_ext, dm_ext, f_ext, df_ext, a_ext, da_ext, as_ext, Kr, dKr, thetaW);
877 evaluateCoefficients(a_rowptr.data(), a_colind.data(), rho, rho_ext, beta, gravity.data(), alpha.data()[elementMaterialTypes.data()[eN]], n.data()[elementMaterialTypes.data()[eN]], thetaR.data()[elementMaterialTypes.data()[eN]],
878 thetaSR.data()[elementMaterialTypes.data()[eN]], &KWs.data()[elementMaterialTypes.data()[eN] * nnz], bc_u_ext, bc_m_ext, bc_dm_ext, bc_f_ext, bc_df_ext, bc_a_ext, bc_da_ext, bc_as_ext, Kr, dKr, thetaW_bc);
879 //
880 //calculate the flux jacobian
881 //
882 for (int j = 0; j < nDOF_trial_element; j++) {
883 exteriorNumericalFluxJacobian(a_rowptr.data(), a_colind.data(), isDOFBoundary_u.data()[ebNE_kb], normal, a_ext, da_ext, grad_u_ext, &u_grad_trial_trace[j * nSpace], df_ext, u_trial_trace_ref.data()[ebN_local_kb * nDOF_test_element + j],
884 ebqe_penalty_ext.data()[ebNE_kb], //penalty,
885 fluxJacobian_u_u[j]);
886 } //j
887 //
888 //update the global Jacobian from the flux Jacobian
889 //
890 for (int i = 0; i < nDOF_test_element; i++) {
891 int eN_i = eN * nDOF_test_element + i;
892 for (int j = 0; j < nDOF_trial_element; j++) {
893 int ebN_i_j = ebN * 4 * nDOF_test_X_trial_element + i * nDOF_trial_element + j;
894 globalJacobian.data()[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] += fluxJacobian_u_u[j] * u_test_dS[i];
895 } //j
896 } //i
897 } //kb
898 } //ebNE
899 } //computeJacobian
900
901 void FCTStep(arguments_dict &args)
902{
903 xt::pyarray<double> &bc_mask = args.array<double>("bc_mask");
904 int NNZ = args.scalar<int>("NNZ"); // number of non-zero entries
905 int numDOFs = args.scalar<int>("numDOFs"); // number of DOFs
906 double dt = args.scalar<double>("dt");
907 xt::pyarray<double> &ML = args.array<double>("ML"); // lumped mass matrix (as vector)
908 xt::pyarray<double> &mn = args.array<double>("mn"); // DOFs at time tn
909 xt::pyarray<double> &mHigh = args.array<double>("mHigh"); // high-order mass at t^{n+1}
910 xt::pyarray<double> &mLow = args.array<double>("mLow"); // low-order mass at t^{n+1}
911 xt::pyarray<double> &mDotLow = args.array<double>("mDotLow");
912 xt::pyarray<double> &limited_solution = args.array<double>("limited_solution");
913 xt::pyarray<int> &csrRowIndeces_DofLoops = args.array<int>("csrRowIndeces_DofLoops");
914 xt::pyarray<int> &csrColumnOffsets_DofLoops = args.array<int>("csrColumnOffsets_DofLoops");
915 xt::pyarray<double> &MC = args.array<double>("MC"); // consistent mass matrix
916 xt::pyarray<double> &dt_times_fH_minus_fL = args.array<double>("dt_times_fH_minus_fL");
917 xt::pyarray<double> &min_m_bc = args.array<double>("min_m_bc");
918 xt::pyarray<double> &max_m_bc = args.array<double>("max_m_bc");
919 xt::pyarray<double> &fluxCorrection = args.array<double>("fluxCorrection");
920 //owned by Python so they survive between the two passes below
921 xt::pyarray<double> &Rpos = args.array<double>("Rpos");
922 xt::pyarray<double> &Rneg = args.array<double>("Rneg");
923 xt::pyarray<double> &FluxCorrectionMatrix = args.array<double>("FluxCorrectionMatrix");
924 // flags
925 int LUMPED_MASS_MATRIX = args.scalar<int>("LUMPED_MASS_MATRIX");
926 int MONOLITHIC = args.scalar<int>("MONOLITHIC");
927 int fct_pass = args.scalar<int>("fct_pass");
928
929 //PARALLEL. The limiter is L_ij = min(Rpos_i, Rneg_j), so an owned row i reads the ratio of a column
930 //j that may be a ghost, and everything a ratio is built from -- the row's own sparsity, ML, mDotLow,
931 //min/max_m_bc -- is INCOMPLETE at a ghost DOF: with one layer of overlap a ghost's element star and
932 //its boundary faces are only partly on this rank. Both ranks sharing a cut edge would then apply a
933 //different L_ij, f_ij = -f_ji dies and the correction manufactures mass along the partition.
934 //The cure is to let the owners' values cross the cut in the middle of the limiter:
935 // fct_pass = 1 -> FluxCorrectionMatrix + the local ratios Rpos/Rneg
936 // [Python forward-inserts Rpos/Rneg, owner -> ghost; the inputs above are scattered before pass 1]
937 // fct_pass = 2 -> apply the limiter
938 // fct_pass = 0 -> both back to back; this is the serial path and is byte-identical to the old
939 // single-pass function.
940 //mDot is pointwise in (mLow, mn) so it is rebuilt in whichever pass needs it rather than carried.
941 const bool doPass1 = (fct_pass != 2);
942 const bool doPass2 = (fct_pass != 1);
943
944 std::vector<double> mDot(numDOFs, 0.0);
945 for (int i = 0; i < numDOFs; i++)
946 mDot.at(i) = (mLow.at(i) - mn.at(i)) / dt; // local time derivative from low-order mass
947
949 // PASS 1: antidiffusive fluxes and Zalesak's ratios //
951 if (doPass1) {
952 int ij = 0;
953 for (int i = 0; i < numDOFs; i++) {
954
955 // initialize local min/max from BC
956 double mini = min_m_bc.at(i);
957 double maxi = max_m_bc.at(i);
958
959 double Pposi = 0.0, Pnegi = 0.0;
960
961 // LOOP OVER THE SPARSITY PATTERN (j-LOOP)
962 for (int offset = csrRowIndeces_DofLoops.at(i); offset < csrRowIndeces_DofLoops.at(i + 1);offset++)
963 {
964 int j = csrColumnOffsets_DofLoops.at(offset);
965
967 // COMPUTE THE BOUNDS //
969 if (GLOBAL_FCT == 0) {
970 if (MONOLITHIC == 0) {
971 mini = fmin(mini, mLow.at(j));
972 maxi = fmax(maxi, mLow.at(j));
973 } else {
974 mini = fmin(mini, mn.at(j));
975 maxi = fmax(maxi, mn.at(j));
976 }
977 }
978
979 if (MONOLITHIC == 0) {
980 FluxCorrectionMatrix.at(ij) = (LUMPED_MASS_MATRIX == 1 ? 0. : 1.) * dt * MC.at(ij) * (mDotLow.at(i) - mDotLow.at(j)) + dt_times_fH_minus_fL.at(ij);
981 } else {
982 FluxCorrectionMatrix.at(ij) = dt_times_fH_minus_fL.at(ij);
983 }
984
986 // COMPUTE P VECTORS //
988 Pposi += FluxCorrectionMatrix.at(ij) * ((FluxCorrectionMatrix.at(ij) > 0) ? 1. : 0.);
989 Pnegi += FluxCorrectionMatrix.at(ij) * ((FluxCorrectionMatrix.at(ij) < 0) ? 1. : 0.);
990
991 // update ij
992 ij += 1;
993 } // j-loop
994
996 // COMPUTE Q VECTORS //
998 double gamma;
999 double Qposi;
1000 double Qnegi;
1001
1002 if (MONOLITHIC == 0) {
1003 Qposi = ML.at(i) * (maxi - mLow.at(i));
1004 Qnegi = ML.at(i) * (mini - mLow.at(i));
1005 } else {
1006 // cek todo: don't think this is right for Richards
1007 gamma = 10.0 * ML.at(i);
1008 Qposi = fmin(0.5 * ML.at(i) * (1.0 - mn.at(i)), gamma * (maxi - mn.at(i)));
1009 Qnegi = fmax(0.5 * ML.at(i) * (0.0 - mn.at(i)), gamma * (mini - mn.at(i)));
1010 }
1011
1013 // COMPUTE R VECTORS //
1015 Rpos.at(i) = ((Pposi == 0.0) ? 1.0 : fmin(1.0, Qposi / Pposi));
1016 Rneg.at(i) = ((Pnegi == 0.0) ? 1.0 : fmin(1.0, Qnegi / Pnegi));
1017 } // i DOFs
1018 } // pass 1
1019
1021 // COMPUTE LIMITERS //
1023 if (doPass2) {
1024 int ij = 0;
1025 // std::cout << "FCT: entering second DOF loop (applying limiters)...\n";
1026
1027 for (int i = 0; i < numDOFs; i++) {
1028 double ith_Limiter_times_FluxCorrectionMatrix = 0.0;
1029 double alpha_fA, alpha_dot, beta_ij = 1.0;
1030 // LOOP OVER THE SPARSITY PATTERN (j-LOOP)
1031 for (int offset = csrRowIndeces_DofLoops.at(i); offset < csrRowIndeces_DofLoops.at(i + 1); offset++)
1032 {
1033 int j = csrColumnOffsets_DofLoops.at(offset);
1034 alpha_fA = ((FluxCorrectionMatrix.at(ij) > 0.0) ? fmin(Rpos.at(i), Rneg.at(j)) : fmin(Rneg.at(i), Rpos.at(j))) * FluxCorrectionMatrix.at(ij);
1035 alpha_dot = fmin(1.0, beta_ij * fabs(alpha_fA) / MC.at(ij) / fmax(1.0e-8, fabs(mDot.at(i) - mDot.at(j))));
1036
1037 if (MONOLITHIC == 0) {
1038 ith_Limiter_times_FluxCorrectionMatrix += alpha_fA;
1039 } else {
1040 ith_Limiter_times_FluxCorrectionMatrix += alpha_fA + (LUMPED_MASS_MATRIX == 1 ? 0. : 1.) * dt * alpha_dot * MC.at(ij) * (mDot.at(i) - mDot.at(j));
1041 }
1042 ij += 1;
1043 }
1044 fluxCorrection.at(i) = -ith_Limiter_times_FluxCorrectionMatrix * bc_mask.at(i) / dt;
1045 limited_solution.at(i) = mLow.at(i) + 1.0 / ML.at(i) * ith_Limiter_times_FluxCorrectionMatrix * bc_mask.at(i);
1046 }
1047 }
1048}
1049
1050
1051 void kth_FCT_step(arguments_dict &args)
1052 {
1053 int NNZ = args.scalar<int>("NNZ"); //number on non-zero entries on sparsity pattern
1054 int numDOFs = args.scalar<int>("numDOFs"); //number of DOFs
1055 int num_fct_iter = args.scalar<int>("num_fct_iter");
1056 double dt = args.scalar<double>("dt");
1057 xt::pyarray<double> &lumped_mass_matrix = args.array<double>("lumped_mass_matrix"); //lumped mass matrix (as vector)
1058 xt::pyarray<double> &soln = args.array<double>("soln"); //DOFs of solution at time tn
1059 xt::pyarray<double> &pn = args.array<double>("pn"); //DOFs of solution at time tn
1060 xt::pyarray<double> &solH = args.array<double>("solH"); //DOFs of high order solution at tnp1
1061 xt::pyarray<double> &uLow = args.array<double>("uLow");
1062 xt::pyarray<double> &uDotLow = args.array<double>("uDotLow");
1063 xt::pyarray<double> &dLow = args.array<double>("dLow");
1064 xt::pyarray<double> &solLim = args.array<double>("limited_solution");
1065 xt::pyarray<double> &MC = args.array<double>("MC");
1066 xt::pyarray<double> &ML = args.array<double>("ML");
1067 xt::pyarray<double> &FluxMatrix = args.array<double>("FluxMatrix");
1068 xt::pyarray<double> &limitedFlux = args.array<double>("limited_Flux");
1069 xt::pyarray<int> &csrRowIndeces_DofLoops = args.array<int>("csrRowIndeces_DofLoops"); //csr row indeces
1070 xt::pyarray<int> &csrColumnOffsets_DofLoops = args.array<int>("csrColumnOffsets_DofLoops"); //csr column offsets
1071 xt::pyarray<double> &MassMatrix = args.array<double>("MassMatrix"); //mass matrix
1072 xt::pyarray<double> &dt_times_fH_minus_fL = args.array<double>("dt_times_fH_minus_fL"); //low minus high order dissipative matrices
1073 xt::pyarray<double> &min_m_bc = args.array<double>("min_m_bc"); //min/max value at BCs. If DOF is not at boundary then min=1E10, max=-1E10
1074 xt::pyarray<double> &max_m_bc = args.array<double>("max_m_bc");
1075 int LUMPED_MASS_MATRIX = args.scalar<int>("LUMPED_MASS_MATRIX");
1076 int MONOLITHIC = args.scalar<int>("MONOLITHIC");
1077 double Rpos[numDOFs], Rneg[numDOFs];
1078 int ij = 0;
1079
1081 // ********** COMPUTE LOW ORDER SOLUTION ********** //
1083 if (num_fct_iter == 0) { // No FCT for global bounds
1084 for (int i = 0; i < numDOFs; i++) { solLim.data()[i] = uLow.data()[i]; }
1085 } else // do FCT iterations (with global bounds) on low order solution
1086 {
1087 for (int iter = 0; iter < num_fct_iter; iter++) {
1088 ij = 0;
1089 for (int i = 0; i < numDOFs; i++) {
1090 double maxi = 1.0, Pposi = 0;
1091 for (int offset = csrRowIndeces_DofLoops.data()[i]; offset < csrRowIndeces_DofLoops.data()[i + 1]; offset++) {
1092 int j = csrColumnOffsets_DofLoops.data()[offset];
1093 // compute Flux correction
1094 double Fluxij = FluxMatrix.data()[ij] - limitedFlux.data()[ij];
1095 Pposi += Fluxij * ((Fluxij > 0) ? 1. : 0.);
1096 // update ij
1097 ij += 1;
1098 }
1099 // compute Q vectors
1100 double mi = ML.data()[i];
1101 double solLimi = solLim.data()[i];
1102 double Qposi = mi * (maxi - solLimi);
1103 // compute R vectors
1104 Rpos[i] = ((Pposi == 0) ? 1. : fmin(1.0, Qposi / Pposi));
1105 }
1106 ij = 0;
1107 for (int i = 0; i < numDOFs; i++) {
1108 double ith_Limiter_times_FluxCorrectionMatrix = 0.;
1109 double Rposi = Rpos[i];
1110 for (int offset = csrRowIndeces_DofLoops.data()[i]; offset < csrRowIndeces_DofLoops.data()[i + 1]; offset++) {
1111 int j = csrColumnOffsets_DofLoops.data()[offset];
1112 // Flux Correction
1113 double Fluxij = FluxMatrix.data()[ij] - limitedFlux.data()[ij];
1114 // compute limiter
1115 double Lij = 1.0;
1116 Lij = (Fluxij > 0 ? Rposi : Rpos[j]);
1117 // compute limited flux
1118 ith_Limiter_times_FluxCorrectionMatrix += Lij * Fluxij;
1119
1120 // update limited flux
1121 limitedFlux.data()[ij] = Lij * Fluxij;
1122
1123 //update FluxMatrix
1124 FluxMatrix.data()[ij] = Fluxij;
1125
1126 //update ij
1127 ij += 1;
1128 }
1129 //update limited solution
1130 double mi = ML.data()[i];
1131 }
1132 }
1133 }
1134
1135 // ***************************************** //
1136 // ********** HIGH ORDER SOLUTION ********** //
1137 // ***************************************** //
1138 ij = 0;
1139 for (int i = 0; i < numDOFs; i++) {
1140 double mini = soln.data()[i], maxi = soln.data()[i];
1141 double Pposi = 0, Pnegi = 0.;
1142 for (int offset = csrRowIndeces_DofLoops.data()[i]; offset < csrRowIndeces_DofLoops.data()[i + 1]; offset++) {
1143 int j = csrColumnOffsets_DofLoops.data()[offset];
1144 // compute local bounds //
1145 mini = fmin(mini, soln.data()[j]);
1146 maxi = fmax(maxi, soln.data()[j]);
1147 // compute P vectors //
1148 double fij = (MC.data()[ij] * (uDotLow.data()[i] - uDotLow.data()[j]) / dt + dLow.data()[ij] * (uLow.data()[i] - uLow.data()[j]));
1149 Pposi += fij * (fij > 0 ? 1. : 0.);
1150 Pnegi += fij * (fij < 0 ? 1. : 0.);
1151 //update ij
1152 ij += 1;
1153 }
1154 // compute Q vectors //
1155 double mi = ML.data()[i];
1156 double Qposi = mi * (maxi - solLim.data()[i]);
1157 double Qnegi = mi * (mini - solLim.data()[i]);
1158 // compute R vectors //
1159 Rpos[i] = ((Pposi == 0) ? 1. : fmin(1.0, Qposi / Pposi));
1160 Rneg[i] = ((Pnegi == 0) ? 1. : fmin(1.0, Qnegi / Pnegi));
1161 }
1162
1163 // COMPUTE LIMITERS //
1164 ij = 0;
1165 for (int i = 0; i < numDOFs; i++) {
1166 double ith_limited_flux_correction = 0;
1167 double Rposi = Rpos[i];
1168 double Rnegi = Rneg[i];
1169 for (int offset = csrRowIndeces_DofLoops.data()[i]; offset < csrRowIndeces_DofLoops.data()[i + 1]; offset++) {
1170 int j = csrColumnOffsets_DofLoops.data()[offset];
1171 // compute flux correction
1172 double fij = (MC.data()[ij] * (uDotLow.data()[i] - uDotLow.data()[j]) / dt + dLow.data()[ij] * (uLow.data()[i] - uLow.data()[j]));
1173
1174 // compute limiters
1175 double Lij = 1.0;
1176 Lij = fij > 0 ? fmin(Rposi, Rneg[j]) : fmin(Rnegi, Rpos[j]);
1177 // compute ith_limited_flux_correction
1178 ith_limited_flux_correction += Lij * fij;
1179 ij += 1;
1180 }
1181 double mi = ML.data()[i];
1182 solLim[i] += 1. / mi * ith_limited_flux_correction;
1183 }
1184 }
1185
1186 void calculateResidual_entropy_viscosity(arguments_dict &args)
1187 {
1188 xt::pyarray<double> &globalJacobian = args.array<double>("globalJacobian");
1189 double Theta = args.scalar<double>("Theta");
1190 double Theta_h = args.scalar<double>("Theta_h");
1191 xt::pyarray<double> &bc_mask = args.array<double>("bc_mask");
1192 double dt = args.scalar<double>("dt");
1193 xt::pyarray<double> &mesh_trial_ref = args.array<double>("mesh_trial_ref");
1194 xt::pyarray<double> &mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1195 xt::pyarray<double> &mesh_dof = args.array<double>("mesh_dof");
1196 xt::pyarray<double> &mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
1197 double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
1198 xt::pyarray<int> &mesh_l2g = args.array<int>("mesh_l2g");
1199 xt::pyarray<double> &dV_ref = args.array<double>("dV_ref");
1200 xt::pyarray<double> &u_trial_ref = args.array<double>("u_trial_ref");
1201 xt::pyarray<double> &u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1202 xt::pyarray<double> &u_test_ref = args.array<double>("u_test_ref");
1203 xt::pyarray<double> &u_grad_test_ref = args.array<double>("u_grad_test_ref");
1204 xt::pyarray<double> &mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
1205 xt::pyarray<double> &mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
1206 xt::pyarray<double> &dS_ref = args.array<double>("dS_ref");
1207 xt::pyarray<double> &u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
1208
1209 xt::pyarray<double> &u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
1210 xt::pyarray<double> &u_test_trace_ref = args.array<double>("u_test_trace_ref");
1211 xt::pyarray<double> &u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
1212 xt::pyarray<double> &normal_ref = args.array<double>("normal_ref");
1213 xt::pyarray<double> &boundaryJac_ref = args.array<double>("boundaryJac_ref");
1214 int nElements_global = args.scalar<int>("nElements_global");
1215 xt::pyarray<double> &ebqe_penalty_ext = args.array<double>("ebqe_penalty_ext");
1216 xt::pyarray<int> &elementMaterialTypes = args.array<int>("elementMaterialTypes");
1217 xt::pyarray<int> &isSeepageFace = args.array<int>("isSeepageFace");
1218 xt::pyarray<int> &a_rowptr = args.array<int>("a_rowptr");
1219 xt::pyarray<int> &a_colind = args.array<int>("a_colind");
1220 double rho = args.scalar<double>("rho");
1221 double beta = args.scalar<double>("beta");
1223 xt::pyarray<double> &q_rho = args.array<double>("q_rho");
1224 xt::pyarray<double> &ebqe_rho = args.array<double>("ebqe_rho");
1225 xt::pyarray<double> &gravity = args.array<double>("gravity");
1226 xt::pyarray<double> &alpha = args.array<double>("alpha");
1227 xt::pyarray<double> &n = args.array<double>("n");
1228 xt::pyarray<double> &thetaR = args.array<double>("thetaR");
1229 xt::pyarray<double> &thetaSR = args.array<double>("thetaSR");
1230 PSK_TYPE_member = args.scalar<int>("PSK_TYPE");
1231 xt::pyarray<double> &KWs = args.array<double>("KWs");
1232 double useMetrics = args.scalar<double>("useMetrics");
1233 double alphaBDF = args.scalar<double>("alphaBDF");
1234 int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
1235 double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
1236 double sc_uref = args.scalar<double>("sc_uref");
1237 double sc_alpha = args.scalar<double>("sc_alpha");
1238 xt::pyarray<int> &u_l2g = args.array<int>("u_l2g");
1239 xt::pyarray<int> &r_l2g = args.array<int>("r_l2g");
1240 xt::pyarray<double> &elementDiameter = args.array<double>("elementDiameter");
1241 int degree_polynomial = args.scalar<int>("degree_polynomial");
1242 xt::pyarray<double> &u_dof = args.array<double>("u_dof");
1243 xt::pyarray<double> &u_dof_old = args.array<double>("u_dof_old");
1244 xt::pyarray<double> &velocity = args.array<double>("velocity");
1245 xt::pyarray<double> &q_m = args.array<double>("q_m");
1246 xt::pyarray<double> &q_theta = args.array<double>("q_theta");
1247 xt::pyarray<double> &q_u = args.array<double>("q_u");
1248 xt::pyarray<double> &q_dV = args.array<double>("q_dV");
1249 xt::pyarray<double> &q_m_betaBDF = args.array<double>("q_m_betaBDF");
1250 xt::pyarray<double> &cfl = args.array<double>("cfl");
1251 xt::pyarray<double> &q_numDiff_u = args.array<double>("q_numDiff_u");
1252 xt::pyarray<double> &q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
1253 int offset_u = args.scalar<int>("offset_u");
1254 int stride_u = args.scalar<int>("stride_u");
1255 xt::pyarray<double> &globalResidual = args.array<double>("globalResidual");
1256 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
1257 xt::pyarray<int> &exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
1258 xt::pyarray<int> &elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
1259 xt::pyarray<int> &elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
1260 xt::pyarray<double> &ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
1261 xt::pyarray<int> &isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
1262 xt::pyarray<double> &ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
1263 xt::pyarray<int> &isFluxBoundary_u = args.array<int>("isFluxBoundary_u");
1264 xt::pyarray<double> &ebqe_bc_flux_ext = args.array<double>("ebqe_bc_flux_ext");
1265 xt::pyarray<double> &ebqe_phi = args.array<double>("ebqe_phi");
1266 double epsFact = args.scalar<double>("epsFact");
1267 xt::pyarray<double> &ebqe_u = args.array<double>("ebqe_u");
1268 xt::pyarray<double> &ebqe_theta = args.array<double>("ebqe_theta");
1269 xt::pyarray<double> &ebqe_flux = args.array<double>("ebqe_flux");
1270 // PARAMETERS FOR EDGE BASED STABILIZATION
1271 double cE = args.scalar<double>("cE");
1272 double cK = args.scalar<double>("cK");
1273 // PARAMETERS FOR LOG BASED ENTROPY FUNCTION
1274 double uL = args.scalar<double>("uL");
1275 double uR = args.scalar<double>("uR");
1276 // PARAMETERS FOR EDGE VISCOSITY
1277 int numDOFs = args.scalar<int>("numDOFs");
1278 int NNZ = args.scalar<int>("NNZ");
1279 xt::pyarray<int> &csrRowIndeces_DofLoops = args.array<int>("csrRowIndeces_DofLoops");
1280 xt::pyarray<int> &csrColumnOffsets_DofLoops = args.array<int>("csrColumnOffsets_DofLoops");
1281 xt::pyarray<int> &csrRowIndeces_CellLoops = args.array<int>("csrRowIndeces_CellLoops");
1282 xt::pyarray<int> &csrColumnOffsets_CellLoops = args.array<int>("csrColumnOffsets_CellLoops");
1283 xt::pyarray<int> &csrColumnOffsets_eb_CellLoops = args.array<int>("csrColumnOffsets_eb_CellLoops");
1284 // C matrices
1285 xt::pyarray<double> &Cx = args.array<double>("Cx");
1286 xt::pyarray<double> &Cy = args.array<double>("Cy");
1287 xt::pyarray<double> &Cz = args.array<double>("Cz");
1288 xt::pyarray<double> &CTx = args.array<double>("CTx");
1289 xt::pyarray<double> &CTy = args.array<double>("CTy");
1290 xt::pyarray<double> &CTz = args.array<double>("CTz");
1291 xt::pyarray<double> &ML = args.array<double>("ML");
1292 xt::pyarray<double> &MC = args.array<double>("MC");
1293
1294 xt::pyarray<double> &delta_x_ij = args.array<double>("delta_x_ij");
1295 // PARAMETERS FOR 1st or 2nd ORDER MPP METHOD
1296 int LUMPED_MASS_MATRIX = args.scalar<int>("LUMPED_MASS_MATRIX");
1297 STABILIZATION STABILIZATION_TYPE{static_cast<STABILIZATION>(args.scalar<int>("STABILIZATION_TYPE"))};
1298
1299 int ENTROPY_TYPE = args.scalar<int>("ENTROPY_TYPE");
1300 // FOR FCT
1301 xt::pyarray<double> &dLow = args.array<double>("dLow");
1302 xt::pyarray<double> &fluxMatrix = args.array<double>("fluxMatrix");
1303 xt::pyarray<double> &mDotLow = args.array<double>("mDotLow");
1304 xt::pyarray<double> &mLow = args.array<double>("mLow");
1305 xt::pyarray<double> &dt_times_fH_minus_fL = args.array<double>("dt_times_fH_minus_fL");
1306 xt::pyarray<double> &min_m_bc = args.array<double>("min_m_bc");
1307 xt::pyarray<double> &max_m_bc = args.array<double>("max_m_bc");
1308 // AUX QUANTITIES OF INTEREST
1309 xt::pyarray<double> &quantDOFs = args.array<double>("quantDOFs");
1310 xt::pyarray<double> &mn = args.array<double>("mn");
1311 xt::pyarray<double> &fluxCorrection = args.array<double>("fluxCorrection");
1312 xt::pyarray<double> &limited_solution = args.array<double>("limited_solution");
1313 xt::pyarray<int> &freeDOFMaterialTypes = args.array<int>("freeDOFMaterialTypes");
1314
1315 xt::pyarray<double> &velocity_couple = args.array<double>("velocity_couple");
1316 xt::pyarray<double> &ebqe_velocity_ext_couple = args.array<double>("ebqe_velocity_ext_couple");
1317 // xt::pyarray<double> &q_x = args.array<double>("q_x");
1318 // xt::pyarray<double> &ebqe_x = args.array<double>("ebqe_x");
1319
1320 xt::pyarray<double> &anb_seepage_flux_n = args.array<double>("anb_seepage_flux_n");
1321 xt::pyarray<double> &q_velocity = args.array<double>("q_velocity");
1322 double &anb_seepage_flux(args.scalar<double>("anb_seepage_flux"));
1323 anb_seepage_flux = 0.0;
1324 xt::pyarray<int> &csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
1325 xt::pyarray<int> &csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
1326 xt::pyarray<int> &csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
1327 // double Rpos[numDOFs], Rneg[numDOFs];
1328 std::vector<double> Rpos(numDOFs, 0.0), Rneg(numDOFs, 0.0);
1329 std::vector<double> TransportMatrix(NNZ, 0.0),
1330 TransportMatrixConsistent(NNZ, 0.0),
1331 TransportMatrixn(NNZ, 0.0),
1332 TransportMatrixConsistentn(NNZ, 0.0);
1333 //double FluxCorrectionMatrix[NNZ];
1334 // NOTE: This function follows a different (but equivalent) implementation of the smoothness based indicator than NCLS.h
1335 // Allocate space for the transport matrices
1336 // This is used for first order KUZMIN'S METHOD
1337 // double TransportMatrix[NNZ], TransportMatrixConsistent[NNZ];
1338 // double TransportMatrixn[NNZ], TransportMatrixConsistentn[NNZ];
1339 std::valarray<double> u_free_dof(numDOFs);
1340 std::valarray<double> u_free_dof_old(numDOFs);
1341 std::valarray<double> ML2(numDOFs);
1342 // Lumped L2 projection buffers for density
1343 std::vector<double> rho_dof(numDOFs, 0.0);
1344 std::vector<double> ML_rho(numDOFs, 0.0);
1345 std::fill(velocity_couple.data(), velocity_couple.data() + velocity_couple.size(), 0.0);
1346 std::fill(ebqe_velocity_ext_couple.data(), ebqe_velocity_ext_couple.data() + ebqe_velocity_ext_couple.size(), 0.0);
1347
1348 for (int eN = 0; eN < nElements_global; eN++)
1349 for (int j = 0; j < nDOF_trial_element; j++) {
1350 int eN_nDOF_trial_element = eN * nDOF_trial_element;
1351 u_free_dof[r_l2g.data()[eN_nDOF_trial_element + j]] = u_dof.data()[u_l2g.data()[eN_nDOF_trial_element + j]];
1352 u_free_dof_old[r_l2g.data()[eN_nDOF_trial_element + j]] = u_dof_old.data()[u_l2g.data()[eN_nDOF_trial_element + j]];
1353 }
1354 for (int i = 0; i < NNZ; i++) {
1355 TransportMatrix[i] = 0.;
1356 TransportMatrixConsistent[i] = 0.;
1357 TransportMatrixn[i] = 0.;
1358 TransportMatrixConsistentn[i] = 0.;
1359 }
1360
1361 // Project quadrature density to nodal DOFs before constructing nodal
1362 // stabilization potentials so Phi uses a true nodal density field.
1363 for (int eN = 0; eN < nElements_global; eN++) {
1364 const int eN_nDOF_trial_element = eN * nDOF_trial_element;
1365 for (int k = 0; k < nQuadraturePoints_element; k++) {
1366 const int eN_k = eN * nQuadraturePoints_element + k;
1367 double jac[nSpace * nSpace], jacDet, jacInv[nSpace * nSpace], x, y, z;
1368 ck.calculateMapping_element(eN, k, mesh_dof.data(), mesh_l2g.data(),
1369 mesh_trial_ref.data(), mesh_grad_trial_ref.data(),
1370 jac, jacDet, jacInv, x, y, z);
1371 const double dV = fabs(jacDet) * dV_ref.data()[k];
1372 for (int i = 0; i < nDOF_test_element; i++) {
1373 const int eN_i = eN * nDOF_test_element + i;
1374 const int free_gi = r_l2g.data()[eN_i];
1375 const double u_test_dV = u_test_ref.data()[k * nDOF_trial_element + i] * dV;
1376 rho_dof[free_gi] += q_rho.data()[eN_k] * u_test_dV;
1377 ML_rho[free_gi] += u_test_dV;
1378 }
1379 }
1380 }
1381 for (int i = 0; i < numDOFs; ++i) {
1382 if (ML_rho[i] > 0.0) rho_dof[i] /= ML_rho[i];
1383 else rho_dof[i] = rho;
1384 }
1385 // Cache the projected, salinity-coupled density for use in invert() so
1386 // the m -> u inversion is consistent with the forward residual.
1387 rho_dof_member = rho_dof;
1388
1389 // compute entropy and init global_entropy_residual and boundary_integral
1390 double psi[numDOFs], eta[numDOFs], global_entropy_residual[numDOFs], boundary_integral[numDOFs];
1391 for (int i = 0; i < numDOFs; i++) {
1392 // NODAL ENTROPY //
1393 if (STABILIZATION_TYPE == STABILIZATION::EV_Stab) //EV stab
1394 {
1395 double solni = 1.0 * u_free_dof_old[i];
1396 eta[i] = ENTROPY_TYPE == 1 ? ENTROPY(solni, uL, uR) : ENTROPY_LOG(solni, uL, uR);
1397 global_entropy_residual[i] = 0.;
1398 }
1399 boundary_integral[i] = 0.;
1400 ML2[i] = 0.0;
1401 }
1402
1404 // ** LOOP IN CELLS FOR CELL BASED TERMS ** //
1406 // HERE WE COMPUTE:
1407 // * Time derivative term. u_t
1408 // * cell based CFL (for reference)
1409 // * Entropy residual
1410 // * Transport matrices
1411
1412 for (int eN = 0; eN < nElements_global; eN++) {
1413 const int eN_nDOF_trial_element = eN * nDOF_trial_element;
1414 const int eN_nDOF_mesh_trial_element = eN * nDOF_mesh_trial_element;
1415 //declare local storage for local contributions and initialize
1416 double elementResidual_u[nDOF_test_element], element_entropy_residual[nDOF_test_element], Phi[nDOF_trial_element], Phi_n[nDOF_trial_element];
1417 double elementTransport[nDOF_test_element][nDOF_trial_element], elementTransportConsistent[nDOF_test_element][nDOF_trial_element];
1418 double elementTransportn[nDOF_test_element][nDOF_trial_element], elementTransportConsistentn[nDOF_test_element][nDOF_trial_element];
1419 for (int j = 0; j < nDOF_trial_element; j++) {
1420 const int u_gj = u_l2g.data()[eN_nDOF_trial_element + j];
1421 const int free_gj = r_l2g.data()[eN_nDOF_trial_element + j];
1422 const int x_gj = mesh_l2g.data()[eN_nDOF_mesh_trial_element + j];
1423 const double rho_node_j = rho_dof[free_gj];
1424 Phi[j] = u_dof.data()[u_gj];
1425 Phi_n[j] = u_dof_old.data()[u_gj];
1426 for (int I = 0; I < nSpace; I++) {
1427 // Match the main variable-density operator: grad(u) - (rho/rho0) g.
1428 Phi[j] -= (rho_node_j / rho) * mesh_dof.data()[x_gj * 3 + I] * gravity[I];
1429 Phi_n[j] -= (rho_node_j / rho) * mesh_dof.data()[x_gj * 3 + I] * gravity[I];
1430 }
1431 }
1432 for (int i = 0; i < nDOF_test_element; i++) {
1433 elementResidual_u[i] = 0.0;
1434 element_entropy_residual[i] = 0.0;
1435 for (int j = 0; j < nDOF_trial_element; j++) {
1436 elementTransport[i][j] = 0.0;
1437 elementTransportConsistent[i][j] = 0.0;
1438 elementTransportn[i][j] = 0.0;
1439 elementTransportConsistentn[i][j] = 0.0;
1440 }
1441 }
1442 //loop over quadrature points and compute integrands
1443 for (int k = 0; k < nQuadraturePoints_element; k++) {
1444 //compute indeces and declare local storage
1445 int eN_k = eN * nQuadraturePoints_element + k, eN_k_nSpace = eN_k * nSpace;
1446 double
1447 // for entropy residual
1448 aux_entropy_residual = 0.,
1449 DENTROPY_un, DENTROPY_uni,
1450 //for mass matrix contributions
1451 u = 0.0, un = 0.0, grad_phi[nSpace], grad_phi_n[nSpace], grad_u_velocity[nSpace], velocity_loc[nSpace], u_test_dV[nDOF_trial_element], u_grad_trial[nDOF_trial_element * nSpace], u_grad_test_dV[nDOF_test_element * nSpace],
1452 //for general use
1453 jac[nSpace * nSpace], jacDet, jacInv[nSpace * nSpace], dV, x, y, z, xt, yt, zt, m, dm, f[nSpace], df[nSpace], a[nnz], da[nnz], as[nnz], mn, dmn, fn[nSpace], dfn[nSpace], an[nnz], dan[nnz], asn[nnz];
1454 //get the physical integration weight
1455 ck.calculateMapping_element(eN, k, mesh_dof.data(), mesh_l2g.data(), mesh_trial_ref.data(), mesh_grad_trial_ref.data(), jac, jacDet, jacInv, x, y, z);
1456 ck.calculateMappingVelocity_element(eN, k, mesh_velocity_dof.data(), mesh_l2g.data(), mesh_trial_ref.data(), xt, yt, zt);
1457 dV = fabs(jacDet) * dV_ref.data()[k];
1458 //get the solution (of Newton's solver). To compute time derivative term
1459 ck.valFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element], &u_trial_ref.data()[k * nDOF_trial_element], u);
1460 //get the solution at quad point at tn and tnm1 for entropy viscosity
1461 ck.valFromDOF(u_dof_old.data(), &u_l2g.data()[eN_nDOF_trial_element], &u_trial_ref.data()[k * nDOF_trial_element], un);
1462 //get the solution gradients at tn for entropy viscosity
1463 ck.gradTrialFromRef(&u_grad_trial_ref.data()[k * nDOF_trial_element * nSpace], jacInv, u_grad_trial);
1464 ck.gradFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element], u_grad_trial, grad_u_velocity);
1465
1466 //populate q_x
1467 // const int eN_k_3d = eN_k * 3;
1468 // q_x.data()[eN_k_3d + 0] = x;
1469 // q_x.data()[eN_k_3d + 1] = y;
1470 // q_x.data()[eN_k_3d + 2] = z;
1471 //precalculate test function products with integration weights for mass matrix terms
1472 for (int I = 0; I < nSpace; I++) {
1473 grad_phi[I] = 0.0;
1474 grad_phi_n[I] = 0.0;
1475 }
1476 for (int j = 0; j < nDOF_trial_element; j++) {
1477 u_test_dV[j] = u_test_ref.data()[k * nDOF_trial_element + j] * dV;
1478 for (int I = 0; I < nSpace; I++) {
1479 grad_phi_n[I] += Phi_n[j] * u_grad_trial[j * nSpace + I];
1480 grad_phi[I] += Phi[j] * u_grad_trial[j * nSpace + I];
1481 u_grad_test_dV[j * nSpace + I] = u_grad_trial[j * nSpace + I] * dV; //cek warning won't work for Petrov-Galerkin
1482 }
1483 }
1484 //
1485 //calculate pde coefficients at quadrature points
1486 //
1487 double Kr, dKr, Krn, dKrn, thetaW, thetaWn;
1488 const double rho_local = q_rho.data()[eN_k];
1489 const double rho_velocity = std::fabs(rho_local) > 1.0e-12 ? rho_local : rho;
1490
1491 evaluateCoefficients(a_rowptr.data(), a_colind.data(), rho, rho_local, beta, gravity.data(), alpha.data()[elementMaterialTypes[eN]], n.data()[elementMaterialTypes[eN]], thetaR.data()[elementMaterialTypes[eN]], thetaSR.data()[elementMaterialTypes[eN]],
1492 &KWs.data()[elementMaterialTypes[eN] * nnz], un, mn, dmn, fn, dfn, an, dan, asn, Krn, dKrn, thetaWn);
1493 evaluateCoefficients(a_rowptr.data(), a_colind.data(), rho, rho_local, beta, gravity.data(), alpha.data()[elementMaterialTypes[eN]], n.data()[elementMaterialTypes[eN]], thetaR.data()[elementMaterialTypes[eN]], thetaSR.data()[elementMaterialTypes[eN]],
1494 &KWs.data()[elementMaterialTypes[eN] * nnz], u, m, dm, f, df, a, da, as, Kr, dKr, thetaW);
1495 q_theta.data()[eN_k] = thetaW;
1496
1497 // Darcy velocity for coupling should use the direct FE gradient of the
1498 // pressure head. The Phi-based gradients are only for stabilization.
1499 for (int I = 0; I < nSpace; ++I) {
1500 q_velocity.data()[eN_k_nSpace + I] = grad_u_velocity[I];
1501 }
1502
1503 double pressure_gradient[nSpace];
1504 const double rho_ratio = rho_velocity / rho;
1505 for (int J = 0; J < nSpace; ++J)
1506 pressure_gradient[J] = grad_u_velocity[J] - rho_ratio * gravity.data()[J];
1507
1508 for (int I = 0; I < nSpace; ++I) {
1509 double acc = 0.0;
1510 for (int ii = a_rowptr.data()[I]; ii < a_rowptr.data()[I+1]; ++ii) {
1511 const int J = a_colind.data()[ii];
1512 acc += (a[ii] / rho_velocity) * pressure_gradient[J];
1513 }
1514 velocity.data()[eN_k_nSpace + I] = -acc;
1515 velocity_couple.data()[eN_k_nSpace + I] = -acc;
1516 }
1517 //
1518 //moving mesh
1519 //
1520 double mesh_velocity[3];
1521 mesh_velocity[0] = xt;
1522 mesh_velocity[1] = yt;
1523 mesh_velocity[2] = zt;
1524 //relative velocity at tn
1525 for (int I = 0; I < nSpace; I++) {
1526 f[I] -= MOVING_DOMAIN * m * mesh_velocity[I];
1527 velocity_loc[I] = df[I] * (2.0 * dm * dm / (dm * dm + fmax(1.0e-16, dm * dm)));
1528 }
1530 // CALCULATE CELL BASED CFL //
1532 calculateCFL(elementDiameter.data()[eN] / degree_polynomial, velocity_loc, cfl.data()[eN_k]);
1534 // CALCULATE ENTROPY RESIDUAL AT QUAD POINT //
1536 if (STABILIZATION_TYPE == STABILIZATION::EV_Stab) // EV stab
1537 {
1538 for (int I = 0; I < nSpace; I++) aux_entropy_residual += velocity_loc[I] * grad_phi_n[I];
1539 DENTROPY_un = ENTROPY_TYPE == 1 ? DENTROPY(un, uL, uR) : DENTROPY_LOG(un, uL, uR);
1540 }
1542 // ith-LOOP //
1544 for (int i = 0; i < nDOF_test_element; i++) {
1545 // VECTOR OF ENTROPY RESIDUAL //
1546 int eN_i = eN * nDOF_test_element + i;
1547 ML2[u_l2g.data()[eN_i]] += u_test_dV[i];
1548 if (STABILIZATION_TYPE == STABILIZATION::EV_Stab) // EV stab
1549 {
1550 int gi = offset_u + stride_u * u_l2g.data()[eN_i]; //global i-th index
1551 double uni = u_dof_old.data()[gi];
1552 DENTROPY_uni = ENTROPY_TYPE == 1 ? DENTROPY(uni, uL, uR) : DENTROPY_LOG(uni, uL, uR);
1553 element_entropy_residual[i] += (DENTROPY_un - DENTROPY_uni) * aux_entropy_residual * u_test_dV[i];
1554 }
1555
1556 elementResidual_u[i] += m * u_test_dV[i];
1558 // j-th LOOP // To construct transport matrices
1560
1561 for (int j = 0; j < nDOF_trial_element; j++) {
1562 int j_nSpace = j * nSpace;
1563 int i_nSpace = i * nSpace;
1564 elementTransport[i][j] += ck.SimpleDiffusionJacobian_weak(a_rowptr.data(), a_colind.data(), as, &u_grad_trial[j_nSpace], &u_grad_test_dV[i_nSpace]);
1565 elementTransportConsistent[i][j] += ck.SimpleDiffusionJacobian_weak(a_rowptr.data(), a_colind.data(), a, &u_grad_trial[j_nSpace], &u_grad_test_dV[i_nSpace]);
1566 elementTransportn[i][j] += ck.SimpleDiffusionJacobian_weak(a_rowptr.data(), a_colind.data(), asn, &u_grad_trial[j_nSpace], &u_grad_test_dV[i_nSpace]);
1567 elementTransportConsistentn[i][j] += ck.SimpleDiffusionJacobian_weak(a_rowptr.data(), a_colind.data(), an, &u_grad_trial[j_nSpace], &u_grad_test_dV[i_nSpace]);
1568 }
1569 } //i
1570 //save solution for other models
1571 q_u.data()[eN_k] = u;
1572 q_m.data()[eN_k] = m;
1573 }
1575 // DISTRIBUTE // load cell based element into global residual
1577 for (int i = 0; i < nDOF_test_element; i++) {
1578 int eN_i = eN * nDOF_test_element + i;
1579 int gi = offset_u + stride_u * r_l2g.data()[eN_i]; //global i-th index
1580 // distribute entropy_residual
1581 if (STABILIZATION_TYPE == STABILIZATION::EV_Stab) // EV Stab
1582 global_entropy_residual[gi] += element_entropy_residual[i];
1583 // distribute transport matrices
1584 for (int j = 0; j < nDOF_trial_element; j++) {
1585 int eN_i_j = eN_i * nDOF_trial_element + j;
1586 TransportMatrix[csrRowIndeces_CellLoops.data()[eN_i] + csrColumnOffsets_CellLoops.data()[eN_i_j]] += elementTransport[i][j];
1587 TransportMatrixConsistent[csrRowIndeces_CellLoops.data()[eN_i] + csrColumnOffsets_CellLoops.data()[eN_i_j]] += elementTransportConsistent[i][j];
1588 TransportMatrixn[csrRowIndeces_CellLoops.data()[eN_i] + csrColumnOffsets_CellLoops.data()[eN_i_j]] += elementTransportn[i][j];
1589 TransportMatrixConsistentn[csrRowIndeces_CellLoops.data()[eN_i] + csrColumnOffsets_CellLoops.data()[eN_i_j]] += elementTransportConsistentn[i][j];
1590 } //j
1591 } //i
1592
1593 }
1594
1595 //loop over exterior element boundaries to calculate surface integrals and load into element and global residuals
1596 //
1597 //ebNE is the Exterior element boundary INdex
1598 //ebN is the element boundary INdex
1599 //eN is the element index
1600 for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++) {
1601 int ebN = exteriorElementBoundariesArray.data()[ebNE], eN = elementBoundaryElementsArray.data()[ebN * 2 + 0], ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN * 2 + 0], eN_nDOF_trial_element = eN * nDOF_trial_element;
1602 double elementResidual_u[nDOF_test_element];
1603 for (int i = 0; i < nDOF_test_element; i++) { elementResidual_u[i] = 0.0; }
1604 for (int kb = 0; kb < nQuadraturePoints_elementBoundary; kb++) {
1605 int ebNE_kb = ebNE * nQuadraturePoints_elementBoundary + kb, ebNE_kb_nSpace = ebNE_kb * nSpace, ebN_local_kb = ebN_local * nQuadraturePoints_elementBoundary + kb, ebN_local_kb_nSpace = ebN_local_kb * nSpace;
1606 double u_ext = 0.0, un_ext, grad_u_ext[nSpace], m_ext = 0.0, dm_ext = 0.0, f_ext[nSpace], df_ext[nSpace], a_ext[nnz], da_ext[nnz], as_ext[nnz],
1607 mn_ext = 0.0, dmn_ext = 0.0, fn_ext[nSpace], dfn_ext[nSpace], an_ext[nnz], dan_ext[nnz], asn_ext[nnz], flux_ext = 0.0, bflux_ext = 0.0,
1608 //anb_seepage_flux=0.0, // for flux calculation
1609 bc_u_ext = 0.0, bc_grad_u_ext[nSpace], bc_m_ext = 0.0, bc_dm_ext = 0.0, bc_f_ext[nSpace], bc_df_ext[nSpace], bc_a_ext[nnz], bc_da_ext[nnz], bc_as_ext[nnz], jac_ext[nSpace * nSpace], jacDet_ext, jacInv_ext[nSpace * nSpace], boundaryJac[nSpace * (nSpace - 1)], metricTensor[(nSpace - 1) * (nSpace - 1)], metricTensorDetSqrt, dS, u_test_dS[nDOF_test_element], u_grad_trial_trace[nDOF_trial_element * nSpace], normal[3], x_ext, y_ext, z_ext, xt_ext, yt_ext, zt_ext, integralScaling, G[nSpace * nSpace], G_dd_G, tr_G, fluxJacobian_u_u[nDOF_trial_element], bfluxJacobian_u_u[nDOF_trial_element], fluxJacobian_un_un[nDOF_trial_element];
1610 for (int j = 0; j < nDOF_trial_element; j++) {
1611 fluxJacobian_u_u[j] = 0.0; bfluxJacobian_u_u[j] = 0.0; fluxJacobian_un_un[j] = 0.0;
1612 }
1613 //
1614 //calculate the solution and gradients at quadrature points
1615 //
1616 //compute information about mapping from reference element to physical element
1617 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,
1618 normal_ref.data(), normal, x_ext, y_ext, z_ext);
1619 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);
1620 dS = ((1.0 - MOVING_DOMAIN) * metricTensorDetSqrt + MOVING_DOMAIN * integralScaling) * dS_ref.data()[kb];
1621 //get the metric tensor
1622 //cek todo use symmetry
1623 ck.calculateG(jacInv_ext, G, G_dd_G, tr_G);
1624 //compute shape and solution information
1625 //shape
1626 ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace * nDOF_trial_element], jacInv_ext, u_grad_trial_trace);
1627 //solution and gradient
1628 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);
1629 ck.valFromDOF(u_dof_old.data(), &u_l2g.data()[eN_nDOF_trial_element], &u_trial_trace_ref.data()[ebN_local_kb * nDOF_test_element], un_ext);
1630 ck.gradFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element], u_grad_trial_trace, grad_u_ext);
1631
1632 //populate ebqe_x
1633 // const int ebNE_kb_3d = ebNE_kb * 3;
1634 // ebqe_x.data()[ebNE_kb_3d + 0] = x_ext;
1635 // ebqe_x.data()[ebNE_kb_3d + 1] = y_ext;
1636 // ebqe_x.data()[ebNE_kb_3d + 2] = z_ext;
1637
1638
1639 //precalculate test function products with integration weights
1640 for (int j = 0; j < nDOF_trial_element; j++) { u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb * nDOF_test_element + j] * dS; }
1641 //
1642 //load the boundary values
1643 //
1644 bc_u_ext = isDOFBoundary_u.data()[ebNE_kb] * ebqe_bc_u_ext.data()[ebNE_kb] + (1 - isDOFBoundary_u.data()[ebNE_kb]) * u_ext;
1645 //
1646 //calculate the pde coefficients using the solution and the boundary values for the solution
1647 //
1648 double bc_Kr, bc_dKr,bc_Kr_ext, bc_dKr_ext, bc_Krn, bc_dKrn, thetaW_ext, thetaWn_ext, thetaW_bc_ext;
1649 const double rho_ext = ebqe_rho.data()[ebNE_kb];
1650 const double rho_velocity_ext = std::fabs(rho_ext) > 1.0e-12 ? rho_ext : rho;
1651
1652 evaluateCoefficients(a_rowptr.data(), a_colind.data(), rho, rho_ext, beta, gravity.data(), alpha.data()[elementMaterialTypes.data()[eN]], n.data()[elementMaterialTypes.data()[eN]], thetaR.data()[elementMaterialTypes.data()[eN]],
1653 thetaSR.data()[elementMaterialTypes.data()[eN]], &KWs.data()[elementMaterialTypes.data()[eN] * nnz], u_ext, m_ext, dm_ext, f_ext, df_ext, a_ext, da_ext, as_ext, bc_Kr, bc_dKr, thetaW_ext);
1654 evaluateCoefficients(a_rowptr.data(), a_colind.data(), rho, rho_ext, beta, gravity.data(), alpha.data()[elementMaterialTypes.data()[eN]], n.data()[elementMaterialTypes.data()[eN]], thetaR.data()[elementMaterialTypes.data()[eN]],
1655 thetaSR.data()[elementMaterialTypes.data()[eN]], &KWs.data()[elementMaterialTypes.data()[eN] * nnz], un_ext, mn_ext, dmn_ext, fn_ext, dfn_ext, an_ext, dan_ext, asn_ext, bc_Krn, bc_dKrn, thetaWn_ext);
1656 evaluateCoefficients(a_rowptr.data(), a_colind.data(), rho, rho_ext, beta, gravity.data(), alpha.data()[elementMaterialTypes.data()[eN]], n.data()[elementMaterialTypes.data()[eN]], thetaR.data()[elementMaterialTypes.data()[eN]],
1657 thetaSR.data()[elementMaterialTypes.data()[eN]], &KWs.data()[elementMaterialTypes.data()[eN] * nnz], bc_u_ext, bc_m_ext, bc_dm_ext, bc_f_ext, bc_df_ext, bc_a_ext, bc_da_ext, bc_as_ext, bc_Kr_ext,bc_dKr_ext, thetaW_bc_ext);
1658 ebqe_theta.data()[ebNE_kb] = thetaW_ext;
1659 //
1660 //Calculate Darcy Velocity at external faces
1661 //
1662
1663 // double darcy_velocity_loc_ext[nSpace];
1664 // for (int I = 0; I < nSpace; I++) { darcy_velocity_loc_ext[I] = 0.0; }
1665
1666 // for (int I = 0; I < nSpace; I++) {
1667 // for (int J = 0; J < nSpace; J++) { darcy_velocity_loc_ext[I] -= bc_Kr * KWs.data()[elementMaterialTypes[eN] * nSpace * nSpace + I * nSpace + J] * (grad_u_ext[J]+ gravity.data()[J]); }
1668 // }
1669 // for (int I = 0; I < nSpace; I++) { ebqe_velocity_ext_couple.data()[ebNE_kb_nSpace + I] = darcy_velocity_loc_ext[I] ; }
1670
1671 double ext_pressure_gradient[nSpace];
1672 const double rho_ratio_ext = rho_velocity_ext / rho;
1673 for (int J = 0; J < nSpace; ++J)
1674 ext_pressure_gradient[J] = grad_u_ext[J] - rho_ratio_ext * gravity.data()[J];
1675
1676 for (int I = 0; I < nSpace; ++I) {
1677 double acc = 0.0;
1678 for (int ii = a_rowptr.data()[I]; ii < a_rowptr.data()[I+1]; ++ii) {
1679 const int J = a_colind.data()[ii];
1680 acc += (a_ext[ii] / rho_velocity_ext) * ext_pressure_gradient[J];
1681 }
1682 ebqe_velocity_ext.data()[ebNE_kb_nSpace + I] = -acc;
1683 ebqe_velocity_ext_couple.data()[ebNE_kb_nSpace + I] = -acc;
1684 }
1685
1686
1687 //
1688 //calculate the numerical fluxes
1689 //
1690 bool useConsistentFlux=false;
1691 if (useConsistentFlux) {
1692 exteriorNumericalFlux(ebqe_bc_flux_ext[ebNE_kb], a_rowptr.data(), a_colind.data(),
1693 isSeepageFace.data()[ebNE], //tricky, this is a face flag not face quad
1694 isDOFBoundary_u.data()[ebNE_kb], normal, bc_u_ext, a_ext, grad_u_ext, u_ext, f_ext,
1695 ebqe_penalty_ext.data()[ebNE_kb], // penalty,
1696 flux_ext);
1697 } else {
1698 exteriorNumericalFlux2(ebqe_bc_flux_ext[ebNE_kb], a_rowptr.data(), a_colind.data(),
1699 isSeepageFace.data()[ebNE], //tricky, this is a face flag not face quad
1700 isDOFBoundary_u.data()[ebNE_kb], normal, bc_u_ext, a_ext, grad_u_ext, u_ext, f_ext,
1701 ebqe_penalty_ext.data()[ebNE_kb], // penalty,
1702 flux_ext, bflux_ext);
1703 }
1704
1705 ebqe_flux.data()[ebNE_kb] = flux_ext;
1706
1707 anb_seepage_flux = seepagefluxcalculator(anb_seepage_flux, isSeepageFace.data()[ebNE], dS, flux_ext);
1708 anb_seepage_flux_n.data()[0] = anb_seepage_flux;
1709 ebqe_u.data()[ebNE_kb] = u_ext;
1710 //seed the FCT bounds with the imposed value. min_m_bc/max_m_bc are otherwise left at +-1e10, so
1711 //Zalesak's bounds at a weakly imposed Dirichlet DOF come only from the neighbours' low order masses
1712 //and the limiter clamps the boundary layer. theta(bc_u) uses the DOF's own soil and density so the
1713 //bound lies on the same retention curve as mLow in the edge loop below. done after the flux call so
1714 //the seepage active set is current
1715 if (isDOFBoundary_u.data()[ebNE_kb]) {
1716 double bc_u_dof = isSeepageFace.data()[ebNE] ? 0.0 : ebqe_bc_u_ext.data()[ebNE_kb];
1717 for (int i = 0; i < nDOF_test_element; i++) {
1718 if (u_test_trace_ref.data()[ebN_local_kb * nDOF_test_element + i] <= 1.0e-12) continue; //DOF not on this face
1719 int free_gi = r_l2g.data()[eN * nDOF_test_element + i], mat_gi = freeDOFMaterialTypes.data()[free_gi];
1720 double m_bc, dm_bc, f_bc[nSpace], df_bc[nSpace], a_bc[nnz], da_bc[nnz], as_bc[nnz], Kr_bc, dKr_bc, thetaW_bc;
1721 evaluateCoefficients(a_rowptr.data(), a_colind.data(), rho, rho_dof[free_gi], beta, gravity.data(), alpha.data()[mat_gi],
1722 n.data()[mat_gi], thetaR.data()[mat_gi], thetaSR.data()[mat_gi], &KWs.data()[mat_gi * nnz], bc_u_dof,
1723 m_bc, dm_bc, f_bc, df_bc, a_bc, da_bc, as_bc, Kr_bc, dKr_bc, thetaW_bc);
1724 min_m_bc.data()[free_gi] = fmin(min_m_bc.data()[free_gi], m_bc);
1725 max_m_bc.data()[free_gi] = fmax(max_m_bc.data()[free_gi], m_bc);
1726 }
1727 }
1728 //
1729 //update residuals
1730 //
1731 for (int i = 0; i < nDOF_test_element; i++) {
1732 if (useConsistentFlux) {
1733 elementResidual_u[i] += ck.ExteriorElementBoundaryFlux(flux_ext, u_test_dS[i]);
1734 } else {
1735 elementResidual_u[i] += ck.ExteriorElementBoundaryFlux(bflux_ext, u_test_dS[i]);
1736 }
1737 } //i
1738 for (int j = 0; j < nDOF_trial_element; j++) {
1739 if (useConsistentFlux) {
1740 exteriorNumericalFluxJacobian(a_rowptr.data(), a_colind.data(), isDOFBoundary_u.data()[ebNE_kb], normal, a_ext, da_ext, grad_u_ext, &u_grad_trial_trace[j * nSpace], df_ext, u_trial_trace_ref.data()[ebN_local_kb * nDOF_test_element + j],
1741 ebqe_penalty_ext.data()[ebNE_kb], //penalty,
1742 fluxJacobian_u_u[j]);
1743 } else {
1744 exteriorNumericalFluxJacobian2(a_rowptr.data(), a_colind.data(), isDOFBoundary_u.data()[ebNE_kb], normal, as_ext, a_ext, da_ext, grad_u_ext, &u_grad_trial_trace[j * nSpace], df_ext, u_trial_trace_ref.data()[ebN_local_kb * nDOF_test_element + j],
1745 ebqe_penalty_ext.data()[ebNE_kb], //penalty,
1746 fluxJacobian_u_u[j],bfluxJacobian_u_u[j]);
1747 }
1748 //probably need isDOFBoundary_un here
1749 //exteriorNumericalFluxJacobian(a_rowptr.data(), a_colind.data(), isDOFBoundary_u.data()[ebNE_kb], normal, asn_ext, dan_ext, grad_u_ext, &u_grad_trial_trace[j * nSpace], dfn_ext, u_trial_trace_ref.data()[ebN_local_kb * nDOF_test_element + j],
1750 // ebqe_penalty_ext.data()[ebNE_kb], //penalty,
1751 // fluxJacobian_un_un[j]);
1752 } //j
1753 //
1754 //update the element and global residual storage
1755 //
1756 for (int i = 0; i < nDOF_test_element; i++) {
1757 int eN_i = eN * nDOF_test_element + i;
1758 for (int j = 0; j < nDOF_trial_element; j++) {
1759 int ebN_i_j = ebN * 4 * nDOF_test_X_trial_element + i * nDOF_trial_element + j;
1760 if (useConsistentFlux) {
1761 globalJacobian.data()[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] += fluxJacobian_u_u[j] * u_test_dS[i];
1762 } else {
1763 globalJacobian.data()[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] += bfluxJacobian_u_u[j] * u_test_dS[i];
1764 TransportMatrix[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] += fluxJacobian_u_u[j] * u_test_dS[i];
1765 TransportMatrixConsistent[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] += fluxJacobian_u_u[j] * u_test_dS[i];
1766 TransportMatrixn[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] += fluxJacobian_un_un[j] * u_test_dS[i];
1767 TransportMatrixConsistentn[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] += fluxJacobian_un_un[j] * u_test_dS[i];
1768 }
1769 } //j
1770 } //i
1771 } //kb
1772 for (int i = 0; i < nDOF_test_element; i++) {
1773 int eN_i = eN * nDOF_test_element + i;
1774 globalResidual.data()[offset_u + stride_u * u_l2g.data()[eN_i]] += elementResidual_u[i];
1775 }//i
1776 } //ebNE
1778 // COMPUTE SMOOTHNESS INDICATOR and NORMALIZE ENTROPY RESIDUAL //
1780 // NOTE: see NCLS.h for a different but equivalent implementation of this.
1781 int ij = 0;
1782 double cflux[numDOFs];
1783 for (int i = 0; i < numDOFs; i++) {
1784 double gi[nSpace], Cij[nSpace], xi[nSpace], etaMaxi, etaMini;
1785 double solni = u_free_dof_old[i];
1786 for (int I = 0; I < nSpace; I++) {
1787 solni -= (rho_dof[i] / rho) * gravity.data()[I] * mesh_dof.data()[i * 3 + I];
1788 }
1789 if (STABILIZATION_TYPE == STABILIZATION::EV_Stab) //EV Stabilization
1790 {
1791 // For eta min and max
1792 etaMaxi = fabs(eta[i]);
1793 etaMini = fabs(eta[i]);
1794 }
1795 // initialize gi and compute xi
1796 for (int I = 0; I < nSpace; I++) {
1797 gi[I] = 0.;
1798 xi[I] = mesh_dof.data()[i * 3 + I];
1799 }
1800 // for smoothness indicator //
1801 double alpha_numerator_pos = 0., alpha_numerator_neg = 0., alpha_denominator_pos = 0., alpha_denominator_neg = 0.;
1802 for (int offset = csrRowIndeces_DofLoops.data()[i]; offset < csrRowIndeces_DofLoops.data()[i + 1]; offset++) { // First loop in j (sparsity pattern)
1803 int j = csrColumnOffsets_DofLoops.data()[offset];
1804 if (STABILIZATION_TYPE == STABILIZATION::EV_Stab) //EV Stabilization
1805 {
1806 // COMPUTE ETA MIN AND ETA MAX //
1807 etaMaxi = fmax(etaMaxi, fabs(eta[j]));
1808 etaMini = fmin(etaMini, fabs(eta[j]));
1809 }
1810 double solnj = u_free_dof_old[j];
1811 for (int I = 0; I < nSpace; I++) {
1812 solnj -= (rho_dof[j] / rho) * gravity.data()[I] * mesh_dof.data()[j * 3 + I];
1813 }
1814 // Update Cij matrices
1815 Cij[0] = Cx[ij];
1816#if nSpace == 2
1817 Cij[1] = Cy[ij];
1818#endif
1819#if nSpace == 3
1820 Cij[2] = Cz[ij];
1821#endif
1822 // COMPUTE gi VECTOR. gi=1/mi*sum_j(Cij*solj)
1823 for (int I = 0; I < nSpace; I++) gi[I] += Cij[I] * solnj;
1824
1825 // COMPUTE numerator and denominator of smoothness indicator
1826 double alpha_num = solni - solnj;
1827 if (alpha_num >= 0.) {
1828 alpha_numerator_pos += alpha_num;
1829 alpha_denominator_pos += alpha_num;
1830 } else {
1831 alpha_numerator_neg += alpha_num;
1832 alpha_denominator_neg += fabs(alpha_num);
1833 }
1834 //update ij
1835 ij += 1;
1836 }
1837 // scale g vector by lumped mass matrix
1838 //double mass_matrix_error = abs(ML.data()[i] - ML2[i]);
1839 //if (mass_matrix_error > 1.0e-16) std::cout << mass_matrix_error<<" ML " << ML.data()[i] << '\t' << ML2[i] << std::endl;
1840 for (int I = 0; I < nSpace; I++) gi[I] /= ML.data()[i];
1841 if (STABILIZATION_TYPE == STABILIZATION::EV_Stab) //EV Stab
1842 {
1843 // Normalizae entropy residual
1844 global_entropy_residual[i] *= etaMini == etaMaxi ? 0. : 2 * cE / (etaMaxi - etaMini);
1845 quantDOFs.data()[i] = fabs(global_entropy_residual[i]);
1846 }
1847
1848 // Now that I have the gi vectors, I can use them for the current i-th DOF
1849 double SumPos = 0., SumNeg = 0.;
1850 for (int offset = csrRowIndeces_DofLoops.data()[i]; offset < csrRowIndeces_DofLoops.data()[i + 1]; offset++) { // second loop in j (sparsity pattern)
1851 int j = csrColumnOffsets_DofLoops.data()[offset];
1852 // compute xj
1853 double xj[nSpace];
1854 for (int I = 0; I < nSpace; I++) xj[I] = mesh_dof.data()[j * 3 + I];
1855 // compute gi*(xi-xj)
1856 double gi_times_x = 0.;
1857 for (int I = 0; I < nSpace; I++) {
1858 gi_times_x += gi[I] * delta_x_ij.data()[offset * 3 + I];
1859 }
1860 // compute the positive and negative part of gi*(xi-xj)
1861 SumPos += gi_times_x > 0 ? gi_times_x : 0;
1862 SumNeg += gi_times_x < 0 ? gi_times_x : 0;
1863 }
1864 double sigmaPosi = fmin(1., (fabs(SumNeg) + 1E-15) / (SumPos + 1E-15));
1865 double sigmaNegi = fmin(1., (SumPos + 1E-15) / (fabs(SumNeg) + 1E-15));
1866 double alpha_numi = fabs(sigmaPosi * alpha_numerator_pos + sigmaNegi * alpha_numerator_neg);
1867 double alpha_deni = sigmaPosi * alpha_denominator_pos + sigmaNegi * alpha_denominator_neg;
1868 if (IS_BETAij_ONE == 1) {
1869 alpha_numi = fabs(alpha_numerator_pos + alpha_numerator_neg);
1870 alpha_deni = alpha_denominator_pos + alpha_denominator_neg;
1871 }
1872 double alphai = alpha_numi / (alpha_deni + 1E-15);
1873 quantDOFs.data()[i] = alphai;
1874
1875 if (POWER_SMOOTHNESS_INDICATOR == 0) psi[i] = 1.0;
1876 else psi[i] = std::pow(alphai, POWER_SMOOTHNESS_INDICATOR); //NOTE: they use alpha^2 in the paper
1877 }
1879 // ** LOOP IN DOFs FOR EDGE BASED TERMS ** //
1881 ij = 0;
1882 for (int i = 0; i < numDOFs; i++) {
1883 int ii;
1884 double sum_abs_dt_times_fH_minus_fL = 0.0, MLi = ML.data()[i];
1885 double Kr, dKr, Krn, dKrn;
1886 double J_ii = 0.0;
1887 double ith_dissipative_term = 0;
1888 double ith_low_order_dissipative_term = 0;
1889 double ith_flux_term = 0;
1890 double ith_consistent_flux_term = 0;
1891 double dLii = 0.;
1892 double m, dm, f[nSpace], df[nSpace], a[nnz], da[nnz], as[nnz];
1893 double dmn, fn[nSpace], dfn[nSpace], an[nnz], dan[nnz], asn[nnz];
1894
1895 const double rho_i = rho_dof[i];
1896 // Rock region of each DOF. The nodal m, mn and upwind k_rw below are
1897 // evaluated at a DOF, so they take that DOF's soil; a single
1898 // elementMaterialTypes[0] here would put the whole mesh on one retention
1899 // curve and, since invert() maps m back with the per-DOF soil, would make
1900 // psi -> m -> psi something other than the identity.
1901 const int mat_i = freeDOFMaterialTypes.data()[i];
1902
1903 double thetaW_tmp = 0.0;
1904 // loop over the sparsity pattern of the i-th DOF
1905 for (int offset = csrRowIndeces_DofLoops.data()[i]; offset < csrRowIndeces_DofLoops.data()[i + 1]; offset++) {
1906 int j = csrColumnOffsets_DofLoops.data()[offset];
1907 if (i == j) ii = ij;
1908 const double rho_j = rho_dof[j];
1909 const int mat_j = freeDOFMaterialTypes.data()[j];
1910 const double rho_edge = 0.5 * (rho_i + rho_j);
1911 double delta_phi = u_free_dof[j] - u_free_dof[i];
1912 double delta_phin = u_free_dof_old[j] - u_free_dof_old[i];
1913 // Match evaluateCoefficients(): use an edge-based hydrostatic jump
1914 // with local density scaling rho/rho0, rather than grad(rho g·x),
1915 // which would introduce a spurious (g·x) grad(rho) term.
1916 for (int I = 0; I < nSpace; I++) {
1917 const double delta_x = mesh_dof.data()[j * 3 + I] - mesh_dof.data()[i * 3 + I];
1918 const double hydrostatic_jump = (rho_edge / rho) * gravity.data()[I] * delta_x;
1919 delta_phi -= hydrostatic_jump;
1920 delta_phin -= hydrostatic_jump;
1921 }
1922 double dLowij, dLij, dEVij, dHij, fH, fL, fA=0.0;
1923 double fL_CN =0.0, fA_CN=0.0 ;
1924 fH = -Theta * TransportMatrixConsistent[ij] * delta_phi - (1 - Theta) * TransportMatrixConsistentn[ij] * delta_phin;
1925 // fH = -Theta_h * TransportMatrixConsistent[ij] * delta_phi - (1 - Theta_h) * TransportMatrixConsistentn[ij] * delta_phin; //previous: Theta
1926 ith_consistent_flux_term += fH;
1927 fA = fH;
1928 fA_CN = fH;
1929
1930 if (-TransportMatrix[ij] * delta_phi <= 0.0) {
1931 evaluateCoefficients(a_rowptr.data(), a_colind.data(), rho, rho_i, beta, gravity.data(),
1932 alpha.data()[mat_i],
1933 n.data()[mat_i], thetaR.data()[mat_i], thetaSR.data()[mat_i], &KWs.data()[mat_i * nnz], u_free_dof[i], m, dm, f, df, a, da, as, Kr, dKr, thetaW_tmp);
1934 fL = Theta * Kr * fmax(0.0, -TransportMatrix[ij]) * delta_phi;
1935 fL_CN = Theta_h * Kr * fmax(0.0, -TransportMatrix[ij]) * delta_phi;
1936
1937 if (i != j) {
1938 globalJacobian.data()[ij] -= Theta * Kr * fmax(0.0, -TransportMatrix[ij]);
1939 J_ii -= -Theta * Kr * fmax(0.0, -TransportMatrix[ij]) + Theta * dKr * fmax(0.0, -TransportMatrix[ij]) * delta_phi;
1940 }
1941 ith_flux_term += fL;
1942 fA -= fL;
1943 } else {
1944 evaluateCoefficients(a_rowptr.data(), a_colind.data(), rho, rho_j, beta, gravity.data(),
1945 alpha.data()[mat_j],
1946 n.data()[mat_j], thetaR.data()[mat_j], thetaSR.data()[mat_j], &KWs.data()[mat_j * nnz], u_free_dof[j], m, dm, f, df, a, da, as, Kr, dKr, thetaW_tmp);
1947 fL = Theta * Kr * fmax(0.0, -TransportMatrix[ij]) * delta_phi;
1948 fL_CN = Theta_h * Kr * fmax(0.0, -TransportMatrix[ij]) * delta_phi;
1949
1950 if (i != j) {
1951 globalJacobian.data()[ij] -= Theta * Kr * fmax(0.0, -TransportMatrix[ij]) + Theta * dKr * fmax(0.0, -TransportMatrix[ij]) * delta_phi;
1952 J_ii -= -Theta * Kr * fmax(0.0, -TransportMatrix[ij]);
1953 }
1954 ith_flux_term += fL;
1955 fA -= fL;
1956 }
1957 if (-TransportMatrixn[ij] * delta_phin <= 0.0) {
1958 evaluateCoefficients(a_rowptr.data(), a_colind.data(), rho, rho_i, beta, gravity.data(),
1959 alpha.data()[mat_i],
1960 n.data()[mat_i], thetaR.data()[mat_i], thetaSR.data()[mat_i], &KWs.data()[mat_i * nnz], u_free_dof_old[i], m, dm, f, df, a, da, as, Kr, dKr, thetaW_tmp);
1961 fL = (1 - Theta) * Kr * fmax(0.0, -TransportMatrixn[ij]) * delta_phin;
1962 fL_CN += (1 - Theta_h) * Kr * fmax(0.0, -TransportMatrixn[ij]) * delta_phin;
1963 ith_flux_term += fL;
1964 fA -= fL;
1965 fA_CN -= fL_CN;
1966 } else {
1967 evaluateCoefficients(a_rowptr.data(), a_colind.data(), rho, rho_j, beta, gravity.data(),
1968 alpha.data()[mat_j],
1969 n.data()[mat_j], thetaR.data()[mat_j], thetaSR.data()[mat_j], &KWs.data()[mat_j * nnz], u_free_dof_old[j], m, dm, f, df, a, da, as, Kr, dKr, thetaW_tmp);
1970 fL = (1 - Theta) * Kr * fmax(0.0, -TransportMatrixn[ij]) * delta_phin;
1971 fL_CN += (1 - Theta_h) * Kr * fmax(0.0, -TransportMatrixn[ij]) * delta_phin;
1972 ith_flux_term += fL;
1973 fA -= fL;
1974 fA_CN -= fL_CN;
1975 }
1976 dt_times_fH_minus_fL.data()[ij] = dt * fA;
1977 //dt_times_fH_minus_fL.data()[ij] = dt * fA_CN;
1978 ij += 1;
1979 }
1980 mDotLow.data()[i] = ith_flux_term/MLi;
1981 cflux[i] = ith_consistent_flux_term;
1982 evaluateCoefficients(a_rowptr.data(), a_colind.data(), rho, rho_i, beta, gravity.data(),
1983 alpha.data()[mat_i],
1984 n.data()[mat_i], thetaR.data()[mat_i], thetaSR.data()[mat_i], &KWs.data()[mat_i * nnz], u_free_dof[i], m, dm, f, df, a, da, as, Kr, dKr, thetaW_tmp);
1985 evaluateCoefficients(a_rowptr.data(), a_colind.data(), rho, rho_i, beta, gravity.data(),
1986 alpha.data()[mat_i],
1987 n.data()[mat_i], thetaR.data()[mat_i], thetaSR.data()[mat_i], &KWs.data()[mat_i * nnz], u_free_dof_old[i], mn.data()[i], dmn, fn, dfn, an, dan, asn, Krn, dKrn, thetaW_tmp);
1988 mLow.data()[i] = m;
1989 globalResidual.data()[i] += bc_mask.data()[i] * (MLi * (m - mn.data()[i]) / dt - ith_flux_term);
1990 globalJacobian.data()[ii] += bc_mask.data()[i] * (MLi * dm / dt + J_ii) + (1.0 - bc_mask.data()[i]);
1991 }
1992 if (STABILIZATION_TYPE == STABILIZATION::Implicit_FCT) {
1993 //FCTStep(args);
1994 for (int i = 0; i < numDOFs; i++) {
1995 globalResidual.data()[i] += fluxCorrection.data()[i];
1996 }
1997 }
1998
1999 }
2000
2001 void invert(arguments_dict &args)
2002 {
2003 xt::pyarray<int> &a_rowptr = args.array<int>("a_rowptr");
2004 xt::pyarray<int> &a_colind = args.array<int>("a_colind");
2005 double rho = args.scalar<double>("rho"); // freshwater reference (fallback)
2006 double beta = args.scalar<double>("beta");
2007 xt::pyarray<double> &gravity = args.array<double>("gravity");
2008 xt::pyarray<double> &alpha = args.array<double>("alpha");
2009 xt::pyarray<double> &n = args.array<double>("n");
2010 xt::pyarray<double> &thetaR = args.array<double>("thetaR");
2011 xt::pyarray<double> &thetaSR = args.array<double>("thetaSR");
2012 PSK_TYPE_member = args.scalar<int>("PSK_TYPE");
2013 xt::pyarray<double> &KWs = args.array<double>("KWs");
2014 xt::pyarray<int> &elementMaterialTypes = args.array<int>("elementMaterialTypes");
2015 xt::pyarray<int> &freeDOFMaterialTypes = args.array<int>("freeDOFMaterialTypes");
2016 int numDOFs = args.scalar<int>("numDOFs");
2017 xt::pyarray<double> &mIn = args.array<double>("limited_solution");
2018 xt::pyarray<double> &pOut = args.array<double>("u_dof");
2019 // Per-DOF density projected from q_rho in
2020 // calculateResidual_entropy_viscosity (cached as a class member there).
2021 // If the cache is empty (e.g. invert called before the first residual
2022 // evaluation, or coupling disabled) fall back to the scalar `rho`.
2023 const bool have_rho_dof = (rho_dof_member.size() == static_cast<std::size_t>(numDOFs));
2024 int USE_NEWTON_INVERT = args.scalar<int>("USE_NEWTON_INVERT");
2025
2026 for (int i = 0; i < numDOFs; i++) {
2027 const int material_i = freeDOFMaterialTypes.data()[i];
2028 // Use the salinity-coupled density at this DOF so the m -> u inversion
2029 // is consistent with the variable-density forward residual
2030 // m_forward = theta * rho_local * (...).
2031 const double rho_i = have_rho_dof ? rho_dof_member[i] : rho;
2032 double dm, f[nSpace], df[nSpace], a[nnz], da[nnz];
2033 if (USE_NEWTON_INVERT){
2034 evaluateInverseCoefficients_Newton(a_rowptr.data(), a_colind.data(), rho_i, beta, gravity.data(), alpha.data()[material_i], n.data()[material_i], thetaR.data()[material_i],
2035 thetaSR.data()[material_i], &KWs.data()[material_i * nnz],
2036 pOut.data()[i], mIn.data()[i],
2037 dm, f, df, a, da);
2038 }
2039 else{
2040 evaluateInverseCoefficients(a_rowptr.data(), a_colind.data(), rho_i, beta, gravity.data(), alpha.data()[material_i], n.data()[material_i], thetaR.data()[material_i],
2041 thetaSR.data()[material_i], &KWs.data()[material_i * nnz],
2042 pOut.data()[i], mIn.data()[i],
2043 dm, f, df, a, da);
2044 }
2045 }
2046 }
2047
2048 void calculateMassMatrix(arguments_dict &args)
2049 {
2050 //element
2051 double dt = args.scalar<double>("dt");
2052 xt::pyarray<double> &mesh_trial_ref = args.array<double>("mesh_trial_ref");
2053 xt::pyarray<double> &mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
2054 xt::pyarray<double> &mesh_dof = args.array<double>("mesh_dof");
2055 xt::pyarray<double> &mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
2056 double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
2057 xt::pyarray<int> &mesh_l2g = args.array<int>("mesh_l2g");
2058 xt::pyarray<double> &dV_ref = args.array<double>("dV_ref");
2059 xt::pyarray<double> &u_trial_ref = args.array<double>("u_trial_ref");
2060 xt::pyarray<double> &u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
2061 xt::pyarray<double> &u_test_ref = args.array<double>("u_test_ref");
2062 xt::pyarray<double> &u_grad_test_ref = args.array<double>("u_grad_test_ref");
2063 //element boundary
2064 xt::pyarray<double> &mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
2065 xt::pyarray<double> &mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
2066 xt::pyarray<double> &dS_ref = args.array<double>("dS_ref");
2067 xt::pyarray<double> &u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
2068 xt::pyarray<double> &u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
2069 xt::pyarray<double> &u_test_trace_ref = args.array<double>("u_test_trace_ref");
2070 xt::pyarray<double> &u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
2071 xt::pyarray<double> &normal_ref = args.array<double>("normal_ref");
2072 xt::pyarray<double> &boundaryJac_ref = args.array<double>("boundaryJac_ref");
2073 //physics
2074 int nElements_global = args.scalar<int>("nElements_global");
2075 //new
2076 xt::pyarray<double> &ebqe_penalty_ext = args.array<double>("ebqe_penalty_ext");
2077 xt::pyarray<int> &elementMaterialTypes = args.array<int>("elementMaterialTypes");
2078 xt::pyarray<int> &isSeepageFace = args.array<int>("isSeepageFace");
2079 xt::pyarray<int> &a_rowptr = args.array<int>("a_rowptr");
2080 xt::pyarray<int> &a_colind = args.array<int>("a_colind");
2081 double rho = args.scalar<double>("rho");
2082 double beta = args.scalar<double>("beta");
2083
2084 xt::pyarray<double> &q_rho = args.array<double>("q_rho");
2085
2086 xt::pyarray<double> &gravity = args.array<double>("gravity");
2087 xt::pyarray<double> &alpha = args.array<double>("alpha");
2088 xt::pyarray<double> &n = args.array<double>("n");
2089 xt::pyarray<double> &thetaR = args.array<double>("thetaR");
2090 xt::pyarray<double> &thetaSR = args.array<double>("thetaSR");
2091 PSK_TYPE_member = args.scalar<int>("PSK_TYPE");
2092 xt::pyarray<double> &KWs = args.array<double>("KWs");
2093 //end new
2094 double useMetrics = args.scalar<double>("useMetrics");
2095 double alphaBDF = args.scalar<double>("alphaBDF");
2096 int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
2097 double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
2098 xt::pyarray<int> &u_l2g = args.array<int>("u_l2g");
2099 xt::pyarray<int> &r_l2g = args.array<int>("r_l2g");
2100 xt::pyarray<double> &elementDiameter = args.array<double>("elementDiameter");
2101 int degree_polynomial = args.scalar<int>("degree_polynomial");
2102 xt::pyarray<double> &u_dof = args.array<double>("u_dof");
2103 xt::pyarray<double> &velocity = args.array<double>("velocity");
2104 xt::pyarray<double> &q_m_betaBDF = args.array<double>("q_m_betaBDF");
2105 xt::pyarray<double> &cfl = args.array<double>("cfl");
2106 xt::pyarray<double> &q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
2107 xt::pyarray<int> &csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
2108 xt::pyarray<int> &csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
2109 xt::pyarray<double> &globalJacobian = args.array<double>("globalJacobian");
2110 xt::pyarray<double> &delta_x_ij = args.array<double>("delta_x_ij");
2111 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
2112 xt::pyarray<int> &exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
2113 xt::pyarray<int> &elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
2114 xt::pyarray<int> &elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
2115 xt::pyarray<double> &ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
2116 xt::pyarray<int> &isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
2117 xt::pyarray<double> &ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
2118 xt::pyarray<int> &isFluxBoundary_u = args.array<int>("isFluxBoundary_u");
2119 xt::pyarray<double> &ebqe_bc_flux_ext = args.array<double>("ebqe_bc_flux_ext");
2120 xt::pyarray<int> &csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
2121 int LUMPED_MASS_MATRIX = args.scalar<int>("LUMPED_MASS_MATRIX");
2122 double Ct_sge = 4.0;
2123 //
2124 //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
2125 //
2126 for (int eN = 0; eN < nElements_global; eN++) {
2127 double elementJacobian_u_u[nDOF_test_element][nDOF_trial_element];
2128 for (int i = 0; i < nDOF_test_element; i++)
2129 for (int j = 0; j < nDOF_trial_element; j++) { elementJacobian_u_u[i][j] = 0.0; }
2130 for (int k = 0; k < nQuadraturePoints_element; k++) {
2131 int eN_k = eN * nQuadraturePoints_element + k, //index to a scalar at a quadrature point
2132 eN_k_nSpace = eN_k * nSpace,
2133 eN_nDOF_trial_element = eN * nDOF_trial_element; //index to a vector at a quadrature point
2134 //declare local storage
2135 double u = 0.0, grad_u[nSpace], m = 0.0, dm = 0.0, f[nSpace], df[nSpace], a[nnz], da[nnz], as[nnz], m_t = 0.0, dm_t = 0.0, dpdeResidual_u_u[nDOF_trial_element], Lstar_u[nDOF_test_element], dsubgridError_u_u[nDOF_trial_element], tau = 0.0, tau0 = 0.0, tau1 = 0.0, jac[nSpace * nSpace], jacDet, jacInv[nSpace * nSpace], u_grad_trial[nDOF_trial_element * nSpace], dV, u_test_dV[nDOF_test_element], u_grad_test_dV[nDOF_test_element * nSpace], x, y, z, xt, yt, zt,
2136 G[nSpace * nSpace], G_dd_G, tr_G;
2137
2138 //get jacobian, etc for mapping reference element
2139 ck.calculateMapping_element(eN, k, mesh_dof.data(), mesh_l2g.data(), mesh_trial_ref.data(), mesh_grad_trial_ref.data(), jac, jacDet, jacInv, x, y, z);
2140 ck.calculateMappingVelocity_element(eN, k, mesh_velocity_dof.data(), mesh_l2g.data(), mesh_trial_ref.data(), xt, yt, zt);
2141 //get the physical integration weight
2142 dV = fabs(jacDet) * dV_ref.data()[k];
2143 ck.calculateG(jacInv, G, G_dd_G, tr_G);
2144 //get the trial function gradients
2145 ck.gradTrialFromRef(&u_grad_trial_ref.data()[k * nDOF_trial_element * nSpace], jacInv, u_grad_trial);
2146 //get the solution
2147 ck.valFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element], &u_trial_ref.data()[k * nDOF_trial_element], u);
2148 //get the solution gradients
2149 ck.gradFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element], u_grad_trial, grad_u);
2150 //precalculate test function products with integration weights
2151 for (int j = 0; j < nDOF_trial_element; j++) {
2152 u_test_dV[j] = u_test_ref.data()[k * nDOF_trial_element + j] * dV;
2153 for (int I = 0; I < nSpace; I++) {
2154 u_grad_test_dV[j * nSpace + I] = u_grad_trial[j * nSpace + I] * dV; //cek warning won't work for Petrov-Galerkin
2155 }
2156 }
2157 //
2158 //calculate pde coefficients and derivatives at quadrature points
2159 //
2160 double Kr, dKr, thetaW;
2161 //const double rho_local = q_rho.data()[eN_k];
2162 evaluateCoefficients(a_rowptr.data(), a_colind.data(), rho, q_rho.data()[eN_k], beta, gravity.data(), alpha.data()[elementMaterialTypes.data()[eN]], n.data()[elementMaterialTypes.data()[eN]], thetaR.data()[elementMaterialTypes.data()[eN]],
2163 thetaSR.data()[elementMaterialTypes.data()[eN]], &KWs.data()[elementMaterialTypes.data()[eN] * nnz], u, m, dm, f, df, a, da, as, Kr, dKr, thetaW);
2164 //
2165 //moving mesh
2166 //
2167 double mesh_velocity[3];
2168 mesh_velocity[0] = xt;
2169 mesh_velocity[1] = yt;
2170 mesh_velocity[2] = zt;
2171 for (int I = 0; I < nSpace; I++) {
2172 f[I] -= MOVING_DOMAIN * m * mesh_velocity[I];
2173 df[I] -= MOVING_DOMAIN * dm * mesh_velocity[I];
2174 }
2175 //
2176 //calculate time derivatives
2177 //
2178 //cek hack
2179 dm = 1.0;
2180 ck.bdf(alphaBDF,
2181 q_m_betaBDF.data()[eN_k], //since m_t isn't used, we don't have to correct mass
2182 m, dm, m_t, dm_t);
2183 //
2184 //calculate subgrid error contribution to the Jacobian (strong residual, adjoint, jacobian of strong residual)
2185 //
2186 //calculate the adjoint times the test functions
2187 for (int i = 0; i < nDOF_test_element; i++) {
2188 int i_nSpace = i * nSpace;
2189 Lstar_u[i] = ck.Advection_adjoint(df, &u_grad_test_dV[i_nSpace]);
2190 }
2191 //calculate the Jacobian of strong residual
2192 for (int j = 0; j < nDOF_trial_element; j++) {
2193 int j_nSpace = j * nSpace;
2194 dpdeResidual_u_u[j] = ck.MassJacobian_strong(dm_t, u_trial_ref.data()[k * nDOF_trial_element + j]) + ck.AdvectionJacobian_strong(df, &u_grad_trial[j_nSpace]);
2195 }
2196 //tau and tau*Res
2197 calculateSubgridError_tau(elementDiameter.data()[eN], dm_t, df, cfl.data()[eN_k], tau0);
2198
2199 calculateSubgridError_tau(Ct_sge, G, dm_t, df, tau1, cfl.data()[eN_k]);
2200 tau = useMetrics * tau1 + (1.0 - useMetrics) * tau0;
2201
2202 for (int j = 0; j < nDOF_trial_element; j++) dsubgridError_u_u[j] = -tau * dpdeResidual_u_u[j];
2203 for (int i = 0; i < nDOF_test_element; i++) {
2204 for (int j = 0; j < nDOF_trial_element; j++) {
2205 if (LUMPED_MASS_MATRIX == 1) {
2206 if (i == j) elementJacobian_u_u[i][j] += u_test_dV[i];
2207 } else {
2208 int j_nSpace = j * nSpace;
2209 int i_nSpace = i * nSpace;
2210 dm_t = 1.0; //we are solving for continuum density explicitly
2211 elementJacobian_u_u[i][j] += ck.MassJacobian_weak(dm_t, u_trial_ref.data()[k * nDOF_trial_element + j], u_test_dV[i]);
2212 }
2213 } //j
2214 } //i
2215 } //k
2216 //
2217 //load into element Jacobian into global Jacobian
2218 //
2219 for (int i = 0; i < nDOF_test_element; i++) {
2220 int eN_i = eN * nDOF_test_element + i;
2221 int I = u_l2g.data()[eN_i];
2222 for (int j = 0; j < nDOF_trial_element; j++) {
2223 int eN_i_j = eN_i * nDOF_trial_element + j;
2224 int J = u_l2g.data()[eN * nDOF_trial_element + j];
2225 //globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_u_u.data()[eN_i_j]] += elementJacobian_u_u[i][j];
2226 delta_x_ij.data()[3 * (csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_u_u.data()[eN_i_j]) + 0] = mesh_dof.data()[I * 3 + 0] - mesh_dof.data()[J * 3 + 0];
2227 delta_x_ij.data()[3 * (csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_u_u.data()[eN_i_j]) + 1] = mesh_dof.data()[I * 3 + 1] - mesh_dof.data()[J * 3 + 1];
2228 delta_x_ij.data()[3 * (csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_u_u.data()[eN_i_j]) + 2] = mesh_dof.data()[I * 3 + 2] - mesh_dof.data()[J * 3 + 2];
2229 } //j
2230 } //i
2231 } //elements
2232 } //computeMassMatrix
2233}; //Richards
2234
2235inline Richards_base *newRichards(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
2236{
2237 if (nSpaceIn == 1)
2238 return proteus::chooseAndAllocateDiscretization1D<Richards_base, Richards, CompKernel>(nSpaceIn, nQuadraturePoints_elementIn, nDOF_mesh_trial_elementIn, nDOF_trial_elementIn, nDOF_test_elementIn, nQuadraturePoints_elementBoundaryIn, CompKernelFlag);
2239 else if (nSpaceIn == 2)
2240 return proteus::chooseAndAllocateDiscretization2D<Richards_base, Richards, CompKernel>(nSpaceIn, nQuadraturePoints_elementIn, nDOF_mesh_trial_elementIn, nDOF_trial_elementIn, nDOF_test_elementIn, nQuadraturePoints_elementBoundaryIn, CompKernelFlag);
2241 else {
2242 assert(nSpaceIn == 3);
2243 return proteus::chooseAndAllocateDiscretization<Richards_base, Richards, CompKernel>(nSpaceIn, nQuadraturePoints_elementIn, nDOF_mesh_trial_elementIn, nDOF_trial_elementIn, nDOF_test_elementIn, nQuadraturePoints_elementBoundaryIn, CompKernelFlag);
2244 }
2245}
2246} // namespace richards
2247} // namespace proteus
2248#endif
n
Int n
Definition Headers.h:28
u
Double u
Definition Headers.h:89
num
Int num
Definition Headers.h:32
z
Double * z
Definition Headers.h:49
v
Double v
Definition Headers.h:95
psi
Double psi
Definition Headers.h:78
cE
#define cE
Definition NCLS3P.h:10
proteus::richards::Richards_base
Definition Richards.h:53
proteus::richards::Richards_base::FCTStep
virtual void FCTStep(arguments_dict &args)=0
proteus::richards::Richards_base::kth_FCT_step
virtual void kth_FCT_step(arguments_dict &args)=0
proteus::richards::Richards_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)=0
proteus::richards::Richards_base::calculateResidual_entropy_viscosity
virtual void calculateResidual_entropy_viscosity(arguments_dict &args)=0
proteus::richards::Richards_base::invert
virtual void invert(arguments_dict &args)=0
proteus::richards::Richards_base::~Richards_base
virtual ~Richards_base()
Definition Richards.h:56
proteus::richards::Richards_base::calculateMassMatrix
virtual void calculateMassMatrix(arguments_dict &args)=0
proteus::richards::Richards_base::calculateResidual
virtual void calculateResidual(arguments_dict &args)=0
proteus::richards::Richards::exteriorNumericalFlux
void exteriorNumericalFlux(const double &bc_flux, int rowptr[nSpace], int colind[nnz], int isSeepageFace, int &isDOFBoundary, double n[nSpace], double bc_u, double K[nnz], double grad_psi[nSpace], double u, double K_rho_g[nSpace], double penalty, double &flux)
Definition Richards.h:216
proteus::richards::Richards::calculateNumericalDiffusion
void calculateNumericalDiffusion(const double &shockCapturingDiffusion, const double &elementDiameter, const double &strong_residual, const double grad_u[nSpace], double &numDiff)
Definition Richards.h:205
proteus::richards::Richards::calculateResidual_entropy_viscosity
void calculateResidual_entropy_viscosity(arguments_dict &args)
Definition Richards.h:1186
proteus::richards::Richards::seepagefluxcalculator
double seepagefluxcalculator(double anb_seepage_flux, int isSeepageFace, double dS, double flux_ext)
Definition Richards.h:307
proteus::richards::Richards::invert
void invert(arguments_dict &args)
Definition Richards.h:2001
proteus::richards::Richards::kth_FCT_step
void kth_FCT_step(arguments_dict &args)
Definition Richards.h:1051
proteus::richards::Richards::nDOF_test_X_trial_element
const int nDOF_test_X_trial_element
Definition Richards.h:69
proteus::richards::Richards::evaluateInverseCoefficients
void evaluateInverseCoefficients(const int rowptr[nSpace], const int colind[nnz], const double rho, const double beta, const double gravity[nSpace], const double alpha, const double n_vg, const double thetaR, const double thetaSR, const double KWs[nnz], double &u, const double &m, const double &dm, const double f[nSpace], const double df[nSpace], const double a[nnz], const double da[nnz])
Definition Richards.h:126
proteus::richards::Richards::evaluateInverseCoefficients_Newton
void evaluateInverseCoefficients_Newton(const int rowptr[nSpace], const int colind[nnz], const double rho, const double beta, const double gravity[nSpace], const double alpha, const double n_vg, const double thetaR, const double thetaSR, const double KWs[nnz], double &u, const double &m, const double &dm, const double f[nSpace], const double df[nSpace], const double a[nnz], const double da[nnz])
Definition Richards.h:142
proteus::richards::Richards::evaluateCoefficients
void evaluateCoefficients(const int rowptr[nSpace], const int colind[nnz], const double rho0, const double rho_transport, const double beta, const double gravity[nSpace], const double alpha, const double n_vg, const double thetaR, const double thetaSR, const double KWs[nnz], const double &u, double &m, double &dm, double f[nSpace], double df[nSpace], double a[nnz], double da[nnz], double as[nnz], double &kr, double &dkr, double &thetaW_out)
Definition Richards.h:81
proteus::richards::Richards::exteriorNumericalFlux2
void exteriorNumericalFlux2(const double &bc_flux, int rowptr[nSpace], int colind[nnz], int isSeepageFace, int &isDOFBoundary, double n[nSpace], double bc_u, double K[nnz], double grad_psi[nSpace], double u, double K_rho_g[nSpace], double penalty, double &flux, double &bflux)
Definition Richards.h:258
proteus::richards::Richards::FCTStep
void FCTStep(arguments_dict &args)
Definition Richards.h:901
proteus::richards::Richards::calculateMassMatrix
void calculateMassMatrix(arguments_dict &args)
Definition Richards.h:2048
proteus::richards::Richards::calculateSubgridError_tau
void calculateSubgridError_tau(const double &elementDiameter, const double &dmt, const double dH[nSpace], double &cfl, double &tau)
Definition Richards.h:185
proteus::richards::Richards::exteriorNumericalFluxJacobian2
void exteriorNumericalFluxJacobian2(const int rowptr[nSpace], const int colind[nnz], const int isDOFBoundary, const double n[nSpace], const double Ks[nnz], const double K[nnz], const double dK[nnz], const double grad_psi[nSpace], const double grad_v[nSpace], const double dK_rho_g[nSpace], const double v, const double penalty, double &fluxJacobian, double &bfluxJacobian)
Definition Richards.h:289
proteus::richards::Richards::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 Richards.h:197
proteus::richards::Richards::calculateCFL
void calculateCFL(const double &elementDiameter, const double df[nSpace], double &cfl)
Definition Richards.h:175
proteus::richards::Richards::rho_dof_member
std::vector< double > rho_dof_member
Definition Richards.h:74
proteus::richards::Richards::exteriorNumericalFluxJacobian
void exteriorNumericalFluxJacobian(const int rowptr[nSpace], const int colind[nnz], const int isDOFBoundary, const double n[nSpace], const double K[nnz], const double dK[nnz], const double grad_psi[nSpace], const double grad_v[nSpace], const double dK_rho_g[nSpace], const double v, const double penalty, double &fluxJacobian)
Definition Richards.h:243
proteus::richards::Richards::calculateResidual
void calculateResidual(arguments_dict &args)
Definition Richards.h:313
proteus::richards::Richards::calculateJacobian
void calculateJacobian(arguments_dict &args)
Definition Richards.h:675
df
double df(double C, double b, double a, int q, int r)
Definition analyticalSolutions.c:2209
POWER_SMOOTHNESS_INDICATOR
#define POWER_SMOOTHNESS_INDICATOR
Definition m_comp_co2.h:22
nnz
#define nnz
Definition m_comp_co2.h:19
GLOBAL_FCT
#define GLOBAL_FCT
Definition m_comp_co2.h:24
IS_BETAij_ONE
#define IS_BETAij_ONE
Definition m_comp_co2.h:23
proteus::richards::psk::vgm_invert_analytic
void vgm_invert_analytic(const double m, const double rho, const double alpha, const double n_vg, const double thetaR, const double thetaSR, double &u)
Definition pskRelations.h:918
proteus::richards::psk::gardner_invert_analytic
void gardner_invert_analytic(const double m, const double rho, const double alpha, const double n_vg, const double thetaR, const double thetaSR, double &u)
Definition pskRelations.h:1376
proteus::richards::psk::bc_wetting
void bc_wetting(const double psiC, const double alpha, const double lam, const double thetaR, const double thetaSR, double &thetaW, double &DthetaW_DpsiC, double &KWr, double &DKWr_DpsiC, const bc_kr kr_model=bc_kr::burdine)
Definition pskRelations.h:1096
proteus::richards::psk::vgm_wetting
void vgm_wetting(const double psiC, const double alpha, const double n_vg, const double thetaR, const double thetaSR, double &thetaW, double &DthetaW_DpsiC, double &KWr, double &DKWr_DpsiC)
Definition pskRelations.h:858
proteus::richards::psk::bc_invert_newton
void bc_invert_newton(const double m, const double rho, const double beta, const double alpha, const double lam, const double thetaR, const double thetaSR, double &u)
Definition pskRelations.h:1194
proteus::richards::psk::gardner_invert_newton
void gardner_invert_newton(const double m, const double rho, const double beta, const double alpha, const double n_vg, const double thetaR, const double thetaSR, double &u)
Definition pskRelations.h:1412
proteus::richards::psk::vgm_invert_newton
void vgm_invert_newton(const double m, const double rho, const double beta, const double alpha, const double n_vg, const double thetaR, const double thetaSR, double &u)
Definition pskRelations.h:954
proteus::richards::psk::gardner_wetting
void gardner_wetting(const double psiC, const double alpha, const double n_vg, const double thetaR, const double thetaSR, double &thetaW, double &DthetaW_DpsiC, double &KWr, double &DKWr_DpsiC)
Definition pskRelations.h:1343
proteus::richards::psk::bc_invert_analytic
void bc_invert_analytic(const double m, const double rho, const double alpha, const double lam, const double thetaR, const double thetaSR, double &u)
Definition pskRelations.h:1136
proteus::richards::DENTROPY_LOG
double DENTROPY_LOG(const double &phi, const double &phiL, const double &phiR)
Definition Richards.h:43
proteus::richards::ENTROPY
double ENTROPY(const double &phi, const double &phiL, const double &phiR)
Definition Richards.h:30
proteus::richards::DENTROPY
double DENTROPY(const double &phi, const double &phiL, const double &phiR)
Definition Richards.h:34
proteus::richards::ENTROPY_LOG
double ENTROPY_LOG(const double &phi, const double &phiL, const double &phiR)
Definition Richards.h:39
proteus::richards::newRichards
Richards_base * newRichards(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition Richards.h:2235
proteus
Definition ADR.h:19
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::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::f
double f(const double &g, const double &h, const double &hZ)
Definition SW2DCV.h:58
pskRelations.h
Richards free-function PSK closures (moved here from proteus/richards/psk_models.h).
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