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

Make the flux actually a velocity, also for the new callpath.

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen 2022-07-07 00:30:05 +02:00
parent 6c0d5ea8c5
commit b6986a24f7
2 changed files with 53 additions and 35 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
opm/models/blackoil/blackoillocalresidualtpfa.hh
View File

@ -176,26 +176,28 @@ public:
} }
/*! /*!
* \copydoc FvBaseLocalResidual::computeFlux * This function works like the ElementContext-based version with [two] one
* main differences:
* - The darcy flux is calculated here, not read from the extensive quantities of the element context.
* [- The flux is not per area (a velocity), i.e. it is in [m^3/s], not [m/s].]
*/ */
static void computeFlux(RateVector& flux, static void computeFlux(RateVector& flux,
const Problem& problem, const Problem& problem,
const unsigned globalIndexIn, const unsigned globalIndexIn,
const unsigned globalIndexEx, const unsigned globalIndexEx,
const IntensiveQuantities& intQuantsIn, const IntensiveQuantities& intQuantsIn,
const IntensiveQuantities& intQuantsEx, const IntensiveQuantities& intQuantsEx,
const unsigned timeIdx) const unsigned timeIdx,
const Scalar trans,
const Scalar faceArea)
{ {
assert(timeIdx == 0); assert(timeIdx == 0);
flux = 0.0; flux = 0.0;
Scalar Vin = problem.model().dofTotalVolume(globalIndexIn); Scalar Vin = problem.model().dofTotalVolume(globalIndexIn);
Scalar Vex = problem.model().dofTotalVolume(globalIndexEx); Scalar Vex = problem.model().dofTotalVolume(globalIndexEx);
//Scalar faceArea = 1.0; // This makes the 'flux' a volume flux not a velocity.
Scalar trans = 1.0; // problem.transmissibility(globalIndexIn,globalIndexEx);
// Scalar faceArea = problem.area(globalIndexIn,globalIndexEx);
Scalar faceArea = 1.0; // NB need correct calculation local residual
Scalar thpres = problem.thresholdPressure(globalIndexIn, globalIndexEx); Scalar thpres = problem.thresholdPressure(globalIndexIn, globalIndexEx);
// estimate the gravity correction: for performance reasons we use a simplified // estimate the gravity correction: for performance reasons we use a simplified
@ -215,8 +217,6 @@ public:
// exterior DOF) // exterior DOF)
Scalar distZ = zIn - zEx; // NB could be precalculated Scalar distZ = zIn - zEx; // NB could be precalculated
//
// const ExtensiveQuantities& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
calculateFluxes_(flux, calculateFluxes_(flux,
problem, // should be removed problem, // should be removed
intQuantsIn, intQuantsIn,

72
opm/models/discretization/common/tpfalinearizer.hh
View File

@ -311,6 +311,10 @@ private:
using NeighborSet = std::set< unsigned >; using NeighborSet = std::set< unsigned >;
std::vector<NeighborSet> sparsityPattern(model.numTotalDof()); std::vector<NeighborSet> sparsityPattern(model.numTotalDof());
unsigned numCells = model.numTotalDof();
neighborInfo_.reserve(numCells, 6 * numCells);
std::vector<NeighborInfo> loc_nbinfo;
ElementIterator elemIt = gridView_().template begin<0>(); ElementIterator elemIt = gridView_().template begin<0>();
const ElementIterator elemEndIt = gridView_().template end<0>(); const ElementIterator elemEndIt = gridView_().template end<0>();
for (; elemIt != elemEndIt; ++elemIt) { for (; elemIt != elemEndIt; ++elemIt) {
@ -319,16 +323,23 @@ private:
for (unsigned primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++primaryDofIdx) { for (unsigned primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++primaryDofIdx) {
unsigned myIdx = stencil.globalSpaceIndex(primaryDofIdx); unsigned myIdx = stencil.globalSpaceIndex(primaryDofIdx);
loc_nbinfo.resize(stencil.numDof() - 1); // Do not include the primary dof.
for (unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) { for (unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx); unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
sparsityPattern[myIdx].insert(neighborIdx); sparsityPattern[myIdx].insert(neighborIdx);
if (dofIdx > 0) {
const double trans = problem_().transmissibility(myIdx, neighborIdx);
const double area = stencil.interiorFace(dofIdx - 1).area();
loc_nbinfo[dofIdx - 1] = NeighborInfo{neighborIdx, trans, area};
}
} }
neighborInfo_.appendRow(loc_nbinfo.begin(), loc_nbinfo.end());
} }
} }
// Do not include auxiliary connections in the linearization sparsity pattern. // Do not include auxiliary connections in the linearization sparsity pattern.
auto reservoirSparsityPattern = sparsityPattern; // auto reservoirSparsityPattern = sparsityPattern;
// add the additional neighbors and degrees of freedom caused by the auxiliary // add the additional neighbors and degrees of freedom caused by the auxiliary
// equations // equations
@ -344,25 +355,26 @@ private:
// Now generate the neighbours_ and trans_ structures for the linearization loop. // Now generate the neighbours_ and trans_ structures for the linearization loop.
// Those should not include an entry for the cell itself. // Those should not include an entry for the cell itself.
for (unsigned globI = 0; globI < model.numTotalDof(); globI++) { // for (unsigned globI = 0; globI < model.numTotalDof(); globI++) {
reservoirSparsityPattern[globI].erase(globI); // reservoirSparsityPattern[globI].erase(globI);
} // }
unsigned numCells = model.numTotalDof(); // unsigned numCells = model.numTotalDof();
neighbours_.reserve(numCells, 6 * numCells); // neighbours_.reserve(numCells, 6 * numCells);
trans_.reserve(numCells, 6 * numCells); // trans_.reserve(numCells, 6 * numCells);
std::vector<double> loctrans; // faceArea_.reserve(numCells, 6 * numCells);
for (unsigned globI = 0; globI < numCells; globI++) { // std::vector<double> loctrans;
const auto& cells = reservoirSparsityPattern[globI]; // for (unsigned globI = 0; globI < numCells; globI++) {
neighbours_.appendRow(cells.begin(), cells.end()); // const auto& cells = reservoirSparsityPattern[globI];
unsigned n = cells.size(); // neighbours_.appendRow(cells.begin(), cells.end());
loctrans.resize(n); // unsigned n = cells.size();
short loc = 0; // loctrans.resize(n);
for (const int& cell : cells) { // short loc = 0;
loctrans[loc] = problem_().transmissibility(globI, cell); // for (const int& cell : cells) {
loc++; // loctrans[loc] = problem_().transmissibility(globI, cell);
} // loc++;
trans_.appendRow(loctrans.begin(), loctrans.end()); // }
} // trans_.appendRow(loctrans.begin(), loctrans.end());
// }
} }
// reset the global linear system of equations. // reset the global linear system of equations.
@ -400,7 +412,7 @@ private:
#pragma omp parallel for #pragma omp parallel for
#endif #endif
for (unsigned globI = 0; globI < numCells; globI++) { for (unsigned globI = 0; globI < numCells; globI++) {
const auto& neighbours = neighbours_[globI]; // this is a set but should maybe be changed const auto& nbInfos = neighborInfo_[globI]; // this is a set but should maybe be changed
// accumulation term // accumulation term
double dt = simulator_().timeStepSize(); double dt = simulator_().timeStepSize();
double volume = model_().dofTotalVolume(globI); double volume = model_().dofTotalVolume(globI);
@ -435,7 +447,8 @@ private:
jacobian_->addToBlock(globI, globI, bMat); jacobian_->addToBlock(globI, globI, bMat);
} }
short loc = 0; short loc = 0;
for (const auto& globJ : neighbours) { for (const auto& nbInfo : nbInfos) {
unsigned globJ = nbInfo.neighbor;
assert(globJ != globI); assert(globJ != globI);
res = 0.0; res = 0.0;
bMat = 0.0; bMat = 0.0;
@ -444,14 +457,14 @@ private:
assert(intQuantsExP); assert(intQuantsExP);
const IntensiveQuantities& intQuantsEx = *intQuantsExP; const IntensiveQuantities& intQuantsEx = *intQuantsExP;
LocalResidual::computeFlux( LocalResidual::computeFlux(
adres, problem_(), globI, globJ, intQuantsIn, intQuantsEx, 0); adres, problem_(), globI, globJ, intQuantsIn, intQuantsEx, 0, nbInfo.trans, nbInfo.faceArea);
adres *= trans_[globI][loc]; adres *= nbInfo.faceArea;
setResAndJacobi(res, bMat, adres); setResAndJacobi(res, bMat, adres);
residual_[globI] += res; residual_[globI] += res;
jacobian_->addToBlock(globI, globI, bMat); jacobian_->addToBlock(globI, globI, bMat);
bMat *= -1.0; bMat *= -1.0;
jacobian_->addToBlock(globJ, globI, bMat); jacobian_->addToBlock(globJ, globI, bMat);
loc++; ++loc;
} }
} }
if (not(well_local)) { if (not(well_local)) {
@ -472,8 +485,13 @@ private:
LinearizationType linearizationType_; LinearizationType linearizationType_;
SparseTable<unsigned> neighbours_; struct NeighborInfo
SparseTable<double> trans_; {
unsigned int neighbor;
double trans;
double faceArea;
};
SparseTable<NeighborInfo> neighborInfo_;
}; };
} // namespace Opm } // namespace Opm