From d9eb3ea30f1799633dd10fc9436e17b87ca39f15 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 12 Sep 2014 13:39:16 +0200 Subject: [PATCH 01/14] Chaning using the new interpolation for viscosity --- opm/core/props/pvt/PvtDead.cpp | 21 ++++++++++++++++++++- opm/core/props/pvt/PvtDead.hpp | 1 + 2 files changed, 21 insertions(+), 1 deletion(-) diff --git a/opm/core/props/pvt/PvtDead.cpp b/opm/core/props/pvt/PvtDead.cpp index 889620bf..6ff187fd 100644 --- a/opm/core/props/pvt/PvtDead.cpp +++ b/opm/core/props/pvt/PvtDead.cpp @@ -41,6 +41,7 @@ namespace Opm // resize the attributes of the object b_.resize(numRegions); viscosity_.resize(numRegions); + inverseBV_.resize(numRegions); for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { Opm::PvdoTable pvdoTable(pvdoKeyword, regionIdx); @@ -54,8 +55,15 @@ namespace Opm for (int i = 0; i < sz; ++i) { bInv[i] = 1.0 / b[i]; } + + std::vector bvInv(sz); + for (int i = 0; i < sz; ++i) { + bvInv[i] = 1.0 / (b[i] * visc[i]); + } + b_[regionIdx] = NonuniformTableLinear(press, bInv); viscosity_[regionIdx] = NonuniformTableLinear(press, visc); + inverseBV_[regionIdx] = NonuniformTableLinear(press, bvInv); } } @@ -67,6 +75,7 @@ namespace Opm // resize the attributes of the object b_.resize(numRegions); viscosity_.resize(numRegions); + inverseBV_.resize(numRegions); for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { Opm::PvdgTable pvdgTable(pvdgKeyword, regionIdx); @@ -81,8 +90,15 @@ namespace Opm for (int i = 0; i < sz; ++i) { bInv[i] = 1.0 / b[i]; } + + std::vector bvInv(sz); + for (int i = 0; i < sz; ++i) { + bvInv[i] = 1.0 / (b[i] * visc[i]); + } + b_[regionIdx] = NonuniformTableLinear(press, bInv); viscosity_[regionIdx] = NonuniformTableLinear(press, visc); + inverseBV_[regionIdx] = NonuniformTableLinear(press, bvInv); } } @@ -136,7 +152,10 @@ namespace Opm // #pragma omp parallel for for (int i = 0; i < n; ++i) { int regionIdx = getTableIndex_(pvtTableIdx, i); - output_mu[i] = viscosity_[regionIdx](p[i]); + double tempInvB = b_[regionIdx](p[i]); + double tempInvBV = inverseBV_[regionIdx](p[i]); + // output_mu[i] = viscosity_[regionIdx](p[i]); + output_mu[i] = tempInvB / tempInvBV; output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]); } std::fill(output_dmudr, output_dmudr + n, 0.0); diff --git a/opm/core/props/pvt/PvtDead.hpp b/opm/core/props/pvt/PvtDead.hpp index 3d01eb2a..42559d5e 100644 --- a/opm/core/props/pvt/PvtDead.hpp +++ b/opm/core/props/pvt/PvtDead.hpp @@ -157,6 +157,7 @@ namespace Opm // table per PVT region. std::vector > b_; std::vector > viscosity_; + std::vector > inverseBV_; }; } From d9acd0f3b83c84f75bf6d5a8c996c6512ef09194 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Sun, 14 Sep 2014 16:24:37 +0200 Subject: [PATCH 02/14] Rewriting the d_mu_d_p for deadGas. --- opm/core/props/pvt/PvtDead.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/opm/core/props/pvt/PvtDead.cpp b/opm/core/props/pvt/PvtDead.cpp index 6ff187fd..33f125f3 100644 --- a/opm/core/props/pvt/PvtDead.cpp +++ b/opm/core/props/pvt/PvtDead.cpp @@ -156,7 +156,10 @@ namespace Opm double tempInvBV = inverseBV_[regionIdx](p[i]); // output_mu[i] = viscosity_[regionIdx](p[i]); output_mu[i] = tempInvB / tempInvBV; - output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]); + // output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]); + // output_dmudp[i] = tempInvB / (tempInvBV * tempInvBV) * inverseBV_[regionIdx].derivative(p[i]); + output_dmudp[i] = (tempInvBV * b_[regionIdx].derivative(p[i]) + - tempInvB * inverseBV_[regionIdx].derivative(p[i])) / (tempInvBV * tempInvBV); } std::fill(output_dmudr, output_dmudr + n, 0.0); From d393d147a423add3430f4d3a240882ee0ed35a82 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Sun, 14 Sep 2014 17:33:43 +0200 Subject: [PATCH 03/14] Updating the calculation of mu and dmudp for PvtDead. --- opm/core/props/pvt/PvtDead.cpp | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/opm/core/props/pvt/PvtDead.cpp b/opm/core/props/pvt/PvtDead.cpp index 33f125f3..6d3d606f 100644 --- a/opm/core/props/pvt/PvtDead.cpp +++ b/opm/core/props/pvt/PvtDead.cpp @@ -56,6 +56,8 @@ namespace Opm bInv[i] = 1.0 / b[i]; } + // TODO: should we change the name of b so that we know it is the + // inverse more explicitly? std::vector bvInv(sz); for (int i = 0; i < sz; ++i) { bvInv[i] = 1.0 / (b[i] * visc[i]); @@ -118,7 +120,10 @@ namespace Opm // #pragma omp parallel for for (int i = 0; i < n; ++i) { int regionIdx = getTableIndex_(pvtTableIdx, i); - output_mu[i] = viscosity_[regionIdx](p[i]); + // output_mu[i] = viscosity_[regionIdx](p[i]); + double tempInvB = b_[regionIdx](p[i]); + double tempInvBV = inverseBV_[regionIdx](p[i]); + output_mu[i] = tempInvB / tempInvBV; } } @@ -133,8 +138,13 @@ namespace Opm // #pragma omp parallel for for (int i = 0; i < n; ++i) { int regionIdx = getTableIndex_(pvtTableIdx, i); - output_mu[i] = viscosity_[regionIdx](p[i]); - output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]); + // output_mu[i] = viscosity_[regionIdx](p[i]); + // output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]); + double tempInvB = b_[regionIdx](p[i]); + double tempInvBV = inverseBV_[regionIdx](p[i]); + output_mu[i] = tempInvB / tempInvBV; + output_dmudp[i] = (tempInvBV * b_[regionIdx].derivative(p[i]) + - tempInvB * inverseBV_[regionIdx].derivative(p[i])) / (tempInvBV * tempInvBV); } std::fill(output_dmudr, output_dmudr + n, 0.0); From 0b69961d9e0139b01ed1132c2e0271b83f5b8b4a Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Wed, 17 Sep 2014 11:15:58 +0200 Subject: [PATCH 04/14] A little more comment about B and b. --- opm/core/props/pvt/PvtDead.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/opm/core/props/pvt/PvtDead.cpp b/opm/core/props/pvt/PvtDead.cpp index 6d3d606f..dc12fa6c 100644 --- a/opm/core/props/pvt/PvtDead.cpp +++ b/opm/core/props/pvt/PvtDead.cpp @@ -58,6 +58,7 @@ namespace Opm // TODO: should we change the name of b so that we know it is the // inverse more explicitly? + // or use the captial B in a exlicit way to show the difference. std::vector bvInv(sz); for (int i = 0; i < sz; ++i) { bvInv[i] = 1.0 / (b[i] * visc[i]); From eadbffa99b10442aca91ed2ba40017661a46913f Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Thu, 18 Sep 2014 12:46:45 +0200 Subject: [PATCH 05/14] Changing the PvtLiveOil to use 1/BV for viscosity. Pssing the simple 2D test. --- opm/core/props/pvt/PvtLiveOil.cpp | 181 ++++++++++++++++++++---------- 1 file changed, 121 insertions(+), 60 deletions(-) diff --git a/opm/core/props/pvt/PvtLiveOil.cpp b/opm/core/props/pvt/PvtLiveOil.cpp index e73a1fdc..1621f1dc 100644 --- a/opm/core/props/pvt/PvtLiveOil.cpp +++ b/opm/core/props/pvt/PvtLiveOil.cpp @@ -46,36 +46,47 @@ namespace Opm const auto saturatedPvto = pvtoTable.getOuterTable(); // OIL, PVTO - saturated_oil_table_[pvtTableIdx].resize(4); + // saturated_oil_table_[pvtTableIdx].resize(4); + // adding a extra colummn to the PVTO to store 1/(B*mu) + saturated_oil_table_[pvtTableIdx].resize(5); const int sz = saturatedPvto->numRows(); - for (int k=0; k<4; ++k) { + // for (int k=0; k<4; ++k) { + for (int k=0; k<5; ++k) { saturated_oil_table_[pvtTableIdx][k].resize(sz); } for (int i=0; igetPressureColumn()[i]; // p saturated_oil_table_[pvtTableIdx][1][i] = 1.0/saturatedPvto->getOilFormationFactorColumn()[i]; // 1/Bo saturated_oil_table_[pvtTableIdx][2][i] = saturatedPvto->getOilViscosityColumn()[i]; // mu_o - saturated_oil_table_[pvtTableIdx][3][i] = saturatedPvto->getGasSolubilityColumn()[i]; // Rs + saturated_oil_table_[pvtTableIdx][3][i] = 1.0 / (saturatedPvto->getOilFormationFactorColumn()[i] * + saturatedPvto->getOilViscosityColumn()[i]); // 1/(Bo*mu_o) + saturated_oil_table_[pvtTableIdx][4][i] = saturatedPvto->getGasSolubilityColumn()[i]; // Rs } undersat_oil_tables_[pvtTableIdx].resize(sz); for (int i=0; inumRows(); undersat_oil_tables_[pvtTableIdx][i][0].resize(tsize); undersat_oil_tables_[pvtTableIdx][i][1].resize(tsize); undersat_oil_tables_[pvtTableIdx][i][2].resize(tsize); + undersat_oil_tables_[pvtTableIdx][i][3].resize(tsize); for (int j=0; jgetPressureColumn()[j]; // p undersat_oil_tables_[pvtTableIdx][i][1][j] = 1.0/undersaturatedPvto->getOilFormationFactorColumn()[j]; // 1/Bo undersat_oil_tables_[pvtTableIdx][i][2][j] = undersaturatedPvto->getOilViscosityColumn()[j]; // mu_o + undersat_oil_tables_[pvtTableIdx][i][3][j] = 1.0 / (undersaturatedPvto->getOilFormationFactorColumn()[j] * + undersaturatedPvto->getOilViscosityColumn()[j]); // 1/(Bo*mu_o) } } // Complete undersaturated tables by extrapolating from existing data // as is done in Eclipse and Mrst + // TODO: check if the following formulations applying to 1/(Bo*mu_o) int iNext = -1; for (int i=0; i > >::size_type sz_t; for (sz_t j=1; j= 2); assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2); double val1 = @@ -408,9 +467,9 @@ namespace Opm press); } else { // Undersaturated case // Interpolate between table sections - int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], maxR); - double w = (maxR - saturated_oil_table_[pvtTableIdx][3][is]) / - (saturated_oil_table_[pvtTableIdx][3][is+1] - saturated_oil_table_[pvtTableIdx][3][is]); + int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], maxR); + double w = (maxR - saturated_oil_table_[pvtTableIdx][4][is]) / + (saturated_oil_table_[pvtTableIdx][4][is+1] - saturated_oil_table_[pvtTableIdx][4][is]); assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2); assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2); double val1 = @@ -435,7 +494,7 @@ namespace Opm { int section; double Rval = linearInterpolation(saturated_oil_table_[pvtTableIdx][0], - saturated_oil_table_[pvtTableIdx][3], + saturated_oil_table_[pvtTableIdx][4], press, section); // derivative with respect to frist component (pressure) if (deriv == 1) { @@ -444,9 +503,9 @@ namespace Opm saturated_oil_table_[pvtTableIdx][item], press); } else { // Undersaturated case - int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r); - double w = (r - saturated_oil_table_[pvtTableIdx][3][is]) / - (saturated_oil_table_[pvtTableIdx][3][is+1] - saturated_oil_table_[pvtTableIdx][3][is]); + int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], r); + double w = (r - saturated_oil_table_[pvtTableIdx][4][is]) / + (saturated_oil_table_[pvtTableIdx][4][is+1] - saturated_oil_table_[pvtTableIdx][4][is]); assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2); assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2); double val1 = @@ -465,7 +524,7 @@ namespace Opm if (Rval <= r ) { // Saturated case return 0; } else { // Undersaturated case - int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r); + int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], r); assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2); assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2); double val1 = @@ -477,7 +536,7 @@ namespace Opm undersat_oil_tables_[pvtTableIdx][is+1][item], press); - double val = (val2 - val1)/(saturated_oil_table_[pvtTableIdx][3][is+1]-saturated_oil_table_[pvtTableIdx][3][is]); + double val = (val2 - val1)/(saturated_oil_table_[pvtTableIdx][4][is+1]-saturated_oil_table_[pvtTableIdx][4][is]); return val; } @@ -489,9 +548,9 @@ namespace Opm press); } else { // Undersaturated case // Interpolate between table sections - int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r); - double w = (r - saturated_oil_table_[pvtTableIdx][3][is]) / - (saturated_oil_table_[pvtTableIdx][3][is+1] - saturated_oil_table_[pvtTableIdx][3][is]); + int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], r); + double w = (r - saturated_oil_table_[pvtTableIdx][4][is]) / + (saturated_oil_table_[pvtTableIdx][4][is+1] - saturated_oil_table_[pvtTableIdx][4][is]); assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2); assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2); double val1 = @@ -524,9 +583,9 @@ namespace Opm saturated_oil_table_[pvtTableIdx][item], press); } else { // Undersaturated case - int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r); - double w = (r - saturated_oil_table_[pvtTableIdx][3][is]) / - (saturated_oil_table_[pvtTableIdx][3][is+1] - saturated_oil_table_[pvtTableIdx][3][is]); + int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], r); + double w = (r - saturated_oil_table_[pvtTableIdx][4][is]) / + (saturated_oil_table_[pvtTableIdx][4][is+1] - saturated_oil_table_[pvtTableIdx][4][is]); assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2); assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2); double val1 = @@ -545,7 +604,7 @@ namespace Opm if (isSat) { // Saturated case return 0; } else { // Undersaturated case - int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r); + int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], r); assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2); assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2); double val1 = @@ -557,35 +616,37 @@ namespace Opm undersat_oil_tables_[pvtTableIdx][is+1][item], press); - double val = (val2 - val1)/(saturated_oil_table_[pvtTableIdx][3][is+1]-saturated_oil_table_[pvtTableIdx][3][is]); + double val = (val2 - val1)/(saturated_oil_table_[pvtTableIdx][4][is+1]-saturated_oil_table_[pvtTableIdx][4][is]); return val; } - } else { - if (isSat) { // Saturated case - return linearInterpolation(saturated_oil_table_[pvtTableIdx][0], - saturated_oil_table_[pvtTableIdx][item], - press); - } else { // Undersaturated case - // Interpolate between table sections - int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r); - double w = (r - saturated_oil_table_[pvtTableIdx][3][is]) / - (saturated_oil_table_[pvtTableIdx][3][is+1] - saturated_oil_table_[pvtTableIdx][3][is]); - assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2); - assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2); - double val1 = - linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0], - undersat_oil_tables_[pvtTableIdx][is][item], - press); - double val2 = - linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0], - undersat_oil_tables_[pvtTableIdx][is+1][item], - press); - double val = val1 + w*(val2 - val1); - return val; + } + // no derivative + else { + if (isSat) { // Saturated case + return linearInterpolation(saturated_oil_table_[pvtTableIdx][0], + saturated_oil_table_[pvtTableIdx][item], + press); + } else { // Undersaturated case + // Interpolate between table sections + int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], r); + double w = (r - saturated_oil_table_[pvtTableIdx][4][is]) / + (saturated_oil_table_[pvtTableIdx][4][is+1] - saturated_oil_table_[pvtTableIdx][4][is]); + assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2); + assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2); + double val1 = + linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0], + undersat_oil_tables_[pvtTableIdx][is][item], + press); + double val2 = + linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0], + undersat_oil_tables_[pvtTableIdx][is+1][item], + press); + double val = val1 + w*(val2 - val1); + return val; + } } } - } } // namespace Opm From 2b18b8f27fd5959489594c832a2e0579fdad2b73 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Thu, 18 Sep 2014 15:28:38 +0200 Subject: [PATCH 06/14] Finshing changing Pvt_liveGas Remains to be verfied. --- opm/core/props/pvt/PvtLiveGas.cpp | 73 ++++++++++++++++++++++--------- 1 file changed, 52 insertions(+), 21 deletions(-) diff --git a/opm/core/props/pvt/PvtLiveGas.cpp b/opm/core/props/pvt/PvtLiveGas.cpp index dd6b437c..aa6d92c1 100644 --- a/opm/core/props/pvt/PvtLiveGas.cpp +++ b/opm/core/props/pvt/PvtLiveGas.cpp @@ -56,28 +56,37 @@ namespace Opm Opm::PvtgTable pvtgTable(pvtgKeyword, pvtTableIdx); // GAS, PVTG - saturated_gas_table_[pvtTableIdx].resize(4); - saturated_gas_table_[pvtTableIdx][0] = pvtgTable.getOuterTable()->getPressureColumn(); - saturated_gas_table_[pvtTableIdx][1] = pvtgTable.getOuterTable()->getGasFormationFactorColumn(); - saturated_gas_table_[pvtTableIdx][2] = pvtgTable.getOuterTable()->getGasViscosityColumn(); - saturated_gas_table_[pvtTableIdx][3] = pvtgTable.getOuterTable()->getOilSolubilityColumn(); + // saturated_gas_table_[pvtTableIdx].resize(4); + // Adding one extra line to PVTG to store 1./(Bg*mu_g) + saturated_gas_table_[pvtTableIdx].resize(5); + saturated_gas_table_[pvtTableIdx][0] = pvtgTable.getOuterTable()->getPressureColumn(); // Pressure + saturated_gas_table_[pvtTableIdx][1] = pvtgTable.getOuterTable()->getGasFormationFactorColumn(); // Bg + saturated_gas_table_[pvtTableIdx][2] = pvtgTable.getOuterTable()->getGasViscosityColumn(); // mu_g + saturated_gas_table_[pvtTableIdx][3] = pvtgTable.getOuterTable()->getGasFormationFactorColumn(); // TODO + saturated_gas_table_[pvtTableIdx][4] = pvtgTable.getOuterTable()->getOilSolubilityColumn(); // Rv int sz = pvtgTable.getOuterTable()->numRows(); undersat_gas_tables_[pvtTableIdx].resize(sz); for (int i=0; i 1/Bg for (int i=0; i Date: Mon, 22 Sep 2014 10:28:23 +0200 Subject: [PATCH 07/14] Using resize to allocate memory for 1/BMu. --- opm/core/props/pvt/PvtLiveGas.cpp | 9 +++++++-- opm/core/props/pvt/PvtLiveOil.cpp | 4 ++-- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/opm/core/props/pvt/PvtLiveGas.cpp b/opm/core/props/pvt/PvtLiveGas.cpp index aa6d92c1..28cec7bf 100644 --- a/opm/core/props/pvt/PvtLiveGas.cpp +++ b/opm/core/props/pvt/PvtLiveGas.cpp @@ -62,7 +62,10 @@ namespace Opm saturated_gas_table_[pvtTableIdx][0] = pvtgTable.getOuterTable()->getPressureColumn(); // Pressure saturated_gas_table_[pvtTableIdx][1] = pvtgTable.getOuterTable()->getGasFormationFactorColumn(); // Bg saturated_gas_table_[pvtTableIdx][2] = pvtgTable.getOuterTable()->getGasViscosityColumn(); // mu_g - saturated_gas_table_[pvtTableIdx][3] = pvtgTable.getOuterTable()->getGasFormationFactorColumn(); // TODO + // The number of the columns + int nColumns = saturated_gas_table_[pvtTableIdx][2].size(); + saturated_gas_table_[pvtTableIdx][3].resize(nColumns); // allocate memory for 1/(Bg*mu_g) + // saturated_gas_table_[pvtTableIdx][3] = pvtgTable.getOuterTable()->getGasFormationFactorColumn(); // TODO saturated_gas_table_[pvtTableIdx][4] = pvtgTable.getOuterTable()->getOilSolubilityColumn(); // Rv int sz = pvtgTable.getOuterTable()->numRows(); @@ -75,7 +78,9 @@ namespace Opm undersat_gas_tables_[pvtTableIdx][i][0] = undersatTable.getOilSolubilityColumn(); // Rv undersat_gas_tables_[pvtTableIdx][i][1] = undersatTable.getGasFormationFactorColumn(); // Bg undersat_gas_tables_[pvtTableIdx][i][2] = undersatTable.getGasViscosityColumn(); // mu_g - undersat_gas_tables_[pvtTableIdx][i][3] = undersatTable.getGasViscosityColumn(); // TODO + int nColumns = undersat_gas_tables_[pvtTableIdx][i][2].size(); + undersat_gas_tables_[pvtTableIdx][i][3].resize(nColumns); // allocate memory for 1/(Bg*mu_g) + // undersat_gas_tables_[pvtTableIdx][i][3] = undersatTable.getGasViscosityColumn(); // TODO } // Bg -> 1/Bg diff --git a/opm/core/props/pvt/PvtLiveOil.cpp b/opm/core/props/pvt/PvtLiveOil.cpp index 1621f1dc..49549d38 100644 --- a/opm/core/props/pvt/PvtLiveOil.cpp +++ b/opm/core/props/pvt/PvtLiveOil.cpp @@ -68,7 +68,7 @@ namespace Opm const auto undersaturatedPvto = pvtoTable.getInnerTable(i); // undersat_oil_tables_[pvtTableIdx][i].resize(3); - // adding a extra colummn to the PVTO to store 1/BV + // adding a extra colummn to the PVTO to store 1/(B*mu) undersat_oil_tables_[pvtTableIdx][i].resize(4); int tsize = undersaturatedPvto->numRows(); undersat_oil_tables_[pvtTableIdx][i][0].resize(tsize); @@ -117,7 +117,7 @@ namespace Opm double mu = (undersat_oil_tables_[pvtTableIdx][i][2].back())*(1.0+0.5*visc)/(1.0-0.5*visc); undersat_oil_tables_[pvtTableIdx][i][2].push_back(mu); - // a try to expolate the 1/BMu + // A try to expolate the 1/BMu with the expolated mu and B double inverseBMu = 1.0 / (B*mu); undersat_oil_tables_[pvtTableIdx][i][3].push_back(inverseBMu); From be26c0af18f423891314240e7fbde2c15f5501d0 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Mon, 22 Sep 2014 12:29:40 +0200 Subject: [PATCH 08/14] Fixing the blackoilfluid test. --- tests/test_blackoilfluid.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/test_blackoilfluid.cpp b/tests/test_blackoilfluid.cpp index 225e069d..3abf5293 100644 --- a/tests/test_blackoilfluid.cpp +++ b/tests/test_blackoilfluid.cpp @@ -251,7 +251,8 @@ BOOST_AUTO_TEST_CASE(test_liveoil) std::vector pvtRegionIdx(n, 0); // the tolerance for acceptable difference in values - const double reltol = 1e-9; + // const double reltol = 1e-9; + const double reltol = 1e-1; std::vector p(n); std::vector r(n); From ad54f1aa6208e2e24bdd074b2ddbb4449d986eec Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Mon, 22 Sep 2014 13:19:35 +0200 Subject: [PATCH 09/14] Fixing test_wetgas in test_blackoilfluid. --- tests/test_blackoilfluid.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/test_blackoilfluid.cpp b/tests/test_blackoilfluid.cpp index f495bff8..eb7c4e11 100644 --- a/tests/test_blackoilfluid.cpp +++ b/tests/test_blackoilfluid.cpp @@ -326,7 +326,8 @@ BOOST_AUTO_TEST_CASE(test_wetgas) std::vector pvtRegionIdx(n, 0); // the tolerance for acceptable difference in values - const double reltol = 1e-9; + // const double reltol = 1e-9; + const double reltol = 1e-1; std::vector p(n); std::vector r(n); From 09404202dafb8ca8faa7db915239132dd08218bb Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Mon, 22 Sep 2014 14:21:33 +0200 Subject: [PATCH 10/14] Changing the naming in PvtDead. To be consistent with the naming int PvtLive. --- opm/core/props/pvt/PvtDead.cpp | 48 +++++++++++++++++----------------- opm/core/props/pvt/PvtDead.hpp | 2 +- 2 files changed, 25 insertions(+), 25 deletions(-) diff --git a/opm/core/props/pvt/PvtDead.cpp b/opm/core/props/pvt/PvtDead.cpp index b7035ef8..fdb57a07 100644 --- a/opm/core/props/pvt/PvtDead.cpp +++ b/opm/core/props/pvt/PvtDead.cpp @@ -41,7 +41,7 @@ namespace Opm // resize the attributes of the object b_.resize(numRegions); viscosity_.resize(numRegions); - inverseBV_.resize(numRegions); + inverseBmu_.resize(numRegions); for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { const Opm::PvdoTable& pvdoTable = pvdoTables[regionIdx]; @@ -52,22 +52,22 @@ namespace Opm const std::vector& visc = pvdoTable.getViscosityColumn(); const int sz = b.size(); - std::vector bInv(sz); + std::vector inverseB(sz); for (int i = 0; i < sz; ++i) { - bInv[i] = 1.0 / b[i]; + inverseB[i] = 1.0 / b[i]; } // TODO: should we change the name of b so that we know it is the // inverse more explicitly? // or use the captial B in a exlicit way to show the difference. - std::vector bvInv(sz); + std::vector inverseBmu(sz); for (int i = 0; i < sz; ++i) { - bvInv[i] = 1.0 / (b[i] * visc[i]); + inverseBmu[i] = 1.0 / (b[i] * visc[i]); } - b_[regionIdx] = NonuniformTableLinear(press, bInv); + b_[regionIdx] = NonuniformTableLinear(press, inverseB); viscosity_[regionIdx] = NonuniformTableLinear(press, visc); - inverseBV_[regionIdx] = NonuniformTableLinear(press, bvInv); + inverseBmu_[regionIdx] = NonuniformTableLinear(press, inverseBmu); } } @@ -79,7 +79,7 @@ namespace Opm // resize the attributes of the object b_.resize(numRegions); viscosity_.resize(numRegions); - inverseBV_.resize(numRegions); + inverseBmu_.resize(numRegions); for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { const Opm::PvdgTable& pvdgTable = pvdgTables[regionIdx]; @@ -90,19 +90,19 @@ namespace Opm const std::vector& visc = pvdgTable.getViscosityColumn(); const int sz = b.size(); - std::vector bInv(sz); + std::vector inverseB(sz); for (int i = 0; i < sz; ++i) { - bInv[i] = 1.0 / b[i]; + inverseB[i] = 1.0 / b[i]; } - std::vector bvInv(sz); + std::vector inverseBmu(sz); for (int i = 0; i < sz; ++i) { - bvInv[i] = 1.0 / (b[i] * visc[i]); + inverseBmu[i] = 1.0 / (b[i] * visc[i]); } - b_[regionIdx] = NonuniformTableLinear(press, bInv); + b_[regionIdx] = NonuniformTableLinear(press, inverseB); viscosity_[regionIdx] = NonuniformTableLinear(press, visc); - inverseBV_[regionIdx] = NonuniformTableLinear(press, bvInv); + inverseBmu_[regionIdx] = NonuniformTableLinear(press, inverseBmu); } } @@ -124,8 +124,8 @@ namespace Opm int regionIdx = getTableIndex_(pvtTableIdx, i); // output_mu[i] = viscosity_[regionIdx](p[i]); double tempInvB = b_[regionIdx](p[i]); - double tempInvBV = inverseBV_[regionIdx](p[i]); - output_mu[i] = tempInvB / tempInvBV; + double tempInvBmu = inverseBmu_[regionIdx](p[i]); + output_mu[i] = tempInvB / tempInvBmu; } } @@ -143,10 +143,10 @@ namespace Opm // output_mu[i] = viscosity_[regionIdx](p[i]); // output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]); double tempInvB = b_[regionIdx](p[i]); - double tempInvBV = inverseBV_[regionIdx](p[i]); - output_mu[i] = tempInvB / tempInvBV; - output_dmudp[i] = (tempInvBV * b_[regionIdx].derivative(p[i]) - - tempInvB * inverseBV_[regionIdx].derivative(p[i])) / (tempInvBV * tempInvBV); + double tempInvBmu = inverseBmu_[regionIdx](p[i]); + output_mu[i] = tempInvB / tempInvBmu; + output_dmudp[i] = (tempInvBmu * b_[regionIdx].derivative(p[i]) + - tempInvB * inverseBmu_[regionIdx].derivative(p[i])) / (tempInvBmu * tempInvBmu); } std::fill(output_dmudr, output_dmudr + n, 0.0); @@ -165,13 +165,13 @@ namespace Opm for (int i = 0; i < n; ++i) { int regionIdx = getTableIndex_(pvtTableIdx, i); double tempInvB = b_[regionIdx](p[i]); - double tempInvBV = inverseBV_[regionIdx](p[i]); + double tempInvBmu = inverseBmu_[regionIdx](p[i]); // output_mu[i] = viscosity_[regionIdx](p[i]); - output_mu[i] = tempInvB / tempInvBV; + output_mu[i] = tempInvB / tempInvBmu; // output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]); // output_dmudp[i] = tempInvB / (tempInvBV * tempInvBV) * inverseBV_[regionIdx].derivative(p[i]); - output_dmudp[i] = (tempInvBV * b_[regionIdx].derivative(p[i]) - - tempInvB * inverseBV_[regionIdx].derivative(p[i])) / (tempInvBV * tempInvBV); + output_dmudp[i] = (tempInvBmu * b_[regionIdx].derivative(p[i]) + - tempInvB * inverseBmu_[regionIdx].derivative(p[i])) / (tempInvBmu * tempInvBmu); } std::fill(output_dmudr, output_dmudr + n, 0.0); diff --git a/opm/core/props/pvt/PvtDead.hpp b/opm/core/props/pvt/PvtDead.hpp index 2a663e2d..ccc9acac 100644 --- a/opm/core/props/pvt/PvtDead.hpp +++ b/opm/core/props/pvt/PvtDead.hpp @@ -155,7 +155,7 @@ namespace Opm // table per PVT region. std::vector > b_; std::vector > viscosity_; - std::vector > inverseBV_; + std::vector > inverseBmu_; }; } From 6a43935ea3910ef23ae50c76027fd9fc3f699274 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Mon, 22 Sep 2014 14:26:00 +0200 Subject: [PATCH 11/14] Deleting some commented lines. --- opm/core/props/pvt/PvtDead.cpp | 12 ++---------- opm/core/props/pvt/PvtLiveGas.cpp | 6 ------ opm/core/props/pvt/PvtLiveOil.cpp | 10 ---------- 3 files changed, 2 insertions(+), 26 deletions(-) diff --git a/opm/core/props/pvt/PvtDead.cpp b/opm/core/props/pvt/PvtDead.cpp index fdb57a07..db9cdb39 100644 --- a/opm/core/props/pvt/PvtDead.cpp +++ b/opm/core/props/pvt/PvtDead.cpp @@ -57,9 +57,6 @@ namespace Opm inverseB[i] = 1.0 / b[i]; } - // TODO: should we change the name of b so that we know it is the - // inverse more explicitly? - // or use the captial B in a exlicit way to show the difference. std::vector inverseBmu(sz); for (int i = 0; i < sz; ++i) { inverseBmu[i] = 1.0 / (b[i] * visc[i]); @@ -122,7 +119,6 @@ namespace Opm // #pragma omp parallel for for (int i = 0; i < n; ++i) { int regionIdx = getTableIndex_(pvtTableIdx, i); - // output_mu[i] = viscosity_[regionIdx](p[i]); double tempInvB = b_[regionIdx](p[i]); double tempInvBmu = inverseBmu_[regionIdx](p[i]); output_mu[i] = tempInvB / tempInvBmu; @@ -140,8 +136,6 @@ namespace Opm // #pragma omp parallel for for (int i = 0; i < n; ++i) { int regionIdx = getTableIndex_(pvtTableIdx, i); - // output_mu[i] = viscosity_[regionIdx](p[i]); - // output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]); double tempInvB = b_[regionIdx](p[i]); double tempInvBmu = inverseBmu_[regionIdx](p[i]); output_mu[i] = tempInvB / tempInvBmu; @@ -166,12 +160,10 @@ namespace Opm int regionIdx = getTableIndex_(pvtTableIdx, i); double tempInvB = b_[regionIdx](p[i]); double tempInvBmu = inverseBmu_[regionIdx](p[i]); - // output_mu[i] = viscosity_[regionIdx](p[i]); output_mu[i] = tempInvB / tempInvBmu; - // output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]); - // output_dmudp[i] = tempInvB / (tempInvBV * tempInvBV) * inverseBV_[regionIdx].derivative(p[i]); output_dmudp[i] = (tempInvBmu * b_[regionIdx].derivative(p[i]) - - tempInvB * inverseBmu_[regionIdx].derivative(p[i])) / (tempInvBmu * tempInvBmu); + - tempInvB * inverseBmu_[regionIdx].derivative(p[i])) + / (tempInvBmu * tempInvBmu); } std::fill(output_dmudr, output_dmudr + n, 0.0); diff --git a/opm/core/props/pvt/PvtLiveGas.cpp b/opm/core/props/pvt/PvtLiveGas.cpp index 76c962f8..fb8e5224 100644 --- a/opm/core/props/pvt/PvtLiveGas.cpp +++ b/opm/core/props/pvt/PvtLiveGas.cpp @@ -65,7 +65,6 @@ namespace Opm // The number of the columns int nColumns = saturated_gas_table_[pvtTableIdx][2].size(); saturated_gas_table_[pvtTableIdx][3].resize(nColumns); // allocate memory for 1/(Bg*mu_g) - // saturated_gas_table_[pvtTableIdx][3] = pvtgTable.getOuterTable()->getGasFormationFactorColumn(); // TODO saturated_gas_table_[pvtTableIdx][4] = pvtgTable.getOuterTable()->getOilSolubilityColumn(); // Rv int sz = pvtgTable.getOuterTable()->numRows(); @@ -80,7 +79,6 @@ namespace Opm undersat_gas_tables_[pvtTableIdx][i][2] = undersatTable.getGasViscosityColumn(); // mu_g int nColumns = undersat_gas_tables_[pvtTableIdx][i][2].size(); undersat_gas_tables_[pvtTableIdx][i][3].resize(nColumns); // allocate memory for 1/(Bg*mu_g) - // undersat_gas_tables_[pvtTableIdx][i][3] = undersatTable.getGasViscosityColumn(); // TODO } // Bg -> 1/Bg @@ -115,7 +113,6 @@ namespace Opm double inverseBMu = miscible_gas(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i), 3, false); output_mu[i] = inverseB / inverseBMu; - // output_mu[i] = miscible_gas(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i), 2, false); } } @@ -162,9 +159,6 @@ namespace Opm output_dmudr[i] = (inverseBMu * dinverseBdr - inverseB * dinverseBmudr) / (inverseBMu * inverseBMu); - // output_mu[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 2, 0); - // output_dmudp[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 2, 1); - // output_dmudr[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 2, 2); } } diff --git a/opm/core/props/pvt/PvtLiveOil.cpp b/opm/core/props/pvt/PvtLiveOil.cpp index 7d0421d8..a9be4e7b 100644 --- a/opm/core/props/pvt/PvtLiveOil.cpp +++ b/opm/core/props/pvt/PvtLiveOil.cpp @@ -147,9 +147,6 @@ namespace Opm double inverseBMu = miscible_oil(p[i], z + num_phases_*i, tableIdx, 3, false); output_mu[i] = inverseB / inverseBMu; - - // output_mu[i] = miscible_oil(p[i], z + num_phases_*i, tableIdx, 2, false); - // output_mu[i] = miscible_oil(p[i], z + num_phases_*i, tableIdx, 2, false); } } @@ -184,10 +181,6 @@ namespace Opm output_dmudr[i] = (inverseBMu * dinverseBdr - inverseB * dinverseBmudr) / (inverseBMu * inverseBMu); - // output_mu[i] = miscible_oil(p[i], r[i], tableIdx, 2, 0); - // output_dmudp[i] = miscible_oil(p[i], r[i], tableIdx, 2, 1); - // output_dmudr[i] = miscible_oil(p[i], r[i], tableIdx, 2, 2); - } } @@ -223,9 +216,6 @@ namespace Opm output_dmudr[i] = (inverseBMu * dinverseBdr - inverseB * dinverseBmudr) / (inverseBMu * inverseBMu); - // output_mu[i] = miscible_oil(p[i], r[i], cnd, tableIdx, 2, 0); - // output_dmudp[i] = miscible_oil(p[i], r[i], cnd, tableIdx, 2, 1); - // output_dmudr[i] = miscible_oil(p[i], r[i], cnd, tableIdx, 2, 2); } } From 175da5450fbf10975aa7f7f1dada2b8bf0baac87 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Mon, 22 Sep 2014 14:29:18 +0200 Subject: [PATCH 12/14] Changing the location of * in long formula. --- opm/core/props/pvt/PvtLiveGas.cpp | 8 ++++---- opm/core/props/pvt/PvtLiveOil.cpp | 4 ++-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/opm/core/props/pvt/PvtLiveGas.cpp b/opm/core/props/pvt/PvtLiveGas.cpp index fb8e5224..91dbf812 100644 --- a/opm/core/props/pvt/PvtLiveGas.cpp +++ b/opm/core/props/pvt/PvtLiveGas.cpp @@ -83,12 +83,12 @@ namespace Opm // Bg -> 1/Bg for (int i=0; igetPressureColumn()[i]; // p saturated_oil_table_[pvtTableIdx][1][i] = 1.0/saturatedPvto->getOilFormationFactorColumn()[i]; // 1/Bo saturated_oil_table_[pvtTableIdx][2][i] = saturatedPvto->getOilViscosityColumn()[i]; // mu_o - saturated_oil_table_[pvtTableIdx][3][i] = 1.0 / (saturatedPvto->getOilFormationFactorColumn()[i] * - saturatedPvto->getOilViscosityColumn()[i]); // 1/(Bo*mu_o) + saturated_oil_table_[pvtTableIdx][3][i] = 1.0 / (saturatedPvto->getOilFormationFactorColumn()[i] + * saturatedPvto->getOilViscosityColumn()[i]); // 1/(Bo*mu_o) saturated_oil_table_[pvtTableIdx][4][i] = saturatedPvto->getGasSolubilityColumn()[i]; // Rs } From ef3e524c309035c766aef2d984747d1643f02cb8 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Tue, 23 Sep 2014 13:17:47 +0200 Subject: [PATCH 13/14] Changing nColumns to nRows. --- opm/core/props/pvt/PvtLiveGas.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/opm/core/props/pvt/PvtLiveGas.cpp b/opm/core/props/pvt/PvtLiveGas.cpp index 91dbf812..63446b3c 100644 --- a/opm/core/props/pvt/PvtLiveGas.cpp +++ b/opm/core/props/pvt/PvtLiveGas.cpp @@ -63,8 +63,8 @@ namespace Opm saturated_gas_table_[pvtTableIdx][1] = pvtgTable.getOuterTable()->getGasFormationFactorColumn(); // Bg saturated_gas_table_[pvtTableIdx][2] = pvtgTable.getOuterTable()->getGasViscosityColumn(); // mu_g // The number of the columns - int nColumns = saturated_gas_table_[pvtTableIdx][2].size(); - saturated_gas_table_[pvtTableIdx][3].resize(nColumns); // allocate memory for 1/(Bg*mu_g) + int nRows = saturated_gas_table_[pvtTableIdx][2].size(); + saturated_gas_table_[pvtTableIdx][3].resize(nRows); // allocate memory for 1/(Bg*mu_g) saturated_gas_table_[pvtTableIdx][4] = pvtgTable.getOuterTable()->getOilSolubilityColumn(); // Rv int sz = pvtgTable.getOuterTable()->numRows(); @@ -77,8 +77,8 @@ namespace Opm undersat_gas_tables_[pvtTableIdx][i][0] = undersatTable.getOilSolubilityColumn(); // Rv undersat_gas_tables_[pvtTableIdx][i][1] = undersatTable.getGasFormationFactorColumn(); // Bg undersat_gas_tables_[pvtTableIdx][i][2] = undersatTable.getGasViscosityColumn(); // mu_g - int nColumns = undersat_gas_tables_[pvtTableIdx][i][2].size(); - undersat_gas_tables_[pvtTableIdx][i][3].resize(nColumns); // allocate memory for 1/(Bg*mu_g) + int nRows = undersat_gas_tables_[pvtTableIdx][i][2].size(); + undersat_gas_tables_[pvtTableIdx][i][3].resize(nRows); // allocate memory for 1/(Bg*mu_g) } // Bg -> 1/Bg From ee265f5cb515d7038056bec060de38b7109fb31b Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Tue, 23 Sep 2014 14:38:32 +0200 Subject: [PATCH 14/14] Changing the tolerance for blackoilfluid tests. Using the original strict relative tolerance for other properties than viscosity. --- tests/test_blackoilfluid.cpp | 30 ++++++++++++++++-------------- 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/tests/test_blackoilfluid.cpp b/tests/test_blackoilfluid.cpp index eb7c4e11..add7d495 100644 --- a/tests/test_blackoilfluid.cpp +++ b/tests/test_blackoilfluid.cpp @@ -247,9 +247,10 @@ BOOST_AUTO_TEST_CASE(test_liveoil) const int np = phase_usage_.num_phases; std::vector pvtRegionIdx(n, 0); - // the tolerance for acceptable difference in values - // const double reltol = 1e-9; - const double reltol = 1e-1; + // the relative tolerance in percentage for acceptable difference in values + const double reltolper = 1e-9; + // the relative tolerance in percentage for acceptable difference in values for viscosity + const double reltolpermu = 1e-1; std::vector p(n); std::vector r(n); @@ -291,13 +292,13 @@ BOOST_AUTO_TEST_CASE(test_liveoil) } - testmu(reltol, n, np, pvtRegionIdx, p, r,z, props_, condition); + testmu(reltolpermu, n, np, pvtRegionIdx, p, r,z, props_, condition); - testb(reltol,n,np,pvtRegionIdx,p,r,z,props_,condition); + testb(reltolper,n,np,pvtRegionIdx,p,r,z,props_,condition); - testrsSat(reltol,n,np,pvtRegionIdx,p,props_); + testrsSat(reltolper,n,np,pvtRegionIdx,p,props_); - testrvSat(reltol,n,np,pvtRegionIdx,p,props_); + testrvSat(reltolper,n,np,pvtRegionIdx,p,props_); } @@ -325,9 +326,10 @@ BOOST_AUTO_TEST_CASE(test_wetgas) const int np = phase_usage_.num_phases; std::vector pvtRegionIdx(n, 0); - // the tolerance for acceptable difference in values - // const double reltol = 1e-9; - const double reltol = 1e-1; + // the relative tolerance in percentage for acceptable difference in values + const double reltolper = 1e-9; + // the relative tolerance in percentage for acceptable difference in values for viscosity + const double reltolpermu = 1e-1; std::vector p(n); std::vector r(n); @@ -369,12 +371,12 @@ BOOST_AUTO_TEST_CASE(test_wetgas) } - testmu(reltol, n, np, pvtRegionIdx, p, r,z, props_, condition); + testmu(reltolpermu, n, np, pvtRegionIdx, p, r,z, props_, condition); - testb(reltol,n,np,pvtRegionIdx,p,r,z,props_,condition); + testb(reltolper,n,np,pvtRegionIdx,p,r,z,props_,condition); - testrsSat(reltol,n,np,pvtRegionIdx,p,props_); + testrsSat(reltolper,n,np,pvtRegionIdx,p,props_); - testrvSat(reltol,n,np,pvtRegionIdx,p,props_); + testrvSat(reltolper,n,np,pvtRegionIdx,p,props_); }