From 8f5e0940fe02b2567dead0c6e713767a890294f2 Mon Sep 17 00:00:00 2001 From: hnil Date: Fri, 17 Jun 2022 11:15:15 +0200 Subject: [PATCH] restructuring to be able to call without local indices --- ebos/eclfluxmoduletpfa.hh | 115 ++++++++++++++++-- ebos/eclthresholdpressure.hh | 1 - opm/simulators/flow/BlackoilModelEbos.hpp | 54 ++++++-- opm/simulators/wells/BlackoilWellModel.hpp | 6 + .../wells/BlackoilWellModel_impl.hpp | 16 +++ 5 files changed, 171 insertions(+), 21 deletions(-) diff --git a/ebos/eclfluxmoduletpfa.hh b/ebos/eclfluxmoduletpfa.hh index e4b3cc52a..9fbf40059 100644 --- a/ebos/eclfluxmoduletpfa.hh +++ b/ebos/eclfluxmoduletpfa.hh @@ -219,7 +219,6 @@ public: EvalType& pressureDifference, const IntensiveQuantities& intQuantsIn, const IntensiveQuantities& intQuantsEx, - const unsigned scvfIdx, const unsigned timeIdx, const unsigned phaseIdx, const unsigned interiorDofIdx, @@ -315,6 +314,106 @@ public: } } + // static void volumeAndPhasePressureDifferences(short (&upIdx)[numPhases] , + // Evaluation (&volumeFlux)[numPhases], + // Evaluation (&pressureDifferences)[numPhases], + // const Problem& problem, + // const unsigned globalIndexIn, + // const unsinged globalIndexOut, + // const IntensiveQuantities& intQuantsIn, + // const IntensiveQuantities& intQuantsIn, + // const unsinged timeIdx) + // { + + // //Valgrind::SetUndefined(*this); + + // //const auto& problem = elemCtx.problem(); + // //const auto& stencil = elemCtx.stencil(timeIdx); + // //const auto& scvf = stencil.interiorFace(scvfIdx); + // //unsigned interiorDofIdx = scvf.interiorIndex(); + // //unsigned exteriorDofIdx = scvf.exteriorIndex(); + // //const auto& globalIndexIn = stencil.globalSpaceIndex(interiorDofIdx); + // //const auto& globalIndexEx = stencil.globalSpaceIndex(exteriorDofIdx); + + + // assert(interiorDofIdx != exteriorDofIdx); + + // //unsigned I = stencil.globalSpaceIndex(interiorDofIdx_); + // //unsigned J = stencil.globalSpaceIndex(exteriorDofIdx_); + // Scalar Vin = model.dofTotalVolume(globalIndexIn, /*timeIdx=*/0); + // Scalar Vex = model.dofTotalVolume(globalIndexOut, /*timeIdx=*/0); + + + // Scalar trans = problem.transmissibility(globalIndexIn,globalIndexOut); + // Scalar faceArea = problem.area(globalIndexIn,globalIndexOut); + // Scalar thpres = problem.thresholdPressure(globalIndexIn, globalIndexEx); + + // // estimate the gravity correction: for performance reasons we use a simplified + // // approach for this flux module that assumes that gravity is constant and always + // // acts into the downwards direction. (i.e., no centrifuge experiments, sorry.) + // constexpr Scalar g = 9.8; + + + // // this is quite hacky because the dune grid interface does not provide a + // // cellCenterDepth() method (so we ask the problem to provide it). The "good" + // // solution would be to take the Z coordinate of the element centroids, but since + // // ECL seems to like to be inconsistent on that front, it needs to be done like + // // here... + // Scalar zIn = problem.dofCenterDepth(globalIndexIn, timeIdx); + // Scalar zEx = problem.dofCenterDepth(globalIndexOut, timeIdx); + + // // the distances from the DOF's depths. (i.e., the additional depth of the + // // exterior DOF) + // Scalar distZ = zIn - zEx; + + // for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) { + // if (!FluidSystem::phaseIsActive(phaseIdx)) + // continue; + // short dnIdx; + // calculatePhasePressureDiff_(upIdx[phaseIdx], + // dnIdx, + // pressureDifferences[phaseIdx], + // intQuantsIn, + // intQuantsEx, + // timeIdx,//input + // phaseIdx,//input + // interiorDofIdx,//input + // exteriorDofIdx,//intput + // Vin, + // Vex, + // globalIndexIn, + // globalIndexEx, + // distZ*g, + // thpres); + // if(pressureDifferences[phaseIdx] == 0){ + // volumeFlux[phaseIdx] = 0.0; + // continue; + // } + // IntensiveQuantities up; + // unsigned globalIndex; + // if(upIdx[phaseIdx] == interiorDofIdx){ + // up = intQuantsIn; + // globalIndex = globalIndexIn; + // }else{ + // up = intQuantsEx; + // globalIndex = globalIndexEx; + // } + // // TODO: should the rock compaction transmissibility multiplier be upstreamed + // // or averaged? all fluids should see the same compaction?! + // //const auto& globalIndex = stencil.globalSpaceIndex(upstreamIdx); + // const Evaluation& transMult = + // problem.template rockCompTransMultiplier(up, globalIndex); + + // if (upIdx[phaseIdx] == interiorDofIdx) + // volumeFlux[phaseIdx] = + // pressureDifferences[phaseIdx]*up.mobility(phaseIdx)*transMult*(-trans/faceArea); + // else + // volumeFlux[phaseIdx] = + // pressureDifferences[phaseIdx]*(Toolbox::value(up.mobility(phaseIdx))*Toolbox::value(transMult)*(-trans/faceArea)); + + // } + // } + static void volumeAndPhasePressureDifferences(short (&upIdx)[numPhases] , Evaluation (&volumeFlux)[numPhases], Evaluation (&pressureDifferences)[numPhases], @@ -326,17 +425,20 @@ public: const auto& problem = elemCtx.problem(); const auto& stencil = elemCtx.stencil(timeIdx); const auto& scvf = stencil.interiorFace(scvfIdx); - unsigned interiorDofIdx = scvf.interiorIndex(); unsigned exteriorDofIdx = scvf.exteriorIndex(); + const auto& globalIndexIn = stencil.globalSpaceIndex(interiorDofIdx); + const auto& globalIndexEx = stencil.globalSpaceIndex(exteriorDofIdx); + + assert(interiorDofIdx != exteriorDofIdx); //unsigned I = stencil.globalSpaceIndex(interiorDofIdx_); //unsigned J = stencil.globalSpaceIndex(exteriorDofIdx_); Scalar Vin = elemCtx.dofVolume(interiorDofIdx, /*timeIdx=*/0); Scalar Vex = elemCtx.dofVolume(exteriorDofIdx, /*timeIdx=*/0); - const auto& globalIndexIn = stencil.globalSpaceIndex(interiorDofIdx); - const auto& globalIndexEx = stencil.globalSpaceIndex(exteriorDofIdx); + + Scalar trans = problem.transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx); Scalar faceArea = scvf.area(); Scalar thpres = problem.thresholdPressure(globalIndexIn, globalIndexEx); @@ -344,7 +446,7 @@ public: // estimate the gravity correction: for performance reasons we use a simplified // approach for this flux module that assumes that gravity is constant and always // acts into the downwards direction. (i.e., no centrifuge experiments, sorry.) - Scalar g = elemCtx.problem().gravity()[dimWorld - 1]; + constexpr Scalar g = 9.8; const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx, timeIdx); const auto& intQuantsEx = elemCtx.intensiveQuantities(exteriorDofIdx, timeIdx); @@ -370,7 +472,6 @@ public: pressureDifferences[phaseIdx], intQuantsIn, intQuantsEx, - scvfIdx,//input timeIdx,//input phaseIdx,//input interiorDofIdx,//input @@ -465,7 +566,6 @@ protected: pressureDifference_[phaseIdx], intQuantsIn, intQuantsEx, - scvfIdx,//input timeIdx,//input phaseIdx,//input interiorDofIdx_,//input @@ -492,6 +592,7 @@ protected: // TODO: should the rock compaction transmissibility multiplier be upstreamed // or averaged? all fluids should see the same compaction?! //const auto& globalIndex = stencil.globalSpaceIndex(upstreamIdx); + //NB as long as this is upwinded it could be an intensive quantity const Evaluation& transMult = problem.template rockCompTransMultiplier(up, globalIndex); diff --git a/ebos/eclthresholdpressure.hh b/ebos/eclthresholdpressure.hh index d86a0c19e..8154f3d3c 100644 --- a/ebos/eclthresholdpressure.hh +++ b/ebos/eclthresholdpressure.hh @@ -170,7 +170,6 @@ private: pressureDifference, intQuantsIn, intQuantsEx, - scvfIdx,//input /*timeIdx*/0,//input phaseIdx,//input i,//input diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index 21d58aaf8..d01e11f65 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -574,20 +574,54 @@ namespace Opm { template void updateIntensiveQuantity(const Problem& /*problem*/, - const SolutionVector& /*solution*/, + SolutionVector& solution, + const BVector& dx, const IntensiveQuantity*) { + ebosSimulator_.model().newtonMethod().update_(/*nextSolution=*/solution, + /*curSolution=*/solution, + /*update=*/dx, + /*resid=*/dx); // the update routines of the black + // oil model do not care about the + // residual ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0); } void updateIntensiveQuantity(const Problem& problem, - const SolutionVector& solution, + SolutionVector& solution, + const BVector& dx, const BlackOilIntensiveQuantitiesSimple* /*intensive*/) { - ebosSimulator_.model().invalidateAndUpdateIntensiveQuantitiesSimple(problem, - solution, - /*timeIdx*/0); + auto& model = ebosSimulator_.model(); + auto& ebosNewtonMethod = model.newtonMethod(); + if(false){ + if (!std::isfinite(dx.one_norm())) + throw NumericalIssue("Non-finite update!"); + + size_t numGridDof = model.numGridDof(); + for (unsigned dofIdx = 0; dofIdx < numGridDof; ++dofIdx) { + ebosNewtonMethod.updatePrimaryVariables_(dofIdx, + solution[dofIdx], + solution[dofIdx], + dx[dofIdx], + dx[dofIdx]); + model.invalidateAndUpdateIntensiveSingleQuantitiesSimple(problem, + solution[dofIdx], + dofIdx, + /*timeIdx*/0); + } + }else{ + ebosNewtonMethod.update_(/*nextSolution*/solution, + /*curSolution=*/solution, + /*update=*/dx, + /*resid=*/dx); // the update routines of the black + // oil model do not care about the + // residual + model.invalidateAndUpdateIntensiveQuantitiesSimple(problem, + solution, + /*timeIdx*/0); + } } /// Apply an update to the primary variables. @@ -597,17 +631,11 @@ namespace Opm { auto& ebosNewtonMethod = ebosSimulator_.model().newtonMethod(); SolutionVector& solution = ebosSimulator_.model().solution(/*timeIdx=*/0); - ebosNewtonMethod.update_(/*nextSolution=*/solution, - /*curSolution=*/solution, - /*update=*/dx, - /*resid=*/dx); // the update routines of the black - // oil model do not care about the - // residual - + // if the solution is updated, the intensive quantities need to be recalculated //ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0); IntensiveQuantities* dummy = NULL;//only for template spesialization - this->updateIntensiveQuantity(problem,solution, dummy); + this->updateIntensiveQuantity(problem,solution, dx, dummy); //ebosSimulator_.model().invalidateAndUpdateIntensiveQuantitiesSimple(problem, //solution, // /*timeIdx=*/0); diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index 13cd0d42a..ed443600e 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -210,7 +210,13 @@ namespace Opm { { endReportStep(); } + + + void computeTotalRatesForDof(RateVector& rate, + unsigned globalIdx, + unsigned timeIdx) const; + template void computeTotalRatesForDof(RateVector& rate, const Context& context, diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index ebde65f63..b28fa609e 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -499,7 +499,23 @@ namespace Opm { this->computeWellTemperature(); } + template + void + BlackoilWellModel:: + computeTotalRatesForDof(RateVector& rate, + unsigned elemIdx, + unsigned timeIdx) const + { + rate = 0; + + if (!is_cell_perforated_[elemIdx]) + return; + for (const auto& well : well_container_) + well->addCellRates(rate, elemIdx); + } + + template template void