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