Add some more high-level tracy hooks.

Mostly for looking at TemperatureModel performance.
No hooks in inner loops (cell or connection loops), so runtime with tracy
active should be very close to without.

opm/simulators/flow/GenericOutputBlackoilModule.cpp

@@ -1050,6 +1050,7 @@ template<class FluidSystem>
 Inplace GenericOutputBlackoilModule<FluidSystem>::
 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<std::string, double>& miscSummaryData,
                           std::map<std::string, std::vector<double>>& regionData) const
 {
+    OPM_TIMEBLOCK(GenericOutputBlackoilModule_updateSummaryRegionValues);
     // The field summary vectors should only use the FIPNUM based region sum.
     {
         for (const auto& phase : Inplace::phases()) {

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<Parameters::RelaxedMaxPvFraction<Scalar>>();
         const bool relax = (sum_pv_not_converged / sum_pv) <  relaxed_max_pv_fraction;
         const auto tolerance_energy_balance = relax? Parameters::Get<Parameters::ToleranceEnergyBalanceRelaxed<Scalar>>():
@@ -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<Parameters::EnableDriftCompensationTemp>();
-        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<Parameters::EnableDriftCompensationTemp>();
+            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"

opm/simulators/linalg/ParallelOverlappingILU0_impl.hpp

@@ -42,6 +42,7 @@ namespace detail
 template<class M>
 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<std::size_t> 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<Matrix>(*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