tpfalinearizer.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
  This file is part of the Open Porous Media project (OPM).

  OPM is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 2 of the License, or
  (at your option) any later version.

  OPM is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with OPM.  If not, see <http://www.gnu.org/licenses/>.

  Consult the COPYING file in the top-level source directory of this
  module for the precise wording of the license and the list of
  copyright holders.
*/
#ifndef TPFA_LINEARIZER_HH
#define TPFA_LINEARIZER_HH

#include <dune/common/version.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>

#include <opm/common/Exceptions.hpp>
#include <opm/common/TimingMacros.hpp>

#include <opm/grid/utility/SparseTable.hpp>

#include <opm/material/common/ConditionalStorage.hpp>

#include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>
#include <opm/input/eclipse/Schedule/BCProp.hpp>

#include <opm/models/blackoil/blackoilproperties.hh>
#include <opm/models/common/multiphasebaseproperties.hh>
#include <opm/models/discretization/common/baseauxiliarymodule.hh>
#include <opm/models/discretization/common/fvbaseproperties.hh>
#include <opm/models/discretization/common/linearizationtype.hh>
#include <opm/models/discretization/common/tpfalinearizerstructs.hh>

#include <opm/simulators/linalg/exportSystem.hpp>

#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
#include <opm/material/fluidsystems/BlackOilFluidSystemNonStatic.hpp>

#include <cassert>
#include <cstddef>
#include <exception>   // current_exception, rethrow_exception
#include <iostream>
#include <map>
#include <memory>
#include <numeric>
#include <set>
#include <stdexcept>
#include <unordered_map>
#include <utility>
#include <vector>

#if HAVE_OPENMP
#include <omp.h>
#endif

#include <fmt/format.h>

#include <opm/common/utility/gpuDecorators.hpp>
#include <opm/common/utility/gpuistl_if_available.hpp>
#include <opm/common/utility/pointerArithmetic.hpp>
#if HAVE_CUDA
#include <opm/models/discretization/common/tpfalinearizergpukernels.hh>
#include <opm/simulators/flow/SimpleFIBlackOilModel.hpp>
#include <opm/simulators/flow/ThermalGasWaterFlowProblem.hpp>
#include <opm/models/discretization/common/tpfalinearizergpuparams.hh>
#if USE_HIP
#include <opm/simulators/linalg/gpuistl_hip/GpuSparseMatrixWrapper.hpp>
#include <opm/simulators/linalg/gpuistl_hip/MiniMatrix.hpp>
#include <opm/simulators/linalg/gpuistl_hip/MiniVector.hpp>
#include <opm/simulators/linalg/gpuistl_hip/detail/gpusparse_matrix_operations.hpp>
#else
#include <opm/simulators/linalg/gpuistl/GpuSparseMatrixWrapper.hpp>
#include <opm/simulators/linalg/gpuistl/MiniMatrix.hpp>
#include <opm/simulators/linalg/gpuistl/MiniVector.hpp>
#include <opm/simulators/linalg/gpuistl/detail/gpusparse_matrix_operations.hpp>
#endif
#endif

namespace Opm::Parameters {

struct SeparateSparseSourceTerms { static constexpr bool value = false; };

} // namespace Opm::Parameters

namespace Opm {

// forward declarations
template<class TypeTag>
class EcfvDiscretization;

template<class TypeTag>
class TpfaLinearizer
{
    using Model = GetPropType<TypeTag, Properties::Model>;
    using Problem = GetPropType<TypeTag, Properties::Problem>;
    using Simulator = GetPropType<TypeTag, Properties::Simulator>;
    using GridView = GetPropType<TypeTag, Properties::GridView>;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;

    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
    using GlobalEqVector = GetPropType<TypeTag, Properties::GlobalEqVector>;
    using SparseMatrixAdapter = GetPropType<TypeTag, Properties::SparseMatrixAdapter>;
    using EqVector = GetPropType<TypeTag, Properties::EqVector>;
    using Constraints = GetPropType<TypeTag, Properties::Constraints>;
    using Stencil = GetPropType<TypeTag, Properties::Stencil>;
    using LocalResidual = GetPropType<TypeTag, Properties::LocalResidual>;
    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
    using Indices = GetPropType<TypeTag, Properties::Indices>;

    using Element = typename GridView::template Codim<0>::Entity;
    using ElementIterator = typename GridView::template Codim<0>::Iterator;

    using Vector = GlobalEqVector;

    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;

    enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
    enum { historySize = getPropValue<TypeTag, Properties::TimeDiscHistorySize>() };
    enum { dimWorld = GridView::dimensionworld };

    constexpr static bool runAssemblyOnGpu = getPropValue<TypeTag, Properties::RunAssemblyOnGpu>();

    using MatrixBlockCPU = typename SparseMatrixAdapter::MatrixBlock;
    using VectorBlockCPU = Dune::FieldVector<Scalar, numEq>;
    using ADVectorBlockCPU = GetPropType<TypeTag, Properties::RateVector>;

#if HAVE_CUDA
    using MatrixBlockGPU = gpuistl::MiniMatrix<Scalar, numEq>;
    using VectorBlockGPU = gpuistl::MiniVector<Scalar, numEq>;
    using ADVectorBlockGPU = gpuistl::MiniVector<Evaluation, numEq>;
    using MatrixBlock = std::conditional_t<runAssemblyOnGpu, MatrixBlockGPU, MatrixBlockCPU>;
    using VectorBlock = std::conditional_t<runAssemblyOnGpu, VectorBlockGPU, VectorBlockCPU>;
    using ADVectorBlock = std::conditional_t<runAssemblyOnGpu, ADVectorBlockGPU, ADVectorBlockCPU>;
#else
    using MatrixBlock = MatrixBlockCPU;
    using VectorBlock = VectorBlockCPU;
    using ADVectorBlock = ADVectorBlockCPU;
#endif

    static constexpr bool linearizeNonLocalElements =
        getPropValue<TypeTag, Properties::LinearizeNonLocalElements>();
    static constexpr bool enableFullyImplicitThermal = (getPropValue<TypeTag, Properties::EnergyModuleType>() == EnergyModules::FullyImplicitThermal);
    static constexpr bool enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>();
    static constexpr bool enableDispersion = getPropValue<TypeTag, Properties::EnableDispersion>();
    static const bool enableBioeffects = getPropValue<TypeTag, Properties::EnableBioeffects>();

    // copying the linearizer is not a good idea
    TpfaLinearizer(const TpfaLinearizer&) = delete;

public:
    TpfaLinearizer()
    {
        simulatorPtr_ = nullptr;
        separateSparseSourceTerms_ = Parameters::Get<Parameters::SeparateSparseSourceTerms>();
        exportIndex_=-1;
        exportCount_=-1;
    }

    static void registerParameters()
    {
        Parameters::Register<Parameters::SeparateSparseSourceTerms>
            ("Treat well source terms all in one go, instead of on a cell by cell basis.");
    }

    void init(Simulator& simulator)
    {
        simulatorPtr_ = &simulator;
        eraseMatrix();
    }

    void eraseMatrix()
    {
        jacobian_.reset();
    }

    void linearize()
    {
        linearizeDomain();
        linearizeAuxiliaryEquations();
    }

    void linearizeDomain()
    {
        int succeeded;
        try {
            linearizeDomain(fullDomain_);
            succeeded = 1;
        }
        catch (const std::exception& e) {
            std::cout << "rank " << simulator_().gridView().comm().rank()
                      << " caught an exception while linearizing:" << e.what()
                      << "\n"  << std::flush;
            succeeded = 0;
        }
        catch (...) {
            std::cout << "rank " << simulator_().gridView().comm().rank()
                      << " caught an exception while linearizing"
                      << "\n"  << std::flush;
            succeeded = 0;
        }
        OPM_TIMEBLOCK(linearizationSynch);
        succeeded = simulator_().gridView().comm().min(succeeded);

        if (!succeeded) {
            throw NumericalProblem("A process did not succeed in linearizing the system");
        }
    }

    template <class SubDomainType>
    void linearizeDomain(const SubDomainType& domain)
    {
        OPM_TIMEBLOCK(linearizeDomain);
        // we defer the initialization of the Jacobian matrix until here because the
        // auxiliary modules usually assume the problem, model and grid to be fully
        // initialized...
        if (!jacobian_) {
            initFirstIteration_();
        }

        // Called here because it is no longer called from linearize_().
        if (problem_().iterationContext().inLocalSolve()) {
            resetSystem_(domain);
        }
        else {
            resetSystem_();
        }

        linearize_(domain);
    }

    void finalize()
    { jacobian_->finalize(); }

    void linearizeAuxiliaryEquations()
    {
        OPM_TIMEBLOCK(linearizeAuxilaryEquations);
        // flush possible local caches into matrix structure
        jacobian_->commit();

        auto& model = model_();
        const auto& comm = simulator_().gridView().comm();
        for (unsigned auxModIdx = 0; auxModIdx < model.numAuxiliaryModules(); ++auxModIdx) {
            bool succeeded = true;
            try {
                model.auxiliaryModule(auxModIdx)->linearize(*jacobian_, residual_);
            }
            catch (const std::exception& e) {
                succeeded = false;

                std::cout << "rank " << simulator_().gridView().comm().rank()
                          << " caught an exception while linearizing:" << e.what()
                          << "\n"  << std::flush;
            }

            succeeded = comm.min(succeeded);

            if (!succeeded) {
                throw NumericalProblem("linearization of an auxiliary equation failed");
            }
        }
    }

    const SparseMatrixAdapter& jacobian() const
    { return *jacobian_; }

    SparseMatrixAdapter& jacobian()
    { return *jacobian_; }

    const GlobalEqVector& residual() const
    { return residual_; }

    GlobalEqVector& residual()
    { return residual_; }

    void exportSystem(const int idx, std::string& tag, const char *path="export")
    {
        const bool export_sparsity = exportIndex_ == -1;

        // increment indices and generate tag
        exportCount_ = exportIndex_ == idx ? ++exportCount_ : 0;
        exportIndex_ = idx;
        tag = fmt::format(fmt::runtime("_{:03d}_{:02d}"), exportIndex_, exportCount_);

        fmt::print(fmt::runtime("index = {:d}\n"), exportIndex_);
        fmt::print(fmt::runtime("count = {:d}\n"), exportCount_);

        Opm::exportSystem(jacobian_->istlMatrix(), residual_, export_sparsity, tag.c_str(), path);
    }

    void setLinearizationType(LinearizationType linearizationType)
    { linearizationType_ = linearizationType; }

    const LinearizationType& getLinearizationType() const
    { return linearizationType_; }

    const auto& getFlowsInfo() const
    { return flowsInfo_; }

    const auto& getFloresInfo() const
    { return floresInfo_; }

    const auto& getVelocityInfo() const
    { return velocityInfo_; }

    const auto& getNeighborInfo() const {
        return neighborInfo_;
    }


    void updateDiscretizationParameters()
    {
        updateStoredTransmissibilities();
    }

    void updateBoundaryConditionData()
    {
        for (auto& bdyInfo : boundaryInfo_) {
            const auto [type, massrateAD] = problem_().boundaryCondition(bdyInfo.cell, bdyInfo.dir);

            // Strip the unnecessary (and zero anyway) derivatives off massrate.
            VectorBlockCPU massrate(0.0);
            for (std::size_t ii = 0; ii < massrate.size(); ++ii) {
                massrate[ii] = massrateAD[ii].value();
            }
            if (type != BCType::NONE) {
                const auto& exFluidState = problem_().boundaryFluidState(bdyInfo.cell, bdyInfo.dir);
                bdyInfo.bcdata.type = type;
                bdyInfo.bcdata.massRate = massrate;
                bdyInfo.bcdata.exFluidState = exFluidState;
            }
        }
    }

    std::map<unsigned, Constraints> constraintsMap() const
    { return {}; }

    template <class SubDomainType>
    void resetSystem_(const SubDomainType& domain)
    {
        if (!jacobian_) {
            initFirstIteration_();
        }
        for (int globI : domain.cells) {
            residual_[globI] = 0.0;
            jacobian_->clearRow(globI, 0.0);
        }
    }

private:
    Simulator& simulator_()
    { return *simulatorPtr_; }

    const Simulator& simulator_() const
    { return *simulatorPtr_; }

    Problem& problem_()
    { return simulator_().problem(); }

    const Problem& problem_() const
    { return simulator_().problem(); }

    Model& model_()
    { return simulator_().model(); }

    const Model& model_() const
    { return simulator_().model(); }

    const GridView& gridView_() const
    { return problem_().gridView(); }

    void initFirstIteration_()
    {
        // initialize the BCRS matrix for the Jacobian of the residual function
        createMatrix_();

        // initialize the Jacobian matrix and the vector for the residual function
        residual_.resize(model_().numTotalDof());
        resetSystem_();

        // initialize the sparse tables for Flows and Flores
        createFlows_();
    }

    // Construct the BCRS matrix for the Jacobian of the residual function
    void createMatrix_()
    {
        OPM_TIMEBLOCK(createMatrix);
        if (!neighborInfo_.empty()) {
            // It is ok to call this function multiple times, but it
            // should not do anything if already called.
            return;
        }
        const auto& model = model_();
        Stencil stencil(gridView_(), model_().dofMapper());

        // for the main model, find out the global indices of the neighboring degrees of
        // freedom of each primary degree of freedom
        using NeighborSet = std::set<unsigned>;
        std::vector<NeighborSet> sparsityPattern(model.numTotalDof());
        const Scalar gravity = problem_().gravity()[dimWorld - 1];
        unsigned numCells = model.numTotalDof();
        neighborInfo_.reserve(numCells, 6 * numCells); // Expect ~6 neighbors per cell
        std::vector<NeighborInfoCPU> loc_nbinfo;
        for (const auto& elem : elements(gridView_())) {
            stencil.update(elem);

            for (unsigned primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++primaryDofIdx) {
                const unsigned myIdx = stencil.globalSpaceIndex(primaryDofIdx);
                loc_nbinfo.resize(stencil.numDof() - 1); // Do not include the primary dof in neighborInfo_

                for (unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
                    const unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
                    sparsityPattern[myIdx].insert(neighborIdx);
                    if (dofIdx > 0) {
                        const Scalar trans = problem_().transmissibility(myIdx, neighborIdx);
                        const auto scvfIdx = dofIdx - 1;
                        const auto& scvf = stencil.interiorFace(scvfIdx);
                        const Scalar area = scvf.area();
                        const Scalar Vin = problem_().model().dofTotalVolume(myIdx);
                        const Scalar Vex = problem_().model().dofTotalVolume(neighborIdx);
                        const Scalar zIn = problem_().dofCenterDepth(myIdx);
                        const Scalar zEx = problem_().dofCenterDepth(neighborIdx);
                        const Scalar dZg = (zIn - zEx)*gravity;
                        const Scalar thpresInToEx = problem_().thresholdPressure(myIdx, neighborIdx);
                        const Scalar thpresExToIn = problem_().thresholdPressure(neighborIdx, myIdx);
                        const auto dirId = scvf.dirId();
                        auto faceDir = dirId < 0 ? FaceDir::DirEnum::Unknown
                                                 : FaceDir::FromIntersectionIndex(dirId);
                        ResidualNBInfo nbinfo{trans,
                                              area,
                                              thpresInToEx,
                                              thpresExToIn,
                                              dZg,
                                              faceDir,
                                              Vin,
                                              Vex,
                                              {},
                                              {},
                                              {},
                                              {}};
                        if constexpr (enableFullyImplicitThermal) {
                            nbinfo.inAlpha = problem_().thermalHalfTransmissibility(myIdx, neighborIdx);
                            nbinfo.outAlpha = problem_().thermalHalfTransmissibility(neighborIdx, myIdx);
                        }
                        if constexpr (enableDiffusion) {
                            nbinfo.diffusivity = problem_().diffusivity(myIdx, neighborIdx);
                        }
                        if constexpr (enableDispersion) {
                            nbinfo.dispersivity = problem_().dispersivity(myIdx, neighborIdx);
                        }
                        loc_nbinfo[dofIdx - 1] = NeighborInfoCPU{neighborIdx, nbinfo, nullptr};
                    }
                }
                neighborInfo_.appendRow(loc_nbinfo.begin(), loc_nbinfo.end());
                if (problem_().nonTrivialBoundaryConditions()) {
                    for (unsigned bfIndex = 0; bfIndex < stencil.numBoundaryFaces(); ++bfIndex) {
                        const auto& bf = stencil.boundaryFace(bfIndex);
                        const int dir_id = bf.dirId();
                        // not for NNCs
                        if (dir_id < 0) {
                            continue;
                        }
                        const auto [type, massrateAD] = problem_().boundaryCondition(myIdx, dir_id);
                        // Strip the unnecessary (and zero anyway) derivatives off massrate.
                        VectorBlockCPU massrate(0.0);
                        for (std::size_t ii = 0; ii < massrate.size(); ++ii) {
                            massrate[ii] = massrateAD[ii].value();
                        }
                        const auto& exFluidState = problem_().boundaryFluidState(myIdx, dir_id);
                        BoundaryConditionDataCPU bcdata {type,
                                                         massrate,
                                                         exFluidState.pvtRegionIndex(),
                                                         bfIndex,
                                                         bf.area(),
                                                         bf.integrationPos()[dimWorld - 1],
                                                         exFluidState};
                        boundaryInfo_.push_back({myIdx, dir_id, bfIndex, bcdata});
                    }
                }
            }
        }

        // add the additional neighbors and degrees of freedom caused by the auxiliary
        // equations
        const std::size_t numAuxMod = model.numAuxiliaryModules();
        for (unsigned auxModIdx = 0; auxModIdx < numAuxMod; ++auxModIdx) {
            model.auxiliaryModule(auxModIdx)->addNeighbors(sparsityPattern);
        }

        // allocate raw matrix
        jacobian_ = std::make_unique<SparseMatrixAdapter>(simulator_());
        diagMatAddress_.resize(numCells);
        // create matrix structure based on sparsity pattern
        jacobian_->reserve(sparsityPattern);
        for (unsigned globI = 0; globI < numCells; globI++) {
            const auto& nbInfos = neighborInfo_[globI];
            diagMatAddress_[globI] = jacobian_->blockAddress(globI, globI);
            for (auto& nbInfo : nbInfos) {
                nbInfo.matBlockAddress = jacobian_->blockAddress(nbInfo.neighbor, globI);
            }
        }

#if HAVE_CUDA && OPM_IS_COMPILING_WITH_GPU_COMPILER
        gpuJacobian_.reset(new gpuistl::GpuSparseMatrixWrapper<Scalar>(
            gpuistl::GpuSparseMatrixWrapper<Scalar>::fromMatrix(jacobian_->istlMatrix())));
        gpuBufferDiagMatAddress_.reset(new gpuistl::GpuBuffer<MatrixBlockGPU*>(
            gpuistl::detail::getDiagPtrsTyped<MatrixBlockGPU>(*gpuJacobian_)));
#endif

        // Create dummy full domain.
        fullDomain_.cells.resize(numCells);
        std::iota(fullDomain_.cells.begin(), fullDomain_.cells.end(), 0);
    }

    // reset the global linear system of equations.
    void resetSystem_()
    {
        residual_ = 0.0;
        // zero all matrix entries
        jacobian_->clear();

#if HAVE_CUDA && OPM_IS_COMPILING_WITH_GPU_COMPILER
        gpuJacobian_->setToZero();
#endif
    }

    // Initialize the flows, flores, and velocity sparse tables
    void createFlows_()
    {
        OPM_TIMEBLOCK(createFlows);
        // If FLOWS/FLORES is set in any RPTRST in the schedule, then we initializate the sparse tables.
        // If DISPERC is in the deck, we initialize the sparse table here as well.
        const bool anyFlows = simulator_().problem().eclWriter().outputModule().getFlows().anyFlows();
        const auto& blockFlows = simulator_().problem().eclWriter().outputModule().getFlows().blockFlows();
        const auto& blockVelocity = simulator_().problem().eclWriter().outputModule().getFlows().blockVelocity();
        const bool isTemp = simulator_().vanguard().eclState().getSimulationConfig().isTemp();
        const bool anyFlores = simulator_().problem().eclWriter().outputModule().getFlows().anyFlores() || isTemp;
        const bool dispersionActive = simulator_().vanguard().eclState().getSimulationConfig().rock_config().dispersion();
        if (!dispersionActive && !enableBioeffects && blockVelocity.empty()
            && !((anyFlows || !blockFlows.empty()) && flowsInfo_.empty())
            && !(anyFlores && floresInfo_.empty())) {
            return;
        }
        const auto& model = model_();
        const auto& nncOutput = simulator_().problem().eclWriter().getOutputNnc().front();
        Stencil stencil(gridView_(), model_().dofMapper());
        const unsigned numCells = model.numTotalDof();
        std::unordered_multimap<int, std::pair<int, int>> nncIndices;
        std::vector<FlowInfo> loc_flinfo;
        std::vector<VelocityInfo> loc_vlinfo;
        unsigned int nncId = 0;
        VectorBlock flow(0.0);

        // Create a nnc structure to use fast lookup
        for (unsigned nncIdx = 0; nncIdx < nncOutput.size(); ++nncIdx) {
            const int ci1 = nncOutput[nncIdx].cell1;
            const int ci2 = nncOutput[nncIdx].cell2;
            nncIndices.emplace(ci1, std::make_pair(ci2, nncIdx));
        }

        if (anyFlows) {
            flowsInfo_.reserve(numCells, 6 * numCells);
        }
        else if (!blockFlows.empty()) {
            flowsInfo_.reserve(numCells, 6 * blockFlows.size());
        }
        if (anyFlores) {
            floresInfo_.reserve(numCells, 6 * numCells);
        }
        if (dispersionActive || enableBioeffects) {
            velocityInfo_.reserve(numCells, 6 * numCells);
        }
        else if (!blockVelocity.empty()) {
            velocityInfo_.reserve(numCells, 6 * blockVelocity.size());
        }

        for (const auto& elem : elements(gridView_())) {
            stencil.update(elem);
            for (unsigned primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++primaryDofIdx) {
                const unsigned myIdx = stencil.globalSpaceIndex(primaryDofIdx);
                bool blockFlowFound = false;
                bool blockVelocityFound = false;
                if (!blockFlows.empty()) {
                    if (std::ranges::binary_search(blockFlows,
                                                   simulator_().vanguard().cartesianIndex(myIdx))) {
                        blockFlowFound = true;
                    }
                    else {
                        flowsInfo_.appendRow(loc_flinfo.begin(), loc_flinfo.begin());
                        if (!dispersionActive && !enableBioeffects && !anyFlores && blockVelocity.empty()) {
                            continue;
                        }
                    }
                }
                if (!blockVelocity.empty() && !(dispersionActive || enableBioeffects)) {
                    if (std::ranges::binary_search(blockVelocity,
                                                   simulator_().vanguard().cartesianIndex(myIdx))) {
                        blockVelocityFound = true;
                    }
                    else {
                        velocityInfo_.appendRow(loc_vlinfo.begin(), loc_vlinfo.begin());
                        if (!anyFlows && blockFlows.empty() && !anyFlores) {
                            continue;
                        }
                    }
                }
                const int numFaces = stencil.numBoundaryFaces() + stencil.numInteriorFaces();
                loc_flinfo.resize(numFaces);
                loc_vlinfo.resize(stencil.numDof() - 1);

                for (unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
                    const unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
                    if (dofIdx > 0) {
                        const auto scvfIdx = dofIdx - 1;
                        const auto& scvf = stencil.interiorFace(scvfIdx);
                        int faceId = scvf.dirId();
                        const int cartMyIdx = simulator_().vanguard().cartesianIndex(myIdx);
                        const int cartNeighborIdx = simulator_().vanguard().cartesianIndex(neighborIdx);
                        const auto& range = nncIndices.equal_range(cartMyIdx);
                        for (auto it = range.first; it != range.second; ++it) {
                            if (it->second.first == cartNeighborIdx){
                                // -1 gives problem since is used for the nncInput from the deck
                                faceId = -2;
                                // the index is stored to be used for writting the outputs
                                nncId = it->second.second;
                            }
                        }
                        loc_flinfo[dofIdx - 1] = FlowInfo{faceId, flow, nncId};
                        loc_vlinfo[dofIdx - 1] = VelocityInfo{faceId, flow};
                    }
                }

                for (unsigned bdfIdx = 0; bdfIdx < stencil.numBoundaryFaces(); ++bdfIdx) {
                    const auto& scvf = stencil.boundaryFace(bdfIdx);
                    const int faceId = scvf.dirId();
                    loc_flinfo[stencil.numInteriorFaces() + bdfIdx] = FlowInfo{faceId, flow, nncId};
                }

                if (anyFlows || blockFlowFound) {
                    flowsInfo_.appendRow(loc_flinfo.begin(), loc_flinfo.end());
                }
                if (anyFlores) {
                    floresInfo_.appendRow(loc_flinfo.begin(), loc_flinfo.end());
                }
                if (dispersionActive || enableBioeffects || blockVelocityFound) {
                    velocityInfo_.appendRow(loc_vlinfo.begin(), loc_vlinfo.end());
                }
            }
        }
    }

public:
    template <class VectorBlockType, class MatrixBlockType, class ADVectorBlockType>
    OPM_HOST_DEVICE static void
    setResAndJacobi(VectorBlockType& res, MatrixBlockType& bMat, const ADVectorBlockType& resid)
    {
        for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
            res[eqIdx] = resid[eqIdx].value();
        }

        for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
            for (unsigned pvIdx = 0; pvIdx < numEq; ++pvIdx) {
                // A[dofIdx][focusDofIdx][eqIdx][pvIdx] is the partial derivative of
                // the residual function 'eqIdx' for the degree of freedom 'dofIdx'
                // with regard to the focus variable 'pvIdx' of the degree of freedom
                // 'focusDofIdx'
                bMat[eqIdx][pvIdx] = resid[eqIdx].derivative(pvIdx);
            }
        }
    }

    void updateFlowsInfo()
    {
        OPM_TIMEBLOCK(updateFlows);
        const bool enableFlows = simulator_().problem().eclWriter().outputModule().getFlows().hasFlows();
        const auto& blockFlows = simulator_().problem().eclWriter().outputModule().getFlows().blockFlows();
        // We reuse the fluxes in the TEMP option
        const bool isTemp = simulator_().vanguard().eclState().getSimulationConfig().isTemp();
        const bool enableFlores = simulator_().problem().eclWriter().outputModule().getFlows().hasFlores() || isTemp;
        if (!enableFlows && !enableFlores && blockFlows.empty()) {
            return;
        }
        const unsigned int numCells = model_().numTotalDof();
#ifdef _OPENMP
#pragma omp parallel for
#endif
        for (unsigned globI = 0; globI < numCells; ++globI) {
            OPM_TIMEBLOCK_LOCAL(linearizationForEachCell, Subsystem::Assembly);
            const auto& nbInfos = neighborInfo_[globI];
            ADVectorBlock adres(0.0);
            ADVectorBlock darcyFlux(0.0);
            const IntensiveQuantities& intQuantsIn = model_().intensiveQuantities(globI, /*timeIdx*/ 0);
            // Flux term.
            {
                OPM_TIMEBLOCK_LOCAL(fluxCalculationForEachCell, Subsystem::Assembly);
                short loc = 0;
                for (const auto& nbInfo : nbInfos) {
                    OPM_TIMEBLOCK_LOCAL(fluxCalculationForEachFace, Subsystem::Assembly);
                    const unsigned globJ = nbInfo.neighbor;
                    assert(globJ != globI);
                    adres = 0.0;
                    darcyFlux = 0.0;
                    const IntensiveQuantities& intQuantsEx = model_().intensiveQuantities(globJ, /*timeIdx*/ 0);
                    LocalResidual::computeFlux(adres, darcyFlux, globI, globJ, intQuantsIn,
                                               intQuantsEx, nbInfo.res_nbinfo, problem_().moduleParams());
                    adres *= nbInfo.res_nbinfo.faceArea;
                    if (!blockFlows.empty()) {
                        if (std::ranges::binary_search(blockFlows,
                                                       simulator_().vanguard().cartesianIndex(globI))) {
                            for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
                                flowsInfo_[globI][loc].flow[eqIdx] = adres[eqIdx].value();
                            }
                        }
                    }
                    else if (enableFlows) {
                        for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
                            flowsInfo_[globI][loc].flow[eqIdx] = adres[eqIdx].value();
                        }
                    }
                    if (enableFlores) {
                        for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
                            floresInfo_[globI][loc].flow[eqIdx] = darcyFlux[eqIdx].value();
                        }
                    }
                    ++loc;
                }
            }
        }

        // Boundary terms. Only looping over cells with nontrivial bcs.
        for (const auto& bdyInfo : boundaryInfo_) {
            if (bdyInfo.bcdata.type == BCType::NONE) {
                continue;
            }

            ADVectorBlockCPU adres(0.0);
            const unsigned globI = bdyInfo.cell;
            const auto& nbInfos = neighborInfo_[globI];
            const IntensiveQuantities& insideIntQuants = model_().intensiveQuantities(globI, /*timeIdx*/ 0);
            LocalResidual::computeBoundaryFlux(adres, problem_(), bdyInfo.bcdata, insideIntQuants, globI);
            adres *= bdyInfo.bcdata.faceArea;
            const unsigned bfIndex = bdyInfo.bfIndex;
            if (enableFlows) {
                for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
                    flowsInfo_[globI][nbInfos.size() + bfIndex].flow[eqIdx] = adres[eqIdx].value();
                }
            }
            // TODO also store Flores?
        }
    }

private:
    template <class SubDomainType>
    void linearize_(const SubDomainType& domain)
    {
        constexpr bool run_assembly_on_gpu = getPropValue<TypeTag, Properties::RunAssemblyOnGpu>();

        // This check should be removed once this is addressed by
        // for example storing the previous timesteps' values for
        // rsmax (for DRSDT) and similar.
        if (!problem_().recycleFirstIterationStorage()) {
            if (!model_().storeIntensiveQuantities() && !model_().enableStorageCache()) {
                OPM_THROW(std::runtime_error, "Must have cached either IQs or storage when we cannot recycle.");
            }
        }

        const double dt = simulator_().timeStepSize();

        OPM_TIMEBLOCK(linearize);

        // We do not call resetSystem_() here, since that will set
        // the full system to zero, not just our part.
        // Instead, that must be called before starting the linearization.
        const bool dispersionActive = simulator_().vanguard().eclState().getSimulationConfig().rock_config().dispersion();
        const unsigned int numCells = domain.cells.size();

        if constexpr (!run_assembly_on_gpu) {
            linearize_parallelization_wrapper<run_assembly_on_gpu, LocalResidual>(
                numCells /*numCells*/,
                domain,
                neighborInfo_,
                diagMatAddress_,
                residual_,
                model_(),
                dt,
                dispersionActive,
                problem_());

            linearize_bc<IntensiveQuantities, Model, LocalResidual>(
                diagMatAddress_, residual_, boundaryInfo_);
        } else {
#if HAVE_CUDA && OPM_IS_COMPILING_WITH_GPU_COMPILER
            // Make sure we have can have the domain on the GPU.
            if constexpr (std::is_same_v<SubDomainType, FullDomain<>>) {

                using GpuParams = TpfaLinearizerGpuParams<TypeTag>;

                int constexpr blockSize = 256; // Experimentally this is a good value for multiple
                                               // GPUs. Autotune this later.

                GpuParams gpuParams(domain,
                                    neighborInfo_,
                                    *gpuJacobian_,
                                    *gpuBufferDiagMatAddress_,
                                    jacobian_->istlMatrix(),
                                    residual_,
                                    boundaryInfo_,
                                    model_(),
                                    problem_(),
                                    numCells);

                linearize_parallelization_wrapper<run_assembly_on_gpu,
                                                  typename GpuParams::LocalResidualGPU>(
                    numCells,
                    gpuParams.domainView(),
                    gpuParams.neighborInfoView(),
                    gpuParams.diagMatAddressView(),
                    gpuParams.residualView(),
                    gpuParams.modelView(),
                    dt,
                    dispersionActive,
                    gpuParams.flowProblemView());

                if (gpuParams.boundaryInfoSize() > 0) {
                    auto boundaryInfoView = gpuParams.boundaryInfoView();
                    linearize_bc_threadsafe<TpfaLinearizer<TypeTag>,
                                            typename GpuParams::GPUBOIQ,
                                            decltype(gpuParams.modelView()),
                                            typename GpuParams::LocalResidualGPU>
                        <<<((gpuParams.boundaryInfoSize() + blockSize - 1) / blockSize),
                           blockSize>>>(gpuParams.diagMatAddressView(),
                                        gpuParams.residualView(),
                                        boundaryInfoView,
                                        gpuParams.modelView(),
                                        gpuParams.flowProblemView());
                }

                // The memory copies here are synchronous and in the default stream,
                // guaranteeing that the GPU kernels have completed.
                gpuParams.copyResidualToHost(residual_, numCells);
                gpuParams.copyJacobianToHost(*jacobian_, *gpuJacobian_);
            } else {
                OPM_THROW(std::logic_error, "Only FullDomain is supported on GPU");
            }
#else
            OPM_THROW(std::logic_error,
                      "Trying to run GPU assembly without compiling with GPU support");
#endif
        }

        // Handle source terms separately as we want the functionality for CPU and GPU cases
        // but for now we cannot handle this inside gpu kernels due to the use of problem_()
        linearize_source_terms(numCells, domain);
    }

    template <bool useGPU,
              class LocalResidualT,
              class ModelClass,
              class DiagPtrType,
              class DomainType,
              class NeighborSparseTable,
              class ResidualType,
              class ProblemT>
    void linearize_parallelization_wrapper(const unsigned int numCells,
                                           const DomainType& domain,
                                           const NeighborSparseTable& neighborInfo,
                                           DiagPtrType& diagMatAddress,
                                           ResidualType& residual,
                                           const ModelClass& model,
                                           Scalar dt,
                                           [[maybe_unused]] bool dispersionActive,
                                           const ProblemT& problem)
    {
        if constexpr (useGPU) {
            static_assert(!enableBioeffects && "Bioeffects not yet supported on GPU");
            assert(!dispersionActive && "Dispersion not yet supported on GPU");
#if HAVE_CUDA && OPM_IS_COMPILING_WITH_GPU_COMPILER
            int constexpr blockSize = 256;
            kernel_linearize<TpfaLinearizer<TypeTag>,
                             ModelClass,
                             LocalResidualT,
                             DiagPtrType,
                             Scalar,
                             DomainType,
                             NeighborSparseTable,
                             ResidualType,
                             ProblemT>
                <<<((numCells + blockSize - 1) / blockSize), blockSize>>>(numCells,
                                                                          domain,
                                                                          neighborInfo,
                                                                          diagMatAddress,
                                                                          residual,
                                                                          model,
                                                                          dt,
                                                                          problem);
#else
            OPM_THROW(std::runtime_error, "Trying to run GPU code without GPU support");
#endif
        } else {
#ifdef _OPENMP
#pragma omp parallel for
#endif
            for (unsigned ii = 0; ii < numCells; ++ii) {
                linearize_cell<false, LocalResidual>(ii,
                                                     domain,
                                                     neighborInfo,
                                                     diagMatAddress,
                                                     residual,
                                                     model,
                                                     velocityInfo_,
                                                     dt,
                                                     problem);
            }
        }
    }

    template <class SubDomainType>
    void linearize_source_terms(unsigned int numCells, const SubDomainType& domain)
    {
#ifdef _OPENMP
#pragma omp parallel for
#endif
        for (unsigned ii = 0; ii < numCells; ++ii) {
            OPM_TIMEBLOCK_LOCAL(linearizationForEachCell, Subsystem::Assembly);
            const unsigned globI = domain.cells[ii];
            VectorBlockCPU res(0.0);
            MatrixBlockCPU bMat(0.0);
            ADVectorBlockCPU adres(0.0);
            const IntensiveQuantities& intQuantsIn = model_().intensiveQuantities(globI, /*timeIdx*/ 0);
            const double volume = model_().dofTotalVolume(globI);

            // Cell-wise source terms.
            // This will include well sources if SeparateSparseSourceTerms is false.
            res = 0.0;
            bMat = 0.0;
            adres = 0.0;
            if (separateSparseSourceTerms_) {
                LocalResidual::computeSourceDense(adres, problem_(), intQuantsIn, globI, 0);
            } else {
                LocalResidual::computeSource(adres, problem_(), intQuantsIn, globI, 0);
            }
            adres *= -volume;
            setResAndJacobi(res, bMat, adres);
            residual_[globI] += res;
            // SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat);
            *diagMatAddress_[globI] += bMat;
        } // end of loop for cell globI.

        // Add sparse source terms. For now only wells.
        if (separateSparseSourceTerms_) {
            problem_().wellModel().addReservoirSourceTerms(residual_, diagMatAddress_);
        }
    }

public:
    template <bool useGPU,
              class LocalResidualT,
              class ProblemType,
              class VelocityInfoType,
              class ModelClass,
              class DiagPtrType,
              class DomainType,
              class NeighborSparseTable,
              class ResidualType>
    OPM_HOST_DEVICE static void linearize_cell(const unsigned int ii,
                                               const DomainType& domain,
                                               const NeighborSparseTable& neighborInfo,
                                               const DiagPtrType& diagMatAddress,
                                               ResidualType& residual,
                                               const ModelClass& model,
                                               VelocityInfoType& velocityInfo,
                                               const Scalar dt,
                                               const ProblemType& problem)
    {
#if OPM_IS_INSIDE_HOST_FUNCTION
        OPM_TIMEBLOCK_LOCAL(linearizationForEachCell, Subsystem::Assembly);
#endif
        const unsigned globI = domain.cells[ii];
        const auto& nbInfos = neighborInfo[globI];
        VectorBlock res(0.0);
        MatrixBlock bMat(0.0);
        ADVectorBlock adres(0.0);
        ADVectorBlock darcyFlux(0.0);
        const auto& intQuantsIn = model.intensiveQuantities(globI, /*timeIdx*/ 0);

        // Flux term.
        {
#if OPM_IS_INSIDE_HOST_FUNCTION
            OPM_TIMEBLOCK_LOCAL(fluxCalculationForEachCell, Subsystem::Assembly);
#endif
            short loc = 0;
            for (const auto& nbInfo : nbInfos) {
                const unsigned globJ = nbInfo.neighbor;
                res = 0.0;
                bMat = 0.0;
                adres = 0.0;
                darcyFlux = 0.0;

                const auto& intQuantsEx = model.intensiveQuantities(globJ, /*timeIdx*/ 0);

                LocalResidualT::computeFlux(adres,
                                            darcyFlux,
                                            globI,
                                            globJ,
                                            intQuantsIn,
                                            intQuantsEx,
                                            nbInfo.res_nbinfo,
                                            problem.moduleParams());

                adres *= nbInfo.res_nbinfo.faceArea;

                // currently not supported on GPU
                if constexpr (!useGPU) {
                    // Use std::cmp_less because sparseTable.size() returns signed integer, and globI is unsigned
                    if (std::cmp_less(globI, velocityInfo.size())) {
                        if (velocityInfo.rowSize(globI) > 0) {
                            for (unsigned phaseIdx = 0; phaseIdx < numEq; ++phaseIdx) {
                                velocityInfo[globI][loc].velocity[phaseIdx]
                                    = darcyFlux[phaseIdx].value() / nbInfo.res_nbinfo.faceArea;
                            }
                        }
                    }
                }

                setResAndJacobi(res, bMat, adres);
                residual[globI] += res;
                // SparseAdapter syntax:  jacobian_->addToBlock(globI, globI, bMat);
                *diagMatAddress[globI] += bMat;
                bMat *= -1.0;
                // SparseAdapter syntax: jacobian_->addToBlock(globJ, globI, bMat);
                *nbInfo.matBlockAddress += bMat;
                ++loc;
            }
        }

        // Accumulation term.
        adres = 0.0;
        {
            LocalResidualT::template computeStorage<Evaluation>(adres, intQuantsIn);
        }
        setResAndJacobi(res, bMat, adres);

        if constexpr (!useGPU) { // Cached storage not enabled for GPU
            // Either use cached storage term, or compute it on the fly.
            if (model.enableStorageCache()) {
                // The cached storage for timeIdx 0 (current time) is not
                // used, but after storage cache is shifted at the end of the
                // timestep, it will become cached storage for timeIdx 1.
                model.updateCachedStorage(globI, /*timeIdx=*/0, res);
                // We should not update the storage cache here for NLDD local solves.
                // This will reset the start-of-step storage to incorrect numbers when
                // we do local solves, where the iteration number will start from 0,
                // but the starting state may not be identical to the start-of-step state.
                // Note that a full assembly must be done before local solves
                // otherwise this will be left un-updated.
                if (problem.iterationContext().isFirstGlobalIteration()) {
                    // Need to update the storage cache.
                    if (problem.recycleFirstIterationStorage()) {
                        // Assumes nothing have changed in the system which
                        // affects masses calculated from primary variables.
                        model.updateCachedStorage(globI, 1, res);
                    } else {
                        VectorBlock tmp;
                        const auto& intQuantOld = model.intensiveQuantities(globI, 1);
                        LocalResidualT::template computeStorage<Scalar>(tmp, intQuantOld);
                        model.updateCachedStorage(globI, 1, tmp);
                    }
                }
                res -= model.cachedStorage(globI, 1);
            } else {
#if OPM_IS_INSIDE_HOST_FUNCTION
                OPM_TIMEBLOCK_LOCAL(computeStorage0, Subsystem::Assembly);
#endif
                VectorBlock tmp;
                const auto& intQuantOld = model.intensiveQuantities(globI, 1);
                LocalResidualT::template computeStorage<Scalar>(tmp, intQuantOld);
                // assume volume do not change
                res -= tmp;
            }
        } else {
            VectorBlock tmp;
            const auto& intQuantOld = model.intensiveQuantities(globI, 1);
            LocalResidualT::template computeStorage<Scalar>(tmp, intQuantOld);
            // assume volume do not change
            res -= tmp;
        }

        const Scalar volume = model.dofTotalVolume(globI);
        const Scalar storefac = volume / dt;
        res *= storefac;
        bMat *= storefac;
        residual[globI] += res;
        // SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat);
        *diagMatAddress[globI] += bMat;
    }

#if HAVE_CUDA && OPM_IS_COMPILING_WITH_GPU_COMPILER
    template <class LocalIntensiveQuantities,
              class LocalResidualT,
              class ModelT,
              class DiagPtrType,
              class ResidualT,
              class BoundaryInfoT,
              class ProblemT>
    OPM_HOST_DEVICE static void linearize_bc_threadsafe_single_cell(DiagPtrType diagMatAdress,
                                                               ResidualT residual,
                                                               const BoundaryInfoT boundaryInfoElement,
                                                               ModelT model,
                                                               ProblemT problem)
    {
        using Scalar = GetPropType<TypeTag, Properties::Scalar>;
        constexpr int numEq = getPropValue<TypeTag, Properties::NumEq>();
        if (boundaryInfoElement.bcdata.type != BCType::NONE) {
            VectorBlock res(0.0);
            MatrixBlock bMat(0.0);
            ADVectorBlock adres(0.0);
            const unsigned globI = boundaryInfoElement.cell;
            const LocalIntensiveQuantities& insideIntQuants
                = model.intensiveQuantities(globI, /*timeIdx*/ 0);
            if constexpr (!std::is_empty_v<GetPropType<TypeTag, Properties::FluidSystem>>) {
                LocalResidualT::computeBoundaryFlux(
                    adres, problem, boundaryInfoElement.bcdata, insideIntQuants, globI);
            }
            adres *= boundaryInfoElement.bcdata.faceArea;
            TpfaLinearizer<TypeTag>::setResAndJacobi(res, bMat, adres);
            auto* residualPtr = &(residual.data()[globI]);
            for (int i = 0; i < numEq; ++i) {
                atomicAdd(&((*residualPtr)[i]), res[i]);
            }
            auto* matPtr = diagMatAdress[globI];
            for (int row = 0; row < bMat.size(); ++row) {
                for (int col = 0; col < bMat.size(); ++col) {
                    Scalar* elemPtr = &((*matPtr)[row][col]);
                    atomicAdd(elemPtr, bMat[row][col]);
                }
            }
        }
    }
#endif

private:
    template <class LocalIntensiveQuantities,
              class ModelClass,
              class LocalResidualT,
              class DiagPtrType,
              class ResidualType,
              class BoundaryInfoT>
    void linearize_bc(DiagPtrType& diagMatAdress,
                      ResidualType& residual,
                      const BoundaryInfoT& boundaryInfo)
    {
        // Boundary terms. Only looping over cells with nontrivial bcs.
        for (const auto& bdyInfo : boundaryInfo) {
            if (bdyInfo.bcdata.type == BCType::NONE) {
                continue;
            }
            VectorBlock res(0.0);
            MatrixBlock bMat(0.0);
            ADVectorBlock adres(0.0);
            const unsigned globI = bdyInfo.cell;
            const LocalIntensiveQuantities& insideIntQuants
                = model_().intensiveQuantities(globI, /*timeIdx*/ 0);
            LocalResidual::computeBoundaryFlux(
                adres, problem_(), bdyInfo.bcdata, insideIntQuants, globI);
            adres *= bdyInfo.bcdata.faceArea;
            setResAndJacobi(res, bMat, adres);
            residual[globI] += res;
            *diagMatAdress[globI] += bMat;
        }
    }

    void updateStoredTransmissibilities()
    {
        if (neighborInfo_.empty()) {
            // This function was called before createMatrix_() was called.
            // We call initFirstIteration_(), not createMatrix_(), because
            // that will also initialize the residual consistently.
            initFirstIteration_();
        }

        const unsigned numCells = model_().numTotalDof();
#ifdef _OPENMP
#pragma omp parallel for
#endif
        for (unsigned globI = 0; globI < numCells; globI++) {
            auto nbInfos = neighborInfo_[globI]; // nbInfos will be a SparseTable<...>::mutable_iterator_range.
            for (auto& nbInfo : nbInfos) {
                const unsigned globJ = nbInfo.neighbor;
                nbInfo.res_nbinfo.trans = problem_().transmissibility(globI, globJ);
            }
        }
    }

    Simulator* simulatorPtr_{};

    // the jacobian matrix
    std::unique_ptr<SparseMatrixAdapter> jacobian_{};
#if HAVE_CUDA
    std::unique_ptr<gpuistl::GpuSparseMatrixWrapper<Scalar>> gpuJacobian_;
    std::unique_ptr<gpuistl::GpuBuffer<MatrixBlockGPU*>> gpuBufferDiagMatAddress_;
#endif

    // the right-hand side
    GlobalEqVector residual_;

    LinearizationType linearizationType_{};

    using ResidualNBInfo = typename LocalResidual::ResidualNBInfo;
    using NeighborInfoCPU = NeighborInfoStruct<ResidualNBInfo, MatrixBlockCPU>;

    SparseTable<NeighborInfoCPU> neighborInfo_{};
    std::vector<MatrixBlockCPU*> diagMatAddress_ {};

    struct FlowInfo
    {
        int faceId;
        VectorBlock flow;
        unsigned int nncId;
    };
    SparseTable<FlowInfo> flowsInfo_;
    SparseTable<FlowInfo> floresInfo_;

    struct VelocityInfo
    {
        int faceId;
        VectorBlock velocity;
    };
    SparseTable<VelocityInfo> velocityInfo_;

    using ScalarFluidState = typename IntensiveQuantities::ScalarFluidState;

    using BoundaryConditionDataCPU = BoundaryConditionData<VectorBlockCPU, ScalarFluidState>;

    using BoundaryInfoCPU = BoundaryInfo<BoundaryConditionDataCPU>;

    std::vector<BoundaryInfoCPU> boundaryInfo_;

    bool separateSparseSourceTerms_ = false;

    FullDomain<> fullDomain_;

    int exportIndex_;
    int exportCount_;
};
} // namespace Opm

#endif // TPFA_LINEARIZER