diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.cpp b/opm/simulators/flow/GenericOutputBlackoilModule.cpp index 272130b03b..8757467007 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.cpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.cpp @@ -1050,6 +1050,7 @@ template Inplace GenericOutputBlackoilModule:: accumulateRegionSums(const Parallel::Communication& comm) { + OPM_TIMEBLOCK(GenericOutputBlackoilModule_accumulateRegionSums); Inplace inplace; for (const auto& region : this->regions_) { @@ -1073,6 +1074,7 @@ updateSummaryRegionValues(const Inplace& inplace, std::map& miscSummaryData, std::map>& regionData) const { + OPM_TIMEBLOCK(GenericOutputBlackoilModule_updateSummaryRegionValues); // The field summary vectors should only use the FIPNUM based region sum. { for (const auto& phase : Inplace::phases()) { diff --git a/opm/simulators/flow/TemperatureModel.hpp b/opm/simulators/flow/TemperatureModel.hpp index 704a7cd5c0..2b5a469f34 100644 --- a/opm/simulators/flow/TemperatureModel.hpp +++ b/opm/simulators/flow/TemperatureModel.hpp @@ -333,6 +333,7 @@ public: if (!this->doTemp()) { return; } + OPM_TIMEBLOCK(TemperatureModel_endTimeStep); // We use the specialized intensive quantities here with only the temperature derivative const unsigned int numCells = simulator_.model().numTotalDof(); @@ -396,6 +397,7 @@ protected: void advanceTemperatureFields() { + OPM_TIMEBLOCK(TemperatureModel_advanceTemperatureFields); const int max_iter = 20; const int min_iter = 1; bool is_converged = false; @@ -419,6 +421,7 @@ protected: void solveAndUpdate() { + OPM_TIMEBLOCK(TemperatureModel_solveAndUpdate); const unsigned int numCells = simulator_.model().numTotalDof(); EnergyVector dx(numCells); bool conv = this->linearSolve_(this->energyMatrix_->istlMatrix(), dx, this->energyVector_); @@ -428,6 +431,7 @@ protected: } } else { + OPM_TIMEBLOCK(TemperatureModel_solveAndUpdate_update); #ifdef _OPENMP #pragma omp parallel for #endif @@ -441,6 +445,7 @@ protected: bool converged(const int iter) { + OPM_TIMEBLOCK(TemperatureModel_converged); Scalar dt = simulator_.timeStepSize(); Scalar maxNorm = 0.0; Scalar sumNorm = 0.0; @@ -464,13 +469,17 @@ protected: sum_pv += pvValue; } } - maxNorm = simulator_.gridView().comm().max(maxNorm); - sumNorm = simulator_.gridView().comm().sum(sumNorm); - sum_pv = simulator_.gridView().comm().sum(sum_pv); - sumNorm /= sum_pv; + { + OPM_TIMEBLOCK(TemperatureModel_converged_communicate); - // Use relaxed tolerance if the fraction of unconverged cells porevolume is less than relaxed_max_pv_fraction - sum_pv_not_converged = simulator_.gridView().comm().sum(sum_pv_not_converged); + maxNorm = simulator_.gridView().comm().max(maxNorm); + sumNorm = simulator_.gridView().comm().sum(sumNorm); + sum_pv = simulator_.gridView().comm().sum(sum_pv); + sumNorm /= sum_pv; + + // Use relaxed tolerance if the fraction of unconverged cells porevolume is less than relaxed_max_pv_fraction + sum_pv_not_converged = simulator_.gridView().comm().sum(sum_pv_not_converged); + } Scalar relaxed_max_pv_fraction = Parameters::Get>(); const bool relax = (sum_pv_not_converged / sum_pv) < relaxed_max_pv_fraction; const auto tolerance_energy_balance = relax? Parameters::Get>(): @@ -571,76 +580,90 @@ protected: void assembleEquations() { + OPM_TIMEBLOCK(TemperatureModel_assembleEquations); + const unsigned int numCells = simulator_.model().numTotalDof(); for (unsigned globI = 0; globI < numCells; ++globI) { this->energyVector_[globI] = 0.0; energyMatrix_->clearRow(globI, 0.0); } - Scalar dt = simulator_.timeStepSize(); - #ifdef _OPENMP - #pragma omp parallel for - #endif - for (unsigned globI = 0; globI < numCells; ++globI) { - MatrixBlockTemp bMat; - Scalar volume = simulator_.model().dofTotalVolume(globI); - Scalar storefac = volume / dt; - Evaluation storage = 0.0; - computeStorageTerm(globI, storage); - this->energyVector_[globI] += storefac * ( getValue(storage) - storage1_[globI] ); - bMat[0][0] = storefac * storage.derivative(temperatureIdx); - *diagMatAddress_[globI] += bMat; + const Scalar dt = simulator_.timeStepSize(); + + // Storage term + { + OPM_TIMEBLOCK(TemperatureModel_assembleEquations_storage); + #ifdef _OPENMP + #pragma omp parallel for + #endif + for (unsigned globI = 0; globI < numCells; ++globI) { + MatrixBlockTemp bMat; + Scalar volume = simulator_.model().dofTotalVolume(globI); + Scalar storefac = volume / dt; + Evaluation storage = 0.0; + computeStorageTerm(globI, storage); + this->energyVector_[globI] += storefac * ( getValue(storage) - storage1_[globI] ); + bMat[0][0] = storefac * storage.derivative(temperatureIdx); + *diagMatAddress_[globI] += bMat; + } } - const auto& floresInfo = this->simulator_.problem().model().linearizer().getFloresInfo(); - const bool enableDriftCompensation = Parameters::Get(); - const auto& problem = simulator_.problem(); + // Flux term + { + OPM_TIMEBLOCK(TemperatureModel_assembleEquations_flux); + const auto& floresInfo = this->simulator_.problem().model().linearizer().getFloresInfo(); + const bool enableDriftCompensation = Parameters::Get(); + const auto& problem = simulator_.problem(); - #ifdef _OPENMP - #pragma omp parallel for - #endif - for (unsigned globI = 0; globI < numCells; ++globI) { - const auto& nbInfos = neighborInfo_[globI]; - const auto& floresInfos = floresInfo[globI]; - int loc = 0; - const auto& intQuantsIn = intQuants_[globI]; - MatrixBlockTemp bMat; - for (const auto& nbInfo : nbInfos) { - unsigned globJ = nbInfo.neighbor; - const auto& intQuantsEx = intQuants_[globJ]; - assert(globJ != globI); - const auto& darcyflux = floresInfos[loc].flow; - // compute convective flux - Evaluation flux = 0.0; - computeFluxTerm(intQuantsIn.fluidStateTemp(), intQuantsEx.fluidStateTemp(), darcyflux, flux); - // compute conductive flux - Evaluation heatFlux = 0.0; - computeHeatFluxTerm(intQuantsIn, intQuantsEx, nbInfo.res_nbinfo, heatFlux); - heatFlux += flux; - this->energyVector_[globI] += getValue(heatFlux); - bMat[0][0] = heatFlux.derivative(temperatureIdx); - *diagMatAddress_[globI] += bMat; - bMat *= -1.0; - //SparseAdapter syntax: jacobian_->addToBlock(globJ, globI, bMat); - *nbInfo.matBlockAddress += bMat; - loc++; - } + #ifdef _OPENMP + #pragma omp parallel for + #endif + for (unsigned globI = 0; globI < numCells; ++globI) { + const auto& nbInfos = neighborInfo_[globI]; + const auto& floresInfos = floresInfo[globI]; + int loc = 0; + const auto& intQuantsIn = intQuants_[globI]; + MatrixBlockTemp bMat; + for (const auto& nbInfo : nbInfos) { + unsigned globJ = nbInfo.neighbor; + const auto& intQuantsEx = intQuants_[globJ]; + assert(globJ != globI); + const auto& darcyflux = floresInfos[loc].flow; + // compute convective flux + Evaluation flux = 0.0; + computeFluxTerm(intQuantsIn.fluidStateTemp(), intQuantsEx.fluidStateTemp(), darcyflux, flux); + // compute conductive flux + Evaluation heatFlux = 0.0; + computeHeatFluxTerm(intQuantsIn, intQuantsEx, nbInfo.res_nbinfo, heatFlux); + heatFlux += flux; + this->energyVector_[globI] += getValue(heatFlux); + bMat[0][0] = heatFlux.derivative(temperatureIdx); + *diagMatAddress_[globI] += bMat; + bMat *= -1.0; + //SparseAdapter syntax: jacobian_->addToBlock(globJ, globI, bMat); + *nbInfo.matBlockAddress += bMat; + loc++; + } - if (enableDriftCompensation) { - auto dofDriftRate = problem.drift()[globI]/dt; - const auto& fs = intQuantsIn.fluidStateTemp(); - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { - const unsigned activeCompIdx = - FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx)); - auto drift_hrate = dofDriftRate[activeCompIdx]*getValue(fs.enthalpy(phaseIdx)) * getValue(fs.density(phaseIdx)) / getValue(fs.invB(phaseIdx)); - this->energyVector_[globI] -= drift_hrate*scalingFactor_; + if (enableDriftCompensation) { + auto dofDriftRate = problem.drift()[globI]/dt; + const auto& fs = intQuantsIn.fluidStateTemp(); + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + const unsigned activeCompIdx = + FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx)); + auto drift_hrate = dofDriftRate[activeCompIdx]*getValue(fs.enthalpy(phaseIdx)) * getValue(fs.density(phaseIdx)) / getValue(fs.invB(phaseIdx)); + this->energyVector_[globI] -= drift_hrate*scalingFactor_; + } } } } // Well terms - const auto& wellPtrs = simulator_.problem().wellModel().localNonshutWells(); - for (const auto& wellPtr : wellPtrs) { - this->assembleEquationWell(*wellPtr); + { + OPM_TIMEBLOCK(TemperatureModel_assembleEquations_wells); + const auto& wellPtrs = simulator_.problem().wellModel().localNonshutWells(); + for (const auto& wellPtr : wellPtrs) { + this->assembleEquationWell(*wellPtr); + } } // For standard ilu0 we need to set the overlapping cells to identity. But since we use "paroverilu0" diff --git a/opm/simulators/linalg/ParallelOverlappingILU0_impl.hpp b/opm/simulators/linalg/ParallelOverlappingILU0_impl.hpp index 6c86be406c..4e1f0fa779 100644 --- a/opm/simulators/linalg/ParallelOverlappingILU0_impl.hpp +++ b/opm/simulators/linalg/ParallelOverlappingILU0_impl.hpp @@ -42,6 +42,7 @@ namespace detail template void ghost_last_bilu0_decomposition (M& A, std::size_t interiorSize) { + OPM_TIMEBLOCK(GhostLastBlockILU0Decomp); // iterator types assert(interiorSize <= A.N()); using rowiterator = typename M::RowIterator; @@ -434,10 +435,13 @@ update() } std::vector inverseOrdering(ordering_.size()); - std::size_t index = 0; - for (const auto newIndex : ordering_) { - inverseOrdering[newIndex] = index++; + OPM_TIMEBLOCK(createInverseOrdering); + std::size_t index = 0; + for (const auto newIndex : ordering_) + { + inverseOrdering[newIndex] = index++; + } } try @@ -453,8 +457,9 @@ update() // create ILU-0 decomposition if (ordering_.empty()) { + OPM_TIMEBLOCK(iluDecompositionUpdateMatrix); if (ILU_) { - OPM_TIMEBLOCK(iluDecompositionMakeMatrix); + OPM_TIMEBLOCK(iluDecompositionCopyEntries); // The ILU_ matrix is already a copy with the same // sparse structure as A_, but the values of A_ may // have changed, so we must copy all elements. @@ -468,6 +473,7 @@ update() } } } else { + OPM_TIMEBLOCK(iluDecompositionDuplicateMatrix); // First call, must duplicate matrix. ILU_ = std::make_unique(*A_); } @@ -550,11 +556,14 @@ update() } // Check whether there was a problem on some process - const bool parallel_failure = comm_ && comm_->communicator().min(ilu_setup_successful) == 0; - const bool local_failure = ilu_setup_successful == 0; - if (local_failure || parallel_failure) { - throw Dune::MatrixBlockError(); + OPM_TIMEBLOCK(checkParallelFailure); + const bool parallel_failure = comm_ && comm_->communicator().min(ilu_setup_successful) == 0; + const bool local_failure = ilu_setup_successful == 0; + if (local_failure || parallel_failure) + { + throw Dune::MatrixBlockError(); + } } // store ILU in simple CRS format