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>;
121 separateSparseSourceTerms_ = Parameters::Get<Parameters::SeparateSparseSourceTerms>();
133 Parameters::Register<Parameters::SeparateSparseSourceTerms>
134 (
"Treat well source terms all in one go, instead of on a cell by cell basis.");
148 simulatorPtr_ = &simulator;
196 catch (
const std::exception&
e)
198 std::cout <<
"rank " << simulator_().
gridView().comm().rank()
199 <<
" caught an exception while linearizing:" <<
e.what()
200 <<
"\n" << std::flush;
205 std::cout <<
"rank " << simulator_().
gridView().comm().rank()
206 <<
" caught an exception while linearizing"
207 <<
"\n" << std::flush;
213 throw NumericalProblem(
"A process did not succeed in linearizing the system");
226 template <
class SubDomainType>
234 initFirstIteration_();
237 if (
domain.cells.size() == model_().numTotalDof()) {
248 { jacobian_->finalize(); }
260 auto& model = model_();
261 const auto& comm = simulator_().
gridView().comm();
265 model.auxiliaryModule(
auxModIdx)->linearize(*jacobian_, residual_);
267 catch (
const std::exception&
e) {
270 std::cout <<
"rank " << simulator_().
gridView().comm().rank()
271 <<
" caught an exception while linearizing:" <<
e.what()
272 <<
"\n" << std::flush;
286 {
return *jacobian_; }
289 {
return *jacobian_; }
295 {
return residual_; }
298 {
return residual_; }
300 void setLinearizationType(LinearizationType linearizationType){
301 linearizationType_ = linearizationType;
304 const LinearizationType& getLinearizationType()
const{
305 return linearizationType_;
335 return velocityInfo_;
338 void updateDiscretizationParameters()
340 updateStoredTransmissibilities();
343 void updateBoundaryConditionData() {
344 for (
auto&
bdyInfo : boundaryInfo_) {
352 if (type != BCType::NONE) {
353 const auto& exFluidState = problem_().boundaryFluidState(
bdyInfo.cell,
bdyInfo.dir);
356 bdyInfo.bcdata.exFluidState = exFluidState;
369 template <
class SubDomainType>
373 initFirstIteration_();
376 residual_[
globI] = 0.0;
377 jacobian_->clearRow(
globI, 0.0);
382 Simulator& simulator_()
383 {
return *simulatorPtr_; }
384 const Simulator& simulator_()
const
385 {
return *simulatorPtr_; }
388 {
return simulator_().
problem(); }
389 const Problem& problem_()
const
390 {
return simulator_().
problem(); }
393 {
return simulator_().
model(); }
394 const Model& model_()
const
395 {
return simulator_().
model(); }
397 const GridView& gridView_()
const
398 {
return problem_().gridView(); }
400 void initFirstIteration_()
406 residual_.resize(model_().numTotalDof());
417 if (!neighborInfo_.empty()) {
422 const auto& model = model_();
423 Stencil stencil(gridView_(), model_().dofMapper());
427 using NeighborSet = std::set< unsigned >;
429 const Scalar gravity = problem_().gravity()[dimWorld - 1];
430 unsigned numCells = model.numTotalDof();
431 neighborInfo_.reserve(numCells, 6 * numCells);
434 stencil.update(
elem);
440 for (
unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
441 unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
445 const auto scvfIdx = dofIdx - 1;
447 const Scalar area =
scvf.area();
448 const Scalar Vin = problem_().model().dofTotalVolume(
myIdx);
449 const Scalar Vex = problem_().model().dofTotalVolume(
neighborIdx);
450 const Scalar
zIn = problem_().dofCenterDepth(
myIdx);
452 const Scalar dZg = (
zIn -
zEx)*gravity;
455 Scalar outAlpha {0.};
456 Scalar diffusivity {0.};
457 Scalar dispersivity {0.};
458 if constexpr(enableEnergy){
462 if constexpr(enableDiffusion){
468 const auto dirId =
scvf.dirId();
469 auto faceDir = dirId < 0 ? FaceDir::DirEnum::Unknown
470 : FaceDir::FromIntersectionIndex(dirId);
471 loc_nbinfo[dofIdx - 1] = NeighborInfo{
neighborIdx, {trans, area, thpres, dZg, faceDir, Vin, Vex, inAlpha, outAlpha, diffusivity, dispersivity},
nullptr};
476 if (problem_().nonTrivialBoundaryConditions()) {
477 for (
unsigned bfIndex = 0; bfIndex < stencil.numBoundaryFaces(); ++bfIndex) {
478 const auto&
bf = stencil.boundaryFace(bfIndex);
489 const auto& exFluidState = problem_().boundaryFluidState(
myIdx,
dir_id);
490 BoundaryConditionData bcdata{type,
492 exFluidState.pvtRegionIndex(),
495 bf.integrationPos()[dimWorld - 1],
497 boundaryInfo_.push_back({
myIdx,
dir_id, bfIndex, bcdata});
505 size_t numAuxMod = model.numAuxiliaryModules();
510 jacobian_.reset(
new SparseMatrixAdapter(simulator_()));
511 diagMatAddress_.resize(numCells);
523 fullDomain_.cells.resize(numCells);
524 std::iota(fullDomain_.cells.begin(), fullDomain_.cells.end(), 0);
542 const bool anyFlows = simulator_().
problem().eclWriter().outputModule().getFlows().anyFlows() ||
543 simulator_().
problem().eclWriter().outputModule().getFlows().hasBlockFlows();
544 const bool anyFlores = simulator_().
problem().eclWriter().outputModule().getFlows().anyFlores();
545 const bool enableDispersion = simulator_().
vanguard().eclState().getSimulationConfig().rock_config().dispersion();
546 if (((!anyFlows || !flowsInfo_.empty()) && (!anyFlores || !floresInfo_.empty())) && !enableDispersion) {
549 const auto& model = model_();
551 Stencil stencil(gridView_(), model_().dofMapper());
552 unsigned numCells = model.numTotalDof();
553 std::unordered_multimap<int, std::pair<int, int>>
nncIndices;
556 unsigned int nncId = 0;
557 VectorBlock flow(0.0);
567 flowsInfo_.reserve(numCells, 6 * numCells);
570 floresInfo_.reserve(numCells, 6 * numCells);
572 if (enableDispersion) {
573 velocityInfo_.reserve(numCells, 6 * numCells);
577 stencil.update(
elem);
580 int numFaces = stencil.numBoundaryFaces() + stencil.numInteriorFaces();
584 for (
unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
585 unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
587 const auto scvfIdx = dofIdx - 1;
589 int faceId =
scvf.dirId();
593 for (
auto it =
range.first; it !=
range.second; ++it) {
598 nncId = it->second.second;
601 loc_flinfo[dofIdx - 1] = FlowInfo{faceId, flow, nncId};
607 const auto&
scvf = stencil.boundaryFace(
bdfIdx);
608 int faceId =
scvf.dirId();
609 loc_flinfo[stencil.numInteriorFaces() +
bdfIdx] = FlowInfo{faceId, flow, nncId};
619 if (enableDispersion) {
643 void updateFlowsInfo() {
645 const bool enableFlows = simulator_().
problem().eclWriter().outputModule().getFlows().hasFlows() ||
646 simulator_().
problem().eclWriter().outputModule().getFlows().hasBlockFlows();
647 const bool enableFlores = simulator_().
problem().eclWriter().outputModule().getFlows().hasFlores();
651 const unsigned int numCells = model_().numTotalDof();
653#pragma omp parallel for
660 const IntensiveQuantities&
intQuantsIn = model_().intensiveQuantities(
globI, 0);
671 const IntensiveQuantities&
intQuantsEx = model_().intensiveQuantities(
globJ, 0);
690 for (
const auto&
bdyInfo : boundaryInfo_) {
691 if (
bdyInfo.bcdata.type == BCType::NONE)
700 const unsigned bfIndex =
bdyInfo.bfIndex;
711 template <
class SubDomainType>
717 if (!problem_().recycleFirstIterationStorage()) {
718 if (!model_().storeIntensiveQuantities() && !model_().enableStorageCache()) {
719 OPM_THROW(std::runtime_error,
"Must have cached either IQs or storage when we cannot recycle.");
728 const bool enableDispersion = simulator_().
vanguard().eclState().getSimulationConfig().rock_config().dispersion();
729 const unsigned int numCells =
domain.cells.size();
733#pragma omp parallel for
735 for (
unsigned ii = 0;
ii < numCells; ++
ii) {
739 VectorBlock
res(0.0);
740 MatrixBlock
bMat(0.0);
743 const IntensiveQuantities&
intQuantsIn = model_().intensiveQuantities(
globI, 0);
757 const IntensiveQuantities&
intQuantsEx = model_().intensiveQuantities(
globJ, 0);
760 if (enableDispersion) {
778 double volume = model_().dofTotalVolume(
globI);
787 if (model_().enableStorageCache()) {
791 model_().updateCachedStorage(
globI, 0,
res);
792 if (model_().newtonMethod().numIterations() == 0) {
794 if (problem_().recycleFirstIterationStorage()) {
804 model_().updateCachedStorage(
globI, 1,
res);
807 Dune::FieldVector<Scalar, numEq> tmp;
810 model_().updateCachedStorage(
globI, 1, tmp);
813 res -= model_().cachedStorage(
globI, 1);
816 Dune::FieldVector<Scalar, numEq> tmp;
833 if (separateSparseSourceTerms_) {
834 LocalResidual::computeSourceDense(
adres, problem_(),
globI, 0);
836 LocalResidual::computeSource(
adres, problem_(),
globI, 0);
846 if (separateSparseSourceTerms_) {
847 problem_().wellModel().addReservoirSourceTerms(residual_, diagMatAddress_);
851 for (
const auto&
bdyInfo : boundaryInfo_) {
852 if (
bdyInfo.bcdata.type == BCType::NONE)
855 VectorBlock
res(0.0);
856 MatrixBlock
bMat(0.0);
869 void updateStoredTransmissibilities()
871 if (neighborInfo_.empty()) {
875 initFirstIteration_();
877 unsigned numCells = model_().numTotalDof();
879#pragma omp parallel for
891 Simulator *simulatorPtr_;
894 std::unique_ptr<SparseMatrixAdapter> jacobian_;
897 GlobalEqVector residual_;
899 LinearizationType linearizationType_;
901 using ResidualNBInfo =
typename LocalResidual::ResidualNBInfo;
904 unsigned int neighbor;
905 ResidualNBInfo res_nbinfo;
906 MatrixBlock* matBlockAddress;
908 SparseTable<NeighborInfo> neighborInfo_;
909 std::vector<MatrixBlock*> diagMatAddress_;
917 SparseTable<FlowInfo> flowsInfo_;
918 SparseTable<FlowInfo> floresInfo_;
922 VectorBlock velocity;
924 SparseTable<VelocityInfo> velocityInfo_;
926 using ScalarFluidState =
typename IntensiveQuantities::ScalarFluidState;
927 struct BoundaryConditionData
930 VectorBlock massRate;
931 unsigned pvtRegionIdx;
932 unsigned boundaryFaceIndex;
935 ScalarFluidState exFluidState;
941 unsigned int bfIndex;
942 BoundaryConditionData bcdata;
944 std::vector<BoundaryInfo> boundaryInfo_;
945 bool separateSparseSourceTerms_ =
false;
948 std::vector<int> cells;
949 std::vector<bool> interior;
951 FullDomain fullDomain_;