proteus 1.9.0
C/C++/Fortran libraries
Loading...
Searching...
No Matches
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*/
11template<class PSK, class DENSITY_W, class DENSITY_N>
12inline 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}
128template<class PSK, class DENSITY_W, class DENSITY_N>
129inline 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
293inline 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*/
322template<class PSK, class DENSITY_W, class DENSITY_N>
323inline 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 */
477template<class PSK>
478static 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 */
565template<class PSK, class DENSITY_W>
566static 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*/
663template<class PSK>
664static 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*/
727template<class PSK, class DENSITY_W, class DENSITY_N>
728static 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 */
812template<class PSK, class DENSITY_N>
813static 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*/
909template<class PSK, class DENSITY_N>
910static 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
984template<class PSK>
985static 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
1090template <class PSK>
1091inline 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
f
Double f
Definition Headers.h:64
u
Double u
Definition Headers.h:89
phi
Double phi
Definition Headers.h:76
df
double df(double C, double b, double a, int q, int r)
Definition analyticalSolutions.c:2209
nnz
#define nnz
Definition m_comp_co2.h:19
pskRelations.h
Richards free-function PSK closures (moved here from proteus/richards/psk_models.h).
CompressibleN_FractionalFlowVariables
Definition pskRelations.h:640
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
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
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
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
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
Generated on for proteus by doxygen 1.18.0