proteus  1.8.1
C/C++/Fortran libraries
pskRelations.h
Go to the documentation of this file.
1 #ifndef PSKRELATIONS_H
2 #define PSKRELATIONS_H
3 #include <cmath>
4 #include <iostream>
5 #include <cassert>
6 #include "densityRelations.h"
7 #include "SubsurfaceTransportCoefficients.h"
16 /* jcc two phase flow, modified by cek and mwf*/
17 using namespace std;
18 
19 class PskRelation
20 {
21 public:
22  double Se,
23  dSe_dSw,
24  Sw_min,
25  Sw_max,
26  krw,
27  dkrw,
28  krn,
29  dkrn,
30  psic,
31  dpsic,
32  //in case need psic--> Se form
33  dSe_dpsic;
34  PskRelation(const double* rwork, const int* iwork = 0):
35  Se(1.0),dSe_dSw(1.0),
36  Sw_min(rwork[0]),
37  Sw_max(rwork[1]),
38  krw(1.0),dkrw(0.0),
39  krn(0.0),dkrn(0.0),psic(0.0),
40  dpsic(0.0),dSe_dpsic(0.0)
41  {}
42  PskRelation():
43  Se(1.0),dSe_dSw(1.0),
44  Sw_min(0.0),
45  Sw_max(1),
46  krw(1.0),dkrw(0.0),
47  krn(0.0),dkrn(0.0),psic(0.0),
48  dpsic(0.0),dSe_dpsic(0.0)
49  {}
50  virtual inline void setParams(const double* rwork, const int* iwork = 0)
51  {
52  Sw_min = rwork[0];
53  Sw_max = rwork[1];
54  }
55  /*for fudge factors aka tolerances in various models*/
56  virtual inline void setTolerances(const double* rwork_tol)
57  {
58  }
59  virtual ~PskRelation(){}
60 
61  /*linear coefficients*/
62  inline void calc(const double& Sw)
63  {
64  calc_Se(Sw);
65 
66  krw = Se;
67  dkrw = dSe_dSw;
68 
69  krn = (1.0-Se);
70  dkrn = -dSe_dSw;
71 
72  psic = Se;
73  dpsic = dSe_dSw;
74  }
75 
76  inline void calc_Se(const double& Sw)
77  {
78  if (Sw < Sw_min)
79  {
80  Se = 0.0;
81  dSe_dSw = 0.0;
82  }
83  else if (Sw > Sw_max)
84  {
85  Se = 1.0;
86  dSe_dSw =0.0;
87  }
88  else
89  {
90  Se = (Sw - Sw_min)/(Sw_max-Sw_min);
91  dSe_dSw = 1.0/(Sw_max - Sw_min);
92  }
93  }
94  //for Sw as a function of capillary pressure head
95  //note dkrw,dkrn are set to be wrt psic
96  virtual inline void calc_from_psic(const double& psicIn)
97  {
98  bool implemented = false;
99  assert(implemented);
100  }
101 
102 };
103 
104 /* quadratic kr */
105 class SimplePSK : public PskRelation
106 {
107 public:
108  SimplePSK(const double* rwork, const int* iwork = 0):
109  PskRelation(rwork)
110  {}
111 
112  inline void calc(const double& Sw)
113  {
114  calc_Se(Sw);
115 
116  krw = Se*Se;
117  dkrw = 2.0*Se*dSe_dSw;
118 
119  krn = (1.0-Se)*(1.0-Se);
120  dkrn = 2.0*(Se-1.0)*dSe_dSw;
121 
122  psic = (1.0-Se)/Sw_min; /*mwf change to test het -Se;*/
123  dpsic = -dSe_dSw;
124  }
125  //
126  virtual inline void calc_from_psic(const double& psicIn)
127  {
128  psic = max(0.0,psicIn);
129  dpsic= 1.0;
130  Se = 1.0-psic*Sw_min;
131  Se = max(0.0,min(1.0,Se));
132 
133  krw = Se*Se;
134  dkrw = 2.0*Se*dSe_dSw;
135 
136  krn = (1.0-Se)*(1.0-Se);
137  dkrn = 2.0*(Se-1.0)*dSe_dSw;
138 
139  }
140 };
141 
142 /* Van Genuchten-Mualem */
143 class VGMorig: public PskRelation
144 {
145 public:
146  double alpha,
147  m,
148  n,
149  Se_eps;
150  double Se_eps_const;
151  VGMorig()
152  {}
153  VGMorig(const double* rwork, const int* iwork = 0):
154  PskRelation(rwork),
155  alpha(rwork[2]),
156  m(rwork[3])
157  {
158  n = (1.0)/(1.0-m);
159  Se_eps_const=1.0e-4;
160  //mwf debug
161  //std::cout<<"VGMorig ctor rwork[2]= "<<rwork[2]<<" alpha= "<<alpha<<" rwork[3]= "<<rwork[3]<<" m= "<<m<<" n= "<<n<<std::endl;
162  }
163  /*for fudge factors aka tolerances in various models*/
164  virtual inline void setTolerances(const double* rwork_tol)
165  {
166  Se_eps_const = rwork_tol[0];
167  }
168 
169  inline void setParams(const double* rwork, const int* iwork = 0)
170  {
171  Sw_min = rwork[0];
172  Sw_max = rwork[1];
173  alpha = rwork[2];
174  m = rwork[3];
175  n = (1.0)/(1.0-m);
176  //mwf debug
177  //std::cout<<"VGMorig setParams rwork[2]= "<<rwork[2]<<" alpha= "<<alpha<<" rwork[3]= "<<rwork[3]<<" m= "<<m<<" n= "<<n<<std::endl;
178  }
179  inline void calc_Se_eps(const double& Se)
180  {
181  /* cek todo, work out when you need to stay away from Se=0 and 1 and what to assign in these cases */
182  if(Se <= Se_eps_const)
183  Se_eps=Se_eps_const;
184  else if(Se >= (1.0-Se_eps_const))
185  Se_eps=1.0-Se_eps_const;
186  else
187  Se_eps=Se;
188  }
189 
190  inline void calc(const double& Sw)
191  {
192  double Seovmmo,Seovm,Seovmmh,S,Smmo,S2mmo,Sm,S2m,pcesub1,pcesub2;
193  calc_Se(Sw);
194  calc_Se_eps(Se);
195 
196  Seovmmo = pow(Se_eps,((1.0/m)-1.0));
197  Seovm = Seovmmo*Se_eps;
198  Seovmmh = pow(Se_eps,((1.0/m)-0.5));
199  S = 1.0 - Seovm;
200  Smmo = pow(S,m-1.0);
201  S2mmo = pow(S,2.0*m-1.0);
202  Sm = Smmo*S;
203  S2m = S2mmo*S;
204 
205  pcesub1 = pow(((1.0/Seovm)-1.0),((1.0/n)-1.0));
206  pcesub2 = pcesub1*((1.0/Seovm)-1.0);
207 
208  krw = sqrt(Se)*(1.0-Sm)*(1.0-Sm);
209  dkrw = (0.5*(1.0/sqrt(Se_eps))*(1.0-Sm)*(1.0-Sm) + 2.0*(1-Sm)*Smmo*Seovmmh)*dSe_dSw;
210 
211  krn = sqrt(1.0-Se)*S2m;
212  dkrn = (-0.5*(1.0/sqrt(1.0-Se_eps))*S2m - 2.0*sqrt(1.0-Se)*(S2mmo)*Seovmmo)*dSe_dSw;
213 
214  psic = pow((pow(Se_eps,(-1.0/m)) - 1.0),(1.0/n))/alpha;
215  dpsic = ((-1.0/(alpha*n*m))*(pow((pow(Se_eps,-1.0/m)-1.0),-m))*(pow(Se_eps,((-1.0/m)-1.0))))*dSe_dSw;
216  }
217 };
218 
219 class VGM : public VGMorig
220 {
221  public:
222  double ns_del,eps_small,sqrt_eps_small;
223  VGM()
224  {}
225  VGM(const double* rwork, const int* iwork = 0):
226  VGMorig(rwork,iwork),
227  ns_del(1.0e-8),
228  eps_small(1.0e-16),
229  sqrt_eps_small(1.0e-8)
230  {}
231  VGM(const VGM& r)
232  {
233  Se=r.Se;
234  dSe_dSw=r.dSe_dSw;
235  Sw_min=r.Sw_min;
236  Sw_max=r.Sw_max;
237  krw=r.krw;
238  dkrw=r.dkrw;
239  krn=r.krn;
240  dkrn=r.dkrn;
241  psic=r.psic;
242  dpsic=r.dpsic;
243  dSe_dpsic=r.dSe_dpsic;
244  alpha=r.alpha;
245  m=r.m;
246  n=r.n;
247  Se_eps_const=r.Se_eps_const;
248  ns_del=r.ns_del;
249  eps_small=r.eps_small;
250  sqrt_eps_small=r.sqrt_eps_small;
251  }
252  /*for fudge factors aka tolerances in various models*/
253  virtual inline void setTolerances(const double* rwork_tol)
254  {
255  eps_small = rwork_tol[0]; //mwf 072110 don't tie to eps_small? sqrt_eps_small = sqrt(eps_small);
256  ns_del = rwork_tol[1];
257  }
258 
259  inline void calc_Se(const double& Sw)
260  {
261  Se = (Sw - Sw_min)/(Sw_max-Sw_min);
262  dSe_dSw = 1.0/(Sw_max - Sw_min);
263  Se = max(eps_small,min(Se,1.0-eps_small));
264  }
265 
266  inline void calc(const double& Sw)
267  {
268  calc_Se(Sw);
269  //taken from MualemVanGenuchten2p in pdetk
270  double sBar,psiC,DsBar_DpC,DDsBar_DDpC,DkrW_DpC,DkrN_DpC;
271  double vBar,uBar,
272  alphaPsiC, alphaPsiC_n, alphaPsiC_nM1, alphaPsiC_nM2,
273  onePlus_alphaPsiC_n,
274  sqrt_sBar, sqrt_1minusSbar,
275  sBarByOnePlus_alphaPsiC_n, sBarBy_onePlus_alphaPsiC_n_2;
276 
277  sBar = Se;
278  //begin MualemVanGenuchten2p setVFraction
279  onePlus_alphaPsiC_n = pow(sBar,1.0/-m);
280  alphaPsiC_n = onePlus_alphaPsiC_n - 1.0;
281  alphaPsiC = pow(alphaPsiC_n,1.0/n);
282  psiC = alphaPsiC/alpha;
283 
284  alphaPsiC_nM1 = alphaPsiC_n/alphaPsiC;
285  sBarByOnePlus_alphaPsiC_n = sBar/onePlus_alphaPsiC_n;
286  sqrt_sBar = sqrt(sBar);
287  sqrt_1minusSbar = sqrt(1.0 - sBar);
288 
289  DsBar_DpC = -alpha*(n-1.0)*alphaPsiC_nM1
290  *sBarByOnePlus_alphaPsiC_n;
291  //DthetaW_DpC = thetaSR[i] * DsBar_DpC;
292 
293  vBar = 1.0-alphaPsiC_nM1*sBar;
294  uBar = alphaPsiC_nM1*sBar;
295 
296  //change names krW--> krw, krN--> krn
297  krw = sqrt_sBar*vBar*vBar;
298  krn = sqrt_1minusSbar*uBar*uBar;
299  psic= psiC;
300  if(psiC<=0.0)
301  {
302  DsBar_DpC = 0.0;
303  //DthetaW_DpC = 0.0;
304  krw = 1.0;
305  krn = 0.0;
306  }
307 
308  //begin MualemVanGenuchten2p calculateDerivatives
309  alphaPsiC_nM2 = alphaPsiC_nM1/alphaPsiC;
310 
311  sBarBy_onePlus_alphaPsiC_n_2 = sBarByOnePlus_alphaPsiC_n
312  /onePlus_alphaPsiC_n;
313  DDsBar_DDpC = alpha*alpha*(n-1.)
314  *((2*n-1.)*alphaPsiC_nM1*alphaPsiC_nM1
315  *sBarBy_onePlus_alphaPsiC_n_2
316  -
317  (n-1.)*alphaPsiC_nM2
318  *sBarByOnePlus_alphaPsiC_n);
319 
320  //DDthetaW_DDpC = thetaSR[i]*DDsBar_DDpC;
321 
322  DkrW_DpC = (0.5/sqrt_sBar)*DsBar_DpC*vBar*vBar
323  -
324  2.0*sqrt_sBar*vBar*
325  (alpha*(n-1.0)*alphaPsiC_nM2*sBar
326  + alphaPsiC_nM1 * DsBar_DpC);
327 
328  //DKW_DpC = KWs[i]*DkrW_DpC;
329 
330 
331  //recalculate if necessary
332  if (sqrt_1minusSbar >= sqrt_eps_small)//SQRT_MACHINE_EPSILON)
333  {
334  DkrN_DpC = -(0.5/sqrt_1minusSbar)*DsBar_DpC*uBar*uBar
335  +
336  2.0*sqrt_1minusSbar*uBar*
337  (alpha*(n-1.0)*alphaPsiC_nM2*sBar
338  + alphaPsiC_nM1 * DsBar_DpC);
339  }
340  else
341  {
342  DkrN_DpC =((1.0 - sBar)/eps_small)*2.0*sqrt_eps_small*uBar*
343  (alpha*(n-1.0)*alphaPsiC_nM2*sBar
344  + alphaPsiC_nM1 * DsBar_DpC)
345  - (DsBar_DpC/eps_small)*sqrt_eps_small*uBar*uBar;
346  }
347 
348  //if we're in the nonsmooth regime
349  if (psiC < ns_del && psiC > 0.0 )
350  {
351  DkrW_DpC = 0.0;
352  }
353 
354  if (psiC <= 0.0)
355  {
356  DDsBar_DDpC = 0.0;
357  //DDthetaW_DDpC = 0.0;
358  DkrW_DpC = 0.0;
359  DkrN_DpC = 0.0;
360  }
361  //end calculateDerivatives
362  double DpC_Dse = 0.0;
363  if (fabs(DsBar_DpC) > 0.0)
364  DpC_Dse = 1.0/DsBar_DpC;
365  double DpC_Dsw = DpC_Dse*dSe_dSw;
366  dkrw = DkrW_DpC*DpC_Dsw;
367  dkrn = DkrN_DpC*DpC_Dsw;
368  dpsic= DpC_Dsw;
369 
370  }
371  //TODO add this to other classes
372  virtual inline void calc_from_psic(const double& psicIn)
373  {
374  //taken from MualemVanGenuchten2p in pdetk
375  double sBar,psiC,DsBar_DpC,DDsBar_DDpC,DkrW_DpC,DkrN_DpC;
376  double vBar,uBar,
377  alphaPsiC, alphaPsiC_n, alphaPsiC_nM1, alphaPsiC_nM2,
378  onePlus_alphaPsiC_n,
379  sqrt_sBar, sqrt_1minusSbar,
380  sBarByOnePlus_alphaPsiC_n, sBarBy_onePlus_alphaPsiC_n_2;
381 
382 
383  psiC = max(0.0,psicIn);
384  alphaPsiC = alpha*psiC;
385  alphaPsiC_n = pow(alphaPsiC,n);
386  alphaPsiC_nM1 = alphaPsiC_n/alphaPsiC;
387  onePlus_alphaPsiC_n = 1.0 + alphaPsiC_n;
388  sBar = pow(onePlus_alphaPsiC_n,-m);
389  sBarByOnePlus_alphaPsiC_n = sBar/onePlus_alphaPsiC_n;
390  sqrt_sBar = sqrt(sBar);
391  sqrt_1minusSbar = sqrt(1.0 - sBar);
392  //thetaW = thetaSR[i]*sBar + thetaR[i];
393  DsBar_DpC = -alpha*(n-1.0)*alphaPsiC_nM1
394  *sBarByOnePlus_alphaPsiC_n;
395  //DthetaW_DpC = thetaSR[i] * DsBar_DpC;
396  vBar = 1.0-alphaPsiC_nM1*sBar;
397  uBar = alphaPsiC_nM1*sBar;
398 
399  //change names krW--> krw, krN--> krn
400  krw = sqrt_sBar*vBar*vBar;
401  krn = sqrt_1minusSbar*uBar*uBar;
402  Se = sBar;
403  psic=psiC;
404  if(psiC<=0.0)
405  {
406  sBar = 1.0;
407  Se = sBar;
408  //thetaW = thetaS[i];
409  DsBar_DpC = 0.0;
410  //DthetaW_DpC = 0.0;
411  krw = 1.0;
412  krn = 0.0;
413  }
414  //mwf debug
415  //std::cout<<"vgm_calc_from_psic alpha= "<<alpha<<" m= "<<m<<" psic= "<<psic<<" sBar= "<<sBar<<" krw= "<<krw<<" krn= "<<krn
416  // <<std::endl;
417  //begin MualemVanGenuchten2p calculateDerivatives
418  alphaPsiC_nM2 = alphaPsiC_nM1/alphaPsiC;
419 
420  sBarBy_onePlus_alphaPsiC_n_2 = sBarByOnePlus_alphaPsiC_n
421  /onePlus_alphaPsiC_n;
422  DDsBar_DDpC = alpha*alpha*(n-1.)
423  *((2*n-1.)*alphaPsiC_nM1*alphaPsiC_nM1
424  *sBarBy_onePlus_alphaPsiC_n_2
425  -
426  (n-1.)*alphaPsiC_nM2
427  *sBarByOnePlus_alphaPsiC_n);
428 
429  //DDthetaW_DDpC = thetaSR[i]*DDsBar_DDpC;
430 
431  DkrW_DpC = (0.5/sqrt_sBar)*DsBar_DpC*vBar*vBar
432  -
433  2.0*sqrt_sBar*vBar*
434  (alpha*(n-1.0)*alphaPsiC_nM2*sBar
435  + alphaPsiC_nM1 * DsBar_DpC);
436 
437  //DKW_DpC = KWs[i]*DkrW_DpC;
438 
439 
440  //recalculate if necessary
441  if (sqrt_1minusSbar >= sqrt_eps_small)//SQRT_MACHINE_EPSILON)
442  {
443  DkrN_DpC = -(0.5/sqrt_1minusSbar)*DsBar_DpC*uBar*uBar
444  +
445  2.0*sqrt_1minusSbar*uBar*
446  (alpha*(n-1.0)*alphaPsiC_nM2*sBar
447  + alphaPsiC_nM1 * DsBar_DpC);
448  }
449  else
450  {
451  DkrN_DpC =((1.0 - sBar)/eps_small)*2.0*sqrt_eps_small*uBar*
452  (alpha*(n-1.0)*alphaPsiC_nM2*sBar
453  + alphaPsiC_nM1 * DsBar_DpC)
454  - (DsBar_DpC/eps_small)*sqrt_eps_small*uBar*uBar;
455  }
456 
457  //if we're in the nonsmooth regime
458  if (psiC < ns_del && psiC > 0.0 )
459  {
460  DkrW_DpC = 0.0;
461  }
462 
463  if (psiC <= 0.0)
464  {
465  DDsBar_DDpC = 0.0;
466  //DDthetaW_DDpC = 0.0;
467  DkrW_DpC = 0.0;
468  DkrN_DpC = 0.0;
469  }
470  //end calculateDerivatives
471  dkrw = DkrW_DpC; //note: \od{k_{rw}}{\psi_c} not \od{k_{rw}}{S_w}
472  dkrn = DkrN_DpC;
473  dpsic= 1.0;
474  dSe_dpsic=DsBar_DpC;
475  }
476 };
477 /* Van Genuchten-Burdine */
478 class VGB : public VGM
479 {
480 public:
481  VGB(const double* rwork, const int* iwork = 0):
482  VGM(rwork,iwork)
483  {}
484 
485  inline void calc(const double& Sw)
486  {
487  double S,Smmo,Sm,alpha,Se1ovMmo;
488  calc_Se(Sw);
489  calc_Se_eps(Se);
490 
491  Se1ovMmo = pow(Se_eps,((1.0/m)-1.0));
492  S = 1.0-Se1ovMmo*Se_eps;
493  Smmo = pow(S,(m-1.0));
494  Sm = Smmo*S;
495 
496  //cek from matlab, needs optimizing
497 
498  krw = Se_eps*Se_eps*(1.0-pow(1.0-pow(Se_eps,1.0/m),1.0/m));
499  dkrw = (2.0*Se_eps*(1.0-pow(1.0-pow(Se_eps,1.0/m),1.0/m))+Se_eps*pow(1.0-pow(Se_eps,1.0/m),1.0/m)/(m*m)*pow(Se_eps,1.0/m)/(1.0-pow(Se_eps,1.0/m)))*dSe_dSw;
500 
501  krn = pow(1.0-Se_eps,2.0)*pow(1.0-pow(Se_eps,1.0/m),1.0*m);
502  dkrn = (-2.0*(1.0-Se_eps)*pow(1.0-pow(Se_eps,1.0/m),1.0*m)-pow(1.0-Se_eps,2.0)*pow(1.0-pow(Se_eps,1.0/m),1.0*m)*pow(Se_eps,1.0/m)/Se_eps/(1.0-pow(Se_eps,1.0/m)))*dSe_dSw;
503 
504 // krw = (Se*Se)*(1.0-Sm);
505 // dkrw = ((2.0*Se)*(1.0-Sm)+(Se*Se)*(Se1ovMmo)*Smmo)*dSe_dSw;
506 
507 // krn = (1.0-Se)*(1.0-Se)*(Sm);
508 // dkrn = (2.0*(Se-1.0)*Sm-((1.0-Se)*(1.0-Se))*Smmo*Se1ovMmo)*dSe_dSw;
509  psic = pow((pow(Se_eps,(-1.0/m)) - 1.0),(1.0/n))/alpha;
510  dpsic = ((-1.0/(alpha*n*m))*(pow((pow(Se_eps,-1.0/m)-1.0),-m))*(pow(Se_eps,((-1.0/m)-1.0))))*dSe_dSw;
511  }
512 };
513 
514 /* Brooks-Corey-Mualem */
515 class BCM : public PskRelation
516 {
517  public:
518  double pd,lambda;
519 
520  BCM(const double* rwork, const int* iwork = 0):
521  PskRelation(rwork,iwork),
522  pd(rwork[2]),
523  lambda(rwork[3])
524  {}
525 
526  inline void setParams(const double* rwork, const int* iwork = 0)
527  {
528  Sw_min = rwork[0];
529  Sw_max = rwork[1];
530  pd = rwork[2];
531  lambda = rwork[3];
532  }
533  inline void calc(const double& Sw)
534  {
535  double Value,Expon,krwovSe,Oovbclpo,Oovbcl,X,sqrt1mu;
536  calc_Se(Sw);
537 
538  Oovbcl = 1.0/lambda;
539  Oovbclpo = Oovbcl+1.0;
540  X = pow(Se,Oovbcl);
541  Value = 1.0-X*Se;
542  Expon = ((4.0+5.0*lambda)/(2.0*lambda));
543  krwovSe = pow(Se,(Expon-1.0));
544  sqrt1mu = sqrt(1.0-Se);
545 
546  krw = krwovSe*Se;
547  dkrw = (Expon*krwovSe)*dSe_dSw;
548 
549  krn = sqrt1mu*Value*Value;
550  dkrn = (-0.5*(1.0/sqrt1mu)*Value*Value - 2.0*sqrt1mu*Value*Oovbclpo*X )*dSe_dSw;
551 
552  psic = 1.0/X;
553  dpsic = (-Oovbcl/(X*Se))*dSe_dSw;
554  }
555 };
556 
557 /* Brooks-Corey-Burdine */
558 class BCB : public BCM
559 {
560 public:
561  BCB(const double* rwork, const int* iwork = 0):BCM(rwork,iwork)
562  {}
563 
564  inline void calc(const double& Sw)
565  {
566  double Se2ovL,Se2,Se3,omSe,Expon,Semoovlmo,Se_cutOff;
567  calc_Se(Sw);
568 
569  Se2ovL = pow(Se,(2.0/lambda));
570  Se2 = Se*Se;
571  Se3 = Se*Se2;
572  omSe = 1.0-Se;
573  Expon = (2.0+3.0*lambda)/lambda;
574  Se_cutOff = max(1.0e-4,Se);
575  Semoovlmo = pow(Se_cutOff,((-1.0/lambda)-1.0));
576 
577  krw = Se2ovL*Se3;
578  dkrw = Expon*Se2ovL*Se2*dSe_dSw;
579 
580  krn = (omSe*omSe)*(1.0-Se2ovL*Se);
581  dkrn = (2.*omSe*(1.0-Se2ovL*Se) - (omSe*omSe)*(Expon-2.0)*Se2ovL*dSe_dSw)*dSe_dSw;
582 
583  /* cek debug */
584  krw = pow(Se,(2.0+3.0*lambda)/lambda);
585  dkrw = (((2.0+3.0*lambda)/lambda)*pow(Se,(2.0+3.0*lambda)/lambda - 1.0))*dSe_dSw;
586 
587  krn = (1.0-Se)*(1.0-Se)*(1.0-pow(Se,(2.0+lambda)/lambda));
588  dkrn = (-2.0*(1.0-Se)*(1.0-pow(Se,(2.0+lambda)/lambda))
589  -((2.0+lambda)/lambda)*(1.0-Se)*(1.0-Se)*pow(Se,(2.0+lambda)/lambda-1.0))*dSe_dSw;
590 
591  psic = pd*Semoovlmo*Se;
592  dpsic = pd*(-1.0/lambda)*Semoovlmo*dSe_dSw;
593  }
594 };
595 
596 class FractionalFlowVariables
597 {
598 public:
599  double muw,
600  mun,
601  lambdaw,
602  dlambdaw,
603  lambdan,
604  dlambdan,
605  lambdat,
606  dlambdat,
607  fw,
608  dfw,
609  fn,
610  dfn;
611 
612  FractionalFlowVariables(double muwIn,double munIn):
613  muw(muwIn),
614  mun(munIn)
615  {}
616 
617  inline void calc(const PskRelation& psk,
618  const DensityRelation& density_w,
619  const DensityRelation& density_n)
620  {
621  lambdaw = density_w.rho*psk.krw/muw;
622  dlambdaw =(density_w.rho/muw)*psk.dkrw;
623 
624  lambdan =(density_n.rho*psk.krn)/mun;
625  dlambdan =(density_n.rho/mun)*psk.dkrn;
626 
627  lambdat = lambdaw + lambdan;
628  dlambdat = dlambdaw + dlambdan;
629 
630  fw = lambdaw/lambdat;
631  dfw = (dlambdaw*lambdat - lambdaw*dlambdat)/(lambdat*lambdat);
632 
633  fn = lambdan/lambdat;
634  dfn = (dlambdan*lambdat - lambdan*dlambdat)/(lambdat*lambdat);
635  }
636 };
637 
638 struct CompressibleN_FractionalFlowVariables : public FractionalFlowVariables
639 {
640 public:
641  CompressibleN_FractionalFlowVariables(double muwIn,double munIn):
642  FractionalFlowVariables(muwIn,munIn)
643  {}
644 
645  double drhon,
646  dlambdaw_psiw,
647  drhon_psiw,
648  dlambdan_psiw,
649  dlambdat_psiw,
650  dfw_psiw,
651  dfn_psiw;
652 
653  inline void calc(const PskRelation& psk,
654  const DensityRelation& density_w,
655  const DensityRelation& density_n)
656  {
657  lambdaw = density_w.rho*psk.krw/muw;
658  dlambdaw = (density_w.rho/muw)*psk.dkrw;
659  dlambdaw_psiw = 0.0;
660 
661  drhon = density_n.drho*psk.dpsic;
662  drhon_psiw = density_n.drho;
663 
664  lambdan = (density_n.rho*psk.krn)/mun;
665  dlambdan = (1.0/mun)*(psk.dkrn*density_n.rho + drhon*psk.krn);
666  dlambdan_psiw = drhon_psiw*(psk.krn/mun);
667 
668  lambdat = lambdaw + lambdan;
669  dlambdat = dlambdaw + dlambdan;
670  dlambdat_psiw = dlambdaw_psiw + dlambdan_psiw;
671 
672  fw = lambdaw/lambdat;
673  dfw = (dlambdaw*lambdat - lambdaw*dlambdat)/(lambdat*lambdat);
674  dfw_psiw = (dlambdaw_psiw*lambdat - lambdaw*dlambdat_psiw)/(lambdat*lambdat);
675 
676  fn = lambdan/lambdat;
677  dfn = (dlambdan*lambdat - lambdan*dlambdat)/(lambdat*lambdat);
678  dfn_psiw = (dlambdan_psiw*lambdat - lambdan*dlambdat_psiw)/(lambdat*lambdat);
679  }
680 };
681 
682 
683 class PskSpline: public PskRelation
684 {
685  /*************************************************************
686  spline psk relations. tables of values are held (externally)
687  in splineArray in the order
688 
689  u (sw, or psic), u^-1 (either sw or psic), krw,krn
690 
691  force dSe_dSw = 1.0, so that will be consistent with analytical evaluations
692  where dkrw, etc are wrt to se
693 
694  ************************************************************/
695 public:
696  PskSpline(const double* rworkIn, const int* iworkIn=0):
697  PskRelation(),
698  nknots(2),
699  lastIndex(0),
700  uinvOffset(1),
701  krwOffset(2),
702  krnOffset(3),
703  splineArray(rworkIn)
704  {
705  assert(iworkIn);
706  nknots = iworkIn[0];
707  }
708  virtual ~PskSpline() {}
709  virtual inline void calc(const double& Sw)
710  {
711  assert(splineArray);
712  //
713  Se = Sw; //note force Se = Sw since splines are evaluated directly
714  dSe_dSw = 1.0;
715  piecewiseLinearTableLookup(Sw,
716  nknots,
717  &lastIndex,
718  &psic,
719  &dpsic,
720  splineArray,
721  splineArray+uinvOffset*nknots);
722 
723  //
724  piecewiseLinearTableLookup(Sw,
725  nknots,
726  &lastIndex,
727  &krw,
728  &dkrw,
729  splineArray,
730  splineArray+krwOffset*nknots);
731  //
732  piecewiseLinearTableLookup(Sw,
733  nknots,
734  &lastIndex,
735  &krn,
736  &dkrn,
737  splineArray,
738  splineArray+krnOffset*nknots);
739 
740  }
741  virtual inline void calc_from_psic(const double& psicIn)
742  {
743  assert(splineArray);
744  //
745  psic = psicIn;
746  piecewiseLinearTableLookup(psic,
747  nknots,
748  &lastIndex,
749  &Se, //same as Sw
750  &dSe_dpsic, //same as dSw_dpsic
751  splineArray,
752  splineArray+uinvOffset*nknots);
753 
754  //
755  piecewiseLinearTableLookup(psic,
756  nknots,
757  &lastIndex,
758  &krw,
759  &dkrw,
760  splineArray,
761  splineArray+krwOffset*nknots);
762  //
763  piecewiseLinearTableLookup(psic,
764  nknots,
765  &lastIndex,
766  &krn,
767  &dkrn,
768  splineArray,
769  splineArray+krnOffset*nknots);
770 
771  }
772  virtual inline void setParams(const double* rwork, const int* iwork = 0)
773  {
774  splineArray = rwork;
775  if (iwork)
776  nknots = iwork[0];
777  }
778 public:
779  int nknots,lastIndex,uinvOffset,krwOffset,krnOffset;
780  const double* splineArray;
781 };
783 #endif
PskRelation::~PskRelation
virtual ~PskRelation()
Definition: pskRelations.h:59
FractionalFlowVariables::fn
double fn
Definition: pskRelations.h:609
CompressibleN_FractionalFlowVariables::dfw_psiw
double dfw_psiw
Definition: pskRelations.h:650
PskRelation::dpsic
double dpsic
Definition: pskRelations.h:31
CompressibleN_FractionalFlowVariables::drhon
double drhon
Definition: pskRelations.h:645
densityRelations.h
PskRelation::dSe_dpsic
double dSe_dpsic
Definition: pskRelations.h:33
PskRelation::dkrn
double dkrn
Definition: pskRelations.h:29
VGM::VGM
VGM(const double *rwork, const int *iwork=0)
Definition: pskRelations.h:225
FractionalFlowVariables::dfw
double dfw
Definition: pskRelations.h:608
PskRelation::Sw_min
double Sw_min
Definition: pskRelations.h:24
FractionalFlowVariables::dlambdat
double dlambdat
Definition: pskRelations.h:606
SubsurfaceTransportCoefficients.h
PskRelation::krn
double krn
Definition: pskRelations.h:28
VGM::setTolerances
virtual void setTolerances(const double *rwork_tol)
Definition: pskRelations.h:253
VGMorig::calc
void calc(const double &Sw)
Definition: pskRelations.h:190
CompressibleN_FractionalFlowVariables::dlambdat_psiw
double dlambdat_psiw
Definition: pskRelations.h:649
VGMorig::n
double n
Definition: pskRelations.h:148
VGMorig::m
double m
Definition: pskRelations.h:147
VGB
Definition: pskRelations.h:479
VGMorig::VGMorig
VGMorig()
Definition: pskRelations.h:151
BCM::BCM
BCM(const double *rwork, const int *iwork=0)
Definition: pskRelations.h:520
PskSpline::~PskSpline
virtual ~PskSpline()
Definition: pskRelations.h:708
VGMorig::calc_Se_eps
void calc_Se_eps(const double &Se)
Definition: pskRelations.h:179
DensityRelation::rho
double rho
Definition: densityRelations.h:23
CompressibleN_FractionalFlowVariables::calc
void calc(const PskRelation &psk, const DensityRelation &density_w, const DensityRelation &density_n)
Definition: pskRelations.h:653
CompressibleN_FractionalFlowVariables::dfn_psiw
double dfn_psiw
Definition: pskRelations.h:651
FractionalFlowVariables
Definition: pskRelations.h:597
VGMorig::setTolerances
virtual void setTolerances(const double *rwork_tol)
Definition: pskRelations.h:164
PskRelation::PskRelation
PskRelation(const double *rwork, const int *iwork=0)
Definition: pskRelations.h:34
DensityRelation::drho
double drho
Definition: densityRelations.h:23
CompressibleN_FractionalFlowVariables::drhon_psiw
double drhon_psiw
Definition: pskRelations.h:647
PskRelation::PskRelation
PskRelation()
Definition: pskRelations.h:42
PskRelation::dkrw
double dkrw
Definition: pskRelations.h:27
n
Int n
Definition: Headers.h:28
FractionalFlowVariables::dlambdan
double dlambdan
Definition: pskRelations.h:604
PskSpline::splineArray
const double * splineArray
Definition: pskRelations.h:780
PskSpline::setParams
virtual void setParams(const double *rwork, const int *iwork=0)
Definition: pskRelations.h:772
SimplePSK::calc
void calc(const double &Sw)
Definition: pskRelations.h:112
SimplePSK
Definition: pskRelations.h:106
PskRelation::calc
void calc(const double &Sw)
Definition: pskRelations.h:62
VGMorig::Se_eps
double Se_eps
Definition: pskRelations.h:149
SimplePSK::SimplePSK
SimplePSK(const double *rwork, const int *iwork=0)
Definition: pskRelations.h:108
min
#define min(a, b)
Definition: jf.h:71
VGM::calc
void calc(const double &Sw)
Definition: pskRelations.h:266
FractionalFlowVariables::lambdaw
double lambdaw
Definition: pskRelations.h:601
FractionalFlowVariables::lambdan
double lambdan
Definition: pskRelations.h:603
VGM::VGM
VGM(const VGM &r)
Definition: pskRelations.h:231
VGM::sqrt_eps_small
double sqrt_eps_small
Definition: pskRelations.h:222
FractionalFlowVariables::dlambdaw
double dlambdaw
Definition: pskRelations.h:602
PskSpline::calc
virtual void calc(const double &Sw)
Definition: pskRelations.h:709
VGM::calc_from_psic
virtual void calc_from_psic(const double &psicIn)
Definition: pskRelations.h:372
VGM::VGM
VGM()
Definition: pskRelations.h:223
SimplePSK::calc_from_psic
virtual void calc_from_psic(const double &psicIn)
Definition: pskRelations.h:126
FractionalFlowVariables::fw
double fw
Definition: pskRelations.h:607
VGM::calc_Se
void calc_Se(const double &Sw)
Definition: pskRelations.h:259
BCB::BCB
BCB(const double *rwork, const int *iwork=0)
Definition: pskRelations.h:561
BCM
Definition: pskRelations.h:516
FractionalFlowVariables::lambdat
double lambdat
Definition: pskRelations.h:605
PskSpline::uinvOffset
int uinvOffset
Definition: pskRelations.h:779
FractionalFlowVariables::muw
double muw
Definition: pskRelations.h:599
PskRelation::Se
double Se
Definition: pskRelations.h:22
PskSpline
Definition: pskRelations.h:684
CompressibleN_FractionalFlowVariables
Definition: pskRelations.h:639
PskRelation::setParams
virtual void setParams(const double *rwork, const int *iwork=0)
Definition: pskRelations.h:50
VGMorig::VGMorig
VGMorig(const double *rwork, const int *iwork=0)
Definition: pskRelations.h:153
VGB::VGB
VGB(const double *rwork, const int *iwork=0)
Definition: pskRelations.h:481
VGMorig::setParams
void setParams(const double *rwork, const int *iwork=0)
Definition: pskRelations.h:169
PskRelation::calc_Se
void calc_Se(const double &Sw)
Definition: pskRelations.h:76
PskRelation::dSe_dSw
double dSe_dSw
Definition: pskRelations.h:23
PskRelation::Sw_max
double Sw_max
Definition: pskRelations.h:25
BCB::calc
void calc(const double &Sw)
Definition: pskRelations.h:564
CompressibleN_FractionalFlowVariables::dlambdaw_psiw
double dlambdaw_psiw
Definition: pskRelations.h:646
pd
#define pd(x)
Definition: jf.h:24
VGB::calc
void calc(const double &Sw)
Definition: pskRelations.h:485
PskRelation::setTolerances
virtual void setTolerances(const double *rwork_tol)
Definition: pskRelations.h:56
VGMorig::alpha
double alpha
Definition: pskRelations.h:146
CompressibleN_FractionalFlowVariables::dlambdan_psiw
double dlambdan_psiw
Definition: pskRelations.h:648
VGMorig::Se_eps_const
double Se_eps_const
Definition: pskRelations.h:150
PskRelation::psic
double psic
Definition: pskRelations.h:30
PskRelation::calc_from_psic
virtual void calc_from_psic(const double &psicIn)
Definition: pskRelations.h:96
r
Double r
Definition: Headers.h:83
FractionalFlowVariables::mun
double mun
Definition: pskRelations.h:600
FractionalFlowVariables::dfn
double dfn
Definition: pskRelations.h:610
CompressibleN_FractionalFlowVariables::CompressibleN_FractionalFlowVariables
CompressibleN_FractionalFlowVariables(double muwIn, double munIn)
Definition: pskRelations.h:641
FractionalFlowVariables::calc
void calc(const PskRelation &psk, const DensityRelation &density_w, const DensityRelation &density_n)
Definition: pskRelations.h:617
VGM
Definition: pskRelations.h:220
BCM::calc
void calc(const double &Sw)
Definition: pskRelations.h:533
DensityRelation
Definition: densityRelations.h:21
BCM::pd
double pd
Definition: pskRelations.h:518
max
#define max(a, b)
Definition: analyticalSolutions.h:14
BCB
Definition: pskRelations.h:559
BCM::setParams
void setParams(const double *rwork, const int *iwork=0)
Definition: pskRelations.h:526
piecewiseLinearTableLookup
void piecewiseLinearTableLookup(double x, int nv, int *start, double *y, double *dy, const double *xv, const double *yv)
Definition: SubsurfaceTransportCoefficients.cpp:62
VGMorig
Definition: pskRelations.h:144
PskSpline::PskSpline
PskSpline(const double *rworkIn, const int *iworkIn=0)
Definition: pskRelations.h:696
PskRelation::krw
double krw
Definition: pskRelations.h:26
PskRelation
Definition: pskRelations.h:20
PskSpline::calc_from_psic
virtual void calc_from_psic(const double &psicIn)
Definition: pskRelations.h:741
FractionalFlowVariables::FractionalFlowVariables
FractionalFlowVariables(double muwIn, double munIn)
Definition: pskRelations.h:612
Generated on Fri Jul 1 2022 11:17:40 for proteus by doxygen 1.8.20