proteus 1.9.0
C/C++/Fortran libraries
Loading...
Searching...
No Matches
richards/ADR.h
Go to the documentation of this file.
1#ifndef ADR_H
2#define ADR_H
3#include <cmath>
4#include <iostream>
5#include <valarray>
6#include "CompKernel.h"
7#include "ModelFactory.h"
8#include "../mprans/ArgumentsDict.h"
9#include "xtensor-python/pyarray.hpp"
10#define nnz nSpace
11
12namespace py = pybind11;
13#define POWER_SMOOTHNESS_INDICATOR 2
14#define IS_BETAij_ONE 0
15#define GLOBAL_FCT 0
16
17namespace proteus
18{
19namespace richards
20{
21 class ADR_base
22 {
23 //The base class defining the interface
24 public:
25 virtual ~ADR_base(){}
26 virtual void calculateResidual(arguments_dict& args)=0;
27 virtual void calculateJacobian(arguments_dict& args)=0;
28 };
29
30 template<class CompKernelType,
31 int nSpace,
32 int nQuadraturePoints_element,
33 int nDOF_mesh_trial_element,
34 int nDOF_trial_element,
35 int nDOF_test_element,
36 int nQuadraturePoints_elementBoundary>
37
38 class ADR : public ADR_base
39 {
40 public:
41 const int nDOF_test_X_trial_element;
42 CompKernelType ck;
43 ADR():
44 nDOF_test_X_trial_element(nDOF_test_element*nDOF_trial_element),
45 ck()
46 {}
47
48 inline
49
50 void evaluateCoefficients(const int rowptr[nSpace],
51 const int colind[nnz],
52 const double rho,
53 const double beta,
54 const double gravity[nSpace],
55 const double alpha,
56 const double n_vg,
57 const double thetaR,
58 const double thetaSR,
59 const double KWs[nnz],
60 const double& u,
61 double& m,
62 double& dm,
63 double f[nSpace],
64 double df[nSpace],
65 double a[nnz],
66 double da[nnz],
67 double as[nnz],
68 double& kr,
69 double& dkr)
70 {
71 const int nSpace2 = nSpace * nSpace;
72 double psiC;
73 double pcBar;
74 double pcBar_n;
75 double pcBar_nM1;
76 double pcBar_nM2;
77 double onePlus_pcBar_n;
78 double sBar;
79 double sqrt_sBar;
80 double DsBar_DpsiC;
81 double thetaW;
82 double DthetaW_DpsiC;
83 double vBar;
84 double vBar2;
85 double DvBar_DpsiC;
86 double KWr;
87 double DKWr_DpsiC;
88 double rho2 = rho * rho;
89 double thetaS;
90 double rhom;
91 double drhom;
92 double m_vg;
93 double pcBarStar;
94 double sqrt_sBarStar;
95 m = u; //rhom*thetaW; //u
96 dm = 1.0; //-rhom*DthetaW_DpsiC+drhom*thetaW;// 1.0;//
97 for (int I=0;I<nSpace;I++)
98 {
99 f[I] = 0.0;
100 df[I] = 0.0;
101
102 for (int ii=rowptr[I]; ii < rowptr[I+1]; ii++)
103 {
104 double velocity = 5.0;
105 double D= 0.02;
106 f[I] = velocity*u; //rho2*KWr*KWs[ii]*gravity[colind[ii]];//velocity *u;//
107 df[I] = velocity; //-rho2*DKWr_DpsiC*KWs[ii]*gravity[colind[ii]];//velocity
108 a[ii] = D ; //rho*KWs[ii]; //rho*KWr*KWs[ii];
109 da[ii] = 0.0; //-rho*DKWr_DpsiC*KWs[ii]; //0.0;//
110 //as[ii] = 1.0; //rho*KWs[ii];//0.0;//rho*KWs[ii];
111 //kr = 1.0;//KWr;// 0.0
112 //dkr=0.0; //mod picard DKWr_DpsiC;
113 }
114 }
115 }
116 void calculateResidual(arguments_dict& args)
117 {
118 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
119 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
120 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
121 xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
122 double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
123 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
124 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
125 xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
126 xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
127 xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
128 xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
129 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
130 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
131 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
132 xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
133 xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
134 xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
135 xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
136 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
137 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
138 int nElements_global = args.scalar<int>("nElements_global");
139 xt::pyarray<double>& ebqe_penalty_ext = args.array<double>("ebqe_penalty_ext");
140 xt::pyarray<int>& elementMaterialTypes = args.array<int>("elementMaterialTypes");
141 xt::pyarray<int>& isSeepageFace = args.array<int>("isSeepageFace");
142 xt::pyarray<int>& a_rowptr = args.array<int>("a_rowptr");
143 xt::pyarray<int>& a_colind = args.array<int>("a_colind");
144 double rho = args.scalar<double>("rho");
145 double beta = args.scalar<double>("beta");
146 xt::pyarray<double>& gravity = args.array<double>("gravity");
147 xt::pyarray<double>& alpha = args.array<double>("alpha");
148 xt::pyarray<double>& n = args.array<double>("n");
149 xt::pyarray<double>& thetaR = args.array<double>("thetaR");
150 xt::pyarray<double>& thetaSR = args.array<double>("thetaSR");
151 xt::pyarray<double>& KWs = args.array<double>("KWs");
152 double useMetrics = args.scalar<double>("useMetrics");
153 double alphaBDF = args.scalar<double>("alphaBDF");
154 int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
155 double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
156 double sc_uref = args.scalar<double>("sc_uref");
157 double sc_alpha = args.scalar<double>("sc_alpha");
158 xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
159 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
160 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
161 xt::pyarray<double>& u_dof_old = args.array<double>("u_dof_old");
162 xt::pyarray<double>& velocity = args.array<double>("velocity");
163 xt::pyarray<double>& q_m = args.array<double>("q_m");
164 xt::pyarray<double>& q_u = args.array<double>("q_u");
165 xt::pyarray<double>& q_dV = args.array<double>("q_dV");
166 xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
167 xt::pyarray<double>& cfl = args.array<double>("cfl");
168 xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
169 xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
170 int offset_u = args.scalar<int>("offset_u");
171 int stride_u = args.scalar<int>("stride_u");
172 xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
173 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
174 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
175 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
176 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
177 xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
178 xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
179 xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
180 xt::pyarray<int>& isFluxBoundary_u = args.array<int>("isFluxBoundary_u");
181 xt::pyarray<double>& ebqe_bc_flux_ext = args.array<double>("ebqe_bc_flux_ext");
182 xt::pyarray<double>& ebqe_phi = args.array<double>("ebqe_phi");
183 double epsFact = args.scalar<double>("epsFact");
184 xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
185 xt::pyarray<double>& ebqe_flux = args.array<double>("ebqe_flux");
186 // PARAMETERS FOR EDGE BASED STABILIZATION
187 double cE = args.scalar<double>("cE");
188 double cK = args.scalar<double>("cK");
189 // PARAMETERS FOR LOG BASED ENTROPY FUNCTION
190 double uL = args.scalar<double>("uL");
191 double uR = args.scalar<double>("uR");
192 // PARAMETERS FOR EDGE VISCOSITY
193 int numDOFs = args.scalar<int>("numDOFs");
194 int NNZ = args.scalar<int>("NNZ");
195 xt::pyarray<int>& csrRowIndeces_DofLoops = args.array<int>("csrRowIndeces_DofLoops");
196 xt::pyarray<int>& csrColumnOffsets_DofLoops = args.array<int>("csrColumnOffsets_DofLoops");
197 xt::pyarray<int>& csrRowIndeces_CellLoops = args.array<int>("csrRowIndeces_CellLoops");
198 xt::pyarray<int>& csrColumnOffsets_CellLoops = args.array<int>("csrColumnOffsets_CellLoops");
199 xt::pyarray<int>& csrColumnOffsets_eb_CellLoops = args.array<int>("csrColumnOffsets_eb_CellLoops");
200 // C matrices
201 xt::pyarray<double>& Cx = args.array<double>("Cx");
202 xt::pyarray<double>& Cy = args.array<double>("Cy");
203 xt::pyarray<double>& Cz = args.array<double>("Cz");
204 xt::pyarray<double>& CTx = args.array<double>("CTx");
205 xt::pyarray<double>& CTy = args.array<double>("CTy");
206 xt::pyarray<double>& CTz = args.array<double>("CTz");
207 xt::pyarray<double>& ML = args.array<double>("ML");
208 xt::pyarray<double>& delta_x_ij = args.array<double>("delta_x_ij");
209 // PARAMETERS FOR 1st or 2nd ORDER MPP METHOD
210 int LUMPED_MASS_MATRIX = args.scalar<int>("LUMPED_MASS_MATRIX");
211 int STABILIZATION_TYPE = args.scalar<int>("STABILIZATION_TYPE");
212 int ENTROPY_TYPE = args.scalar<int>("ENTROPY_TYPE");
213 // FOR FCT
214 xt::pyarray<double>& dLow = args.array<double>("dLow");
215 xt::pyarray<double>& fluxMatrix = args.array<double>("fluxMatrix");
216 xt::pyarray<double>& uDotLow = args.array<double>("uDotLow");
217 xt::pyarray<double>& uLow = args.array<double>("uLow");
218 xt::pyarray<double>& dt_times_fH_minus_fL = args.array<double>("dt_times_fH_minus_fL");
219 xt::pyarray<double>& min_s_bc = args.array<double>("min_s_bc");
220 xt::pyarray<double>& max_s_bc = args.array<double>("max_s_bc");
221 // AUX QUANTITIES OF INTEREST
222 xt::pyarray<double>& quantDOFs = args.array<double>("quantDOFs");
223 xt::pyarray<double>& sLow = args.array<double>("sLow");
224 xt::pyarray<double>& sn = args.array<double>("sn");
225
226 assert(a_rowptr.data()[nSpace] == nnz);
227 assert(a_rowptr.data()[nSpace] == nSpace);
228 //cek should this be read in?
229 double Ct_sge = 4.0;
230 // if (LUMPED_MASS_MATRIX ==1)
231 // {
232 // mass_lumping(nElements_global, nDOF_test_element, nQuadraturePoints_element, u_test_ref, dV_ref, ML, u_l2g);
233 // }
234
235 //loop over elements to compute volume integrals and load them into element and global residual
236 //
237 //eN is the element index
238 //eN_k is the quadrature point index for a scalar
239 //eN_k_nSpace is the quadrature point index for a vector
240 //eN_i is the element test function index
241 //eN_j is the element trial function index
242 //eN_k_j is the quadrature point index for a trial function
243 //eN_k_i is the quadrature point index for a trial function
244 for(int eN=0;eN<nElements_global;eN++)
245 {
246 //declare local storage for element residual and initialize
247 double elementResidual_u[nDOF_test_element];
248 for (int i=0;i<nDOF_test_element;i++)
249 {
250 elementResidual_u[i]=0.0;
251 }//i
252 //loop over quadrature points and compute integrands
253 for (int k=0;k<nQuadraturePoints_element;k++)
254 {
255 //compute indeces and declare local storage
256 int eN_k = eN*nQuadraturePoints_element+k,
257 eN_k_nSpace = eN_k*nSpace,
258 eN_nDOF_trial_element = eN*nDOF_trial_element;
259 double u=0.0,grad_u[nSpace],grad_u_old[nSpace],
260 m=0.0,dm=0.0,
261 f[nSpace],df[nSpace],
262 a[nnz],da[nnz],as[nnz],
263 m_t=0.0,dm_t=0.0,
264 pdeResidual_u=0.0,
265 Lstar_u[nDOF_test_element],
266 subgridError_u=0.0,
267 tau=0.0,tau0=0.0,tau1=0.0,
268 numDiff0=0.0,numDiff1=0.0,
269 jac[nSpace*nSpace],
270 jacDet,
271 jacInv[nSpace*nSpace],
272 u_grad_trial[nDOF_trial_element*nSpace],
273 u_test_dV[nDOF_trial_element],
274 u_grad_test_dV[nDOF_test_element*nSpace],
275 dV,x,y,z,xt,yt,zt,
276 G[nSpace*nSpace],G_dd_G,tr_G,norm_Rv;
277 //
278 //compute solution and gradients at quadrature points
279 //
280 ck.calculateMapping_element(eN,
281 k,
282 mesh_dof.data(),
283 mesh_l2g.data(),
284 mesh_trial_ref.data(),
285 mesh_grad_trial_ref.data(),
286 jac,
287 jacDet,
288 jacInv,
289 x,y,z);
290 ck.calculateMappingVelocity_element(eN,
291 k,
292 mesh_velocity_dof.data(),
293 mesh_l2g.data(),
294 mesh_trial_ref.data(),
295 xt,yt,zt);
296 //get the physical integration weight
297 dV = fabs(jacDet)*dV_ref.data()[k];
298 q_dV.data()[eN_k] = dV;
299 ck.calculateG(jacInv,G,G_dd_G,tr_G);
300 //get the trial function gradients
301 ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
302 //get the solution
303 ck.valFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],u);
304 //get the solution gradients
305 ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial,grad_u);
306 //precalculate test function products with integration weights
307 for (int j=0;j<nDOF_trial_element;j++)
308 {
309 u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
310 for (int I=0;I<nSpace;I++)
311 {
312 u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
313 }
314 }
315 //
316 // //calculate pde coefficients at quadrature points
317 // //
318 double Kr,dKr;
319 evaluateCoefficients(a_rowptr.data(),
320 a_colind.data(),
321 rho,
322 beta,
323 gravity.data(),
324 alpha.data()[elementMaterialTypes.data()[eN]],
325 n.data()[elementMaterialTypes.data()[eN]],
326 thetaR.data()[elementMaterialTypes.data()[eN]],
327 thetaSR.data()[elementMaterialTypes.data()[eN]],
328 &KWs.data()[elementMaterialTypes.data()[eN]*nnz],
329 u,
330 m,
331 dm,
332 f,
333 df,
334 a,
335 da,
336 as,
337 Kr,
338 dKr);
339 //
340 // //calculate time derivative at quadrature points
341 // //
342 ck.bdf(alphaBDF,
343 q_m_betaBDF.data()[eN_k],
344 m,
345 dm,
346 m_t,
347 dm_t);
348
349
350 // //update element residual
351 // //
352 for(int i=0;i<nDOF_test_element;i++)
353 {
354 int eN_k_i=eN_k*nDOF_test_element+i,
355 eN_k_i_nSpace = eN_k_i*nSpace,
356 i_nSpace=i*nSpace;
357 if (LUMPED_MASS_MATRIX==1)
358 {
359 // Lumped mass matrix contribution
360 globalResidual.data()[offset_u + stride_u*u_l2g.data()[eN*nDOF_test_element + i]] += u_test_dV[i] * m_t;
361 }
362 else
363 {
364 elementResidual_u[i] += ck.Mass_weak(m_t, u_test_dV[i]);
365 }
366
367 elementResidual_u[i] += ck.Advection_weak(f,&u_grad_test_dV[i_nSpace]) +
368 ck.Diffusion_weak(a_rowptr.data(),a_colind.data(),a,grad_u,&u_grad_test_dV[i_nSpace]);
369 //+ ck.Mass_weak(m_t,u_test_dV[i])
370 /* + */
371 // /* ck.SubgridError(subgridError_u,Lstar_u[i]) + */
372 // /* ck.NumericalDiffusion(q_numDiff_u_last[eN_k],grad_u,&u_grad_test_dV[i_nSpace]); */
373 }//i
374 // //
375 q_m.data()[eN_k] = m;
376 q_u.data()[eN_k] = u;
377 }
378 // //
379 //load element into global residual and save element residual
380 //
381 for(int i=0;i<nDOF_test_element;i++)
382 {
383 int eN_i=eN*nDOF_test_element+i;
384
385 globalResidual.data()[offset_u+stride_u*u_l2g.data()[eN_i]] += elementResidual_u[i];
386 }//i
387 }//elements
388 }
389
390 void calculateJacobian(arguments_dict& args)
391 {
392 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
393 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
394 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
395 xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
396 double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
397 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
398 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
399 xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
400 xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
401 xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
402 xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
403 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
404 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
405 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
406 xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
407 xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
408 xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
409 xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
410 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
411 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
412 int nElements_global = args.scalar<int>("nElements_global");
413 xt::pyarray<double>& ebqe_penalty_ext = args.array<double>("ebqe_penalty_ext");
414 xt::pyarray<int>& elementMaterialTypes = args.array<int>("elementMaterialTypes");
415 xt::pyarray<int>& isSeepageFace = args.array<int>("isSeepageFace");
416 xt::pyarray<int>& a_rowptr = args.array<int>("a_rowptr");
417 xt::pyarray<int>& a_colind = args.array<int>("a_colind");
418 double rho = args.scalar<double>("rho");
419 double beta = args.scalar<double>("beta");
420 xt::pyarray<double>& gravity = args.array<double>("gravity");
421 xt::pyarray<double>& alpha = args.array<double>("alpha");
422 xt::pyarray<double>& n = args.array<double>("n");
423 xt::pyarray<double>& thetaR = args.array<double>("thetaR");
424 xt::pyarray<double>& thetaSR = args.array<double>("thetaSR");
425 xt::pyarray<double>& KWs = args.array<double>("KWs");
426 double useMetrics = args.scalar<double>("useMetrics");
427 double alphaBDF = args.scalar<double>("alphaBDF");
428 int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
429 double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
430 xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
431 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
432 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
433 xt::pyarray<double>& velocity = args.array<double>("velocity");
434 xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
435 xt::pyarray<double>& cfl = args.array<double>("cfl");
436 xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
437 xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
438 xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
439 xt::pyarray<double>& globalJacobian = args.array<double>("globalJacobian");
440 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
441 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
442 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
443 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
444 xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
445 xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
446 xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
447 xt::pyarray<int>& isFluxBoundary_u = args.array<int>("isFluxBoundary_u");
448 xt::pyarray<double>& ebqe_bc_flux_ext = args.array<double>("ebqe_bc_flux_ext");
449 xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
450 int LUMPED_MASS_MATRIX = args.scalar<int>("LUMPED_MASS_MATRIX");
451 assert(a_rowptr.data()[nSpace] == nnz);
452 assert(a_rowptr.data()[nSpace] == nSpace);
453 double Ct_sge = 4.0;
454
455 //
456 //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
457 //
458 for(int eN=0;eN<nElements_global;eN++)
459 {
460 double elementJacobian_u_u[nDOF_test_element][nDOF_trial_element];
461 for (int i=0;i<nDOF_test_element;i++)
462 {
463 for (int j=0;j<nDOF_trial_element;j++)
464 {
465 elementJacobian_u_u[i][j]=0.0;
466 }
467 }
468 for (int k=0;k<nQuadraturePoints_element;k++)
469 {
470 int eN_k = eN*nQuadraturePoints_element+k, //index to a scalar at a quadrature point
471 eN_k_nSpace = eN_k*nSpace,
472 eN_nDOF_trial_element = eN*nDOF_trial_element; //index to a vector at a quadrature point
473
474 //declare local storage
475 double u=0.0,
476 grad_u[nSpace],
477 m=0.0,dm=0.0,
478 f[nSpace],df[nSpace],
479 a[nnz],da[nnz],as[nnz],
480 m_t=0.0,dm_t=0.0,
481 dpdeResidual_u_u[nDOF_trial_element],
482 Lstar_u[nDOF_test_element],
483 dsubgridError_u_u[nDOF_trial_element],
484 tau=0.0,tau0=0.0,tau1=0.0,
485 jac[nSpace*nSpace],
486 jacDet,
487 jacInv[nSpace*nSpace],
488 u_grad_trial[nDOF_trial_element*nSpace],
489 dV,
490 u_test_dV[nDOF_test_element],
491 u_grad_test_dV[nDOF_test_element*nSpace],
492 x,y,z,xt,yt,zt,
493 G[nSpace*nSpace],G_dd_G,tr_G;
494 //
495 //calculate solution and gradients at quadrature points
496 //
497 //get jacobian, etc for mapping reference element
498 ck.calculateMapping_element(eN,
499 k,
500 mesh_dof.data(),
501 mesh_l2g.data(),
502 mesh_trial_ref.data(),
503 mesh_grad_trial_ref.data(),
504 jac,
505 jacDet,
506 jacInv,
507 x,y,z);
508 ck.calculateMappingVelocity_element(eN,
509 k,
510 mesh_velocity_dof.data(),
511 mesh_l2g.data(),
512 mesh_trial_ref.data(),
513 xt,yt,zt);
514 //get the physical integration weight
515 dV = fabs(jacDet)*dV_ref.data()[k];
516 ck.calculateG(jacInv,G,G_dd_G,tr_G);
517 //get the trial function gradients
518 ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
519 //get the solution
520 ck.valFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],u);
521 //get the solution gradients
522 ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial,grad_u);
523 //precalculate test function products with integration weights
524 for (int j=0;j<nDOF_trial_element;j++)
525 {
526 u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
527 for (int I=0;I<nSpace;I++)
528 {
529 u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
530 }
531 }
532 //
533 //calculate pde coefficients and derivatives at quadrature points
534 //
535 double Kr,dKr;
536 evaluateCoefficients(a_rowptr.data(),
537 a_colind.data(),
538 rho,
539 beta,
540 gravity.data(),
541 alpha.data()[elementMaterialTypes.data()[eN]],
542 n.data()[elementMaterialTypes.data()[eN]],
543 thetaR.data()[elementMaterialTypes.data()[eN]],
544 thetaSR.data()[elementMaterialTypes.data()[eN]],
545 &KWs.data()[elementMaterialTypes.data()[eN]*nnz],
546 u,
547 m,
548 dm,
549 f,
550 df,
551 a,
552 da,
553 as,
554 Kr,
555 dKr);
556 //
557 //calculate time derivatives
558
559 ck.bdf(alphaBDF,
560 q_m_betaBDF.data()[eN_k],
561 m,
562 dm,
563 m_t,
564 dm_t);
565
566 for(int i=0;i<nDOF_test_element;i++)
567 {
568 //int eN_k_i=eN_k*nDOF_test_element+i;
569 //int eN_k_i_nSpace=eN_k_i*nSpace;
570 for(int j=0;j<nDOF_trial_element;j++)
571 {
572 //int eN_k_j=eN_k*nDOF_trial_element+j;
573 //int eN_k_j_nSpace = eN_k_j*nSpace;
574 int j_nSpace = j*nSpace;
575 int i_nSpace = i*nSpace;
576
577 elementJacobian_u_u[i][j] += ck.MassJacobian_weak(dm_t,u_trial_ref.data()[k*nDOF_trial_element+j],u_test_dV[i]) +
578 ck.AdvectionJacobian_weak(df,u_trial_ref.data()[k*nDOF_trial_element+j],&u_grad_test_dV[i_nSpace]) +
579 ck.DiffusionJacobian_weak(a_rowptr.data(),a_colind.data(),a,da,
580 grad_u,&u_grad_test_dV[i_nSpace],1.0,
581 u_trial_ref.data()[k*nDOF_trial_element+j],&u_grad_trial[j_nSpace]);
582 //+
583 // +
584 // ck.SubgridErrorJacobian(dsubgridError_u_u[j],Lstar_u[i]) +
585 // ck.NumericalDiffusionJacobian(q_numDiff_u_last[eN_k],&u_grad_trial[j_nSpace],&u_grad_test_dV[i_nSpace]);
586 }//j
587 }//i
588 }//k
589 // //
590 //load into element Jacobian into global Jacobian
591 //
592 for (int i=0;i<nDOF_test_element;i++)
593 {
594 int eN_i = eN*nDOF_test_element+i;
595 for (int j=0;j<nDOF_trial_element;j++)
596 {
597 int eN_i_j = eN_i*nDOF_trial_element+j;
598 globalJacobian.data()[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_u_u[eN_i_j]] += elementJacobian_u_u[i][j];
599 }//j
600 }//i
601 }//elements
602
603 }//computeJacobian
604
605 };//ADR
606
607
608
609 inline ADR_base* newADR(int nSpaceIn,
610 int nQuadraturePoints_elementIn,
611 int nDOF_mesh_trial_elementIn,
612 int nDOF_trial_elementIn,
613 int nDOF_test_elementIn,
614 int nQuadraturePoints_elementBoundaryIn,
615 int CompKernelFlag)
616 {
617 if (nSpaceIn == 1)
618 {
619 assert(nSpaceIn == 1);
620 return proteus::chooseAndAllocateDiscretization1D<ADR_base,ADR,CompKernel>(nSpaceIn,
621 nQuadraturePoints_elementIn,
622 nDOF_mesh_trial_elementIn,
623 nDOF_trial_elementIn,
624 nDOF_test_elementIn,
625 nQuadraturePoints_elementBoundaryIn,
626 CompKernelFlag);
627 }
628 else if (nSpaceIn == 2)
629 {
630 assert(nSpaceIn == 3);
631 return proteus::chooseAndAllocateDiscretization2D<ADR_base,ADR,CompKernel>(nSpaceIn,
632 nQuadraturePoints_elementIn,
633 nDOF_mesh_trial_elementIn,
634 nDOF_trial_elementIn,
635 nDOF_test_elementIn,
636 nQuadraturePoints_elementBoundaryIn,
637 CompKernelFlag);
638 }
639 else
640 {
641 assert(nSpaceIn == 3);
642 return proteus::chooseAndAllocateDiscretization<ADR_base,ADR,CompKernel>(nSpaceIn,
643 nQuadraturePoints_elementIn,
644 nDOF_mesh_trial_elementIn,
645 nDOF_trial_elementIn,
646 nDOF_test_elementIn,
647 nQuadraturePoints_elementBoundaryIn,
648 CompKernelFlag);
649 }
650 }
651}//richards
652}//proteus
653#endif
n
Int n
Definition Headers.h:28
u
Double u
Definition Headers.h:89
z
Double * z
Definition Headers.h:49
cE
#define cE
Definition NCLS3P.h:10
proteus::richards::ADR_base
Definition richards/ADR.h:22
proteus::richards::ADR_base::calculateResidual
virtual void calculateResidual(arguments_dict &args)=0
proteus::richards::ADR_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)=0
proteus::richards::ADR_base::~ADR_base
virtual ~ADR_base()
Definition richards/ADR.h:25
proteus::richards::ADR::nDOF_test_X_trial_element
const int nDOF_test_X_trial_element
Definition richards/ADR.h:41
proteus::richards::ADR::calculateResidual
void calculateResidual(arguments_dict &args)
Definition richards/ADR.h:116
proteus::richards::ADR::calculateJacobian
void calculateJacobian(arguments_dict &args)
Definition richards/ADR.h:390
proteus::richards::ADR::evaluateCoefficients
void evaluateCoefficients(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], 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)
Definition richards/ADR.h:50
df
double df(double C, double b, double a, int q, int r)
Definition analyticalSolutions.c:2209
D
#define D(x, y)
Definition jf.h:9
nnz
#define nnz
Definition m_comp_co2.h:19
proteus::richards::newADR
ADR_base * newADR(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/ADR.h:609
proteus
Definition ADR.h:19
proteus::chooseAndAllocateDiscretization
Model_Base * chooseAndAllocateDiscretization(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nDOF_v_trial_elementIn, int nDOF_v_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition MixedModelFactory.h:48
proteus::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
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