tpsamodel.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_MODEL_HPP
#define TPSA_MODEL_HPP
 
#include <dune/grid/common/gridenums.hh>
 
#include <dune/common/dynmatrix.hh>
#include <dune/common/dynvector.hh>
 
#include <opm/grid/utility/ElementChunks.hpp>
 
#include <opm/common/utility/SymmTensor.hpp>
#include <opm/common/utility/VoigtArray.hpp>
#include <opm/material/common/MathToolbox.hpp>
 
#include <opm/models/parallel/threadmanager.hpp>
#include <opm/models/tpsa/tpsabaseproperties.hpp>
#include <opm/models/utils/linearleastsquares.hpp>
#include <opm/models/utils/propertysystem.hh>
 
#include <array>
#include <memory>
#include <unordered_map>
 
 
namespace Opm {
 
template <class TypeTag>
class TpsaModel
{
    using DofMapper = GetPropType<TypeTag, Properties::DofMapper>;
    using Evaluation = GetPropType<TypeTag, Properties::EvaluationTPSA>;
    using GlobalEqVector = GetPropType<TypeTag, Properties::GlobalEqVectorTPSA>;
    using GridView = GetPropType<TypeTag, Properties::GridView>;
    using Indices = GetPropType<TypeTag, Properties::IndicesTPSA>;
    using Linearizer = GetPropType<TypeTag, Properties::LinearizerTPSA>;
    using NewtonMethod = GetPropType<TypeTag, Properties::NewtonMethodTPSA>;
    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariablesTPSA>;
    using Simulator = GetPropType<TypeTag, Properties::Simulator>;
    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVectorTPSA>;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
 
    enum { dimWorld = GridView::dimensionworld };
    enum { historySize = getPropValue<TypeTag, Properties::SolutionHistorySizeTPSA>() };
    enum { numEq = getPropValue<TypeTag, Properties::NumEqTPSA>() };
 
    enum { disp0Idx = Indices::disp0Idx };
    enum { rot0Idx = Indices::rot0Idx };
    enum { solidPres0Idx = Indices::solidPres0Idx };
 
    using MaterialState = MaterialStateTPSA<Evaluation>;
 
    using DimVector = Dune::FieldVector<Scalar, dimWorld>;
    using SymTensor = SymmTensor<Scalar>;
    using PotForceVector = Dune::BlockVector<Scalar>;
 
public:
    class TpsaBlockVectorWrapper
    {
    protected:
        SolutionVector blockVector_;
    public:
        TpsaBlockVectorWrapper(const std::string&, const std::size_t size)
            : blockVector_(size)
        {}
 
        TpsaBlockVectorWrapper() = default;
 
        static TpsaBlockVectorWrapper serializationTestObject()
        {
            TpsaBlockVectorWrapper result("dummy", 3);
            result.blockVector_[0] = 1.0;
            result.blockVector_[1] = 2.0;
            result.blockVector_[2] = 3.0;
 
            return result;
        }
 
        SolutionVector& blockVector()
        { return blockVector_; }
 
        const SolutionVector& blockVector() const
        { return blockVector_; }
 
        bool operator==(const TpsaBlockVectorWrapper& wrapper) const
        {
            return std::equal(this->blockVector_.begin(), this->blockVector_.end(),
                              wrapper.blockVector_.begin(), wrapper.blockVector_.end());
        }
 
        template<class Serializer>
        void serializeOp(Serializer& serializer)
        {
            serializer(blockVector_);
        }
    };
 
    // ///
    // Public functions
    // ///
    explicit TpsaModel(Simulator& simulator)
        : linearizer_(std::make_unique<Linearizer>())
        , newtonMethod_(simulator)
        , simulator_(simulator)
        , element_chunks_(simulator.gridView(), Dune::Partitions::all, ThreadManager::maxThreads())
    {
        // Initialize equation weights to 1.0
        eqWeights_.resize(numEq, 1.0);
 
        // Initialize historic solution vectors
        // OBS: need at least history size = 2, due to time-derivative of solid-pressure in Flow coupling term
        const std::size_t numDof = simulator_.model().numGridDof();
        for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
            solution_[timeIdx] = std::make_unique<TpsaBlockVectorWrapper>("solution", numDof);
        }
 
        // Initialize potential force vectors
        mechPotPresForce_.resize(numDof);
        mechPotTempForce_.resize(numDof);
    }
 
    void finishInit()
    {
        // Initialize the linearizer
        linearizer_->init(simulator_);
 
        // Resize material state vector
        resizeMaterialState_();
    }
 
    static void registerParameters()
    {
        // Newton method parameters
        NewtonMethod::registerParameters();
    }
 
    void prepareTPSA()
    {
        // Update historic solution
        solution(/*timeIdx=*/1) = solution(/*timeIdx=*/0);
    }
 
    void syncOverlap()
    {
        // Syncronize the solution on the ghost and overlap elements
        using GhostSyncHandle = GridCommHandleGhostSync<PrimaryVariables,
                                                        SolutionVector,
                                                        DofMapper,
                                                        /*commCodim=*/0>;
 
        auto ghostSync = GhostSyncHandle(solution(/*timeIdx=*/0),
                                         simulator_.model().dofMapper());
 
        simulator_.gridView().communicate(ghostSync,
                                          Dune::InteriorBorder_All_Interface,
                                          Dune::ForwardCommunication);
    }
 
    void updateMaterialState(const unsigned /*timeIdx*/)
    {
        // Loop over all elements chuncks and update material state from current solution
        const auto& elementMapper = simulator_.model().elementMapper();
#ifdef _OPENMP
#pragma omp parallel for
#endif
        for (const auto& chunk : element_chunks_) {
            for (const auto& elem : chunk) {
                const unsigned globalIdx = elementMapper.index(elem);
                auto& currSol = solution(/*timeIdx=*/0)[globalIdx];
                setMaterialState_(globalIdx, /*timeIdx=*/0, currSol);
            }
        }
    }
 
    // ///
    // Public get and set functions
    // ///
    const Linearizer& linearizer() const
    {
        return *linearizer_;
    }
 
    Linearizer& linearizer()
    {
        return *linearizer_;
    }
 
    const NewtonMethod& newtonMethod() const
    {
        return newtonMethod_;
    }
 
    NewtonMethod& newtonMethod()
    {
        return newtonMethod_;
    }
 
    const SolutionVector& solution(unsigned timeIdx) const
    {
        return solution_[timeIdx]->blockVector();
    }
 
    SolutionVector& solution(unsigned timeIdx)
    {
        return solution_[timeIdx]->blockVector();
    }
 
    std::size_t numGridDof() const
    {
        return simulator_.model().numGridDof();
    }
 
    std::size_t numTotalDof() const
    {
        return numGridDof() + numAuxiliaryDof();
    }
 
    Scalar dofTotalVolume(unsigned globalIdx) const
    {
        return simulator_.model().dofTotalVolume(globalIdx);
    }
 
    Scalar eqWeight(unsigned /*dofIdx*/, unsigned eqIdx) const
    {
        return eqWeights_[eqIdx];
    }
 
    void setEqWeight(unsigned eqIdx, Scalar value)
    {
        eqWeights_[eqIdx] = value;
    }
 
    std::size_t numAuxiliaryModules() const
    {
        return 0;
    }
 
    std::size_t numAuxiliaryDof() const
    {
        return 0;
    }
 
    const MaterialState& materialState(const unsigned globalIdx, unsigned /*timeIdx*/) const
    {
        return materialState_[globalIdx];
    }
 
    DimVector disp(const unsigned globalIdx, const bool /*with_fracture*/) const
    {
        DimVector d;
        for (std::size_t i = 0; i < 3; ++i) {
            d[i] = decay<Scalar>(materialState_[globalIdx].displacement(i));
        }
        return d;
    }
 
    DimVector rotation(const unsigned globalIdx) const
    {
        DimVector rot;
        for (std::size_t i = 0; i < 3; ++i) {
            rot[i] = decay<Scalar>(materialState_[globalIdx].rotation(i));
        }
        return rot;
    }
 
    Scalar solidPressure(const unsigned globalIdx) const
    {
        return decay<Scalar>(materialState_[globalIdx].solidPressure());
    }
 
    SymTensor delstress(const unsigned globalIdx) const
    {
        return stress(globalIdx, false);
    }
 
    SymTensor fractureStress(const unsigned /*globalIdx*/) const
    {
        SymTensor val;
        return val;
    }
 
    SymTensor linstress(const unsigned globalIdx) const
    {
        // Subtract the potential forces from pressure and temperature from the main diagonal of
        // the total stress
        SymTensor linStressTensor = stress(globalIdx, false);
        const auto potForce = mechPotentialForce(globalIdx);
        for (const auto& dirIdx : SymTensor::diag_indices) {
            linStressTensor[dirIdx] -= potForce;
        }
        return -1.0 * linStressTensor;
    }
 
    SymTensor stress(const unsigned globalIdx, const bool /*with_fracture*/) const
    {
        SymTensor stressOutput;
        const auto& stressInfo = linearizer_->getStressInfo();
        if (!stressInfo.empty()) {
            const auto& stressInfoGlobI = stressInfo[globalIdx];
 
            // Setup least-squares matrix of face normals and right-hand side of traction vectors
            // per faces. Matrix shape = 3 rows per face x 6 columns for each (symmetric) stress
            // tensor component.
            std::unordered_multimap<int, std::size_t> faceIdMap;
            std::size_t mapSize = 0;
            for (std::size_t sInfoIdx = 0; sInfoIdx < stressInfoGlobI.size(); ++sInfoIdx) {
                const auto faceId  = stressInfoGlobI[sInfoIdx].faceId;
 
                // OBS: NNC faces are not added to least squares system, since TPSA does not handle
                // NNC connections, like numerical aquifers, yet!
                if (faceId < 0 || stressInfoGlobI[sInfoIdx].faceArea == 0.0) {
                    continue;
                }
                if (!faceIdMap.contains(faceId)) {
                    ++mapSize;
                }
                faceIdMap.insert({faceId, sInfoIdx});
            }
 
            // If no valid faces exist for cell, return zero SymTensor
            if (mapSize == 0) {
                return stressOutput;
            }
 
            Dune::DynamicMatrix<Scalar> mat(3 * mapSize, 6);
            Dune::DynamicVector<Scalar> rhs(3 * mapSize);
            std::size_t rowIdx = 0;
            for (auto it = faceIdMap.begin(); it != faceIdMap.end();) {
                auto [first, last] = faceIdMap.equal_range(it->first);
 
                // Loop over possible multiple of the same face
                Scalar sumFaceArea = 0.0;
                for (auto inner = first; inner != last; ++inner) {
                    const auto sInfoIdx = inner->second;
                    const auto& faceStressInfo = stressInfoGlobI[sInfoIdx];
 
                    // One normal per face
                    if (inner == first) {
                        const auto& fNormal = faceStressInfo.faceNormal;
                        mat[3*rowIdx][0] = fNormal[0];
                        mat[3*rowIdx][3] = fNormal[1];
                        mat[3*rowIdx][4] = fNormal[2];
 
                        mat[3*rowIdx + 1][1] = fNormal[1];
                        mat[3*rowIdx + 1][3] = fNormal[0];
                        mat[3*rowIdx + 1][5] = fNormal[2];
 
                        mat[3*rowIdx + 2][2] = fNormal[2];
                        mat[3*rowIdx + 2][4] = fNormal[0];
                        mat[3*rowIdx + 2][5] = fNormal[1];
                    }
 
                    // Add traction vectors for same faces
                    const auto& fTraction = faceStressInfo.traction;
                    rhs[3*rowIdx] += fTraction[0];
                    rhs[3*rowIdx + 1] += fTraction[1];
                    rhs[3*rowIdx + 2] += fTraction[2];
 
                    // Sum face area
                    sumFaceArea += faceStressInfo.faceArea;
                }
                // Divide rhs for (unique) face by sum of face areas
                rhs[3*rowIdx] /= sumFaceArea;
                rhs[3*rowIdx + 1] /= sumFaceArea;
                rhs[3*rowIdx + 2] /= sumFaceArea;
 
                // Move to next unique face
                ++rowIdx;
                it = last;
            }
 
            // Use least squares solver for underdetermined system
            LinearLeastSquares<Scalar> lsq(mat, rhs);
            lsq.solve();
 
            // Reconstructed stress tensor at cell center
            const auto& stress = lsq.x();
            stressOutput[VoigtIndex::XX] = stress[0]; // XX
            stressOutput[VoigtIndex::YY] = stress[1]; // YY
            stressOutput[VoigtIndex::ZZ] = stress[2]; // ZZ
            stressOutput[VoigtIndex::YZ] = stress[5]; // YZ
            stressOutput[VoigtIndex::XZ] = stress[4]; // XZ
            stressOutput[VoigtIndex::XY] = stress[3]; // XY
        }
        return stressOutput;
    }
 
    std::array<DimVector, 6> traction(const unsigned globalIdx) const
    {
        std::array<DimVector, 6> tractionOutput{};
        const auto& stressInfo = linearizer_->getStressInfo();
        if (stressInfo.empty()) {
            return tractionOutput;
        }
 
        // Average traction (force per area) over all face entries sharing the same face
        // direction, consistent with the face-area weighting used in stress().
        std::array<Scalar, 6> sumFaceArea{};
        for (const auto& faceStressInfo : stressInfo[globalIdx]) {
            const auto faceId = faceStressInfo.faceId;
 
            // OBS: NNC faces are not added, since TPSA does not handle NNC connections yet!
            if (faceId < 0 || faceStressInfo.faceArea == 0.0) {
                continue;
            }
 
            tractionOutput[faceId] += faceStressInfo.traction;
            sumFaceArea[faceId] += faceStressInfo.faceArea;
        }
 
        for (std::size_t faceId = 0; faceId < 6; ++faceId) {
            if (sumFaceArea[faceId] > 0.0) {
                tractionOutput[faceId] /= sumFaceArea[faceId];
            }
        }
 
        return tractionOutput;
    }
 
    SymTensor strain(const unsigned globalIdx, const bool /*with_fracture*/) const
    {
        // Deviatoric stress
        auto stressDev = this->linstress(globalIdx);
        Scalar traceStress = 1.0 / 3.0 * stressDev.trace();
        stressDev[VoigtIndex::XX] -= traceStress;
        stressDev[VoigtIndex::YY] -= traceStress;
        stressDev[VoigtIndex::ZZ] -= traceStress;
 
        // Deviatoric strain
        auto& problem = simulator_.problem();
        const auto sMod = problem.shearModulus(globalIdx);
        SymTensor strainDev = 1.0 / (2.0 * sMod) * stressDev;
 
        // Volumetric strain
        const auto lameParam = problem.lame(globalIdx);
        Scalar strainVolTerm = traceStress / (3.0 * lameParam + 2 * sMod);
 
        SymTensor strainVol;
        strainVol[VoigtIndex::XX] = strainVolTerm;
        strainVol[VoigtIndex::YY] = strainVolTerm;
        strainVol[VoigtIndex::ZZ] = strainVolTerm;
 
        // Total
        SymTensor strainOutput = strainDev + strainVol;
        return strainOutput;
    }
 
    Scalar mechPotentialForce(unsigned globalIdx) const
    {
        return mechPotentialPressForce(globalIdx) + mechPotentialTempForce(globalIdx);
    }
 
    Scalar mechPotentialPressForce(unsigned globalIdx) const
    {
        return mechPotPresForce_[globalIdx];
    }
 
    void setMechPotentialPressForce(unsigned globalIdx, Scalar val)
    {
        mechPotPresForce_[globalIdx] = val;
    }
 
    Scalar mechPotentialTempForce(unsigned globalIdx) const
    {
        return mechPotTempForce_[globalIdx];
    }
 
    void setMechPotentialTempForce(unsigned globalIdx, Scalar val)
    {
        mechPotTempForce_[globalIdx] = val;
    }
 
protected:
    // ///
    // Protected functions
    // ///
    void resizeMaterialState_()
    {
        const std::size_t numDof = simulator_.model().numGridDof();
        materialState_.resize(numDof);
    }
 
private:
    // ///
    // Private functions
    // ///
    void setMaterialState_(const unsigned globalIdx, const unsigned /*timeIdx*/, PrimaryVariables& values)
    {
        auto& dofMaterialState = materialState_[globalIdx];
            for (unsigned dirIdx = 0; dirIdx < 3; ++dirIdx) {
                dofMaterialState.setDisplacement(dirIdx, values.makeEvaluation(disp0Idx + dirIdx, 0));
                dofMaterialState.setRotation(dirIdx, values.makeEvaluation(rot0Idx + dirIdx, 0));
            }
            dofMaterialState.setSolidPressure(values.makeEvaluation(solidPres0Idx, 0));
    }
 
    std::unique_ptr<Linearizer> linearizer_;
    NewtonMethod newtonMethod_;
    Simulator& simulator_;
    ElementChunks<GridView, Dune::Partitions::All> element_chunks_;
 
    std::array<std::unique_ptr<TpsaBlockVectorWrapper>, historySize> solution_;
    std::vector<Scalar> eqWeights_;
    std::vector<MaterialState> materialState_;
    PotForceVector mechPotPresForce_;
    PotForceVector mechPotTempForce_;
};  // class TpsaModel
 
}  // namespace Opm
 
 
#endif