95 using Element =
typename GridView::template Codim<0>::Entity;
96 using ElementIterator =
typename GridView::template Codim<0>::Iterator;
98 using Vector = GlobalEqVector;
102 enum { dimWorld = GridView::dimensionworld };
104 using MatrixBlock =
typename SparseMatrixAdapter::MatrixBlock;
105 using VectorBlock = Dune::FieldVector<Scalar, numEq>;
122 separateSparseSourceTerms_ = Parameters::Get<Parameters::SeparateSparseSourceTerms>();
134 Parameters::Register<Parameters::SeparateSparseSourceTerms>
135 (
"Treat well source terms all in one go, instead of on a cell by cell basis.");
149 simulatorPtr_ = &simulator;
197 catch (
const std::exception&
e)
199 std::cout <<
"rank " << simulator_().
gridView().comm().rank()
200 <<
" caught an exception while linearizing:" <<
e.what()
201 <<
"\n" << std::flush;
206 std::cout <<
"rank " << simulator_().
gridView().comm().rank()
207 <<
" caught an exception while linearizing"
208 <<
"\n" << std::flush;
214 throw NumericalProblem(
"A process did not succeed in linearizing the system");
227 template <
class SubDomainType>
235 initFirstIteration_();
238 if (
domain.cells.size() == model_().numTotalDof()) {
249 { jacobian_->finalize(); }
261 auto& model = model_();
262 const auto& comm = simulator_().
gridView().comm();
266 model.auxiliaryModule(
auxModIdx)->linearize(*jacobian_, residual_);
268 catch (
const std::exception&
e) {
271 std::cout <<
"rank " << simulator_().
gridView().comm().rank()
272 <<
" caught an exception while linearizing:" <<
e.what()
273 <<
"\n" << std::flush;
287 {
return *jacobian_; }
290 {
return *jacobian_; }
296 {
return residual_; }
299 {
return residual_; }
301 void setLinearizationType(LinearizationType linearizationType){
302 linearizationType_ = linearizationType;
305 const LinearizationType& getLinearizationType()
const{
306 return linearizationType_;
336 return velocityInfo_;
339 void updateDiscretizationParameters()
341 updateStoredTransmissibilities();
344 void updateBoundaryConditionData() {
345 for (
auto&
bdyInfo : boundaryInfo_) {
353 if (type != BCType::NONE) {
354 const auto& exFluidState = problem_().boundaryFluidState(
bdyInfo.cell,
bdyInfo.dir);
357 bdyInfo.bcdata.exFluidState = exFluidState;
370 template <
class SubDomainType>
374 initFirstIteration_();
377 residual_[
globI] = 0.0;
378 jacobian_->clearRow(
globI, 0.0);
383 Simulator& simulator_()
384 {
return *simulatorPtr_; }
385 const Simulator& simulator_()
const
386 {
return *simulatorPtr_; }
389 {
return simulator_().
problem(); }
390 const Problem& problem_()
const
391 {
return simulator_().
problem(); }
394 {
return simulator_().
model(); }
395 const Model& model_()
const
396 {
return simulator_().
model(); }
398 const GridView& gridView_()
const
399 {
return problem_().gridView(); }
401 void initFirstIteration_()
407 residual_.resize(model_().numTotalDof());
418 if (!neighborInfo_.empty()) {
423 const auto& model = model_();
424 Stencil stencil(gridView_(), model_().dofMapper());
428 using NeighborSet = std::set< unsigned >;
430 const Scalar gravity = problem_().gravity()[dimWorld - 1];
431 unsigned numCells = model.numTotalDof();
432 neighborInfo_.reserve(numCells, 6 * numCells);
435 stencil.update(
elem);
441 for (
unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
442 unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
446 const auto scvfIdx = dofIdx - 1;
448 const Scalar area =
scvf.area();
449 const Scalar Vin = problem_().model().dofTotalVolume(
myIdx);
450 const Scalar Vex = problem_().model().dofTotalVolume(
neighborIdx);
451 const Scalar
zIn = problem_().dofCenterDepth(
myIdx);
453 const Scalar dZg = (
zIn -
zEx)*gravity;
456 Scalar outAlpha {0.};
457 Scalar diffusivity {0.};
458 Scalar dispersivity {0.};
459 if constexpr(enableEnergy){
463 if constexpr(enableDiffusion){
469 const auto dirId =
scvf.dirId();
470 auto faceDir = dirId < 0 ? FaceDir::DirEnum::Unknown
471 : FaceDir::FromIntersectionIndex(dirId);
472 loc_nbinfo[dofIdx - 1] = NeighborInfo{
neighborIdx, {trans, area, thpres, dZg, faceDir, Vin, Vex, inAlpha, outAlpha, diffusivity, dispersivity},
nullptr};
477 if (problem_().nonTrivialBoundaryConditions()) {
478 for (
unsigned bfIndex = 0; bfIndex < stencil.numBoundaryFaces(); ++bfIndex) {
479 const auto&
bf = stencil.boundaryFace(bfIndex);
490 const auto& exFluidState = problem_().boundaryFluidState(
myIdx,
dir_id);
491 BoundaryConditionData bcdata{type,
493 exFluidState.pvtRegionIndex(),
496 bf.integrationPos()[dimWorld - 1],
498 boundaryInfo_.push_back({
myIdx,
dir_id, bfIndex, bcdata});
506 size_t numAuxMod = model.numAuxiliaryModules();
511 jacobian_.reset(
new SparseMatrixAdapter(simulator_()));
512 diagMatAddress_.resize(numCells);
524 fullDomain_.cells.resize(numCells);
525 std::iota(fullDomain_.cells.begin(), fullDomain_.cells.end(), 0);
543 const bool anyFlows = simulator_().
problem().eclWriter().outputModule().getFlows().anyFlows() ||
544 simulator_().
problem().eclWriter().outputModule().getFlows().hasBlockFlows();
545 const bool anyFlores = simulator_().
problem().eclWriter().outputModule().getFlows().anyFlores();
546 const bool enableDispersion = simulator_().
vanguard().eclState().getSimulationConfig().rock_config().dispersion();
547 if (((!anyFlows || !flowsInfo_.empty()) && (!anyFlores || !floresInfo_.empty())) && (!enableDispersion && !enableMICP)) {
550 const auto& model = model_();
552 Stencil stencil(gridView_(), model_().dofMapper());
553 unsigned numCells = model.numTotalDof();
554 std::unordered_multimap<int, std::pair<int, int>>
nncIndices;
557 unsigned int nncId = 0;
558 VectorBlock flow(0.0);
568 flowsInfo_.reserve(numCells, 6 * numCells);
571 floresInfo_.reserve(numCells, 6 * numCells);
573 if (enableDispersion || enableMICP) {
574 velocityInfo_.reserve(numCells, 6 * numCells);
578 stencil.update(
elem);
581 int numFaces = stencil.numBoundaryFaces() + stencil.numInteriorFaces();
585 for (
unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
586 unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
588 const auto scvfIdx = dofIdx - 1;
590 int faceId =
scvf.dirId();
594 for (
auto it =
range.first; it !=
range.second; ++it) {
599 nncId = it->second.second;
602 loc_flinfo[dofIdx - 1] = FlowInfo{faceId, flow, nncId};
608 const auto&
scvf = stencil.boundaryFace(
bdfIdx);
609 int faceId =
scvf.dirId();
610 loc_flinfo[stencil.numInteriorFaces() +
bdfIdx] = FlowInfo{faceId, flow, nncId};
620 if (enableDispersion || enableMICP) {
644 void updateFlowsInfo() {
646 const bool enableFlows = simulator_().
problem().eclWriter().outputModule().getFlows().hasFlows() ||
647 simulator_().
problem().eclWriter().outputModule().getFlows().hasBlockFlows();
648 const bool enableFlores = simulator_().
problem().eclWriter().outputModule().getFlows().hasFlores();
652 const unsigned int numCells = model_().numTotalDof();
654#pragma omp parallel for
661 const IntensiveQuantities&
intQuantsIn = model_().intensiveQuantities(
globI, 0);
672 const IntensiveQuantities&
intQuantsEx = model_().intensiveQuantities(
globJ, 0);
691 for (
const auto&
bdyInfo : boundaryInfo_) {
692 if (
bdyInfo.bcdata.type == BCType::NONE)
701 const unsigned bfIndex =
bdyInfo.bfIndex;
712 template <
class SubDomainType>
718 if (!problem_().recycleFirstIterationStorage()) {
719 if (!model_().storeIntensiveQuantities() && !model_().enableStorageCache()) {
720 OPM_THROW(std::runtime_error,
"Must have cached either IQs or storage when we cannot recycle.");
729 const bool enableDispersion = simulator_().
vanguard().eclState().getSimulationConfig().rock_config().dispersion();
730 const unsigned int numCells =
domain.cells.size();
734#pragma omp parallel for
736 for (
unsigned ii = 0;
ii < numCells; ++
ii) {
740 VectorBlock
res(0.0);
741 MatrixBlock
bMat(0.0);
744 const IntensiveQuantities&
intQuantsIn = model_().intensiveQuantities(
globI, 0);
758 const IntensiveQuantities&
intQuantsEx = model_().intensiveQuantities(
globJ, 0);
761 if (enableDispersion || enableMICP) {
779 double volume = model_().dofTotalVolume(
globI);
788 if (model_().enableStorageCache()) {
792 model_().updateCachedStorage(
globI, 0,
res);
793 if (model_().newtonMethod().numIterations() == 0) {
795 if (problem_().recycleFirstIterationStorage()) {
805 model_().updateCachedStorage(
globI, 1,
res);
808 Dune::FieldVector<Scalar, numEq> tmp;
811 model_().updateCachedStorage(
globI, 1, tmp);
814 res -= model_().cachedStorage(
globI, 1);
817 Dune::FieldVector<Scalar, numEq> tmp;
834 if (separateSparseSourceTerms_) {
847 if (separateSparseSourceTerms_) {
848 problem_().wellModel().addReservoirSourceTerms(residual_, diagMatAddress_);
852 for (
const auto&
bdyInfo : boundaryInfo_) {
853 if (
bdyInfo.bcdata.type == BCType::NONE)
856 VectorBlock
res(0.0);
857 MatrixBlock
bMat(0.0);
870 void updateStoredTransmissibilities()
872 if (neighborInfo_.empty()) {
876 initFirstIteration_();
878 unsigned numCells = model_().numTotalDof();
880#pragma omp parallel for
892 Simulator *simulatorPtr_;
895 std::unique_ptr<SparseMatrixAdapter> jacobian_;
898 GlobalEqVector residual_;
900 LinearizationType linearizationType_;
902 using ResidualNBInfo =
typename LocalResidual::ResidualNBInfo;
905 unsigned int neighbor;
906 ResidualNBInfo res_nbinfo;
907 MatrixBlock* matBlockAddress;
909 SparseTable<NeighborInfo> neighborInfo_;
910 std::vector<MatrixBlock*> diagMatAddress_;
918 SparseTable<FlowInfo> flowsInfo_;
919 SparseTable<FlowInfo> floresInfo_;
923 VectorBlock velocity;
925 SparseTable<VelocityInfo> velocityInfo_;
927 using ScalarFluidState =
typename IntensiveQuantities::ScalarFluidState;
928 struct BoundaryConditionData
931 VectorBlock massRate;
932 unsigned pvtRegionIdx;
933 unsigned boundaryFaceIndex;
936 ScalarFluidState exFluidState;
942 unsigned int bfIndex;
943 BoundaryConditionData bcdata;
945 std::vector<BoundaryInfo> boundaryInfo_;
946 bool separateSparseSourceTerms_ =
false;
949 std::vector<int> cells;
950 std::vector<bool> interior;
952 FullDomain fullDomain_;