proteus  1.8.1
C/C++/Fortran libraries
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 "CompKernel.h"
6 #include "ModelFactory.h"
7 #include "../mprans/ArgumentsDict.h"
8 #include "xtensor-python/pyarray.hpp"
9 #define nnz nSpace
10 
11 namespace py = pybind11;
12 
13 namespace proteus
14 {
15  class Richards_base
16  {
17  //The base class defining the interface
18  public:
19  virtual ~Richards_base() {}
20  virtual void calculateResidual(arguments_dict& args)=0;
21  virtual void calculateJacobian(arguments_dict& args)=0;
22  };
23 
24  template<class CompKernelType,
25  int nSpace,
26  int nQuadraturePoints_element,
27  int nDOF_mesh_trial_element,
28  int nDOF_trial_element,
29  int nDOF_test_element,
30  int nQuadraturePoints_elementBoundary>
31  class Richards : public Richards_base
32  {
33  public:
34  const int nDOF_test_X_trial_element;
35  CompKernelType ck;
36  Richards():
37  nDOF_test_X_trial_element(nDOF_test_element*nDOF_trial_element),
38  ck()
39  {}
40  inline
41  void evaluateCoefficients(const int rowptr[nSpace],
42  const int colind[nnz],
43  const double rho,
44  const double beta,
45  const double gravity[nSpace],
46  const double alpha,
47  const double n_vg,
48  const double thetaR,
49  const double thetaSR,
50  const double KWs[nnz],
51  const double& u,
52  double& m,
53  double& dm,
54  double f[nSpace],
55  double df[nSpace],
56  double a[nnz],
57  double da[nnz])
58  {
59  const int nSpace2 = nSpace * nSpace;
60  double psiC;
61  double pcBar;
62  double pcBar_n;
63  double pcBar_nM1;
64  double pcBar_nM2;
65  double onePlus_pcBar_n;
66  double sBar;
67  double sqrt_sBar;
68  double DsBar_DpsiC;
69  double thetaW;
70  double DthetaW_DpsiC;
71  double vBar;
72  double vBar2;
73  double DvBar_DpsiC;
74  double KWr;
75  double DKWr_DpsiC;
76  double rho2 = rho * rho;
77  double thetaS;
78  double rhom;
79  double drhom;
80  double m_vg;
81  double pcBarStar;
82  double sqrt_sBarStar;
83 
84  psiC = -u;
85  m_vg = 1.0 - 1.0 / n_vg;
86  thetaS = thetaR + thetaSR;
87  if (psiC > 0.0)
88  {
89  pcBar = alpha * psiC;
90  pcBarStar = pcBar;
91  if (pcBar < 1.0e-8)
92  pcBarStar = 1.0e-8;
93  pcBar_nM2 = pow(pcBarStar, n_vg - 2);
94  pcBar_nM1 = pcBar_nM2 * pcBar;
95  pcBar_n = pcBar_nM1 * pcBar;
96  onePlus_pcBar_n = 1.0 + pcBar_n;
97 
98  sBar = pow(onePlus_pcBar_n, -m_vg);
99  /* using -mn = 1-n */
100  DsBar_DpsiC = alpha * (1.0 - n_vg) * (sBar/onePlus_pcBar_n)*pcBar_nM1;
101 
102  vBar = 1.0-pcBar_nM1*sBar;
103  vBar2 = vBar*vBar;
104  DvBar_DpsiC = -alpha*(n_vg-1.0)*pcBar_nM2*sBar - pcBar_nM1*DsBar_DpsiC;
105 
106  thetaW = thetaSR*sBar + thetaR;
107  DthetaW_DpsiC = thetaSR * DsBar_DpsiC;
108 
109  sqrt_sBar = sqrt(sBar);
110  sqrt_sBarStar = sqrt_sBar;
111  if (sqrt_sBar < 1.0e-8)
112  sqrt_sBarStar = 1.0e-8;
113  KWr= sqrt_sBar*vBar2;
114  DKWr_DpsiC= ((0.5/sqrt_sBarStar)*DsBar_DpsiC*vBar2
115  +
116  2.0*sqrt_sBar*vBar*DvBar_DpsiC);
117  }
118  else
119  {
120  thetaW = thetaS;
121  DthetaW_DpsiC = 0.0;
122  KWr = 1.0;
123  DKWr_DpsiC = 0.0;
124  }
125  //slight compressibility
126  rhom = rho*exp(beta*u);
127  drhom = beta*rhom;
128  m = rhom*thetaW;
129  dm = -rhom*DthetaW_DpsiC+drhom*thetaW;
130  for (int I=0;I<nSpace;I++)
131  {
132  f[I] = 0.0;
133  df[I] = 0.0;
134  for (int ii=rowptr[I]; ii < rowptr[I+1]; ii++)
135  {
136  f[I] += rho2*KWr*KWs[ii]*gravity[colind[ii]];
137  df[I] += -rho2*DKWr_DpsiC*KWs[ii]*gravity[colind[ii]];
138  a[ii] = rho*KWr*KWs[ii];
139  da[ii] = -rho*DKWr_DpsiC*KWs[ii];
140  }
141  }
142  }
143 
144  inline
145  void calculateSubgridError_tau(const double& elementDiameter,
146  const double& dmt,
147  const double dH[nSpace],
148  double& cfl,
149  double& tau)
150  {
151  double h,nrm_v,oneByAbsdt;
152  h = elementDiameter;
153  nrm_v=0.0;
154  for(int I=0;I<nSpace;I++)
155  nrm_v+=dH[I]*dH[I];
156  nrm_v = sqrt(nrm_v);
157  cfl = nrm_v/h;
158  oneByAbsdt = fabs(dmt);
159  tau = 1.0/(2.0*nrm_v/h + oneByAbsdt + 1.0e-8);
160  }
161 
162  inline
163  void calculateSubgridError_tau(const double& Ct_sge,
164  const double G[nSpace*nSpace],
165  const double& A0,
166  const double Ai[nSpace],
167  double& tau_v,
168  double& cfl)
169  {
170  double v_d_Gv=0.0;
171  for(int I=0;I<nSpace;I++)
172  for (int J=0;J<nSpace;J++)
173  v_d_Gv += Ai[I]*G[I*nSpace+J]*Ai[J];
174 
175  tau_v = 1.0/sqrt(Ct_sge*A0*A0 + v_d_Gv);
176  }
177 
178  inline
179  void calculateNumericalDiffusion(const double& shockCapturingDiffusion,
180  const double& elementDiameter,
181  const double& strong_residual,
182  const double grad_u[nSpace],
183  double& numDiff)
184  {
185  double h,
186  num,
187  den,
188  n_grad_u;
189  h = elementDiameter;
190  n_grad_u = 0.0;
191  for (int I=0;I<nSpace;I++)
192  n_grad_u += grad_u[I]*grad_u[I];
193  num = shockCapturingDiffusion*0.5*h*fabs(strong_residual);
194  den = sqrt(n_grad_u) + 1.0e-8;
195  numDiff = num/den;
196  }
197 
198  inline
199  void exteriorNumericalFlux(const double& bc_flux,
200  int rowptr[nSpace],
201  int colind[nnz],
202  int isSeepageFace,
203  int& isDOFBoundary,
204  double n[nSpace],
205  double bc_u,
206  double K[nnz],
207  double grad_psi[nSpace],
208  double u,
209  double K_rho_g[nSpace],
210  double penalty,
211  double& flux)
212  {
213  double v_I,bc_u_seepage=0.0;
214  if (isSeepageFace || isDOFBoundary)
215  {
216  flux = 0.0;
217  for(int I=0;I<nSpace;I++)
218  {
219  //gravity
220  v_I = K_rho_g[I];
221  //pressure head
222  for(int m=rowptr[I];m<rowptr[I+1];m++)
223  {
224  v_I -= K[m]*grad_psi[colind[m]];
225  }
226  flux += v_I*n[I];
227  }
228  if (isSeepageFace)
229  bc_u = bc_u_seepage;
230  flux += penalty*(u-bc_u);
231  if (isSeepageFace)
232  {
233  if (flux > 0.0)
234  {
235  isDOFBoundary = 1;
236  bc_u = bc_u_seepage;
237  }
238  else
239  {
240  isDOFBoundary = 0;
241  flux = 0.0;
242  }
243  }
244  /* //set DOF flag and flux correctly if seepage face */
245  /* if (isSeepageFace) */
246  /* { */
247  /* if (flux < 0.0 || u < bc_u_seepage) */
248  /* { */
249  /* isDOFBoundary = 0; */
250  /* flux = 0.0; */
251  /* } */
252  /* else */
253  /* { */
254  /* isDOFBoundary = 1; */
255  /* bc_u = bc_u_seepage; */
256  /* } */
257  /* } */
258  /* //Dirichlet penalty */
259  /* if (isDOFBoundary) */
260  /* flux += penalty*(u-bc_u); */
261  }
262  else
263  flux = bc_flux;
264  }
265 
266  void exteriorNumericalFluxJacobian(const int rowptr[nSpace],
267  const int colind[nnz],
268  const int isDOFBoundary,
269  const double n[nSpace],
270  const double K[nnz],
271  const double dK[nnz],
272  const double grad_psi[nSpace],
273  const double grad_v[nSpace],
274  const double dK_rho_g[nSpace],
275  const double v,
276  const double penalty,
277  double& fluxJacobian)
278  {
279  if (isDOFBoundary)
280  {
281  fluxJacobian = 0.0;
282  for(int I=0;I<nSpace;I++)
283  {
284  //gravity
285  fluxJacobian += dK_rho_g[I]*v*n[I];
286  //pressure head
287  for(int m=rowptr[I]; m<rowptr[I+1]; m++)
288  {
289  fluxJacobian -= (K[m]*grad_v[colind[m]] + dK[m]*v*grad_psi[colind[m]])*n[I];
290  }
291  }
292  //Dirichlet penalty
293  fluxJacobian += penalty*v;
294  }
295  else
296  fluxJacobian = 0.0;
297  }
298 
299  void calculateResidual(arguments_dict& args)
300  {
301  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
302  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
303  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
304  xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
305  double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
306  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
307  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
308  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
309  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
310  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
311  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
312  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
313  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
314  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
315  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
316  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
317  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
318  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
319  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
320  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
321  int nElements_global = args.scalar<int>("nElements_global");
322  xt::pyarray<double>& ebqe_penalty_ext = args.array<double>("ebqe_penalty_ext");
323  xt::pyarray<int>& elementMaterialTypes = args.array<int>("elementMaterialTypes");
324  xt::pyarray<int>& isSeepageFace = args.array<int>("isSeepageFace");
325  xt::pyarray<int>& a_rowptr = args.array<int>("a_rowptr");
326  xt::pyarray<int>& a_colind = args.array<int>("a_colind");
327  double rho = args.scalar<double>("rho");
328  double beta = args.scalar<double>("beta");
329  xt::pyarray<double>& gravity = args.array<double>("gravity");
330  xt::pyarray<double>& alpha = args.array<double>("alpha");
331  xt::pyarray<double>& n = args.array<double>("n");
332  xt::pyarray<double>& thetaR = args.array<double>("thetaR");
333  xt::pyarray<double>& thetaSR = args.array<double>("thetaSR");
334  xt::pyarray<double>& KWs = args.array<double>("KWs");
335  double useMetrics = args.scalar<double>("useMetrics");
336  double alphaBDF = args.scalar<double>("alphaBDF");
337  int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
338  double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
339  double sc_uref = args.scalar<double>("sc_uref");
340  double sc_alpha = args.scalar<double>("sc_alpha");
341  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
342  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
343  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
344  xt::pyarray<double>& u_dof_old = args.array<double>("u_dof_old");
345  xt::pyarray<double>& velocity = args.array<double>("velocity");
346  xt::pyarray<double>& q_m = args.array<double>("q_m");
347  xt::pyarray<double>& q_u = args.array<double>("q_u");
348  xt::pyarray<double>& q_dV = args.array<double>("q_dV");
349  xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
350  xt::pyarray<double>& cfl = args.array<double>("cfl");
351  xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
352  xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
353  int offset_u = args.scalar<int>("offset_u");
354  int stride_u = args.scalar<int>("stride_u");
355  xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
356  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
357  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
358  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
359  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
360  xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
361  xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
362  xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
363  xt::pyarray<int>& isFluxBoundary_u = args.array<int>("isFluxBoundary_u");
364  xt::pyarray<double>& ebqe_bc_flux_ext = args.array<double>("ebqe_bc_flux_ext");
365  xt::pyarray<double>& ebqe_phi = args.array<double>("ebqe_phi");
366  double epsFact = args.scalar<double>("epsFact");
367  xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
368  xt::pyarray<double>& ebqe_flux = args.array<double>("ebqe_flux");
369  assert(a_rowptr.data()[nSpace] == nnz);
370  assert(a_rowptr.data()[nSpace] == nSpace);
371  //cek should this be read in?
372  double Ct_sge = 4.0;
373 
374  //loop over elements to compute volume integrals and load them into element and global residual
375  //
376  //eN is the element index
377  //eN_k is the quadrature point index for a scalar
378  //eN_k_nSpace is the quadrature point index for a vector
379  //eN_i is the element test function index
380  //eN_j is the element trial function index
381  //eN_k_j is the quadrature point index for a trial function
382  //eN_k_i is the quadrature point index for a trial function
383  for(int eN=0;eN<nElements_global;eN++)
384  {
385  //declare local storage for element residual and initialize
386  register double elementResidual_u[nDOF_test_element];
387  for (int i=0;i<nDOF_test_element;i++)
388  {
389  elementResidual_u[i]=0.0;
390  }//i
391  //loop over quadrature points and compute integrands
392  for (int k=0;k<nQuadraturePoints_element;k++)
393  {
394  //compute indeces and declare local storage
395  register int eN_k = eN*nQuadraturePoints_element+k,
396  eN_k_nSpace = eN_k*nSpace,
397  eN_nDOF_trial_element = eN*nDOF_trial_element;
398  register double u=0.0,grad_u[nSpace],grad_u_old[nSpace],
399  m=0.0,dm=0.0,
400  f[nSpace],df[nSpace],
401  a[nnz],da[nnz],
402  m_t=0.0,dm_t=0.0,
403  pdeResidual_u=0.0,
404  Lstar_u[nDOF_test_element],
405  subgridError_u=0.0,
406  tau=0.0,tau0=0.0,tau1=0.0,
407  numDiff0=0.0,numDiff1=0.0,
408  jac[nSpace*nSpace],
409  jacDet,
410  jacInv[nSpace*nSpace],
411  u_grad_trial[nDOF_trial_element*nSpace],
412  u_test_dV[nDOF_trial_element],
413  u_grad_test_dV[nDOF_test_element*nSpace],
414  dV,x,y,z,xt,yt,zt,
415  G[nSpace*nSpace],G_dd_G,tr_G,norm_Rv;
416  //
417  //compute solution and gradients at quadrature points
418  //
419  ck.calculateMapping_element(eN,
420  k,
421  mesh_dof.data(),
422  mesh_l2g.data(),
423  mesh_trial_ref.data(),
424  mesh_grad_trial_ref.data(),
425  jac,
426  jacDet,
427  jacInv,
428  x,y,z);
429  ck.calculateMappingVelocity_element(eN,
430  k,
431  mesh_velocity_dof.data(),
432  mesh_l2g.data(),
433  mesh_trial_ref.data(),
434  xt,yt,zt);
435  //get the physical integration weight
436  dV = fabs(jacDet)*dV_ref.data()[k];
437  q_dV.data()[eN_k] = dV;
438  ck.calculateG(jacInv,G,G_dd_G,tr_G);
439  //get the trial function gradients
440  ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
441  //get the solution
442  ck.valFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],u);
443  //get the solution gradients
444  ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial,grad_u);
445  //precalculate test function products with integration weights
446  for (int j=0;j<nDOF_trial_element;j++)
447  {
448  u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
449  for (int I=0;I<nSpace;I++)
450  {
451  u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
452  }
453  }
454  //
455  //calculate pde coefficients at quadrature points
456  //
457  evaluateCoefficients(a_rowptr.data(),
458  a_colind.data(),
459  rho,
460  beta,
461  gravity.data(),
462  alpha.data()[elementMaterialTypes.data()[eN]],
463  n.data()[elementMaterialTypes.data()[eN]],
464  thetaR.data()[elementMaterialTypes.data()[eN]],
465  thetaSR.data()[elementMaterialTypes.data()[eN]],
466  &KWs.data()[elementMaterialTypes.data()[eN]*nnz],
467  u,
468  m,
469  dm,
470  f,
471  df,
472  a,
473  da);
474  //
475  //calculate time derivative at quadrature points
476  //
477  ck.bdf(alphaBDF,
478  q_m_betaBDF.data()[eN_k],
479  m,
480  dm,
481  m_t,
482  dm_t);
483  /* // */
484  /* //calculate subgrid error (strong residual and adjoint) */
485  /* // */
486  /* //calculate strong residual */
487  /* pdeResidual_u = ck.Mass_strong(m_t) + ck.Advection_strong(df,grad_u); */
488  /* //calculate adjoint */
489  /* for (int i=0;i<nDOF_test_element;i++) */
490  /* { */
491  /* // register int eN_k_i_nSpace = (eN_k*nDOF_trial_element+i)*nSpace; */
492  /* // Lstar_u[i] = ck.Advection_adjoint(df,&u_grad_test_dV[eN_k_i_nSpace]); */
493  /* register int i_nSpace = i*nSpace; */
494  /* Lstar_u[i] = ck.Advection_adjoint(df,&u_grad_test_dV[i_nSpace]); */
495  /* } */
496  /* //calculate tau and tau*Res */
497  /* calculateSubgridError_tau(elementDiameter[eN],dm_t,df,cfl[eN_k],tau0); */
498  /* calculateSubgridError_tau(Ct_sge, */
499  /* G, */
500  /* dm_t, */
501  /* df, */
502  /* tau1, */
503  /* cfl[eN_k]); */
504 
505  /* tau = useMetrics*tau1+(1.0-useMetrics)*tau0; */
506 
507  /* subgridError_u = -tau*pdeResidual_u; */
508  /* // */
509  /* //calculate shock capturing diffusion */
510  /* // */
511 
512 
513  /* ck.calculateNumericalDiffusion(shockCapturingDiffusion,elementDiameter[eN],pdeResidual_u,grad_u,numDiff0); */
514  /* //ck.calculateNumericalDiffusion(shockCapturingDiffusion,G,pdeResidual_u,grad_u_old,numDiff1); */
515  /* ck.calculateNumericalDiffusion(shockCapturingDiffusion,sc_uref, sc_alpha,G,G_dd_G,pdeResidual_u,grad_u,numDiff1); */
516  /* q_numDiff_u[eN_k] = useMetrics*numDiff1+(1.0-useMetrics)*numDiff0; */
517  //std::cout<<tau<<" "<<q_numDiff_u[eN_k]<<std::endl;
518  //
519  //update element residual
520  //
521  for(int i=0;i<nDOF_test_element;i++)
522  {
523  register int eN_k_i=eN_k*nDOF_test_element+i,
524  eN_k_i_nSpace = eN_k_i*nSpace,
525  i_nSpace=i*nSpace;
526 
527  elementResidual_u[i] += ck.Mass_weak(m_t,u_test_dV[i]) +
528  ck.Advection_weak(f,&u_grad_test_dV[i_nSpace]) +
529  ck.Diffusion_weak(a_rowptr.data(),a_colind.data(),a,grad_u,&u_grad_test_dV[i_nSpace]);
530  /* + */
531  /* ck.SubgridError(subgridError_u,Lstar_u[i]) + */
532  /* ck.NumericalDiffusion(q_numDiff_u_last[eN_k],grad_u,&u_grad_test_dV[i_nSpace]); */
533  }//i
534  //
535  q_m.data()[eN_k] = m;
536  q_u.data()[eN_k] = u;
537  }
538  //
539  //load element into global residual and save element residual
540  //
541  for(int i=0;i<nDOF_test_element;i++)
542  {
543  register int eN_i=eN*nDOF_test_element+i;
544 
545  globalResidual.data()[offset_u+stride_u*u_l2g.data()[eN_i]] += elementResidual_u[i];
546  }//i
547  }//elements
548  //
549  //loop over exterior element boundaries to calculate surface integrals and load into element and global residuals
550  //
551  //ebNE is the Exterior element boundary INdex
552  //ebN is the element boundary INdex
553  //eN is the element index
554  for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
555  {
556  register int ebN = exteriorElementBoundariesArray.data()[ebNE],
557  eN = elementBoundaryElementsArray.data()[ebN*2+0],
558  ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
559  eN_nDOF_trial_element = eN*nDOF_trial_element;
560  register double elementResidual_u[nDOF_test_element];
561  for (int i=0;i<nDOF_test_element;i++)
562  {
563  elementResidual_u[i]=0.0;
564  }
565  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
566  {
567  register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
568  ebNE_kb_nSpace = ebNE_kb*nSpace,
569  ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
570  ebN_local_kb_nSpace = ebN_local_kb*nSpace;
571  register double u_ext=0.0,
572  grad_u_ext[nSpace],
573  m_ext=0.0,
574  dm_ext=0.0,
575  f_ext[nSpace],
576  df_ext[nSpace],
577  a_ext[nnz],
578  da_ext[nnz],
579  flux_ext=0.0,
580  bc_u_ext=0.0,
581  bc_grad_u_ext[nSpace],
582  bc_m_ext=0.0,
583  bc_dm_ext=0.0,
584  bc_f_ext[nSpace],
585  bc_df_ext[nSpace],
586  bc_a_ext[nnz],
587  bc_da_ext[nnz],
588  jac_ext[nSpace*nSpace],
589  jacDet_ext,
590  jacInv_ext[nSpace*nSpace],
591  boundaryJac[nSpace*(nSpace-1)],
592  metricTensor[(nSpace-1)*(nSpace-1)],
593  metricTensorDetSqrt,
594  dS,
595  u_test_dS[nDOF_test_element],
596  u_grad_trial_trace[nDOF_trial_element*nSpace],
597  normal[3],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
598  G[nSpace*nSpace],G_dd_G,tr_G;
599  //
600  //calculate the solution and gradients at quadrature points
601  //
602  //compute information about mapping from reference element to physical element
603  ck.calculateMapping_elementBoundary(eN,
604  ebN_local,
605  kb,
606  ebN_local_kb,
607  mesh_dof.data(),
608  mesh_l2g.data(),
609  mesh_trial_trace_ref.data(),
610  mesh_grad_trial_trace_ref.data(),
611  boundaryJac_ref.data(),
612  jac_ext,
613  jacDet_ext,
614  jacInv_ext,
615  boundaryJac,
616  metricTensor,
617  metricTensorDetSqrt,
618  normal_ref.data(),
619  normal,
620  x_ext,y_ext,z_ext);
621  ck.calculateMappingVelocity_elementBoundary(eN,
622  ebN_local,
623  kb,
624  ebN_local_kb,
625  mesh_velocity_dof.data(),
626  mesh_l2g.data(),
627  mesh_trial_trace_ref.data(),
628  xt_ext,yt_ext,zt_ext,
629  normal,
630  boundaryJac,
631  metricTensor,
632  integralScaling);
633  dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref.data()[kb];
634  //get the metric tensor
635  //cek todo use symmetry
636  ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
637  //compute shape and solution information
638  //shape
639  ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,u_grad_trial_trace);
640  //solution and gradient
641  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);
642  ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial_trace,grad_u_ext);
643  //precalculate test function products with integration weights
644  for (int j=0;j<nDOF_trial_element;j++)
645  {
646  u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
647  }
648  //
649  //load the boundary values
650  //
651  bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+(1-isDOFBoundary_u.data()[ebNE_kb])*u_ext;
652  //
653  //calculate the pde coefficients using the solution and the boundary values for the solution
654  //
655  evaluateCoefficients(a_rowptr.data(),
656  a_colind.data(),
657  rho,
658  beta,
659  gravity.data(),
660  alpha.data()[elementMaterialTypes.data()[eN]],
661  n.data()[elementMaterialTypes.data()[eN]],
662  thetaR.data()[elementMaterialTypes.data()[eN]],
663  thetaSR.data()[elementMaterialTypes.data()[eN]],
664  &KWs.data()[elementMaterialTypes.data()[eN]*nnz],
665  u_ext,
666  m_ext,
667  dm_ext,
668  f_ext,
669  df_ext,
670  a_ext,
671  da_ext);
672  evaluateCoefficients(a_rowptr.data(),
673  a_colind.data(),
674  rho,
675  beta,
676  gravity.data(),
677  alpha.data()[elementMaterialTypes.data()[eN]],
678  n.data()[elementMaterialTypes.data()[eN]],
679  thetaR.data()[elementMaterialTypes.data()[eN]],
680  thetaSR.data()[elementMaterialTypes.data()[eN]],
681  &KWs.data()[elementMaterialTypes.data()[eN]*nnz],
682  bc_u_ext,
683  bc_m_ext,
684  bc_dm_ext,
685  bc_f_ext,
686  bc_df_ext,
687  bc_a_ext,
688  bc_da_ext);
689  //
690  //calculate the numerical fluxes
691  //
692  exteriorNumericalFlux(ebqe_bc_flux_ext[ebNE_kb],
693  a_rowptr.data(),
694  a_colind.data(),
695  isSeepageFace.data()[ebNE],//tricky, this is a face flag not face quad
696  isDOFBoundary_u.data()[ebNE_kb],
697  normal,
698  bc_u_ext,
699  a_ext,
700  grad_u_ext,
701  u_ext,
702  f_ext,
703  ebqe_penalty_ext.data()[ebNE_kb],// penalty,
704  flux_ext);
705  ebqe_flux.data()[ebNE_kb] = flux_ext;
706  ebqe_u.data()[ebNE_kb] = u_ext;
707  //
708  //update residuals
709  //
710  for (int i=0;i<nDOF_test_element;i++)
711  {
712  int ebNE_kb_i = ebNE_kb*nDOF_test_element+i;
713 
714  elementResidual_u[i] += ck.ExteriorElementBoundaryFlux(flux_ext,u_test_dS[i]);
715  }//i
716  }//kb
717  //
718  //update the element and global residual storage
719  //
720  for (int i=0;i<nDOF_test_element;i++)
721  {
722  int eN_i = eN*nDOF_test_element+i;
723 
724  globalResidual.data()[offset_u+stride_u*u_l2g.data()[eN_i]] += elementResidual_u[i];
725  }//i
726  }//ebNE
727  }
728 
729  void calculateJacobian(arguments_dict& args)
730  {
731  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
732  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
733  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
734  xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
735  double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
736  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
737  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
738  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
739  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
740  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
741  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
742  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
743  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
744  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
745  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
746  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
747  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
748  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
749  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
750  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
751  int nElements_global = args.scalar<int>("nElements_global");
752  xt::pyarray<double>& ebqe_penalty_ext = args.array<double>("ebqe_penalty_ext");
753  xt::pyarray<int>& elementMaterialTypes = args.array<int>("elementMaterialTypes");
754  xt::pyarray<int>& isSeepageFace = args.array<int>("isSeepageFace");
755  xt::pyarray<int>& a_rowptr = args.array<int>("a_rowptr");
756  xt::pyarray<int>& a_colind = args.array<int>("a_colind");
757  double rho = args.scalar<double>("rho");
758  double beta = args.scalar<double>("beta");
759  xt::pyarray<double>& gravity = args.array<double>("gravity");
760  xt::pyarray<double>& alpha = args.array<double>("alpha");
761  xt::pyarray<double>& n = args.array<double>("n");
762  xt::pyarray<double>& thetaR = args.array<double>("thetaR");
763  xt::pyarray<double>& thetaSR = args.array<double>("thetaSR");
764  xt::pyarray<double>& KWs = args.array<double>("KWs");
765  double useMetrics = args.scalar<double>("useMetrics");
766  double alphaBDF = args.scalar<double>("alphaBDF");
767  int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
768  double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
769  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
770  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
771  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
772  xt::pyarray<double>& velocity = args.array<double>("velocity");
773  xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
774  xt::pyarray<double>& cfl = args.array<double>("cfl");
775  xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
776  xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
777  xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
778  xt::pyarray<double>& globalJacobian = args.array<double>("globalJacobian");
779  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
780  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
781  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
782  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
783  xt::pyarray<double>& ebqe_velocity_ext = args.array<double>("ebqe_velocity_ext");
784  xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
785  xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
786  xt::pyarray<int>& isFluxBoundary_u = args.array<int>("isFluxBoundary_u");
787  xt::pyarray<double>& ebqe_bc_flux_ext = args.array<double>("ebqe_bc_flux_ext");
788  xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
789  assert(a_rowptr.data()[nSpace] == nnz);
790  assert(a_rowptr.data()[nSpace] == nSpace);
791  double Ct_sge = 4.0;
792 
793  //
794  //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
795  //
796  for(int eN=0;eN<nElements_global;eN++)
797  {
798  register double elementJacobian_u_u[nDOF_test_element][nDOF_trial_element];
799  for (int i=0;i<nDOF_test_element;i++)
800  {
801  for (int j=0;j<nDOF_trial_element;j++)
802  {
803  elementJacobian_u_u[i][j]=0.0;
804  }
805  }
806  for (int k=0;k<nQuadraturePoints_element;k++)
807  {
808  int eN_k = eN*nQuadraturePoints_element+k, //index to a scalar at a quadrature point
809  eN_k_nSpace = eN_k*nSpace,
810  eN_nDOF_trial_element = eN*nDOF_trial_element; //index to a vector at a quadrature point
811 
812  //declare local storage
813  register double u=0.0,
814  grad_u[nSpace],
815  m=0.0,dm=0.0,
816  f[nSpace],df[nSpace],
817  a[nnz],da[nnz],
818  m_t=0.0,dm_t=0.0,
819  dpdeResidual_u_u[nDOF_trial_element],
820  Lstar_u[nDOF_test_element],
821  dsubgridError_u_u[nDOF_trial_element],
822  tau=0.0,tau0=0.0,tau1=0.0,
823  jac[nSpace*nSpace],
824  jacDet,
825  jacInv[nSpace*nSpace],
826  u_grad_trial[nDOF_trial_element*nSpace],
827  dV,
828  u_test_dV[nDOF_test_element],
829  u_grad_test_dV[nDOF_test_element*nSpace],
830  x,y,z,xt,yt,zt,
831  G[nSpace*nSpace],G_dd_G,tr_G;
832  //
833  //calculate solution and gradients at quadrature points
834  //
835  //get jacobian, etc for mapping reference element
836  ck.calculateMapping_element(eN,
837  k,
838  mesh_dof.data(),
839  mesh_l2g.data(),
840  mesh_trial_ref.data(),
841  mesh_grad_trial_ref.data(),
842  jac,
843  jacDet,
844  jacInv,
845  x,y,z);
846  ck.calculateMappingVelocity_element(eN,
847  k,
848  mesh_velocity_dof.data(),
849  mesh_l2g.data(),
850  mesh_trial_ref.data(),
851  xt,yt,zt);
852  //get the physical integration weight
853  dV = fabs(jacDet)*dV_ref.data()[k];
854  ck.calculateG(jacInv,G,G_dd_G,tr_G);
855  //get the trial function gradients
856  ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
857  //get the solution
858  ck.valFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],u);
859  //get the solution gradients
860  ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial,grad_u);
861  //precalculate test function products with integration weights
862  for (int j=0;j<nDOF_trial_element;j++)
863  {
864  u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
865  for (int I=0;I<nSpace;I++)
866  {
867  u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
868  }
869  }
870  //
871  //calculate pde coefficients and derivatives at quadrature points
872  //
873  evaluateCoefficients(a_rowptr.data(),
874  a_colind.data(),
875  rho,
876  beta,
877  gravity.data(),
878  alpha.data()[elementMaterialTypes.data()[eN]],
879  n.data()[elementMaterialTypes.data()[eN]],
880  thetaR.data()[elementMaterialTypes.data()[eN]],
881  thetaSR.data()[elementMaterialTypes.data()[eN]],
882  &KWs.data()[elementMaterialTypes.data()[eN]*nnz],
883  u,
884  m,
885  dm,
886  f,
887  df,
888  a,
889  da);
890  //
891  //calculate time derivatives
892  //
893  ck.bdf(alphaBDF,
894  q_m_betaBDF.data()[eN_k],
895  m,
896  dm,
897  m_t,
898  dm_t);
899  // //
900  // //calculate subgrid error contribution to the Jacobian (strong residual, adjoint, jacobian of strong residual)
901  // //
902  // //calculate the adjoint times the test functions
903  // for (int i=0;i<nDOF_test_element;i++)
904  // {
905  // // int eN_k_i_nSpace = (eN_k*nDOF_trial_element+i)*nSpace;
906  // // Lstar_u[i]=ck.Advection_adjoint(df,&u_grad_test_dV[eN_k_i_nSpace]);
907  // register int i_nSpace = i*nSpace;
908  // Lstar_u[i]=ck.Advection_adjoint(df,&u_grad_test_dV[i_nSpace]);
909  // }
910  // //calculate the Jacobian of strong residual
911  // for (int j=0;j<nDOF_trial_element;j++)
912  // {
913  // //int eN_k_j=eN_k*nDOF_trial_element+j;
914  // //int eN_k_j_nSpace = eN_k_j*nSpace;
915  // int j_nSpace = j*nSpace;
916  // dpdeResidual_u_u[j]= ck.MassJacobian_strong(dm_t,u_trial_ref[k*nDOF_trial_element+j]) +
917  // ck.AdvectionJacobian_strong(df,&u_grad_trial[j_nSpace]);
918  // }
919  // //tau and tau*Res
920  // calculateSubgridError_tau(elementDiameter[eN],
921  // dm_t,
922  // df,
923  // cfl[eN_k],
924  // tau0);
925  // calculateSubgridError_tau(Ct_sge,
926  // G,
927  // dm_t,
928  // df,
929  // tau1,
930  // cfl[eN_k]);
931  // tau = useMetrics*tau1+(1.0-useMetrics)*tau0;
932  // for(int j=0;j<nDOF_trial_element;j++)
933  // dsubgridError_u_u[j] = -tau*dpdeResidual_u_u[j];
934  for(int i=0;i<nDOF_test_element;i++)
935  {
936  //int eN_k_i=eN_k*nDOF_test_element+i;
937  //int eN_k_i_nSpace=eN_k_i*nSpace;
938  for(int j=0;j<nDOF_trial_element;j++)
939  {
940  //int eN_k_j=eN_k*nDOF_trial_element+j;
941  //int eN_k_j_nSpace = eN_k_j*nSpace;
942  int j_nSpace = j*nSpace;
943  int i_nSpace = i*nSpace;
944  elementJacobian_u_u[i][j] += ck.MassJacobian_weak(dm_t,u_trial_ref.data()[k*nDOF_trial_element+j],u_test_dV[i]) +
945  ck.AdvectionJacobian_weak(df,u_trial_ref.data()[k*nDOF_trial_element+j],&u_grad_test_dV[i_nSpace]) +
946  ck.DiffusionJacobian_weak(a_rowptr.data(),a_colind.data(),a,da,
947  grad_u,&u_grad_test_dV[i_nSpace],1.0,
948  u_trial_ref.data()[k*nDOF_trial_element+j],&u_grad_trial[j_nSpace]);
949  // +
950  // ck.SubgridErrorJacobian(dsubgridError_u_u[j],Lstar_u[i]) +
951  // ck.NumericalDiffusionJacobian(q_numDiff_u_last[eN_k],&u_grad_trial[j_nSpace],&u_grad_test_dV[i_nSpace]);
952  }//j
953  }//i
954  }//k
955  //
956  //load into element Jacobian into global Jacobian
957  //
958  for (int i=0;i<nDOF_test_element;i++)
959  {
960  int eN_i = eN*nDOF_test_element+i;
961  for (int j=0;j<nDOF_trial_element;j++)
962  {
963  int eN_i_j = eN_i*nDOF_trial_element+j;
964  globalJacobian.data()[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_u_u[eN_i_j]] += elementJacobian_u_u[i][j];
965  }//j
966  }//i
967  }//elements
968  //
969  //loop over exterior element boundaries to compute the surface integrals and load them into the global Jacobian
970  //
971  for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
972  {
973  register int ebN = exteriorElementBoundariesArray.data()[ebNE];
974  register int eN = elementBoundaryElementsArray.data()[ebN*2+0],
975  ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
976  eN_nDOF_trial_element = eN*nDOF_trial_element;
977  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
978  {
979  register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
980  ebNE_kb_nSpace = ebNE_kb*nSpace,
981  ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
982  ebN_local_kb_nSpace = ebN_local_kb*nSpace;
983 
984  register double u_ext=0.0,
985  grad_u_ext[nSpace],
986  m_ext=0.0,
987  dm_ext=0.0,
988  f_ext[nSpace],
989  df_ext[nSpace],
990  a_ext[nnz],
991  da_ext[nnz],
992  dflux_u_u_ext=0.0,
993  bc_u_ext=0.0,
994  //bc_grad_u_ext[nSpace],
995  bc_m_ext=0.0,
996  bc_dm_ext=0.0,
997  bc_f_ext[nSpace],
998  bc_df_ext[nSpace],
999  bc_a_ext[nnz],
1000  bc_da_ext[nnz],
1001  fluxJacobian_u_u[nDOF_trial_element],
1002  jac_ext[nSpace*nSpace],
1003  jacDet_ext,
1004  jacInv_ext[nSpace*nSpace],
1005  boundaryJac[nSpace*(nSpace-1)],
1006  metricTensor[(nSpace-1)*(nSpace-1)],
1007  metricTensorDetSqrt,
1008  dS,
1009  u_test_dS[nDOF_test_element],
1010  u_grad_trial_trace[nDOF_trial_element*nSpace],
1011  normal[3],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
1012  G[nSpace*nSpace],G_dd_G,tr_G;
1013  //
1014  //calculate the solution and gradients at quadrature points
1015  //
1016  ck.calculateMapping_elementBoundary(eN,
1017  ebN_local,
1018  kb,
1019  ebN_local_kb,
1020  mesh_dof.data(),
1021  mesh_l2g.data(),
1022  mesh_trial_trace_ref.data(),
1023  mesh_grad_trial_trace_ref.data(),
1024  boundaryJac_ref.data(),
1025  jac_ext,
1026  jacDet_ext,
1027  jacInv_ext,
1028  boundaryJac,
1029  metricTensor,
1030  metricTensorDetSqrt,
1031  normal_ref.data(),
1032  normal,
1033  x_ext,y_ext,z_ext);
1034  ck.calculateMappingVelocity_elementBoundary(eN,
1035  ebN_local,
1036  kb,
1037  ebN_local_kb,
1038  mesh_velocity_dof.data(),
1039  mesh_l2g.data(),
1040  mesh_trial_trace_ref.data(),
1041  xt_ext,yt_ext,zt_ext,
1042  normal,
1043  boundaryJac,
1044  metricTensor,
1045  integralScaling);
1046  dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref.data()[kb];
1047  //dS = metricTensorDetSqrt*dS_ref.data()[kb];
1048  ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
1049  //compute shape and solution information
1050  //shape
1051  ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,u_grad_trial_trace);
1052  //solution and gradients
1053  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);
1054  ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial_trace,grad_u_ext);
1055  //precalculate test function products with integration weights
1056  for (int j=0;j<nDOF_trial_element;j++)
1057  {
1058  u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
1059  }
1060  //
1061  //load the boundary values
1062  //
1063  bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+(1-isDOFBoundary_u.data()[ebNE_kb])*u_ext;
1064  //
1065  //calculate the internal and external trace of the pde coefficients
1066  //
1067  evaluateCoefficients(a_rowptr.data(),
1068  a_colind.data(),
1069  rho,
1070  beta,
1071  gravity.data(),
1072  alpha.data()[elementMaterialTypes.data()[eN]],
1073  n.data()[elementMaterialTypes.data()[eN]],
1074  thetaR.data()[elementMaterialTypes.data()[eN]],
1075  thetaSR.data()[elementMaterialTypes.data()[eN]],
1076  &KWs.data()[elementMaterialTypes.data()[eN]*nnz],
1077  u_ext,
1078  m_ext,
1079  dm_ext,
1080  f_ext,
1081  df_ext,
1082  a_ext,
1083  da_ext);
1084  evaluateCoefficients(a_rowptr.data(),
1085  a_colind.data(),
1086  rho,
1087  beta,
1088  gravity.data(),
1089  alpha.data()[elementMaterialTypes.data()[eN]],
1090  n.data()[elementMaterialTypes.data()[eN]],
1091  thetaR.data()[elementMaterialTypes.data()[eN]],
1092  thetaSR.data()[elementMaterialTypes.data()[eN]],
1093  &KWs.data()[elementMaterialTypes.data()[eN]*nnz],
1094  bc_u_ext,
1095  bc_m_ext,
1096  bc_dm_ext,
1097  bc_f_ext,
1098  bc_df_ext,
1099  bc_a_ext,
1100  bc_da_ext);
1101  //
1102  //calculate the flux jacobian
1103  //
1104  for (int j=0;j<nDOF_trial_element;j++)
1105  {
1106  //register int ebNE_kb_j = ebNE_kb*nDOF_trial_element+j;
1107  register int ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j;
1108  exteriorNumericalFluxJacobian(a_rowptr.data(),
1109  a_colind.data(),
1110  isDOFBoundary_u.data()[ebNE_kb],
1111  normal,
1112  a_ext,
1113  da_ext,
1114  grad_u_ext,
1115  &u_grad_trial_trace[j*nSpace],
1116  df_ext,
1117  u_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element+j],
1118  ebqe_penalty_ext.data()[ebNE_kb],//penalty,
1119  fluxJacobian_u_u[j]);
1120  }//j
1121  //
1122  //update the global Jacobian from the flux Jacobian
1123  //
1124  for (int i=0;i<nDOF_test_element;i++)
1125  {
1126  register int eN_i = eN*nDOF_test_element+i;
1127  //register int ebNE_kb_i = ebNE_kb*nDOF_test_element+i;
1128  for (int j=0;j<nDOF_trial_element;j++)
1129  {
1130  register int ebN_i_j = ebN*4*nDOF_test_X_trial_element + i*nDOF_trial_element + j;
1131  globalJacobian.data()[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] += fluxJacobian_u_u[j]*u_test_dS[i];
1132  }//j
1133  }//i
1134  }//kb
1135  }//ebNE
1136  }//computeJacobian
1137  };//Richards
1138 
1139  inline Richards_base* newRichards(int nSpaceIn,
1140  int nQuadraturePoints_elementIn,
1141  int nDOF_mesh_trial_elementIn,
1142  int nDOF_trial_elementIn,
1143  int nDOF_test_elementIn,
1144  int nQuadraturePoints_elementBoundaryIn,
1145  int CompKernelFlag)
1146  {
1147  return proteus::chooseAndAllocateDiscretization<Richards_base,Richards,CompKernel>(nSpaceIn,
1148  nQuadraturePoints_elementIn,
1149  nDOF_mesh_trial_elementIn,
1150  nDOF_trial_elementIn,
1151  nDOF_test_elementIn,
1152  nQuadraturePoints_elementBoundaryIn,
1153  CompKernelFlag);
1154  }
1155 }//proteus
1156 #endif
proteus::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:163
proteus::Richards_base
Definition: Richards.h:16
proteus::Richards::Richards
Richards()
Definition: Richards.h:36
proteus::Richards::calculateNumericalDiffusion
void calculateNumericalDiffusion(const double &shockCapturingDiffusion, const double &elementDiameter, const double &strong_residual, const double grad_u[nSpace], double &numDiff)
Definition: Richards.h:179
proteus::Richards::ck
CompKernelType ck
Definition: Richards.h:35
proteus::Richards::calculateJacobian
void calculateJacobian(arguments_dict &args)
Definition: Richards.h:729
proteus::Richards_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)=0
n
Int n
Definition: Headers.h:28
proteus::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:199
df
double df(double C, double b, double a, int q, int r)
Definition: analyticalSolutions.c:2209
ModelFactory.h
CompKernel.h
proteus::arguments_dict::scalar
T & scalar(const std::string &key)
proteus::Richards_base::calculateResidual
virtual void calculateResidual(arguments_dict &args)=0
proteus::arguments_dict::array
xt::pyarray< T > & array(const std::string &key)
num
Int num
Definition: Headers.h:32
proteus::Richards::nDOF_test_X_trial_element
const int nDOF_test_X_trial_element
Definition: Richards.h:34
v
Double v
Definition: Headers.h:95
z
Double * z
Definition: Headers.h:49
proteus::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:1139
u
Double u
Definition: Headers.h:89
xt
Definition: AddedMass.cpp:7
proteus::Richards_base::~Richards_base
virtual ~Richards_base()
Definition: Richards.h:19
proteus::Richards
Definition: Richards.h:32
proteus
Definition: ADR.h:17
proteus::Richards::calculateSubgridError_tau
void calculateSubgridError_tau(const double &elementDiameter, const double &dmt, const double dH[nSpace], double &cfl, double &tau)
Definition: Richards.h:145
proteus::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:266
proteus::f
double f(const double &g, const double &h, const double &hZ)
Definition: SW2DCV.h:58
proteus::arguments_dict
Definition: ArgumentsDict.h:70
nnz
#define nnz
Definition: Richards.h:9
proteus::Richards::calculateResidual
void calculateResidual(arguments_dict &args)
Definition: Richards.h:299
proteus::Richards::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])
Definition: Richards.h:41
Generated on Fri Jul 1 2022 11:17:40 for proteus by doxygen 1.8.20