mumpssolver.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:
//
// SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
#ifndef DUMUX_LINEAR_MUMPS_SOLVER_HH
#define DUMUX_LINEAR_MUMPS_SOLVER_HH
 
#include <string>
#include <vector>
#include <tuple>
#include <numeric>
#include <limits>
#include <algorithm>
#include <functional>
#include <chrono>
#include <iostream>
#include <cstring>
#include <cmath>
#include <memory>
#include <optional>
 
#include <dune/common/exceptions.hh>
#include <dune/common/indices.hh>
#include <dune/common/hybridutilities.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/istl/bvector.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/solver.hh>
#include <dune/istl/solvercategory.hh>
 
#include <dumux/common/parameters.hh>
#include <dumux/linear/solver.hh>
#include <dumux/linear/solvercategory.hh>
#include <dumux/linear/parallelhelpers.hh>
#include <dumux/linear/linearsolvertraits.hh>
#include <dumux/linear/matrixconverter.hh>
#include <dumux/common/typetraits/matrix.hh>
#include <dumux/common/typetraits/vector.hh>
#include <dumux/parallel/vectorcommdatahandle.hh>
 
#if DUMUX_HAVE_MUMPS
 
#include <mpi.h>
#include <dmumps_c.h>
 
namespace Dumux::Detail {
 
// Provides GridView and DofMapper from LSTraits if the traits belong to a grid,
// or inert dummy types otherwise (SeqLinearSolverTraits).
template<class LSTraits>
struct MumpsGridTypes {
    static constexpr bool hasGrid = false;
    struct GridView {};
    struct DofMapper { std::size_t size() const { return 0; } };
};
 
template<class LSTraits>
requires requires { typename LSTraits::GridView; typename LSTraits::DofMapper; }
struct MumpsGridTypes<LSTraits> {
    static constexpr bool hasGrid = true;
    using GridView = typename LSTraits::GridView;
    using DofMapper = typename LSTraits::DofMapper;
};
 
// Dereference pointer-like objects (e.g. std::shared_ptr<GridGeometry>) and pass through
// references/values unchanged. Used so the multidomain constructor accepts a tuple of either
// (shared) pointers to grid geometries or the grid geometries themselves.
template<class T>
decltype(auto) mumpsDeref(const T& t)
{
    if constexpr (requires (const T& x) { *x; })
        return (*t);
    else
        return (t);
}
 
} // end namespace Dumux::Detail
 
namespace Dumux {

template<class LSTraits, class LATraits>
class DirectSolverMumps : public LinearSolver
{
    using Matrix = typename LATraits::Matrix;
    using XVector = typename LATraits::Vector;
    using BVector = typename LATraits::Vector;
 
    using Scalar = double;
    // Global scalar-DOF index. MUMPS itself uses 32-bit MUMPS_INT for the matrix
    // entries; we keep a wider type internally and cast (with a check) at the boundary.
    using GlobalIndex = long long;
 
    using GridTypes_ = Detail::MumpsGridTypes<LSTraits>;
    using GridView = typename GridTypes_::GridView;
    using DofMapper = typename GridTypes_::DofMapper;
    static constexpr bool hasGrid = GridTypes_::hasGrid;
 
    // True when XVector is a MultiTypeBlockVector (multidomain assembler).
    static constexpr bool isMultiType = isMultiTypeBlockVector<XVector>::value;
    // Scalar DOFs per block DOF. MultiType is always converted to scalar (blockSize=1).
    static constexpr std::size_t blockSize = []() constexpr -> std::size_t {
        if constexpr (isMultiTypeBlockVector<XVector>::value) return 1;
        else return XVector::block_type::size();
    }();
 
    // Flattened scalar block vector type produced by VectorConverter (1 scalar per block).
    // The per-subdomain RHS/ghost closures operate on slices of a vector of this type.
    using ScalarBlockVector = Dune::BlockVector<Dune::FieldVector<Scalar, 1>>;
 
    // Per-subdomain parallel bookkeeping (MultiType only). One entry per subdomain, in the
    // same order as the MultiTypeBlockVector subvectors / the flattened scalar layout.
    struct SubdomainInfo
    {
        std::size_t scalarBS = 1; // scalar DOFs per block DOF in this subdomain
        std::size_t nBlockLocal = 0; // local block DOFs of this subdomain on this rank
        std::size_t scalarOffset = 0; // start index of this subdomain in the flattened scalar vector
        bool rowsComplete = true; // owner border rows already complete after local assembly?
        // Sum partial border-DOF RHS contributions onto the owner (non-overlapping incomplete rows).
        // Empty when not needed (overlapping or complete rows). Operates on its scalar slice.
        std::function<void(ScalarBlockVector&)> rhsPreprocess;
        // Broadcast the owner's solution value to all copies (zero non-owned beforehand).
        // Operates on its scalar slice.
        std::function<void(ScalarBlockVector&)> ghostBroadcast;
    };
 
public:
    struct SolverResult : public Dune::InverseOperatorResult
    {
        SolverResult() = default;
        SolverResult(const Dune::InverseOperatorResult& o) : InverseOperatorResult(o) {}
        operator bool() const { return this->converged; }
    };

    explicit DirectSolverMumps(const std::string& paramGroup = "")
    : LinearSolver(paramGroup)
    , isParallel_(false)
    {
        if constexpr (hasGrid)
            if (Dune::MPIHelper::getCommunication().size() > 1)
                DUNE_THROW(Dune::InvalidStateException, "Using sequential constructor for parallel run. Use signature with gridView and dofMapper!");
 
        mpiComm_ = MPI_COMM_SELF;
        initMumps_();
        readSolverParams_();
    }

    DirectSolverMumps(const GridView& gridView,
                      const DofMapper& dofMapper,
                      const std::string& paramGroup = "")
    requires (hasGrid)
    : LinearSolver(paramGroup)
    , gridView_(gridView)
    , mapper_(&dofMapper)
    , isParallel_(Detail::solverCategory<LSTraits>(gridView) != Dune::SolverCategory::sequential)
    {
        if constexpr (LSTraits::canCommunicate)
            mpiComm_ = static_cast<MPI_Comm>(gridView.comm());
        else
            mpiComm_ = MPI_COMM_SELF;
        initMumps_();
        readSolverParams_();
 
        buildGlobalDofIndices_(gridView, dofMapper);
    }

    template<class GGTuple>
    explicit DirectSolverMumps(const GGTuple& ggTuple, const std::string& paramGroup = "")
    requires (isMultiType && LSTraits::canCommunicate && requires { std::tuple_size<GGTuple>::value; })
    : LinearSolver(paramGroup)
    , gridView_(Detail::mumpsDeref(std::get<0>(ggTuple)).gridView())
    , isParallel_(Detail::mumpsDeref(std::get<0>(ggTuple)).gridView().comm().size() > 1)
    {
        constexpr std::size_t numSubDomains = std::tuple_size_v<GGTuple>;
        static_assert(numSubDomains == std::tuple_size_v<XVector>,
                      "Number of grid geometries must match the number of MultiTypeBlockVector subdomains");
 
        mpiComm_ = static_cast<MPI_Comm>(Detail::mumpsDeref(std::get<0>(ggTuple)).gridView().comm());
        initMumps_();
        readSolverParams_();
 
        buildMultiTypeMaps_(ggTuple);
    }
 
    ~DirectSolverMumps()
    {
        if (initialized_)
        {
            id_.job = -2; // terminate the MUMPS instance
            dmumps_c(&id_);
        }
    }
 
    // Non-copyable: owns a MUMPS instance and raw pointers into our buffers.
    DirectSolverMumps(const DirectSolverMumps&) = delete;
    DirectSolverMumps& operator=(const DirectSolverMumps&) = delete;

    void updateAfterGridAdaption(const GridView& gridView, const DofMapper& dofMapper)
    requires (hasGrid)
    {
        gridView_ = gridView;
        mapper_ = &dofMapper;
        buildGlobalDofIndices_(gridView, dofMapper);
        analyzed_ = false;
        factorized_ = false;
    }

    template<class GGTuple>
    void updateAfterGridAdaption(const GGTuple& ggTuple)
    requires (isMultiType && LSTraits::canCommunicate && requires { std::tuple_size<GGTuple>::value; })
    {
        gridView_ = Detail::mumpsDeref(std::get<0>(ggTuple)).gridView();
        buildMultiTypeMaps_(ggTuple);
        analyzed_ = false;
        factorized_ = false;
    }

    SolverResult solve(const Matrix& A, XVector& x, const BVector& b)
    {
        setMatrix(A);
        return solve(x, b);
    }

    void setMatrix(const Matrix& A)
    {
        factorized_ = false;
        prepareMaps_(A);
        factorize_(A);
    }

    void setMatrix(std::shared_ptr<Matrix> A)
    { setMatrix(*A); }

    SolverResult solve(XVector& x, const BVector& b)
    { return solveColumns_({&x}, {&b}); }

    SolverResult solve(std::vector<XVector>& x, const std::vector<BVector>& b)
    {
        std::vector<XVector*> xColumns;
        std::vector<const BVector*> bColumns;
        for (auto& xj : x)
            xColumns.push_back(&xj);
        for (const auto& bj : b)
            bColumns.push_back(&bj);
        return solveColumns_(xColumns, bColumns);
    }
 
    const Dune::InverseOperatorResult& result() const
    { return result_; }
 
    Scalar norm(const XVector& x) const
    {
        if constexpr (isMultiType)
        {
            if (!isParallel_ || localToGlobal_.empty())
            {
                auto y = VectorConverter<XVector>::multiTypeToBlockVector(x);
                return y.two_norm();
            }
            // Parallel: sum border-DOF partial contributions before the owned-DOF global norm.
            auto v = VectorConverter<XVector>::multiTypeToBlockVector(x);
            for (const auto& sub : subs_)
                if (sub.rhsPreprocess)
                    sub.rhsPreprocess(v);
 
            double localSumSq = 0.0;
            for (std::size_t r = 0; r < v.size(); ++r)
                if (isOwned_[r])
                    localSumSq += v[r][0] * v[r][0];
 
            double globalSumSq = 0.0;
            MPI_Allreduce(&localSumSq, &globalSumSq, 1, MPI_DOUBLE, MPI_SUM, mpiComm_);
            using std::sqrt;
            return sqrt(globalSumSq);
        }
        else
        {
            if (!isParallel_ || localToGlobal_.empty())
                return x.two_norm();
 
            // Parallel single-type: true global residual norm over owned DOFs.
            if constexpr (LSTraits::canCommunicate)
            {
                if (LSTraits::isNonOverlapping(*gridView_))
                {
                    auto xcopy = x;
                    if constexpr (requires { LSTraits::dofCodims; })
                    {
                        MultiCodimParallelVectorHelper<GridView, DofMapper> vecHelper(*gridView_, *mapper_);
                        vecHelper.makeNonOverlappingConsistent(xcopy, LSTraits::dofCodims);
                    }
                    else
                    {
                        ParallelVectorHelper<GridView, DofMapper, LSTraits::dofCodim> vecHelper(*gridView_, *mapper_);
                        vecHelper.makeNonOverlappingConsistent(xcopy);
                    }
                    return ownedGlobalNorm_(xcopy);
                }
            }
            return ownedGlobalNorm_(x);
        }
    }
 
    std::string name() const { return "DuMux direct MUMPS solver"; }
 
private:
    // Accessor for MUMPS control parameters using the manual's 1-based ICNTL(i) numbering.
    int& icntl_(int i) { return id_.icntl[i - 1]; }
    const int& infog_(int i) const { return id_.infog[i - 1]; }

    void initMumps_()
    {
        // Zero the whole control struct so distributed-RHS/solution fields MUMPS may read
        // (nz_rhs, nloc_rhs, pointers, ...) start well-defined rather than holding stack garbage.
        std::memset(&id_, 0, sizeof(id_));
        id_.comm_fortran = static_cast<int>(MPI_Comm_c2f(mpiComm_));
        id_.par = 1;
        id_.sym = 0;
        id_.job = -1;
        dmumps_c(&id_);
        initialized_ = true;
        checkError_("initialization");
 
        icntl_(1) = -1; // suppress error output stream
        icntl_(2) = -1; // suppress diagnostic/warning output stream
        icntl_(3) = -1; // suppress global information output stream
        icntl_(4) =  1; // verbosity: errors only
        icntl_(5) =  0; // assembled matrix (not elemental)
        icntl_(18) = 3; // distributed assembled input: each rank supplies its local triplets
        icntl_(20) = 0; // dense, centralized RHS on the host
        icntl_(21) = 0; // centralized solution on the host
    }

    void readSolverParams_()
    {
        const int ordering = getParamFromGroup<int>(this->paramGroup(), "LinearSolver.MumpsOrdering", -1);
        if (ordering >= 0)
            icntl_(7) = ordering;
 
        // ICNTL(14): percentage increase of the estimated workspace, needed if numerical pivoting
        // creates more fill-in than the symbolic analysis predicted
        const int workspaceIncrease = getParamFromGroup<int>(this->paramGroup(), "LinearSolver.MumpsWorkspaceIncrease", -1);
        if (workspaceIncrease >= 0)
            icntl_(14) = workspaceIncrease;
 
        // ICNTL(28): symbolic analysis mode. 0=automatic, 1=sequential (default), 2=parallel.
        // Parallel analysis requires MUMPS built with a parallel ordering tool (PT-SCOTCH or
        // ParMETIS); without one MUMPS silently falls back to sequential analysis. ICNTL(29)
        // selects that tool (0=auto, 1=PT-SCOTCH, 2=ParMETIS).
        const int parAnalysis = getParamFromGroup<int>(this->paramGroup(), "LinearSolver.MumpsParallelAnalysis", -1);
        if (parAnalysis >= 0)
            icntl_(28) = parAnalysis;
        const int parOrderingTool = getParamFromGroup<int>(this->paramGroup(), "LinearSolver.MumpsParallelOrderingTool", -1);
        if (parOrderingTool >= 0)
            icntl_(29) = parOrderingTool;
    }

    void buildGlobalDofIndices_(const GridView& gridView, const DofMapper& mapper)
    {
        const std::size_t nLocal = mapper.size();
        localToGlobal_.assign(nLocal, std::numeric_limits<GlobalIndex>::max());
        isOwned_.assign(nLocal, false);
        isGhost_.assign(nLocal, false);
 
        if constexpr (LSTraits::canCommunicate)
        {
            if (isParallel_)
            {
                ParallelISTLHelper<LSTraits> pHelper(gridView, mapper);
                for (std::size_t i = 0; i < nLocal; ++i)
                {
                    isOwned_[i] = pHelper.isOwned(i);
                    isGhost_[i] = pHelper.isGhost(i);
                }
 
                GlobalIndex numOwned = static_cast<GlobalIndex>(
                    std::count(isOwned_.begin(), isOwned_.end(), true));
                const int commSize = gridView.comm().size();
                const int myRank = gridView.comm().rank();
                std::vector<GlobalIndex> counts(commSize);
                gridView.comm().allgather(&numOwned, 1, counts.data());
 
                GlobalIndex offset = 0;
                for (int p = 0; p < myRank; ++p)
                    offset += counts[p];
                numBlocksGlobal_ = 0;
                for (auto c : counts) numBlocksGlobal_ += static_cast<std::size_t>(c);
 
                for (std::size_t i = 0; i < nLocal; ++i)
                    if (isOwned_[i])
                        localToGlobal_[i] = offset++;
 
                // Propagate global IDs from owners to ghost DOFs via min communication.
                if constexpr (requires { LSTraits::dofCodims; }) {
                    MultiCodimVectorCommDataHandleMin<DofMapper, std::vector<GlobalIndex>, GridView::dimension, GlobalIndex>
                        handle(mapper, localToGlobal_, LSTraits::dofCodims);
                    gridView.communicate(handle, Dune::InteriorBorder_All_Interface, Dune::ForwardCommunication);
                } else {
                    VectorCommDataHandleMin<DofMapper, std::vector<GlobalIndex>, LSTraits::dofCodim, GlobalIndex>
                        handle(mapper, localToGlobal_);
                    gridView.communicate(handle, Dune::InteriorBorder_All_Interface, Dune::ForwardCommunication);
                }
                return;
            }
        }
 
        // Sequential: all DOFs owned with consecutive IDs.
        numBlocksGlobal_ = nLocal;
        std::fill(isOwned_.begin(), isOwned_.end(), true);
        for (std::size_t i = 0; i < nLocal; ++i)
            localToGlobal_[i] = static_cast<GlobalIndex>(i);
    }
 
    // Obtain a subdomain's DOF mapper discretization-generically: most grid geometries expose
    // dofMapper(); FCStaggered only provides one through its linear-solver traits (a lightweight
    // by-value wrapper around the grid view).
    template<class SubTraits, class SubGG>
    static decltype(auto) subDofMapper_(const SubGG& gg)
    {
        if constexpr (requires { gg.dofMapper(); })
            return (gg.dofMapper());
        else
            return SubTraits::dofMapper(gg);
    }
 
    // True if a subdomain's DOF mapper can drive the generic codim-based grid communication
    // (it must expose indices(entity)); the minimal FCStaggered mapper wrapper cannot, so
    // FCStaggered subdomains are only supported sequentially.
    template<class SubTraits, class GV>
    static constexpr bool subCanComm_()
    {
        using Mapper = typename SubTraits::DofMapper;
        using Element = typename GV::template Codim<0>::Entity;
        return SubTraits::canCommunicate
            && requires (const Mapper& m, const Element& e) { m.indices(e); };
    }
 
    // Build global block indices for one subdomain (MultiType). Returns the subdomain's
    // total number of global block DOFs. gids[i] = global block ID (offset to be globally unique).
    // CanComm gates the parallel branch so it is not instantiated for mappers that cannot
    // communicate (e.g. FCStaggered); such subdomains fall back to sequential numbering.
    template<class SubLSTraits, bool CanComm, class GV, class Mapper>
    GlobalIndex buildSubdomainGIDs_(const GV& gridView, const Mapper& mapper,
                                    GlobalIndex globalOffset,
                                    std::vector<GlobalIndex>& gids,
                                    std::vector<bool>& owned,
                                    std::vector<bool>& ghost)
    {
        const std::size_t nLocal = mapper.size();
        gids.assign(nLocal, std::numeric_limits<GlobalIndex>::max());
        owned.assign(nLocal, false);
        ghost.assign(nLocal, false);
 
        if constexpr (CanComm)
        {
            if (isParallel_)
            {
                ParallelISTLHelper<SubLSTraits> pHelper(gridView, mapper);
                for (std::size_t i = 0; i < nLocal; ++i)
                {
                    owned[i] = pHelper.isOwned(i);
                    ghost[i] = pHelper.isGhost(i);
                }
 
                GlobalIndex numOwned = static_cast<GlobalIndex>(
                    std::count(owned.begin(), owned.end(), true));
                const int commSize = gridView.comm().size();
                const int myRank = gridView.comm().rank();
                std::vector<GlobalIndex> counts(commSize);
                gridView.comm().allgather(&numOwned, 1, counts.data());
 
                GlobalIndex localOffset = globalOffset;
                for (int p = 0; p < myRank; ++p)
                    localOffset += counts[p];
 
                GlobalIndex nGlobal = 0;
                for (auto c : counts) nGlobal += c;
 
                for (std::size_t i = 0; i < nLocal; ++i)
                    if (owned[i])
                        gids[i] = localOffset++;
 
                if constexpr (requires { SubLSTraits::dofCodims; }) {
                    MultiCodimVectorCommDataHandleMin<Mapper, std::vector<GlobalIndex>, GV::dimension, GlobalIndex>
                        handle(mapper, gids, SubLSTraits::dofCodims);
                    gridView.communicate(handle, Dune::InteriorBorder_All_Interface, Dune::ForwardCommunication);
                } else {
                    constexpr int codim = SubLSTraits::dofCodim;
                    VectorCommDataHandleMin<Mapper, std::vector<GlobalIndex>, codim, GlobalIndex>
                        handle(mapper, gids);
                    gridView.communicate(handle, Dune::InteriorBorder_All_Interface, Dune::ForwardCommunication);
                }
                return nGlobal;
            }
        }
 
        // Sequential subdomain: all DOFs owned with consecutive IDs.
        std::fill(owned.begin(), owned.end(), true);
        for (std::size_t i = 0; i < nLocal; ++i)
            gids[i] = globalOffset + static_cast<GlobalIndex>(i);
        return static_cast<GlobalIndex>(nLocal);
    }

    template<class GGTuple>
    void buildMultiTypeMaps_(const GGTuple& ggTuple)
    {
        constexpr std::size_t numSubDomains = std::tuple_size_v<GGTuple>;
        localToGlobal_.clear();
        isOwned_.clear();
        isGhost_.clear();
        rowComplete_.clear();
        subs_.assign(numSubDomains, SubdomainInfo{});
 
        // Running offsets while we append the subdomains' scalar DOFs in order.
        GlobalIndex scalarGlobalOffset = 0; // global scalar index where the current subdomain starts
        std::size_t scalarLocalOffset = 0; // local scalar index where the current subdomain starts
 
        using Dune::Hybrid::forEach;
        forEach(std::make_index_sequence<numSubDomains>{}, [&](auto I)
        {
            constexpr std::size_t i = decltype(I)::value;
            const auto& gg = Detail::mumpsDeref(std::get<i>(ggTuple));
            using SubGG = std::decay_t<decltype(gg)>;
            using SubTraits = LinearSolverTraits<SubGG>;
            using SubBlockVec = std::tuple_element_t<i, XVector>;
            constexpr std::size_t scalarBS = SubBlockVec::block_type::size();
            using GV = typename SubGG::GridView;
            constexpr bool canComm = subCanComm_<SubTraits, GV>();
 
            // Parallel runs require every subdomain's DOF mapper to support the generic grid
            // communication; bail out with a clear message otherwise (e.g. FCStaggered momentum).
            if (isParallel_ && !canComm)
                DUNE_THROW(Dune::NotImplemented,
                    "DirectSolverMumps: parallel multidomain solve is not supported for subdomain "
                    << i << " — its discretization's DOF mapper (e.g. FaceCenteredStaggered) cannot "
                    "drive the generic grid communication. This combination works sequentially; for "
                    "parallel runs use a momentum discretization with a full mapper (CVFE/FCDiamond).");
 
            // Block-level global IDs for this subdomain (numbered within the subdomain, [0, nGlobalBlock)).
            std::vector<GlobalIndex> gids;
            std::vector<bool> owned, ghost;
            const GlobalIndex nGlobalBlock = buildSubdomainGIDs_<SubTraits, canComm>(
                gg.gridView(), subDofMapper_<SubTraits>(gg), 0, gids, owned, ghost);
            const std::size_t nBlockLocal = gids.size();
 
            // In overlapping decompositions the overlap supplies the neighbour contributions to the
            // owner's border rows. In non-overlapping ones every process assembles only its own
            // elements, so a border row is split among the processes sharing the dof.
            const bool rowsComplete = !SubTraits::isNonOverlapping(gg.gridView());
 
            // Append this subdomain's scalar DOFs to the global bookkeeping (in subvector order,
            // matching MatrixConverter/VectorConverter's flattened scalar layout).
            for (std::size_t b = 0; b < nBlockLocal; ++b)
                for (std::size_t k = 0; k < scalarBS; ++k)
                {
                    localToGlobal_.push_back(scalarGlobalOffset + scalarBS * gids[b] + static_cast<GlobalIndex>(k));
                    isOwned_.push_back(owned[b]);
                    isGhost_.push_back(ghost[b]);
                    rowComplete_.push_back(rowsComplete);
                }
 
            subs_[i].scalarBS = scalarBS;
            subs_[i].nBlockLocal = nBlockLocal;
            subs_[i].scalarOffset = scalarLocalOffset;
            subs_[i].rowsComplete = rowsComplete;
            if (isParallel_)
                if constexpr (canComm)
                    setupSubdomainClosures_<i>(gg);
 
            scalarLocalOffset += scalarBS * nBlockLocal;
            scalarGlobalOffset += static_cast<GlobalIndex>(scalarBS) * nGlobalBlock;
        });
 
        numBlocksGlobal_ = static_cast<std::size_t>(scalarGlobalOffset);
    }

    template<std::size_t i, class SubGG>
    void setupSubdomainClosures_(const SubGG& gg)
    {
        using SubTraits = LinearSolverTraits<SubGG>;
        using SubBlockVec = std::tuple_element_t<i, XVector>;
        using GV = typename SubGG::GridView;
        using Mapper = typename SubTraits::DofMapper;
        constexpr std::size_t scalarBS = SubBlockVec::block_type::size();
 
        const GV gv = gg.gridView();
        // Store our own copy of the DOF mapper so the captured closures own a stable instance
        // regardless of whether the grid geometry hands it out by reference or by value.
        auto mapper = std::make_shared<Mapper>(subDofMapper_<SubTraits>(gg));
        const std::size_t offset = subs_[i].scalarOffset;
        const std::size_t nBlock = subs_[i].nBlockLocal;
 
        // Pack/unpack helpers between the scalar slice and a typed subdomain block vector.
        auto unpack = [offset, nBlock](const ScalarBlockVector& v, SubBlockVec& block) {
            block.resize(nBlock);
            for (std::size_t b = 0; b < nBlock; ++b)
                for (std::size_t k = 0; k < scalarBS; ++k)
                    block[b][k] = v[offset + scalarBS * b + k][0];
        };
        auto pack = [offset, nBlock](const SubBlockVec& block, ScalarBlockVector& v) {
            for (std::size_t b = 0; b < nBlock; ++b)
                for (std::size_t k = 0; k < scalarBS; ++k)
                    v[offset + scalarBS * b + k] = block[b][k];
        };
 
        // Ghost broadcast: sum each DOF's copies to the owner's value (callers zero non-owned first).
        subs_[i].ghostBroadcast = [gv, mapper, unpack, pack](ScalarBlockVector& v) {
            SubBlockVec block;
            unpack(v, block);
            if constexpr (requires { SubTraits::dofCodims; }) {
                MultiCodimVectorCommDataHandleSum<Mapper, SubBlockVec, GV::dimension>
                    handle(*mapper, block, SubTraits::dofCodims);
                gv.communicate(handle, Dune::All_All_Interface, Dune::ForwardCommunication);
            } else {
                VectorCommDataHandleSum<Mapper, SubBlockVec, SubTraits::dofCodim>
                    handle(*mapper, block);
                gv.communicate(handle, Dune::All_All_Interface, Dune::ForwardCommunication);
            }
            pack(block, v);
        };
 
        // RHS preprocessing: only needed for incomplete owner rows (sum partial border
        // contributions onto the owner before the centralized gather).
        if (!subs_[i].rowsComplete) {
            subs_[i].rhsPreprocess = [gv, mapper, unpack, pack](ScalarBlockVector& v) {
                SubBlockVec block;
                unpack(v, block);
                if constexpr (requires { SubTraits::dofCodims; })
                {
                    MultiCodimParallelVectorHelper<GV, Mapper> vecHelper(gv, *mapper);
                    vecHelper.makeNonOverlappingConsistent(block, SubTraits::dofCodims);
                }
                else
                {
                    ParallelVectorHelper<GV, Mapper, SubTraits::dofCodim> vecHelper(gv, *mapper);
                    vecHelper.makeNonOverlappingConsistent(block);
                }
                pack(block, v);
            };
        }
    }

    template<class MTMatrix>
    void initSeqDomainOffsets_(const MTMatrix& A)
    {
        constexpr std::size_t numSubDomains = MTMatrix::N();
        if (subs_.size() < numSubDomains)
            subs_.resize(numSubDomains);
        std::size_t offset = 0;
        using namespace Dune::Hybrid;
        forEach(std::make_index_sequence<numSubDomains>{}, [&](auto I)
        {
            constexpr std::size_t i = decltype(I)::value;
            using Block = typename std::decay_t<decltype(A[Dune::index_constant<i>{}][Dune::index_constant<i>{}])>::block_type;
            subs_[i].scalarBS = Block::rows;
            subs_[i].nBlockLocal = A[Dune::index_constant<i>{}][Dune::index_constant<i>{}].N();
            subs_[i].scalarOffset = offset;
            offset += subs_[i].scalarBS * subs_[i].nBlockLocal;
        });
    }
 
    template<class MTMatrix>
    bool subdomainSizesMatch_(const MTMatrix& A) const
    {
        bool match = true;
        using namespace Dune::Hybrid;
        forEach(std::make_index_sequence<MTMatrix::N()>{}, [&](auto I)
        {
            constexpr std::size_t i = decltype(I)::value;
            match = match && A[I][I].N() == subs_[i].nBlockLocal;
        });
        return match;
    }
 
    // Initialize DOF maps for sequential use from a known scalar DOF count.
    void initSeqMaps_(std::size_t nDofs)
    {
        numBlocksGlobal_ = nDofs;
        localToGlobal_.resize(nDofs);
        isOwned_.assign(nDofs, true);
        isGhost_.assign(nDofs, false);
        rowComplete_.assign(nDofs, true);
        for (std::size_t i = 0; i < nDofs; ++i)
            localToGlobal_[i] = static_cast<GlobalIndex>(i);
    }

    void buildIsSrcRow_()
    {
        const std::size_t n = isOwned_.size();
        isSrcRow_.assign(n, false);
        if constexpr (isMultiType)
        {
            for (std::size_t r = 0; r < n; ++r)
                isSrcRow_[r] = rowComplete_[r] ? isOwned_[r] : !isGhost_[r];
            return;
        }
        else if constexpr (LSTraits::canCommunicate)
        {
            if (isParallel_)
            {
                const bool nonOverlapping = LSTraits::isNonOverlapping(*gridView_);
                for (std::size_t r = 0; r < n; ++r)
                    isSrcRow_[r] = nonOverlapping ? !isGhost_[r] : isOwned_[r];
                return;
            }
        }
        isSrcRow_ = isOwned_;
    }

    auto tripletWriter_(bool firstCall, std::size_t& k, bool& samePattern)
    {
        if (firstCall)
        {
            irnLoc_.clear();
            jcnLoc_.clear();
            aLoc_.clear();
        }
 
        return [this, firstCall, &k, &samePattern](GlobalIndex row, GlobalIndex col, double value)
        {
            if (firstCall)
            {
                irnLoc_.push_back(toMumpsInt_(row + 1));
                jcnLoc_.push_back(toMumpsInt_(col + 1));
                aLoc_.push_back(value);
            }
            else if (k >= aLoc_.size() || irnLoc_[k] != row + 1 || jcnLoc_[k] != col + 1)
                samePattern = false;
            else if (samePattern)
                aLoc_[k] = value;
            ++k;
        };
    }

    template<class AMatrix>
    bool assembleTriplets_(const AMatrix& A, bool firstCall)
    requires (isMultiTypeBlockMatrix<AMatrix>::value)
    {
        using namespace Dune::Hybrid;
        constexpr std::size_t numSubDomains = AMatrix::N();
        std::size_t k = 0;
        bool samePattern = true;
        auto addTriplet = tripletWriter_(firstCall, k, samePattern);
        forEach(std::make_index_sequence<numSubDomains>{}, [&](auto I)
        {
            forEach(std::make_index_sequence<numSubDomains>{}, [&](auto J)
            {
                const auto& block = A[I][J];
                using BT = typename std::decay_t<decltype(block)>::block_type;
                const std::size_t offI = subs_[decltype(I)::value].scalarOffset;
                const std::size_t offJ = subs_[decltype(J)::value].scalarOffset;
                for (auto rowIt = block.begin(); rowIt != block.end(); ++rowIt)
                {
                    for (int ki = 0; ki < BT::rows; ++ki)
                    {
                        const std::size_t r = offI + BT::rows*rowIt.index() + ki;
                        if (!isSrcRow_[r]) continue;
                        const GlobalIndex gRow = localToGlobal_[r];
                        for (auto colIt = rowIt->begin(); colIt != rowIt->end(); ++colIt)
                            for (int kj = 0; kj < BT::cols; ++kj)
                                addTriplet(gRow, localToGlobal_[offJ + BT::cols*colIt.index() + kj], (*colIt)[ki][kj]);
                    }
                }
            });
        });
        return samePattern && k == aLoc_.size();
    }
 
    template<class AMatrix>
    bool assembleTriplets_(const AMatrix& A, bool firstCall)
    requires (!isMultiTypeBlockMatrix<AMatrix>::value)
    {
        using BlockType = typename AMatrix::block_type;
        constexpr int BS_I = BlockType::rows;
        constexpr int BS_J = BlockType::cols;
 
        std::size_t k = 0;
        bool samePattern = true;
        auto addTriplet = tripletWriter_(firstCall, k, samePattern);
        for (auto rowIt = A.begin(); rowIt != A.end(); ++rowIt)
        {
            const std::size_t r = rowIt.index();
            if (!isSrcRow_[r]) continue;
            const GlobalIndex gBlockRow = localToGlobal_[r];
 
            for (int ki = 0; ki < BS_I; ++ki)
            {
                const GlobalIndex gRow = gBlockRow * BS_I + ki;
                for (auto colIt = rowIt->begin(); colIt != rowIt->end(); ++colIt)
                {
                    const GlobalIndex gBlockCol = localToGlobal_[colIt.index()];
                    for (int kj = 0; kj < BS_J; ++kj)
                        addTriplet(gRow, gBlockCol * BS_J + kj, (*colIt)[ki][kj]);
                }
            }
        }
        return samePattern && k == aLoc_.size();
    }

    void setupRhsGather_()
    {
        // MPI counts and displacements are int
        const std::size_t nScalar = numBlocksGlobal_ * blockSize;
        if (nScalar > static_cast<std::size_t>(std::numeric_limits<int>::max()))
            DUNE_THROW(Dune::NotImplemented, "DirectSolverMumps: gathering the right-hand side of "
                       << nScalar << " unknowns");
 
        // Owned scalar global indices in the (ascending-local) packing order.
        std::vector<int> myIdx;
        for (std::size_t i = 0; i < localToGlobal_.size(); ++i)
            if (isOwned_[i])
                for (std::size_t k = 0; k < blockSize; ++k)
                    myIdx.push_back(static_cast<int>(localToGlobal_[i] * blockSize + k));
        nOwnedScalar_ = static_cast<int>(myIdx.size());
 
        int commSize = 1, myRank = 0;
        MPI_Comm_size(mpiComm_, &commSize);
        MPI_Comm_rank(mpiComm_, &myRank);
        isHost_ = (myRank == 0);
 
        recvCounts_.assign(isHost_ ? commSize : 0, 0);
        MPI_Gather(&nOwnedScalar_, 1, MPI_INT,
                   isHost_ ? recvCounts_.data() : nullptr, 1, MPI_INT, 0, mpiComm_);
 
        if (isHost_)
        {
            displs_.assign(commSize, 0);
            std::partial_sum(recvCounts_.begin(), recvCounts_.end() - 1, displs_.begin() + 1);
            allGlobalIdx_.assign(nScalar, 0);
            rhsHost_.assign(nScalar, 0.0); // natural (MUMPS) order, length n per right-hand side
            tmpHost_.assign(nScalar, 0.0); // gather order
        }
        // Gather the owned global indices (constant across solves).
        MPI_Gatherv(myIdx.data(), nOwnedScalar_, MPI_INT,
                    isHost_ ? allGlobalIdx_.data() : nullptr,
                    isHost_ ? recvCounts_.data() : nullptr,
                    isHost_ ? displs_.data() : nullptr,
                    MPI_INT, 0, mpiComm_);
        sendBuf_.resize(static_cast<std::size_t>(nOwnedScalar_));
    }

    template<class BVec, class XVec>
    void solveCentralized_(const std::vector<const BVec*>& b, const std::vector<XVec*>& x)
    {
        const auto start = std::chrono::steady_clock::now();
        const std::size_t numColumns = b.size();
        if (numColumns == 0)
            return;
 
        const std::size_t n = tmpHost_.size();
        if (isHost_)
            rhsHost_.resize(n*numColumns);
 
        for (std::size_t j = 0; j < numColumns; ++j)
        {
            // Pack owned scalar RHS values in the same order as setupRhsGather_'s indices.
            const auto& bj = *b[j];
            std::size_t s = 0;
            for (std::size_t i = 0; i < bj.size(); ++i)
                if (isOwned_[i])
                    for (std::size_t k = 0; k < blockSize; ++k)
                        sendBuf_[s++] = bj[i][k];
 
            MPI_Gatherv(sendBuf_.data(), nOwnedScalar_, MPI_DOUBLE,
                        isHost_ ? tmpHost_.data() : nullptr,
                        isHost_ ? recvCounts_.data() : nullptr,
                        isHost_ ? displs_.data() : nullptr,
                        MPI_DOUBLE, 0, mpiComm_);
 
            if (isHost_)
                for (std::size_t p = 0; p < n; ++p)
                    rhsHost_[j*n + static_cast<std::size_t>(allGlobalIdx_[p])] = tmpHost_[p];
        }
 
        if (isHost_)
        {
            id_.nrhs = static_cast<MUMPS_INT>(numColumns);
            id_.lrhs = id_.n;
            id_.rhs = rhsHost_.data();
        }
        id_.job = 3; // solve
        dmumps_c(&id_);
        checkError_("solve");
 
        for (std::size_t j = 0; j < numColumns; ++j)
        {
            // The solutions are centralized on the host (natural order); permute back to gather order.
            if (isHost_)
                for (std::size_t p = 0; p < n; ++p)
                    tmpHost_[p] = rhsHost_[j*n + static_cast<std::size_t>(allGlobalIdx_[p])];
 
            MPI_Scatterv(isHost_ ? tmpHost_.data() : nullptr,
                         isHost_ ? recvCounts_.data() : nullptr,
                         isHost_ ? displs_.data() : nullptr,
                         MPI_DOUBLE, sendBuf_.data(), nOwnedScalar_, MPI_DOUBLE, 0, mpiComm_);
 
            auto& xj = *x[j];
            std::size_t s = 0;
            for (std::size_t i = 0; i < xj.size(); ++i)
                if (isOwned_[i])
                    for (std::size_t k = 0; k < blockSize; ++k)
                        xj[i][k] = sendBuf_[s++];
        }
 
        if (this->verbosity() > 0 && isHost_)
            std::cout << "[MUMPS] solve+gather "
                      << std::chrono::duration<double, std::milli>(std::chrono::steady_clock::now() - start).count()
                      << " ms" << (numColumns > 1 ? " for " + std::to_string(numColumns) + " right-hand sides" : std::string{})
                      << std::endl;
    }
 
    // Symbolic analysis (repeated only if the matrix structure changed) and numeric factorization
    void factorize_(const Matrix& A)
    {
        using Clock = std::chrono::steady_clock;
        const bool time = this->verbosity() > 0;
        auto ms = [](Clock::time_point a, Clock::time_point z) {
            return std::chrono::duration<double, std::milli>(z - a).count();
        };
 
        double tAssemble = 0.0, tAnalysis = 0.0;
        auto t0 = Clock::now();
        // the analysis is collective, so all processes redo it if one process' structure changed
        const bool analyze = !analyzed_ || !onAllProcesses_(assembleTriplets_(A, /*firstCall=*/false));
        if (analyze)
        {
            buildIsSrcRow_();
            setupRhsGather_();
 
            assembleTriplets_(A, /*firstCall=*/true);
            id_.n = toMumpsInt_(static_cast<GlobalIndex>(numBlocksGlobal_ * blockSize));
            id_.nz_loc = toMumpsInt_(static_cast<GlobalIndex>(aLoc_.size()));
            id_.irn_loc = irnLoc_.data();
            id_.jcn_loc = jcnLoc_.data();
            id_.a_loc = aLoc_.data();
            auto t1 = Clock::now(); tAssemble = ms(t0, t1);
 
            id_.job = 1; // symbolic analysis
            dmumps_c(&id_);
            checkError_("symbolic analysis");
            analyzed_ = true;
            tAnalysis = ms(t1, Clock::now());
        }
        else
            tAssemble = ms(t0, Clock::now());
 
        auto t2 = Clock::now();
        id_.job = 2; // numeric factorization
        dmumps_c(&id_);
        checkError_("numeric factorization");
        factorized_ = true;
        const double tNumeric = ms(t2, Clock::now());
 
        if (time && isHost_)
            std::cout << "[MUMPS] assemble " << tAssemble << " ms"
                      << (tAnalysis > 0.0 ? (", analysis " + std::to_string(tAnalysis) + " ms") : std::string{})
                      << ", numeric factorization " << tNumeric << " ms" << std::endl;
    }
 
    // Sequential dof maps follow the size of the system of A, parallel ones are checked against it.
    void prepareMaps_(const Matrix& A)
    {
        if constexpr (isMultiType)
        {
            if (isParallel_)
                updateMapsToSystemSize_(numScalarRows_(A), subdomainSizesMatch_(A));
            else
            {
                updateMapsToSystemSize_(numScalarRows_(A));
                initSeqDomainOffsets_(A);
            }
        }
        else
            updateMapsToSystemSize_(A.N());
    }
 
    template<class MTMatrix>
    static std::size_t numScalarRows_(const MTMatrix& A)
    {
        std::size_t numRows = 0;
        using namespace Dune::Hybrid;
        forEach(std::make_index_sequence<MTMatrix::N()>{}, [&](auto I)
        {
            const auto& block = A[I][I];
            numRows += block.N()*std::decay_t<decltype(block)>::block_type::rows;
        });
        return numRows;
    }
 
    SolverResult solveColumns_(const std::vector<XVector*>& x, const std::vector<const BVector*>& b)
    {
        if (!factorized_)
            DUNE_THROW(Dune::InvalidStateException, "Called solve(x, b) but no matrix has been set");
 
        checkColumns_(x, b);
        result_.clear();
        if constexpr (isMultiType)
            solveMultiType_(x, b);
        else
            solveSingleType_(x, b);
        result_.converged = true;
        return result_;
    }
 
    // The vectors have to match the factorized system and the number of right-hand sides has to be
    // the same on all processes, which take part in every gather; otherwise all processes throw.
    void checkColumns_(const std::vector<XVector*>& x, const std::vector<const BVector*>& b) const
    {
        bool match = x.size() == b.size();
        for (std::size_t j = 0; match && j < x.size(); ++j)
            match = sizeInDofMaps_(*x[j]) == localToGlobal_.size() && sizeInDofMaps_(*b[j]) == localToGlobal_.size();
 
        if (isParallel_)
        {
            long long counts[2] = { static_cast<long long>(b.size()), -static_cast<long long>(b.size()) };
            MPI_Allreduce(MPI_IN_PLACE, counts, 2, MPI_LONG_LONG, MPI_MIN, mpiComm_);
            match = match && counts[0] == -counts[1];
        }
 
        if (!onAllProcesses_(match))
            DUNE_THROW(Dune::InvalidStateException,
                       "DirectSolverMumps: the solutions and right-hand sides do not match each other or the "
                       "factorized system on at least one process (here " << x.size() << " solutions and "
                       << b.size() << " right-hand sides, the system has " << localToGlobal_.size() << " rows)");
    }
 
    // The size of a vector in the entries of the dof maps (blocks for single-type, scalars for MultiType)
    template<class V>
    static std::size_t sizeInDofMaps_(const V& v)
    {
        if constexpr (isMultiTypeBlockVector<V>::value)
        {
            std::size_t size = 0;
            Dune::Hybrid::forEach(std::make_index_sequence<V::N()>{}, [&](auto i)
            {
                using Block = typename std::decay_t<decltype(v[i])>::block_type;
                size += v[i].size()*Block::size();
            });
            return size;
        }
        else
            return v.size();
    }
 
    bool onAllProcesses_(bool value) const
    {
        if (!isParallel_)
            return value;
        int local = value ? 1 : 0;
        int global = 0;
        MPI_Allreduce(&local, &global, 1, MPI_INT, MPI_MIN, mpiComm_);
        return global == 1;
    }
 
    // The sequential dof maps follow the size of the system; the parallel ones are rebuilt by
    // updateAfterGridAdaption. In parallel, a mismatch on one process makes all processes throw.
    void updateMapsToSystemSize_(std::size_t numRows, bool subdomainSizesMatch = true)
    {
        const bool sizesMatch = numRows == localToGlobal_.size() && subdomainSizesMatch;
        if (isParallel_)
        {
            if (!onAllProcesses_(sizesMatch))
                DUNE_THROW(Dune::InvalidStateException,
                           "DirectSolverMumps: the dof distribution of the system differs from the one of "
                           "the solver on at least one process (here " << numRows << " rows, the solver "
                           "expects " << localToGlobal_.size() << "). Call updateAfterGridAdaption after "
                           "the grid changed.");
        }
        else if (!sizesMatch)
        {
            initSeqMaps_(numRows);
            analyzed_ = false;
        }
    }
 
    void solveSingleType_(const std::vector<XVector*>& x, const std::vector<const BVector*>& b)
    {
        std::vector<BVector> bConsistent;
        bConsistent.reserve(b.size());
        for (const auto* bj : b)
        {
            bConsistent.push_back(*bj);
            makeRhsConsistent_(bConsistent.back());
        }
 
        std::vector<const BVector*> bColumns;
        for (const auto& bj : bConsistent)
            bColumns.push_back(&bj);
        solveCentralized_(bColumns, x);
 
        for (auto* xj : x)
            distributeSolution_(*xj);
    }
 
    // Non-overlapping parallel: accumulate partial border RHS onto owners before the gather.
    void makeRhsConsistent_(BVector& b) const
    {
        if constexpr (LSTraits::canCommunicate)
        {
            if (isParallel_ && LSTraits::isNonOverlapping(*gridView_))
            {
                if constexpr (requires { LSTraits::dofCodims; })
                {
                    MultiCodimParallelVectorHelper<GridView, DofMapper> vecHelper(*gridView_, *mapper_);
                    vecHelper.makeNonOverlappingConsistent(b, LSTraits::dofCodims);
                }
                else
                {
                    ParallelVectorHelper<GridView, DofMapper, LSTraits::dofCodim> vecHelper(*gridView_, *mapper_);
                    vecHelper.makeNonOverlappingConsistent(b);
                }
            }
        }
    }
 
    // Distribute the solution: zero non-owned, then sum-communicate to the owner's value.
    void distributeSolution_(XVector& x) const
    {
        if constexpr (LSTraits::canCommunicate)
        {
            if (isParallel_)
            {
                for (std::size_t i = 0; i < x.size(); ++i)
                    if (!isOwned_[i])
                        x[i] = 0;
                if constexpr (requires { LSTraits::dofCodims; }) {
                    MultiCodimVectorCommDataHandleSum<DofMapper, XVector, GridView::dimension> handle(*mapper_, x, LSTraits::dofCodims);
                    gridView_->communicate(handle, Dune::All_All_Interface, Dune::ForwardCommunication);
                } else {
                    VectorCommDataHandleSum<DofMapper, XVector, LSTraits::dofCodim> handle(*mapper_, x);
                    gridView_->communicate(handle, Dune::All_All_Interface, Dune::ForwardCommunication);
                }
            }
        }
    }
 
    void solveMultiType_(const std::vector<XVector*>& x, const std::vector<const BVector*>& b)
    {
        // The flattened vectors have the scalar layout of the system factorized by MUMPS.
        std::vector<ScalarBlockVector> bScalar;
        bScalar.reserve(b.size());
        for (const auto* bj : b)
        {
            bScalar.push_back(VectorConverter<BVector>::multiTypeToBlockVector(*bj));
 
            // RHS preprocessing per subdomain: sum partial border-DOF contributions onto the owner.
            if (isParallel_)
                for (const auto& sub : subs_)
                    if (sub.rhsPreprocess)
                        sub.rhsPreprocess(bScalar.back());
        }
 
        auto xScalar = bScalar; // same layout; owned entries overwritten by the solve
        std::vector<const ScalarBlockVector*> bColumns;
        std::vector<ScalarBlockVector*> xColumns;
        for (std::size_t j = 0; j < bScalar.size(); ++j)
        {
            bColumns.push_back(&bScalar[j]);
            xColumns.push_back(&xScalar[j]);
        }
        solveCentralized_(bColumns, xColumns);
 
        for (std::size_t j = 0; j < x.size(); ++j)
        {
            // Ghost broadcast per subdomain: zero non-owned DOFs, sum-communicate to the owner's value.
            if (isParallel_)
            {
                for (std::size_t r = 0; r < xScalar[j].size(); ++r)
                    if (!isOwned_[r]) xScalar[j][r] = 0.0;
                for (const auto& sub : subs_)
                    if (sub.ghostBroadcast)
                        sub.ghostBroadcast(xScalar[j]);
            }
 
            VectorConverter<XVector>::retrieveValues(*x[j], xScalar[j]);
        }
    }
 
    Scalar ownedGlobalNorm_(const XVector& x) const
    {
        double localSumSq = 0.0;
        for (std::size_t i = 0; i < x.size(); ++i)
            if (isOwned_[i])
                for (std::size_t k = 0; k < blockSize; ++k)
                    localSumSq += x[i][k] * x[i][k];
        double globalSumSq = 0.0;
        MPI_Allreduce(&localSumSq, &globalSumSq, 1, MPI_DOUBLE, MPI_SUM, mpiComm_);
        using std::sqrt;
        return sqrt(globalSumSq);
    }
 
    static MUMPS_INT toMumpsInt_(GlobalIndex v)
    {
        if (v > static_cast<GlobalIndex>(std::numeric_limits<MUMPS_INT>::max()))
            DUNE_THROW(Dune::RangeError, "Index " << v << " exceeds MUMPS_INT range; "
                       "rebuild MUMPS with 64-bit integers (-DINTSIZE64)");
        return static_cast<MUMPS_INT>(v);
    }
 
    void checkError_(const char* phase) const
    {
        const int error = infog_(1);
        if (error < 0)
        {
            // the errors that the MUMPS user guide resolves by increasing ICNTL(14)
            const bool workspaceTooSmall = error == -8 || error == -9 || error == -14
                                           || error == -15 || error == -17 || error == -20;
            DUNE_THROW(Dune::Exception, "MUMPS " << phase << " failed: INFOG(1)=" << error
                       << ", INFOG(2)=" << infog_(2)
                       << (workspaceTooSmall ? ". Increase LinearSolver.MumpsWorkspaceIncrease (ICNTL(14))" : ""));
        }
    }
 
    std::optional<GridView> gridView_;
    const DofMapper* mapper_ = nullptr;
    bool isParallel_;
 
    MPI_Comm mpiComm_;
    DMUMPS_STRUC_C id_;
    bool initialized_ = false;
    bool analyzed_ = false;
    bool factorized_ = false;
    SolverResult result_;
 
    // Global DOF bookkeeping (block level for single-type, scalar level for MultiType).
    std::size_t numBlocksGlobal_ = 0;
    std::vector<GlobalIndex> localToGlobal_;
    std::vector<bool> isOwned_;
    std::vector<bool> isGhost_;
    std::vector<bool> isSrcRow_; // block rows this rank contributes to the distributed assembly
    std::vector<bool> rowComplete_; // per (scalar) row: owner row already complete? (MultiType)
 
    // Distributed assembled input (ICNTL(18)=3). irn_loc/jcn_loc hold the structure of the last
    // symbolic analysis; a_loc is refilled every solve.
    std::vector<MUMPS_INT> irnLoc_;
    std::vector<MUMPS_INT> jcnLoc_;
    std::vector<double> aLoc_;
 
    // Centralized dense RHS/solution gather (host = rank 0). The assembled RHS is gathered to the
    // host, solved, and the solution scattered back.
    bool isHost_ = true;
    int nOwnedScalar_ = 0;
    std::vector<int> recvCounts_; // host only
    std::vector<int> displs_; // host only
    std::vector<int> allGlobalIdx_; // host only, gather-order -> natural global index
    std::vector<double> rhsHost_; // host only, natural (MUMPS) order, length n per right-hand side
    std::vector<double> tmpHost_; // host only, gather order
    std::vector<double> sendBuf_; // per-rank owned scalar values / solution
 
    // MultiType (multidomain) bookkeeping: one entry per subdomain, in subvector order.
    std::vector<SubdomainInfo> subs_;
};
 
} // end namespace Dumux
 
#endif // DUMUX_HAVE_MUMPS
 
#endif