OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #1206 from GitPaean/switching_well_primary_variables

Browse Source 
Not switching the order of well primary variables and well equations
This commit is contained in:
Atgeirr Flø Rasmussen 2017-06-22 19:37:40 +02:00 committed by GitHub
commit 2ceca48f04
2 changed files with 17 additions and 34 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
opm/autodiff/StandardWellsDense.hpp
View File

@ -174,8 +174,6 @@ enum WellVariablePositions {
int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const;
int ebosCompToFlowPhaseIdx( const int compIdx ) const;
std::vector<double>
extractPerfData(const std::vector<double>& in) const;

49
opm/autodiff/StandardWellsDense_impl.hpp
View File

@ -209,23 +209,23 @@ namespace Opm {
}
// subtract sum of phase fluxes in the well equations.
resWell_[w][flowPhaseToEbosCompIdx(componentIdx)] -= cq_s[componentIdx].value();
resWell_[w][componentIdx] -= cq_s[componentIdx].value();
// assemble the jacobians
for (int pvIdx = 0; pvIdx < numWellEq; ++pvIdx) {
if (!only_wells) {
// also need to consider the efficiency factor when manipulating the jacobians.
ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(componentIdx)][flowToEbosPvIdx(pvIdx)] -= cq_s_effective.derivative(pvIdx);
duneB_[w][cell_idx][flowToEbosPvIdx(pvIdx)][flowPhaseToEbosCompIdx(componentIdx)] -= cq_s_effective.derivative(pvIdx+numEq); // intput in transformed matrix
duneC_[w][cell_idx][flowPhaseToEbosCompIdx(componentIdx)][flowToEbosPvIdx(pvIdx)] -= cq_s_effective.derivative(pvIdx);
duneB_[w][cell_idx][pvIdx][flowPhaseToEbosCompIdx(componentIdx)] -= cq_s_effective.derivative(pvIdx+numEq); // intput in transformed matrix
duneC_[w][cell_idx][componentIdx][flowToEbosPvIdx(pvIdx)] -= cq_s_effective.derivative(pvIdx);
}
invDuneD_[w][w][flowPhaseToEbosCompIdx(componentIdx)][flowToEbosPvIdx(pvIdx)] -= cq_s[componentIdx].derivative(pvIdx+numEq);
invDuneD_[w][w][componentIdx][pvIdx] -= cq_s[componentIdx].derivative(pvIdx+numEq);
}
// add trivial equation for 2p cases (Only support water + oil)
if (numComp == 2) {
assert(!active_[ Gas ]);
invDuneD_[w][w][flowPhaseToEbosCompIdx(Gas)][flowToEbosPvIdx(Gas)] = 1.0;
invDuneD_[w][w][Gas][Gas] = 1.0;
}
// Store the perforation phase flux for later usage.
@ -245,9 +245,9 @@ namespace Opm {
EvalWell resWell_loc = (wellSurfaceVolumeFraction(w, componentIdx) - F0_[w + nw*componentIdx]) * volume / dt;
resWell_loc += getQs(w, componentIdx);
for (int pvIdx = 0; pvIdx < numWellEq; ++pvIdx) {
invDuneD_[w][w][flowPhaseToEbosCompIdx(componentIdx)][flowToEbosPvIdx(pvIdx)] += resWell_loc.derivative(pvIdx+numEq);
invDuneD_[w][w][componentIdx][pvIdx] += resWell_loc.derivative(pvIdx+numEq);
}
resWell_[w][flowPhaseToEbosCompIdx(componentIdx)] += resWell_loc.value();
resWell_[w][componentIdx] += resWell_loc.value();
}
}
@ -501,20 +501,6 @@ namespace Opm {
template<typename TypeTag>
int
StandardWellsDense<TypeTag>::
ebosCompToFlowPhaseIdx( const int compIdx ) const
{
assert(compIdx < 3);
const int compToPhase[ 3 ] = { Oil, Water, Gas };
return compToPhase[ compIdx ];
}
template<typename TypeTag>
std::vector<double>
StandardWellsDense<TypeTag>::
@ -968,10 +954,9 @@ namespace Opm {
const int numComp = numComponents();
std::vector<double> res(numComp*nw);
for( int compIdx = 0; compIdx < numComp; ++compIdx) {
const int ebosCompIdx = flowPhaseToEbosCompIdx(compIdx);
for (int wellIdx = 0; wellIdx < nw; ++wellIdx) {
int idx = wellIdx + nw*compIdx;
res[idx] = resWell_[ wellIdx ][ ebosCompIdx ];
res[idx] = resWell_[ wellIdx ][ compIdx ];
}
}
return res;
@ -1250,20 +1235,20 @@ namespace Opm {
// update the second and third well variable (The flux fractions)
std::vector<double> F(np,0.0);
if (active_[ Water ]) {
const int sign2 = dwells[w][flowPhaseToEbosCompIdx(WFrac)] > 0 ? 1: -1;
const double dx2_limited = sign2 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(WFrac)]),dFLimit);
const int sign2 = dwells[w][WFrac] > 0 ? 1: -1;
const double dx2_limited = sign2 * std::min(std::abs(dwells[w][WFrac]),dFLimit);
well_state.wellSolutions()[WFrac*nw + w] = xvar_well_old[WFrac*nw + w] - dx2_limited;
}
if (active_[ Gas ]) {
const int sign3 = dwells[w][flowPhaseToEbosCompIdx(GFrac)] > 0 ? 1: -1;
const double dx3_limited = sign3 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(GFrac)]),dFLimit);
const int sign3 = dwells[w][GFrac] > 0 ? 1: -1;
const double dx3_limited = sign3 * std::min(std::abs(dwells[w][GFrac]),dFLimit);
well_state.wellSolutions()[GFrac*nw + w] = xvar_well_old[GFrac*nw + w] - dx3_limited;
}
if (has_solvent_) {
const int sign4 = dwells[w][flowPhaseToEbosCompIdx(SFrac)] > 0 ? 1: -1;
const double dx4_limited = sign4 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(SFrac)]),dFLimit);
const int sign4 = dwells[w][SFrac] > 0 ? 1: -1;
const double dx4_limited = sign4 * std::min(std::abs(dwells[w][SFrac]),dFLimit);
well_state.wellSolutions()[SFrac*nw + w] = xvar_well_old[SFrac*nw + w] - dx4_limited;
}
@ -1367,7 +1352,7 @@ namespace Opm {
case THP: // The BHP and THP both uses the total rate as first well variable.
case BHP:
{
well_state.wellSolutions()[nw*XvarWell + w] = xvar_well_old[nw*XvarWell + w] - dwells[w][flowPhaseToEbosCompIdx(XvarWell)];
well_state.wellSolutions()[nw*XvarWell + w] = xvar_well_old[nw*XvarWell + w] - dwells[w][XvarWell];
switch (wells().type[w]) {
case INJECTOR:
@ -1435,8 +1420,8 @@ namespace Opm {
case SURFACE_RATE: // Both rate controls use bhp as first well variable
case RESERVOIR_RATE:
{
const int sign1 = dwells[w][flowPhaseToEbosCompIdx(XvarWell)] > 0 ? 1: -1;
const double dx1_limited = sign1 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(XvarWell)]),std::abs(xvar_well_old[nw*XvarWell + w])*dBHPLimit);
const int sign1 = dwells[w][XvarWell] > 0 ? 1: -1;
const double dx1_limited = sign1 * std::min(std::abs(dwells[w][XvarWell]),std::abs(xvar_well_old[nw*XvarWell + w])*dBHPLimit);
well_state.wellSolutions()[nw*XvarWell + w] = std::max(xvar_well_old[nw*XvarWell + w] - dx1_limited,1e5);
well_state.bhp()[w] = well_state.wellSolutions()[nw*XvarWell + w];