OpenCAEPoro  v0.5.0
A simulator for multicomponent porous media flow
MixtureComp.cpp
Go to the documentation of this file.
1 
12 #include "MixtureComp.hpp"
13 
14 COMP::COMP(const vector<string>& comp)
15 {
16  name = comp[0];
17  Pc = stod(comp[1]);
18  Tc = stod(comp[2]);
19  acf = stod(comp[3]);
20  MW = stod(comp[4]);
21  VcMW = stod(comp[5]);
22  Vc = MW * VcMW;
23  OmegaA = stod(comp[6]);
24  OmegaB = stod(comp[7]);
25  Vshift = stod(comp[8]);
26 }
27 
28 MixtureComp::MixtureComp(const ComponentParam& param, const USI& tarId)
29 {
30  mixtureType = EOS_PVTW;
31  // if Water don't exist?
32  // for Mixture class
33  // Now, only one case is considered: oil, gas, water could exist
34  numPhase = param.numPhase + 1;
35  numCom = param.numCom + 1;
36  Allocate();
37  vjp.resize(numPhase, 0);
38  vji.resize(numPhase);
39  for (auto& v : vji) {
40  v.resize(numCom, 0);
41  }
42  rhoN.resize(numPhase * numCom);
43  xiN.resize(numPhase * numCom);
44  muN.resize(numPhase * numCom);
45 
46  // for MixtureComp class
47  NC = param.numCom;
48  NPmax = param.numPhase;
49 
50  zi.resize(NC);
51 
52  Cname = param.Cname;
53  if (param.Tc.activity)
54  Tc = param.Tc.data[tarId];
55  else
56  OCP_ABORT("TCRIT hasn't been input!");
57  if (param.Pc.activity)
58  Pc = param.Pc.data[tarId];
59  else
60  OCP_ABORT("PCRIT hasn't been input!");
61 
62  if (param.Vc.activity)
63  Vc = param.Vc.data[tarId];
64  else if (param.Zc.activity) {
65  Zc = param.Zc.data[tarId];
66  Vc.resize(NC);
67  for (USI i = 0; i < NC; i++) {
68  Vc[i] = 10.73159 * Zc[i] * Tc[i] / Pc[i];
69  }
70  } else
71  OCP_ABORT("VCRIT or ZCRIT hasn't been input!");
72 
73  if (param.MW.activity)
74  MWC = param.MW.data[tarId];
75  else
76  OCP_ABORT("MW hasn't been input!");
77  if (param.Acf.activity)
78  Acf = param.Acf.data[tarId];
79  else
80  OCP_ABORT("ACF hasn't been input!");
81  if (param.OmegaA.activity)
82  OmegaA = param.OmegaA.data[tarId];
83  else
84  OmegaA.resize(NC, 0.457235529);
85  if (param.OmegaB.activity)
86  OmegaB = param.OmegaB.data[tarId];
87  else
88  OmegaB.resize(NC, 0.077796074);
89 
90  if (param.Vshift.activity) {
91  Vshift = param.Vshift.data[tarId];
92  for (USI i = 0; i < NC; i++)
93  Vshift[i] *= (GAS_CONSTANT * OmegaB[i] * Tc[i] / Pc[i]);
94  } else
95  Vshift.resize(NC, 0);
96 
97  if (param.Vcvis.activity)
98  Vcvis = param.Vcvis.data[tarId];
99  else if (param.Zcvis.activity) {
100  Zcvis = param.Zcvis.data[tarId];
101  Vcvis.resize(NC);
102  for (USI i = 0; i < NC; i++) {
103  Vcvis[i] = GAS_CONSTANT * Zcvis[i] * Tc[i] / Pc[i];
104  }
105  } else
106  Vcvis = Vc;
107 
108  LBCcoef = param.LBCcoef;
109  for (auto& lbc : LBCcoef) {
110  lbc *= 10;
111  }
112 
113  InputMiscibleParam(param, tarId);
114 
115  CallId();
116 
117  USI len = NC * NC;
118  BIC.resize(len, 0);
119 
120  if (param.BIC[tarId].size() != len) {
121  USI iter = 0;
122  for (USI i = 1; i < NC; i++) {
123  for (USI j = 0; j < i; j++) {
124  BIC[i * NC + j] = param.BIC[tarId][iter];
125  BIC[j * NC + i] = BIC[i * NC + j];
126  iter++;
127  }
128  }
129  } else {
130  BIC = param.BIC[tarId];
131  }
132 
133  for (USI i = 0; i < NC; i++) {
134  for (USI j = 0; j < NC; j++) {
135  cout << setw(10) << BIC[i * NC + j] << " ";
136  }
137  cout << endl;
138  }
139 
140  EoSctrl.SSMsta.maxIt = stoi(param.SSMparamSTA[0]);
141  EoSctrl.SSMsta.tol = stod(param.SSMparamSTA[1]);
142  EoSctrl.SSMsta.eYt = stod(param.SSMparamSTA[2]);
143  EoSctrl.SSMsta.tol2 = EoSctrl.SSMsta.tol * EoSctrl.SSMsta.tol;
144 
145  EoSctrl.NRsta.maxIt = stoi(param.NRparamSTA[0]);
146  EoSctrl.NRsta.tol = stod(param.NRparamSTA[1]);
147  EoSctrl.NRsta.tol2 = EoSctrl.NRsta.tol * EoSctrl.NRsta.tol;
148 
149  EoSctrl.SSMsp.maxIt = stoi(param.SSMparamSP[0]);
150  EoSctrl.SSMsp.tol = stod(param.SSMparamSP[1]);
151  EoSctrl.SSMsp.tol2 = EoSctrl.SSMsp.tol * EoSctrl.SSMsp.tol;
152 
153  EoSctrl.NRsp.maxIt = stoi(param.NRparamSP[0]);
154  EoSctrl.NRsp.tol = stod(param.NRparamSP[1]);
155  EoSctrl.NRsp.tol2 = EoSctrl.NRsp.tol * EoSctrl.NRsp.tol;
156 
157  EoSctrl.RR.maxIt = stoi(param.RRparam[0]);
158  EoSctrl.RR.tol = stod(param.RRparam[1]);
159  EoSctrl.RR.tol2 = EoSctrl.RR.tol * EoSctrl.RR.tol;
160 
161  AllocateEoS();
162  AllocatePhase();
163  AllocateMethod();
164  AllocateOthers();
165 }
166 
167 void MixtureComp::Flash(const OCP_DBL& Pin, const OCP_DBL& Tin, const OCP_DBL* Niin)
168 {
169  ftype = 0;
170  lNP = 0;
171  InitPTN(Pin, Tin + CONV5, Niin);
172  CalFlash();
173 
174  // Water Properties
175  const USI Wpid = numPhase - 1;
176  const USI Wcid = numCom - 1;
177  phaseExist[Wpid] = OCP_TRUE;
178  xij[Wpid * numCom + Wcid] = 1.0;
179  Nt = Nh + Ni[Wcid];
180  PVTW.Eval_All(0, P, data, cdata);
181  OCP_DBL Pw0 = data[0];
182  OCP_DBL bw0 = data[1];
183  OCP_DBL cbw = data[2];
184  OCP_DBL bw = bw0 * (1 - cbw * (P - Pw0));
185  mu[Wpid] = data[3];
186  xi[Wpid] = 1 / (CONV1 * bw);
187  rho[Wpid] = std_RhoW / bw;
188  vj[Wpid] = CONV1 * Ni[Wcid] * bw;
189  vf += vj[Wpid];
190 
191  // Calculate Sj
192  CalSaturation();
193 }
194 
195 void MixtureComp::InitFlashIMPEC(const OCP_DBL& Pin,
196  const OCP_DBL& Pbbin,
197  const OCP_DBL& Tin,
198  const OCP_DBL* Sjin,
199  const OCP_DBL& Vpore,
200  const OCP_DBL* Ziin,
201  const OCP_USI& bId)
202 {
203  // Attention: zi[numCom - 1] = 0 here, that's Zw = 0;
204  SetBulkId(bId);
205  ftype = 0;
206  lNP = 0;
207 
208  InitPTZ(Pin, Tin + CONV5, Ziin);
209  PhaseEquilibrium();
210  // Attention Nt = 1 now
211  CalMW();
212  CalVfXiRho();
213  CalViscosity();
214  CalSurfaceTension();
215  IdentifyPhase();
216  CopyPhase();
217 
218  // Calulate Nt, water is exclued
219  OCP_DBL Sw = Sjin[numPhase - 1];
220  Nh = Vpore * (1 - Sw) / vf;
221 
222  // Next, nu represents moles of phase instead of molar fraction of phase
223  Dscalar(NP, Nh, &nu[0]);
224  // correct vj, vf with new Nt
225  Dscalar(NPmax, Nh, &vj[0]);
226  vf *= Nh;
227  // CalVfiVfp_full01();
228  CalVfiVfp_full02();
229  // Calculate Ni
230  for (USI i = 0; i < NC; i++) {
231  Ni[i] = zi[i] * Nh;
232  }
233 
234  // Water Properties
235  USI Wpid = numPhase - 1;
236  USI Wcid = numCom - 1;
237  phaseExist[Wpid] = OCP_TRUE;
238  xij[Wpid * numCom + Wcid] = 1.0;
239  PVTW.Eval_All(0, P, data, cdata);
240  OCP_DBL Pw0 = data[0];
241  OCP_DBL bw0 = data[1];
242  OCP_DBL cbw = data[2];
243  OCP_DBL bw = bw0 * (1 - cbw * (P - Pw0));
244  OCP_DBL bwp = -cbw * bw0;
245  mu[Wpid] = data[3];
246  xi[Wpid] = 1 / (CONV1 * bw);
247  rho[Wpid] = std_RhoW / bw;
248  Ni[Wcid] = Vpore * Sw * xi[Wpid];
249  Nt = Nh + Ni[Wcid];
250  vj[Wpid] = CONV1 * Ni[Wcid] * bw;
251  vf += vj[Wpid];
252  vfi[Wcid] = CONV1 * bw;
253  vfP += CONV1 * Ni[Wcid] * bwp;
254 
255  // Calculate Sj
256  CalSaturation();
257 
258  CalSkipForNextStep();
259 }
260 
261 void MixtureComp::InitFlashFIM(const OCP_DBL& Pin,
262  const OCP_DBL& Pbbin,
263  const OCP_DBL& Tin,
264  const OCP_DBL* Sjin,
265  const OCP_DBL& Vpore,
266  const OCP_DBL* Ziin,
267  const OCP_USI& bId)
268 {
269  // Attention: zi[numCom - 1] = 0 here, that's Zw = 0;
270  SetBulkId(bId);
271  ftype = 0;
272  lNP = 0;
273 
274  InitPTZ(Pin, Tin + CONV5, Ziin);
275  PhaseEquilibrium();
276  // Attention Nt = 1 now
277  CalMW();
278  CalVfXiRho();
279  CalViscosity();
280  CalSurfaceTension();
281  IdentifyPhase();
282  CopyPhase();
283 
284  // Calulate Nt, water is exclued
285  OCP_DBL Sw = Sjin[numPhase - 1];
286  Nh = Vpore * (1 - Sw) / vf;
287 
288  // Next, nu represents moles of phase instead of molar fraction of phase
289  Dscalar(NP, Nh, &nu[0]);
290  // correct vj, vf with new Nt
291  Dscalar(NPmax, Nh, &vj[0]);
292  vf *= Nh;
293  // CalVfiVfp_full01();
294  CalVfiVfp_full02();
295  // Calculate Ni
296  for (USI i = 0; i < NC; i++) {
297  Ni[i] = zi[i] * Nh;
298  }
299 
300  // Water Properties
301  USI Wpid = numPhase - 1;
302  USI Wcid = numCom - 1;
303  phaseExist[Wpid] = OCP_TRUE;
304  xij[Wpid * numCom + Wcid] = 1.0;
305  PVTW.Eval_All(0, P, data, cdata);
306  OCP_DBL Pw0 = data[0];
307  OCP_DBL bw0 = data[1];
308  OCP_DBL cbw = data[2];
309  OCP_DBL bw = bw0 * (1 - cbw * (P - Pw0));
310  OCP_DBL bwp = -cbw * bw0;
311  mu[Wpid] = data[3];
312  xi[Wpid] = 1 / (CONV1 * bw);
313  rho[Wpid] = std_RhoW / bw;
314  muP[Wpid] = cdata[3];
315  xiP[Wpid] = -bwp / (bw * bw * CONV1);
316  rhoP[Wpid] = CONV1 * xiP[Wpid] * std_RhoW;
317  Ni[Wcid] = Vpore * Sw * xi[Wpid];
318  nj[Wpid] = Ni[Wcid];
319  Nt = Nh + Ni[Wcid];
320  vj[Wpid] = CONV1 * Ni[Wcid] * bw;
321  vfi[Wcid] = CONV1 * bw;
322  const OCP_DBL vwp = CONV1 * Ni[Wcid] * bwp;
323  vf += vj[Wpid];
324  vfP += vwp;
325  vji[numPhase - 1][numCom - 1] = vfi[Wcid];
326  vjp[numPhase - 1] = vwp;
327 
328  CalSaturation();
329 
330 #ifdef OCP_OLD_FIM
331  CaldXsdXpAPI02();
332 #else
333  CaldXsdXpAPI02p();
334 #endif // OCP_OLD_FIM
335  CalXiPNX_partial();
336  CalRhoPX_partial();
337  CalMuPX_partial();
338 
339  // Calculate pSderExist and pVnumCom
340  fill(pSderExist.begin(), pSderExist.end(), OCP_FALSE);
341  fill(pVnumCom.begin(), pVnumCom.end(), 0.0);
342  if (phaseExist[0]) {
343  pSderExist[0] = OCP_TRUE;
344  pVnumCom[0] = NC;
345  }
346  if (phaseExist[1]) {
347  pSderExist[1] = OCP_TRUE;
348  pVnumCom[1] = NC;
349  }
350  pSderExist[2] = OCP_TRUE;
351 
352  CalSkipForNextStep();
353 }
354 
355 void MixtureComp::InitFlashFIMn(const OCP_DBL& Pin,
356  const OCP_DBL& Pbbin,
357  const OCP_DBL& Tin,
358  const OCP_DBL* Sjin,
359  const OCP_DBL& Vpore,
360  const OCP_DBL* Ziin,
361  const OCP_USI& bId)
362 {
363  // Attention: zi[numCom - 1] = 0 here, that's Zw = 0;
364  SetBulkId(bId);
365  ftype = 0;
366  lNP = 0;
367 
368  InitPTZ(Pin, Tin + CONV5, Ziin);
369  PhaseEquilibrium();
370  // Attention Nt = 1 now
371  CalMW();
372  CalVfXiRho();
373  CalViscosity();
374  CalSurfaceTension();
375  IdentifyPhase();
376  CopyPhase();
377 
378  // Calulate Nt, water is exclued
379  OCP_DBL Sw = Sjin[numPhase - 1];
380  Nh = Vpore * (1 - Sw) / vf;
381 
382  // Next, nu represents moles of phase instead of molar fraction of phase
383  Dscalar(NP, Nh, &nu[0]);
384  // correct vj, vf with new Nt
385  Dscalar(NPmax, Nh, &vj[0]);
386  vf *= Nh;
387  // Calculate Ni
388  for (USI i = 0; i < NC; i++) {
389  Ni[i] = zi[i] * Nh;
390  }
391 
392  CalVjpVfpVfn_partial();
393  // Water Properties
394  USI Wpid = numPhase - 1;
395  USI Wcid = numCom - 1;
396  phaseExist[Wpid] = OCP_TRUE;
397  xij[Wpid * numCom + Wcid] = 1.0;
398  PVTW.Eval_All(0, P, data, cdata);
399  OCP_DBL Pw0 = data[0];
400  OCP_DBL bw0 = data[1];
401  OCP_DBL cbw = data[2];
402  OCP_DBL bw = bw0 * (1 - cbw * (P - Pw0));
403  OCP_DBL bwp = -cbw * bw0;
404  mu[Wpid] = data[3];
405  xi[Wpid] = 1 / (CONV1 * bw);
406  rho[Wpid] = std_RhoW / bw;
407  muP[Wpid] = cdata[3];
408  xiP[Wpid] = -bwp / (bw * bw * CONV1);
409  rhoP[Wpid] = CONV1 * xiP[Wpid] * std_RhoW;
410  Ni[Wcid] = Vpore * Sw * xi[Wpid];
411  nj[Wpid] = Ni[Wcid];
412  Nt = Nh + Ni[Wcid];
413  vj[Wpid] = CONV1 * Ni[Wcid] * bw;
414  vfi[Wcid] = CONV1 * bw;
415  const OCP_DBL vwp = CONV1 * Ni[Wcid] * bwp;
416  vf += vj[Wpid];
417  vfP += vwp;
418  vji[numPhase - 1][numCom - 1] = vfi[Wcid];
419  vjp[numPhase - 1] = vwp;
420 
421  CalSaturation();
422  CaldXsdXpAPI03();
423  CalXiPn_partial();
424  CalRhoPn_partial();
425  CalMuPn_partial();
426  CalVfiVfp_full03();
427 
428  // Calculate pSderExist and pVnumCom
429  fill(pSderExist.begin(), pSderExist.end(), OCP_FALSE);
430  fill(pVnumCom.begin(), pVnumCom.end(), 0.0);
431  if (phaseExist[0]) {
432  pSderExist[0] = OCP_TRUE;
433  pVnumCom[0] = NC;
434  }
435  if (phaseExist[1]) {
436  pSderExist[1] = OCP_TRUE;
437  pVnumCom[1] = NC;
438  }
439  pSderExist[2] = OCP_TRUE;
440 
441  CalSkipForNextStep();
442 }
443 
444 void MixtureComp::FlashIMPEC(const OCP_DBL& Pin,
445  const OCP_DBL& Tin,
446  const OCP_DBL* Niin,
447  const USI& lastNP,
448  const OCP_DBL* xijin,
449  const OCP_USI& bId)
450 {
451  SetBulkId(bId);
452  // Hydroncarbon phase, if lNp = 0, then strict stability analysis will be used
453  lNP = lastNP > 0 ? lastNP - 1 : 0;
454  if (lNP == 2) {
455  for (USI i = 0; i < NC; i++) {
456  lKs[i] = xijin[i] / xijin[i + numCom];
457  }
458  }
459 
460  InitPTN(Pin, Tin + CONV5, Niin);
461  CalFtypeIMPEC();
462  CalFlash();
463  // Calculate derivates for hydrocarbon phase and components
464  // d vf / d Ni, d vf / d P
465  // CalVfiVfp_full01();
466  CalVfiVfp_full02();
467 
468  // Water Properties
469  const USI Wpid = numPhase - 1;
470  const USI Wcid = numCom - 1;
471  phaseExist[Wpid] = OCP_TRUE;
472  xij[Wpid * numCom + Wcid] = 1.0;
473  Nt = Nh + Ni[Wcid];
474  PVTW.Eval_All(0, P, data, cdata);
475  OCP_DBL Pw0 = data[0];
476  OCP_DBL bw0 = data[1];
477  OCP_DBL cbw = data[2];
478  OCP_DBL bw = bw0 * (1 - cbw * (P - Pw0));
479  OCP_DBL bwp = -cbw * bw0;
480  mu[Wpid] = data[3];
481  xi[Wpid] = 1 / (CONV1 * bw);
482  rho[Wpid] = std_RhoW / bw;
483  vj[Wpid] = CONV1 * Ni[Wcid] * bw;
484  vf += vj[Wpid];
485  vfi[Wcid] = CONV1 * bw;
486  vfP += CONV1 * Ni[Wcid] * bwp;
487 
488  // Calculate Sj
489  CalSaturation();
490 
491  CalSkipForNextStep();
492 }
493 
494 void MixtureComp::FlashFIM(const OCP_DBL& Pin,
495  const OCP_DBL& Tin,
496  const OCP_DBL* Niin,
497  const OCP_DBL* Sjin,
498  const USI& lastNP,
499  const OCP_DBL* xijin,
500  const OCP_USI& bId)
501 {
502 
503  SetBulkId(bId);
504 
505  // Hydroncarbon phase, if lNp = 0, then strict stability analysis will be used
506  lNP = lastNP > 0 ? lastNP - 1 : 0;
507  if (lNP == 2) {
508  for (USI i = 0; i < NC; i++) {
509  lKs[i] = xijin[i] / xijin[i + numCom];
510  }
511  }
512 
513  InitPTN(Pin, Tin + CONV5, Niin);
514  CalFtypeFIM(Sjin);
515  CalFlash();
516 
517  // Calculate derivates for hydrocarbon phase and components
518  // d vf / d Ni, d vf / d P
519  CalVfiVfp_full02();
520 
521  // Water Properties
522  USI Wpid = numPhase - 1;
523  USI Wcid = numCom - 1;
524  phaseExist[Wpid] = OCP_TRUE;
525  xij[Wpid * numCom + Wcid] = 1.0;
526  nj[Wpid] = Ni[Wcid];
527  Nt = Nh + Ni[Wcid];
528  PVTW.Eval_All(0, P, data, cdata);
529  OCP_DBL Pw0 = data[0];
530  OCP_DBL bw0 = data[1];
531  OCP_DBL cbw = data[2];
532  OCP_DBL bw = bw0 * (1 - cbw * (P - Pw0));
533  OCP_DBL bwp = -cbw * bw0;
534  mu[Wpid] = data[3];
535  xi[Wpid] = 1 / (CONV1 * bw);
536  rho[Wpid] = std_RhoW / bw;
537  muP[Wpid] = cdata[3];
538  xiP[Wpid] = -bwp / (bw * bw * CONV1);
539  rhoP[Wpid] = CONV1 * xiP[Wpid] * std_RhoW;
540  vj[Wpid] = CONV1 * Ni[Wcid] * bw;
541  vfi[Wcid] = CONV1 * bw;
542  const OCP_DBL vwp = CONV1 * Ni[Wcid] * bwp;
543  vf += vj[Wpid];
544  vfP += vwp;
545  vji[numPhase - 1][numCom - 1] = vfi[Wcid];
546  vjp[numPhase - 1] = vwp;
547 
548  CalSaturation();
549 
550 #ifdef OCP_OLD_FIM
551  CaldXsdXpAPI02();
552 #else
553  CaldXsdXpAPI02p();
554 #endif // OCP_OLD_FIM
555  CalXiPNX_partial();
556  CalRhoPX_partial();
557  CalMuPX_partial();
558 
559  // Calculate pSderExist and pVnumCom
560  fill(pSderExist.begin(), pSderExist.end(), OCP_FALSE);
561  fill(pVnumCom.begin(), pVnumCom.end(), 0);
562  if (phaseExist[0]) {
563  pSderExist[0] = OCP_TRUE;
564  pVnumCom[0] = NC;
565  }
566  if (phaseExist[1]) {
567  pSderExist[1] = OCP_TRUE;
568  pVnumCom[1] = NC;
569  }
570  pSderExist[2] = OCP_TRUE;
571 
572  CalSkipForNextStep();
573 }
574 
575 void MixtureComp::FlashFIMn(const OCP_DBL& Pin,
576  const OCP_DBL& Tin,
577  const OCP_DBL* Niin,
578  const OCP_DBL* Sjin,
579  const OCP_DBL* xijin,
580  const OCP_DBL* njin,
581  const USI* phaseExistin,
582  const USI& lastNP,
583  const OCP_USI& bId)
584 {
585  SetBulkId(bId);
586  inputNP = 0;
587  for (USI j = 0; j < numPhase; j++) {
588  if (phaseExistin[j] == 1) inputNP++;
589  }
590  inputNP--;
591 
592  if (inputNP == 1 || OCP_TRUE) {
593  // Hydroncarbon phase, if lNp = 0, then strict stability analysis will be used
594  lNP = lastNP > 0 ? lastNP - 1 : 0;
595  if (lNP == 2) {
596  for (USI i = 0; i < NC; i++) {
597  lKs[i] = xijin[i] / xijin[i + numCom];
598  }
599  }
600  InitPTN(Pin, Tin + CONV5, Niin);
601  CalFtypeFIM(Sjin);
602  CalFlash();
603  } else {
605  NP = inputNP;
606  InitPTN(Pin, Tin + CONV5, Niin);
607  CalAiBi();
608  for (USI j = 0; j < NP; j++) {
609  phaseExist[j] = phaseExistin[j];
610  S[j] = Sjin[j];
611  nu[j] = njin[j];
612  nj[j] = njin[j];
613  Dcopy(NC, &x[j][0], &xijin[j * numCom]);
614  Dcopy(NC, &xij[j * numCom], &xijin[j * numCom]);
615  }
616  CalFugPhiAll();
617  S[numPhase - 1] = Sjin[numPhase - 1];
618  phaseLabel[0] = OIL;
619  phaseLabel[1] = GAS;
620  CalMW();
621  CalVfXiRho();
622  CalViscosity();
623  CopyPhase();
624  CalSurfaceTension();
625  }
626 
627  CalVjpVfpVfn_partial();
628 
629  // Water Properties
630  USI Wpid = numPhase - 1;
631  USI Wcid = numCom - 1;
632  phaseExist[Wpid] = OCP_TRUE;
633  xij[Wpid * numCom + Wcid] = 1.0;
634  nj[Wpid] = Ni[Wcid];
635  Nt = Nh + Ni[Wcid];
636  PVTW.Eval_All(0, P, data, cdata);
637  OCP_DBL Pw0 = data[0];
638  OCP_DBL bw0 = data[1];
639  OCP_DBL cbw = data[2];
640  OCP_DBL bw = bw0 * (1 - cbw * (P - Pw0));
641  OCP_DBL bwp = -cbw * bw0;
642  mu[Wpid] = data[3];
643  xi[Wpid] = 1 / (CONV1 * bw);
644  rho[Wpid] = std_RhoW / bw;
645  muP[Wpid] = cdata[3];
646  xiP[Wpid] = -bwp / (bw * bw * CONV1);
647  rhoP[Wpid] = CONV1 * xiP[Wpid] * std_RhoW;
648  vj[Wpid] = CONV1 * Ni[Wcid] * bw;
649  vfi[Wcid] = CONV1 * bw;
650  const OCP_DBL vwp = CONV1 * Ni[Wcid] * bwp;
651  vf += vj[Wpid];
652  vfP += vwp;
653  vji[numPhase - 1][numCom - 1] = vfi[Wcid];
654  vjp[numPhase - 1] = vwp;
655 
656  CalSaturation();
657 
658  // calculate diff
659  if (inputNP == NP && NP == 2 && OCP_FALSE) {
660  cout << scientific << setprecision(3);
661  cout << "--------- " << NP << " --------- " << inputNP << endl;
662  for (USI j = 0; j < NP; j++) {
663  USI j1 = phaseLabel[j];
664  cout << setw(10) << S[j1] - Sjin[j1] << " " << setw(10)
665  << (nj[j1] - njin[j1]) / njin[j1] << " ";
666  for (USI i = 0; i < NC; i++) {
667  cout << setw(10) << xij[j1 * numCom + i] - xijin[j1 * numCom + i]
668  << " ";
669  }
670  cout << endl;
671  }
672  cout << S[2] - Sjin[2] << endl;
673  }
674 
675  CaldXsdXpAPI03();
676  CalXiPn_partial();
677  CalRhoPn_partial();
678  CalMuPn_partial();
679  CalVfiVfp_full03();
680 
681  // Calculate pSderExist and pVnumCom
682  fill(pSderExist.begin(), pSderExist.end(), OCP_FALSE);
683  fill(pVnumCom.begin(), pVnumCom.end(), 0.0);
684  if (phaseExist[0]) {
685  pSderExist[0] = OCP_TRUE;
686  pVnumCom[0] = NC;
687  }
688  if (phaseExist[1]) {
689  pSderExist[1] = OCP_TRUE;
690  pVnumCom[1] = NC;
691  }
692  pSderExist[2] = OCP_TRUE;
693 
694  CalSkipForNextStep();
695 }
696 
697 void MixtureComp::CalFlash()
698 {
699  PhaseEquilibrium();
700  // Next, nu represents moles of phase instead of molar fraction of phase
701  Dscalar(NP, Nh, &nu[0]);
702  CalMW();
703  CalVfXiRho();
704  CalViscosity();
705  CalSurfaceTension();
706  IdentifyPhase();
707  CopyPhase();
708 }
709 
710 OCP_DBL
711 MixtureComp::XiPhase(const OCP_DBL& Pin,
712  const OCP_DBL& Tin,
713  const OCP_DBL* Ziin,
714  const USI& tarPhase)
715 {
716  // assume that only single phase exists here
717  if (tarPhase == WATER) {
718  // water phase
719  PVTW.Eval_All(0, Pin, data, cdata);
720  OCP_DBL Pw0 = data[0];
721  OCP_DBL bw0 = data[1];
722  OCP_DBL cbw = data[2];
723  OCP_DBL bw = bw0 * (1 - cbw * (Pin - Pw0));
724  OCP_DBL xitmp = 1 / (CONV1 * bw);
725  return xitmp;
726  } else {
727  // hydrocarbon phase
728  InitPTZ(Pin, Tin + CONV5, Ziin); // P, T has been Set !!!
729  NP = 1;
730  CalAiBi();
731  CalAjBj(Aj[0], Bj[0], zi);
732  SolEoS(Zj[0], Aj[0], Bj[0]);
733  OCP_DBL vtmp = Zj[0] * GAS_CONSTANT * T / P;
734  for (USI i = 0; i < NC; i++) {
735  vtmp -= zi[i] * Vshift[i];
736  }
737  OCP_DBL xitmp = 1 / vtmp;
738  return xitmp;
739  }
740 }
741 
742 OCP_DBL
743 MixtureComp::RhoPhase(const OCP_DBL& Pin,
744  const OCP_DBL& Pbb,
745  const OCP_DBL& Tin,
746  const OCP_DBL* Ziin,
747  const USI& tarPhase)
748 {
749 
750  // assume that only single phase exists here
751  if (tarPhase == WATER) {
752  // water phase
753  PVTW.Eval_All(0, Pin, data, cdata);
754  OCP_DBL Pw0 = data[0];
755  OCP_DBL bw0 = data[1];
756  OCP_DBL cbw = data[2];
757  OCP_DBL bw = (bw0 * (1 - cbw * (Pin - Pw0)));
758 
759  return std_RhoW / bw;
760  } else {
761  // hydrocarbon phase
762  OCP_DBL xitmp = XiPhase(Pin, Tin, Ziin, tarPhase);
763  NP = 1;
764  x[0] = zi;
765  CalMW();
766  return MW[0] * xitmp;
767  }
768 }
769 
770 void MixtureComp::SetupWellOpt(WellOpt& opt,
771  const vector<SolventINJ>& sols,
772  const OCP_DBL& Psurf,
773  const OCP_DBL& Tsurf)
774 {
775  const USI wellType = opt.WellType();
776  if (wellType == INJ) {
777  const string fluidName = opt.InjFluidType();
778  vector<OCP_DBL> tmpZi(numCom, 0);
779  if (fluidName == "WAT") {
780  tmpZi.back() = 1;
781  opt.SetInjProdPhase(WATER);
782  } else {
783  // inj phase is gas
784  opt.SetInjProdPhase(GAS);
785  const USI len = sols.size();
786  for (USI i = 0; i < len; i++) {
787  if (fluidName == sols[i].name) {
788  tmpZi = sols[i].data;
789  tmpZi.resize(numCom);
790  // Convert volume units Mscf/stb to molar units lbmoles for
791  // injfluid Use flash in Bulk in surface condition
792  OCP_DBL tmp = 1000 * XiPhase(Psurf, Tsurf, &tmpZi[0], GAS);
793  opt.SetInjFactor(tmp);
794  break;
795  }
796  if (i == len - 1) {
797  OCP_ABORT("Wrong FluidType!");
798  }
799  }
800  }
801  opt.SetInjZi(tmpZi);
802  } else if (wellType == PROD) {
803  vector<OCP_DBL> tmpWght(numPhase, 0);
804  switch (opt.OptMode()) {
805  case BHP_MODE:
806  break;
807  case ORATE_MODE:
808  tmpWght[0] = 1;
809  break;
810  case GRATE_MODE:
811  tmpWght[1] = 1;
812  break;
813  case WRATE_MODE:
814  tmpWght[2] = 1;
815  break;
816  case LRATE_MODE:
817  tmpWght[0] = 1;
818  tmpWght[2] = 1;
819  break;
820  default:
821  OCP_ABORT("WRONG optMode");
822  break;
823  }
824  opt.SetProdPhaseWeight(tmpWght);
825  } else {
826  OCP_ABORT("Wrong Well Type!");
827  }
828 }
829 
830 void MixtureComp::CalProdWeight(const OCP_DBL& Pin,
831  const OCP_DBL& Tin,
832  const OCP_DBL* Niin,
833  const vector<OCP_DBL>& prodPhase,
834  vector<OCP_DBL>& prodWeight)
835 {
836  Flash(Pin, Tin, Niin);
837 
838  OCP_DBL qt = Nt;
839  vector<OCP_DBL> factor(numPhase, 0);
840 
841  factor[0] = vj[0] / qt / CONV1; // stb / lbmol
842  factor[1] = vj[1] / qt / 1000; // Mscf / lbmol
843  factor[2] = xi[2] * vj[2] / qt; // stb / stb
844 
845  OCP_DBL tmp = 0;
846  for (USI i = 0; i < 3; i++) {
847  tmp += factor[i] * prodPhase[i];
848  }
849  if (tmp < 1E-12 || !isfinite(tmp)) {
850  OCP_ABORT("Wrong Condition!");
851  }
852  fill(prodWeight.begin(), prodWeight.end(), tmp);
853 }
854 
855 void MixtureComp::CalProdRate(const OCP_DBL& Pin,
856  const OCP_DBL& Tin,
857  const OCP_DBL* Niin,
858  vector<OCP_DBL>& prodRate)
859 {
860  Flash(Pin, Tin, Niin);
861 
862  prodRate[0] = vj[0] / CONV1; // stb
863  prodRate[1] = vj[1] / 1000; // Mscf
864  prodRate[2] = vj[2] * xi[2]; // stb
865 }
866 
867 void MixtureComp::CallId()
868 {
869  lId = 0;
870  for (USI i = 1; i < NC; i++) {
871  if (MWC[i] < MWC[lId]) lId = i;
872  }
873 }
874 
875 void MixtureComp::AllocateEoS()
876 {
877  // Allocate Memoery for EoS variables
878  Ai.resize(NC);
879  Bi.resize(NC);
880  Aj.resize(NPmax);
881  Bj.resize(NPmax);
882  Zj.resize(NPmax);
883  Ztmp.resize(3);
884  delta1P2 = delta1 + delta2;
885  delta1M2 = delta1 - delta2;
886  delta1T2 = delta1 * delta2;
887 }
888 
889 void MixtureComp::SolEoS(OCP_DBL& ZjT, const OCP_DBL& AjT, const OCP_DBL& BjT) const
890 {
891  const OCP_DBL aj = AjT;
892  const OCP_DBL bj = BjT;
893 
894  const OCP_DBL a = (delta1 + delta2 - 1) * bj - 1;
895  const OCP_DBL b =
896  (aj + delta1 * delta2 * bj * bj - (delta1 + delta2) * bj * (bj + 1));
897  const OCP_DBL c = -(aj * bj + delta1 * delta2 * bj * bj * (bj + 1));
898 
899  USI flag = CubicRoot(a, b, c, OCP_TRUE); // True with NT
900  if (flag == 1) {
901  ZjT = Ztmp[0];
902  } else {
903  OCP_DBL zj1 = Ztmp[0];
904  OCP_DBL zj2 = Ztmp[2];
905  OCP_DBL dG = (zj2 - zj1) + log((zj1 - bj) / (zj2 - bj)) -
906  aj / (bj * (delta2 - delta1)) *
907  log((zj1 + delta1 * bj) * (zj2 + delta2 * bj) /
908  ((zj1 + delta2 * bj) * (zj2 + delta1 * bj)));
909  if (dG > 0)
910  ZjT = zj1;
911  else
912  ZjT = zj2;
913  }
914 }
915 
916 void MixtureComp::CalAiBi()
917 {
918  // Calculate Ai, Bi
919  OCP_DBL acf, mwi;
920  OCP_DBL Pri, Tri;
921 
922  for (USI i = 0; i < NC; i++) {
923  acf = Acf[i];
924  // PR
925  if (acf <= 0.49) {
926  mwi = 0.37464 + 1.54226 * acf - 0.26992 * pow(acf, 2);
927  } else {
928  mwi = 0.379642 + 1.48503 * acf - 0.164423 * pow(acf, 2) +
929  0.016667 * pow(acf, 3);
930  }
931 
932  Pri = P / Pc[i];
933  Tri = T / Tc[i];
934  Ai[i] = OmegaA[i] * Pri / pow(Tri, 2) * pow((1 + mwi * (1 - sqrt(Tri))), 2);
935  Bi[i] = OmegaB[i] * Pri / Tri;
936  }
937 }
938 
939 void MixtureComp::CalAjBj(OCP_DBL& AjT, OCP_DBL& BjT, const vector<OCP_DBL>& xj) const
940 {
941  AjT = 0;
942  BjT = 0;
943 
944  for (USI i1 = 0; i1 < NC; i1++) {
945  BjT += Bi[i1] * xj[i1];
946  AjT += xj[i1] * xj[i1] * Ai[i1] * (1 - BIC[i1 * NC + i1]);
947 
948  for (USI i2 = 0; i2 < i1; i2++) {
949  AjT +=
950  2 * xj[i1] * xj[i2] * sqrt(Ai[i1] * Ai[i2]) * (1 - BIC[i1 * NC + i2]);
951  }
952  }
953 }
954 
955 void MixtureComp::CalAjBj(OCP_DBL& AjT, OCP_DBL& BjT, const OCP_DBL* xj) const
956 {
957  AjT = 0;
958  BjT = 0;
959 
960  for (USI i1 = 0; i1 < NC; i1++) {
961  BjT += Bi[i1] * xj[i1];
962  AjT += xj[i1] * xj[i1] * Ai[i1] * (1 - BIC[i1 * NC + i1]);
963 
964  for (USI i2 = 0; i2 < i1; i2++) {
965  AjT +=
966  2 * xj[i1] * xj[i2] * sqrt(Ai[i1] * Ai[i2]) * (1 - BIC[i1 * NC + i2]);
967  }
968  }
969 }
970 
971 void MixtureComp::AllocatePhase()
972 {
973  // Allocate Memoery for Phase variables
974  vC.resize(NPmax);
975  nu.resize(NPmax);
976  xiC.resize(NPmax);
977  rhoC.resize(NPmax);
978  MW.resize(NPmax);
979  phaseLabel.resize(NPmax);
980  x.resize(NPmax);
981  phi.resize(NPmax);
982  fug.resize(NPmax);
983  n.resize(NPmax);
984  for (USI j = 0; j < NPmax; j++) {
985  x[j].resize(NC);
986  phi[j].resize(NC);
987  fug[j].resize(NC);
988  n[j].resize(NC);
989  }
990  ln = n;
991 }
992 
993 void MixtureComp::CalFugPhi(vector<OCP_DBL>& phiT,
994  vector<OCP_DBL>& fugT,
995  const vector<OCP_DBL>& xj)
996 {
997  OCP_DBL aj, bj, zj;
998  CalAjBj(aj, bj, xj);
999  SolEoS(zj, aj, bj);
1000 
1001  const OCP_DBL m1 = delta1;
1002  const OCP_DBL m2 = delta2;
1003  // const OCP_DBL m1Mm2 = delta1M2;
1004 
1005  OCP_DBL tmp;
1006  for (USI i = 0; i < NC; i++) {
1007  tmp = 0;
1008  for (USI k = 0; k < NC; k++) {
1009  tmp += 2 * (1 - BIC[i * NC + k]) * sqrt(Ai[i] * Ai[k]) * xj[k];
1010  }
1011  phiT[i] = exp(Bi[i] / bj * (zj - 1) - log(zj - bj) -
1012  aj / (m1 - m2) / bj * (tmp / aj - Bi[i] / bj) *
1013  log((zj + m1 * bj) / (zj + m2 * bj)));
1014  fugT[i] = phiT[i] * xj[i] * P;
1015  }
1016 
1017  // OCP_DBL tmp01 = log(zj - bj);
1018  // OCP_DBL tmp02 = aj / (m1Mm2) / bj * log((zj + m1 * bj) / (zj + m2 * bj));
1019 
1020  // for (USI i = 0; i < NC; i++) {
1021  // tmp = 0;
1022  // for (int k = 0; k < NC; k++) {
1023  // tmp += 2 * (1 - BIC[i * NC + k]) * sqrt(Ai[i] * Ai[k]) * xj[k];
1024  // }
1025  // phiT[i] = exp(Bi[i] / bj * (zj - 1) - tmp01 - tmp02 * (tmp / aj - Bi[i] /
1026  // bj));
1027  // fugT[i] = phiT[i] * xj[i] * P;
1028  //}
1029 
1030  Asta = aj;
1031  Bsta = bj;
1032  Zsta = zj;
1033 }
1034 
1035 void MixtureComp::CalFugPhi(OCP_DBL* phiT, OCP_DBL* fugT, const OCP_DBL* xj)
1036 {
1037  OCP_DBL aj, bj, zj;
1038  CalAjBj(aj, bj, xj);
1039  SolEoS(zj, aj, bj);
1040 
1041  const OCP_DBL m1 = delta1;
1042  const OCP_DBL m2 = delta2;
1043  const OCP_DBL m1Mm2 = delta1M2;
1044 
1045  // OCP_DBL tmp01 = log(zj - bj);
1046  // OCP_DBL tmp02 = aj / (m1Mm2) / bj * log((zj + m1 * bj) / (zj + m2 * bj));
1047 
1048  OCP_DBL tmp;
1049  for (USI i = 0; i < NC; i++) {
1050  tmp = 0;
1051  for (USI k = 0; k < NC; k++) {
1052  tmp += 2 * (1 - BIC[i * NC + k]) * sqrt(Ai[i] * Ai[k]) * xj[k];
1053  }
1054  phiT[i] = exp(Bi[i] / bj * (zj - 1) - log(zj - bj) -
1055  aj / (m1Mm2) / bj * (tmp / aj - Bi[i] / bj) *
1056  log((zj + m1 * bj) / (zj + m2 * bj)));
1057  fugT[i] = phiT[i] * xj[i] * P;
1058  }
1059  // for (USI i = 0; i < NC; i++) {
1060  // tmp = 0;
1061  // for (int k = 0; k < NC; k++) {
1062  // tmp += 2 * (1 - BIC[i * NC + k]) * sqrt(Ai[i] * Ai[k]) * xj[k];
1063  // }
1064  // phiT[i] = exp(Bi[i] / bj * (zj - 1) - tmp01 - tmp02 * (tmp / aj - Bi[i] /
1065  // bj)); fugT[i] = phiT[i] * xj[i] * P;
1066  //}
1067 
1068  Asta = aj;
1069  Bsta = bj;
1070  Zsta = zj;
1071 }
1072 
1073 void MixtureComp::CalFugPhi(OCP_DBL* fugT, const OCP_DBL* xj)
1074 {
1075  OCP_DBL aj, bj, zj;
1076  CalAjBj(aj, bj, xj);
1077  SolEoS(zj, aj, bj);
1078 
1079  const OCP_DBL m1 = delta1;
1080  const OCP_DBL m2 = delta2;
1081  const OCP_DBL m1Mm2 = delta1M2;
1082 
1083  // OCP_DBL tmp01 = log(zj - bj);
1084  // OCP_DBL tmp02 = aj / (m1Mm2) / bj * log((zj + m1 * bj) / (zj + m2 * bj));
1085 
1086  OCP_DBL tmp;
1087  for (USI i = 0; i < NC; i++) {
1088  tmp = 0;
1089  for (USI k = 0; k < NC; k++) {
1090  tmp += 2 * (1 - BIC[i * NC + k]) * sqrt(Ai[i] * Ai[k]) * xj[k];
1091  }
1092  tmp = exp(Bi[i] / bj * (zj - 1) - log(zj - bj) -
1093  aj / (m1Mm2) / bj * (tmp / aj - Bi[i] / bj) *
1094  log((zj + m1 * bj) / (zj + m2 * bj)));
1095  fugT[i] = tmp * xj[i] * P;
1096  }
1097 
1098  Asta = aj;
1099  Bsta = bj;
1100  Zsta = zj;
1101 }
1102 
1103 void MixtureComp::CalFugPhiAll()
1104 {
1105  OCP_DBL tmp; // , tmp01, tmp02;
1106  const OCP_DBL m1 = delta1;
1107  const OCP_DBL m2 = delta2;
1108  const OCP_DBL m1Mm2 = delta1M2;
1109 
1110  for (USI j = 0; j < NP; j++) {
1111  const vector<OCP_DBL>& xj = x[j];
1112  vector<OCP_DBL>& phiT = phi[j];
1113  vector<OCP_DBL>& fugT = fug[j];
1114  OCP_DBL& aj = Aj[j];
1115  OCP_DBL& bj = Bj[j];
1116  OCP_DBL& zj = Zj[j];
1117 
1118  CalAjBj(aj, bj, xj);
1119  SolEoS(zj, aj, bj);
1120 
1121  for (USI i = 0; i < NC; i++) {
1122  tmp = 0;
1123  for (USI k = 0; k < NC; k++) {
1124  tmp += 2 * (1 - BIC[i * NC + k]) * sqrt(Ai[i] * Ai[k]) * xj[k];
1125  }
1126  phiT[i] = exp(Bi[i] / bj * (zj - 1) - log(zj - bj) -
1127  aj / (m1Mm2) / bj * (tmp / aj - Bi[i] / bj) *
1128  log((zj + m1 * bj) / (zj + m2 * bj)));
1129  fugT[i] = phiT[i] * xj[i] * P;
1130  }
1131 
1132  // tmp01 = log(zj - bj);
1133  // tmp02 = aj / (m1Mm2) / bj * log((zj + m1 * bj) / (zj + m2 * bj));
1134  // for (USI i = 0; i < NC; i++) {
1135  // tmp = 0;
1136  // for (int k = 0; k < NC; k++) {
1137  // tmp += 2 * (1 - BIC[i * NC + k]) * sqrt(Ai[i] * Ai[k]) * xj[k];
1138  // }
1139  // phiT[i] =
1140  // exp(Bi[i] / bj * (zj - 1) - tmp01 - tmp02 * (tmp / aj - Bi[i] /
1141  // bj));
1142  // fugT[i] = phiT[i] * xj[i] * P;
1143  //}
1144  }
1145 }
1146 
1147 void MixtureComp::CalMW()
1148 {
1149  // Calculate Molecular Weight of phase
1150  for (USI j = 0; j < NP; j++) {
1151  MW[j] = 0;
1152  for (USI i = 0; i < NC; i++) {
1153  MW[j] += x[j][i] * MWC[i];
1154  }
1155  }
1156 }
1157 
1158 void MixtureComp::CalVfXiRho()
1159 {
1160  // Attention that nu is moles of phase instead of molar fraction of phase now
1161  vf = 0;
1162  OCP_DBL tmp;
1163  for (USI j = 0; j < NP; j++) {
1164 
1165  vector<OCP_DBL>& xj = x[j];
1166  tmp = Zj[j] * GAS_CONSTANT * T / P;
1167  for (USI i = 0; i < NC; i++) {
1168  tmp -= xj[i] * Vshift[i];
1169  }
1170 
1171  vC[j] = tmp * nu[j];
1172  vf += vC[j];
1173  xiC[j] = 1 / tmp;
1174  rhoC[j] = MW[j] * xiC[j];
1175  }
1176 }
1177 
1178 void MixtureComp::CalSaturation()
1179 {
1180  for (USI j = 0; j < numPhase; j++) {
1181  S[j] = 0;
1182  if (phaseExist[j]) {
1183  S[j] = vj[j] / vf;
1184  }
1185  }
1186 }
1187 
1188 USI MixtureComp::FindMWmax()
1189 {
1190  // find the phase id with the highest Molecular Weight
1191  USI tmpId = 0;
1192  OCP_DBL tmpMW = MW[0];
1193  for (USI j = 1; j < NP; j++) {
1194  if (tmpMW < MW[j]) {
1195  tmpMW = MW[j];
1196  tmpId = j;
1197  }
1198  }
1199  return tmpId;
1200 }
1201 
1202 void MixtureComp::x2n()
1203 {
1204  // Total moles are supposed to be 1
1205  for (USI j = 0; j < NP; j++) {
1206 
1207  // nu[j] = fabs(nu[j]);
1208  for (USI i = 0; i < NC; i++) {
1209  n[j][i] = nu[j] * x[j][i];
1210  }
1211  }
1212 }
1213 
1214 void MixtureComp::PrintX()
1215 {
1216  for (USI j = 0; j < NP; j++) {
1217  for (USI i = 0; i < NC; i++) {
1218  cout << x[j][i] << " ";
1219  }
1220  cout << endl;
1221  }
1222  cout << "----------------------" << endl;
1223 }
1224 
1225 void MixtureComp::AllocateMethod()
1226 {
1227 
1228  Kw.resize(4);
1229  for (USI i = 0; i < 4; i++) {
1230  Kw[i].resize(NC);
1231  }
1232  Ks.resize(NPmax - 1);
1233  for (USI i = 0; i < NPmax - 1; i++) {
1234  Ks[i].resize(NC);
1235  }
1236  phiSta.resize(NC);
1237  fugSta.resize(NC);
1238 
1239  Y.resize(NC);
1240  di.resize(NC);
1241  resSTA.resize(NC);
1242 
1243  JmatSTA.resize(NC * NC);
1244  Ax.resize(NC);
1245  Bx.resize(NC);
1246  Zx.resize(NC);
1247  lKs.resize(NC);
1248 
1249  resRR.resize(NPmax - 1);
1250  resSP.resize(static_cast<size_t>(NC) * NPmax);
1251  JmatSP.resize(static_cast<size_t>(NC * NC) * NPmax * NPmax);
1252  fugX.resize(NPmax);
1253  fugN.resize(NPmax);
1254  Zn.resize(NPmax);
1255  for (USI j = 0; j < NPmax; j++) {
1256  fugX[j].resize(NC * NC);
1257  fugN[j].resize(NC * NC);
1258  Zn[j].resize(NC);
1259  }
1260  An.resize(NC);
1261  Bn.resize(NC);
1262  pivot.resize(static_cast<size_t>(NC + 1) * NPmax, 1);
1263  lJmatWork = NC * (NPmax - 1);
1264  JmatWork.resize(lJmatWork);
1265 
1266  pivot.resize(NPmax * static_cast<size_t>(NC) + numPhase, 1);
1267 }
1268 
1269 void MixtureComp::CalKwilson()
1270 {
1271  for (USI i = 0; i < NC; i++) {
1272  Kw[0][i] = (Pc[i] / P) * exp(5.373 * (1 + Acf[i]) * (1 - Tc[i] / T));
1273  Kw[1][i] = 1 / Kw[0][i];
1274  Kw[2][i] = pow(Kw[0][i], 1.0 / 3);
1275  Kw[3][i] = pow(Kw[1][i], 1.0 / 3);
1276  }
1277 }
1278 
1279 void MixtureComp::PhaseEquilibrium()
1280 {
1281  // Attention: sum of components' moles equals 1
1282  switch (ftype) {
1283  case 0:
1284  // flash from single phase
1285  flagSkip = OCP_TRUE;
1286  NP = 1;
1287  nu[0] = 1;
1288  x[0] = zi;
1289  CalAiBi();
1290  CalAjBj(Aj[0], Bj[0], x[0]);
1291  SolEoS(Zj[0], Aj[0], Bj[0]);
1292  CalKwilson();
1293  while (!PhaseStable()) {
1294  NP++;
1295  PhaseSplit();
1296  if (NP == NPmax || NP == 1) break;
1297  }
1298  if (NP > 1) {
1299  flagSkip = OCP_FALSE;
1300  }
1301  // record error
1302  if (NP == 1) {
1303  if (EoSctrl.NRsta.conflag)
1304  ePEC = EoSctrl.NRsta.realTol;
1305  else if (EoSctrl.SSMsta.conflag)
1306  ePEC = EoSctrl.SSMsta.realTol;
1307  else
1308  ePEC = 1E8;
1309  } else {
1310  if (EoSctrl.NRsp.conflag)
1311  ePEC = EoSctrl.NRsp.realTol;
1312  else if (EoSctrl.SSMsp.conflag)
1313  ePEC = EoSctrl.SSMsp.realTol;
1314  else
1315  ePEC = 1E8;
1316  }
1317 
1318  break;
1319  case 1:
1320  // Skip Phase Stability analysis, only single phase exists
1321  flagSkip = OCP_TRUE;
1322  NP = 1;
1323  nu[0] = 1;
1324  x[0] = zi;
1325  CalAiBi();
1326  CalAjBj(Aj[0], Bj[0], x[0]);
1327  SolEoS(Zj[0], Aj[0], Bj[0]);
1328  // record error
1329  ePEC = 0.0;
1330  break;
1331 
1332  case 2:
1333  // Skip Phase Stability analysis, two phases exist
1334  flagSkip = OCP_FALSE;
1335  NP = 2;
1336  Yt = 1.01;
1337  CalAiBi();
1338  CalKwilson();
1339  PhaseSplit();
1340 
1341  if (EoSctrl.NRsp.conflag)
1342  ePEC = EoSctrl.NRsp.realTol;
1343  else if (EoSctrl.SSMsp.conflag)
1344  ePEC = EoSctrl.SSMsp.realTol;
1345  else
1346  ePEC = 1E8;
1347  break;
1348 
1349  default:
1350  OCP_ABORT("Wrong flash type!");
1351  break;
1352  }
1353 }
1354 
1355 OCP_BOOL MixtureComp::PhaseStable()
1356 {
1357  if (NP == 1) {
1358  testPId = 0;
1359  } else {
1360  CalMW();
1361  testPId = FindMWmax();
1362  }
1363 
1364  EoSctrl.SSMsta.conflag = OCP_FALSE;
1365  EoSctrl.SSMsta.curIt = 0;
1366  EoSctrl.NRsta.conflag = OCP_FALSE;
1367  EoSctrl.NRsta.curIt = 0;
1368 
1369  // Test if a phase is stable, if stable return OCP_TRUE, else return OCP_FALSE
1370  OCP_BOOL flag;
1371  USI tmpNP = NP;
1372 
1373  if (lNP == 0) {
1374  // strict stability ananlysis
1375  flag = StableSSM(testPId);
1376  } else {
1377  flag = StableSSM01(testPId);
1378  if (!flag) tmpNP++;
1379 
1380  if (tmpNP != lNP) {
1381  flag = StableSSM(testPId);
1382  flagSkip = OCP_FALSE;
1383  }
1384  }
1385  itersSSMSTA += EoSctrl.SSMsta.curIt;
1386  itersNRSTA += EoSctrl.NRsta.curIt;
1387  countsSSMSTA++;
1388  countsNRSTA++;
1389  return flag;
1390 }
1391 
1392 OCP_BOOL MixtureComp::StableSSM01(const USI& Id)
1393 {
1394  // if unsatble, return OCP_FALSE
1395  // if stable, return OCP_TRUE
1396 
1397  OCP_DBL Stol = EoSctrl.SSMsta.tol2;
1398  USI maxIt = EoSctrl.SSMsta.maxIt;
1399  OCP_DBL eYt = EoSctrl.SSMsta.eYt;
1400  OCP_DBL Ktol = EoSctrl.SSMsta.Ktol;
1401  OCP_DBL dYtol = EoSctrl.SSMsta.dYtol;
1402  // OCP_DBL& Sk = EoSctrl.SSMsta.curSk;
1403  OCP_DBL Se, Sk, dY;
1404 
1405  OCP_BOOL flag, Tsol; // Tsol, trivial solution
1406  USI iter, k;
1407 
1408  const vector<OCP_DBL>& xj = x[Id];
1409  CalFugPhi(phi[Id], fug[Id], xj);
1410  const vector<OCP_DBL>& fugId = fug[Id];
1411 
1412  vector<OCP_DBL>& ks = Ks[0];
1413  for (k = 0; k < 2; k++) {
1414 
1415  ks = Kw[k];
1416  iter = 0;
1417  flag = OCP_FALSE;
1418  Tsol = OCP_FALSE;
1419  while (OCP_TRUE) {
1420  Yt = 0;
1421  for (USI i = 0; i < NC; i++) {
1422  Y[i] = xj[i] * ks[i];
1423  Yt += Y[i];
1424  }
1425  Dscalar(NC, 1 / Yt, &Y[0]);
1426 
1427  if (iter > 0) {
1428  dY = 0;
1429  for (USI i = 0; i < NC; i++) {
1430  dY += pow((Y[i] - di[i]), 2);
1431  }
1432  if (dY < dYtol) {
1433  // converges
1434  flag = OCP_TRUE;
1435  break;
1436  }
1437  }
1438 
1439  CalFugPhi(&fugSta[0], &Y[0]);
1440  Se = 0;
1441  Sk = 0;
1442  for (USI i = 0; i < NC; i++) {
1443  di[i] = fugId[i] / (fugSta[i] * Yt);
1444  ks[i] *= di[i];
1445  Se += pow(di[i] - 1, 2);
1446  // Sk += pow(ks[i] - 1, 2);
1447  Sk += pow(log(ks[i]), 2);
1448  }
1449 
1450  iter++;
1451  if (Se < Stol) {
1452  flag = OCP_TRUE;
1453  break;
1454  }
1455  if (Sk < Ktol) {
1456  // Sk < Ktol -> trivial solution
1457  flag = OCP_TRUE;
1458  Tsol = OCP_TRUE;
1459  break;
1460  }
1461  if (iter > maxIt) {
1462  break;
1463  }
1464 
1465  // Record last Y with di
1466  di = Y;
1467  }
1468 
1469  if (!Tsol) {
1470  // flag = StableNR(Id);
1471  }
1472 
1473  // cout << "Yt = " << setprecision(8) << scientific << Yt << " " << setw(2)
1474  // << "Sk = " << setprecision(3) << scientific << Sk << " " << setw(2)
1475  // << iter << " ";
1476  EoSctrl.SSMsta.curIt += iter;
1477 
1478  if (flag && Yt > 1 - 0.1 && Sk > 1) {
1479  // close to phase boundary, or more than 1 phase, So don't skip at next step
1480  flagSkip = OCP_FALSE;
1481  }
1482  if (flag && Yt > 1 + eYt) {
1483  EoSctrl.SSMsta.conflag = OCP_TRUE;
1484  EoSctrl.SSMsta.realTol = sqrt(Se);
1485  return OCP_FALSE;
1486  }
1487  }
1488  EoSctrl.SSMsta.realTol = sqrt(Se);
1489  return OCP_TRUE;
1490 }
1491 
1492 OCP_BOOL MixtureComp::StableSSM(const USI& Id)
1493 {
1494  // if unsatble, return OCP_FALSE
1495  // if stable, return OCP_TRUE
1496 
1497  const vector<OCP_DBL>& xj = x[Id];
1498  CalFugPhi(phi[Id], fug[Id], xj);
1499  const vector<OCP_DBL>& fugId = fug[Id];
1500 
1501  for (USI i = 0; i < NC; i++) {
1502  di[i] = phi[Id][i] * xj[i];
1503  }
1504 
1505  OCP_DBL Stol = EoSctrl.SSMsta.tol2;
1506  USI maxIt = EoSctrl.SSMsta.maxIt;
1507  OCP_DBL eYt = EoSctrl.SSMsta.eYt;
1508  OCP_DBL Se;
1509  OCP_BOOL flag;
1510  USI iter;
1511  USI k;
1512 
1513  // cout << "SSMBEGINS" << endl;
1514 
1515  for (k = 0; k < 4; k++) {
1516 
1517  Yt = 0;
1518  for (USI i = 0; i < NC; i++) {
1519  Y[i] = xj[i] / Kw[k][i];
1520  Yt += Y[i];
1521  }
1522  Dscalar(NC, 1 / Yt, &Y[0]);
1523  // CalFugPhi(phiSta, fugSta, Y);
1524  CalFugPhi(&phiSta[0], &fugSta[0], &Y[0]);
1525 
1526  Se = 0;
1527  for (USI i = 0; i < NC; i++) {
1528  Se += pow(log(fugSta[i] / fugId[i] * Yt), 2);
1529  }
1530 
1531  flag = OCP_TRUE;
1532  iter = 0;
1533 
1534  // cout << "ssmbegins" << endl;
1535 
1536  while (Se > Stol) {
1537 
1538  // cout << setprecision(6) << scientific << Se;
1539  // cout << " ";
1540  // cout << setprecision(12) << scientific << Yt;
1541  // cout << " " << iter << endl;
1542  // PrintDX(NC, &xj[0]);
1543  // PrintDX(NC, &Y[0]);
1544 
1545  Yt = 0;
1546  for (USI i = 0; i < NC; i++) {
1547  Y[i] = di[i] / phiSta[i];
1548  Yt += Y[i];
1549  }
1550  Dscalar(NC, 1 / Yt, &Y[0]);
1551  // CalFugPhi(phiSta, fugSta, Y);
1552  CalFugPhi(&phiSta[0], &fugSta[0], &Y[0]);
1553  Se = 0;
1554  for (USI i = 0; i < NC; i++) {
1555  Se += pow(log(fugSta[i] / fugId[i] * Yt), 2);
1556  }
1557 
1558  iter++;
1559  if (iter > maxIt) {
1560  flag = OCP_FALSE;
1561  break;
1562  }
1563  }
1564  EoSctrl.SSMsta.curIt += iter;
1565  flag = StableNR(Id);
1566  // here, a relaxation is needed, on the one hand it can prevent the influence
1567  // of rounding error, on the other hand, if Yt is too close to 1, phase
1568  // splitting calculation may get into trouble and single phase is indentified
1569  // finally
1570  if (flag && Yt > 1 + eYt) {
1571  // cout << "Yt = " << setprecision(12) << scientific << Yt << " " <<
1572  // setw(3)
1573  // << EoSctrl.SSMsta.curIt << " " << setw(3) << EoSctrl.NRsta.curIt
1574  // << " " << flag << " " << k << " " << 2 << " ";
1575  EoSctrl.SSMsta.conflag = OCP_TRUE;
1576  EoSctrl.SSMsta.realTol = sqrt(Se);
1577  return OCP_FALSE;
1578  }
1579  }
1580  // cout << "Yt = " << setprecision(12) << scientific << Yt << " " << setw(3)
1581  // << EoSctrl.SSMsta.curIt << " " << setw(3) << EoSctrl.NRsta.curIt << " "
1582  // << flag << " " << k << " " << 1 << " ";
1583  /*if (!flag) {
1584  OCP_WARNING("SSM not converged in Stability Analysis");
1585  }*/
1586  EoSctrl.SSMsta.realTol = sqrt(Se);
1587  return OCP_TRUE;
1588 }
1589 
1590 OCP_BOOL MixtureComp::StableNR(const USI& Id)
1591 {
1592 
1593  for (USI i = 0; i < NC; i++) {
1594  resSTA[i] = log(fug[Id][i] / (fugSta[i] * Yt));
1595  }
1596 
1597  USI maxIt = EoSctrl.NRsta.maxIt;
1598  OCP_DBL Stol = EoSctrl.NRsta.tol;
1599  OCP_DBL Se = Dnorm2(NC, &resSTA[0]);
1600  OCP_DBL alpha = 1;
1601  USI iter = 0;
1602  // OCP_DBL Se0 = Se;
1603 
1604  while (Se > Stol) {
1605 
1606  CalFugXSTA();
1607  AssembleJmatSTA();
1608  // LUSolve(1, NC, &JmatSTA[0], &resSTA[0], &pivot[0]);
1609  SYSSolve(1, &uplo, NC, &JmatSTA[0], &resSTA[0], &pivot[0], &JmatWork[0],
1610  lJmatWork);
1611  Dscalar(NC, Yt, &Y[0]);
1612  Daxpy(NC, alpha, &resSTA[0], &Y[0]);
1613  Yt = 0;
1614  for (USI i = 0; i < NC; i++) {
1615  Yt += Y[i];
1616  }
1617  Dscalar(NC, 1 / Yt, &Y[0]);
1618 
1619  CalFugPhi(&fugSta[0], &Y[0]);
1620  for (USI i = 0; i < NC; i++) {
1621  resSTA[i] = log(fug[Id][i] / (fugSta[i] * Yt));
1622  }
1623  Se = Dnorm2(NC, &resSTA[0]);
1624  iter++;
1625  if (iter > maxIt) {
1626  EoSctrl.NRsta.curIt += iter;
1627  EoSctrl.NRsta.conflag = OCP_FALSE;
1628  EoSctrl.NRsta.realTol = Se;
1629  return OCP_FALSE;
1630  }
1631  }
1632  EoSctrl.NRsta.curIt += iter;
1633  EoSctrl.NRsta.conflag = OCP_TRUE;
1634  EoSctrl.NRsta.realTol = Se;
1635  return OCP_TRUE;
1636 }
1637 
1638 void MixtureComp::CalFugXSTA()
1639 {
1640  // Y sums to be 1 now, it's actually the mole fraction of spliting phase
1641  // for stability analysis
1642  vector<OCP_DBL>& fugx = fugX[0];
1643  OCP_DBL aj = Asta;
1644  OCP_DBL bj = Bsta;
1645  OCP_DBL zj = Zsta;
1646  OCP_DBL tmp = 0;
1647 
1648  const OCP_DBL m1 = delta1;
1649  const OCP_DBL m2 = delta2;
1650  const OCP_DBL m1Pm2 = delta1P2;
1651  const OCP_DBL m1Mm2 = delta1M2;
1652  const OCP_DBL m1Tm2 = delta1T2;
1653 
1654  Bx = Bi;
1655  for (USI i = 0; i < NC; i++) {
1656  tmp = 0;
1657  for (USI k = 0; k < NC; k++) {
1658  tmp += Y[k] * (1 - BIC[i * NC + k]) * sqrt(Ai[i] * Ai[k]);
1659  }
1660  Ax[i] = 2 * tmp;
1661  Zx[i] = ((bj - zj) * Ax[i] + ((aj + m1Tm2 * (3 * bj * bj + 2 * bj)) +
1662  ((m1Pm2) * (2 * bj + 1) - 2 * m1Tm2 * bj) * zj -
1663  (m1Pm2 - 1) * zj * zj) *
1664  Bx[i]) /
1665  (3 * zj * zj + 2 * ((m1Pm2 - 1) * bj - 1) * zj +
1666  (aj + m1Tm2 * bj * bj - (m1Pm2)*bj * (bj + 1)));
1667  }
1668 
1669  OCP_DBL C, E, G;
1670  OCP_DBL Cxk, Dxk, Exk, Gxk;
1671  OCP_DBL aik;
1672  G = (zj + m1 * bj) / (zj + m2 * bj);
1673 
1674  for (USI i = 0; i < NC; i++) {
1675 
1676  C = Y[i] * P / (zj - bj);
1677  // C = 1 / (zj - bj);
1678  // D = Bx[i] * (zj - 1) / bj;
1679  E = -aj / ((m1Mm2)*bj) * (Ax[i] / aj - Bx[i] / bj);
1680 
1681  for (USI k = 0; k < NC; k++) {
1682  aik = (1 - BIC[i * NC + k]) * sqrt(Ai[i] * Ai[k]);
1683 
1684  // Cxk = -Y[i] * (Zx[k] - Bx[k]) / ((zj - bj) * (zj - bj));
1685  Cxk = ((zj - bj) * delta(i, k) - Y[i] * (Zx[k] - Bx[k])) * P /
1686  ((zj - bj) * (zj - bj));
1687  Dxk = Bx[i] / bj * (Zx[k] - Bx[k] * (zj - 1) / bj);
1688  /*Exk = (Ax[k] * bj - aj * Bx[k]) / (bj * bj) * (Ax[i] / aj - Bx[i] / bj) +
1689  aj / bj * (2 * (1 - BIC[i * NC + k]) * sqrt(Ai[i] * Ai[k]) / aj -
1690  Ax[k] * Ax[i] / (aj * aj) + Bx[i] * Bx[k] / (bj * bj));*/
1691  Exk = (2 * (aj / bj * Bx[k] * Bx[i] + bj * aik) - Ax[i] * Bx[k] -
1692  Ax[k] * Bi[i]) /
1693  (bj * bj);
1694  Exk /= -(m1Mm2);
1695  Gxk = (m1Mm2) / (zj + m2 * bj) / (zj + m2 * bj) * (zj * Bx[k] - Zx[k] * bj);
1696  fugx[i * NC + k] = 1 / C * Cxk + Dxk + Exk * log(G) + E / G * Gxk;
1697  }
1698  }
1699  // cout << "PhiX" << endl;
1700  // for (USI i = 0; i < NC; i++) {
1701  // for (USI k = 0; k < NC; k++) {
1702  // cout << fugx[i * NC + k] << " ";
1703  // }
1704  // cout << endl;
1705  //}
1706 
1707 #ifdef OCP_NANCHECK
1708  for (USI j = 0; j < NP; j++) {
1709  if (!CheckNan(fugX[j].size(), &fugX[j][0])) {
1710  OCP_ABORT("INF or NAN in fugX !");
1711  }
1712  }
1713 #endif // NANCHECK
1714 }
1715 
1716 void MixtureComp::AssembleJmatSTA()
1717 {
1718  vector<OCP_DBL>& fugx = fugX[0];
1719  fill(JmatSTA.begin(), JmatSTA.end(), 0.0);
1720  OCP_DBL tmp;
1721  for (USI i = 0; i < NC; i++) {
1722 
1723  tmp = 0;
1724  for (USI k = 0; k < NC; k++) {
1725  tmp += Y[k] * fugx[i * NC + k];
1726  }
1727 
1728  for (USI j = 0; j < NC; j++) {
1729  // Symmetric
1730  // JmatSTA[i * NC + j] = (fugx[i * NC + j] - tmp + delta(i, j) / Y[i]) / Yt;
1731  JmatSTA[i * NC + j] = (fugx[i * NC + j] - tmp + 1) / Yt;
1732  }
1733  }
1734  // cout << "JmatSTA" << endl;
1735  // for (USI i = 0; i < NC; i++) {
1736  // for (USI k = 0; k < NC; k++) {
1737  // cout << JmatSTA[k * NC + i] << " ";
1738  // }
1739  // cout << endl;
1740  //}
1741 }
1742 
1743 void MixtureComp::PhaseSplit()
1744 {
1745  EoSctrl.SSMsp.conflag = OCP_FALSE;
1746  EoSctrl.SSMsp.curIt = 0;
1747  EoSctrl.NRsp.conflag = OCP_FALSE;
1748  EoSctrl.NRsp.curIt = 0;
1749  EoSctrl.RR.curIt = 0;
1750 
1751  // cout << "begin" << endl;
1752  SplitSSM(OCP_FALSE);
1753  SplitNR();
1754  while (!EoSctrl.NRsp.conflag) {
1755  SplitSSM(OCP_TRUE);
1756  SplitNR();
1757  if (!CheckSplit()) break;
1758  if (EoSctrl.SSMsp.conflag) break;
1759  }
1760  CheckSplit();
1761 
1762  itersSSMSP += EoSctrl.SSMsp.curIt;
1763  itersNRSP += EoSctrl.NRsp.curIt;
1764  itersRR += EoSctrl.RR.curIt;
1765  countsSSMSP++;
1766  countsNRSP++;
1767  countsRR++;
1768 
1769  // cout << scientific << setprecision(8);
1770  // for (USI i = 0; i < NC; i++) {
1771  // cout << x[0][i] / x[1][i] << " ";
1772  // }
1773  // cout << endl;
1774  // cout << "Yt = " << scientific << setprecision(12) << Yt << " "
1775  // << setw(3) << "SSMtol = " << setprecision(6) << sqrt(EoSctrl.SSMsp.realTol)
1776  // << " "
1777  // << setw(3) << EoSctrl.SSMsp.curIt << " "
1778  // << setw(3) << "NRtol = " << setprecision(6) << EoSctrl.NRsp.realTol << " "
1779  // << setw(3) << EoSctrl.NRsp.curIt << " "
1780  // << setw(2) << lNP << " " << setw(2) << NP << " "
1781  // << (lNP == NP ? "N" : "Y") << " ";
1782 }
1783 
1784 OCP_BOOL MixtureComp::CheckSplit()
1785 {
1786  if (NP == 2) {
1787 
1788  OCP_DBL eX = 0;
1789  OCP_DBL nuMax = max(nu[0], nu[1]);
1790 
1791  for (USI i = 0; i < NC; i++) {
1792  eX += (x[0][i] - x[1][i]) * (x[0][i] - x[1][i]);
1793  }
1794 
1795  if (OCP_TRUE) {
1796  // Calculate Gibbs Energy
1797 
1798  CalFugPhi(phiSta, fugSta, zi);
1799  GibbsEnergyB = 0;
1800  GibbsEnergyE = 0;
1801  for (USI i = 0; i < NC; i++) {
1802  GibbsEnergyB += zi[i] * log(fugSta[i]);
1803  GibbsEnergyE += (n[0][i] * log(fug[0][i]) + n[1][i] * log(fug[1][i]));
1804  }
1805 
1806  cout << scientific << setprecision(6);
1807  // cout << GibbsEnergyB << " " << GibbsEnergyE << endl;
1808  if (GibbsEnergyE > GibbsEnergyB) {
1809  cout << ftype << " ";
1810  cout << GibbsEnergyB << " ";
1811  cout << GibbsEnergyE << " ";
1812  cout << nuMax << " ";
1813  cout << eX << " ";
1814  cout << EoSctrl.NRsp.conflag << " ";
1815  cout << bulkId << endl;
1816  }
1817  }
1818 
1819  if (nuMax < 1 && EoSctrl.NRsp.conflag && isfinite(eX)) {
1820  // accept this result
1821  } else {
1822  if (!isfinite(eX) || 1 - nuMax < 1E-3) {
1823  NP = 1;
1824  x[0] = zi;
1825  nu[0] = 1;
1826  CalAjBj(Aj[0], Bj[0], x[0]);
1827  SolEoS(Zj[0], Aj[0], Bj[0]);
1828 
1829  EoSctrl.SSMsta.conflag = OCP_FALSE;
1830  EoSctrl.NRsta.conflag = OCP_FALSE;
1831  return OCP_FALSE;
1832  }
1833  }
1834  }
1835  return OCP_TRUE;
1836 }
1837 
1838 void MixtureComp::SplitSSM(const OCP_BOOL& flag)
1839 {
1840  if (NP == 2) {
1841  SplitSSM2(flag);
1842  } else {
1843  SplitSSM3(flag);
1844  }
1845 }
1846 
1847 void MixtureComp::SplitSSM2(const OCP_BOOL& flag)
1848 {
1849  // NP = 2 in this case
1850  // Ks is very IMPORTANT!
1851  // flag = OCP_TRUE : Restart SSM
1852  // flag = OCP_FALSE : New SSM
1853  EoSctrl.SSMsp.conflag = OCP_TRUE;
1854  OCP_DBL Se = 1;
1855  OCP_DBL Stol = EoSctrl.SSMsp.tol2;
1856  USI maxIt = EoSctrl.SSMsp.maxIt;
1857 
1858  if (!flag) {
1859  if (lNP == 2) {
1860  Ks[NP - 2] = lKs;
1861  } else {
1862  if (Yt < 1.1 || OCP_TRUE) {
1863  Ks[NP - 2] = Kw[0];
1864  } else {
1865  for (USI i = 0; i < NC; i++) {
1866  // Ks[NP - 2][i] = phi[testPId][i] / phiSta[i];
1867  Ks[NP - 2][i] = Y[i] / x[testPId][i];
1868  }
1869  }
1870  }
1871  } else {
1872  Stol *= 1E-1;
1873  maxIt *= 2;
1874  }
1875 
1876  if (Yt < 1.1) {
1877  Stol *= 1E-1;
1878  maxIt *= 2;
1879  }
1880 
1881  USI iter = 0;
1882  while (Se > Stol) {
1883 
1884  RachfordRice2();
1885  UpdateXRR();
1886  CalFugPhiAll();
1887  Se = 0;
1888  for (USI i = 0; i < NC; i++) {
1889  Se += pow(fug[1][i] / fug[0][i] - 1, 2);
1890  Ks[0][i] = phi[1][i] / phi[0][i];
1891  }
1892 
1893  iter++;
1894  if (iter > maxIt) {
1895  // OCP_WARNING("SSM not converged in Phase Spliting!");
1896  EoSctrl.SSMsp.conflag = OCP_FALSE;
1897  break;
1898  }
1899  }
1900 
1901  EoSctrl.SSMsp.realTol = sqrt(Se);
1902  EoSctrl.SSMsp.curIt += iter;
1903 }
1904 
1905 void MixtureComp::SplitSSM3(const OCP_BOOL& flag) {}
1906 
1907 void MixtureComp::RachfordRice2()
1908 {
1909  const vector<OCP_DBL>& Ktmp = Ks[0];
1910  OCP_DBL Kmin = Ktmp[0];
1911  OCP_DBL Kmax = Ktmp[0];
1912 
1913  for (USI i = 1; i < NC; i++) {
1914  if (Ktmp[i] < Kmin) Kmin = Ktmp[i];
1915  if (Ktmp[i] > Kmax) Kmax = Ktmp[i];
1916  }
1917 
1918  const OCP_DBL numin = 1 / (1 - Kmax);
1919  const OCP_DBL numax = 1 / (1 - Kmin);
1920 
1921  nu[0] = 0.5 * (numin + numax);
1922 
1923  // Solve RR with NR
1924  OCP_DBL tmp, rj, J, dnuj, tmpnu;
1925 
1926  USI iter = 0;
1927  const OCP_DBL tol = EoSctrl.RR.tol;
1928  const OCP_DBL maxIt = EoSctrl.RR.maxIt;
1929  while (OCP_TRUE) {
1930 
1931  rj = 0;
1932  J = 0;
1933  for (USI i = 0; i < NC; i++) {
1934  tmp = 1 + nu[0] * (Ktmp[i] - 1);
1935  rj += zi[i] * (Ktmp[i] - 1) / tmp;
1936  J -= zi[i] * (Ktmp[i] - 1) * (Ktmp[i] - 1) / (tmp * tmp);
1937  }
1938 
1939  if (fabs(rj) < tol || iter > maxIt) break;
1940 
1941  dnuj = -rj / J;
1942  tmpnu = nu[0] + dnuj;
1943  if (tmpnu < numax && tmpnu > numin) {
1944  nu[0] = tmpnu;
1945  } else {
1946  if (dnuj > 0) {
1947  nu[0] = (nu[0] + numax) / 2;
1948  } else {
1949  nu[0] = (nu[0] + numin) / 2;
1950  }
1951  }
1952  iter++;
1953  }
1954 
1955  EoSctrl.RR.curIt += iter;
1956  nu[1] = 1 - nu[0];
1957 
1958  // cout << scientific << setprecision(6) << nu[0] << " " << nu[1] << endl;
1959 }
1960 
1961 void MixtureComp::RachfordRice2P()
1962 {
1963  // modified RachfordRice equations
1964  // less iterations but more divergence --- unstable!
1965 
1966  const vector<OCP_DBL>& Ktmp = Ks[0];
1967  OCP_DBL Kmin = Ktmp[0];
1968  OCP_DBL Kmax = Ktmp[0];
1969 
1970  for (USI i = 1; i < NC; i++) {
1971  if (Ktmp[i] < Kmin) Kmin = Ktmp[i];
1972  if (Ktmp[i] > Kmax) Kmax = Ktmp[i];
1973  }
1974 
1975  const OCP_DBL numin = 1 / (1 - Kmax);
1976  const OCP_DBL numax = 1 / (1 - Kmin);
1977 
1978  nu[0] = 0.5 * (numin + numax);
1979 
1980  // Solve RR with NR
1981  OCP_DBL tmp, rj, J, dnuj, tmpnu;
1982  OCP_DBL f, df;
1983 
1984  USI iter = 0;
1985  const OCP_DBL tol = EoSctrl.RR.tol;
1986  const OCP_DBL maxIt = EoSctrl.RR.maxIt;
1987  while (OCP_TRUE) {
1988 
1989  rj = 0;
1990  J = 0;
1991  for (USI i = 0; i < NC; i++) {
1992  tmp = 1 + nu[0] * (Ktmp[i] - 1);
1993  rj += zi[i] * (Ktmp[i] - 1) / tmp;
1994  J -= zi[i] * (Ktmp[i] - 1) * (Ktmp[i] - 1) / (tmp * tmp);
1995  }
1996  f = (nu[0] - numin) * (numax - nu[0]);
1997  df = -2 * nu[0] + (numax + numin);
1998  J *= f;
1999  J += df * rj;
2000  rj *= f;
2001 
2002  if (fabs(rj) < tol || iter > maxIt) break;
2003 
2004  dnuj = -rj / J;
2005  tmpnu = nu[0] + dnuj;
2006  if (tmpnu < numax && tmpnu > numin) {
2007  nu[0] = tmpnu;
2008  } else {
2009  if (dnuj > 0) {
2010  nu[0] = (nu[0] + numax) / 2;
2011  } else {
2012  nu[0] = (nu[0] + numin) / 2;
2013  }
2014  }
2015  iter++;
2016  }
2017 
2018  EoSctrl.RR.curIt += iter;
2019 
2020  cout << iter << " " << scientific << setprecision(3) << fabs(rj) << " "
2021  << nu[0] << " " << numin << " " << numax << endl;
2022 
2023  nu[1] = 1 - nu[0];
2024 }
2025 
2026 void MixtureComp::RachfordRice3()
2027 {
2028 }
2029 
2030 void MixtureComp::UpdateXRR()
2031 {
2032  OCP_DBL tmp = 0;
2033  for (USI i = 0; i < NC; i++) {
2034  tmp = 1;
2035  for (USI j = 0; j < NP - 1; j++) {
2036  tmp += nu[j] * (Ks[j][i] - 1);
2037  }
2038  x[NP - 1][i] = zi[i] / tmp;
2039  for (USI j = 0; j < NP - 1; j++) {
2040  x[j][i] = Ks[j][i] * x[NP - 1][i];
2041  }
2042  }
2043 }
2044 
2045 void MixtureComp::SplitBFGS()
2046 {
2047  // Use BFGS to calculate phase splitting
2048  // JmatSP will store the BFGS mat
2049  // resSP will store the resSP - lresSP if necessary
2050  // n will store the n - ln if necessary
2051 
2052  // get initial value, n and ln, resSP and lresSP, H0
2053 
2054  // begin BFGS
2055 }
2056 
2057 void MixtureComp::SplitNR()
2058 {
2059  EoSctrl.NRsp.conflag = OCP_FALSE;
2060  // for (USI j = 0; j < NP; j++) {
2061  // nu[j] = fabs(nu[j]);
2062  // }
2063 
2064  USI len = NC * (NP - 1);
2065  x2n();
2066  CalResSP();
2067  OCP_DBL eNR0;
2068  OCP_DBL eNR = Dnorm2(len, &resSP[0]);
2069  OCP_DBL NRtol = EoSctrl.NRsp.tol;
2070  OCP_DBL alpha;
2071 
2072  OCP_DBL en;
2073  USI iter = 0;
2074  eNR0 = eNR;
2075  while (eNR > NRtol) {
2076 
2077  // eNR0 = eNR;
2078  ln = n;
2079  CalFugNAll();
2080  AssembleJmatSP();
2081 
2082  // LUSolve(1, len, &JmatSP[0], &resSP[0], &pivot[0]);
2083  // PrintDX(NC, &resSP[0]);
2084 
2085  int info = SYSSolve(1, &uplo, len, &JmatSP[0], &resSP[0], &pivot[0],
2086  &JmatWork[0], lJmatWork);
2087  if (info > 0 && false) {
2088  for (USI i = 0; i < len; i++) {
2089  for (USI j = 0; j < len; j++) {
2090  cout << scientific << setprecision(9) << JmatSP[i * len + j]
2091  << " ";
2092  }
2093  cout << ";\n";
2094  }
2095  }
2096  // PrintDX(NC, &resSP[0]);
2097 
2098  alpha = CalStepNRsp();
2099 
2100  n[NP - 1] = zi;
2101  for (USI j = 0; j < NP - 1; j++) {
2102  Daxpy(NC, alpha, &resSP[j * NC], &n[j][0]);
2103  Daxpy(NC, -1, &n[j][0], &n[NP - 1][0]);
2104 
2105  // nu[j] = Dnorm1(NC, &n[j][0]);
2106  nu[j] = 0;
2107  for (USI i = 0; i < NC; i++) {
2108  nu[j] += n[j][i];
2109  }
2110 
2111  for (USI i = 0; i < NC; i++) {
2112  // x[j][i] = n[j][i] / nu[j];
2113  x[j][i] = fabs(n[j][i] / nu[j]);
2114  }
2115  }
2116 
2117  // for (USI i = 0; i < NC; i++) {
2118  // n[NP - 1][i] = fabs(n[NP - 1][i]);
2119  // }
2120  // nu[NP - 1] = Dnorm1(NC, &n[NP - 1][0]);
2121  // for (USI i = 0; i < NC; i++) {
2122  // x[NP - 1][i] = n[NP - 1][i] / nu[NP - 1];
2123  // }
2124  nu[NP - 1] = 0;
2125  for (USI i = 0; i < NC; i++) {
2126  nu[NP - 1] += n[NP - 1][i];
2127  }
2128  for (USI i = 0; i < NC; i++) {
2129  x[NP - 1][i] = fabs(n[NP - 1][i] / nu[NP - 1]);
2130  }
2131 
2132  CalFugPhiAll();
2133  CalResSP();
2134  eNR = Dnorm2(len, &resSP[0]);
2135  iter++;
2136  if (eNR > eNR0 || iter > EoSctrl.NRsp.maxIt) {
2137  break;
2138  }
2139 
2140  // Maybe it should be execute before "eNR > eNR0 || iter > EoSctrl.NRsp.maxIt"
2141  en = 0;
2142  for (USI j = 0; j < NP; j++) {
2143  Daxpy(NC, -1, &n[j][0], &ln[j][0]);
2144  en += Dnorm2(NC, &ln[j][0]);
2145  }
2146  if (en / (NP * NC) < 1E-8) {
2147  EoSctrl.NRsp.conflag = OCP_TRUE;
2148  break;
2149  }
2150  }
2151  EoSctrl.NRsp.realTol = eNR;
2152  if (eNR < NRtol) EoSctrl.NRsp.conflag = OCP_TRUE;
2153  EoSctrl.NRsp.curIt += iter;
2154 
2155  // cout << iter << " " << scientific << setprecision(3) << eNR << endl;
2156 }
2157 
2158 void MixtureComp::CalResSP()
2159 {
2160  // So it equals -res
2161  for (USI j = 0; j < NP - 1; j++) {
2162  for (USI i = 0; i < NC; i++) {
2163  resSP[j * NC + i] = log(fug[NP - 1][i] / fug[j][i]);
2164  }
2165  }
2166 }
2167 
2168 void MixtureComp::CalFugNAll(const OCP_BOOL& Znflag)
2169 {
2170  OCP_DBL C, E, G;
2171  OCP_DBL Cnk, Dnk, Enk, Gnk;
2172  OCP_DBL tmp, aik;
2173 
2174  for (USI j = 0; j < NP; j++) {
2175  // j th phase
2176  vector<OCP_DBL>& fugn = fugN[j];
2177  const OCP_DBL& aj = Aj[j];
2178  const OCP_DBL& bj = Bj[j];
2179  const OCP_DBL& zj = Zj[j];
2180  const vector<OCP_DBL>& xj = x[j];
2181  vector<OCP_DBL>& Znj = Zn[j];
2182 
2183  for (USI i = 0; i < NC; i++) {
2184  tmp = 0;
2185  for (USI m = 0; m < NC; m++) {
2186  tmp += (1 - BIC[i * NC + m]) * sqrt(Ai[i] * Ai[m]) * xj[m];
2187  }
2188  An[i] = 2 / nu[j] * (tmp - aj);
2189  Bn[i] = 1 / nu[j] * (Bi[i] - bj);
2190  if (Znflag) {
2191  Znj[i] =
2192  ((bj - zj) * An[i] +
2193  ((aj + delta1 * delta2 * (3 * bj * bj + 2 * bj)) +
2194  ((delta1 + delta2) * (2 * bj + 1) - 2 * delta1 * delta2 * bj) *
2195  zj -
2196  (delta1 + delta2 - 1) * zj * zj) *
2197  Bn[i]) /
2198  (3 * zj * zj + 2 * ((delta1 + delta2 - 1) * bj - 1) * zj +
2199  (aj + delta1 * delta2 * bj * bj -
2200  (delta1 + delta2) * bj * (bj + 1)));
2201  }
2202  }
2203 
2204  G = (zj + delta1 * bj) / (zj + delta2 * bj);
2205 
2206  for (USI i = 0; i < NC; i++) {
2207  // i th fugacity
2208  C = xj[i] * P / (zj - bj);
2209  // D = Bi[i] / bj * (zj - 1);
2210  tmp = 0;
2211  for (USI k = 0; k < NC; k++) {
2212  tmp += (1 - BIC[i * NC + k]) * sqrt(Ai[i] * Ai[k]) * xj[k];
2213  }
2214  E = -aj / ((delta1 - delta2) * bj) * (2 * tmp / aj - Bi[i] / bj);
2215 
2216  for (USI k = 0; k <= i; k++) {
2217  // k th components
2218 
2219  aik = (1 - BIC[i * NC + k]) * sqrt(Ai[i] * Ai[k]);
2220 
2221  Cnk = P / (zj - bj) / (zj - bj) *
2222  ((zj - bj) / nu[j] * (delta(i, k) - xj[i]) -
2223  xj[i] * (Znj[k] - Bn[k]));
2224  Dnk = Bi[i] / bj * (Znj[k] - (Bi[k] - bj) * (zj - 1) / (nu[j] * bj));
2225  Gnk = (delta1 - delta2) / ((zj + delta2 * bj) * (zj + delta2 * bj)) *
2226  (Bn[k] * zj - Znj[k] * bj);
2227  /*Enk = 1 / ((delta1 - delta2) * bj * bj) * (An[k] * bj - Bn[k] * aj) *
2228  (Bi[i] / bj - 2 * tmp / aj)
2229  + aj / ((delta1 - delta2) * bj) * (-Bi[i] / (bj * bj) * Bn[k] - 2 /
2230  (aj * aj) * (aj * (aik - tmp) / nu[j] - An[k] * tmp));*/
2231  Enk = -1 / (delta1 - delta2) / (bj * bj) *
2232  (2 * (bj * aik / nu[j] + Bn[k] * (Bi[i] * aj / bj - tmp)) -
2233  An[k] * Bi[i] - aj * Bi[i] / nu[j]);
2234  Enk -= E / nu[j];
2235  fugn[i * NC + k] = 1 / C * Cnk + Dnk + Enk * log(G) + E / G * Gnk;
2236  fugn[k * NC + i] = fugn[i * NC + k];
2237  /*cout << "fnn[" << j << "][" << i << "][" << k << "] = " << fugn[i * NC
2238  + k]; cout << endl;*/
2239  }
2240  }
2241  }
2242  // PrintFugN();
2243 #ifdef OCP_NANCHECK
2244  for (USI j = 0; j < NP; j++) {
2245  if (!CheckNan(fugN[j].size(), &fugN[j][0])) {
2246  OCP_ABORT("INF or NAN in fugN !");
2247  }
2248  }
2249 #endif // NANCHECK
2250 }
2251 
2252 void MixtureComp::PrintFugN()
2253 {
2254  for (USI j = 0; j < NP; j++) {
2255  for (USI i = 0; i < NC; i++) {
2256  for (USI k = 0; k < NC; k++) {
2257  cout << fugN[j][i * NC + k] << " ";
2258  }
2259  cout << endl;
2260  }
2261  cout << "---------------" << endl;
2262  }
2263 }
2264 
2265 void MixtureComp::AssembleJmatSP()
2266 {
2267  // Dim: (NP-1)*NC
2268  // Attention that fugNj is sysmetric
2269  fill(JmatSP.begin(), JmatSP.end(), 0);
2270 
2271  OCP_DBL* Jtmp = &JmatSP[0];
2272  const OCP_DBL* fugNp;
2273  const OCP_DBL* fugNj;
2274 
2275  for (USI j = 0; j < NP - 1; j++) {
2276  fugNp = &fugN[NP - 1][0];
2277  fugNj = &fugN[j][0];
2278 
2279  for (USI i = 0; i < NC; i++) {
2280  // ith components
2281  for (USI k = 0; k < NC; k++) {
2282  // kth fugacity
2283  Jtmp[k] = fugNj[k] + fugNp[k];
2284  }
2285  Jtmp += NC * (NP - 1);
2286  fugNp += NC;
2287  fugNj += NC;
2288  }
2289  Jtmp += NC;
2290  }
2291 
2292  // cout << endl << "Jmat" << endl;
2293  // for (USI i = 0; i < NC; i++) {
2294  // for (USI j = 0; j < NC; j++) {
2295  // cout << scientific << setprecision(6) << JmatSP[i * NC + j] << " ";
2296  // }
2297  // cout << endl;
2298  // }
2299 }
2300 
2301 OCP_DBL MixtureComp::CalStepNRsp()
2302 {
2303  OCP_DBL alpha = 1;
2304  OCP_DBL tmp;
2305 
2306  for (USI j = 0; j < NP - 1; j++) {
2307 
2308  const OCP_DBL* nj = &n[j][0];
2309  const OCP_DBL* r = &resSP[j * NC];
2310 
2311  for (USI i = 0; i < NC; i++) {
2312  tmp = nj[i] + alpha * r[i];
2313  if (tmp < 0) {
2314  alpha *= 0.9 * fabs(nj[i] / r[i]);
2315  }
2316  }
2317  }
2318  return alpha;
2319 }
2320 
2321 void MixtureComp::AllocateOthers()
2322 {
2323  sqrtMWi.resize(NC);
2324  for (USI i = 0; i < NC; i++) sqrtMWi[i] = sqrt(MWC[i]);
2325  muC.resize(NPmax);
2326  muAux.resize(NPmax);
2327  for (USI i = 0; i < NPmax; i++) {
2328  muAux[i].resize(5);
2329  }
2330  muAux1I.resize(NC);
2331  for (USI i = 0; i < NC; i++) {
2332  muAux1I[i] = 5.4402 * pow(Tc[i], 1.0 / 6) / pow(Pc[i], 2.0 / 3);
2333  }
2334  fugP.resize(NPmax);
2335  for (USI j = 0; j < NPmax; j++) {
2336  fugP[j].resize(NC);
2337  }
2338  Zp.resize(NPmax);
2339  JmatDer.resize(NPmax * NPmax * (NC + 1) * (NC + 1));
2340  JmatTmp = JmatDer;
2341  rhsDer.resize(NPmax * (NC + 1) * (NC + 1));
2342  xixC.resize(NPmax * NC);
2343  xiPC.resize(NPmax);
2344  xiNC.resize(NPmax * NC);
2345 
2346  // new
2347  JmatDer.resize((numPhase + NPmax * NC) * (numPhase + NPmax * NC));
2348  JmatTmp = JmatDer;
2349  rhsDer.resize((numPhase + NPmax * NC) * (numCom + 1 + 1));
2350 }
2351 
2352 void MixtureComp::IdentifyPhase()
2353 {
2354  phaseExist[0] = OCP_FALSE;
2355  phaseExist[1] = OCP_FALSE;
2356  if (NP == 1) {
2357  // Critical Temperature Method
2358  OCP_DBL A = 0;
2359  OCP_DBL B = 0;
2360  for (USI i = 0; i < NC; i++) {
2361  A += x[0][i] * Vc[i] * Tc[i];
2362  B += x[0][i] * Vc[i];
2363  }
2364  OCP_DBL Tc = A / B;
2365  if (T > Tc) {
2366  phaseLabel[0] = GAS;
2367  phaseExist[1] = OCP_TRUE;
2368  } else {
2369  phaseLabel[0] = OIL;
2370  phaseExist[0] = OCP_TRUE;
2371  }
2372  } else {
2373  // Compare MW
2374  if (MW[0] > MW[1]) {
2375  phaseLabel[0] = OIL;
2376  phaseLabel[1] = GAS;
2377  } else {
2378  phaseLabel[0] = GAS;
2379  phaseLabel[1] = OIL;
2380  }
2381  phaseExist[0] = OCP_TRUE;
2382  phaseExist[1] = OCP_TRUE;
2383  }
2384 }
2385 
2386 void MixtureComp::CopyPhase()
2387 {
2388  // copy vj, x, mu, xi, rho
2389  for (USI j = 0; j < NP; j++) {
2390  const USI j1 = phaseLabel[j];
2391  vj[j1] = vC[j];
2392  Dcopy(NC, &xij[j1 * numCom], &x[j][0]);
2393  mu[j1] = muC[j];
2394  xi[j1] = xiC[j];
2395  rho[j1] = rhoC[j];
2396  nj[j1] = nu[j];
2397  }
2398 }
2399 
2400 void MixtureComp::CalViscosity() { CalViscoLBC(); }
2401 
2402 void MixtureComp::CalViscoLBC()
2403 {
2404  OCP_DBL tmp;
2405  OCP_DBL Tri;
2406  OCP_DBL xijT;
2407  OCP_DBL xijP;
2408  OCP_DBL xijV;
2409 
2410  for (USI j = 0; j < NP; j++) {
2411  const vector<OCP_DBL>& xj = x[j];
2412  vector<OCP_DBL>& muA = muAux[j];
2413  fill(muA.begin(), muA.end(), 0.0);
2414  xijT = 0;
2415  xijP = 0;
2416  xijV = 0;
2417 
2418  for (USI i = 0; i < NC; i++) {
2419  tmp = 5.4402 * pow(Tc[i], 1.0 / 6) / sqrt(MW[j]) / pow(Pc[i], 2.0 / 3);
2420  // tmp = muAux1I[i] / sqrt(MW[j]);
2421  Tri = T / Tc[i];
2422  if (Tri <= 1.5) {
2423  tmp = 34 * 1E-5 * pow(Tri, 0.94) / tmp;
2424  } else {
2425  tmp = 17.78 * 1E-5 * pow((4.58 * Tri - 1.67), 0.625) / tmp;
2426  }
2427  muA[0] += xj[i] * sqrt(MWC[i]) * tmp;
2428  muA[1] += xj[i] * sqrt(MWC[i]);
2429  xijT += xj[i] * Tc[i];
2430  xijP += xj[i] * Pc[i];
2431  xijV += xj[i] * Vcvis[i];
2432  }
2433  muA[2] = 5.4402 * pow(xijT, 1.0 / 6) / sqrt(MW[j]) / pow(xijP, 2.0 / 3);
2434  muA[3] = xiC[j] * xijV;
2435 
2436  if (muA[3] <= 0.18 && OCP_FALSE) {
2437  muC[j] = muA[0] / muA[1] + 2.05 * 1E-4 * muA[3] / muA[2];
2438  } else {
2439  muA[4] = muA[3] * (muA[3] * (muA[3] * (LBCcoef[4] * muA[3] + LBCcoef[3]) +
2440  LBCcoef[2]) +
2441  LBCcoef[1]) +
2442  LBCcoef[0];
2443  muC[j] = muA[0] / muA[1] + 1E-4 * (pow(muA[4], 4) - 1) / muA[2];
2444  }
2445  }
2446 }
2447 
2448 void MixtureComp::CalViscoHZYT() {}
2449 
2450 void MixtureComp::CalFugXAll()
2451 {
2452  for (USI j = 0; j < NP; j++) {
2453 
2454  vector<OCP_DBL>& fugx = fugX[j];
2455  vector<OCP_DBL>& xj = x[j];
2456  OCP_DBL aj = Aj[j];
2457  OCP_DBL bj = Bj[j];
2458  OCP_DBL zj = Zj[j];
2459  OCP_DBL tmp = 0;
2460 
2461  Bx = Bi;
2462  for (USI i = 0; i < NC; i++) {
2463  tmp = 0;
2464  for (USI k = 0; k < NC; k++) {
2465  tmp += xj[k] * (1 - BIC[i * NC + k]) * sqrt(Ai[i] * Ai[k]);
2466  }
2467  Ax[i] = 2 * tmp;
2468  Zx[i] =
2469  ((bj - zj) * Ax[i] +
2470  ((aj + delta1 * delta2 * (3 * bj * bj + 2 * bj)) +
2471  ((delta1 + delta2) * (2 * bj + 1) - 2 * delta1 * delta2 * bj) * zj -
2472  (delta1 + delta2 - 1) * zj * zj) *
2473  Bx[i]) /
2474  (3 * zj * zj + 2 * ((delta1 + delta2 - 1) * bj - 1) * zj +
2475  (aj + delta1 * delta2 * bj * bj - (delta1 + delta2) * bj * (bj + 1)));
2476  }
2477 
2478  OCP_DBL C, E, G;
2479  OCP_DBL Cxk, Dxk, Exk, Gxk;
2480  OCP_DBL aik;
2481  G = (zj + delta1 * bj) / (zj + delta2 * bj);
2482 
2483  for (USI i = 0; i < NC; i++) {
2484  // ith fugacity
2485  C = xj[i] * P / (zj - bj);
2486  // C = 1 / (zj - bj);
2487  // D = Bx[i] * (zj - 1) / bj;
2488  E = -aj / ((delta1 - delta2) * bj) * (Ax[i] / aj - Bx[i] / bj);
2489 
2490  for (USI k = 0; k < NC; k++) {
2491  aik = (1 - BIC[i * NC + k]) * sqrt(Ai[i] * Ai[k]);
2492 
2493  // kth components
2494  Cxk = ((zj - bj) * delta(i, k) - xj[i] * (Zx[k] - Bx[k])) * P /
2495  ((zj - bj) * (zj - bj));
2496  Dxk = Bx[i] / bj * (Zx[k] - Bx[k] * (zj - 1) / bj);
2497  /*Exk = (Ax[k] * bj - aj * Bx[k]) / (bj * bj) * (Ax[i] / aj - Bx[i] /
2498  bj) + aj / bj * (2 * (1 - BIC[i * NC + k]) * sqrt(Ai[i] * Ai[k]) / aj
2499  - Ax[k] * Ax[i] / (aj * aj) + Bx[i] * Bx[k] / (bj * bj));*/
2500  Exk = (2 * (aj / bj * Bx[k] * Bx[i] + bj * aik) - Ax[i] * Bx[k] -
2501  Ax[k] * Bi[i]) /
2502  (bj * bj);
2503  Exk /= -(delta1 - delta2);
2504  Gxk = (delta1 - delta2) / (zj + delta2 * bj) / (zj + delta2 * bj) *
2505  (zj * Bx[k] - Zx[k] * bj);
2506  fugx[i * NC + k] = 1 / C * Cxk + Dxk + Exk * log(G) + E / G * Gxk;
2507  }
2508  }
2509  }
2510 #ifdef OCP_NANCHECK
2511  for (USI j = 0; j < NP; j++) {
2512  if (!CheckNan(fugX[j].size(), &fugX[j][0])) {
2513  OCP_ABORT("INF or NAN in fugX !");
2514  }
2515  }
2516 #endif // NANCHECK
2517 }
2518 
2519 void MixtureComp::CalFugPAll(const OCP_BOOL& Zpflag)
2520 {
2521 
2522  OCP_DBL C, E, G;
2523  OCP_DBL Cp, Dp, Gp;
2524  OCP_DBL tmp;
2525 
2526  for (USI j = 0; j < NP; j++) {
2527 
2528  vector<OCP_DBL>& fugp = fugP[j];
2529  vector<OCP_DBL>& xj = x[j];
2530  OCP_DBL& aj = Aj[j];
2531  OCP_DBL& bj = Bj[j];
2532  OCP_DBL& zj = Zj[j];
2533 
2534  OCP_DBL Ap = aj / P;
2535  OCP_DBL Bp = bj / P;
2536  if (Zpflag) {
2537  Zp[j] =
2538  ((bj - zj) * Ap +
2539  ((aj + delta1 * delta2 * (3 * bj * bj + 2 * bj)) +
2540  ((delta1 + delta2) * (2 * bj + 1) - 2 * delta1 * delta2 * bj) * zj -
2541  (delta1 + delta2 - 1) * zj * zj) *
2542  Bp) /
2543  (3 * zj * zj + 2 * ((delta1 + delta2 - 1) * bj - 1) * zj +
2544  (aj + delta1 * delta2 * bj * bj - (delta1 + delta2) * bj * (bj + 1)));
2545  }
2546 
2547  G = (zj + delta1 * bj) / (zj + delta2 * bj);
2548  Gp = (delta1 - delta2) / ((zj + delta2 * bj) * (zj + delta2 * bj)) *
2549  (Bp * zj - Zp[j] * bj);
2550  for (USI i = 0; i < NC; i++) {
2551 
2552  C = P / (zj - bj);
2553  // D = Bi[i] / bj * (zj - 1);
2554 
2555  tmp = 0;
2556  for (USI m = 0; m < NC; m++) {
2557  tmp += (1 - BIC[i * NC + m]) * sqrt(Ai[i] * Ai[m]) * xj[m];
2558  }
2559 
2560  E = -aj / ((delta1 - delta2) * bj) * (2 * tmp / aj - Bi[i] / bj);
2561 
2562  Cp = ((zj - bj) - P * (Zp[j] - Bp)) / ((zj - bj) * (zj - bj));
2563  Dp = Bi[i] / bj * Zp[j];
2564  // Ep = 0;
2565 
2566  // Attention that if xj[i] = 0, then fugp[i] = nan
2567  // but Cp also contains xj[i], so eliminate it first
2568  fugp[i] = Cp / C + Dp + E / G * Gp;
2569  }
2570  }
2571 
2572 #ifdef OCP_NANCHECK
2573  for (USI j = 0; j < NP; j++) {
2574  if (!CheckNan(fugP[j].size(), &fugP[j][0])) {
2575  OCP_ABORT("INF or NAN in fugP !");
2576  }
2577  }
2578 #endif // NANCHECK
2579 }
2580 
2581 void MixtureComp::CalVjpVfpVfn_partial()
2582 {
2583  // Vfi = 0
2584  // Vjp, Vfp
2585  // Vjn(dVj / dnij) in Vji
2586  OCP_DBL CgTP = GAS_CONSTANT * T / P;
2587  fill(vfi.begin(), vfi.end(), 0.0);
2588  vfP = 0;
2589 
2590  for (USI j = 0; j < NP; j++) {
2591  const OCP_DBL& aj = Aj[j];
2592  const OCP_DBL& bj = Bj[j];
2593  const OCP_DBL& zj = Zj[j];
2594  const vector<OCP_DBL>& xj = x[j];
2595  vector<OCP_DBL>& Znj = Zn[j];
2596  OCP_DBL tmp;
2597  const USI j1 = phaseLabel[j];
2598  for (USI i = 0; i < NC; i++) {
2599  tmp = 0;
2600  for (USI m = 0; m < NC; m++) {
2601  tmp += (1 - BIC[i * NC + m]) * sqrt(Ai[i] * Ai[m]) * xj[m];
2602  }
2603  An[i] = 2 / nu[j] * (tmp - aj);
2604  Bn[i] = 1 / nu[j] * (Bi[i] - bj);
2605  Znj[i] =
2606  ((bj - zj) * An[i] +
2607  ((aj + delta1 * delta2 * (3 * bj * bj + 2 * bj)) +
2608  ((delta1 + delta2) * (2 * bj + 1) - 2 * delta1 * delta2 * bj) * zj -
2609  (delta1 + delta2 - 1) * zj * zj) *
2610  Bn[i]) /
2611  (3 * zj * zj + 2 * ((delta1 + delta2 - 1) * bj - 1) * zj +
2612  (aj + delta1 * delta2 * bj * bj - (delta1 + delta2) * bj * (bj + 1)));
2613  vji[j1][i] = CgTP * (zj + nu[j] * Znj[i]) - Vshift[i];
2614  }
2615  Zp[j] = ((bj - zj) * aj +
2616  ((aj + delta1 * delta2 * (3 * bj * bj + 2 * bj)) +
2617  ((delta1 + delta2) * (2 * bj + 1) - 2 * delta1 * delta2 * bj) * zj -
2618  (delta1 + delta2 - 1) * zj * zj) *
2619  bj) /
2620  P /
2621  (3 * zj * zj + 2 * ((delta1 + delta2 - 1) * bj - 1) * zj +
2622  (aj + delta1 * delta2 * bj * bj - (delta1 + delta2) * bj * (bj + 1)));
2623  vjp[j1] = CgTP * nu[j] * (Zp[j] - zj / P);
2624  vfP += vjp[j1];
2625  }
2626 }
2627 
2628 void MixtureComp::CalXiPn_partial()
2629 {
2630  // xiN = 0
2631  // xiP
2632  // xin(dxi_j / d nij) in xix
2633  OCP_DBL CgTP = GAS_CONSTANT * T / P;
2634  OCP_DBL tmp;
2635  fill(xiN.begin(), xiN.end() - numCom, 0.0);
2636 
2637  for (USI j = 0; j < NP; j++) {
2638  const USI j1 = phaseLabel[j];
2639 
2640  tmp = -xi[j1] * xi[j1] * CgTP;
2641  for (USI i = 0; i < NC; i++) {
2642  xix[j1 * numCom + i] = tmp * Zn[j][i];
2643  }
2644  xiP[j1] = tmp * (Zp[j] - Zj[j] / P);
2645  }
2646 }
2647 
2648 void MixtureComp::CalRhoPn_partial()
2649 {
2650  // rhoN = 0
2651  // rhoP
2652  // rhon(drhoj / d nij) in rhox
2653  fill(rhoN.begin(), rhoN.end() - numCom, 0.0);
2654  USI j1;
2655  for (USI j = 0; j < NP; j++) {
2656  j1 = phaseLabel[j];
2657  rhoP[j1] = xiP[j1] * MW[j];
2658  for (USI m = 0; m < NC; m++) {
2659  rhox[j1 * numCom + m] =
2660  xix[j1 * numCom + m] * MW[j] + xi[j1] / nu[j] * (MWC[m] - MW[j]);
2661  }
2662  }
2663 }
2664 
2665 void MixtureComp::CalMuPn_partial() { CalMuPnLBC_partial(); }
2666 
2667 void MixtureComp::CalMuPnLBC_partial()
2668 {
2669  // muN = 0
2670  // muP
2671  // mun in mux
2672 
2673  fill(muN.begin(), muN.end() - numCom, 0.0);
2674 
2675  OCP_DBL val1IJ, val2IJ;
2676  OCP_DBL der1IJ, der2IJ, der3J, der4J, der6J, der7J, der8J;
2677  OCP_DBL Tri, tmp;
2678  OCP_DBL xTj, xPj, xVj;
2679  OCP_DBL derxTj, derxPj, derMWj;
2680 
2681  // dxij / dP = 0, d MJ / dP = 0
2682  // Calculate dmuj / dP
2683  // der2IJ = der3J = der4J = der6J = 0;
2684 
2685  for (USI j = 0; j < NP; j++) {
2686  const USI j1 = phaseLabel[j];
2687  const vector<OCP_DBL>& xj = x[j];
2688  const vector<OCP_DBL>& muAuxj = muAux[j];
2689 
2690  xTj = xPj = xVj = 0;
2691  for (USI i = 0; i < NC; i++) {
2692  xTj += xj[i] * Tc[i];
2693  xPj += xj[i] * Pc[i];
2694  xVj += xj[i] * Vcvis[i];
2695  }
2696  der7J = xVj * xiP[j1];
2697  if (muAuxj[3] <= 0.18 && OCP_FALSE) {
2698  muP[j1] = (2.05 * 1E-4) * der7J / muAuxj[2];
2699  } else {
2700  der8J = der7J * (LBCcoef[1] +
2701  muAuxj[3] * (2 * LBCcoef[2] +
2702  muAuxj[3] * (3 * LBCcoef[3] +
2703  muAuxj[3] * 4 * LBCcoef[4])));
2704  muP[j1] = (4 * 1E-4) * pow(muAuxj[4], 3) * der8J / muAuxj[2];
2705  }
2706 
2707  // Calculate dmuj / nkj
2708  const USI bId = numCom * j1;
2709  for (USI k = 0; k < NC; k++) {
2710  derxTj = (Tc[k] - xTj) / nu[j];
2711  derxPj = (Pc[k] - xPj) / nu[j];
2712  derMWj = (MWC[k] - MW[j]) / nu[j];
2713  der3J = 0;
2714  der4J = sqrtMWi[k];
2715  for (USI i = 0; i < NC; i++) {
2716  val1IJ = muAux1I[i] / sqrt(MW[j]);
2717  der1IJ = -(1 / 2) * muAux1I[i] * pow(MW[j], -1.5) * derMWj;
2718  Tri = T / Tc[i];
2719  if (Tri <= 1.5) {
2720  tmp = 34 * 1E-5 * pow(Tri, 0.94);
2721  } else {
2722  tmp = 17.78 * 1E-5 * pow(4.58 * Tri - 1.67, 0.625);
2723  }
2724  val2IJ = tmp / val1IJ;
2725  der2IJ = -tmp * der1IJ / (val1IJ * val1IJ);
2726  der3J += sqrtMWi[i] *
2727  (xj[i] * der2IJ + (delta(i, k) - xj[i]) * val2IJ / nu[j]);
2728  der4J -= xj[i] * sqrtMWi[i];
2729  }
2730  der4J /= nu[j];
2731  der6J =
2732  5.4402 *
2733  (1.0 / 6 * pow(xTj, -5.0 / 6) * derxTj -
2734  pow(xTj, 1.0 / 6) * (0.5 / MW[j] * derMWj + 2.0 / 3 / xPj * derxPj)) /
2735  (sqrt(MW[j]) * pow(xPj, 2.0 / 3));
2736  der7J = xix[j1 * numCom + k] * xVj + (Vcvis[k] - xVj) * xi[j1] / nu[j];
2737  if (muAuxj[3] <= 0.18 && OCP_FALSE) {
2738  mux[bId + k] =
2739  (der3J * muAuxj[1] - muAuxj[0] * der4J) / (muAuxj[1] * muAuxj[1]) +
2740  2.05 * 1E-4 * (der7J * muAuxj[2] - muAuxj[3] * der6J) /
2741  (muAuxj[2] * muAuxj[2]);
2742  } else {
2743  der8J =
2744  der7J * (LBCcoef[1] +
2745  muAuxj[3] * (2 * LBCcoef[2] +
2746  muAuxj[3] * (3 * LBCcoef[3] +
2747  muAuxj[3] * 4 * LBCcoef[4])));
2748  mux[bId + k] =
2749  (der3J * muAuxj[1] - muAuxj[0] * der4J) / (muAuxj[1] * muAuxj[1]) +
2750  1E-4 *
2751  (4 * pow(muAuxj[4], 3) * der8J * muAuxj[2] -
2752  (pow(muAuxj[4], 4) - 1) * der6J) /
2753  (muAuxj[2] * muAuxj[2]);
2754  }
2755  }
2756  }
2757 }
2758 
2759 void MixtureComp::CalVjpVfpVfx_partial()
2760 {
2761  // Wrong Now !
2762 
2763  // Vfi = 0
2764  // Vjp, Vfp
2765  // Vjn(dVj / dxij) in Vji
2766  OCP_DBL CgTP = GAS_CONSTANT * T / P;
2767  fill(vfi.begin(), vfi.end(), 0.0);
2768  vfP = 0;
2769 
2770  for (USI j = 0; j < NP; j++) {
2771  const USI j1 = phaseLabel[j];
2772  const vector<OCP_DBL>& xj = x[j];
2773  const OCP_DBL aj = Aj[j];
2774  const OCP_DBL bj = Bj[j];
2775  const OCP_DBL zj = Zj[j];
2776  OCP_DBL tmp = 0;
2777 
2778  Bx = Bi;
2779  for (USI i = 0; i < NC; i++) {
2780  tmp = 0;
2781  for (USI k = 0; k < NC; k++) {
2782  tmp += xj[k] * (1 - BIC[i * NC + k]) * sqrt(Ai[i] * Ai[k]);
2783  }
2784  Ax[i] = 2 * tmp;
2785  Zx[i] =
2786  ((bj - zj) * Ax[i] +
2787  ((aj + delta1 * delta2 * (3 * bj * bj + 2 * bj)) +
2788  ((delta1 + delta2) * (2 * bj + 1) - 2 * delta1 * delta2 * bj) * zj -
2789  (delta1 + delta2 - 1) * zj * zj) *
2790  Bx[i]) /
2791  (3 * zj * zj + 2 * ((delta1 + delta2 - 1) * bj - 1) * zj +
2792  (aj + delta1 * delta2 * bj * bj - (delta1 + delta2) * bj * (bj + 1)));
2793 
2794  vji[j1][i] = CgTP * nu[j] * Zx[i];
2795  }
2796  Zp[j] = ((bj - zj) * aj +
2797  ((aj + delta1 * delta2 * (3 * bj * bj + 2 * bj)) +
2798  ((delta1 + delta2) * (2 * bj + 1) - 2 * delta1 * delta2 * bj) * zj -
2799  (delta1 + delta2 - 1) * zj * zj) *
2800  bj) /
2801  P /
2802  (3 * zj * zj + 2 * ((delta1 + delta2 - 1) * bj - 1) * zj +
2803  (aj + delta1 * delta2 * bj * bj - (delta1 + delta2) * bj * (bj + 1)));
2804  vjp[j1] = CgTP * nu[j] * (Zp[j] - zj / P);
2805  vfP += vjp[j1];
2806  }
2807 }
2808 
2809 void MixtureComp::CaldXsdXpAPI04()
2810 {
2811  // Wrong Now !
2812 
2813  // water not in oil or gas; hydroncarbon not in water
2814  // inexist phase will bot be stored
2815  // only hydrocarbon phases, hydrocarbon components, water phase will be stored
2816 
2817  // dS / dP
2818  // S = Sj, xij
2819  // P = P, Ni
2820  // water is included
2821  fill(dXsdXp.begin(), dXsdXp.end(), 0.0);
2822  fill(res.begin(), res.end(), 0.0);
2823  resPc = 0;
2824 
2825  const USI ncol = numCom + 1;
2826  const OCP_DBL vf2 = vf * vf;
2827  const OCP_DBL vwp = vjp[numPhase - 1];
2828  const OCP_DBL vw = vj[numPhase - 1];
2829 
2830  if (NP == 1) {
2831  // when NP = 1, vfi = vhi, vfP = vhp + vwp, h = hydrocarbon
2832  const USI j1 = phaseLabel[0];
2833  OCP_DBL* bId = &dXsdXp[0];
2834  // dS / dP
2835  bId[0] = ((vfP - vwp) * vf - vfP * vj[j1]) / vf2;
2836  bId++;
2837  // dSh / dNm
2838  for (USI m = 0; m < NC; m++) {
2839  bId[m] = vfi[m] * vw / vf2;
2840  }
2841  bId += NC;
2842  // dSh / dNw
2843  bId[0] = -vfi[numCom - 1] * vj[j1] / vf2;
2844  bId++;
2845  // dSw / dP, dNm
2846  for (USI m = 0; m < ncol; m++) {
2847  bId[m] = -dXsdXp[m];
2848  }
2849  bId += ncol + 1;
2850  // dxij / dNm
2851  for (USI i = 0; i < NC; i++) {
2852  for (USI m = 0; m < NC; m++) {
2853  bId[m] = (delta(i, m) * Nh - Ni[i]) / (Nh * Nh);
2854  }
2855  bId += ncol;
2856  }
2857  // res = 0
2858  } else {
2859  CalFugXAll();
2860  CalFugPAll(OCP_FALSE);
2861  CaldXsdXp04();
2862  }
2863 }
2864 
2865 void MixtureComp::CaldXsdXp04()
2866 {
2867  // Wrong Now !
2868 
2869  fill(JmatTmp.begin(), JmatTmp.end(), 0.0);
2870  OCP_DBL* MTmp = &JmatTmp[0];
2871 
2872  // -dF / dXs
2873  // -dFn / dXs
2874  const USI dim = NP + 1 + NP * NC;
2875  for (USI i = 0; i < NC; i++) {
2876  // -dFn / dSh
2877  for (USI m = 0; m < NP; m++) {
2878  MTmp[m] = vf * xi[phaseLabel[m]] * x[m][i];
2879  }
2880  MTmp += NP;
2881  // -dFn / dSw
2882  MTmp++;
2883  // -dFn / dxnm
2884  for (USI m = 0; m < NP; m++) {
2885  const USI m1 = phaseLabel[m];
2886  for (USI n = 0; n < NC; n++) {
2887  MTmp[n] = vf * S[m1] *
2888  (xix[m1 * numCom + n] * x[m][i] + delta(i, n) * xi[m1]);
2889  for (USI j = 0; j < NP; j++) {
2890  MTmp[n] +=
2891  vji[m1][n] * x[m][i] * S[phaseLabel[j]] * xi[phaseLabel[j]];
2892  }
2893  }
2894  MTmp += NC;
2895  }
2896  }
2897  // -dFnw / dFs
2898  // -dFnw / dSw
2899  MTmp[numPhase - 1] = vf * xi[numPhase - 1];
2900  MTmp += (NP + 1);
2901  // -dFnw / dxnm
2902  for (USI m = 0; m < NP; m++) {
2903  for (USI n = 0; n < NC; n++) {
2904  MTmp[n] = vji[phaseLabel[m]][n] * xi[numPhase - 1] * S[numPhase - 1];
2905  }
2906  MTmp += NC;
2907  }
2908 
2909  // -dFx / dXs
2910  for (USI j = 0; j < NP; j++) {
2911  // -dFx / dSm
2912  MTmp += (NP + 1);
2913  // -dFx / dxnm, m = j
2914  MTmp += j * NC;
2915  for (USI n = 0; n < NC; n++) {
2916  MTmp[n] = -1.0;
2917  }
2918  MTmp += (NP - j) * NC;
2919  }
2920 
2921  // -dFf / dXs
2922  for (USI j = 0; j < NP - 1; j++) {
2923  const OCP_DBL* fugXJ = &fugX[j][0];
2924  const OCP_DBL* fugXNP = &fugX[NP - 1][0];
2925  for (USI i = 0; i < NC; i++) {
2926  // -dFf / dSm
2927  MTmp += (NP + 1);
2928  // -dFf / dxnm, m = j or m = np-1
2929  // m = j
2930  MTmp += j * NC;
2931  for (USI n = 0; n < NC; n++) {
2932  MTmp[n] = -fugXJ[n];
2933  }
2934  fugXJ += NC;
2935  // m = np-1
2936  MTmp += (NP - 1 - j) * NC;
2937  for (USI n = 0; n < NC; n++) {
2938  MTmp[n] = fugXNP[n];
2939  }
2940  fugXNP += NC;
2941  MTmp += NC;
2942  }
2943  }
2944 
2945  cout << endl << "dFdS" << endl;
2946  cout << scientific << setprecision(1);
2947  for (USI i = 0; i < dim; i++) {
2948  for (USI j = 0; j < dim; j++) cout << setw(8) << JmatTmp[i * dim + j] << " ";
2949  cout << endl;
2950  }
2951 
2952  // Transpose JmatTmp
2953 
2954  for (USI i = 0; i < dim; i++) {
2955  for (USI j = 0; j < dim; j++) {
2956  JmatDer[i * dim + j] = JmatTmp[j * dim + i];
2957  }
2958  }
2959 
2960  // dF / dXp
2961  // d fij / dP have been calculated in CalVfiVfp_full01() before
2962  fill(JmatTmp.begin(), JmatTmp.end(), 0.0);
2963  MTmp = &JmatTmp[0];
2964  OCP_DBL tmp01;
2965  OCP_DBL tmp02;
2966  USI j1;
2967  const USI nrhs = numCom + 1;
2968  const USI nrow = dim;
2969  // dFn / dXp
2970  for (USI i = 0; i < NC; i++) {
2971  // dFn / dP
2972  tmp01 = 0;
2973  tmp02 = 0;
2974  for (USI j = 0; j < NP; j++) {
2975  j1 = phaseLabel[j];
2976  tmp01 += S[j1] * x[j][i] * xiP[j1];
2977  tmp02 += S[j1] * x[j][i] * xi[j1];
2978  }
2979  MTmp[0] = -(vf * tmp01 + vfP * tmp02);
2980  MTmp += 1;
2981  // dFn / dNk
2982  MTmp[i] = 1;
2983  MTmp += NC;
2984  // dFn / dNw
2985  MTmp[0] = -vfi[numCom - 1] * tmp02;
2986  MTmp++;
2987  }
2988  // dFnw / dXp
2989  // dFnw / dP
2990  MTmp[0] = -S[numPhase - 1] * (vfP * xi[numPhase - 1] + vf * xiP[numPhase - 1]);
2991  // dFnw / dNw
2992  MTmp[numCom] = 1 - vfi[numCom - 1] * S[numPhase - 1] * xi[numPhase - 1];
2993  MTmp += nrhs;
2994 
2995  // dFx / dXp
2996  MTmp += NP * nrhs;
2997 
2998  // dFf / dXp
2999  const vector<OCP_DBL>& fugPNP = fugP[NP - 1];
3000  for (USI j = 0; j < NP - 1; j++) {
3001  const vector<OCP_DBL>& fugPj = fugP[j];
3002  for (USI i = 0; i < NC; i++) {
3003  // dFf / dP
3004  MTmp[0] = fugPj[i] - fugPNP[i];
3005  MTmp += nrhs;
3006  }
3007  }
3008 
3009  cout << endl << "dFdp" << endl;
3010  cout << scientific << setprecision(3);
3011  for (USI i = 0; i < dim; i++) {
3012  for (USI j = 0; j < numCom + 1; j++)
3013  cout << setw(10) << JmatTmp[i * (numCom + 1) + j] << " ";
3014  cout << endl;
3015  }
3016 
3017  // Transpose JmatTmp
3018  for (USI i = 0; i < nrhs; i++) {
3019  for (USI j = 0; j < nrow; j++) {
3020  rhsDer[i * nrow + j] = JmatTmp[j * nrhs + i];
3021  }
3022  }
3023 
3024  LUSolve(nrhs, dim, &JmatDer[0], &rhsDer[0], &pivot[0]);
3025 
3026  cout << setprecision(6);
3027  cout << endl << "dXsdXp" << endl;
3028  for (USI i = 0; i < nrow; i++) {
3029  for (USI j = 0; j < nrhs; j++) {
3030  cout << setw(13) << rhsDer[j * nrow + i] << " ";
3031  }
3032  cout << endl;
3033  }
3034  cout << endl;
3035 }
3036 
3037 void MixtureComp::CalXiPNX_partial()
3038 {
3039  // Here!
3040  // dxi/dP = dxi/dP
3041  // dxi/dNk = dxi/dNk = 0
3042  // dxij / dP = dxij / dNk = 0
3043  // See MixtureComp::CalXiPNX_full01()
3044 
3045  // Calculate xix and xiP
3046  OCP_DBL CgTP = GAS_CONSTANT * T / P;
3047  OCP_DBL tmp;
3048 
3049  for (USI j = 0; j < NP; j++) {
3050  const vector<OCP_DBL>& xj = x[j];
3051  OCP_DBL aj = Aj[j];
3052  OCP_DBL bj = Bj[j];
3053  OCP_DBL zj = Zj[j];
3054 
3055  Bx = Bi;
3056  for (USI i = 0; i < NC; i++) {
3057  tmp = 0;
3058  for (USI k = 0; k < NC; k++) {
3059  tmp += xj[k] * (1 - BIC[i * NC + k]) * sqrt(Ai[i] * Ai[k]);
3060  }
3061  Ax[i] = 2 * tmp;
3062  Zx[i] =
3063  ((bj - zj) * Ax[i] +
3064  ((aj + delta1 * delta2 * (3 * bj * bj + 2 * bj)) +
3065  ((delta1 + delta2) * (2 * bj + 1) - 2 * delta1 * delta2 * bj) * zj -
3066  (delta1 + delta2 - 1) * zj * zj) *
3067  Bx[i]) /
3068  (3 * zj * zj + 2 * ((delta1 + delta2 - 1) * bj - 1) * zj +
3069  (aj + delta1 * delta2 * bj * bj - (delta1 + delta2) * bj * (bj + 1)));
3070  }
3071 
3072  tmp = -xiC[j] * xiC[j];
3073  for (USI i = 0; i < NC; i++) {
3074  xixC[j * NC + i] = tmp * (CgTP * Zx[i] - Vshift[i]);
3075  }
3076  xiPC[j] = tmp * CgTP * (Zp[j] - Zj[j] / P);
3077  }
3078 
3079  // xiNC
3080  fill(xiNC.begin(), xiNC.end(), 0.0);
3081 
3082  // Copoy to xiP, xix
3083  // fill(xiP.begin(), xiP.end(), 0.0);
3084  // fill(xix.begin(), xix.end(), 0.0);
3085  fill(xiN.begin(), xiN.end(), 0.0);
3086  USI j1;
3087  for (USI j = 0; j < NP; j++) {
3088  j1 = phaseLabel[j];
3089  xiP[j1] = xiPC[j];
3090  Dcopy(NC, &xix[j1 * numCom], &xixC[j * NC]);
3091  }
3092 }
3093 
3094 void MixtureComp::CalRhoPX_partial()
3095 {
3096  // Here!
3097  // drho/dP = drho/dP
3098  // drho/dNk = drho/dNk = 0
3099  // dxij / dP = dxij / dNk = 0
3100  // See MixtureComp::CalRhoPNX_full()
3101 
3102  // fill(rhox.begin(), rhox.end() - numCom, 0.0);
3103  fill(rhoN.begin(), rhoN.end() - numCom, 0.0);
3104  // fill(rhoP.begin(), rhoP.end() - 1, 0.0);
3105 
3106  USI j1;
3107  for (USI j = 0; j < NP; j++) {
3108  j1 = phaseLabel[j];
3109  rhoP[j1] = xiP[j1] * MW[j];
3110  for (USI i = 0; i < NC; i++) {
3111  rhox[j1 * numCom + i] = xix[j1 * numCom + i] * MW[j] + xi[j1] * MWC[i];
3112  }
3113  }
3114 }
3115 
3116 void MixtureComp::CalMuPX_partial()
3117 {
3118  // fill(muP.begin(), muP.end() - 1, 0.0);
3119  fill(muN.begin(), muN.end() - numCom, 0.0);
3120  // fill(mux.begin(), mux.end() - numCom, 0.0);
3121 
3122  CalMuPXLBC_partial();
3123 }
3124 
3125 void MixtureComp::CalMuPXLBC_partial()
3126 {
3127  // Here!
3128  // dmu/dP = dmu/dP
3129  // dmu/dNk = dmu/dNk = 0
3130  // // dxij / dP = dxij / dNk = 0
3131  // See MixtureComp::CalMuPXLBC_full01()
3132 
3133  OCP_DBL val1IJ, val2IJ;
3134  OCP_DBL der1IJ, der2IJ, der3J, der4J, der6J, der7J, der8J;
3135  OCP_DBL Tri, tmp;
3136  OCP_DBL xTj, xPj, xVj;
3137  OCP_DBL derxTj, derxPj, derMWj;
3138 
3139  // dxij / dP = 0, d MJ / dP = 0
3140  // Calculate dmuj / dP
3141  // der2IJ = der3J = der4J = der6J = 0;
3142 
3143  for (USI j = 0; j < NP; j++) {
3144  const USI j1 = phaseLabel[j];
3145  const vector<OCP_DBL>& xj = x[j];
3146  const vector<OCP_DBL>& muAuxj = muAux[j];
3147 
3148  xTj = xPj = xVj = 0;
3149  for (USI i = 0; i < NC; i++) {
3150  xTj += xj[i] * Tc[i];
3151  xPj += xj[i] * Pc[i];
3152  xVj += xj[i] * Vcvis[i];
3153  }
3154  der7J = xVj * xiP[j1];
3155  if (muAuxj[3] <= 0.18 && OCP_FALSE) {
3156  muP[j1] = (2.05 * 1E-4) * der7J / muAuxj[2];
3157  } else {
3158  der8J = der7J * (LBCcoef[1] +
3159  muAuxj[3] * (2 * LBCcoef[2] +
3160  muAuxj[3] * (3 * LBCcoef[3] +
3161  muAuxj[3] * 4 * LBCcoef[4])));
3162  muP[j1] = (4 * 1E-4) * pow(muAuxj[4], 3) * der8J / muAuxj[2];
3163  }
3164 
3165  // Calculate dmuj / xkj
3166  const USI bId = numCom * j1;
3167  for (USI k = 0; k < NC; k++) {
3168  derxTj = Tc[k];
3169  derxPj = Pc[k];
3170  derMWj = MWC[k];
3171  der3J = 0;
3172  for (USI i = 0; i < NC; i++) {
3173  val1IJ = muAux1I[i] / sqrt(MW[j]);
3174  der1IJ = -(1 / 2) * muAux1I[i] * pow(MW[j], -1.5) * derMWj;
3175  Tri = T / Tc[i];
3176  if (Tri <= 1.5) {
3177  tmp = 34 * 1E-5 * pow(Tri, 0.94);
3178  } else {
3179  tmp = 17.78 * 1E-5 * pow(4.58 * Tri - 1.67, 0.625);
3180  }
3181  val2IJ = tmp / val1IJ;
3182  der2IJ = -tmp * der1IJ / (val1IJ * val1IJ);
3183  der3J +=
3184  xj[i] * sqrtMWi[i] * der2IJ + delta(i, k) * sqrtMWi[k] * val2IJ;
3185  }
3186  der4J = sqrtMWi[k];
3187  der6J =
3188  5.4402 *
3189  (1.0 / 6 * pow(xTj, -5.0 / 6) * derxTj -
3190  pow(xTj, 1.0 / 6) * (0.5 / MW[j] * derMWj + 2.0 / 3 / xPj * derxPj)) /
3191  (sqrt(MW[j]) * pow(xPj, 2.0 / 3));
3192  der7J = xix[j1 * numCom + k] * xVj + xi[j1] * Vcvis[k];
3193  if (muAuxj[3] <= 0.18 && OCP_FALSE) {
3194  mux[bId + k] =
3195  (der3J * muAuxj[1] - muAuxj[0] * der4J) / (muAuxj[1] * muAuxj[1]) +
3196  2.05 * 1E-4 * (der7J * muAuxj[2] - muAuxj[3] * der6J) /
3197  (muAuxj[2] * muAuxj[2]);
3198  } else {
3199  der8J =
3200  der7J * (LBCcoef[1] +
3201  muAuxj[3] * (2 * LBCcoef[2] +
3202  muAuxj[3] * (3 * LBCcoef[3] +
3203  muAuxj[3] * 4 * LBCcoef[4])));
3204  mux[bId + k] =
3205  (der3J * muAuxj[1] - muAuxj[0] * der4J) / (muAuxj[1] * muAuxj[1]) +
3206  1E-4 *
3207  (4 * pow(muAuxj[4], 3) * der8J * muAuxj[2] -
3208  (pow(muAuxj[4], 4) - 1) * der6J) /
3209  (muAuxj[2] * muAuxj[2]);
3210  }
3211  }
3212  }
3213 }
3214 
3215 void MixtureComp::CalXiPNX_full01()
3216 {
3217  // call CalVfiVfp_full01() before, use rhsDer
3218  // Calculate xiPC, xiNC, xixC
3219  // Calculate xixC first
3220  OCP_DBL CgTP = GAS_CONSTANT * T / P;
3221  OCP_DBL tmp;
3222  for (USI j = 0; j < NP; j++) {
3223  const vector<OCP_DBL>& xj = x[j];
3224  OCP_DBL aj = Aj[j];
3225  OCP_DBL bj = Bj[j];
3226  OCP_DBL zj = Zj[j];
3227 
3228  Bx = Bi;
3229  for (USI i = 0; i < NC; i++) {
3230  tmp = 0;
3231  for (USI k = 0; k < NC; k++) {
3232  tmp += xj[k] * (1 - BIC[i * NC + k]) * sqrt(Ai[i] * Ai[k]);
3233  }
3234  Ax[i] = 2 * tmp;
3235  Zx[i] =
3236  ((bj - zj) * Ax[i] +
3237  ((aj + delta1 * delta2 * (3 * bj * bj + 2 * bj)) +
3238  ((delta1 + delta2) * (2 * bj + 1) - 2 * delta1 * delta2 * bj) * zj -
3239  (delta1 + delta2 - 1) * zj * zj) *
3240  Bx[i]) /
3241  (3 * zj * zj + 2 * ((delta1 + delta2 - 1) * bj - 1) * zj +
3242  (aj + delta1 * delta2 * bj * bj - (delta1 + delta2) * bj * (bj + 1)));
3243  }
3244 
3245  tmp = -xiC[j] * xiC[j] * CgTP;
3246  for (USI i = 0; i < NC; i++) {
3247  xixC[j * NC + i] = tmp * Zx[i];
3248  }
3249  }
3250  // Calculate xiPC and xiNC
3251  if (NP == 1) {
3252  // NP = 1
3253  tmp = -xiC[0] * xiC[0] * CgTP;
3254  xiPC[0] = tmp * (Zp[0] - Zj[0] / P);
3255  for (USI i = 0; i < NC; i++) {
3256  xiNC[i] = tmp * Zn[0][i];
3257  }
3258  } else {
3259  // NP > 1
3260  const OCP_DBL* dnkjdNP = &rhsDer[0];
3261  OCP_DBL dertmp;
3262 
3263  // Calculate xiNC
3264  // Now dnkjdNP reach to dnkj / dNi after CalVfiVfp_full01
3265  for (USI i = 0; i < NC; i++) {
3266  for (USI j = 0; j < NP; j++) {
3267 
3268  dertmp = 0;
3269  tmp = -xiC[j] * xiC[j] * CgTP;
3270  const vector<OCP_DBL>& Znj = Zn[j];
3271 
3272  for (USI k = 0; k < NC; k++) {
3273  dertmp += Znj[k] * dnkjdNP[k];
3274  }
3275  xiNC[j * NC + i] = tmp * dertmp;
3276  dnkjdNP += NC;
3277  }
3278  }
3279  // Calculate xiPC
3280  // Now dnkjdNP reach to dnkj / dP
3281  for (USI j = 0; j < NP; j++) {
3282  tmp = -xiC[j] * xiC[j] * CgTP;
3283  xiPC[j] = Zp[j] - Zj[j] / P;
3284  const vector<OCP_DBL>& Znj = Zn[j];
3285 
3286  // in OCP
3287  for (USI k = 0; k < NC; k++) {
3288  xiPC[j] += Znj[k] * dnkjdNP[k];
3289  }
3290  dnkjdNP += NC;
3291  xiPC[j] *= tmp;
3292  }
3293  }
3294 
3295  // Copoy to xiP, xix, xiN
3296  // fill(xiP.begin(), xiP.end(), 0.0);
3297  // fill(xix.begin(), xix.end(), 0.0);
3298  USI j1;
3299  for (USI j = 0; j < NP; j++) {
3300  j1 = phaseLabel[j];
3301  xiP[j1] = xiPC[j];
3302  Dcopy(NC, &xiN[j1 * numCom], &xiNC[j * NC]);
3303  Dcopy(NC, &xix[j1 * numCom], &xixC[j * NC]);
3304  }
3305 }
3306 
3307 void MixtureComp::CalXiPNX_full02()
3308 {
3309  // Attention!
3310  // NP = 1 or NP = 2
3311 
3312  // Call CalVfiVfp_full02() before, use rhsder
3313  // Calculate xiPC, xiNC, xixC
3314  // Calculate xixC first
3315  OCP_DBL CgTP = GAS_CONSTANT * T / P;
3316  OCP_DBL tmp;
3317  for (USI j = 0; j < NP; j++) {
3318  const vector<OCP_DBL>& xj = x[j];
3319  OCP_DBL aj = Aj[j];
3320  OCP_DBL bj = Bj[j];
3321  OCP_DBL zj = Zj[j];
3322 
3323  Bx = Bi;
3324  for (USI i = 0; i < NC; i++) {
3325  tmp = 0;
3326  for (USI k = 0; k < NC; k++) {
3327  tmp += xj[k] * (1 - BIC[i * NC + k]) * sqrt(Ai[i] * Ai[k]);
3328  }
3329  Ax[i] = 2 * tmp;
3330  Zx[i] =
3331  ((bj - zj) * Ax[i] +
3332  ((aj + delta1 * delta2 * (3 * bj * bj + 2 * bj)) +
3333  ((delta1 + delta2) * (2 * bj + 1) - 2 * delta1 * delta2 * bj) * zj -
3334  (delta1 + delta2 - 1) * zj * zj) *
3335  Bx[i]) /
3336  (3 * zj * zj + 2 * ((delta1 + delta2 - 1) * bj - 1) * zj +
3337  (aj + delta1 * delta2 * bj * bj - (delta1 + delta2) * bj * (bj + 1)));
3338  }
3339 
3340  tmp = -xiC[j] * xiC[j] * CgTP;
3341  for (USI i = 0; i < NC; i++) {
3342  xixC[j * NC + i] = tmp * Zx[i];
3343  }
3344  }
3345 
3346  // Calculate xiPC and xiNC
3347  if (NP == 1) {
3348  // NP = 1
3349  tmp = -xiC[0] * xiC[0] * CgTP;
3350  xiPC[0] = tmp * (Zp[0] - Zj[0] / P);
3351  for (USI i = 0; i < NC; i++) {
3352  xiNC[i] = tmp * Zn[0][i];
3353  }
3354  } else {
3355  // NP = 2
3356  OCP_DBL* nijPN = &rhsDer[0];
3357  // Calculate xiPC
3358  xiPC[0] = (Zp[0] - Zj[0] / P);
3359  xiPC[1] = (Zp[1] - Zj[1] / P);
3360  for (USI i = 0; i < NC; i++) {
3361  xiPC[0] += Zn[0][i] * nijPN[i];
3362  xiPC[1] -= Zn[1][i] * nijPN[i];
3363  }
3364  nijPN += NC;
3365  xiPC[0] *= -xiC[0] * xiC[0] * CgTP;
3366  xiPC[1] *= -xiC[1] * xiC[1] * CgTP;
3367  // Calculate xiNC
3368  OCP_DBL tmp0 = -xiC[0] * xiC[0] * CgTP;
3369  OCP_DBL tmp1 = -xiC[1] * xiC[1] * CgTP;
3370  for (USI i = 0; i < NC; i++) {
3371  xiNC[i] = 0;
3372  xiNC[NC + i] = 0;
3373  for (USI k = 0; k < NC; k++) {
3374  xiNC[i] += Zn[0][k] * nijPN[k];
3375  xiNC[NC + i] -= Zn[1][k] * nijPN[k];
3376  }
3377  xiNC[NC + i] += Zn[1][i];
3378  nijPN += NC;
3379  xiNC[i] *= tmp0;
3380  xiNC[NC + i] *= tmp1;
3381  }
3382  }
3383 
3384  // Copoy to xiP, xix, xiN
3385  // fill(xiP.begin(), xiP.end(), 0.0);
3386  // fill(xix.begin(), xix.end(), 0.0);
3387  USI j1;
3388  for (USI j = 0; j < NP; j++) {
3389  j1 = phaseLabel[j];
3390  xiP[j1] = xiPC[j];
3391  Dcopy(NC, &xiN[j1 * numCom], &xiNC[j * NC]);
3392  Dcopy(NC, &xix[j1 * numCom], &xixC[j * NC]);
3393  }
3394 }
3395 
3396 void MixtureComp::CalRhoPNX_full01()
3397 {
3398  // CaldXsdXp01() shoule be called before
3399  // Using rhsDer
3400  // Cal rhoP, rhoX
3401  // Water has been calculated
3402  // fill(rhox.begin(), rhox.end() - numCom, 0.0);
3403  // fill(rhoP.begin(), rhoP.end() - 1, 0.0);
3404 
3405  USI j1, bId;
3406  for (USI j = 0; j < NP; j++) {
3407  j1 = phaseLabel[j];
3408  bId = j1 * numCom;
3409  rhoP[j1] = xiP[j1] * MW[j];
3410  for (USI m = 0; m < NC; m++) {
3411  rhox[bId + m] = xix[bId + m] * MW[j] + xi[j1] * MWC[m];
3412  }
3413  }
3414  if (NP == 1) {
3415  OCP_DBL tmp;
3416  j1 = phaseLabel[0];
3417  bId = j1 * numCom;
3418  for (USI m = 0; m < NC; m++) {
3419  rhoN[bId + m] = xiN[bId + m] * MW[0];
3420  tmp = MWC[m] * Nh;
3421  for (USI i = 0; i < NC; i++) {
3422  tmp -= MWC[i] * Ni[i];
3423  }
3424  rhoN[bId + m] += tmp * xi[j1] / (Nh * Nh);
3425  }
3426  }
3427  if (NP > 1) {
3428  // correct rhoP
3429  // use rhsDer, see CaldXsdXp01(), attention that it only contains hydrocarbon
3430  OCP_DBL tmp;
3431  const OCP_DBL* xijP = &rhsDer[NP];
3432  for (USI j = 0; j < NP; j++) {
3433  j1 = phaseLabel[j];
3434  tmp = 0;
3435  for (USI i = 0; i < NC; i++) {
3436  tmp += MWC[i] * xijP[i];
3437  }
3438  rhoP[j1] += xi[j1] * tmp;
3439  xijP += NC;
3440  }
3441  // Calculate rhoN
3442  const OCP_DBL* xijN = xijP;
3443  for (USI m = 0; m < NC; m++) {
3444  xijN += NP;
3445  for (USI j = 0; j < NP; j++) {
3446  j1 = phaseLabel[j];
3447  bId = j1 * numCom;
3448  rhoN[bId + m] = xiN[bId + m] * MW[j];
3449  tmp = 0;
3450  for (USI i = 0; i < NC; i++) {
3451  tmp += MWC[i] * xijN[i];
3452  }
3453  rhoN[bId + m] += tmp * xi[j1];
3454  xijN += NC;
3455  }
3456  }
3457  }
3458 }
3459 
3460 void MixtureComp::CalRhoPNX_full()
3461 {
3462  // CaldXsdXpAPI01() or CaldXsdXpAPI02() shoule be called before
3463  // Using dXsdXp
3464  // Cal rhoP, rhoX
3465  // Water has been calculated
3466  // fill(rhox.begin(), rhox.end() - numCom, 0.0);
3467  // fill(rhoP.begin(), rhoP.end() - 1, 0.0);
3468 
3469  USI j1, bId;
3470  OCP_DBL tmp = 0;
3471  if (NP == 1) {
3472  j1 = phaseLabel[0];
3473  bId = j1 * numCom;
3474  rhoP[j1] = xiP[j1] * MW[0];
3475  for (USI m = 0; m < NC; m++) {
3476  rhox[bId + m] = xix[bId + m] * MW[0] + xi[j1] * MWC[m];
3477  rhoN[bId + m] = xiN[bId + m] * MW[0];
3478  tmp = MWC[m] * Nh;
3479  for (USI i = 0; i < NC; i++) {
3480  tmp -= MWC[i] * Ni[i];
3481  }
3482  rhoN[bId + m] += tmp * xi[j1] / (Nh * Nh);
3483  }
3484  } else {
3485  // Using dXsdXp
3486  fill(rhoP.begin(), rhoP.end() - 1, 0.0);
3487  fill(rhoN.begin(), rhoN.end() - numCom, 0.0);
3488  OCP_DBL ncol = numCom + 1;
3489  for (USI j = 0; j < NP; j++) {
3490  j1 = phaseLabel[j];
3491  bId = j1 * numCom;
3492  const OCP_DBL* xijPN = &dXsdXp[(numPhase + bId) * ncol];
3493 
3494  for (USI i = 0; i < NC; i++) {
3495  rhoP[j1] += MWC[i] * xijPN[0];
3496  xijPN++;
3497  for (USI m = 0; m < NC; m++) {
3498  rhoN[bId + m] += MWC[i] * xijPN[m];
3499  }
3500  xijPN += numCom;
3501  }
3502  rhoP[j1] *= xi[j1];
3503  rhoP[j1] += xiP[j1] * MW[j];
3504  for (USI m = 0; m < NC; m++) {
3505  rhoN[bId + m] *= xi[j1];
3506  rhoN[bId + m] += xiN[bId + m] * MW[j];
3507  rhox[bId + m] = xix[bId + m] * MW[j] + xi[j1] * MWC[m];
3508  }
3509  }
3510  }
3511 }
3512 
3513 void MixtureComp::CalMuPX_full01()
3514 {
3515  fill(muP.begin(), muP.end() - 1, 0.0);
3516  fill(muN.begin(), muN.end() - numCom, 0.0);
3517  fill(mux.begin(), mux.end() - numCom, 0.0);
3518 
3519  CalMuPXLBC_full01();
3520 }
3521 
3522 void MixtureComp::CalMuPXLBC_full01()
3523 {
3524  // CaldXsdXp01() before
3525  // use rhsDer
3526  OCP_DBL val1IJ, val2IJ;
3527  OCP_DBL der1IJ, der2IJ, der3J, der4J, der6J, der7J, der8J;
3528  OCP_DBL Tri, tmp;
3529  OCP_DBL xTj, xPj, xVj;
3530  OCP_DBL derxTj, derxPj, derMWj, derxVj;
3531 
3532  // Calculate dmuj / dxkj
3533  for (USI j = 0; j < NP; j++) {
3534  const vector<OCP_DBL>& xj = x[j];
3535  const vector<OCP_DBL>& muAuxj = muAux[j];
3536  const USI j1 = phaseLabel[j];
3537  const USI bId = numCom * j1;
3538  xTj = xPj = xVj = 0;
3539  for (USI i = 0; i < NC; i++) {
3540  xTj += xj[i] * Tc[i];
3541  xPj += xj[i] * Pc[i];
3542  xVj += xj[i] * Vcvis[i];
3543  }
3544  for (USI k = 0; k < NC; k++) {
3545  derxTj = Tc[k];
3546  derxPj = Pc[k];
3547  derMWj = MWC[k];
3548  derxVj = Vcvis[k];
3549  der3J = 0;
3550  for (USI i = 0; i < NC; i++) {
3551  val1IJ = muAux1I[i] / sqrt(MW[j]);
3552  der1IJ = -(1 / 2) * muAux1I[i] * pow(MW[j], -1.5) * derMWj;
3553  Tri = T / Tc[i];
3554  if (Tri <= 1.5) {
3555  tmp = 34 * 1E-5 * pow(Tri, 0.94);
3556  } else {
3557  tmp = 17.78 * 1E-5 * pow(4.58 * Tri - 1.67, 0.625);
3558  }
3559  val2IJ = tmp / val1IJ;
3560  der2IJ = -tmp * der1IJ / (val1IJ * val1IJ);
3561  der3J += sqrtMWi[i] * (delta(i, k) * val2IJ + xj[i] * der2IJ);
3562  }
3563  der4J = sqrtMWi[k];
3564  der6J =
3565  5.4402 *
3566  (1.0 / 6 * pow(xTj, -5.0 / 6) * derxTj -
3567  pow(xTj, 1.0 / 6) * (0.5 / MW[j] * derMWj + 2.0 / 3 / xPj * derxPj)) /
3568  (sqrt(MW[j]) * pow(xPj, 2.0 / 3));
3569  der7J = xix[bId + k] * xVj + xi[j1] * derxVj;
3570  if (muAuxj[3] <= 0.18 && OCP_FALSE) {
3571  mux[bId + k] =
3572  (der3J * muAuxj[1] - muAuxj[0] * der4J) / (muAuxj[1] * muAuxj[1]) +
3573  2.05 * 1E-4 * (der7J * muAuxj[2] - muAuxj[3] * der6J) /
3574  (muAuxj[2] * muAuxj[2]);
3575  } else {
3576  der8J =
3577  der7J * (LBCcoef[1] +
3578  muAuxj[3] * (2 * LBCcoef[2] +
3579  muAuxj[3] * (3 * LBCcoef[3] +
3580  muAuxj[3] * 4 * LBCcoef[4])));
3581  mux[bId + k] =
3582  (der3J * muAuxj[1] - muAuxj[0] * der4J) / (muAuxj[1] * muAuxj[1]) +
3583  1E-4 *
3584  (4 * pow(muAuxj[4], 3) * der8J * muAuxj[2] -
3585  (pow(muAuxj[4], 4) - 1) * der6J) /
3586  (muAuxj[2] * muAuxj[2]);
3587  }
3588  }
3589  }
3590 
3591  if (NP == 1) {
3592  // NP = 1, then dxij / dP = 0, d MJ / dP = 0
3593  // Calculate dmuj / dP
3594  // der2IJ = der3J = der4J = der6J = 0;
3595  const vector<OCP_DBL>& xj = x[0];
3596  const vector<OCP_DBL>& muAuxj = muAux[0];
3597  xTj = xPj = xVj = 0;
3598  for (USI i = 0; i < NC; i++) {
3599  xTj += xj[i] * Tc[i];
3600  xPj += xj[i] * Pc[i];
3601  xVj += xj[i] * Vcvis[i];
3602  }
3603  der7J = xVj * xiPC[0];
3604  if (muAuxj[3] <= 0.18 && OCP_FALSE) {
3605  muP[phaseLabel[0]] = (2.05 * 1E-4) * der7J / muAuxj[2];
3606  } else {
3607  der8J = der7J * (LBCcoef[1] +
3608  muAuxj[3] * (2 * LBCcoef[2] +
3609  muAuxj[3] * (3 * LBCcoef[3] +
3610  muAuxj[3] * 4 * LBCcoef[4])));
3611  muP[phaseLabel[0]] = (4 * 1E-4) * pow(muAuxj[4], 3) * der8J / muAuxj[2];
3612  }
3613  } else {
3614  // NP > 1
3615  // Calculate dmuj / dP
3616  // use rhsDer (after CaldXsdXpAPI01)
3617  for (USI j = 0; j < NP; j++) {
3618  const OCP_DBL* xijP = &rhsDer[NP + j * NC];
3619  const vector<OCP_DBL>& xj = x[j];
3620  const vector<OCP_DBL>& muAuxj = muAux[j];
3621  const USI j1 = phaseLabel[j];
3622  xTj = xPj = xVj = 0;
3623  derxTj = derxPj = derMWj = 0;
3624  for (USI i = 0; i < NC; i++) {
3625  xTj += xj[i] * Tc[i];
3626  xPj += xj[i] * Pc[i];
3627  xVj += xj[i] * Vcvis[i];
3628  derxTj += xijP[i] * Tc[i];
3629  derxPj += xijP[i] * Pc[i];
3630  derMWj += xijP[i] * MWC[i];
3631  }
3632  der3J = der4J = der7J = 0;
3633  for (USI i = 0; i < NC; i++) {
3634  val1IJ = muAux1I[i] / sqrt(MW[j]);
3635  der1IJ = -(1 / 2) * muAux1I[i] * pow(MW[j], -1.5) * derMWj;
3636  Tri = T / Tc[i];
3637  if (Tri <= 1.5) {
3638  tmp = 34 * 1E-5 * pow(Tri, 0.94);
3639  } else {
3640  tmp = 17.78 * 1E-5 * pow(4.58 * Tri - 1.67, 0.625);
3641  }
3642  val2IJ = tmp / val1IJ;
3643  der2IJ = -tmp * der1IJ / (val1IJ * val1IJ);
3644  der3J += sqrtMWi[i] * (xijP[i] * val2IJ + xj[i] * der2IJ);
3645  der4J += sqrtMWi[i] * xijP[i];
3646  der7J += xijP[i] * Vcvis[i];
3647  }
3648  der7J *= xi[j1];
3649  der7J += xiP[j1] * xVj;
3650  der6J =
3651  5.4402 *
3652  (1.0 / 6 * pow(xTj, -5.0 / 6) * derxTj -
3653  pow(xTj, 1.0 / 6) * (0.5 / MW[j] * derMWj + 2.0 / 3 / xPj * derxPj)) /
3654  (sqrt(MW[j]) * pow(xPj, 2.0 / 3));
3655  if (muAuxj[3] <= 0.18 && OCP_FALSE) {
3656  muP[j1] =
3657  (der3J * muAuxj[1] - muAuxj[0] * der4J) / (muAuxj[1] * muAuxj[1]) +
3658  2.05 * 1E-4 * (der7J * muAuxj[2] - muAuxj[3] * der6J) /
3659  (muAuxj[2] * muAuxj[2]);
3660  } else {
3661  der8J =
3662  der7J * (LBCcoef[1] +
3663  muAuxj[3] * (2 * LBCcoef[2] +
3664  muAuxj[3] * (3 * LBCcoef[3] +
3665  muAuxj[3] * 4 * LBCcoef[4])));
3666  muP[j1] =
3667  (der3J * muAuxj[1] - muAuxj[0] * der4J) / (muAuxj[1] * muAuxj[1]) +
3668  1E-4 *
3669  (4 * pow(muAuxj[4], 3) * der8J * muAuxj[2] -
3670  (pow(muAuxj[4], 4) - 1) * der6J) /
3671  (muAuxj[2] * muAuxj[2]);
3672  }
3673  }
3674  }
3675 }
3676 
3677 void MixtureComp::CalXiRhoMuPN_pfullx()
3678 {
3679  // Using dXsdXp
3680  // CalXiPNX_partial(),CalRhoPX_partial(),CalMuPX_partial() have been called
3681  // s: S,x; p: P,N Call CaldXsdXpAPI01 or CaldXsdXpAPI02 before
3682 
3683  const USI ncol = numCom + 1;
3684  const OCP_DBL* xijPN = &dXsdXp[numPhase * ncol];
3685 
3686  if (NP == 1) {
3687  const USI j1 = phaseLabel[0];
3688  const USI bId = j1 * numCom;
3689  xijPN += bId * ncol;
3690 
3691  for (USI i = 0; i < NC; i++) {
3692  xijPN++;
3693  // dN
3694  for (USI m = 0; m < NC; m++) {
3695  xiN[bId + m] += xix[bId + i] * xijPN[m];
3696  rhoN[bId + m] += rhox[bId + i] * xijPN[m];
3697  muN[bId + m] += mux[bId + i] * xijPN[m];
3698  }
3699  xijPN += numCom;
3700  }
3701  } else {
3702  for (USI j = 0; j < numPhase - 1; j++) {
3703 
3704  if (!phaseExist[j]) { // skip
3705  xijPN += numCom * ncol;
3706  continue;
3707  }
3708 
3709  const USI bId = j * numCom;
3710  for (USI i = 0; i < NC; i++) {
3711  // dP
3712  xiP[j] += xix[bId + i] * xijPN[0];
3713  rhoP[j] += rhox[bId + i] * xijPN[0];
3714  muP[j] += mux[bId + i] * xijPN[0];
3715  xijPN++;
3716  // xiN
3717  for (USI m = 0; m < NC; m++) {
3718  xiN[bId + m] += xix[bId + i] * xijPN[m];
3719  rhoN[bId + m] += rhox[bId + i] * xijPN[m];
3720  muN[bId + m] += mux[bId + i] * xijPN[m];
3721  }
3722  xijPN += numCom;
3723  }
3724  // skip water component
3725  xijPN += ncol;
3726  }
3727  }
3728 }
3729 
3730 void MixtureComp::CalXiRhoMuPN_pfullxn(const OCP_BOOL& xflag)
3731 {
3732  // Using dXsdXp
3733  // s: S,n; p: P,N Call CaldXsdXpAPI03, CaldXsdXpAPI02p before
3734  // CalXiPn_partial(); CalRhoPn_partial(); CalMuPn_partial() have been called
3735 
3736  if (NP == 1 && !xflag) {
3737  const USI j1 = phaseLabel[0];
3738  const USI bId = j1 * numCom;
3739 
3740  // dN
3741  for (USI i = 0; i < NC; i++) {
3742  xiN[bId + i] = xix[bId + i];
3743  rhoN[bId + i] = rhox[bId + i];
3744  muN[bId + i] = mux[bId + i];
3745  }
3746  } else {
3747  const USI ncol = numCom + 1;
3748  const OCP_DBL* xijPN = &dXsdXp[(NP + 1) * ncol];
3749 
3750  for (USI j = 0; j < numPhase - 1; j++) {
3751 
3752  if (!phaseExist[j]) // skip
3753  continue;
3754 
3755  const USI bId = j * numCom;
3756  for (USI i = 0; i < NC; i++) {
3757  // dP
3758  xiP[j] += xix[bId + i] * xijPN[0];
3759  rhoP[j] += rhox[bId + i] * xijPN[0];
3760  muP[j] += mux[bId + i] * xijPN[0];
3761  xijPN++;
3762  // xiN
3763  for (USI m = 0; m < NC; m++) {
3764  xiN[bId + m] += xix[bId + i] * xijPN[m];
3765  rhoN[bId + m] += rhox[bId + i] * xijPN[m];
3766  muN[bId + m] += mux[bId + i] * xijPN[m];
3767  }
3768  xijPN += numCom;
3769  }
3770  // don't skip water component
3771  }
3772  }
3773 }
3774 
3775 void MixtureComp::CalVfiVfp_full01()
3776 {
3777  OCP_DBL CgTP = GAS_CONSTANT * T / P;
3778 
3779  if (NP == 1) {
3780  // NP = 1
3781  const OCP_DBL& aj = Aj[0];
3782  const OCP_DBL& bj = Bj[0];
3783  const OCP_DBL& zj = Zj[0];
3784  const vector<OCP_DBL>& xj = x[0];
3785  vector<OCP_DBL>& Znij = Zn[0];
3786  OCP_DBL tmp;
3787 
3788  for (USI i = 0; i < NC; i++) {
3789  tmp = 0;
3790  for (USI m = 0; m < NC; m++) {
3791  tmp += (1 - BIC[i * NC + m]) * sqrt(Ai[i] * Ai[m]) * xj[m];
3792  }
3793  An[i] = 2 / nu[0] * (tmp - aj);
3794  Bn[i] = 1 / nu[0] * (Bi[i] - bj);
3795  Znij[i] =
3796  ((bj - zj) * An[i] +
3797  ((aj + delta1 * delta2 * (3 * bj * bj + 2 * bj)) +
3798  ((delta1 + delta2) * (2 * bj + 1) - 2 * delta1 * delta2 * bj) * zj -
3799  (delta1 + delta2 - 1) * zj * zj) *
3800  Bn[i]) /
3801  (3 * zj * zj + 2 * ((delta1 + delta2 - 1) * bj - 1) * zj +
3802  (aj + delta1 * delta2 * bj * bj - (delta1 + delta2) * bj * (bj + 1)));
3803  vfi[i] = CgTP * (zj + nu[0] * Znij[i]);
3804  }
3805  Zp[0] = ((bj - zj) * aj +
3806  ((aj + delta1 * delta2 * (3 * bj * bj + 2 * bj)) +
3807  ((delta1 + delta2) * (2 * bj + 1) - 2 * delta1 * delta2 * bj) * zj -
3808  (delta1 + delta2 - 1) * zj * zj) *
3809  bj) /
3810  P /
3811  (3 * zj * zj + 2 * ((delta1 + delta2 - 1) * bj - 1) * zj +
3812  (aj + delta1 * delta2 * bj * bj - (delta1 + delta2) * bj * (bj + 1)));
3813  vfP = CgTP * nu[0] * (Zp[0] - zj / P);
3814  } else {
3815  // NP > 1
3816  CalFugNAll();
3817  CalFugPAll();
3818  AssembleMatVfiVfp_full01();
3819  AssembleRhsVfiVfp_full01();
3820  LUSolve(NC + 1, NC * NP, &JmatDer[0], &rhsDer[0], &pivot[0]);
3821  // Calculate Vfi
3822  OCP_DBL* dnkjdNP = &rhsDer[0];
3823  for (USI i = 0; i < NC; i++) {
3824  vfi[i] = 0;
3825  for (USI j = 0; j < NP; j++) {
3826  for (USI k = 0; k < NC; k++) {
3827  vfi[i] += dnkjdNP[j * NC + k] * (Zj[j] + nu[j] * Zn[j][k]);
3828  }
3829  }
3830  vfi[i] *= CgTP;
3831  dnkjdNP += NP * NC;
3832  }
3833  // Calculate Vfp
3834  vfP = 0;
3835  for (USI j = 0; j < NP; j++) {
3836  vfP += nu[j] / P * (Zp[j] * P - Zj[j]);
3837  for (USI k = 0; k < NC; k++) {
3838  vfP += dnkjdNP[j * NC + k] * (Zj[j] + nu[j] * Zn[j][k]);
3839  }
3840  }
3841  vfP *= CgTP;
3842  }
3843 
3844 #ifdef OCP_NANCHECK
3845  if (!CheckNan(vfi.size(), &vfi[0])) {
3846  OCP_ABORT("INF or NAN in vfi !");
3847  }
3848  if (!CheckNan(1, &vfP)) {
3849  OCP_ABORT("INF or NAN in vfP !");
3850  }
3851 #endif // NANCHECK
3852 }
3853 
3854 void MixtureComp::AssembleMatVfiVfp_full01()
3855 {
3856  fill(JmatDer.begin(), JmatDer.end(), 0.0);
3857  // Attention 1: JmatDer should be sorted by column
3858  // Attention 2: d ln fij / d nkj is symetric for each j;
3859  OCP_DBL* bId = &JmatDer[0];
3860  for (USI j = 0; j < NP - 1; j++) {
3861  // for jth phase
3862  OCP_DBL* fugNj = &fugN[j][0];
3863  for (USI i = 0; i < NC; i++) {
3864  // for ith components
3865  Dcopy(NC, bId, fugNj);
3866  bId += (NP - 1 - j) * NC;
3867  bId[i] = 1.0;
3868  bId += (1 + j) * NC;
3869  fugNj += NC;
3870  }
3871  bId += NC;
3872  }
3873  // NP - 1 phase
3874  bId = &JmatDer[(NP - 1) * (NP * NC * NC)];
3875  OCP_DBL* fugNj = &fugN[NP - 1][0];
3876  for (USI i = 0; i < NC; i++) {
3877  for (USI j = 0; j < NP - 1; j++) {
3878  Dcopy(NC, bId, fugNj);
3879  bId += NC;
3880  }
3881  Dscalar((NP - 1) * NC, -1.0, bId - (NP - 1) * NC);
3882  fugNj += NC;
3883  bId[i] = 1.0;
3884  bId += NC;
3885  }
3886 }
3887 
3888 void MixtureComp::AssembleRhsVfiVfp_full01()
3889 {
3890  fill(rhsDer.begin(), rhsDer.end(), 0.0);
3891  OCP_DBL* rhstmp = &rhsDer[0];
3892  for (USI k = 0; k < NC; k++) {
3893  // d Nk
3894  rhstmp[NC * (NP - 1) + k] = 1;
3895  rhstmp += NP * NC;
3896  }
3897  // d P
3898  for (USI j = 0; j < NP; j++) {
3899  for (USI i = 0; i < NC; i++) {
3900  rhstmp[j * NC + i] = fugP[NP - 1][i] - fugP[j][i];
3901  }
3902  }
3903 
3904 #ifdef OCP_NANCHECK
3905  if (!CheckNan(rhsDer.size(), &rhsDer[0])) {
3906  OCP_ABORT("INF or NAN in rhsDer !");
3907  }
3908 #endif // NANCHECK
3909 }
3910 
3911 void MixtureComp::CalVfiVfp_full02()
3912 {
3913  // Attention!
3914  // NP = 1 or NP = 2
3915 
3916  OCP_DBL CgTP = GAS_CONSTANT * T / P;
3917 
3918  if (NP == 1) {
3919  // NP = 1
3920  const OCP_DBL& aj = Aj[0];
3921  const OCP_DBL& bj = Bj[0];
3922  const OCP_DBL& zj = Zj[0];
3923  const vector<OCP_DBL>& xj = x[0];
3924  vector<OCP_DBL>& Znij = Zn[0];
3925  OCP_DBL tmp;
3926 
3927  for (USI i = 0; i < NC; i++) {
3928  tmp = 0;
3929  for (USI m = 0; m < NC; m++) {
3930  tmp += (1 - BIC[i * NC + m]) * sqrt(Ai[i] * Ai[m]) * xj[m];
3931  }
3932  An[i] = 2 / nu[0] * (tmp - aj);
3933  Bn[i] = 1 / nu[0] * (Bi[i] - bj);
3934  Znij[i] =
3935  ((bj - zj) * An[i] +
3936  ((aj + delta1 * delta2 * (3 * bj * bj + 2 * bj)) +
3937  ((delta1 + delta2) * (2 * bj + 1) - 2 * delta1 * delta2 * bj) * zj -
3938  (delta1 + delta2 - 1) * zj * zj) *
3939  Bn[i]) /
3940  (3 * zj * zj + 2 * ((delta1 + delta2 - 1) * bj - 1) * zj +
3941  (aj + delta1 * delta2 * bj * bj - (delta1 + delta2) * bj * (bj + 1)));
3942  vfi[i] = CgTP * (zj + nu[0] * Znij[i]) - Vshift[i];
3943  }
3944  Zp[0] = ((bj - zj) * aj +
3945  ((aj + delta1 * delta2 * (3 * bj * bj + 2 * bj)) +
3946  ((delta1 + delta2) * (2 * bj + 1) - 2 * delta1 * delta2 * bj) * zj -
3947  (delta1 + delta2 - 1) * zj * zj) *
3948  bj) /
3949  P /
3950  (3 * zj * zj + 2 * ((delta1 + delta2 - 1) * bj - 1) * zj +
3951  (aj + delta1 * delta2 * bj * bj - (delta1 + delta2) * bj * (bj + 1)));
3952  vfP = CgTP * nu[0] * (Zp[0] - zj / P);
3953  } else {
3954  // NP = 2, IF NP > 2 -> WRONG!
3955  CalFugNAll();
3956  CalFugPAll();
3957  AssembleMatVfiVfp_full02();
3958  AssembleRhsVfiVfp_full02();
3959  LUSolve(NC + 1, NC, &JmatDer[0], &rhsDer[0], &pivot[0]);
3960  // now d nm0 / dP(dNk) has been available
3961  const OCP_DBL* dnkjdNP = &rhsDer[0];
3962  // Calculate Vfp
3963  const USI j0 = phaseLabel[0];
3964  const USI j1 = phaseLabel[1];
3965  vjp[j0] = CgTP * nu[0] * (Zp[0] - Zj[0] / P);
3966  vjp[j1] = CgTP * nu[1] * (Zp[1] - Zj[1] / P);
3967  for (USI k = 0; k < NC; k++) {
3968  vjp[j0] += (CgTP * (Zj[0] + nu[0] * Zn[0][k]) - Vshift[k]) * dnkjdNP[k];
3969  vjp[j1] -= (CgTP * (Zj[1] + nu[1] * Zn[1][k]) - Vshift[k]) * dnkjdNP[k];
3970  }
3971  vfP = vjp[j0] + vjp[j1];
3972  dnkjdNP += NC;
3973 
3974  // Calculate Vfi
3975  for (USI i = 0; i < NC; i++) {
3976  vfi[i] = 0;
3977  vji[j0][i] = 0;
3978  vji[j1][i] = 0;
3979  for (USI k = 0; k < NC; k++) {
3980  vji[j0][i] +=
3981  (CgTP * (Zj[0] + nu[0] * Zn[0][k]) - Vshift[k]) * dnkjdNP[k];
3982  vji[j1][i] += (CgTP * (Zj[1] + nu[1] * Zn[1][k]) - Vshift[k]) *
3983  (delta(i, k) - dnkjdNP[k]);
3984  }
3985  vfi[i] = vji[j0][i] + vji[j1][i];
3986  dnkjdNP += NC;
3987  }
3988  }
3989 }
3990 
3991 void MixtureComp::AssembleMatVfiVfp_full02()
3992 {
3993  // NP = 2
3994  fill(JmatDer.begin(), JmatDer.end(), 0.0);
3995  // Attention 1: JmatDer should be sorted by column
3996  // Attention 2: d ln fij / d nkj is symetric for each j;
3997  OCP_DBL* tmpMat = &JmatDer[0];
3998  for (USI i = 0; i < NC; i++) {
3999  for (USI m = 0; m < NC; m++) {
4000  tmpMat[m] = fugN[0][i * NC + m] + fugN[1][i * NC + m];
4001  }
4002  tmpMat += NC;
4003  }
4004 }
4005 
4006 void MixtureComp::AssembleRhsVfiVfp_full02()
4007 {
4008  // NP = 2
4009  fill(rhsDer.begin(), rhsDer.end(), 0.0);
4010  OCP_DBL* rhstmp = &rhsDer[0];
4011 
4012  // dP
4013  for (USI i = 0; i < NC; i++) {
4014  rhstmp[i] = fugP[1][i] - fugP[0][i];
4015  }
4016  rhstmp += NC;
4017 
4018  // dNk
4019  for (USI k = 0; k < NC; k++) {
4020  for (USI i = 0; i < NC; i++) {
4021  rhstmp[i] = fugN[1][k * NC + i]; // d lnfij / d nkj = d lnfkj / d nij
4022  }
4023  rhstmp += NC;
4024  }
4025 }
4026 
4027 void MixtureComp::CaldXsdXp01()
4028 {
4029 
4030  // NP > 1 in this function
4031  // only hydrocarbon was considered
4032  // if water exists, vf, vfP, and S should be updated
4033  // CalVfiVfp_full01() and CalXiPNX_full01() should be called before
4034 
4035  fill(JmatTmp.begin(), JmatTmp.end(), 0.0);
4036  OCP_DBL* MTmp = &JmatTmp[0];
4037 
4038  // dF / dXs
4039  // dFn / dXs
4040  for (USI i = 0; i < NC; i++) {
4041  // dFn / dSm
4042  for (USI m = 0; m < NP; m++) {
4043  MTmp[m] = vf * xiC[m] * x[m][i];
4044  }
4045  MTmp += NP;
4046  // dFn / dxnm
4047  for (USI m = 0; m < NP; m++) {
4048  for (USI n = 0; n < NC; n++) {
4049  MTmp[n] =
4050  vf * S[m] * (xixC[m * NC + n] * x[m][i] + delta(i, n) * xiC[m]);
4051  }
4052  MTmp += NC;
4053  }
4054  }
4055  // dFx / dXs
4056  for (USI j = 0; j < NP; j++) {
4057  // dFx / dSm
4058  MTmp += NP;
4059  // dFx / dxnm, m = j
4060  MTmp += j * NC;
4061  for (USI n = 0; n < NC; n++) {
4062  MTmp[n] = 1.0;
4063  }
4064  MTmp += (NP - j) * NC;
4065  }
4066  // dFf / dXs
4067  for (USI j = 0; j < NP - 1; j++) {
4068  const OCP_DBL* fugXJ = &fugX[j][0];
4069  const OCP_DBL* fugXNP = &fugX[NP - 1][0];
4070  for (USI i = 0; i < NC; i++) {
4071  // dFf / dSm
4072  MTmp += NP;
4073  // dFf / dxnm, m = j or m = np-1
4074  // m = j
4075  MTmp += j * NC;
4076  for (USI n = 0; n < NC; n++) {
4077  MTmp[n] = fugXJ[n];
4078  }
4079  fugXJ += NC;
4080  // m = np-1
4081  MTmp += (NP - 1 - j) * NC;
4082  for (USI n = 0; n < NC; n++) {
4083  MTmp[n] = -fugXNP[n];
4084  }
4085  fugXNP += NC;
4086  MTmp += NC;
4087  }
4088  }
4089 
4090  // USI row01 = NP * (NC + 1);
4091  // USI col01 = row01;
4092  // cout << "FsXs" << endl;
4093  // for (USI i = 0; i < row01; i++) {
4094  // for (USI j = 0; j < col01; j++) {
4095  // cout << JmatTmp[i * row01 + j] << " ";
4096  // }
4097  // cout << endl;
4098  //}
4099 
4100  // Transpose JmatTmp
4101  USI dim = NP * (NC + 1);
4102  for (USI i = 0; i < dim; i++) {
4103  for (USI j = 0; j < dim; j++) {
4104  JmatDer[i * dim + j] = JmatTmp[j * dim + i];
4105  }
4106  }
4107 
4108  // dF / dXp
4109  // d fij / dP have been calculated in CalVfiVfp_full01() before
4110  fill(JmatTmp.begin(), JmatTmp.end(), 0.0);
4111  MTmp = &JmatTmp[0];
4112  OCP_DBL tmp01;
4113  OCP_DBL tmp02;
4114  // dFn / dXp
4115  for (USI i = 0; i < NC; i++) {
4116  // dFn / dP
4117  tmp01 = 0;
4118  tmp02 = 0;
4119  for (USI j = 0; j < NP; j++) {
4120  tmp01 += S[j] * x[j][i] * xiPC[j];
4121  tmp02 += S[j] * x[j][i] * xiC[j];
4122  }
4123  MTmp[0] = -(vf * tmp01 + vfP * tmp02); // in OCP
4124  // MTmp[0] = -(vf * tmp01); // in PennSim
4125  MTmp += 1;
4126  // dFn / dNk
4127  // in OCP
4128  for (USI k = 0; k < NC; k++) {
4129  tmp01 = 0;
4130  for (USI j = 0; j < NP; j++) {
4131  tmp01 += S[j] * x[j][i] * xiNC[j * NC + k];
4132  }
4133  MTmp[k] = -(vf * tmp01 + vfi[k] * tmp02 - delta(i, k));
4134  }
4135  // in PennSim
4136  // MTmp[i] = 1.0;
4137  MTmp += NC;
4138  }
4139 
4140  // dFx / dXp
4141  MTmp += NP * (NC + 1);
4142 
4143  // dFf / dXp
4144  const vector<OCP_DBL>& fugPNP = fugP[NP - 1];
4145  for (USI j = 0; j < NP - 1; j++) {
4146  const vector<OCP_DBL>& fugPj = fugP[j];
4147  for (USI i = 0; i < NC; i++) {
4148  // dFf / dP
4149  MTmp[0] = -(fugPj[i] - fugPNP[i]);
4150  MTmp += (NC + 1);
4151  }
4152  }
4153 
4154  // USI row02 = NP * (NC + 1);
4155  // USI col02 = NC + 1;
4156  // cout << "FsXp" << endl;
4157  // for (USI i = 0; i < row02; i++) {
4158  // for (USI j = 0; j < col02; j++) {
4159  // cout << JmatTmp[i * col02 + j] << " ";
4160  // }
4161  // cout << endl;
4162  //}
4163 
4164  // Transpose JmatTmp
4165  USI nrhs = NC + 1;
4166  USI nrow = (NC + 1) * NP;
4167  for (USI i = 0; i < nrhs; i++) {
4168  for (USI j = 0; j < nrow; j++) {
4169  rhsDer[i * nrow + j] = JmatTmp[j * nrhs + i];
4170  }
4171  }
4172  LUSolve(NC + 1, (NC + 1) * NP, &JmatDer[0], &rhsDer[0], &pivot[0]);
4173  // now in rhsDer: dXs / dP, dXs / dNk
4174  // print solution
4175  // cout << "Solution" << endl;
4176  // for (USI i = 0; i < nrhs; i++) {
4177  // for (USI j = 0; j < nrow; j++) {
4178  // cout << rhsDer[i * nrow + j] << " ";
4179  // }
4180  // cout << endl;
4181  //}
4182 }
4183 
4184 void MixtureComp::CaldXsdXpAPI01()
4185 {
4186  // Calculate derivates for hydrocarbon phase and components
4187  // if water exists, vf, vfP, and S should be updated
4188  fill(dXsdXp.begin(), dXsdXp.end(), 0.0);
4189 
4190  OCP_DBL vw = vj[numPhase - 1];
4191  OCP_DBL vwp = vjp[numPhase - 1];
4192 
4193  if (NP == 1) {
4194  USI bId = (numCom + 1) * phaseLabel[0];
4195  OCP_DBL tmp = vf * vf;
4196  // only dSj / dP, dSj / dNk ---- hydrocarbon phase
4197  dXsdXp[bId] = (vfP * vw - vf * vwp) / tmp;
4198  for (USI k = 0; k < NC; k++) {
4199  dXsdXp[bId + k + 1] = (vfi[k] * vw) / tmp;
4200  }
4201  // dxij / dNk --- new
4202  bId = (numCom + 1) * (numPhase + numCom * phaseLabel[0]) + 1;
4203  for (USI i = 0; i < NC; i++) {
4204  for (USI k = 0; k < NC; k++) {
4205  dXsdXp[bId + k] = (delta(i, k) * Nh - Ni[i]) / (Nh * Nh);
4206  }
4207  bId += (numCom + 1);
4208  }
4209  } else {
4210  // NP > 1
4211  // Calculate Saturation
4212  for (USI j = 0; j < NP; j++) {
4213  S[j] = vC[j] / vf;
4214  }
4215  CalFugXAll();
4216  CaldXsdXp01();
4217  // copy from rhsDer to dXsdXp
4218  OCP_DBL* DTmp;
4219  const OCP_DBL* STmp;
4220  USI nrow = NP * (NC + 1); // row num of rhsDer
4221  USI ncol2 = numCom + 1; // col num of dXsdXp
4222  // Saturation
4223  for (USI j = 0; j < NP; j++) {
4224  STmp = &rhsDer[j];
4225  DTmp = &dXsdXp[phaseLabel[j] * ncol2];
4226  // dS / dP
4227  DTmp[0] = STmp[0];
4228  // dS / dNk
4229  DTmp++;
4230  for (USI k = 0; k < NC; k++) {
4231  STmp += nrow;
4232  DTmp[k] = STmp[0];
4233  }
4234  // water is excluded
4235  }
4236  // xij ---- mole fraction of component i in phase j
4237  OCP_DBL* DbId = &dXsdXp[numPhase * ncol2];
4238  for (USI j = 0; j < NP; j++) {
4239  DTmp = DbId + phaseLabel[j] * numCom * ncol2;
4240  for (USI i = 0; i < NC; i++) {
4241  STmp = &rhsDer[NP + j * NC + i];
4242  // dxij / dP
4243  DTmp[0] = STmp[0];
4244  // dxij / dNk
4245  DTmp++;
4246  for (USI k = 0; k < NC; k++) {
4247  STmp += nrow;
4248  DTmp[k] = STmp[0];
4249  }
4250  DTmp += numCom;
4251  // water is excluded
4252  }
4253  }
4254  }
4255  // Correct Sj
4256  CalSaturation();
4257 
4258  USI Wpid = numPhase - 1;
4259  USI Wcid = numCom - 1;
4260 
4261  // Calculate water derivatives
4262  // only dSj / dNw, dSw / dP , dSw / dNk needs to be considered
4263  USI ncol = 1 + numCom;
4264  // dSj / dNw
4265  for (USI j = 0; j < NP; j++) {
4266  dXsdXp[(phaseLabel[j] + 1) * ncol - 1] = -S[phaseLabel[j]] * vfi[Wcid] / vf;
4267  }
4268  OCP_DBL* Dtmp = &dXsdXp[Wpid * ncol];
4269  // dSw / dP
4270  OCP_DBL vf2 = vf * vf;
4271  Dtmp[0] = (vwp * vf - vw * vfP) / vf2;
4272  Dtmp++;
4273  // dSw / d Nk
4274  for (USI k = 0; k < NC; k++) {
4275  Dtmp[k] = -vw * vfi[k] / vf2;
4276  }
4277  // dSw / d Nw
4278  Dtmp[NC] = (vfi[Wcid] * vf - vw * vfi[Wcid]) / vf2;
4279 }
4280 
4281 void MixtureComp::CaldXsdXpAPI02()
4282 {
4283  // Attention!
4284  // NP = 1 or NP = 2
4285 
4286  // dS / dP
4287  // S = Sj, xij
4288  // P = P, Ni
4289  // water is included
4290  fill(dXsdXp.begin(), dXsdXp.end(), 0);
4291  USI ncol = numCom + 1;
4292 
4293  OCP_DBL vf2 = vf * vf;
4294  OCP_DBL vwp = vjp[numPhase - 1];
4295  OCP_DBL vw = vj[numPhase - 1];
4296 
4297  if (NP == 1) {
4298  // when NP = 1, vfi = vhi, vfP = vhp + vwp, h = hydrocarbon
4299  USI j0 = phaseLabel[0];
4300  OCP_DBL* bId = &dXsdXp[j0 * ncol];
4301  // dS / dP
4302  bId[0] = ((vfP - vwp) * vf - vfP * vj[j0]) / vf2;
4303  bId++;
4304  // dSh / dNm
4305  for (USI m = 0; m < NC; m++) {
4306  bId[m] = vfi[m] * vw / vf2;
4307  }
4308  bId += NC;
4309  // dSh / dNw
4310  bId[0] = -vfi[numCom - 1] * vj[j0] / vf2;
4311  // dSw / dP, dNm
4312  bId = &dXsdXp[(numPhase - 1) * ncol];
4313  for (USI m = 0; m < ncol; m++) {
4314  bId[m] = -dXsdXp[j0 * ncol + m];
4315  }
4316  // dxij / dNm
4317  bId = &dXsdXp[(numPhase + j0 * numCom) * ncol + 1];
4318  for (USI i = 0; i < NC; i++) {
4319  for (USI m = 0; m < NC; m++) {
4320  bId[m] = (delta(i, m) * Nh - Ni[i]) / (Nh * Nh);
4321  }
4322  bId += ncol;
4323  }
4324  } else {
4325  // NP = 2
4326  // dSj / dP, dSj / dNm
4327  OCP_DBL* bId;
4328  for (USI j = 0; j < 2; j++) {
4329  const USI j1 = phaseLabel[j];
4330  bId = &dXsdXp[j1 * ncol];
4331  // dSj / dP
4332  bId[0] = (vjp[j1] * vf - vfP * vj[j1]) / vf2;
4333  bId++;
4334  // dSj / dNm
4335  for (USI m = 0; m < numCom; m++) {
4336  bId[m] = (vji[j1][m] * vf - vfi[m] * vj[j1]) / vf2;
4337  }
4338  }
4339  // dSw / dP, dSw / dNm
4340  bId = &dXsdXp[(numPhase - 1) * ncol];
4341  // dSw / dP
4342  bId[0] = (vwp * vf - vfP * vw) / vf2;
4343  bId++;
4344  // dSw / dNm
4345  for (USI m = 0; m < numCom; m++) {
4346  bId[m] = (vji[(numPhase - 1)][m] * vf - vfi[m] * vw) / vf2;
4347  }
4348  bId += numCom;
4349 
4350  // dxij / dP, dxij / dNm
4351  const USI j0 = phaseLabel[0];
4352  const USI j1 = phaseLabel[1];
4353  const OCP_DBL* dnkjdNP = &rhsDer[0];
4354  OCP_DBL* bId1 = bId + j0 * numCom * ncol;
4355  OCP_DBL* bId2 = bId + j1 * numCom * ncol;
4356  OCP_DBL njDerSum = 0;
4357  // dxij / dP
4358  for (USI i = 0; i < NC; i++) njDerSum += dnkjdNP[i];
4359  for (USI i = 0; i < NC; i++) {
4360  bId1[0] = (dnkjdNP[i] * nu[0] - njDerSum * Nh * n[0][i]) / (nu[0] * nu[0]);
4361  bId2[0] = (-dnkjdNP[i] * nu[1] + njDerSum * Nh * n[1][i]) / (nu[1] * nu[1]);
4362  bId1 += ncol;
4363  bId2 += ncol;
4364  }
4365  dnkjdNP += NC;
4366  // dxij / dNm
4367  for (USI m = 0; m < NC; m++) {
4368  njDerSum = 0;
4369  for (USI i = 0; i < NC; i++) njDerSum += dnkjdNP[i];
4370  bId1 = bId + j0 * numCom * ncol + m + 1;
4371  bId2 = bId + j1 * numCom * ncol + m + 1;
4372  for (USI i = 0; i < NC; i++) {
4373  bId1[0] =
4374  (dnkjdNP[i] * nu[0] - njDerSum * Nh * n[0][i]) / (nu[0] * nu[0]);
4375  bId2[0] = ((delta(i, m) - dnkjdNP[i]) * nu[1] -
4376  (1 - njDerSum) * Nh * n[1][i]) /
4377  (nu[1] * nu[1]);
4378  bId1 += ncol;
4379  bId2 += ncol;
4380  }
4381  dnkjdNP += NC;
4382  }
4383  }
4384 }
4385 
4386 void MixtureComp::CaldXsdXpAPI02p()
4387 {
4388  // Attention!
4389  // NP = 1 or NP = 2
4390 
4391  // water not in oil or gas; hydroncarbon not in water
4392  // inexist phase will bot be stored
4393  // only hydrocarbon phases, hydrocarbon components, water phase will be stored
4394 
4395  // dS / dP
4396  // S = Sj, xij
4397  // P = P, Ni
4398  // water is included
4399  fill(dXsdXp.begin(), dXsdXp.end(), 0);
4400  USI ncol = numCom + 1;
4401 
4402  OCP_DBL vf2 = vf * vf;
4403  OCP_DBL vwp = vjp[numPhase - 1];
4404  OCP_DBL vw = vj[numPhase - 1];
4405 
4406  if (NP == 1) {
4407  // when NP = 1, vfi = vhi, vfP = vhp + vwp, h = hydrocarbon
4408  const USI j1 = phaseLabel[0];
4409  OCP_DBL* bId = &dXsdXp[0];
4410  // dS / dP
4411  bId[0] = ((vfP - vwp) * vf - vfP * vj[j1]) / vf2;
4412  bId++;
4413  // dSh / dNm
4414  for (USI m = 0; m < NC; m++) {
4415  bId[m] = vfi[m] * vw / vf2;
4416  }
4417  bId += NC;
4418  // dSh / dNw
4419  bId[0] = -vfi[numCom - 1] * vj[j1] / vf2;
4420  bId++;
4421  // dSw / dP, dNm
4422  for (USI m = 0; m < ncol; m++) {
4423  bId[m] = -dXsdXp[m];
4424  }
4425  bId += ncol + 1;
4426  // dxij / dNm
4427  for (USI i = 0; i < NC; i++) {
4428  for (USI m = 0; m < NC; m++) {
4429  bId[m] = (delta(i, m) * Nh - Ni[i]) / (Nh * Nh);
4430  }
4431  bId += ncol;
4432  }
4433  } else {
4434  // NP = 2
4435  // dSj / dP, dSj / dNm
4436  OCP_DBL* bId;
4437  for (USI j = 0; j < 2; j++) {
4438  const USI j1 = phaseLabel[j];
4439  bId = &dXsdXp[j1 * ncol];
4440  // dSj / dP
4441  bId[0] = (vjp[j1] * vf - vfP * vj[j1]) / vf2;
4442  bId++;
4443  // dSj / dNm
4444  for (USI m = 0; m < numCom; m++) {
4445  bId[m] = (vji[j1][m] * vf - vfi[m] * vj[j1]) / vf2;
4446  }
4447  }
4448  // dSw / dP, dSw / dNm
4449  bId = &dXsdXp[(numPhase - 1) * ncol];
4450  // dSw / dP
4451  bId[0] = (vwp * vf - vfP * vw) / vf2;
4452  bId++;
4453  // dSw / dNm
4454  for (USI m = 0; m < numCom; m++) {
4455  bId[m] = (vji[(numPhase - 1)][m] * vf - vfi[m] * vw) / vf2;
4456  }
4457  bId += numCom;
4458 
4459  // dxij / dP, dxij / dNm
4460  const USI j0 = phaseLabel[0];
4461  const USI j1 = phaseLabel[1];
4462  const OCP_DBL* dnkjdNP = &rhsDer[0];
4463  OCP_DBL* bId1 = bId + j0 * NC * ncol;
4464  OCP_DBL* bId2 = bId + j1 * NC * ncol;
4465  OCP_DBL njDerSum = 0;
4466  // dxij / dP
4467  for (USI i = 0; i < NC; i++) njDerSum += dnkjdNP[i];
4468  for (USI i = 0; i < NC; i++) {
4469  bId1[0] = (dnkjdNP[i] * nu[0] - njDerSum * Nh * n[0][i]) / (nu[0] * nu[0]);
4470  bId2[0] = (-dnkjdNP[i] * nu[1] + njDerSum * Nh * n[1][i]) / (nu[1] * nu[1]);
4471  bId1 += ncol;
4472  bId2 += ncol;
4473  }
4474  dnkjdNP += NC;
4475  // dxij / dNm
4476  for (USI m = 0; m < NC; m++) {
4477  njDerSum = 0;
4478  for (USI i = 0; i < NC; i++) njDerSum += dnkjdNP[i];
4479  bId1 = bId + j0 * NC * ncol + m + 1;
4480  bId2 = bId + j1 * NC * ncol + m + 1;
4481  for (USI i = 0; i < NC; i++) {
4482  bId1[0] =
4483  (dnkjdNP[i] * nu[0] - njDerSum * Nh * n[0][i]) / (nu[0] * nu[0]);
4484  bId2[0] = ((delta(i, m) - dnkjdNP[i]) * nu[1] -
4485  (1 - njDerSum) * Nh * n[1][i]) /
4486  (nu[1] * nu[1]);
4487  bId1 += ncol;
4488  bId2 += ncol;
4489  }
4490  dnkjdNP += NC;
4491  }
4492  }
4493 }
4494 
4495 void MixtureComp::CaldXsdXpAPI03()
4496 {
4497  // p : P, Ni
4498  // s : S, nij
4499  // water not in oil or gas; hydroncarbon not in water
4500  // inexist phase will bot be stored
4501  // only hydrocarbon phases, hydrocarbon components, water phase will be stored
4502 
4503  // water is included
4504  fill(dXsdXp.begin(), dXsdXp.end(), 0.0);
4505  fill(res.begin(), res.end(), 0.0);
4506  resPc = 0;
4507  const USI ncol = numCom + 1;
4508  const OCP_DBL vf2 = vf * vf;
4509  const OCP_DBL vwp = vjp[numPhase - 1];
4510  const OCP_DBL vw = vj[numPhase - 1];
4511 
4512  if (NP == 1) {
4513  // when NP = 1, vfi = vhi, vfP = vhp + vwp, h = hydrocarbon
4514  // and vfi = vji
4515  const USI j0 = phaseLabel[0];
4516  OCP_DBL* bId = &dXsdXp[0];
4517  // dS / dP
4518  bId[0] = ((vfP - vwp) * vf - vfP * vj[j0]) / vf2;
4519  bId++;
4520  // dSh / dNm
4521  for (USI m = 0; m < NC; m++) {
4522  bId[m] = vji[j0][m] * vw / vf2;
4523  }
4524  bId += NC;
4525  // dSh / dNw
4526  bId[0] = -vfi[numCom - 1] * vj[j0] / vf2;
4527  bId++;
4528  // dSw / dP, dNm
4529  for (USI m = 0; m < ncol; m++) {
4530  bId[m] = -dXsdXp[m];
4531  }
4532 
4533  // dnij / dNm
4534  bId = &dXsdXp[2 * ncol + 1];
4535  for (USI i = 0; i < NC; i++) {
4536  bId[i] = 1;
4537  bId += ncol;
4538  }
4539  // res = 0
4540 
4541  } else {
4542  CalFugNAll(OCP_FALSE);
4543  CalFugPAll(OCP_FALSE);
4544  CaldXsdXp03();
4545 
4546  // Assemble dXsdXp
4547  // // inexist phase will bot be stored
4548  // no d xiw / dP, d xiw / dNk
4549  // no d xwj / dP, d xwj / dNk
4550 
4551  vector<USI> sIndex(numPhase, 0); // store index
4552  for (USI j = 0; j < numPhase - 1; j++) {
4553  if (phaseExist[j]) {
4554  for (USI j1 = j + 1; j1 < numPhase; j1++) {
4555  sIndex[j1]++;
4556  }
4557  }
4558  }
4559 
4560  const OCP_DBL* STmp = &rhsDer[0];
4561  const USI nrhs = numCom + 1;
4562  const USI wNP = NP + 1;
4563 
4564  USI bId;
4565  for (USI i = 0; i < nrhs; i++) {
4566  // Sh
4567  for (USI j = 0; j < NP; j++) {
4568  dXsdXp[sIndex[phaseLabel[j]] * nrhs + i] = STmp[j];
4569  }
4570  STmp += NP;
4571  // Sw
4572  dXsdXp[NP * nrhs + i] = STmp[0];
4573  STmp++;
4574  // nij
4575  for (USI j = 0; j < NP; j++) {
4576  bId = wNP + sIndex[phaseLabel[j]] * NC;
4577  for (USI k = 0; k < NC; k++) {
4578  dXsdXp[(bId + k) * nrhs + i] = STmp[k];
4579  }
4580  STmp += NC;
4581  }
4582  }
4583 
4584  // res
4585  // Sh
4586  for (USI j = 0; j < NP; j++) {
4587  res[sIndex[phaseLabel[j]]] = STmp[j];
4588  }
4589  STmp += NP;
4590  // Sw
4591  res[NP] = STmp[0];
4592  STmp++;
4593  // nij
4594  for (USI j = 0; j < NP; j++) {
4595  bId = wNP + sIndex[phaseLabel[j]] * NC;
4596  for (USI k = 0; k < NC; k++) {
4597  res[bId + k] = STmp[k];
4598  }
4599  STmp += NC;
4600  }
4601 
4602  OCP_DBL* myres = &res[NP + 1];
4603  // precalculate a value
4604  for (USI j = 0; j < NP; j++) {
4605  const USI j1 = phaseLabel[j];
4606  for (USI i = 0; i < NC; i++) {
4607  resPc += vji[j1][i] * myres[i];
4608  }
4609  myres += NC;
4610  }
4611 
4612  // resPc = 0;
4613  // fill(res.begin(), res.end(), 0.0);
4614 
4615  // cout << "res1" << endl;
4616  // const USI nrow = wNP + NP * NC;
4617  // for (USI i = 0; i < nrow; i++) {
4618  // for (USI j = nrhs; j < nrhs + 1; j++) {
4619  // cout << setw(12) << rhsDer[j * nrow + i] << " ";
4620  // }
4621  // cout << endl;
4622  // }
4623  // cout << endl;
4624  // cout << "resPc = " << resPc << endl;
4625  }
4626 
4627 #ifdef OCP_NANCHECK
4628  if (!CheckNan(dXsdXp.size(), &dXsdXp[0])) {
4629  OCP_ABORT("INF or INF in bmat !");
4630  }
4631 #endif
4632 }
4633 
4634 void MixtureComp::CaldXsdXp03()
4635 {
4636  // p : P, Ni
4637  // s : S, nij
4638 
4639  fill(JmatTmp.begin(), JmatTmp.end(), 0.0);
4640  OCP_DBL* MTmp = &JmatTmp[0];
4641  const USI dim = NP + 1 + NP * NC;
4642  const OCP_DBL vf2 = vf * vf;
4643  // -dF / ds
4644 
4645  // -dFf / ds
4646  for (USI j = 0; j < NP - 1; j++) {
4647  const OCP_DBL* fugNJ = &fugN[j][0];
4648  const OCP_DBL* fugNNP = &fugN[NP - 1][0];
4649  for (USI i = 0; i < NC; i++) {
4650  // -dFf / dS = 0
4651  MTmp += (NP + 1);
4652  // -dFf / dnij
4653  MTmp += j * NC;
4654  for (USI k = 0; k < NC; k++) {
4655  MTmp[k] = -fugNJ[k];
4656  }
4657  fugNJ += NC;
4658  // np-1 phase
4659  MTmp += (NP - 1 - j) * NC;
4660  for (USI k = 0; k < NC; k++) {
4661  MTmp[k] = fugNNP[k];
4662  }
4663  fugNNP += NC;
4664  MTmp += NC;
4665  }
4666  }
4667  // -dFn / ds
4668  for (USI i = 0; i < NC; i++) {
4669  // -dFn / dS = 0
4670  MTmp += (NP + 1);
4671  // -dFn / dnij
4672  for (USI j = 0; j < NP; j++) {
4673  MTmp[i] = 1;
4674  MTmp += NC;
4675  }
4676  }
4677  // -dFs / ds ---- Hydroncarbon
4678  for (USI j = 0; j < NP; j++) {
4679  const USI j1 = phaseLabel[j];
4680  // -dFs / dS
4681  MTmp[j] = -1;
4682  MTmp += (NP + 1);
4683  // -dFs / dnkm
4684  for (USI m = 0; m < NP; m++) {
4685  const USI m1 = phaseLabel[m];
4686  if (m1 == j1) {
4687  for (USI k = 0; k < NC; k++) {
4688  MTmp[k] = vji[j1][k] * (vf - vj[j1]) / vf2;
4689  }
4690  } else {
4691  for (USI k = 0; k < NC; k++) {
4692  MTmp[k] = -vji[m1][k] * vj[j1] / vf2;
4693  }
4694  }
4695  MTmp += NC;
4696  }
4697  }
4698  // -dFsw / ds
4699  // -dFsw / dSw
4700  MTmp[NP] = -1;
4701  MTmp += (NP + 1);
4702  // -dFsw / dnkm
4703  for (USI m = 0; m < NP; m++) {
4704  const USI m1 = phaseLabel[m];
4705  for (USI k = 0; k < NC; k++) {
4706  MTmp[k] = -vji[m1][k] * vj[numPhase - 1] / vf2;
4707  }
4708  MTmp += NC;
4709  }
4710 
4711  // cout << "dFdS" << endl;
4712  // cout << scientific << setprecision(1);
4713  // for (USI i = 0; i < dim; i++) {
4714  // for (USI j = 0; j < dim; j++)
4715  // cout << JmatTmp[i * dim + j] << " ";
4716  // cout << endl;
4717  // if (i == dim - 4)
4718  // cout << endl;
4719  // }
4720  // cout << phaseLabel[0] << " " << phaseLabel[1] << endl;
4721 
4722  // Transpose JmatTmp
4723  for (USI i = 0; i < dim; i++) {
4724  for (USI j = 0; j < dim; j++) {
4725  JmatDer[i * dim + j] = JmatTmp[j * dim + i];
4726  }
4727  }
4728 
4729  fill(JmatTmp.begin(), JmatTmp.end(), 0.0);
4730  MTmp = &JmatTmp[0];
4731 
4732  // dF / dp
4733  const USI ncol2 = numCom + 1;
4734  // dFf / dp
4735  for (USI j = 0; j < NP - 1; j++) {
4736  for (USI i = 0; i < NC; i++) {
4737  MTmp[0] = fugP[j][i] - fugP[NP - 1][i];
4738  MTmp += ncol2;
4739  }
4740  }
4741  // dFn / dp
4742  for (USI i = 0; i < NC; i++) {
4743  MTmp++;
4744  MTmp[i] = 1;
4745  MTmp += numCom;
4746  }
4747  // dFs / dp ---- Hydroncarbon
4748  for (USI j = 0; j < NP; j++) {
4749  const USI j1 = phaseLabel[j];
4750  // dFs / dP
4751  MTmp[0] = (vfP * vj[j1] - vjp[j1] * vf) / vf2;
4752  // dFs / dNh
4753  MTmp += numCom;
4754  // dFs / dNw
4755  MTmp[0] = vfi[numCom - 1] * vj[j1] / vf2;
4756  MTmp++;
4757  }
4758  // dFsw / dp
4759  // dFsw / dP
4760  MTmp[0] = (vfP * vj[numPhase - 1] - vjp[numPhase - 1] * vf) / vf2;
4761  // dFsw / dNh
4762  MTmp += numCom;
4763  // dFsw / dNw
4764  MTmp[0] = vfi[numCom - 1] * (vj[numPhase - 1] - vf) / vf2;
4765  MTmp++;
4766 
4767  // cout << "dFdp" << endl;
4768  // cout << scientific << setprecision(3);
4769  // for (USI i = 0; i < dim; i++) {
4770  // for (USI j = 0; j < numCom+1; j++)
4771  // cout << JmatTmp[i * (numCom + 1) + j] << " ";
4772  // cout << endl;
4773  // }
4774 
4775  // Transpose JmatTmp
4776  const USI nrhs = numCom + 1;
4777  const USI nrow = NP * NC + NP + 1;
4778  for (USI i = 0; i < nrhs; i++) {
4779  for (USI j = 0; j < nrow; j++) {
4780  rhsDer[i * nrow + j] = JmatTmp[j * nrhs + i];
4781  }
4782  }
4783 
4784  // Calculate res
4785  // resF
4786  OCP_DBL* myRes = &rhsDer[nrhs * nrow];
4787  for (USI j = 0; j < NP - 1; j++) {
4788  const OCP_DBL* fugJ = &fug[j][0];
4789  const OCP_DBL* fugNP = &fug[NP - 1][0];
4790  for (USI i = 0; i < NC; i++) {
4791  myRes[i] = fugJ[i] - fugNP[i];
4792  }
4793  myRes += NC;
4794  }
4795  // resN
4796  for (USI i = 0; i < NC; i++) {
4797  myRes[i] = Ni[i];
4798  for (USI j = 0; j < NP; j++) {
4799  myRes[i] -= x[j][i] * nu[j];
4800  }
4801  }
4802  myRes += NC;
4803  // resS
4804  for (USI j = 0; j < NP; j++) {
4805  myRes[j] = S[phaseLabel[j]] - vj[phaseLabel[j]] / vf;
4806  }
4807  myRes[NP] = S[numPhase - 1] - vj[numPhase - 1] / vf;
4808 
4809  // cout << scientific << setprecision(4);
4810  // cout << "res0" << endl;
4811  // for (USI i = 0; i < nrow; i++) {
4812  // for (USI j = nrhs; j < nrhs + 1; j++) {
4813  // cout << setw(12) << rhsDer[j * nrow + i] << " ";
4814  // }
4815  // cout << endl;
4816  // }
4817  // cout << endl;
4818 
4819  LUSolve(numCom + 1 + 1, nrow, &JmatDer[0], &rhsDer[0], &pivot[0]);
4820 }
4821 
4822 void MixtureComp::CalVfiVfp_full03()
4823 {
4824  // Call CaldXsdXpAPI03() before
4825  // use dXsdXp: s = Sj,nij; p = P,Ni
4826  USI j1;
4827  if (NP == 1) {
4828  // vfP = vfP; vfi = vji
4829  j1 = phaseLabel[0];
4830  for (USI i = 0; i < NC; i++) {
4831  vfi[i] = vji[j1][i];
4832  }
4833  } else {
4834  const USI ncol = numCom + 1;
4835  const OCP_DBL* nijPN = &dXsdXp[(NP + 1) * ncol];
4836  for (USI j = 0; j < numPhase - 1; j++) {
4837 
4838  if (!phaseExist[j]) // skip
4839  continue;
4840 
4841  for (USI i = 0; i < NC; i++) {
4842  vfP += vji[j][i] * nijPN[0];
4843  nijPN++;
4844  for (USI m = 0; m < NC; m++) {
4845  vfi[m] += vji[j][i] * nijPN[m];
4846  }
4847  nijPN += numCom;
4848  }
4849  // don't skip water, no water in dXsdXp here.
4850  }
4851  }
4852 }
4853 
4854 void MixtureComp::CalKeyDerx()
4855 {
4856  // Call CaldXsdXpAPI01 or CaldXsdXpAPI02 before
4857  // Calculate d (xij*xi/mu) / dP or dNk
4858  // no water component
4859  // use dXsdXp
4860  // s = S, xij, p = P, Nk
4861 
4862  fill(keyDer.begin(), keyDer.end(), 0.0);
4863  const USI ncol = numCom + 1;
4864  OCP_DBL tmp;
4865 
4866  const OCP_DBL* sTmp = &dXsdXp[numPhase * ncol];
4867  const OCP_DBL* xiNtmp = &xiN[0];
4868  const OCP_DBL* muNtmp = &muN[0];
4869  const OCP_DBL* xijtmp = &xij[0];
4870  OCP_DBL* dTmp = &keyDer[0];
4871 
4872  // hydrocarbon phase
4873  for (USI j = 0; j < numPhase - 1; j++) {
4874  if (!phaseExist[j]) {
4875  sTmp += numCom * ncol;
4876  xiNtmp += numCom;
4877  muNtmp += numCom;
4878  xijtmp += numCom;
4879  continue;
4880  }
4881  tmp = xi[j] / (mu[j] * mu[j]);
4882  for (USI i = 0; i < NC; i++) {
4883  // dP
4884  dTmp[0] = (sTmp[0] * xi[j] + xijtmp[i] * xiP[j]) / mu[j] -
4885  xijtmp[i] * muP[j] * tmp;
4886  dTmp++;
4887  sTmp++;
4888  // dNk
4889  for (USI k = 0; k < NC; k++) {
4890  dTmp[k] = (sTmp[k] * xi[j] + xijtmp[i] * xiNtmp[k]) / mu[j] -
4891  xijtmp[i] * muNtmp[k] * tmp;
4892  }
4893  dTmp += numCom;
4894  sTmp += numCom;
4895  }
4896  xiNtmp += numCom;
4897  muNtmp += numCom;
4898  xijtmp += numCom;
4899  // skip water components
4900  sTmp += ncol;
4901  }
4902 
4903  // water phase
4904  // dP
4905  const USI wpid = numPhase - 1;
4906  keyDer[NP * NC * ncol] =
4907  (xiP[wpid] * mu[wpid] - muP[wpid] * xi[wpid]) / (mu[wpid] * mu[wpid]);
4908 }
4909 
4910 void MixtureComp::CalKeyDern()
4911 {
4912  // Call CaldXsdXpAPI03 before
4913  // Calculate d (xij*xi/mu) / dP or dNk
4914  // no water component
4915  // use dXsdXp
4916  // s = S, nij, p = P, Nk
4917  // xij = nij / nj
4918 
4919  vector<OCP_DBL> njPN(NC + 1, 0);
4920 
4921  fill(keyDer.begin(), keyDer.end(), 0.0);
4922  const USI ncol = numCom + 1;
4923  const OCP_DBL* sTmp = &dXsdXp[(NP + 1) * ncol];
4924  const OCP_DBL* xiNtmp = &xiN[0];
4925  const OCP_DBL* muNtmp = &muN[0];
4926  const OCP_DBL* xijtmp = &xij[0];
4927  OCP_DBL* dTmp = &keyDer[0];
4928  OCP_DBL tmp01, tmp02;
4929 
4930  if (NP == 1) {
4931  const USI j = phaseLabel[0];
4932  xiNtmp += j * numCom;
4933  muNtmp += j * numCom;
4934  xijtmp += j * numCom;
4935  const OCP_DBL tmp = xi[j] / (mu[j] * mu[j]);
4936 
4937  for (USI i = 0; i < NC; i++) {
4938  // dP
4939  dTmp[0] = (xijtmp[i] * xiP[j]) / mu[j] - xijtmp[i] * muP[j] * tmp;
4940  dTmp++;
4941  // dNk
4942  for (USI k = 0; k < NC; k++) {
4943  dTmp[k] = ((delta(i, k) - xijtmp[i]) * xi[j] / nj[j] +
4944  xijtmp[i] * xiNtmp[k]) /
4945  mu[j] -
4946  xijtmp[i] * muNtmp[k] * tmp;
4947  }
4948  dTmp += numCom;
4949  }
4950  } else {
4951  // hydrocarbon phase
4952  for (USI j = 0; j < numPhase - 1; j++) {
4953  if (!phaseExist[j]) {
4954  xiNtmp += numCom;
4955  muNtmp += numCom;
4956  xijtmp += numCom;
4957  continue;
4958  }
4959 
4960  fill(njPN.begin(), njPN.end(), 0.0);
4961  for (USI i = 0; i < NC; i++) {
4962  for (USI k = 0; k < NC + 1; k++) {
4963  njPN[k] += sTmp[k];
4964  }
4965  sTmp += ncol;
4966  }
4967  sTmp -= NC * ncol;
4968 
4969  tmp01 = xi[j] / (mu[j] * mu[j]);
4970  for (USI i = 0; i < NC; i++) {
4971  tmp02 = (sTmp[0] - njPN[0] * xijtmp[i]) / nj[j];
4972  // dP
4973  dTmp[0] = (tmp02 * xi[j] + xijtmp[i] * xiP[j]) / mu[j] -
4974  xijtmp[i] * muP[j] * tmp01;
4975  dTmp++;
4976  sTmp++;
4977  // dNk
4978  for (USI k = 0; k < NC; k++) {
4979  tmp02 = (sTmp[k] - njPN[k + 1] * xijtmp[i]) / nj[j];
4980  dTmp[k] = (tmp02 * xi[j] + xijtmp[i] * xiNtmp[k]) / mu[j] -
4981  xijtmp[i] * muNtmp[k] * tmp01;
4982  }
4983  dTmp += numCom;
4984  sTmp += numCom;
4985  }
4986  xiNtmp += numCom;
4987  muNtmp += numCom;
4988  xijtmp += numCom;
4989  }
4990  }
4991 
4992  // water phase
4993  // dP
4994  const USI wpid = numPhase - 1;
4995  keyDer[NP * NC * ncol] =
4996  (xiP[wpid] * mu[wpid] - muP[wpid] * xi[wpid]) / (mu[wpid] * mu[wpid]);
4997 }
4998 
4999 USI MixtureComp::CubicRoot(const OCP_DBL& a,
5000  const OCP_DBL& b,
5001  const OCP_DBL& c,
5002  const OCP_BOOL& NTflag) const
5003 {
5004 
5005  OCP_DBL Q = (a * a - 3 * b) / 9;
5006  OCP_DBL R = (2 * a * a * a - 9 * a * b + 27 * c) / 54;
5007 
5008  OCP_DBL Q3 = Q * Q * Q;
5009  OCP_DBL M = R * R - Q3;
5010 
5011  if (M <= 0) {
5012  // 3 real roots
5013  OCP_DBL theta = acos(R / sqrt(Q3));
5014  Ztmp[0] = -2 * sqrt(Q) * cos(theta / 3) - a / 3;
5015  Ztmp[1] = -2 * sqrt(Q) * cos((theta + 2 * PI) / 3) - a / 3;
5016  Ztmp[2] = -2 * sqrt(Q) * cos((theta - 2 * PI) / 3) - a / 3;
5017 
5018  if (NTflag) {
5019  NTcubicroot(Ztmp[0], a, b, c);
5020  NTcubicroot(Ztmp[1], a, b, c);
5021  NTcubicroot(Ztmp[2], a, b, c);
5022  }
5023 
5024  sort(Ztmp.begin(), Ztmp.end());
5025 
5026  // vector<OCP_DBL> e(3, 0);
5027  // for (USI i = 0; i < 3; i++) {
5028  // e[i] = Ztmp[i] * (Ztmp[i] * (Ztmp[i] + a) + b) + c;
5029  //}
5030  // for (USI i = 0; i < 3; i++) {
5031  // cout << scientific << e[i] << "\t";
5032  //}
5033 
5034  return 3;
5035  } else {
5036  OCP_DBL tmp1 = -R + sqrt(M);
5037  OCP_DBL tmp2 = R + sqrt(M);
5038  OCP_DBL S = signD(tmp1) * pow(fabs(tmp1), 1.0 / 3);
5039  OCP_DBL T = -signD(tmp2) * pow(fabs(tmp2), 1.0 / 3);
5040  Ztmp[0] = S + T - a / 3;
5041 
5042  if (NTflag) {
5043  NTcubicroot(Ztmp[0], a, b, c);
5044  }
5045 
5046  // vector<OCP_DBL> e(1, 0);
5047  // for (USI i = 0; i < 1; i++) {
5048  // e[i] = Ztmp[i] * (Ztmp[i] * (Ztmp[i] + a) + b) + c;
5049  //}
5050  // for (USI i = 0; i < 1; i++) {
5051  // cout << scientific << e[i] << "\t";
5052  //}
5053 
5054  return 1;
5055  }
5056 }
5057 
5059 OCP_DBL signD(const OCP_DBL& d)
5060 {
5061  if (d > 0) {
5062  return 1.0;
5063  } else if (d < 0) {
5064  return -1.0;
5065  } else {
5066  return 0.0;
5067  }
5068 }
5069 
5070 OCP_DBL delta(const USI& i, const USI& j)
5071 {
5072  if (i == j) {
5073  return 1.0;
5074  }
5075  return 0.0;
5076 }
5077 
5078 void NTcubicroot(OCP_DBL& root, const OCP_DBL& a, const OCP_DBL& b, const OCP_DBL& c)
5079 {
5080  OCP_DBL e = root * (root * (root + a) + b) + c;
5081  OCP_DBL df;
5082  OCP_DBL iter = 0;
5083  OCP_DBL optroot = root;
5084  OCP_DBL opte = fabs(e);
5085 
5086  while (fabs(e) > 1E-8) {
5087 
5088  df = root * (3 * root + 2 * a) + b;
5089  root = root - e / df;
5090  iter++;
5091  if (iter > 10) {
5092  // std::cout << "WARNING: INEXACT ROOT FOR CUBIC EQUATIONS" << std::endl;
5093  break;
5094  }
5095  e = root * (root * (root + a) + b) + c;
5096  if (fabs(e) <= opte) {
5097  opte = fabs(e);
5098  optroot = root;
5099  }
5100  }
5101  root = optroot;
5102 }
5103 
5105 // For Output
5107 
5108 void MixtureComp::OutMixtureIters() const
5109 {
5110  cout << "SSMSTA: " << setw(12) << itersSSMSTA << setw(15)
5111  << itersSSMSTA * 1.0 / countsSSMSTA << endl;
5112  cout << "NRSTA: " << setw(12) << itersNRSTA << setw(15)
5113  << itersNRSTA * 1.0 / countsNRSTA << endl;
5114  cout << "SSMSP: " << setw(12) << itersSSMSP << setw(15)
5115  << itersSSMSP * 1.0 / countsSSMSP << endl;
5116  cout << "NRSP: " << setw(12) << itersNRSP << setw(15)
5117  << itersNRSP * 1.0 / countsNRSP << endl;
5118  cout << "NRRR: " << setw(12) << itersRR << setw(15)
5119  << itersRR * 1.0 / countsRR << endl;
5120 }
5121 
5123 // Optional Features
5125 
5126 void MixtureComp::SetupOptionalFeatures(OptionalFeatures& optFeatures,
5127  const OCP_USI& numBulk)
5128 {
5129  skipSta = &optFeatures.skipStaAnaly;
5130  if (skipSta->IfUseSkip()) {
5131  skipSta->Setup(numBulk, numPhase - 1, numCom - 1);
5132  AllocateSkip();
5133  }
5134  misTerm = &optFeatures.miscible;
5135  if (misTerm->IfUseMiscible() == ifUseMiscible == OCP_TRUE) {
5136  misTerm->Setup(numBulk);
5137  } else if (misTerm->IfUseMiscible() != ifUseMiscible) {
5138  OCP_ABORT(
5139  "Keywords MISCIBLE, PARACHOR, MISCSTR aren't been given simultaneously!");
5140  }
5141 }
5142 
5144 // Accelerate PVT
5146 
5147 void MixtureComp::AllocateSkip()
5148 {
5149  phiN.resize(NC * NC);
5150  skipMatSTA.resize(NC * NC);
5151  eigenSkip.resize(NC);
5152  eigenWork.resize(2 * NC + 1);
5153 }
5154 
5155 void MixtureComp::CalPhiNSkip()
5156 {
5157  OCP_DBL C, E, G;
5158  OCP_DBL Cnk, Dnk, Enk, Gnk;
5159  OCP_DBL tmp, aik;
5160 
5161  // 0 th phase
5162  const OCP_DBL& aj = Aj[0];
5163  const OCP_DBL& bj = Bj[0];
5164  const OCP_DBL& zj = Zj[0];
5165  const vector<OCP_DBL>& xj = x[0];
5166  vector<OCP_DBL>& Znj = Zn[0];
5167 
5168  for (USI i = 0; i < NC; i++) {
5169  tmp = 0;
5170  for (USI m = 0; m < NC; m++) {
5171  tmp += (1 - BIC[i * NC + m]) * sqrt(Ai[i] * Ai[m]) * xj[m];
5172  }
5173  An[i] = 2 / nu[0] * (tmp - aj);
5174  Bn[i] = 1 / nu[0] * (Bi[i] - bj);
5175  Znj[i] = ((bj - zj) * An[i] +
5176  ((aj + delta1 * delta2 * (3 * bj * bj + 2 * bj)) +
5177  ((delta1 + delta2) * (2 * bj + 1) - 2 * delta1 * delta2 * bj) * zj -
5178  (delta1 + delta2 - 1) * zj * zj) *
5179  Bn[i]) /
5180  (3 * zj * zj + 2 * ((delta1 + delta2 - 1) * bj - 1) * zj +
5181  (aj + delta1 * delta2 * bj * bj - (delta1 + delta2) * bj * (bj + 1)));
5182  }
5183 
5184  G = (zj + delta1 * bj) / (zj + delta2 * bj);
5185 
5186  for (USI i = 0; i < NC; i++) {
5187  // i th fugacity
5188  C = 1 / (zj - bj);
5189  // D = Bi[i] / bj * (zj - 1);
5190  tmp = 0;
5191  for (USI k = 0; k < NC; k++) {
5192  tmp += (1 - BIC[i * NC + k]) * sqrt(Ai[i] * Ai[k]) * xj[k];
5193  }
5194  E = -aj / ((delta1 - delta2) * bj) * (2 * tmp / aj - Bi[i] / bj);
5195 
5196  for (USI k = 0; k <= i; k++) {
5197  // k th components
5198 
5199  aik = (1 - BIC[i * NC + k]) * sqrt(Ai[i] * Ai[k]);
5200 
5201  Cnk = (Bn[k] - Znj[k]) / ((zj - bj) * (zj - bj));
5202  Dnk = Bi[i] / bj * (Znj[k] - (Bi[k] - bj) * (zj - 1) / (nu[0] * bj));
5203  Gnk = (delta1 - delta2) / ((zj + delta2 * bj) * (zj + delta2 * bj)) *
5204  (Bn[k] * zj - Znj[k] * bj);
5205  /*Enk = 1 / ((delta1 - delta2) * bj * bj) * (An[k] * bj - Bn[k] * aj) *
5206  (Bi[i] / bj - 2 * tmp / aj)
5207  + aj / ((delta1 - delta2) * bj) * (-Bi[i] / (bj * bj) * Bn[k] - 2 /
5208  (aj * aj) * (aj * (aik - tmp) / nu[j] - An[k] * tmp));*/
5209  Enk = -1 / (delta1 - delta2) / (bj * bj) *
5210  (2 * (bj * aik / nu[0] + Bn[k] * (Bi[i] * aj / bj - tmp)) -
5211  An[k] * Bi[i] - aj * Bi[i] / nu[0]);
5212  Enk -= E / nu[0];
5213  phiN[i * NC + k] = 1 / C * Cnk + Dnk + Enk * log(G) + E / G * Gnk;
5214  phiN[k * NC + i] = phiN[i * NC + k];
5215  }
5216  }
5217 }
5218 
5219 void MixtureComp::AssembleSkipMatSTA()
5220 {
5221  // Sysmetric Matrix
5222  // stored by colum
5223  vector<OCP_DBL>& xj = x[0];
5224 
5225  for (USI i = 0; i < NC; i++) {
5226  for (USI j = 0; j <= i; j++) {
5227  skipMatSTA[i * NC + j] =
5228  delta(i, j) + nu[0] * sqrt(xj[i] * xj[j]) * phiN[i * NC + j];
5229  skipMatSTA[j * NC + i] = skipMatSTA[i * NC + j];
5230  }
5231  }
5232 }
5233 
5234 void MixtureComp::CalSkipForNextStep()
5235 {
5236  if (flagSkip && ftype == 0) {
5237  // 1. Np == 1 is base for Skipping
5238  // 2. If flagSkip == true, then next stablity analysis is possible to be
5239  // skipped, it depends on if conditions are met
5240  // 3. If ftype == 0, then the range should be calculated, which also means last
5241  // skip is unsatisfied
5242  CalPhiNSkip();
5243  AssembleSkipMatSTA();
5244 #ifdef DEBUG
5245  if (!CheckNan(skipMatSTA.size(), &skipMatSTA[0])) {
5246  OCP_WARNING("Nan in skipMatSTA!");
5247  }
5248 #endif // DEBUG
5249 
5250  CalEigenSY(NC, &skipMatSTA[0], &eigenSkip[0], &eigenWork[0], 2 * NC + 1);
5251  skipSta->AssignValue(bulkId, eigenSkip[0], P, T, zi);
5252  }
5253  skipSta->SetFlagSkip(bulkId, flagSkip);
5254 }
5255 
5257 // Miscible
5259 
5260 void MixtureComp::InputMiscibleParam(const ComponentParam& param, const USI& tarId)
5261 {
5262  ifUseMiscible = param.miscible;
5263  if (param.Parachor.activity) parachor = param.Parachor.data[tarId];
5264  if (ifUseMiscible && parachor.empty()) {
5265  OCP_ABORT("PARACHOR has not been Input!");
5266  }
5267 }
5268 
5269 void MixtureComp::CalSurfaceTension()
5270 {
5271  // be careful!
5272  // phase molar densities should be converted into gm-M/cc here
5273  if (ifUseMiscible) {
5274  if (NP == 1)
5275  surTen = 100;
5276  else {
5277  const OCP_DBL b0 = xiC[0] * CONV7;
5278  const OCP_DBL b1 = xiC[1] * CONV7;
5279  surTen = 0;
5280  for (USI i = 0; i < NC; i++)
5281  surTen += parachor[i] * (b0 * x[0][i] - b1 * x[1][i]);
5282  surTen = pow(surTen, 4.0);
5283  }
5284  misTerm->AssignValue(bulkId, surTen);
5285  }
5286 }
5287 
5288 /*----------------------------------------------------------------------------*/
5289 /* Brief Change History of This File */
5290 /*----------------------------------------------------------------------------*/
5291 /* Author Date Actions */
5292 /*----------------------------------------------------------------------------*/
5293 /* Shizhe Li Jan/05/2022 Create file */
5294 /*----------------------------------------------------------------------------*/
Dscalar
void Dscalar(const int &n, const double &alpha, double *x)
Scales a vector by a constant.
Definition: DenseMat.cpp:61
CalEigenSY
void CalEigenSY(const int &N, float *A, float *w, float *work, const int &lwork)
Calculate the minimal eigenvalue for symmetric matrix with mkl lapack.
Definition: DenseMat.cpp:14
LUSolve
void LUSolve(const int &nrhs, const int &N, double *A, double *b, int *pivot)
Calls dgesv to solve the linear system for general matrices.
Definition: DenseMat.cpp:221
CheckNan
bool CheckNan(const int &N, const T *x)
check NaN
Definition: DenseMat.hpp:209
Daxpy
void Daxpy(const int &n, const double &alpha, const double *x, double *y)
Constant times a vector plus a vector.
Definition: DenseMat.cpp:68
SYSSolve
int SYSSolve(const int &nrhs, const char *uplo, const int &N, double *A, double *b, int *pivot, double *work, const int &lwork)
Calls dsysy to solve the linear system for symm matrices.
Definition: DenseMat.cpp:234
Dcopy
void Dcopy(const int &N, double *dst, const double *src)
Calculate the minimal eigenvalue for symmetric matrix with mkl lapack.
Definition: DenseMat.cpp:36
Dnorm2
double Dnorm2(const int &N, double *x)
Computes the L2-norm of a vector.
Definition: DenseMat.cpp:55
signD
OCP_DBL signD(const OCP_DBL &d)
Return the sign of double di.
Definition: MixtureComp.cpp:5059
MixtureComp.hpp
MixtureComp class declaration.
EOS_PVTW
const USI EOS_PVTW
Mixture model = equation-of-state.
Definition: OCPConst.hpp:103
USI
unsigned int USI
Generic unsigned integer.
Definition: OCPConst.hpp:23
GAS
const USI GAS
Fluid type = gas.
Definition: OCPConst.hpp:92
OCP_DBL
double OCP_DBL
Double precision.
Definition: OCPConst.hpp:27
CONV1
const OCP_DBL CONV1
1 bbl = CONV1 ft3
Definition: OCPConst.hpp:63
LRATE_MODE
const USI LRATE_MODE
Well option = fixed fluid rate???
Definition: OCPConst.hpp:134
WRATE_MODE
const USI WRATE_MODE
Well option = fixed water rate.
Definition: OCPConst.hpp:133
WATER
const USI WATER
Fluid type = water.
Definition: OCPConst.hpp:93
ORATE_MODE
const USI ORATE_MODE
Well option = fixed oil rate.
Definition: OCPConst.hpp:131
CONV7
const OCP_DBL CONV7
lbm/ft3 -> gm-M/cc
Definition: OCPConst.hpp:69
OIL
const USI OIL
Fluid type = oil.
Definition: OCPConst.hpp:91
OCP_USI
unsigned int OCP_USI
Long unsigned integer.
Definition: OCPConst.hpp:25
CONV5
const OCP_DBL CONV5
0 F = CONV5 R
Definition: OCPConst.hpp:67
PROD
const USI PROD
Well type = producer.
Definition: OCPConst.hpp:123
GAS_CONSTANT
const OCP_DBL GAS_CONSTANT
Gas Constant.
Definition: OCPConst.hpp:37
PI
const OCP_DBL PI
Pi.
Definition: OCPConst.hpp:39
INJ
const USI INJ
Well type = injector.
Definition: OCPConst.hpp:122
BHP_MODE
const USI BHP_MODE
Well option = fixed bottom-hole-pressure.
Definition: OCPConst.hpp:135
OCP_BOOL
unsigned int OCP_BOOL
OCP_BOOL in OCP.
Definition: OCPConst.hpp:29
GRATE_MODE
const USI GRATE_MODE
Well option = fixed gas rate.
Definition: OCPConst.hpp:132
OCP_WARNING
#define OCP_WARNING(msg)
Log warning messages.
Definition: UtilError.hpp:39
OCP_ABORT
#define OCP_ABORT(msg)
Abort if critical error happens.
Definition: UtilError.hpp:47
COMP::MW
OCP_DBL MW
Molecular Weight.
Definition: MixtureComp.hpp:117
COMP::OmegaA
OCP_DBL OmegaA
Param A of Components.
Definition: MixtureComp.hpp:120
COMP::Vc
OCP_DBL Vc
VcMW * MW.
Definition: MixtureComp.hpp:119
COMP::VcMW
OCP_DBL VcMW
Critical Volume / MW.
Definition: MixtureComp.hpp:118
COMP::acf
OCP_DBL acf
Acentric Factor.
Definition: MixtureComp.hpp:116
COMP::Tc
OCP_DBL Tc
Critical Temperature.
Definition: MixtureComp.hpp:115
COMP::OmegaB
OCP_DBL OmegaB
Param B of Components.
Definition: MixtureComp.hpp:121
COMP::Vshift
OCP_DBL Vshift
shift volume
Definition: MixtureComp.hpp:122
COMP::Pc
OCP_DBL Pc
Critical Pressure.
Definition: MixtureComp.hpp:114
COMP::name
string name
Name of components.
Definition: MixtureComp.hpp:113
ComponentParam
ComponentParam contains information of components.
Definition: ParamReservoir.hpp:88
ComponentParam::Tc
Type_A_r< vector< OCP_DBL > > Tc
Critical temperature of hydrocarbon components.
Definition: ParamReservoir.hpp:122
ComponentParam::Vcvis
Type_A_r< vector< OCP_DBL > > Vcvis
Critical volume used for viscosity calculations only.
Definition: ParamReservoir.hpp:134
ComponentParam::numPhase
USI numPhase
num of phase, water is excluded, constant now.
Definition: ParamReservoir.hpp:120
ComponentParam::Zc
Type_A_r< vector< OCP_DBL > > Zc
Critical Z-factor of hydrocarbon components.
Definition: ParamReservoir.hpp:125
ComponentParam::SSMparamSP
vector< string > SSMparamSP
Params for Solving Phase Spliting with SSM.
Definition: ParamReservoir.hpp:190
ComponentParam::MW
Type_A_r< vector< OCP_DBL > > MW
Molecular Weight of hydrocarbon components.
Definition: ParamReservoir.hpp:126
ComponentParam::SSMparamSTA
vector< string > SSMparamSTA
Params for Solving Phase Spliting with SSM.
Definition: ParamReservoir.hpp:188
ComponentParam::numCom
USI numCom
num of components, water is excluded.
Definition: ParamReservoir.hpp:119
ComponentParam::Cname
vector< string > Cname
Name of hydrocarbon components.
Definition: ParamReservoir.hpp:121
ComponentParam::Acf
Type_A_r< vector< OCP_DBL > > Acf
Acentric factor of hydrocarbon components.
Definition: ParamReservoir.hpp:127
ComponentParam::Vshift
Type_A_r< vector< OCP_DBL > > Vshift
Volume shift of hydrocarbon components.
Definition: ParamReservoir.hpp:130
ComponentParam::NRparamSP
vector< string > NRparamSP
Params for Solving Phase Spliting with NR.
Definition: ParamReservoir.hpp:191
ComponentParam::Zcvis
Type_A_r< vector< OCP_DBL > > Zcvis
Critical Z-factor used for viscosity calculations only.
Definition: ParamReservoir.hpp:136
ComponentParam::OmegaB
Type_A_r< vector< OCP_DBL > > OmegaB
OMEGA_B of hydrocarbon components.
Definition: ParamReservoir.hpp:129
ComponentParam::NRparamSTA
vector< string > NRparamSTA
Params for Solving Phase Spliting with NR.
Definition: ParamReservoir.hpp:189
ComponentParam::RRparam
vector< string > RRparam
Params for Solving Rachford-Rice equations.
Definition: ParamReservoir.hpp:192
ComponentParam::Vc
Type_A_r< vector< OCP_DBL > > Vc
Critical volume of hydrocarbon components.
Definition: ParamReservoir.hpp:124
ComponentParam::Pc
Type_A_r< vector< OCP_DBL > > Pc
Critical pressure of hydrocarbon components.
Definition: ParamReservoir.hpp:123
ComponentParam::OmegaA
Type_A_r< vector< OCP_DBL > > OmegaA
OMEGA_A of hydrocarbon components.
Definition: ParamReservoir.hpp:128
ComponentParam::BIC
vector< vector< OCP_DBL > > BIC
Binary interaction.
Definition: ParamReservoir.hpp:138
ComponentParam::LBCcoef
vector< OCP_DBL > LBCcoef
LBC coefficients for viscosity calculation.
Definition: ParamReservoir.hpp:137
ComponentParam::Parachor
Type_A_r< vector< OCP_DBL > > Parachor
PARACHOR of hydrocarbon components.
Definition: ParamReservoir.hpp:131
Miscible::IfUseMiscible
OCP_BOOL IfUseMiscible() const
Return ifUseMiscible.
Definition: PhasePermeability.hpp:36
Miscible::AssignValue
void AssignValue(const OCP_USI &n, const OCP_DBL &v)
Assign value to surTen.
Definition: PhasePermeability.hpp:38
Miscible::Setup
void Setup(const OCP_USI &numBulk)
Allocate memory for Miscible term.
Definition: PhasePermeability.cpp:32
MixtureComp::CopyPhase
void CopyPhase()
Copy the basic properties from MixtureComp to Mixture.
Definition: MixtureComp.cpp:2386
MixtureComp::StableSSM
OCP_BOOL StableSSM(const USI &Id)
strict SSM
Definition: MixtureComp.cpp:1492
MixtureComp::FlashIMPEC
void FlashIMPEC(const OCP_DBL &Pin, const OCP_DBL &Tin, const OCP_DBL *Niin, const USI &lastNP, const OCP_DBL *xijin, const OCP_USI &bId) override
Flash calculation with moles of components.
Definition: MixtureComp.cpp:444
MixtureComp::CalFtypeIMPEC
void CalFtypeIMPEC()
Calculate Flash type for IMPEC.
Definition: MixtureComp.hpp:568
MixtureComp::RachfordRice2P
void RachfordRice2P()
Used when NP = 2, improved Rachford-Rice2.
Definition: MixtureComp.cpp:1961
MixtureComp::CalProdWeight
void CalProdWeight(const OCP_DBL &Pin, const OCP_DBL &Tin, const OCP_DBL *Niin, const vector< OCP_DBL > &prodPhase, vector< OCP_DBL > &prodWeight) override
Calculate ProdWeight for PROD well.
Definition: MixtureComp.cpp:830
MixtureComp::CalFtypeFIM
void CalFtypeFIM(const OCP_DBL *Sjin)
Calculate Flash type for FIM.
Definition: MixtureComp.hpp:570
MixtureComp::CalProdRate
void CalProdRate(const OCP_DBL &Pin, const OCP_DBL &Tin, const OCP_DBL *Niin, vector< OCP_DBL > &prodRate) override
Calculate Production rate for PROD well.
Definition: MixtureComp.cpp:855
MixtureComp::parachor
vector< OCP_DBL > parachor
Parachor params of hydrocarbon components.
Definition: MixtureComp.hpp:605
MixtureComp::CubicRoot
USI CubicRoot(const OCP_DBL &a, const OCP_DBL &b, const OCP_DBL &c, const OCP_BOOL &NTflag=OCP_FALSE) const
Result is stored in Ztmp.
Definition: MixtureComp.cpp:4999
MixtureComp::XiPhase
OCP_DBL XiPhase(const OCP_DBL &Pin, const OCP_DBL &Tin, const OCP_DBL *Ziin, const USI &tarPhase) override
return mass density of phase
Definition: MixtureComp.cpp:711
MixtureComp::RachfordRice3
void RachfordRice3()
Used when NP > 2.
Definition: MixtureComp.cpp:2026
MixtureComp::FlashFIMn
void FlashFIMn(const OCP_DBL &Pin, const OCP_DBL &Tin, const OCP_DBL *Niin, const OCP_DBL *Sjin, const OCP_DBL *xijin, const OCP_DBL *njin, const USI *phaseExistin, const USI &lastNP, const OCP_USI &bId) override
Definition: MixtureComp.cpp:575
MixtureComp::Flash
void Flash(const OCP_DBL &Pin, const OCP_DBL &Tin, const OCP_DBL *Niin) override
flash calculation with saturation of phases.
Definition: MixtureComp.cpp:167
MixtureComp::AssembleSkipMatSTA
void AssembleSkipMatSTA()
Assemble matrix to Calculated eigen value used for skipping.
Definition: MixtureComp.cpp:5219
MixtureComp::CalPhiNSkip
void CalPhiNSkip()
Calculate d ln phi[i][j] / d n[k][j].
Definition: MixtureComp.cpp:5155
MixtureComp::skipMatSTA
vector< OCP_SIN > skipMatSTA
Definition: MixtureComp.hpp:583
MixtureComp::SplitBFGS
void SplitBFGS()
Use BFGS to calculate phase splitting.
Definition: MixtureComp.cpp:2045
MixtureComp::flagSkip
OCP_BOOL flagSkip
If ture, then skipping could be try,.
Definition: MixtureComp.hpp:581
MixtureComp::AllocateSkip
void AllocateSkip()
Allocate memory for variables used in skipping stability analysis.
Definition: MixtureComp.cpp:5147
MixtureComp::phiN
vector< OCP_DBL > phiN
d ln phi[i][j] / d n[k][j]
Definition: MixtureComp.hpp:582
MixtureComp::misTerm
Miscible * misTerm
Miscible term pointing to OptionalFeature.
Definition: MixtureComp.hpp:598
MixtureComp::CalFugXSTA
void CalFugXSTA()
Calculate d ln(Fug) / dx for Y.
Definition: MixtureComp.cpp:1638
MixtureComp::skipSta
SkipStaAnaly * skipSta
Skip analysis Term pointing to OptionalFeature.
Definition: MixtureComp.hpp:576
MixtureComp::SetupWellOpt
void SetupWellOpt(WellOpt &opt, const vector< SolventINJ > &sols, const OCP_DBL &Psurf, const OCP_DBL &Tsurf) override
Definition: MixtureComp.cpp:770
MixtureComp::ifUseMiscible
OCP_BOOL ifUseMiscible
Definition: MixtureComp.hpp:600
MixtureComp::FlashFIM
void FlashFIM(const OCP_DBL &Pin, const OCP_DBL &Tin, const OCP_DBL *Niin, const OCP_DBL *Sjin, const USI &lastNP, const OCP_DBL *xijin, const OCP_USI &bId) override
Flash calculation with moles of components and Calculate the derivative.
Definition: MixtureComp.cpp:494
MixtureComp::eigenWork
vector< OCP_SIN > eigenWork
work space for computing eigenvalues with ssyevd_
Definition: MixtureComp.hpp:587
MixtureComp::UpdateXRR
void UpdateXRR()
Update X according to RR.
Definition: MixtureComp.cpp:2030
MixtureComp::RhoPhase
OCP_DBL RhoPhase(const OCP_DBL &Pin, const OCP_DBL &Pbb, const OCP_DBL &Tin, const OCP_DBL *Ziin, const USI &tarPhase) override
return mass density of phase
Definition: MixtureComp.cpp:743
MixtureComp::eigenSkip
vector< OCP_SIN > eigenSkip
eigen values of matrix for skipping Skip Stability Analysis.
Definition: MixtureComp.hpp:586
MixtureComp::RachfordRice2
void RachfordRice2()
Used when NP = 2.
Definition: MixtureComp.cpp:1907
MixtureComp::surTen
OCP_DBL surTen
Surface tension between hydrocarbons phases.
Definition: MixtureComp.hpp:603
MixtureComp::StableSSM01
OCP_BOOL StableSSM01(const USI &Id)
relaxed SSM
Definition: MixtureComp.cpp:1392
MixtureComp::x2n
void x2n()
x[j][i] -> n[j][i]
Definition: MixtureComp.cpp:1202
MixtureComp::SplitNR
void SplitNR()
Use NR to calculate phase splitting.
Definition: MixtureComp.cpp:2057
MixtureComp::CalSkipForNextStep
void CalSkipForNextStep()
Calculate skip info for next step.
Definition: MixtureComp.cpp:5234
Mixture::resPc
OCP_DBL resPc
a precalculated value
Definition: Mixture.hpp:290
Mixture::rhoP
vector< OCP_DBL > rhoP
d rho / dP: numphase
Definition: Mixture.hpp:265
Mixture::numCom
USI numCom
num of components.
Definition: Mixture.hpp:242
Mixture::dXsdXp
vector< OCP_DBL > dXsdXp
derivatives of second variables wrt. primary variables
Definition: Mixture.hpp:275
Mixture::rho
vector< OCP_DBL > rho
mass density of phase: numPhase
Definition: Mixture.hpp:254
Mixture::mu
vector< OCP_DBL > mu
viscosity of phase: numPhase
Definition: Mixture.hpp:256
Mixture::pSderExist
vector< OCP_BOOL > pSderExist
Existence of derivative of phase saturation.
Definition: Mixture.hpp:286
Mixture::numPhase
USI numPhase
num of phases.
Definition: Mixture.hpp:241
Mixture::xix
vector< OCP_DBL > xix
d xi[j] / d x[i][j]: numphase * numCom
Definition: Mixture.hpp:270
Mixture::vj
vector< OCP_DBL > vj
volume of phase: numPhase;
Definition: Mixture.hpp:251
Mixture::bulkId
OCP_USI bulkId
index of current bulk
Definition: Mixture.hpp:240
Mixture::vfP
OCP_DBL vfP
Definition: Mixture.hpp:259
Mixture::rhox
vector< OCP_DBL > rhox
d rho[j] / d x[i][j]: numphase * numCom
Definition: Mixture.hpp:267
Mixture::S
vector< OCP_DBL > S
saturation of phase: numPhase
Definition: Mixture.hpp:250
Mixture::xij
vector< OCP_DBL > xij
Nij / nj: numPhase*numCom.
Definition: Mixture.hpp:253
Mixture::Nt
OCP_DBL Nt
Total moles of Components.
Definition: Mixture.hpp:247
Mixture::nj
vector< OCP_DBL > nj
mole number of phase j
Definition: Mixture.hpp:252
Mixture::Ni
vector< OCP_DBL > Ni
moles of component: numCom
Definition: Mixture.hpp:248
Mixture::vfi
vector< OCP_DBL > vfi
Definition: Mixture.hpp:262
Mixture::phaseExist
vector< OCP_BOOL > phaseExist
existence of phase: numPhase
Definition: Mixture.hpp:249
Mixture::keyDer
vector< OCP_DBL > keyDer
d (xij*xi/mu) / dP or dNk
Definition: Mixture.hpp:291
Mixture::xiP
vector< OCP_DBL > xiP
d xi / dP: numphase
Definition: Mixture.hpp:268
Mixture::vf
OCP_DBL vf
volume of total fluids.
Definition: Mixture.hpp:246
Mixture::pVnumCom
vector< USI > pVnumCom
num of variable components in the phase
Definition: Mixture.hpp:287
Mixture::muP
vector< OCP_DBL > muP
d mu / dP: numPhase
Definition: Mixture.hpp:271
Mixture::mux
vector< OCP_DBL > mux
d mu[j] / d x[i][j]: numphase * numCom
Definition: Mixture.hpp:273
Mixture::res
vector< OCP_DBL > res
residual of a set of equations
Definition: Mixture.hpp:289
Mixture::mixtureType
USI mixtureType
Definition: Mixture.hpp:238
Mixture::Allocate
void Allocate()
Allocate memory for common variables for basic class.
Definition: Mixture.hpp:40
Mixture::xi
vector< OCP_DBL > xi
molar density of phase: numPhase
Definition: Mixture.hpp:255
NRparamSP::curIt
USI curIt
current Iters
Definition: MixtureComp.hpp:80
NRparamSP::tol2
OCP_DBL tol2
tol*tol
Definition: MixtureComp.hpp:76
NRparamSP::maxIt
USI maxIt
Max Iteration.
Definition: MixtureComp.hpp:74
NRparamSP::tol
OCP_DBL tol
Tolerance.
Definition: MixtureComp.hpp:75
NRparamSP::conflag
OCP_BOOL conflag
convergence flag, if converges, conflag = OCP_TRUE
Definition: MixtureComp.hpp:78
NRparamSP::realTol
OCP_DBL realTol
Real tol.
Definition: MixtureComp.hpp:77
NRparamSTA::conflag
OCP_BOOL conflag
convergence flag, if converges, conflag = OCP_TRUE
Definition: MixtureComp.hpp:52
NRparamSTA::tol
OCP_DBL tol
Tolerance.
Definition: MixtureComp.hpp:49
NRparamSTA::realTol
OCP_DBL realTol
Real tol.
Definition: MixtureComp.hpp:51
NRparamSTA::curIt
USI curIt
current Iters
Definition: MixtureComp.hpp:54
NRparamSTA::maxIt
USI maxIt
Max Iteration.
Definition: MixtureComp.hpp:48
NRparamSTA::tol2
OCP_DBL tol2
tol*tol
Definition: MixtureComp.hpp:50
OCPTable::Eval_All
USI Eval_All(const USI &j, const OCP_DBL &val, vector< OCP_DBL > &outdata, vector< OCP_DBL > &slope)
Careful: the memory outdata and slope have not be allocated before.
Definition: OCPTable.cpp:47
OptionalFeatures::miscible
Miscible miscible
Miscible term for Compositional Model.
Definition: OptionalFeatures.hpp:54
OptionalFeatures::skipStaAnaly
SkipStaAnaly skipStaAnaly
Skip Stability Analysis term.
Definition: OptionalFeatures.hpp:47
RRparam::tol2
OCP_DBL tol2
tol*tol
Definition: MixtureComp.hpp:89
RRparam::maxIt
USI maxIt
Max Iteration.
Definition: MixtureComp.hpp:87
RRparam::curIt
USI curIt
current Iters
Definition: MixtureComp.hpp:91
RRparam::tol
OCP_DBL tol
Tolerance.
Definition: MixtureComp.hpp:88
SSMparamSP::curIt
USI curIt
current Iters
Definition: MixtureComp.hpp:67
SSMparamSP::tol
OCP_DBL tol
Tolerance.
Definition: MixtureComp.hpp:62
SSMparamSP::conflag
OCP_BOOL conflag
convergence flag, if converges, conflag = OCP_TRUE
Definition: MixtureComp.hpp:65
SSMparamSP::tol2
OCP_DBL tol2
tol*tol
Definition: MixtureComp.hpp:63
SSMparamSP::realTol
OCP_DBL realTol
Real tol.
Definition: MixtureComp.hpp:64
SSMparamSP::maxIt
USI maxIt
Max Iteration.
Definition: MixtureComp.hpp:61
SSMparamSTA::realTol
OCP_DBL realTol
Real tol.
Definition: MixtureComp.hpp:37
SSMparamSTA::curIt
USI curIt
current Iterations
Definition: MixtureComp.hpp:40
SSMparamSTA::eYt
OCP_DBL eYt
if Yt > 1 + eYt, then single phase is unstable
Definition: MixtureComp.hpp:35
SSMparamSTA::Ktol
OCP_DBL Ktol
tolerace^2 for K
Definition: MixtureComp.hpp:33
SSMparamSTA::tol2
OCP_DBL tol2
tol*tol
Definition: MixtureComp.hpp:36
SSMparamSTA::tol
OCP_DBL tol
Tolerance.
Definition: MixtureComp.hpp:32
SSMparamSTA::conflag
OCP_BOOL conflag
convergence flag, if converges, conflag = OCP_TRUE
Definition: MixtureComp.hpp:38
SSMparamSTA::maxIt
USI maxIt
Max Iteration.
Definition: MixtureComp.hpp:31
SkipStaAnaly::IfUseSkip
OCP_BOOL IfUseSkip() const
Return ifUseSkip.
Definition: AcceleratePVT.hpp:32
SkipStaAnaly::Setup
void Setup(const OCP_USI &numBulk, const USI &np, const USI &nc)
Allocate memory for SkipStaAnaly term.
Definition: AcceleratePVT.cpp:18
SkipStaAnaly::SetFlagSkip
void SetFlagSkip(const OCP_USI &n, const OCP_BOOL &flagSkip)
Set flag for skipping.
Definition: AcceleratePVT.hpp:36
SkipStaAnaly::AssignValue
void AssignValue(const OCP_USI &n, const OCP_DBL &minEigenSkip, const OCP_DBL &PSkip, const OCP_DBL &TSkip, const vector< OCP_DBL > &ziSkip)
Update variables used for determine if skipping will happen.
Definition: AcceleratePVT.cpp:40
Type_A_r::activity
OCP_BOOL activity
If OCP_FALSE, this param is not given.
Definition: ParamReservoir.hpp:82
Type_A_r::data
vector< T > data
Data of param.
Definition: ParamReservoir.hpp:83
Generated on Tue Feb 7 2023 01:49:51 for OpenCAEPoro by doxygen 1.9.1