proteus  1.8.1
C/C++/Fortran libraries
twophaseDarcyCoefficients.h
Go to the documentation of this file.
1 #ifndef TWOPHASEDARCYCOEFFICIENTS_H
2 #define TWOPHASEDARCYCOEFFICIENTS_H
3 #include "densityRelations.h"
4 #include "pskRelations.h"
5 
6 /* fully coupled, phase continuity formulation,
7  slight-compressibility assumption is used for densities
8  heterogeneity is represented by material types
9  and diffusion tensors are sparse
10 */
11 template<class PSK, class DENSITY_W, class DENSITY_N>
12 inline int twophaseDarcy_fc_sd_het_matType(int nSimplex,
13  int nPointsPerSimplex,
14  int nSpace,
15  int nParams,
16  const int* rowptr,
17  const int* colind,
18  const int* materialTypes,
19  double muw,
20  double mun,
21  const double* omega,
22  const double* Kbar, /*now has to be tensor*/
23  double b,
24  const double* rwork_psk, const int* iwork_psk,
25  const double* rwork_psk_tol,
26  const double* rwork_density_w,
27  const double* rwork_density_n,
28  const double* g,
29  const double* x,
30  const double* sw,
31  const double* psiw,
32  double* mw,
33  double* dmw_dsw,
34  double* dmw_dpsiw,
35  double* mn,
36  double* dmn_dsw,
37  double* dmn_dpsiw,
38  double* psin,
39  double* dpsin_dsw,
40  double* dpsin_dpsiw,
41  double* phi_psiw,
42  double* dphi_psiw_dpsiw,
43  double* phi_psin,
44  double* dphi_psin_dpsiw,
45  double* dphi_psin_dsw,
46  double* aw,
47  double* daw_dsw,
48  double* daw_dpsiw,
49  double* an,
50  double* dan_dsw,
51  double* dan_dpsiw)
52 {
53  int matID;
54  PSK psk(rwork_psk,iwork_psk); psk.setTolerances(rwork_psk_tol);
55  DENSITY_W density_w(rwork_density_w);
56  DENSITY_N density_n(rwork_density_n);
57  double drhow_dsw,drhon_dsw,drhow_dpsiw,drhon_dpsiw;
58 
59  const int nnz = rowptr[nSpace];
60  for (int eN=0;eN<nSimplex;eN++)
61  {
62  matID = materialTypes[eN];
63  psk.setParams(&rwork_psk[matID*nParams]);
64  for(int pN=0,i;pN<nPointsPerSimplex;pN++)
65  {
66  i = eN*nPointsPerSimplex+pN;
67  psk.calc(sw[i]);
68  /*non-wetting phase pressure head */
69  psin[i] = psiw[i] + psk.psic;
70  dpsin_dsw[i] = psk.dpsic;
71  dpsin_dpsiw[i]= 1.0;
72 
73  /*mwf need to make sure normalized by rhow and rhon respectively!*/
74  density_w.calc(psiw[i]);
75  density_n.calc(psin[i]);
76  /*density of wetting phase just a function of psiw*/
77  drhow_dpsiw = density_w.drho;
78  drhow_dsw = 0.0;
79  /*density of nonwetting phase a function of sw and psiw through psin*/
80  drhon_dpsiw = density_n.drho;
81  drhon_dsw = density_n.drho*psk.dpsic;
82 
83  /* w-phase mass */
84  mw[i] = omega[matID]*density_w.rho*sw[i];
85  dmw_dsw[i] = omega[matID]*(density_w.rho + sw[i]*drhow_dsw);
86  dmw_dpsiw[i]= omega[matID]*sw[i]*drhow_dpsiw;
87 
88  /* n-phase mass */
89  mn[i] = omega[matID]*density_n.rho*(1.0-sw[i]);
90  dmn_dsw[i] = omega[matID]*((1.0-sw[i])*drhon_dsw -density_n.rho);
91  dmn_dpsiw[i]= omega[matID]*(1.0-sw[i])*drhon_dpsiw;
92 
93  /* w-phase potential */
94  phi_psiw[i] = psiw[i];
95  dphi_psiw_dpsiw[i] = 1.0; /*include density dependence below*/
96 
97  /* n-phase potential */
98  phi_psin[i] = psiw[i] + psk.psic;
99  dphi_psin_dpsiw[i]= 1.0;
100  dphi_psin_dsw[i]= psk.dpsic;
101 
102  for (int I=0;I<nSpace;I++)
103  {
104  /* update potentials with gravity */
105  /*slight compressibility, normalized density=1.0*/
106  phi_psiw[i] -= g[I]*x[i*3+I];
107 
108  phi_psin[i] -= b*g[I]*x[i*3+I];
109 
110  for (int m=rowptr[I]; m < rowptr[I+1]; m++)
111  {
112  /* w-phase diffusion */
113 
114  aw[i*nnz+m] = density_w.rho*Kbar[matID*nnz+m]*psk.krw/muw;
115  daw_dsw[i*nnz+m] = density_w.rho*Kbar[matID*nnz+m]*psk.dkrw/muw + drhow_dsw*Kbar[matID*nnz+m]*psk.krw/muw;
116  daw_dpsiw[i*nnz+m]= drhow_dpsiw*Kbar[matID*nnz+m]*psk.krw/muw;
117 
118  /* n-phase diffusion */
119  an[i*nnz+m] = density_n.rho*Kbar[matID*nnz+m]*psk.krn/mun;
120  dan_dsw[i*nnz+m] = density_n.rho*Kbar[matID*nnz+m]*psk.dkrn/mun + drhon_dsw*Kbar[matID*nnz+m]*psk.krn/mun;
121  dan_dpsiw[i*nnz+m] = drhon_dpsiw*Kbar[matID*nnz+m]*psk.krn/mun;
122  }
123  }
124  }
125  }
126  return 0;
127 }
128 template<class PSK, class DENSITY_W, class DENSITY_N>
129 inline int twophaseDarcy_fc_sd_het_matType_nonPotentialForm(int compressibilityFlag,
130  int nSimplex,
131  int nPointsPerSimplex,
132  int nSpace,
133  int nParams,
134  const int* rowptr,
135  const int* colind,
136  const int* materialTypes,
137  double muw,
138  double mun,
139  const double* omega,
140  const double* Kbar, /*now has to be tensor*/
141  double b,
142  const double* rwork_psk, const int* iwork_psk,
143  const double* rwork_psk_tol,
144  const double* rwork_density_w,
145  const double* rwork_density_n,
146  const double* g,
147  const double* x,
148  const double* sw,
149  const double* psiw,
150  double* mw,
151  double* dmw_dsw,
152  double* dmw_dpsiw,
153  double* mn,
154  double* dmn_dsw,
155  double* dmn_dpsiw,
156  double* psin,
157  double* dpsin_dsw,
158  double* dpsin_dpsiw,
159  double* phi_psiw,
160  double* dphi_psiw_dpsiw,
161  double* phi_psin,
162  double* dphi_psin_dpsiw,
163  double* dphi_psin_dsw,
164  double* fw,
165  double* dfw_dsw,
166  double* dfw_dpsiw,
167  double* fn,
168  double* dfn_dsw,
169  double* dfn_dpsiw,
170  double* aw,
171  double* daw_dsw,
172  double* daw_dpsiw,
173  double* an,
174  double* dan_dsw,
175  double* dan_dpsiw)
176 {
177  int matID;
178  PSK psk(rwork_psk,iwork_psk); psk.setTolerances(rwork_psk_tol);
179  DENSITY_W density_w(rwork_density_w);
180  DENSITY_N density_n(rwork_density_n);
181  double drhow_dsw,drhon_dsw,drhow_dpsiw,drhon_dpsiw;
182  double rhow_x,drhow_x_dsw,drhow_x_dpsiw,rhon_x,drhon_x_dsw,drhon_x_dpsiw;
183 
184  const int nnz = rowptr[nSpace];
185  for (int eN=0;eN<nSimplex;eN++)
186  {
187  matID = materialTypes[eN];
188  psk.setParams(&rwork_psk[matID*nParams]);
189  for(int pN=0,i;pN<nPointsPerSimplex;pN++)
190  {
191  i = eN*nPointsPerSimplex+pN;
192  psk.calc(sw[i]);
193  /*non-wetting phase pressure head */
194  psin[i] = psiw[i] + psk.psic;
195  dpsin_dsw[i] = psk.dpsic;
196  dpsin_dpsiw[i]= 1.0;
197 
198  /*mwf need to make sure normalized by rhow and rhon respectively!*/
199  density_w.calc(psiw[i]);
200  density_n.calc(psin[i]);
201  /*density of wetting phase just a function of psiw*/
202  drhow_dpsiw = density_w.drho;
203  drhow_dsw = 0.0;
204  /*density of nonwetting phase a function of sw and psiw through psin*/
205  drhon_dpsiw = density_n.drho;
206  drhon_dsw = density_n.drho*psk.dpsic;
207 
208  /* w-phase mass */
209  mw[i] = omega[matID]*density_w.rho*sw[i];
210  dmw_dsw[i] = omega[matID]*(density_w.rho + sw[i]*drhow_dsw);
211  dmw_dpsiw[i]= omega[matID]*sw[i]*drhow_dpsiw;
212 
213  /* n-phase mass */
214  mn[i] = omega[matID]*density_n.rho*(1.0-sw[i]);
215  dmn_dsw[i] = omega[matID]*((1.0-sw[i])*drhon_dsw -density_n.rho);
216  dmn_dpsiw[i]= omega[matID]*(1.0-sw[i])*drhon_dpsiw;
217 
218  /* w-phase potential */
219  phi_psiw[i] = psiw[i];
220  dphi_psiw_dpsiw[i] = 1.0;
221 
222  /* n-phase potential */
223  phi_psin[i] = psiw[i] + psk.psic;
224  dphi_psin_dpsiw[i]= 1.0;
225  dphi_psin_dsw[i]= psk.dpsic;
226 
227  rhow_x = density_w.rho;
228  drhow_x_dsw = drhow_dsw;
229  drhow_x_dpsiw= drhow_dpsiw;
230  rhon_x = density_n.rho;
231  drhon_x_dsw = drhon_dsw;
232  drhon_x_dpsiw= drhon_dpsiw;
233  if (compressibilityFlag == 0) /*slight compressibility*/
234  {
235  rhow_x = 1.0;
236  drhow_x_dsw = 0.0;
237  drhow_x_dpsiw= 0.0;
238  rhon_x = 1.0;
239  drhon_x_dsw = 0.0;
240  drhon_x_dpsiw= 0.0;
241  }
242 
243  for (int I=0;I<nSpace;I++)
244  {
245  fw[i*nSpace+I] = 0.0;
246  dfw_dsw[i*nSpace+I] = 0.0;
247  dfw_dpsiw[i*nSpace+I]= 0.0;
248 
249  fn[i*nSpace+I] = 0.0;
250  dfn_dsw[i*nSpace+I] = 0.0;
251  dfn_dpsiw[i*nSpace+I]= 0.0;
252 
253  for (int m=rowptr[I]; m < rowptr[I+1]; m++)
254  {
255  /*assuming slight compressibility*/
256  /* w-phase advection (gravity) */
257  fw[i*nSpace+I] += rhow_x*rhow_x*Kbar[matID*nnz+m]*psk.krw/muw*g[colind[m]];
258  dfw_dsw[i*nSpace+I] += rhow_x*rhow_x*Kbar[matID*nnz+m]*psk.dkrw/muw*g[colind[m]] + 2.0*rhow_x*drhow_x_dsw*Kbar[matID*nnz+m]*psk.krw/muw*g[colind[m]];
259  dfw_dpsiw[i*nSpace+I] += 2.0*rhow_x*drhow_x_dpsiw*Kbar[matID*nnz+m]*psk.krw/muw*g[colind[m]];
260  /* n-phase advection (gravity) */
261  fn[i*nSpace+I] += rhon_x*rhon_x*Kbar[matID*nnz+m]*psk.krn/mun*b*g[colind[m]];
262  dfn_dsw[i*nSpace+I] += rhon_x*rhon_x*Kbar[matID*nnz+m]*psk.dkrn/mun*b*g[colind[m]] + 2.0*rhon_x*drhon_x_dsw*Kbar[matID*nnz+m]*psk.krn/mun*b*g[colind[m]];
263  dfn_dpsiw[i*nSpace+I] += 2.0*rhon_x*drhon_x_dpsiw*Kbar[matID*nnz+m]*psk.krn/mun*b*g[colind[m]];
264 
265  /* w-phase diffusion */
266  aw[i*nnz+m] = rhow_x*Kbar[matID*nnz+m]*psk.krw/muw;
267  daw_dsw[i*nnz+m] = rhow_x*Kbar[matID*nnz+m]*psk.dkrw/muw + drhow_x_dsw*Kbar[matID*nnz+m]*psk.krw/muw;
268  daw_dpsiw[i*nnz+m]= drhow_x_dpsiw*Kbar[matID*nnz+m]*psk.krw/muw;
269 
270  /* n-phase diffusion */
271  an[i*nnz+m] = rhon_x*Kbar[matID*nnz+m]*psk.krn/mun;
272  dan_dsw[i*nnz+m] = rhon_x*Kbar[matID*nnz+m]*psk.dkrn/mun + drhon_x_dsw*Kbar[matID*nnz+m]*psk.krn/mun;
273  dan_dpsiw[i*nnz+m] = drhon_x_dpsiw*Kbar[matID*nnz+m]*psk.krn/mun;
274 
275 //mwf original
276 // /* w-phase diffusion */
277 
278 // aw[i*nnz+m] = density_w.rho*Kbar[matID*nnz+m]*psk.krw/muw;
279 // daw_dsw[i*nnz+m] = density_w.rho*Kbar[matID*nnz+m]*psk.dkrw/muw + drhow_dsw*Kbar[matID*nnz+m]*psk.krw/muw;
280 // daw_dpsiw[i*nnz+m]= drhow_dpsiw*Kbar[matID*nnz+m]*psk.krw/muw;
281 
282 // /* n-phase diffusion */
283 // an[i*nnz+m] = density_n.rho*Kbar[matID*nnz+m]*psk.krn/mun;
284 // dan_dsw[i*nnz+m] = density_n.rho*Kbar[matID*nnz+m]*psk.dkrn/mun + drhon_dsw*Kbar[matID*nnz+m]*psk.krn/mun;
285 // dan_dpsiw[i*nnz+m] = drhon_dpsiw*Kbar[matID*nnz+m]*psk.krn/mun;
286  }
287  }
288  }
289  }
290  return 0;
291 }
292 
293 inline int twophaseDarcy_vol_frac(int nSimplex,
294  int nPointsPerSimplex,
295  const int* materialTypes,
296  const double* omega,
297  const double* sw,
298  double * vol_frac_w,
299  double * vol_frac_n)
300 {
301  int matID;
302  for (int eN=0; eN < nSimplex; eN++)
303  {
304  matID = materialTypes[eN];
305  for (int pN=0,i;pN<nPointsPerSimplex;pN++)
306  {
307  i = eN*nPointsPerSimplex+pN;
308  vol_frac_w[i] = omega[matID]*sw[i];
309  vol_frac_n[i] = omega[matID]*(1.0-sw[i]);
310 
311  }
312  }
313  return 0;
314 }
315 
316 /* fully coupled, phase continuity formulation, wetting phase pressure and capillary pressure are primary dependent
317  variables except say u_1 = 1-S for u_1 <= 0
318  slight-compressibility assumption is used for densities
319  heterogeneity is represented by material types
320  and diffusion tensors are sparse
321 */
322 template<class PSK, class DENSITY_W, class DENSITY_N>
323 inline int twophaseDarcy_fc_pp_sd_het_matType(int nSimplex,
324  int nPointsPerSimplex,
325  int nSpace,
326  int nParams,
327  const int* rowptr,
328  const int* colind,
329  const int* materialTypes,
330  double muw,
331  double mun,
332  const double* omega,
333  const double* Kbar, /*now has to be tensor*/
334  double b,
335  const double* rwork_psk, const int* iwork_psk,
336  const double* rwork_psk_tol,
337  const double* rwork_density_w,
338  const double* rwork_density_n,
339  const double* g,
340  const double* x,
341  const double* psiw,
342  const double* psic,
343  double *sw,
344  double* mw,
345  double* dmw_dpsiw,
346  double* dmw_dpsic,
347  double* mn,
348  double* dmn_dpsiw,
349  double* dmn_dpsic,
350  double* phi_psiw,
351  double* dphi_psiw_dpsiw,
352  double* phi_psin,
353  double* dphi_psin_dpsiw,
354  double* dphi_psin_dpsic,
355  double* aw,
356  double* daw_dpsiw,
357  double* daw_dpsic,
358  double* an,
359  double* dan_dpsiw,
360  double* dan_dpsic)
361 {
362  int matID;
363  PSK psk(rwork_psk,iwork_psk); psk.setTolerances(rwork_psk_tol);
364  DENSITY_W density_w(rwork_density_w);
365  DENSITY_N density_n(rwork_density_n);
366  double drhow_dpsiw,drhow_dpsic,drhon_dpsiw,drhon_dpsic,
367  swi,psin,dsw_dpsic,dpsin_dpsic;
368  const int nnz = rowptr[nSpace];
369  for (int eN=0;eN<nSimplex;eN++)
370  {
371  matID = materialTypes[eN];
372  psk.setParams(&rwork_psk[matID*nParams]);
373  for(int pN=0,i;pN<nPointsPerSimplex;pN++)
374  {
375  i = eN*nPointsPerSimplex+pN;
376  psk.calc_from_psic(psic[i]);
377 
378  /*non-wetting phase pressure head */
379  psin = psiw[i] + psic[i];
380  dpsin_dpsic = 1.0;
381  if (psic[i] <= 0.0)
382  {
383  psin = psiw[i];
384  dpsin_dpsic = 0.0;
385  }
386  /*aqueous phase saturation*/
387  swi = psk.Se*(psk.Sw_max-psk.Sw_min) + psk.Sw_min;
388  sw[i] = swi;
389  dsw_dpsic = psk.dSe_dpsic/psk.dSe_dSw;
390 
391  /*mwf need to make sure normalized by rhow and rhon respectively!*/
392  density_w.calc(psiw[i]);
393  density_n.calc(psin);
394  /*density of wetting phase just a function of psiw*/
395  drhow_dpsiw = density_w.drho;
396  drhow_dpsic = 0.0;
397  /*density of nonwetting phase a function of psic and psiw through psin*/
398  drhon_dpsiw = density_n.drho;
399  drhon_dpsic = density_n.drho;
400  if (psic[i] <= 0.0)
401  drhon_dpsic = 0.0;
402  /* w-phase mass */
403  mw[i] = omega[matID]*density_w.rho*swi;
404  dmw_dpsic[i]= omega[matID]*(density_w.rho*dsw_dpsic + swi*drhow_dpsic);
405  dmw_dpsiw[i]= omega[matID]*swi*drhow_dpsiw;
406 
407  /* n-phase mass */
408  mn[i] = omega[matID]*density_n.rho*(1.0-swi);
409  dmn_dpsic[i]= omega[matID]*((1.0-swi)*drhon_dpsic -density_n.rho*dsw_dpsic);
410  dmn_dpsiw[i]= omega[matID]*(1.0-swi)*drhon_dpsiw;
411 
412  /* w-phase potential */
413  phi_psiw[i] = psiw[i];
414  dphi_psiw_dpsiw[i] = 1.0; /*include density dependence below*/
415 
416  /* n-phase potential */
417  phi_psin[i] = psin;
418  dphi_psin_dpsiw[i]= 1.0;
419  dphi_psin_dpsic[i]= dpsin_dpsic;
420 
421  for (int I=0;I<nSpace;I++)
422  {
423  /* update potentials with gravity */
424  /*slight compressibility, normalized density=1.0*/
425  phi_psiw[i] -= g[I]*x[i*3+I];
426 
427  phi_psin[i] -= b*g[I]*x[i*3+I];
428 
429  for (int m=rowptr[I]; m < rowptr[I+1]; m++)
430  {
431  /* w-phase diffusion */
432 
433  aw[i*nnz+m] = density_w.rho*Kbar[matID*nnz+m]*psk.krw/muw;
434  daw_dpsic[i*nnz+m]= density_w.rho*Kbar[matID*nnz+m]*psk.dkrw/muw + drhow_dpsic*Kbar[matID*nnz+m]*psk.krw/muw;
435  daw_dpsiw[i*nnz+m]= drhow_dpsiw*Kbar[matID*nnz+m]*psk.krw/muw;
436 
437  /* n-phase diffusion */
438  an[i*nnz+m] = density_n.rho*Kbar[matID*nnz+m]*psk.krn/mun;
439  dan_dpsic[i*nnz+m] = density_n.rho*Kbar[matID*nnz+m]*psk.dkrn/mun + drhon_dpsic*Kbar[matID*nnz+m]*psk.krn/mun;
440  dan_dpsiw[i*nnz+m] = drhon_dpsiw*Kbar[matID*nnz+m]*psk.krn/mun;
441  }
442  }
443  /*adjust mass terms if psic <= 0, */
444  if (psic[i] <= 0.0)/*0.0)*/
445  {
446  /*mwf debug
447  printf("psic[%d]=%g krn=%g dkrn=%g krw=%g dkrw=%g sw=%g dsw_dpsic=%g drhow_dpsic=%g drhon_dpsic=%g \n",i,psic[i],psk.krn,psk.dkrn,psk.krw,psk.dkrw,swi,
448  dsw_dpsic,drhow_dpsic,drhon_dpsic);
449  */
450 
451  dmn_dpsic[i] = omega[matID]*density_n.rho;
452  dmw_dpsic[i] =-omega[matID]*density_n.rho;//0.0;
453  //dmn_dpsiw[i] = 0.0;
454  phi_psin[i] = 0.0;
455  dphi_psin_dpsiw[i]= 0.0;
456  dphi_psin_dpsic[i]= 0.0;
457  for (int I=0;I<nSpace;I++)
458  {
459  for (int m=rowptr[I]; m < rowptr[I+1]; m++)
460  {
461  an[i*nnz+m] = 0.0;
462  dan_dpsic[i*nnz+m] = 0.0;
463  dan_dpsiw[i*nnz+m] = 0.0;
464  daw_dpsic[i*nnz+m] = 0.0;
465  }
466  }
467 
468  }
469  }
470  }
471  return 0;
472 }
473 
474 /* saturation equation in fractional flow formulation for incompressible flow
475  heterogeneity is represented by materialTypes
476  */
477 template<class PSK>
478 static inline int twophaseDarcy_incompressible_split_sd_saturation_het_matType(int nSimplex,
479  int nPointsPerSimplex,
480  int nSpace,
481  int nParams,
482  const int* rowptr,
483  const int* colind,
484  const int* materialTypes,
485  double muw,
486  double mun,
487  const double* omega,
488  const double* Kbar,
489  double b,
490  double capillaryDiffusionScaling,
491  double advectionScaling,
492  const double* rwork_psk, const int* iwork_psk,
493  const double* rwork_psk_tol,
494  const double* rwork_density_w,
495  const double* rwork_density_n,
496  const double* g,
497  const double* qt,
498  const double* sw,
499  double* m,
500  double* dm,
501  double* phi,
502  double* dphi,
503  double* f,
504  double* df,
505  double* a,
506  double* da)
507 {
508  int matID;
509  PSK psk(rwork_psk,iwork_psk); psk.setTolerances(rwork_psk_tol);
510  FractionalFlowVariables fracFlow(muw,mun);
511  ConstantDensity density_w(rwork_density_w),density_n(rwork_density_n);
512  const int nnz = rowptr[nSpace];
513  //mwf debug
514  //std::cout<<"entering twpffinc rhow= "<<density_w.rho<<" rhon= "<<density_n.rho<<" b= "<<b<<" ";
515  //for (int I=0; I < nSpace; I++)
516  // {
517  // std::cout<<"g["<<I<<"]= "<<g[I]<<" ";
518  // }
519  for (int eN=0;eN<nSimplex;eN++)
520  {
521  matID = materialTypes[eN];
522  psk.setParams(&rwork_psk[matID*nParams]);
523  for(int pN=0,i;pN<nPointsPerSimplex;pN++)
524  {
525  i = eN*nPointsPerSimplex+pN;
526  psk.calc(sw[i]);
527  fracFlow.calc(psk,density_w,density_n);
528 
529  /* wetting phase mass */
530  m[i] = omega[matID]*density_w.rho*sw[i];
531  dm[i] = omega[matID]*density_w.rho;
532 
533  /* capillary potential */
534  phi[i] = psk.psic;
535  dphi[i]= psk.dpsic;
536 
537  for (int I=0;I<nSpace;I++)
538  {
539  /* wetting phase advection */
540  /* todo, remove diagonal assumption on K*/
541  for (int m=rowptr[I]; m < rowptr[I+1]; m++)
542  {
543  const int J = colind[m];
544  if (I==J)
545  {
546  f[i*nSpace+I] = advectionScaling*(qt[i*nSpace+I]*fracFlow.fw
547  - Kbar[matID*nnz+m]*fracFlow.lambdaw*fracFlow.fn*(b*density_n.rho-density_w.rho)*g[I]) ;
548  df[i*nSpace+I] = advectionScaling*(qt[i*nSpace+I]*fracFlow.dfw
549  - (Kbar[matID*nnz+m]*g[I]*(b*density_n.rho-density_w.rho))*(fracFlow.lambdaw*fracFlow.dfn + fracFlow.fn*fracFlow.dlambdaw));
550  }
551  /* wetting phase capillary diffusion */
552  /*include scaling factor in case want to turn of capillary diffusion to test hyperbolic approximations*/
553  a[i*nnz+m] = -capillaryDiffusionScaling*Kbar[matID*nnz+m]*fracFlow.lambdaw*fracFlow.fn;
554  da[i*nnz+m] = -capillaryDiffusionScaling*Kbar[matID*nnz+m]*(fracFlow.dlambdaw*fracFlow.fn + fracFlow.lambdaw*fracFlow.dfn);
555  }
556  }
557  }
558  }
559  return 0;
560 }
561 
562 /* saturation equation in fractional flow formulation for slight compressible flow
563  heterogeneity is represented by materialTypes
564  */
565 template<class PSK, class DENSITY_W>
566 static inline int twophaseDarcy_slightCompressible_split_sd_saturation_het_matType(int nSimplex,
567  int nPointsPerSimplex,
568  int nSpace,
569  int nParams,
570  const int* rowptr,
571  const int* colind,
572  const int* materialTypes,
573  double muw,
574  double mun,
575  const double* omega,
576  const double* Kbar,
577  double b,
578  double capillaryDiffusionScaling,
579  double advectionScaling,
580  const double* rwork_psk, const int* iwork_psk,
581  const double* rwork_psk_tol,
582  const double* rwork_density_w,
583  const double* rwork_density_n,
584  const double* g,
585  const double* qt,
586  const double* psiw,
587  const double* sw,
588  double* m,
589  double* dm,
590  double* phi,
591  double* dphi,
592  double* f,
593  double* df,
594  double* a,
595  double* da)
596 {
597  int matID;
598  PSK psk(rwork_psk,iwork_psk); psk.setTolerances(rwork_psk_tol);
599  FractionalFlowVariables fracFlow(muw,mun);
600  /*normalized densities \rho_{\alpha} = \varrho_{\alpha}/\varrho_{\alpha,0}
601  for spatial term, assuming slight compressiblity so assume \rho_{\alpha} = 1
602  */
603  const double rwork_density_w_x[1] = {1.0}; const double rwork_density_n_x[1] = {1.0};
604  ConstantDensity density_w_x(rwork_density_w_x),density_n_x(rwork_density_n_x);
605  /*for accumulation term*/
606  DENSITY_W density_w(rwork_density_w);
607  const int nnz = rowptr[nSpace];
608  //mwf debug
609  //std::cout<<"entering twpffsl rhow= "<<density_w.rho<<" rhon= "<<density_n.rho<<" b= "<<b<<" ";
610  //for (int I=0; I < nSpace; I++)
611  // {
612  // std::cout<<"g["<<I<<"]= "<<g[I]<<" ";
613  // }
614  for (int eN=0;eN<nSimplex;eN++)
615  {
616  matID = materialTypes[eN];
617  psk.setParams(&rwork_psk[matID*nParams]);
618  for(int pN=0,i;pN<nPointsPerSimplex;pN++)
619  {
620  i = eN*nPointsPerSimplex+pN;
621  psk.calc(sw[i]);
622  density_w.calc(psiw[i]);
623  fracFlow.calc(psk,density_w_x,density_n_x);
624 
625 
626  /* wetting phase mass */
627  m[i] = omega[matID]*density_w.rho*sw[i];
628  /*density of wetting phase just a function of psiw*/
629  dm[i] = omega[matID]*density_w.rho;
630 
631  /* capillary potential */
632  phi[i] = psk.psic;
633  dphi[i]= psk.dpsic;
634 
635  for (int I=0;I<nSpace;I++)
636  {
637  /* wetting phase advection */
638  /* todo remove diagonal K assumption*/
639  for (int m=rowptr[I]; m < rowptr[I+1]; m++)
640  {
641  const int J = colind[m];
642  if (I==J)
643  {
644  f[i*nSpace+I] = advectionScaling*(qt[i*nSpace+I]*fracFlow.fw
645  - Kbar[matID*nnz+m]*fracFlow.lambdaw*fracFlow.fn*(b*density_n_x.rho-density_w_x.rho)*g[I]) ;
646  df[i*nSpace+I] = advectionScaling*(qt[i*nSpace+I]*fracFlow.dfw
647  - (Kbar[matID*nnz+m]*g[I]*(b*density_n_x.rho-density_w_x.rho))*(fracFlow.lambdaw*fracFlow.dfn + fracFlow.fn*fracFlow.dlambdaw));
648  }
649  /* wetting phase capillary diffusion */
650  /*include scaling factor in case want to turn of capillary diffusion to test hyperbolic approximations*/
651  a[i*nnz+m] = -capillaryDiffusionScaling*Kbar[matID*nnz+m]*fracFlow.lambdaw*fracFlow.fn;
652  da[i*nnz+m] = -capillaryDiffusionScaling*Kbar[matID*nnz+m]*(fracFlow.dlambdaw*fracFlow.fn + fracFlow.lambdaw*fracFlow.dfn);
653  }
654  }
655  }
656  }
657  return 0;
658 }
659 
660 /* pressure equation in fractional flow formulation for incompressible flow
661  heterogeneity is represented by materialTypes
662 */
663 template<class PSK>
664 static inline int twophaseDarcy_incompressible_split_sd_pressure_het_matType(int nSimplex,
665  int nPointsPerSimplex,
666  int nSpace,
667  int nParams,
668  const int* rowptr,
669  const int* colind,
670  const int* materialTypes,
671  double muw,
672  double mun,
673  const double* omega,
674  const double* Kbar,
675  double b,
676  double capillaryDiffusionScaling,
677  const double* rwork_psk, const int* iwork_psk,
678  const double* rwork_psk_tol,
679  const double* rwork_density_w,
680  const double* rwork_density_n,
681  const double* g,
682  const double* sw,
683  const double* grad_psic,
684  double* f,
685  double* a)
686 {
687  int matID;
688  PSK psk(rwork_psk,iwork_psk); psk.setTolerances(rwork_psk_tol);
689  FractionalFlowVariables fracFlow(muw,mun);
690  ConstantDensity density_w(rwork_density_w),density_n(rwork_density_n);
691  const int nnz = rowptr[nSpace];
692  for (int eN=0;eN<nSimplex;eN++)
693  {
694  matID = materialTypes[eN];
695  psk.setParams(&rwork_psk[matID*nParams]);
696  for(int pN=0,i;pN<nPointsPerSimplex;pN++)
697  {
698  i = eN*nPointsPerSimplex+pN;
699  psk.calc(sw[i]);
700  fracFlow.calc(psk,density_w,density_n);
701  /*todo remove diagonal assumption on f*/
702  for (int I=0;I<nSpace;I++)
703  {
704  /*include scaling factor in case want to turn of capillary diffusion to test hyperbolic approximations*/
705  for (int m=rowptr[I]; m < rowptr[I+1]; m++)
706  {
707  const int J = colind[m];
708  if (I == J)
709  {
710  //mwf original
711  f[i*nSpace+I] = - capillaryDiffusionScaling*Kbar[matID*nnz+m]*fracFlow.lambdat*(fracFlow.fn*grad_psic[i*nSpace+I]) +
712  Kbar[matID*nnz+m]*fracFlow.lambdat*(density_w.rho + fracFlow.fn*(b*density_n.rho-density_w.rho))*g[I];
713  //f[i*nSpace+I] = - capillaryDiffusionScaling*Kbar[matID*nnz+m]*fracFlow.lambdat*(fracFlow.fn*psk.dpsic*grad_sw[i*nSpace+I]) +
714  // Kbar[matID*nnz+m]*fracFlow.lambdat*(density_w.rho + fracFlow.fn*(b*density_n.rho-density_w.rho))*g[I];
715  }
716  a[i*nnz+m] = Kbar[matID*nnz+m]*fracFlow.lambdat;
717  }
718  }
719  }
720  }
721  return 0.0;
722 }
723 
724 /* pressure equation in fractional flow formulation for incompressible flow
725  heterogeneity is represented by materialTypes
726 */
727 template<class PSK, class DENSITY_W, class DENSITY_N>
728 static inline int twophaseDarcy_slightCompressible_split_sd_pressure_het_matType(int nSimplex,
729  int nPointsPerSimplex,
730  int nSpace,
731  int nParams,
732  const int* rowptr,
733  const int* colind,
734  const int* materialTypes,
735  double muw,
736  double mun,
737  const double* omega,
738  const double* Kbar,
739  double b,
740  double capillaryDiffusionScaling,
741  const double* rwork_psk, const int* iwork_psk,
742  const double* rwork_psk_tol,
743  const double* rwork_density_w,
744  const double* rwork_density_n,
745  const double* g,
746  const double* sw,
747  const double* psiw,
748  const double* psin,
749  const double* grad_psic,
750  double* m,
751  double* dm,
752  double* f,
753  double* a)
754 {
755  int matID;
756  PSK psk(rwork_psk,iwork_psk); psk.setTolerances(rwork_psk_tol);
757  FractionalFlowVariables fracFlow(muw,mun);
758  /*normalized densities \rho_{\alpha} = \varrho_{\alpha}/\varrho_{\alpha,0}
759  for spatial term, assuming slight compressiblity so assume \rho_{\alpha} = 1
760  */
761  const double rwork_density_w_x[1] = {1.0}; const double rwork_density_n_x[1] = {1.0};
762  ConstantDensity density_w_x(rwork_density_w_x),density_n_x(rwork_density_n_x);
763  DENSITY_W density_w(rwork_density_w);
764  DENSITY_N density_n(rwork_density_n);
765  const int nnz = rowptr[nSpace];
766  double drhow_dpsiw,drhon_dpsiw;
767  for (int eN=0;eN<nSimplex;eN++)
768  {
769  matID = materialTypes[eN];
770  psk.setParams(&rwork_psk[matID*nParams]);
771  for(int pN=0,i;pN<nPointsPerSimplex;pN++)
772  {
773  i = eN*nPointsPerSimplex+pN;
774  psk.calc(sw[i]);
775  density_w.calc(psiw[i]);
776  density_n.calc(psin[i]);
777 
778  fracFlow.calc(psk,density_w_x,density_n_x);
779 
780  /*density of wetting phase just a function of psiw*/
781  drhow_dpsiw = density_w.drho;
782  /*density of nonwetting phase a function of sw and psiw through psin but ignore sw derivative because of splitting*/
783  drhon_dpsiw = density_n.drho;
784  /*mwf debug
785  printf("twpff slc psiw=%g psin=%g rhow=%g rhon=%g drhow=%g drhon=%g \n",psiw[i],psin[i],density_w.rho,density_n.rho,drhow_dpsiw,drhon_dpsiw);
786  */
787  m[i] = omega[matID]*density_w.rho*sw[i] + omega[matID]*density_n.rho*(1.0-sw[i]);
788  dm[i]= omega[matID]*drhow_dpsiw*sw[i] + omega[matID]*drhon_dpsiw*(1.0-sw[i]);
789 
790  for (int I=0;I<nSpace;I++)
791  {
792  /*include scaling factor in case want to turn of capillary diffusion to test hyperbolic approximations*/
793  for (int m=rowptr[I]; m < rowptr[I+1]; m++)
794  {
795  const int J = colind[m];
796  if (I == J)
797  {
798  f[i*nSpace+I] = - capillaryDiffusionScaling*Kbar[matID*nnz+m]*fracFlow.lambdat*(fracFlow.fn*grad_psic[i*nSpace+I]) +
799  Kbar[matID*nnz+m]*fracFlow.lambdat*(density_w_x.rho + fracFlow.fn*(b*density_n_x.rho-density_w_x.rho))*g[I];
800  }
801  a[i*nnz+m] = Kbar[matID*nnz+m]*fracFlow.lambdat;
802  }
803  }
804  }
805  }
806  return 0.0;
807 }
808 
809 /* saturation equation in fractional flow formulation for slight compressible flow
810  heterogeneity is represented by materialTypes
811  */
812 template<class PSK, class DENSITY_N>
813 static inline int twophaseDarcy_compressibleN_split_sd_saturation_het_matType(int nSimplex,
814  int nPointsPerSimplex,
815  int nSpace,
816  int nParams,
817  const int* rowptr,
818  const int* colind,
819  const int* materialTypes,
820  double muw,
821  double mun,
822  const double* omega,
823  const double* Kbar,
824  double b,
825  double capillaryDiffusionScaling,
826  double advectionScaling,
827  const double* rwork_psk, const int* iwork_psk,
828  const double* rwork_psk_tol,
829  const double* rwork_density_w,
830  const double* rwork_density_n,
831  const double* g,
832  const double* qt,
833  const double* psiw,
834  const double* sw,
835  double* m,
836  double* dm,
837  double* phi,
838  double* dphi,
839  double* f,
840  double* df,
841  double* a,
842  double* da)
843 {
844  int matID;
845  PSK psk(rwork_psk,iwork_psk); psk.setTolerances(rwork_psk_tol);
846  CompressibleN_FractionalFlowVariables fracFlow(muw,mun);
847  ConstantDensity density_w(rwork_density_w);
848  DENSITY_N density_n(rwork_density_n);
849  const int nnz = rowptr[nSpace];
850  //mwf debug
851  //std::cout<<"entering twpffc rhow= "<<density_w.rho<<" rhon= "<<density_n.rho<<" b= "<<b<<" ";
852  //for (int I=0; I < nSpace; I++)
853  // {
854  // std::cout<<"g["<<I<<"]= "<<g[I]<<" ";
855  // }
856  double psin;
857  for (int eN=0;eN<nSimplex;eN++)
858  {
859  matID = materialTypes[eN];
860  psk.setParams(&rwork_psk[matID*nParams]);
861  for(int pN=0,i;pN<nPointsPerSimplex;pN++)
862  {
863  i = eN*nPointsPerSimplex+pN;
864  psk.calc(sw[i]);
865 
866  psin = psiw[i] + psk.psic;
867 
868  density_w.calc(psiw[i]);
869  density_n.calc(psin);
870 
871  fracFlow.calc(psk,density_w,density_n);
872 
873 
874  /* wetting phase mass */
875  m[i] = omega[matID]*density_w.rho*sw[i];
876  dm[i] = omega[matID]*density_w.rho;
877 
878  /* capillary potential */
879  phi[i] = psk.psic;
880  dphi[i]= psk.dpsic;
881 
882  for (int I=0;I<nSpace;I++)
883  {
884  /* wetting phase advection */
885  for (int m=rowptr[I]; m < rowptr[I+1]; m++)
886  {
887  const int J = colind[m];
888  if (I==J)
889  {
890  f[i*nSpace+I] = advectionScaling*(qt[i*nSpace+I]*fracFlow.fw
891  - Kbar[matID*nnz+m]*fracFlow.lambdaw*fracFlow.fn*(b*density_n.rho-density_w.rho)*g[I]) ;
892  df[i*nSpace+I] = advectionScaling*(qt[i*nSpace+I]*fracFlow.dfw
893  - (Kbar[matID*nnz+m]*g[I]*(b*density_n.rho-density_w.rho))*(fracFlow.lambdaw*fracFlow.dfn + fracFlow.fn*fracFlow.dlambdaw));
894  }
895  /* wetting phase capillary diffusion */
896  /*include scaling factor in case want to turn of capillary diffusion to test hyperbolic approximations*/
897  a[i*nnz+m] = -capillaryDiffusionScaling*Kbar[matID*nnz+m]*fracFlow.lambdaw*fracFlow.fn;
898  da[i*nnz+m] = -capillaryDiffusionScaling*Kbar[matID*nnz+m]*(fracFlow.dlambdaw*fracFlow.fn + fracFlow.lambdaw*fracFlow.dfn);
899  }
900  }
901  }
902  }
903  return 0;
904 }
905 
906 /* pressure equation in fractional flow formulation for incompressible flow
907  heterogeneity is represented by materialTypes
908 */
909 template<class PSK, class DENSITY_N>
910 static inline int twophaseDarcy_compressibleN_split_sd_pressure_het_matType(int nSimplex,
911  int nPointsPerSimplex,
912  int nSpace,
913  int nParams,
914  const int* rowptr,
915  const int* colind,
916  const int* materialTypes,
917  double muw,
918  double mun,
919  const double* omega,
920  const double* Kbar,
921  double b,
922  double capillaryDiffusionScaling,
923  const double* rwork_psk, const int* iwork_psk,
924  const double* rwork_psk_tol,
925  const double* rwork_density_w,
926  const double* rwork_density_n,
927  const double* g,
928  const double* sw,
929  const double* psiw,
930  const double* psin,
931  const double* grad_psic,
932  double* m,
933  double* dm,
934  double* f,
935  double* a)
936 {
937  int matID;
938  PSK psk(rwork_psk,iwork_psk); psk.setTolerances(rwork_psk_tol);
939  CompressibleN_FractionalFlowVariables fracFlow(muw,mun);
940  ConstantDensity density_w(rwork_density_w);
941  DENSITY_N density_n(rwork_density_n);
942  const int nnz = rowptr[nSpace];
943  double drhon_dpsiw;
944  for (int eN=0;eN<nSimplex;eN++)
945  {
946  matID = materialTypes[eN];
947  psk.setParams(&rwork_psk[matID*nParams]);
948  for(int pN=0,i;pN<nPointsPerSimplex;pN++)
949  {
950  i = eN*nPointsPerSimplex+pN;
951  psk.calc(sw[i]);
952  density_w.calc(psiw[i]);
953  density_n.calc(psin[i]);
954 
955  fracFlow.calc(psk,density_w,density_n);
956 
957  /*density of nonwetting phase a function of sw and psiw through psin but ignore sw derivative because of splitting*/
958  drhon_dpsiw = density_n.drho;
959  /*mwf debug
960  printf("twpffc psiw=%g psin=%g rhow=%g rhon=%g drhon=%g \n",psiw[i],psin[i],density_w.rho,density_n.rho,drhon_dpsiw);
961  */
962  m[i] = omega[matID]*density_w.rho*sw[i] + omega[matID]*density_n.rho*(1.0-sw[i]);
963  dm[i]= omega[matID]*drhon_dpsiw*(1.0-sw[i]);
964 
965  for (int I=0;I<nSpace;I++)
966  {
967  /*include scaling factor in case want to turn of capillary diffusion to test hyperbolic approximations*/
968  for (int m=rowptr[I]; m < rowptr[I+1]; m++)
969  {
970  const int J = colind[m];
971  if (I == J)
972  {
973  f[i*nSpace+I] = - capillaryDiffusionScaling*Kbar[matID*nnz+m]*fracFlow.lambdat*(fracFlow.fn*grad_psic[i*nSpace+I]) +
974  Kbar[matID*nnz+m]*fracFlow.lambdat*(density_w.rho + fracFlow.fn*(b*density_n.rho-density_w.rho))*g[I];
975  }
976  a[i*nnz+m] = Kbar[matID*nnz+m]*fracFlow.lambdat;
977  }
978  }
979  }
980  }
981  return 0.0;
982 }
983 
984 template<class PSK>
985 static inline int twophaseDarcy_incompressible_split_pp_sd_saturation_het_matType(int nSimplex,
986  int nPointsPerSimplex,
987  int nSpace,
988  int nParams,
989  const int* rowptr,
990  const int* colind,
991  const int* materialTypes,
992  double muw,
993  double mun,
994  const double* omega,
995  const double* Kbar,
996  double b,
997  double capillaryDiffusionScaling,
998  double advectionScaling,
999  const double* rwork_psk, const int* iwork_psk,
1000  const double* rwork_psk_tol,
1001  const double* rwork_density_w,
1002  const double* rwork_density_n,
1003  const double* g,
1004  const double* qt,
1005  const double* u,//psic u >= 0, sn u < 0
1006  double* sw,
1007  double* m,
1008  double* dm,
1009  double* phi,
1010  double* dphi,
1011  double* f,
1012  double* df,
1013  double* a,
1014  double* da)
1015 {
1016  int matID;
1017  PSK psk(rwork_psk,iwork_psk); psk.setTolerances(rwork_psk_tol);
1018  FractionalFlowVariables fracFlow(muw,mun);
1019  ConstantDensity density_w(rwork_density_w),density_n(rwork_density_n);
1020  const int nnz = rowptr[nSpace];
1021  double psic_eval(0.),dpsic_eval(0.),dsw_du(0.);
1022  //mwf debug
1023  //std::cout<<"entering twpffinc rhow= "<<density_w.rho<<" rhon= "<<density_n.rho<<" b= "<<b<<" ";
1024  //for (int I=0; I < nSpace; I++)
1025  // {
1026  // std::cout<<"g["<<I<<"]= "<<g[I]<<" ";
1027  // }
1028  for (int eN=0;eN<nSimplex;eN++)
1029  {
1030  matID = materialTypes[eN];
1031  psk.setParams(&rwork_psk[matID*nParams]);
1032  for(int pN=0,i;pN<nPointsPerSimplex;pN++)
1033  {
1034  i = eN*nPointsPerSimplex+pN;
1035  if (u[i] < 0.0)
1036  {
1037  sw[i] = psk.Sw_max;
1038  psic_eval = 0.0;
1039  dsw_du = -1.0;
1040  psk.calc(sw[i]); //calculate psk as function of sw=1-u
1041  dpsic_eval=psk.dpsic;
1042  }
1043  else
1044  {
1045  psic_eval = u[i]; dpsic_eval=1.0;
1046  psk.calc_from_psic(psic_eval);
1047  sw[i] = psk.Se*(psk.Sw_max-psk.Sw_min) + psk.Sw_min;
1048  dsw_du = psk.dSe_dpsic/psk.dSe_dSw;
1049  }
1050  fracFlow.calc(psk,density_w,density_n);
1051 
1052  /* wetting phase mass */
1053  m[i] = omega[matID]*density_w.rho*sw[i];
1054  dm[i] = omega[matID]*density_w.rho*dsw_du;
1055 
1056  /* capillary potential */
1057  phi[i] = psic_eval;
1058  dphi[i]= dpsic_eval;
1059 
1060  for (int I=0;I<nSpace;I++)
1061  {
1062  /* wetting phase advection */
1063  /* todo, remove diagonal assumption on K*/
1064  for (int m=rowptr[I]; m < rowptr[I+1]; m++)
1065  {
1066  const int J = colind[m];
1067  if (I==J)
1068  {
1069  f[i*nSpace+I] = advectionScaling*(qt[i*nSpace+I]*fracFlow.fw
1070  - Kbar[matID*nnz+m]*fracFlow.lambdaw*fracFlow.fn*(b*density_n.rho-density_w.rho)*g[I]) ;
1071  df[i*nSpace+I] = advectionScaling*(qt[i*nSpace+I]*fracFlow.dfw
1072  - (Kbar[matID*nnz+m]*g[I]*(b*density_n.rho-density_w.rho))*(fracFlow.lambdaw*fracFlow.dfn + fracFlow.fn*fracFlow.dlambdaw));
1073  }
1074  /* wetting phase capillary diffusion */
1075  /*include scaling factor in case want to turn of capillary diffusion to test hyperbolic approximations*/
1076  a[i*nnz+m] = -capillaryDiffusionScaling*Kbar[matID*nnz+m]*fracFlow.lambdaw*fracFlow.fn;
1077  da[i*nnz+m] = -capillaryDiffusionScaling*Kbar[matID*nnz+m]*(fracFlow.dlambdaw*fracFlow.fn + fracFlow.lambdaw*fracFlow.dfn);
1078  }
1079  }
1080  }
1081  }
1082  return 0;
1083 }
1084 
1085 //loop through and generate a spline table using nknots points given in array domain
1086 //insert values into splineTable in the order
1087 //u_0,..u_{nk-1},uinv_0,...,uinv_{nk-1},krw_0,...,krw_{nk-1},krn_0,...,krn_{nk-1}
1088 //where uinv = psic if u=Sw and vice versa
1089 //assumes all for one media type
1090 template <class PSK>
1091 inline void generateSplineTables(int nknots,
1092  int startIndex,
1093  int calcFlag, //0 --- genate tables for f(S_w), 1, generate tables for f(psi_c)
1094  const double* domain,
1095  const double* rwork_psk,
1096  const int* iwork_psk,
1097  const double* rwork_psk_tol,
1098  double* splineTable)
1099 {
1100  PSK psk(rwork_psk,iwork_psk);
1101  psk.setTolerances(rwork_psk_tol);
1102 
1103  if (calcFlag == 0)
1104  {
1105  for (int i=0; i < nknots; i++)
1106  {
1107  double sw = domain[i];
1108  psk.calc(sw);
1109  //mwf debug
1110  std::cout<<"generate splineTable calcFlag=0 startIndex= "<<startIndex<<" nknots= "<<nknots<<" sw["<<i<<"]= "<<sw<<" psic= "<<psk.psic<<" krw= "<<psk.krw<<" krn= "<<psk.krn<<std::endl;
1111  splineTable[startIndex+i] =sw;
1112  splineTable[startIndex + nknots +i]=psk.psic;
1113  splineTable[startIndex + nknots*2+i]=psk.krw;
1114  splineTable[startIndex + nknots*3+i]=psk.krn;
1115  }
1116  //go ahead and fix enpoints manually
1117  //capillary pressure at dry end
1118  splineTable[startIndex + nknots] = splineTable[startIndex + nknots+1]+1.0;
1119  }
1120  else
1121  {
1122  for (int i=0; i < nknots; i++)
1123  {
1124  double psic = domain[i];
1125  psk.calc_from_psic(psic);
1126  splineTable[startIndex+i] =psic;
1127  splineTable[startIndex + nknots +i]=psk.Se*(psk.Sw_max-psk.Sw_min) + psk.Sw_min;
1128  splineTable[startIndex + nknots*2+i]=psk.krw;
1129  splineTable[startIndex + nknots*3+i]=psk.krn;
1130  }
1131  //fix capillary pressure at wet end
1132  splineTable[startIndex + nknots] = splineTable[startIndex + nknots+1]*(1.0+1.e-8);
1133  //krw at wet end
1134  splineTable[startIndex + nknots*2] = splineTable[startIndex + nknots*2+1];
1135  //krn at wet end
1136  splineTable[startIndex + nknots*3] = splineTable[startIndex + nknots*3+1];
1137 
1138 
1139  }
1140 
1141 }
1142 #endif
densityRelations.h
ConstantDensity
Definition: densityRelations.h:27
twophaseDarcy_fc_sd_het_matType
int twophaseDarcy_fc_sd_het_matType(int nSimplex, int nPointsPerSimplex, int nSpace, int nParams, const int *rowptr, const int *colind, const int *materialTypes, double muw, double mun, const double *omega, const double *Kbar, double b, const double *rwork_psk, const int *iwork_psk, const double *rwork_psk_tol, const double *rwork_density_w, const double *rwork_density_n, const double *g, const double *x, const double *sw, const double *psiw, double *mw, double *dmw_dsw, double *dmw_dpsiw, double *mn, double *dmn_dsw, double *dmn_dpsiw, double *psin, double *dpsin_dsw, double *dpsin_dpsiw, double *phi_psiw, double *dphi_psiw_dpsiw, double *phi_psin, double *dphi_psin_dpsiw, double *dphi_psin_dsw, double *aw, double *daw_dsw, double *daw_dpsiw, double *an, double *dan_dsw, double *dan_dpsiw)
Definition: twophaseDarcyCoefficients.h:12
f
Double f
Definition: Headers.h:64
FractionalFlowVariables
Definition: pskRelations.h:597
twophaseDarcy_fc_pp_sd_het_matType
int twophaseDarcy_fc_pp_sd_het_matType(int nSimplex, int nPointsPerSimplex, int nSpace, int nParams, const int *rowptr, const int *colind, const int *materialTypes, double muw, double mun, const double *omega, const double *Kbar, double b, const double *rwork_psk, const int *iwork_psk, const double *rwork_psk_tol, const double *rwork_density_w, const double *rwork_density_n, const double *g, const double *x, const double *psiw, const double *psic, double *sw, double *mw, double *dmw_dpsiw, double *dmw_dpsic, double *mn, double *dmn_dpsiw, double *dmn_dpsic, double *phi_psiw, double *dphi_psiw_dpsiw, double *phi_psin, double *dphi_psin_dpsiw, double *dphi_psin_dpsic, double *aw, double *daw_dpsiw, double *daw_dpsic, double *an, double *dan_dpsiw, double *dan_dpsic)
Definition: twophaseDarcyCoefficients.h:323
df
double df(double C, double b, double a, int q, int r)
Definition: analyticalSolutions.c:2209
phi
Double phi
Definition: Headers.h:76
twophaseDarcy_fc_sd_het_matType_nonPotentialForm
int twophaseDarcy_fc_sd_het_matType_nonPotentialForm(int compressibilityFlag, int nSimplex, int nPointsPerSimplex, int nSpace, int nParams, const int *rowptr, const int *colind, const int *materialTypes, double muw, double mun, const double *omega, const double *Kbar, double b, const double *rwork_psk, const int *iwork_psk, const double *rwork_psk_tol, const double *rwork_density_w, const double *rwork_density_n, const double *g, const double *x, const double *sw, const double *psiw, double *mw, double *dmw_dsw, double *dmw_dpsiw, double *mn, double *dmn_dsw, double *dmn_dpsiw, double *psin, double *dpsin_dsw, double *dpsin_dpsiw, double *phi_psiw, double *dphi_psiw_dpsiw, double *phi_psin, double *dphi_psin_dpsiw, double *dphi_psin_dsw, double *fw, double *dfw_dsw, double *dfw_dpsiw, double *fn, double *dfn_dsw, double *dfn_dpsiw, double *aw, double *daw_dsw, double *daw_dpsiw, double *an, double *dan_dsw, double *dan_dpsiw)
Definition: twophaseDarcyCoefficients.h:129
u
Double u
Definition: Headers.h:89
CompressibleN_FractionalFlowVariables
Definition: pskRelations.h:639
pskRelations.h
twophaseDarcy_vol_frac
int twophaseDarcy_vol_frac(int nSimplex, int nPointsPerSimplex, const int *materialTypes, const double *omega, const double *sw, double *vol_frac_w, double *vol_frac_n)
Definition: twophaseDarcyCoefficients.h:293
generateSplineTables
void generateSplineTables(int nknots, int startIndex, int calcFlag, const double *domain, const double *rwork_psk, const int *iwork_psk, const double *rwork_psk_tol, double *splineTable)
Definition: twophaseDarcyCoefficients.h:1091
nnz
#define nnz
Definition: Richards.h:9
Generated on Fri Jul 1 2022 11:17:46 for proteus by doxygen 1.8.20