From f6b5e27e5273e9aab3bfe28c150368c736dea401 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Tue, 5 Dec 2023 13:06:11 +0100 Subject: [PATCH] factor out flux output --- .../blackoil/blackoilintensivequantities.hh | 3 + .../blackoil/blackoillocalresidualtpfa.hh | 20 +--- .../discretization/common/fvbaselinearizer.hh | 4 + .../discretization/common/tpfalinearizer.hh | 94 ++++++++++++------- 4 files changed, 71 insertions(+), 50 deletions(-) diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index 8a4e19878..7912b48cc 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -515,10 +515,13 @@ public: using Dir = FaceDir::DirEnum; if (dirMob_) { switch(facedir) { + case Dir::XMinus: case Dir::XPlus: return dirMob_->mobilityX_[phaseIdx]; + case Dir::YMinus: case Dir::YPlus: return dirMob_->mobilityY_[phaseIdx]; + case Dir::ZMinus: case Dir::ZPlus: return dirMob_->mobilityZ_[phaseIdx]; default: diff --git a/opm/models/blackoil/blackoillocalresidualtpfa.hh b/opm/models/blackoil/blackoillocalresidualtpfa.hh index 17654a9c8..e15e455a0 100644 --- a/opm/models/blackoil/blackoillocalresidualtpfa.hh +++ b/opm/models/blackoil/blackoillocalresidualtpfa.hh @@ -272,8 +272,6 @@ public: const auto& globalIndexEx = stencil.globalSpaceIndex(exteriorDofIdx); Scalar trans = problem.transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx); Scalar faceArea = scvf.area(); - const auto& materialLawManager = problem.materialLawManager(); - const auto dirid = scvf.dirId(); Scalar thpres = problem.thresholdPressure(globalIndexIn, globalIndexEx); @@ -789,25 +787,11 @@ public: static FaceDir::DirEnum faceDirFromDirId(const int dirId) { - using Dir = FaceDir::DirEnum; // NNC does not have a direction if (dirId < 0 ) { - return Dir::Unknown; - } - switch(dirId) { - case 0: - case 1: - return Dir::XPlus; - case 2: - case 3: - return Dir::YPlus; - case 4: - case 5: - return Dir::ZPlus; - default: - OPM_THROW(std::runtime_error, - "Unexpected face id" + std::to_string(dirId)); + return FaceDir::DirEnum::Unknown; } + return FaceDir::FromIntersectionIndex(dirId); } }; diff --git a/opm/models/discretization/common/fvbaselinearizer.hh b/opm/models/discretization/common/fvbaselinearizer.hh index e36a7d272..c53035a46 100644 --- a/opm/models/discretization/common/fvbaselinearizer.hh +++ b/opm/models/discretization/common/fvbaselinearizer.hh @@ -311,6 +311,10 @@ public: // This linearizer stores no such data. } + void updateFlowsInfo() { + // This linearizer stores no such data. + } + /*! * \brief Returns the map of constraint degrees of freedom. * diff --git a/opm/models/discretization/common/tpfalinearizer.hh b/opm/models/discretization/common/tpfalinearizer.hh index 56d23ded9..f06c21b4e 100644 --- a/opm/models/discretization/common/tpfalinearizer.hh +++ b/opm/models/discretization/common/tpfalinearizer.hh @@ -434,8 +434,6 @@ private: unsigned numCells = model.numTotalDof(); neighborInfo_.reserve(numCells, 6 * numCells); std::vector loc_nbinfo; - const auto& materialLawManager = problem_().materialLawManager(); - using FaceDirection = FaceDir::DirEnum; for (const auto& elem : elements(gridView_())) { stencil.update(elem); @@ -643,6 +641,68 @@ public: } } + void updateFlowsInfo() { + OPM_TIMEBLOCK(updateFlows); + const bool& enableFlows = simulator_().problem().eclWriter()->eclOutputModule().hasFlows() || + simulator_().problem().eclWriter()->eclOutputModule().hasBlockFlows(); + const bool& enableFlores = simulator_().problem().eclWriter()->eclOutputModule().hasFlores(); + if (!enableFlows && !enableFlores) { + return; + } + const unsigned int numCells = model_().numTotalDof(); + +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (unsigned globI = 0; globI < numCells; ++globI) { + OPM_TIMEBLOCK_LOCAL(linearizationForEachCell); + const auto& nbInfos = neighborInfo_[globI]; + ADVectorBlock adres(0.0); + ADVectorBlock darcyFlux(0.0); + const IntensiveQuantities& intQuantsIn = model_().intensiveQuantities(globI, /*timeIdx*/ 0); + // Flux term. + { + OPM_TIMEBLOCK_LOCAL(fluxCalculationForEachCell); + for (const auto& nbInfo : nbInfos) { + OPM_TIMEBLOCK_LOCAL(fluxCalculationForEachFace); + unsigned globJ = nbInfo.neighbor; + assert(globJ != globI); + adres = 0.0; + darcyFlux = 0.0; + const IntensiveQuantities& intQuantsEx = model_().intensiveQuantities(globJ, /*timeIdx*/ 0); + LocalResidual::computeFlux(adres,darcyFlux, globI, globJ, intQuantsIn, intQuantsEx, nbInfo.res_nbinfo); + adres *= nbInfo.res_nbinfo.faceArea; + if (enableFlows) { + for (unsigned phaseIdx = 0; phaseIdx < numEq; ++ phaseIdx) { + flowsInfo_[globI][nbInfo.res_nbinfo.dirId].flow[phaseIdx] = adres[phaseIdx].value(); + } + } + if (enableFlores) { + for (unsigned phaseIdx = 0; phaseIdx < numEq; ++ phaseIdx) { + floresInfo_[globI][nbInfo.res_nbinfo.dirId].flow[phaseIdx] = darcyFlux[phaseIdx].value(); + } + } + } + } + } + + // Boundary terms. Only looping over cells with nontrivial bcs. + for (const auto& bdyInfo : boundaryInfo_) { + ADVectorBlock adres(0.0); + const unsigned globI = bdyInfo.cell; + const IntensiveQuantities& insideIntQuants = model_().intensiveQuantities(globI, /*timeIdx*/ 0); + LocalResidual::computeBoundaryFlux(adres, problem_(), bdyInfo.bcdata, insideIntQuants, globI); + adres *= bdyInfo.bcdata.faceArea; + + if (enableFlows) { + for (unsigned phaseIdx = 0; phaseIdx < numEq; ++ phaseIdx) { + flowsInfo_[globI][bdyInfo.dir].flow[phaseIdx] = adres[phaseIdx].value(); + } + } + // TODO also store Flores? + } + } + private: template void linearize_(const SubDomainType& domain) @@ -662,9 +722,6 @@ private: // the full system to zero, not just our part. // Instead, that must be called before starting the linearization. const bool& enableDispersion = simulator_().vanguard().eclState().getSimulationConfig().rock_config().dispersion(); - const bool& enableFlows = simulator_().problem().eclWriter()->eclOutputModule().hasFlows() || - simulator_().problem().eclWriter()->eclOutputModule().hasBlockFlows(); - const bool& enableFlores = simulator_().problem().eclWriter()->eclOutputModule().hasFlores(); const unsigned int numCells = domain.cells.size(); const bool on_full_domain = (numCells == model_().numTotalDof()); @@ -701,16 +758,6 @@ private: velocityInfo_[globI][loc].velocity[phaseIdx] = darcyFlux[phaseIdx].value() / nbInfo.res_nbinfo.faceArea; } } - if (enableFlows) { - for (unsigned phaseIdx = 0; phaseIdx < numEq; ++ phaseIdx) { - flowsInfo_[globI][nbInfo.res_nbinfo.dirId].flow[phaseIdx] = adres[phaseIdx].value(); - } - } - if (enableFlores) { - for (unsigned phaseIdx = 0; phaseIdx < numEq; ++ phaseIdx) { - floresInfo_[globI][nbInfo.res_nbinfo.dirId].flow[phaseIdx] = darcyFlux[phaseIdx].value(); - } - } setResAndJacobi(res, bMat, adres); residual_[globI] += res; //SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat); @@ -809,23 +856,6 @@ private: residual_[globI] += res; ////SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat); *diagMatAddress_[globI] += bMat; - - if (enableFlows) { - for (unsigned phaseIdx = 0; phaseIdx < numEq; ++ phaseIdx) { - if (adres[phaseIdx].value() == 0) - continue; - - flowsInfo_[globI][bdyInfo.dir].flow[phaseIdx] = adres[phaseIdx].value(); - } - } - // if (enableFlores) { - // for (unsigned phaseIdx = 0; phaseIdx < numEq; ++ phaseIdx) { - // if (adres[phaseIdx].value() == 0) - // continue; - - // floresInfo_[globI][bdyInfo.dir].flow[phaseIdx] = adres[phaseIdx].value(); - // } - // } } }