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);
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().anyFlows();
543 const bool anyFlores = simulator_().
problem().eclWriter()->outputModule().anyFlores();
544 const bool enableDispersion = simulator_().
vanguard().eclState().getSimulationConfig().rock_config().dispersion();
545 if (((!anyFlows || !flowsInfo_.empty()) && (!anyFlores || !floresInfo_.empty())) && !enableDispersion) {
548 const auto& model = model_();
550 Stencil stencil(gridView_(), model_().dofMapper());
551 unsigned numCells = model.numTotalDof();
552 std::unordered_multimap<int, std::pair<int, int>>
nncIndices;
555 unsigned int nncId = 0;
556 VectorBlock flow(0.0);
566 flowsInfo_.reserve(numCells, 6 * numCells);
569 floresInfo_.reserve(numCells, 6 * numCells);
571 if (enableDispersion) {
572 velocityInfo_.reserve(numCells, 6 * numCells);
576 stencil.update(
elem);
579 int numFaces = stencil.numBoundaryFaces() + stencil.numInteriorFaces();
588 int faceId =
scvf.dirId();
592 for (
auto it =
range.first; it !=
range.second; ++it) {
597 nncId = it->second.second;
606 const auto&
scvf = stencil.boundaryFace(
bdfIdx);
607 int faceId =
scvf.dirId();
608 loc_flinfo[stencil.numInteriorFaces() +
bdfIdx] = FlowInfo{faceId, flow, nncId};
618 if (enableDispersion) {
642 void updateFlowsInfo() {
644 const bool&
enableFlows = simulator_().
problem().eclWriter()->outputModule().hasFlows() ||
645 simulator_().
problem().eclWriter()->outputModule().hasBlockFlows();
650 const unsigned int numCells = model_().numTotalDof();
652#pragma omp parallel for
659 const IntensiveQuantities&
intQuantsIn = model_().intensiveQuantities(
globI, 0);
670 const IntensiveQuantities&
intQuantsEx = model_().intensiveQuantities(
globJ, 0);
689 for (
const auto&
bdyInfo : boundaryInfo_) {
690 if (
bdyInfo.bcdata.type == BCType::NONE)
699 const unsigned bfIndex =
bdyInfo.bfIndex;
710 template <
class SubDomainType>
716 if (!problem_().recycleFirstIterationStorage()) {
717 if (!model_().storeIntensiveQuantities() && !model_().enableStorageCache()) {
718 OPM_THROW(std::runtime_error,
"Must have cached either IQs or storage when we cannot recycle.");
727 const bool& enableDispersion = simulator_().
vanguard().eclState().getSimulationConfig().rock_config().dispersion();
728 const unsigned int numCells =
domain.cells.size();
732#pragma omp parallel for
734 for (
unsigned ii = 0;
ii < numCells; ++
ii) {
738 VectorBlock
res(0.0);
739 MatrixBlock
bMat(0.0);
742 const IntensiveQuantities&
intQuantsIn = model_().intensiveQuantities(
globI, 0);
756 const IntensiveQuantities&
intQuantsEx = model_().intensiveQuantities(
globJ, 0);
759 if (enableDispersion) {
777 double volume = model_().dofTotalVolume(
globI);
786 if (model_().enableStorageCache()) {
790 model_().updateCachedStorage(
globI, 0,
res);
791 if (model_().newtonMethod().numIterations() == 0) {
793 if (problem_().recycleFirstIterationStorage()) {
803 model_().updateCachedStorage(
globI, 1,
res);
806 Dune::FieldVector<Scalar, numEq> tmp;
809 model_().updateCachedStorage(
globI, 1, tmp);
812 res -= model_().cachedStorage(
globI, 1);
815 Dune::FieldVector<Scalar, numEq> tmp;
832 if (separateSparseSourceTerms_) {
833 LocalResidual::computeSourceDense(
adres, problem_(),
globI, 0);
835 LocalResidual::computeSource(
adres, problem_(),
globI, 0);
845 if (separateSparseSourceTerms_) {
846 problem_().wellModel().addReservoirSourceTerms(residual_, diagMatAddress_);
850 for (
const auto&
bdyInfo : boundaryInfo_) {
851 if (
bdyInfo.bcdata.type == BCType::NONE)
854 VectorBlock
res(0.0);
855 MatrixBlock
bMat(0.0);
868 void updateStoredTransmissibilities()
870 if (neighborInfo_.empty()) {
874 initFirstIteration_();
876 unsigned numCells = model_().numTotalDof();
878#pragma omp parallel for
890 Simulator *simulatorPtr_;
893 std::unique_ptr<SparseMatrixAdapter> jacobian_;
896 GlobalEqVector residual_;
898 LinearizationType linearizationType_;
900 using ResidualNBInfo =
typename LocalResidual::ResidualNBInfo;
903 unsigned int neighbor;
904 ResidualNBInfo res_nbinfo;
905 MatrixBlock* matBlockAddress;
908 std::vector<MatrixBlock*> diagMatAddress_;
921 VectorBlock velocity;
925 using ScalarFluidState =
typename IntensiveQuantities::ScalarFluidState;
926 struct BoundaryConditionData
929 VectorBlock massRate;
930 unsigned pvtRegionIdx;
931 unsigned boundaryFaceIndex;
934 ScalarFluidState exFluidState;
940 unsigned int bfIndex;
941 BoundaryConditionData bcdata;
943 std::vector<BoundaryInfo> boundaryInfo_;
944 bool separateSparseSourceTerms_ =
false;
947 std::vector<int> cells;
948 std::vector<bool> interior;
950 FullDomain fullDomain_;