My Project
Loading...
Searching...
No Matches
tpfalinearizer.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*
4 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 2 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
28#ifndef TPFA_LINEARIZER_HH
29#define TPFA_LINEARIZER_HH
30
31#include <dune/common/version.hh>
32#include <dune/common/fvector.hh>
33#include <dune/common/fmatrix.hh>
34
35#include <opm/common/Exceptions.hpp>
36#include <opm/common/TimingMacros.hpp>
37
38#include <opm/grid/utility/SparseTable.hpp>
39
40#include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>
41#include <opm/input/eclipse/Schedule/BCProp.hpp>
42
46
47#include <exception> // current_exception, rethrow_exception
48#include <iostream>
49#include <numeric>
50#include <set>
51#include <type_traits>
52#include <vector>
53
54namespace Opm::Parameters {
55
56struct SeparateSparseSourceTerms { static constexpr bool value = false; };
57
58} // namespace Opm::Parameters
59
60namespace Opm {
61
62// forward declarations
63template<class TypeTag>
65
75template<class TypeTag>
77{
85
94
95 using Element = typename GridView::template Codim<0>::Entity;
96 using ElementIterator = typename GridView::template Codim<0>::Iterator;
97
98 using Vector = GlobalEqVector;
99
102 enum { dimWorld = GridView::dimensionworld };
103
104 using MatrixBlock = typename SparseMatrixAdapter::MatrixBlock;
105 using VectorBlock = Dune::FieldVector<Scalar, numEq>;
107
109 static const bool enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>();
110 static const bool enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>();
111
112 // copying the linearizer is not a good idea
113 TpfaLinearizer(const TpfaLinearizer&) = delete;
115
116public:
118 : jacobian_()
119 {
120 simulatorPtr_ = 0;
121 separateSparseSourceTerms_ = Parameters::Get<Parameters::SeparateSparseSourceTerms>();
122 }
123
125 {
126 }
127
131 static void registerParameters()
132 {
133 Parameters::Register<Parameters::SeparateSparseSourceTerms>
134 ("Treat well source terms all in one go, instead of on a cell by cell basis.");
135 }
136
146 void init(Simulator& simulator)
147 {
148 simulatorPtr_ = &simulator;
149 eraseMatrix();
150 }
151
160 {
161 jacobian_.reset();
162 }
163
174 {
177 }
178
190 {
191 int succeeded;
192 try {
193 linearizeDomain(fullDomain_);
194 succeeded = 1;
195 }
196 catch (const std::exception& e)
197 {
198 std::cout << "rank " << simulator_().gridView().comm().rank()
199 << " caught an exception while linearizing:" << e.what()
200 << "\n" << std::flush;
201 succeeded = 0;
202 }
203 catch (...)
204 {
205 std::cout << "rank " << simulator_().gridView().comm().rank()
206 << " caught an exception while linearizing"
207 << "\n" << std::flush;
208 succeeded = 0;
209 }
210 succeeded = simulator_().gridView().comm().min(succeeded);
211
212 if (!succeeded)
213 throw NumericalProblem("A process did not succeed in linearizing the system");
214 }
215
226 template <class SubDomainType>
228 {
230 // we defer the initialization of the Jacobian matrix until here because the
231 // auxiliary modules usually assume the problem, model and grid to be fully
232 // initialized...
233 if (!jacobian_)
234 initFirstIteration_();
235
236 // Called here because it is no longer called from linearize_().
237 if (domain.cells.size() == model_().numTotalDof()) {
238 // We are on the full domain.
239 resetSystem_();
240 } else {
241 resetSystem_(domain);
242 }
243
244 linearize_(domain);
245 }
246
247 void finalize()
248 { jacobian_->finalize(); }
249
255 {
257 // flush possible local caches into matrix structure
258 jacobian_->commit();
259
260 auto& model = model_();
261 const auto& comm = simulator_().gridView().comm();
262 for (unsigned auxModIdx = 0; auxModIdx < model.numAuxiliaryModules(); ++auxModIdx) {
263 bool succeeded = true;
264 try {
265 model.auxiliaryModule(auxModIdx)->linearize(*jacobian_, residual_);
266 }
267 catch (const std::exception& e) {
268 succeeded = false;
269
270 std::cout << "rank " << simulator_().gridView().comm().rank()
271 << " caught an exception while linearizing:" << e.what()
272 << "\n" << std::flush;
273 }
274
275 succeeded = comm.min(succeeded);
276
277 if (!succeeded)
278 throw NumericalProblem("linearization of an auxiliary equation failed");
279 }
280 }
281
285 const SparseMatrixAdapter& jacobian() const
286 { return *jacobian_; }
287
288 SparseMatrixAdapter& jacobian()
289 { return *jacobian_; }
290
294 const GlobalEqVector& residual() const
295 { return residual_; }
296
297 GlobalEqVector& residual()
298 { return residual_; }
299
300 void setLinearizationType(LinearizationType linearizationType){
301 linearizationType_ = linearizationType;
302 };
303
304 const LinearizationType& getLinearizationType() const{
305 return linearizationType_;
306 };
307
313 const auto& getFlowsInfo() const{
314
315 return flowsInfo_;
316 }
317
323 const auto& getFloresInfo() const{
324
325 return floresInfo_;
326 }
327
333 const auto& getVelocityInfo() const{
334
335 return velocityInfo_;
336 }
337
338 void updateDiscretizationParameters()
339 {
340 updateStoredTransmissibilities();
341 }
342
343 void updateBoundaryConditionData() {
344 for (auto& bdyInfo : boundaryInfo_) {
345 const auto [type, massrateAD] = problem_().boundaryCondition(bdyInfo.cell, bdyInfo.dir);
346
347 // Strip the unnecessary (and zero anyway) derivatives off massrate.
348 VectorBlock massrate(0.0);
349 for (size_t ii = 0; ii < massrate.size(); ++ii) {
350 massrate[ii] = massrateAD[ii].value();
351 }
352 if (type != BCType::NONE) {
353 const auto& exFluidState = problem_().boundaryFluidState(bdyInfo.cell, bdyInfo.dir);
354 bdyInfo.bcdata.type = type;
355 bdyInfo.bcdata.massRate = massrate;
356 bdyInfo.bcdata.exFluidState = exFluidState;
357 }
358 }
359 }
360
366 const std::map<unsigned, Constraints> constraintsMap() const
367 { return {}; }
368
369 template <class SubDomainType>
370 void resetSystem_(const SubDomainType& domain)
371 {
372 if (!jacobian_) {
373 initFirstIteration_();
374 }
375 for (int globI : domain.cells) {
376 residual_[globI] = 0.0;
377 jacobian_->clearRow(globI, 0.0);
378 }
379 }
380
381private:
382 Simulator& simulator_()
383 { return *simulatorPtr_; }
384 const Simulator& simulator_() const
385 { return *simulatorPtr_; }
386
387 Problem& problem_()
388 { return simulator_().problem(); }
389 const Problem& problem_() const
390 { return simulator_().problem(); }
391
392 Model& model_()
393 { return simulator_().model(); }
394 const Model& model_() const
395 { return simulator_().model(); }
396
397 const GridView& gridView_() const
398 { return problem_().gridView(); }
399
400 void initFirstIteration_()
401 {
402 // initialize the BCRS matrix for the Jacobian of the residual function
403 createMatrix_();
404
405 // initialize the Jacobian matrix and the vector for the residual function
406 residual_.resize(model_().numTotalDof());
407 resetSystem_();
408
409 // initialize the sparse tables for Flows and Flores
410 createFlows_();
411 }
412
413 // Construct the BCRS matrix for the Jacobian of the residual function
414 void createMatrix_()
415 {
417 if (!neighborInfo_.empty()) {
418 // It is ok to call this function multiple times, but it
419 // should not do anything if already called.
420 return;
421 }
422 const auto& model = model_();
423 Stencil stencil(gridView_(), model_().dofMapper());
424
425 // for the main model, find out the global indices of the neighboring degrees of
426 // freedom of each primary degree of freedom
427 using NeighborSet = std::set< unsigned >;
428 std::vector<NeighborSet> sparsityPattern(model.numTotalDof());
429 const Scalar gravity = problem_().gravity()[dimWorld - 1];
430 unsigned numCells = model.numTotalDof();
431 neighborInfo_.reserve(numCells, 6 * numCells);
432 std::vector<NeighborInfo> loc_nbinfo;
433 for (const auto& elem : elements(gridView_())) {
434 stencil.update(elem);
435
436 for (unsigned primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++primaryDofIdx) {
437 unsigned myIdx = stencil.globalSpaceIndex(primaryDofIdx);
438 loc_nbinfo.resize(stencil.numDof() - 1); // Do not include the primary dof in neighborInfo_
439
440 for (unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
441 unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
443 if (dofIdx > 0) {
444 const Scalar trans = problem_().transmissibility(myIdx, neighborIdx);
445 const auto scvfIdx = dofIdx - 1;
446 const auto& scvf = stencil.interiorFace(scvfIdx);
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);
451 const Scalar zEx = problem_().dofCenterDepth(neighborIdx);
452 const Scalar dZg = (zIn - zEx)*gravity;
453 const Scalar thpres = problem_().thresholdPressure(myIdx, neighborIdx);
454 Scalar inAlpha {0.};
455 Scalar outAlpha {0.};
456 Scalar diffusivity {0.};
457 Scalar dispersivity {0.};
458 if constexpr(enableEnergy){
459 inAlpha = problem_().thermalHalfTransmissibility(myIdx, neighborIdx);
460 outAlpha = problem_().thermalHalfTransmissibility(neighborIdx, myIdx);
461 }
462 if constexpr(enableDiffusion){
463 diffusivity = problem_().diffusivity(myIdx, neighborIdx);
464 }
465 if (simulator_().vanguard().eclState().getSimulationConfig().rock_config().dispersion()) {
466 dispersivity = problem_().dispersivity(myIdx, neighborIdx);
467 }
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};
472
473 }
474 }
475 neighborInfo_.appendRow(loc_nbinfo.begin(), loc_nbinfo.end());
476 if (problem_().nonTrivialBoundaryConditions()) {
477 for (unsigned bfIndex = 0; bfIndex < stencil.numBoundaryFaces(); ++bfIndex) {
478 const auto& bf = stencil.boundaryFace(bfIndex);
479 const int dir_id = bf.dirId();
480 // not for NNCs
481 if (dir_id < 0)
482 continue;
483 const auto [type, massrateAD] = problem_().boundaryCondition(myIdx, dir_id);
484 // Strip the unnecessary (and zero anyway) derivatives off massrate.
485 VectorBlock massrate(0.0);
486 for (size_t ii = 0; ii < massrate.size(); ++ii) {
487 massrate[ii] = massrateAD[ii].value();
488 }
489 const auto& exFluidState = problem_().boundaryFluidState(myIdx, dir_id);
490 BoundaryConditionData bcdata{type,
491 massrate,
492 exFluidState.pvtRegionIndex(),
493 bfIndex,
494 bf.area(),
495 bf.integrationPos()[dimWorld - 1],
496 exFluidState};
497 boundaryInfo_.push_back({myIdx, dir_id, bfIndex, bcdata});
498 }
499 }
500 }
501 }
502
503 // add the additional neighbors and degrees of freedom caused by the auxiliary
504 // equations
505 size_t numAuxMod = model.numAuxiliaryModules();
506 for (unsigned auxModIdx = 0; auxModIdx < numAuxMod; ++auxModIdx)
507 model.auxiliaryModule(auxModIdx)->addNeighbors(sparsityPattern);
508
509 // allocate raw matrix
510 jacobian_.reset(new SparseMatrixAdapter(simulator_()));
511 diagMatAddress_.resize(numCells);
512 // create matrix structure based on sparsity pattern
513 jacobian_->reserve(sparsityPattern);
514 for (unsigned globI = 0; globI < numCells; globI++) {
515 const auto& nbInfos = neighborInfo_[globI];
516 diagMatAddress_[globI] = jacobian_->blockAddress(globI, globI);
517 for (auto& nbInfo : nbInfos) {
518 nbInfo.matBlockAddress = jacobian_->blockAddress(nbInfo.neighbor, globI);
519 }
520 }
521
522 // Create dummy full domain.
523 fullDomain_.cells.resize(numCells);
524 std::iota(fullDomain_.cells.begin(), fullDomain_.cells.end(), 0);
525 }
526
527 // reset the global linear system of equations.
528 void resetSystem_()
529 {
530 residual_ = 0.0;
531 // zero all matrix entries
532 jacobian_->clear();
533 }
534
535 // Initialize the flows, flores, and velocity sparse tables
536 void createFlows_()
537 {
539 // If FLOWS/FLORES is set in any RPTRST in the schedule, then we initializate the sparse tables
540 // For now, do the same also if any block flows are requested (TODO: only save requested cells...)
541 // If DISPERC is in the deck, we initialize the sparse table here as well.
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) {
547 return;
548 }
549 const auto& model = model_();
550 const auto& nncOutput = simulator_().problem().eclWriter().getOutputNnc();
551 Stencil stencil(gridView_(), model_().dofMapper());
552 unsigned numCells = model.numTotalDof();
553 std::unordered_multimap<int, std::pair<int, int>> nncIndices;
554 std::vector<FlowInfo> loc_flinfo;
555 std::vector<VelocityInfo> loc_vlinfo;
556 unsigned int nncId = 0;
557 VectorBlock flow(0.0);
558
559 // Create a nnc structure to use fast lookup
560 for (unsigned int nncIdx = 0; nncIdx < nncOutput.size(); ++nncIdx) {
561 const int ci1 = nncOutput[nncIdx].cell1;
562 const int ci2 = nncOutput[nncIdx].cell2;
563 nncIndices.emplace(ci1, std::make_pair(ci2, nncIdx));
564 }
565
566 if (anyFlows) {
567 flowsInfo_.reserve(numCells, 6 * numCells);
568 }
569 if (anyFlores) {
570 floresInfo_.reserve(numCells, 6 * numCells);
571 }
572 if (enableDispersion) {
573 velocityInfo_.reserve(numCells, 6 * numCells);
574 }
575
576 for (const auto& elem : elements(gridView_())) {
577 stencil.update(elem);
578 for (unsigned primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++primaryDofIdx) {
579 unsigned myIdx = stencil.globalSpaceIndex(primaryDofIdx);
580 int numFaces = stencil.numBoundaryFaces() + stencil.numInteriorFaces();
581 loc_flinfo.resize(numFaces);
582 loc_vlinfo.resize(stencil.numDof() - 1);
583
584 for (unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
585 unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
586 if (dofIdx > 0) {
587 const auto scvfIdx = dofIdx - 1;
588 const auto& scvf = stencil.interiorFace(scvfIdx);
589 int faceId = scvf.dirId();
590 const int cartMyIdx = simulator_().vanguard().cartesianIndex(myIdx);
591 const int cartNeighborIdx = simulator_().vanguard().cartesianIndex(neighborIdx);
592 const auto& range = nncIndices.equal_range(cartMyIdx);
593 for (auto it = range.first; it != range.second; ++it) {
594 if (it->second.first == cartNeighborIdx){
595 // -1 gives problem since is used for the nncInput from the deck
596 faceId = -2;
597 // the index is stored to be used for writting the outputs
598 nncId = it->second.second;
599 }
600 }
601 loc_flinfo[dofIdx - 1] = FlowInfo{faceId, flow, nncId};
602 loc_vlinfo[dofIdx - 1] = VelocityInfo{flow};
603 }
604 }
605
606 for (unsigned bdfIdx = 0; bdfIdx < stencil.numBoundaryFaces(); ++bdfIdx) {
607 const auto& scvf = stencil.boundaryFace(bdfIdx);
608 int faceId = scvf.dirId();
609 loc_flinfo[stencil.numInteriorFaces() + bdfIdx] = FlowInfo{faceId, flow, nncId};
610 }
611
612
613 if (anyFlows) {
614 flowsInfo_.appendRow(loc_flinfo.begin(), loc_flinfo.end());
615 }
616 if (anyFlores) {
617 floresInfo_.appendRow(loc_flinfo.begin(), loc_flinfo.end());
618 }
619 if (enableDispersion) {
620 velocityInfo_.appendRow(loc_vlinfo.begin(), loc_vlinfo.end());
621 }
622 }
623 }
624 }
625
626public:
627 void setResAndJacobi(VectorBlock& res, MatrixBlock& bMat, const ADVectorBlock& resid) const
628 {
629 for (unsigned eqIdx = 0; eqIdx < numEq; eqIdx++)
630 res[eqIdx] = resid[eqIdx].value();
631
632 for (unsigned eqIdx = 0; eqIdx < numEq; eqIdx++) {
633 for (unsigned pvIdx = 0; pvIdx < numEq; pvIdx++) {
634 // A[dofIdx][focusDofIdx][eqIdx][pvIdx] is the partial derivative of
635 // the residual function 'eqIdx' for the degree of freedom 'dofIdx'
636 // with regard to the focus variable 'pvIdx' of the degree of freedom
637 // 'focusDofIdx'
638 bMat[eqIdx][pvIdx] = resid[eqIdx].derivative(pvIdx);
639 }
640 }
641 }
642
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();
648 if (!enableFlows && !enableFlores) {
649 return;
650 }
651 const unsigned int numCells = model_().numTotalDof();
652#ifdef _OPENMP
653#pragma omp parallel for
654#endif
655 for (unsigned globI = 0; globI < numCells; ++globI) {
657 const auto& nbInfos = neighborInfo_[globI];
658 ADVectorBlock adres(0.0);
660 const IntensiveQuantities& intQuantsIn = model_().intensiveQuantities(globI, /*timeIdx*/ 0);
661 // Flux term.
662 {
664 short loc = 0;
665 for (const auto& nbInfo : nbInfos) {
667 unsigned globJ = nbInfo.neighbor;
668 assert(globJ != globI);
669 adres = 0.0;
670 darcyFlux = 0.0;
671 const IntensiveQuantities& intQuantsEx = model_().intensiveQuantities(globJ, /*timeIdx*/ 0);
672 LocalResidual::computeFlux(adres,darcyFlux, globI, globJ, intQuantsIn, intQuantsEx, nbInfo.res_nbinfo, problem_().moduleParams());
673 adres *= nbInfo.res_nbinfo.faceArea;
674 if (enableFlows) {
675 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
676 flowsInfo_[globI][loc].flow[eqIdx] = adres[eqIdx].value();
677 }
678 }
679 if (enableFlores) {
680 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
681 floresInfo_[globI][loc].flow[eqIdx] = darcyFlux[eqIdx].value();
682 }
683 }
684 ++loc;
685 }
686 }
687 }
688
689 // Boundary terms. Only looping over cells with nontrivial bcs.
690 for (const auto& bdyInfo : boundaryInfo_) {
691 if (bdyInfo.bcdata.type == BCType::NONE)
692 continue;
693
694 ADVectorBlock adres(0.0);
695 const unsigned globI = bdyInfo.cell;
696 const auto& nbInfos = neighborInfo_[globI];
697 const IntensiveQuantities& insideIntQuants = model_().intensiveQuantities(globI, /*timeIdx*/ 0);
698 LocalResidual::computeBoundaryFlux(adres, problem_(), bdyInfo.bcdata, insideIntQuants, globI);
699 adres *= bdyInfo.bcdata.faceArea;
700 const unsigned bfIndex = bdyInfo.bfIndex;
701 if (enableFlows) {
702 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
703 flowsInfo_[globI][nbInfos.size() + bfIndex].flow[eqIdx] = adres[eqIdx].value();
704 }
705 }
706 // TODO also store Flores?
707 }
708 }
709
710private:
711 template <class SubDomainType>
712 void linearize_(const SubDomainType& domain)
713 {
714 // This check should be removed once this is addressed by
715 // for example storing the previous timesteps' values for
716 // rsmax (for DRSDT) and similar.
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.");
720 }
721 }
722
724
725 // We do not call resetSystem_() here, since that will set
726 // the full system to zero, not just our part.
727 // Instead, that must be called before starting the linearization.
728 const bool enableDispersion = simulator_().vanguard().eclState().getSimulationConfig().rock_config().dispersion();
729 const unsigned int numCells = domain.cells.size();
730 const bool on_full_domain = (numCells == model_().numTotalDof());
731
732#ifdef _OPENMP
733#pragma omp parallel for
734#endif
735 for (unsigned ii = 0; ii < numCells; ++ii) {
737 const unsigned globI = domain.cells[ii];
738 const auto& nbInfos = neighborInfo_[globI];
739 VectorBlock res(0.0);
740 MatrixBlock bMat(0.0);
741 ADVectorBlock adres(0.0);
743 const IntensiveQuantities& intQuantsIn = model_().intensiveQuantities(globI, /*timeIdx*/ 0);
744
745 // Flux term.
746 {
748 short loc = 0;
749 for (const auto& nbInfo : nbInfos) {
751 unsigned globJ = nbInfo.neighbor;
752 assert(globJ != globI);
753 res = 0.0;
754 bMat = 0.0;
755 adres = 0.0;
756 darcyFlux = 0.0;
757 const IntensiveQuantities& intQuantsEx = model_().intensiveQuantities(globJ, /*timeIdx*/ 0);
758 LocalResidual::computeFlux(adres,darcyFlux, globI, globJ, intQuantsIn, intQuantsEx, nbInfo.res_nbinfo, problem_().moduleParams());
759 adres *= nbInfo.res_nbinfo.faceArea;
760 if (enableDispersion) {
761 for (unsigned phaseIdx = 0; phaseIdx < numEq; ++ phaseIdx) {
762 velocityInfo_[globI][loc].velocity[phaseIdx] = darcyFlux[phaseIdx].value() / nbInfo.res_nbinfo.faceArea;
763 }
764 }
765 setResAndJacobi(res, bMat, adres);
766 residual_[globI] += res;
767 //SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat);
768 *diagMatAddress_[globI] += bMat;
769 bMat *= -1.0;
770 //SparseAdapter syntax: jacobian_->addToBlock(globJ, globI, bMat);
771 *nbInfo.matBlockAddress += bMat;
772 ++loc;
773 }
774 }
775
776 // Accumulation term.
777 double dt = simulator_().timeStepSize();
778 double volume = model_().dofTotalVolume(globI);
779 Scalar storefac = volume / dt;
780 adres = 0.0;
781 {
782 OPM_TIMEBLOCK_LOCAL(computeStorage);
783 LocalResidual::computeStorage(adres, intQuantsIn);
784 }
785 setResAndJacobi(res, bMat, adres);
786 // Either use cached storage term, or compute it on the fly.
787 if (model_().enableStorageCache()) {
788 // The cached storage for timeIdx 0 (current time) is not
789 // used, but after storage cache is shifted at the end of the
790 // timestep, it will become cached storage for timeIdx 1.
791 model_().updateCachedStorage(globI, /*timeIdx=*/0, res);
792 if (model_().newtonMethod().numIterations() == 0) {
793 // Need to update the storage cache.
794 if (problem_().recycleFirstIterationStorage()) {
795 // Assumes nothing have changed in the system which
796 // affects masses calculated from primary variables.
797 if (on_full_domain) {
798 // This is to avoid resetting the start-of-step storage
799 // to incorrect numbers when we do local solves, where the iteration
800 // number will start from 0, but the starting state may not be identical
801 // to the start-of-step state.
802 // Note that a full assembly must be done before local solves
803 // otherwise this will be left un-updated.
804 model_().updateCachedStorage(globI, /*timeIdx=*/1, res);
805 }
806 } else {
807 Dune::FieldVector<Scalar, numEq> tmp;
808 IntensiveQuantities intQuantOld = model_().intensiveQuantities(globI, 1);
809 LocalResidual::computeStorage(tmp, intQuantOld);
810 model_().updateCachedStorage(globI, /*timeIdx=*/1, tmp);
811 }
812 }
813 res -= model_().cachedStorage(globI, 1);
814 } else {
816 Dune::FieldVector<Scalar, numEq> tmp;
817 IntensiveQuantities intQuantOld = model_().intensiveQuantities(globI, 1);
818 LocalResidual::computeStorage(tmp, intQuantOld);
819 // assume volume do not change
820 res -= tmp;
821 }
822 res *= storefac;
823 bMat *= storefac;
824 residual_[globI] += res;
825 //SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat);
826 *diagMatAddress_[globI] += bMat;
827
828 // Cell-wise source terms.
829 // This will include well sources if SeparateSparseSourceTerms is false.
830 res = 0.0;
831 bMat = 0.0;
832 adres = 0.0;
833 if (separateSparseSourceTerms_) {
834 LocalResidual::computeSourceDense(adres, problem_(), globI, 0);
835 } else {
836 LocalResidual::computeSource(adres, problem_(), globI, 0);
837 }
838 adres *= -volume;
839 setResAndJacobi(res, bMat, adres);
840 residual_[globI] += res;
841 //SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat);
842 *diagMatAddress_[globI] += bMat;
843 } // end of loop for cell globI.
844
845 // Add sparse source terms. For now only wells.
846 if (separateSparseSourceTerms_) {
847 problem_().wellModel().addReservoirSourceTerms(residual_, diagMatAddress_);
848 }
849
850 // Boundary terms. Only looping over cells with nontrivial bcs.
851 for (const auto& bdyInfo : boundaryInfo_) {
852 if (bdyInfo.bcdata.type == BCType::NONE)
853 continue;
854
855 VectorBlock res(0.0);
856 MatrixBlock bMat(0.0);
857 ADVectorBlock adres(0.0);
858 const unsigned globI = bdyInfo.cell;
859 const IntensiveQuantities& insideIntQuants = model_().intensiveQuantities(globI, /*timeIdx*/ 0);
860 LocalResidual::computeBoundaryFlux(adres, problem_(), bdyInfo.bcdata, insideIntQuants, globI);
861 adres *= bdyInfo.bcdata.faceArea;
862 setResAndJacobi(res, bMat, adres);
863 residual_[globI] += res;
865 *diagMatAddress_[globI] += bMat;
866 }
867 }
868
869 void updateStoredTransmissibilities()
870 {
871 if (neighborInfo_.empty()) {
872 // This function was called before createMatrix_() was called.
873 // We call initFirstIteration_(), not createMatrix_(), because
874 // that will also initialize the residual consistently.
875 initFirstIteration_();
876 }
877 unsigned numCells = model_().numTotalDof();
878#ifdef _OPENMP
879#pragma omp parallel for
880#endif
881 for (unsigned globI = 0; globI < numCells; globI++) {
882 auto nbInfos = neighborInfo_[globI]; // nbInfos will be a SparseTable<...>::mutable_iterator_range.
883 for (auto& nbInfo : nbInfos) {
884 unsigned globJ = nbInfo.neighbor;
885 nbInfo.res_nbinfo.trans = problem_().transmissibility(globI, globJ);
886 }
887 }
888 }
889
890
891 Simulator *simulatorPtr_;
892
893 // the jacobian matrix
894 std::unique_ptr<SparseMatrixAdapter> jacobian_;
895
896 // the right-hand side
897 GlobalEqVector residual_;
898
899 LinearizationType linearizationType_;
900
901 using ResidualNBInfo = typename LocalResidual::ResidualNBInfo;
902 struct NeighborInfo
903 {
904 unsigned int neighbor;
905 ResidualNBInfo res_nbinfo;
906 MatrixBlock* matBlockAddress;
907 };
908 SparseTable<NeighborInfo> neighborInfo_;
909 std::vector<MatrixBlock*> diagMatAddress_;
910
911 struct FlowInfo
912 {
913 int faceId;
914 VectorBlock flow;
915 unsigned int nncId;
916 };
917 SparseTable<FlowInfo> flowsInfo_;
918 SparseTable<FlowInfo> floresInfo_;
919
920 struct VelocityInfo
921 {
922 VectorBlock velocity;
923 };
924 SparseTable<VelocityInfo> velocityInfo_;
925
926 using ScalarFluidState = typename IntensiveQuantities::ScalarFluidState;
927 struct BoundaryConditionData
928 {
929 BCType type;
930 VectorBlock massRate;
931 unsigned pvtRegionIdx;
932 unsigned boundaryFaceIndex;
933 double faceArea;
934 double faceZCoord;
935 ScalarFluidState exFluidState;
936 };
937 struct BoundaryInfo
938 {
939 unsigned int cell;
940 int dir;
941 unsigned int bfIndex;
942 BoundaryConditionData bcdata;
943 };
944 std::vector<BoundaryInfo> boundaryInfo_;
945 bool separateSparseSourceTerms_ = false;
946 struct FullDomain
947 {
948 std::vector<int> cells;
949 std::vector<bool> interior;
950 };
951 FullDomain fullDomain_;
952};
953
954} // namespace Opm
955
956#endif // TPFA_LINEARIZER
Base class for specifying auxiliary equations.
The base class for the element-centered finite-volume discretization scheme.
Definition ecfvdiscretization.hh:147
Definition matrixblock.hh:227
Manages the initializing and running of time dependent problems.
Definition simulator.hh:92
Scalar timeStepSize() const
Returns the time step length so that we don't miss the beginning of the next episode or cross the en...
Definition simulator.hh:453
Vanguard & vanguard()
Return a reference to the grid manager of simulation.
Definition simulator.hh:273
Problem & problem()
Return the object which specifies the pysical setup of the simulation.
Definition simulator.hh:304
const GridView & gridView() const
Return the grid view for which the simulation is done.
Definition simulator.hh:285
Model & model()
Return the physical model used in the simulation.
Definition simulator.hh:291
The common code for the linearizers of non-linear systems of equations.
Definition tpfalinearizer.hh:77
const auto & getFloresInfo() const
Return constant reference to the floresInfo.
Definition tpfalinearizer.hh:323
void linearize()
Linearize the full system of non-linear equations.
Definition tpfalinearizer.hh:173
const auto & getFlowsInfo() const
Return constant reference to the flowsInfo.
Definition tpfalinearizer.hh:313
void linearizeDomain(const SubDomainType &domain)
Linearize the part of the non-linear system of equations that is associated with a part of the spatia...
Definition tpfalinearizer.hh:227
void init(Simulator &simulator)
Initialize the linearizer.
Definition tpfalinearizer.hh:146
static void registerParameters()
Register all run-time parameters for the Jacobian linearizer.
Definition tpfalinearizer.hh:131
void linearizeDomain()
Linearize the part of the non-linear system of equations that is associated with the spatial domain.
Definition tpfalinearizer.hh:189
const SparseMatrixAdapter & jacobian() const
Return constant reference to global Jacobian matrix backend.
Definition tpfalinearizer.hh:285
void eraseMatrix()
Causes the Jacobian matrix to be recreated from scratch before the next iteration.
Definition tpfalinearizer.hh:159
const auto & getVelocityInfo() const
Return constant reference to the velocityInfo.
Definition tpfalinearizer.hh:333
const std::map< unsigned, Constraints > constraintsMap() const
Returns the map of constraint degrees of freedom.
Definition tpfalinearizer.hh:366
void linearizeAuxiliaryEquations()
Linearize the part of the non-linear system of equations that is associated with the spatial domain.
Definition tpfalinearizer.hh:254
const GlobalEqVector & residual() const
Return constant reference to global residual vector.
Definition tpfalinearizer.hh:294
Declare the properties used by the infrastructure code of the finite volume discretizations.
The common code for the linearizers of non-linear systems of equations.
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition propertysystem.hh:235
Definition tpfalinearizer.hh:56