diff --git a/opm/models/blackoil/blackoillocalresidualtpfa.hh b/opm/models/blackoil/blackoillocalresidualtpfa.hh index 154843c45..3a8a4c83c 100644 --- a/opm/models/blackoil/blackoillocalresidualtpfa.hh +++ b/opm/models/blackoil/blackoillocalresidualtpfa.hh @@ -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, const Problem& problem, const unsigned globalIndexIn, const unsigned globalIndexEx, const IntensiveQuantities& intQuantsIn, const IntensiveQuantities& intQuantsEx, - const unsigned timeIdx) + const unsigned timeIdx, + const Scalar trans, + const Scalar faceArea) { assert(timeIdx == 0); flux = 0.0; Scalar Vin = problem.model().dofTotalVolume(globalIndexIn); 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); // estimate the gravity correction: for performance reasons we use a simplified @@ -215,8 +217,6 @@ public: // exterior DOF) Scalar distZ = zIn - zEx; // NB could be precalculated - // - // const ExtensiveQuantities& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx); calculateFluxes_(flux, problem, // should be removed intQuantsIn, diff --git a/opm/models/discretization/common/tpfalinearizer.hh b/opm/models/discretization/common/tpfalinearizer.hh index 437721183..fa92b0f70 100644 --- a/opm/models/discretization/common/tpfalinearizer.hh +++ b/opm/models/discretization/common/tpfalinearizer.hh @@ -311,6 +311,10 @@ private: using NeighborSet = std::set< unsigned >; std::vector sparsityPattern(model.numTotalDof()); + unsigned numCells = model.numTotalDof(); + neighborInfo_.reserve(numCells, 6 * numCells); + std::vector loc_nbinfo; + ElementIterator elemIt = gridView_().template begin<0>(); const ElementIterator elemEndIt = gridView_().template end<0>(); for (; elemIt != elemEndIt; ++elemIt) { @@ -319,16 +323,23 @@ private: for (unsigned primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++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) { unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx); 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. - auto reservoirSparsityPattern = sparsityPattern; + // auto reservoirSparsityPattern = sparsityPattern; // add the additional neighbors and degrees of freedom caused by the auxiliary // equations @@ -344,25 +355,26 @@ private: // Now generate the neighbours_ and trans_ structures for the linearization loop. // Those should not include an entry for the cell itself. - for (unsigned globI = 0; globI < model.numTotalDof(); globI++) { - reservoirSparsityPattern[globI].erase(globI); - } - unsigned numCells = model.numTotalDof(); - neighbours_.reserve(numCells, 6 * numCells); - trans_.reserve(numCells, 6 * numCells); - std::vector loctrans; - for (unsigned globI = 0; globI < numCells; globI++) { - const auto& cells = reservoirSparsityPattern[globI]; - neighbours_.appendRow(cells.begin(), cells.end()); - unsigned n = cells.size(); - loctrans.resize(n); - short loc = 0; - for (const int& cell : cells) { - loctrans[loc] = problem_().transmissibility(globI, cell); - loc++; - } - trans_.appendRow(loctrans.begin(), loctrans.end()); - } + // for (unsigned globI = 0; globI < model.numTotalDof(); globI++) { + // reservoirSparsityPattern[globI].erase(globI); + // } + // unsigned numCells = model.numTotalDof(); + // neighbours_.reserve(numCells, 6 * numCells); + // trans_.reserve(numCells, 6 * numCells); + // faceArea_.reserve(numCells, 6 * numCells); + // std::vector loctrans; + // for (unsigned globI = 0; globI < numCells; globI++) { + // const auto& cells = reservoirSparsityPattern[globI]; + // neighbours_.appendRow(cells.begin(), cells.end()); + // unsigned n = cells.size(); + // loctrans.resize(n); + // short loc = 0; + // for (const int& cell : cells) { + // loctrans[loc] = problem_().transmissibility(globI, cell); + // loc++; + // } + // trans_.appendRow(loctrans.begin(), loctrans.end()); + // } } // reset the global linear system of equations. @@ -400,7 +412,7 @@ private: #pragma omp parallel for #endif 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 double dt = simulator_().timeStepSize(); double volume = model_().dofTotalVolume(globI); @@ -435,7 +447,8 @@ private: jacobian_->addToBlock(globI, globI, bMat); } short loc = 0; - for (const auto& globJ : neighbours) { + for (const auto& nbInfo : nbInfos) { + unsigned globJ = nbInfo.neighbor; assert(globJ != globI); res = 0.0; bMat = 0.0; @@ -444,14 +457,14 @@ private: assert(intQuantsExP); const IntensiveQuantities& intQuantsEx = *intQuantsExP; LocalResidual::computeFlux( - adres, problem_(), globI, globJ, intQuantsIn, intQuantsEx, 0); - adres *= trans_[globI][loc]; + adres, problem_(), globI, globJ, intQuantsIn, intQuantsEx, 0, nbInfo.trans, nbInfo.faceArea); + adres *= nbInfo.faceArea; setResAndJacobi(res, bMat, adres); residual_[globI] += res; jacobian_->addToBlock(globI, globI, bMat); bMat *= -1.0; jacobian_->addToBlock(globJ, globI, bMat); - loc++; + ++loc; } } if (not(well_local)) { @@ -472,8 +485,13 @@ private: LinearizationType linearizationType_; - SparseTable neighbours_; - SparseTable trans_; + struct NeighborInfo + { + unsigned int neighbor; + double trans; + double faceArea; + }; + SparseTable neighborInfo_; }; } // namespace Opm