tpsalinearizer.hpp

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:
/*
  Copyright 2025 NORCE AS
 
  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 TPSA_LINEARIZER_HPP
#define TPSA_LINEARIZER_HPP
 
#include <dune/common/fvector.hh>
 
#include <opm/common/TimingMacros.hpp>
 
#include <opm/grid/utility/SparseTable.hpp>
 
#include <opm/input/eclipse/Schedule/BCState.hpp>
 
#include <opm/material/materialstates/MaterialStateTPSA.hpp>
 
#include <opm/models/discretization/common/linearizationtype.hh>
#include <opm/models/tpsa/tpsabaseproperties.hpp>
 
#include <cassert>
#include <set>
#include <vector>
 
 
namespace Opm {
 
template<class TypeTag>
class TpsaLinearizer
{
    using Constraints = GetPropType<TypeTag, Properties::Constraints>;  // TODO: make TPSA constraints
    using Evaluation = GetPropType<TypeTag, Properties::EvaluationTPSA>;
    using FlowModel = GetPropType<TypeTag, Properties::Model>;
    using GeomechModel = GetPropType<TypeTag, Properties::ModelTPSA>;
    using GlobalEqVector = GetPropType<TypeTag, Properties::GlobalEqVectorTPSA>;
    using GridView = GetPropType<TypeTag, Properties::GridView>;
    using LocalResidual = GetPropType<TypeTag, Properties::LocalResidualTPSA>;
    using Problem = GetPropType<TypeTag, Properties::Problem>;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Simulator = GetPropType<TypeTag, Properties::Simulator>;
    using SparseMatrixAdapter = GetPropType<TypeTag, Properties::SparseMatrixAdapterTPSA>;
    using Stencil = GetPropType<TypeTag, Properties::Stencil>;
 
    using MaterialState = MaterialStateTPSA<Evaluation>;
 
    enum { numEq = getPropValue<TypeTag, Properties::NumEqTPSA>() };
 
    using ADVectorBlock = Dune::FieldVector<Evaluation, numEq>;
    using MatrixBlock = typename SparseMatrixAdapter::MatrixBlock;
    using BlockAddress = typename SparseMatrixAdapter::BlockAddress;
    using VectorBlock = Dune::FieldVector<Scalar, numEq>;
 
    using StressInfoVector = Dune::FieldVector<Scalar, 3>;
 
public:
    // ///
    // Public functions
    // ///
    TpsaLinearizer()
    {
        simulatorPtr_ = nullptr;
    }
 
    void init(Simulator& simulator)
    {
        simulatorPtr_ = &simulator;
        eraseMatrix();
    }
 
    void eraseMatrix()
    {
        jacobian_.reset();
    }
 
    void finalize()
    {
        jacobian_->finalize();
    }
 
    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);
        }
    }
 
    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 TPSA system:" << e.what()
                      << "\n"  << std::flush;
            succeeded = 0;
        }
        catch (...) {
            std::cout << "rank " << simulator_().gridView().comm().rank()
                      << " caught an exception while linearizing TPSA system"
                      << "\n"  << std::flush;
            succeeded = 0;
        }
        succeeded = simulator_().gridView().comm().min(succeeded);
 
        if (!succeeded) {
            throw NumericalProblem("A process did not succeed in linearizing the TPSA 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_();
        }
 
        if (domain.cells.size() == flowModel_().numTotalDof()) {
            // We are on the full domain.
            resetSystem_();
        }
        else {
            resetSystem_(domain);
        }
 
        linearize_(domain);
    }
 
    void linearizeAuxiliaryEquations()
    { }
 
    void setResAndJacobi(VectorBlock& res, MatrixBlock& bMat, const ADVectorBlock& resid) const
    {
        // Scalar local residual
        for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
            res[eqIdx] = resid[eqIdx].value();
        }
 
        // 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'
        for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
            for (unsigned pvIdx = 0; pvIdx < numEq; ++pvIdx) {
                bMat[eqIdx][pvIdx] = resid[eqIdx].derivative(pvIdx);
            }
        }
    }
 
    void updateBoundaryConditionData()
    {
        for (auto& bdyInfo : boundaryInfo_) {
            // Get boundary information from problem
            const auto& mechBC = problem_().mechBoundaryCondition(bdyInfo.cell, bdyInfo.dir);
 
            // Strip the unnecessary (and zero anyway) derivatives off displacement
            std::vector<double> displacement(3, 0.0);
            for (std::size_t ii = 0; ii < displacement.size(); ++ii) {
                displacement[ii] = mechBC.displacement[ii].value();
            }
 
            // Update boundary information
            bdyInfo.bcdata.type = mechBC.type;
            bdyInfo.bcdata.displacement = displacement;
            bdyInfo.bcdata.shearModulus = mechBC.shearModulus;
            bdyInfo.bcdata.distance = mechBC.distance;
        }
    }
 
    void updateStressInfo()
    {
        OPM_TIMEBLOCK(updateStressInfoTPSA);
 
        if (stressInfo_.empty()) {
            return;
        }
 
        const auto& geoMechModel = geoMechModel_();
        auto& problem = problem_();
        const unsigned int numCells = fullDomain_.cells.size();
 
#ifdef _OPENMP
#pragma omp parallel for
#endif
        // Interior traction of every cell, one per face to its neighbors
        for (unsigned ii = 0; ii < numCells; ++ii) {
            const unsigned globI = fullDomain_.cells[ii];
            const MaterialState& materialStateIn =
                geoMechModel.materialState(globI, /*timeIdx=*/0);
 
            short loc = 0;
            for (const auto& nbInfo : neighborInfo_[globI]) {
                const unsigned globJ = nbInfo.neighbor;
                const MaterialState& materialStateEx =
                    geoMechModel.materialState(globJ, /*timeIdx=*/0);
 
                // Compute local face term
                ADVectorBlock adres(0.0);
                LocalResidual::computeFaceTerm(adres,
                                               materialStateIn,
                                               materialStateEx,
                                               problem,
                                               globI,
                                               globJ);
                adres *= nbInfo.faceArea;
 
                // Insert interior traction vector
                // OBS: Assume traction vector is the three first entries in residual!
                stressInfo_[globI][loc].faceNormal = problem.cellFaceNormal(globI, globJ);
                stressInfo_[globI][loc].faceArea = nbInfo.faceArea;
                for (unsigned tractionIdx = 0; tractionIdx < 3; ++tractionIdx) {
                    stressInfo_[globI][loc].traction[tractionIdx] = adres[tractionIdx].value();
                }
                ++loc;
            }
        }
 
        // Boundary traction
        for (const auto& bdyInfo : boundaryInfo_) {
            const unsigned globI = bdyInfo.cell;
            const MaterialState& materialStateIn = geoMechModel.materialState(globI, /*timeIdx=*/0);
 
            // Compute local boundary condition
            ADVectorBlock adres(0.0);
            LocalResidual::computeBoundaryTerm(adres,
                                               materialStateIn,
                                               bdyInfo.bcdata,
                                               problem,
                                               globI);
            adres *= bdyInfo.bcdata.faceArea;
 
            // Insert boundary traction vector
            // OBS: Assume traction vector is the three first entries in residual!
            const short loc = neighborInfo_[globI].size() + bdyInfo.bfIndex;
            stressInfo_[globI][loc].faceNormal =
                problem.cellFaceNormalBoundary(globI, bdyInfo.bfIndex);
            stressInfo_[globI][loc].faceArea = bdyInfo.bcdata.faceArea;
            for (unsigned tractionIdx = 0; tractionIdx < 3; ++tractionIdx) {
                stressInfo_[globI][loc].traction[tractionIdx] = adres[tractionIdx].value();
            }
        }
    }
 
    // ///
    // Public get and set functions
    // ///
    const SparseMatrixAdapter& jacobian() const
    { return *jacobian_; }
 
    SparseMatrixAdapter& jacobian()
    { return *jacobian_; }
 
    const GlobalEqVector& residual() const
    { return residual_; }
 
    GlobalEqVector& residual()
    { return residual_; }
 
    const LinearizationType& getLinearizationType() const
    { return linearizationType_; }
 
    const auto& getStressInfo() const
    {
        return stressInfo_;
    }
 
    std::map<unsigned, Constraints> constraintsMap() const
    { return {}; }
 
    void setLinearizationType(LinearizationType linearizationType)
    { linearizationType_ = linearizationType; }
 
private:
    // ///
    // Private functions
    // ///
    void createMatrix_()
    {
        OPM_TIMEBLOCK(createMatrixTPSA);
 
        // If the Jacobian has been initialize before, we jump out
        if (!neighborInfo_.empty()) {
            return;
        }
 
        // Init. the stencil
        const auto& flowModel = flowModel_();
        Stencil stencil(gridView_(), flowModel.dofMapper());
 
        // Build up sparsity patterns and neighboring information for Jacobian and linearization
        std::vector<std::set<unsigned>> sparsityPattern(flowModel.numTotalDof());
        unsigned numCells = flowModel.numTotalDof();
        neighborInfo_.reserve(numCells, 6 * numCells);
        std::vector<NeighborInfo> loc_nbinfo;
 
        stressInfo_.reserve(numCells, 6 * numCells);
        std::vector<StressInfo> loc_stressinfo;
        StressInfoVector nullTraction(0.0);
        StressInfoVector nullFaceNormal(0.0);
        for (const auto& elem : elements(gridView_())) {
            // Loop over primary dofs in the element
            stencil.update(elem);
 
            for (unsigned primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++primaryDofIdx) {
                // Build up neighboring information for curret primary dof
                const unsigned myIdx = stencil.globalSpaceIndex(primaryDofIdx);
                loc_nbinfo.resize(stencil.numDof() - 1);
 
                // Build information for stress output
                const int numFaces = stencil.numBoundaryFaces() + stencil.numInteriorFaces();
                loc_stressinfo.resize(numFaces);
 
                for (unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
                    // NOTE: NeighborInfo could/should be expanded with cell face parameters located in problem_()
                    // needed when computing face terms in LocalResidual
                    const unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
                    sparsityPattern[myIdx].insert(neighborIdx);
                    if (dofIdx > 0) {
                        const auto scvfIdx = dofIdx - 1;
                        const auto& scvf = stencil.interiorFace(scvfIdx);
                        const Scalar area = scvf.area();
                        loc_nbinfo[dofIdx - 1] = NeighborInfo{ neighborIdx, area, BlockAddress{} };
 
                        int faceId = scvf.dirId();
                        loc_stressinfo[dofIdx - 1] =
                            StressInfo{faceId, nullTraction, nullFaceNormal, area};
                    }
                }
 
                // Insert local neighbor info
                neighborInfo_.appendRow(loc_nbinfo.begin(), loc_nbinfo.end());
 
                // Boundary condition information
                unsigned bfIndex = 0;
                for (const auto& intersection : intersections(gridView_(), elem)) {
                    if (intersection.boundary()) {
                        // Get boundary face direction
                        const auto& bf = stencil.boundaryFace(bfIndex);
                        const int dir_id = bf.dirId();
                        const auto bfArea = bf.area();
 
                        // Initialize boundary information for stress output container
                        // OBS: NNCs are not implemented with TPSA, hence traction and face
                        // normal vector will stay zero for NNC connections!
                        loc_stressinfo[stencil.numInteriorFaces() + bfIndex] =
                            StressInfo{dir_id, nullTraction, nullFaceNormal, bfArea};
 
                        // Skip NNCs
                        if (dir_id < 0) {
                            continue;
                        }
 
                        // Get boundary information from problem()
                        const auto& mechBC = problem_().mechBoundaryCondition(myIdx, dir_id);
 
                        // Strip the unnecessary (and zero anyway) derivatives off displacement
                        std::vector<double> displacement(3, 0.0);
                        for (std::size_t ii = 0; ii < displacement.size(); ++ii) {
                            displacement[ii] = mechBC.displacement[ii].value();
                        }
 
                        // Insert boundary condition data in container
                        BoundaryConditionData bcdata{mechBC.type,
                                                     displacement,
                                                     mechBC.shearModulus,
                                                     mechBC.distance,
                                                     bfIndex,
                                                     bfArea};
                        boundaryInfo_.push_back( { myIdx, dir_id, bfIndex, bcdata } );
                        ++bfIndex;
                        continue;
                    }
                    if (!intersection.neighbor()) {
                        ++bfIndex;
                        continue;
                    }
                }
                stressInfo_.appendRow(loc_stressinfo.begin(), loc_stressinfo.end());
            }
        }
 
        // Allocate Jacobian matrix and pointers to its sub-blocks
        jacobian_ = std::make_unique<SparseMatrixAdapter>(simulator_());
        diagMatAddress_.resize(numCells);
        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);
            }
        }
 
        // Create full domain
        fullDomain_.cells.resize(numCells);
        std::iota(fullDomain_.cells.begin(), fullDomain_.cells.end(), 0);
    }
 
    template <class SubDomainType>
    void linearize_(const SubDomainType& domain)
    {
        OPM_TIMEBLOCK(linearizeTPSA);
 
        // Extract misc. variables used in linearization
        const auto& flowModel = flowModel_();
        const auto& geoMechModel = geoMechModel_();
        auto& problem = problem_();
        const unsigned int numCells = domain.cells.size();
 
#ifdef _OPENMP
#pragma omp parallel for
#endif
        // Loop over cells in the domain and compute local residual and jacobian
        for (unsigned ii = 0; ii < numCells; ++ii) {
            OPM_TIMEBLOCK_LOCAL(linearizationForEachCellTPSA, Subsystem::Assembly);
 
            const unsigned globI = domain.cells[ii];
            const auto& nbInfos = neighborInfo_[globI];
            VectorBlock res(0.0);
            MatrixBlock bMat(0.0);
            ADVectorBlock adres(0.0);
            const MaterialState& materialStateIn = geoMechModel.materialState(globI, /*timeIdx=*/0);
 
            // ///
            // Face term
            // ///
            {
                OPM_TIMEBLOCK_LOCAL(faceCalculationForEachCellTPSA, Subsystem::Assembly);
 
                // Loop over neighboring cells
                for (auto& nbInfo : nbInfos) {
                    OPM_TIMEBLOCK_LOCAL(calculationForEachFaceTPSA, Subsystem::Assembly);
 
                    // Reset local residual and Jacobian
                    res = 0.0;
                    bMat = 0.0;
                    adres = 0.0;
 
                    // Neighbor information
                    const unsigned globJ = nbInfo.neighbor;
                    assert(globJ != globI);
                    const MaterialState& materialStateEx = geoMechModel.materialState(globJ, /*timeIdx=*/0);
 
                    // Compute local face term
                    LocalResidual::computeFaceTerm(adres,
                                                   materialStateIn,
                                                   materialStateEx,
                                                   problem,
                                                   globI,
                                                   globJ);
                    adres *= nbInfo.faceArea;
 
                    // Extract residual and sub-block Jacobian entries from computed AD residual
                    setResAndJacobi(res, bMat, adres);
 
                    // Insert into global residual
                    residual_[globI] += res;
 
                    // Insert contribution to (globI, globI) sub-block
                    // SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat);
                    *diagMatAddress_[globI] += bMat;
 
                    // Insert contribution to (globJ, globI) sub-block
                    // Note: since LocalResidual::computeFaceTerm have Evaluation on globI primary variables, it is
                    // natural to insert Jacobian entries for (globJ, globI) here, since we only need to flip the signs
                    // in the calculated face terms
                    // SparseAdapter syntax: jacobian_->addToBlock(globJ, globI, bMat);
                    bMat *= -1.0;
                    *nbInfo.matBlockAddress += bMat;
                }
            }
 
            // ///
            // Volume term
            // //
            adres = 0.0;
            {
                OPM_TIMEBLOCK_LOCAL(computeVolumeTerm, Subsystem::Assembly);
 
                // Compute local volume term
                LocalResidual::computeVolumeTerm(adres,
                                                 materialStateIn,
                                                 problem,
                                                 globI);
            }
            const double volume = flowModel.dofTotalVolume(globI);
            adres *= volume;
 
            // Extract residual and sub-block Jacobian entries from computed AD residual
            setResAndJacobi(res, bMat, adres);
 
            // Insert in global residual
            residual_[globI] += res;
 
            // Insert contribution to (globI, globI) sub-block
            // SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat);
            *diagMatAddress_[globI] += bMat;
 
            // ///
            // Source term
            // ///
            res = 0.0;
            bMat = 0.0;
            adres = 0.0;
 
            // Compute local source term
            LocalResidual::computeSourceTerm(adres,
                                             problem,
                                             globI,
                                             0);
            adres *= -volume;
 
            // Extract residual and sub-block Jacobian entries from computed AD residual
            setResAndJacobi(res, bMat, adres);
 
            // Insert into global residual
            residual_[globI] += res;
 
            // Insert contribution to (globI, globI) sub-block
            // SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat);
            *diagMatAddress_[globI] += bMat;
        }  // globI loop
 
        // ///
        // Boundary term
        // ///
        for (const auto& bdyInfo : boundaryInfo_) {
            VectorBlock res(0.0);
            MatrixBlock bMat(0.0);
            ADVectorBlock adres(0.0);
            const unsigned globI = bdyInfo.cell;
            const MaterialState& materialStateIn = geoMechModel.materialState(globI, /*timeIdx=*/0);
 
            // Compute local boundary condition
            LocalResidual::computeBoundaryTerm(adres,
                                               materialStateIn,
                                               bdyInfo.bcdata,
                                               problem,
                                               globI);
            adres *= bdyInfo.bcdata.faceArea;
 
            // Extract residual and sub-block Jacobian entries from computed AD residual
            setResAndJacobi(res, bMat, adres);
 
            // Insert in global residual
            residual_[globI] += res;
 
            // Insert contribution to (globI, globI) sub-block
            // SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat);
            *diagMatAddress_[globI] += bMat;
        }
    }
 
    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(flowModel_().numTotalDof());
        resetSystem_();
    }
 
    void resetSystem_()
    {
        // Set residual vector entries to zero
        residual_ = 0.0;
 
        // Set all Jacobian matrix entries to zero
        jacobian_->clear();
    }
 
    // ///
    // Private get functions
    // ///
    Simulator& simulator_()
    { return *simulatorPtr_; }
 
    const Simulator& simulator_() const
    { return *simulatorPtr_; }
 
    Problem& problem_()
    { return simulator_().problem(); }
 
    const Problem& problem_() const
    { return simulator_().problem(); }
 
    FlowModel& flowModel_()
    { return simulator_().model(); }
 
    const FlowModel& flowModel_() const
    { return simulator_().model(); }
 
    GeomechModel& geoMechModel_()
    { return problem_().geoMechModel(); }
 
    const GeomechModel& geoMechModel_() const
    { return problem_().geoMechModel(); }
 
    const GridView& gridView_() const
    { return problem_().gridView(); }
 
    Simulator* simulatorPtr_{};
    LinearizationType linearizationType_{};
 
    std::vector<BlockAddress> diagMatAddress_{};
    std::unique_ptr<SparseMatrixAdapter> jacobian_{};
    GlobalEqVector residual_;
 
    //
    // Helper structs
    //
    struct NeighborInfo
    {
        unsigned int neighbor;
        double faceArea;
        BlockAddress matBlockAddress;
    };
    SparseTable<NeighborInfo> neighborInfo_{};
 
    struct StressInfo
    {
        int faceId;
        StressInfoVector traction;
        StressInfoVector faceNormal;
        double faceArea;
    };
 
    SparseTable<StressInfo> stressInfo_{};
 
    struct BoundaryConditionData
    {
        BCMECHType type;
        std::vector<double> displacement;
        Scalar shearModulus;
        Scalar distance;
        unsigned boundaryFaceIndex;
        double faceArea;
    };
 
    struct BoundaryInfo
    {
        unsigned int cell;
        int dir;
        unsigned int bfIndex;
        BoundaryConditionData bcdata;
    };
    std::vector<BoundaryInfo> boundaryInfo_;
 
    struct FullDomain
    {
        std::vector<int> cells;
        std::vector<bool> interior;
    };
    FullDomain fullDomain_;
};  // class TpsaLinearizer
 
}  // namespace Opm
 
#endif