stencil_solver.tpp

/*---------------------------------------------------------------------------*\
 *
 *  bitpit
 *
 *  Copyright (C) 2015-2021 OPTIMAD engineering Srl
 *
 *  -------------------------------------------------------------------------
 *  License
 *  This file is part of bitpit.
 *
 *  bitpit is free software: you can redistribute it and/or modify it
 *  under the terms of the GNU Lesser General Public License v3 (LGPL)
 *  as published by the Free Software Foundation.
 *
 *  bitpit 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 Lesser General Public
 *  License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public License
 *  along with bitpit. If not, see <http://www.gnu.org/licenses/>.
 *
\*---------------------------------------------------------------------------*/
 
#include "stencil_solver.hpp"
 
#include <cmath>
 
namespace bitpit {



template<typename stencil_t, typename stencil_container_t>
DiscretizationStencilProxyBaseStorage<stencil_t, stencil_container_t>::DiscretizationStencilProxyBaseStorage(const stencil_container_t *stencils, int stride)
    : DiscretizationStencilStorageInterface<stencil_t>(),
      m_stencils(stencils), m_stride(stride)
{
}

template<typename stencil_t, typename stencil_container_t>
DiscretizationStencilProxyStorage<stencil_t, stencil_container_t>::DiscretizationStencilProxyStorage(const stencil_container_t *stencils)
    : DiscretizationStencilProxyBaseStorage<stencil_t, stencil_container_t>(stencils, 1)
{
}

 * \result The size of the container.
 */
template<typename stencil_t, typename stencil_container_t>
std::size_t DiscretizationStencilProxyStorage<stencil_t, stencil_container_t>::size() const
{
    return this->m_stencils->size();
}

/*!
 * Get the stencil associated with the specified row.
 *
 * \param rowIndex is the index of the row in the storage
 * \result The stencil associated with the specified row.
 */
template<typename stencil_t, typename stencil_container_t>
const stencil_t & DiscretizationStencilProxyStorage<stencil_t, stencil_container_t>::at(long rowIndex) const
{
    return (*(this->m_stencils))[rowIndex];
}

template<typename stencil_t, typename stencil_container_t>
const stencil_t & DiscretizationStencilProxyStorage<stencil_t, stencil_container_t>::rawAt(std::size_t rowRawIndex) const
{
    return (*(this->m_stencils))[rowRawIndex];
}

template<typename stencil_t, typename stencil_container_t>
const stencil_t & DiscretizationStencilProxyStorage<stencil_t, stencil_container_t>::at(long blockIndex, int componentIdx) const
{
    return (*(this->m_stencils))[blockIndex * this->m_stride + componentIdx];
}

 *
 * \param blockRawIndex is the raw index of the block
 * \param componentIdx is the index of the component inside the block
 * \result The stencil associated with the specified row.
 */
template<typename stencil_t, typename stencil_container_t>
const stencil_t & DiscretizationStencilProxyStorage<stencil_t, stencil_container_t>::rawAt(std::size_t blockRawIndex, int componentIdx) const
{
    return (*(this->m_stencils))[blockRawIndex * this->m_stride + componentIdx];
}

/*!
 * Constructor.
 *
 * \param stencils are the stencils
 */
template<typename stencil_t>
DiscretizationStencilProxyStorage<stencil_t, PiercedStorage<stencil_t>>::DiscretizationStencilProxyStorage(const PiercedStorage<stencil_t> *stencils)
    : DiscretizationStencilProxyBaseStorage<stencil_t, PiercedStorage<stencil_t>>(stencils, stencils->getFieldCount())
{
}

/*!
 * Get the size of the container.
 *
 * \result The size of the container.
 */
template<typename stencil_t>
std::size_t DiscretizationStencilProxyStorage<stencil_t, PiercedStorage<stencil_t>>::size() const
{
    return this->m_stencils->getKernel()->size() * this->m_stride;
}

/*!
 * Get the stencil associated with the specified row.
 *
 * \param rowIndex is the index of the row in the storage
 * \result The stencil associated with the specified row.
 */
template<typename stencil_t>
const stencil_t & DiscretizationStencilProxyStorage<stencil_t, PiercedStorage<stencil_t>>::at(long rowIndex) const
{
    return this->m_stencils->at(rowIndex);
}

template<typename stencil_t>
const stencil_t & DiscretizationStencilProxyStorage<stencil_t, PiercedStorage<stencil_t>>::rawAt(std::size_t rowRawIndex) const
{
    return this->m_stencils->rawAt(rowRawIndex);
}

 *
 * \param blockIndex is the index of the block
 * \param componentIdx is the index of the component inside the block
 * \result The stencil associated with the specified row.
 */
template<typename stencil_t>
const stencil_t & DiscretizationStencilProxyStorage<stencil_t, PiercedStorage<stencil_t>>::at(long blockIndex, int componentIdx) const
{
    return this->m_stencils->at(blockIndex, componentIdx);
}

 * \param blockRawIndex is the raw index of the block
 * \param componentIdx is the index of the component inside the block
 * \result The stencil associated with the specified row.
 */
template<typename stencil_t>
const stencil_t & DiscretizationStencilProxyStorage<stencil_t, PiercedStorage<stencil_t>>::rawAt(std::size_t blockRawIndex, int componentIdx) const
{
    return this->m_stencils->rawAt(blockRawIndex, componentIdx);
}

 * \class DiscretizationStencilSolverAssembler
 * \ingroup discretization
 *
 * \brief The DiscretizationStencilSolverAssembler class defines an assembler
 * for building the stencil solver.
 */
 
#if BITPIT_ENABLE_MPI==1
 *
 * \param stencils are the stencils
 * \param assemblerKernelArgs are the arguments that will be passed to the constructor of the
 * assembler of the solver kernel
 */
template<typename stencil_t, typename solver_kernel_t>
template<typename stencil_container_t, typename... AssemblerKernelArgs>
DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::DiscretizationStencilSolverAssembler(const stencil_container_t *stencils,
                                                                                                       AssemblerKernelArgs&&... assemblerKernelArgs)
    : DiscretizationStencilSolverAssembler(MPI_COMM_SELF, false, stencils, std::forward<AssemblerKernelArgs>(assemblerKernelArgs)...)
{
}

template<typename stencil_t, typename solver_kernel_t>
template<typename stencil_container_t, typename... AssemblerKernelArgs>
DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::DiscretizationStencilSolverAssembler(MPI_Comm communicator, bool partitioned,
                                                                                                       const stencil_container_t *stencils,
                                                                                                       AssemblerKernelArgs&&... assemblerKernelArgs)
    : DiscretizationStencilSolverAssembler(communicator, partitioned,
                                           std::unique_ptr<DiscretizationStencilStorageInterface<stencil_t>>(new DiscretizationStencilProxyStorage<stencil_t, stencil_container_t>(stencils)),
                                           std::forward<AssemblerKernelArgs>(assemblerKernelArgs)...)
{
}
#else
/*!
 * Constructor.
 *
 * \param stencils are the stencils
 * \param assemblerKernelArgs are the arguments that will be passed to the constructor of the
 * assembler of the solver kernel
 */
template<typename stencil_t, typename solver_kernel_t>
template<typename stencil_container_t, typename... AssemblerKernelArgs>
DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::DiscretizationStencilSolverAssembler(const stencil_container_t *stencils,
                                                                                                       AssemblerKernelArgs&&... assemblerKernelArgs)
    : DiscretizationStencilSolverAssembler(std::unique_ptr<DiscretizationStencilStorageInterface<stencil_t>>(new DiscretizationStencilProxyStorage<stencil_t, stencil_container_t>(stencils)),
                                           std::forward<AssemblerKernelArgs>(assemblerKernelArgs)...)
{
}
#endif
 
#if BITPIT_ENABLE_MPI==1
 * assembler of the solver kernel
 */
template<typename stencil_t, typename solver_kernel_t>
template<typename... AssemblerKernelArgs>
DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::DiscretizationStencilSolverAssembler(std::unique_ptr<DiscretizationStencilStorageInterface<stencil_t>> &&stencils,
                                                                                                       AssemblerKernelArgs&&... assemblerKernelArgs)
    : DiscretizationStencilSolverAssembler(MPI_COMM_SELF, false, std::move(stencils), std::forward<AssemblerKernelArgs>(assemblerKernelArgs)...)
{
}

template<typename stencil_t, typename solver_kernel_t>
template<typename... AssemblerKernelArgs>
DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::DiscretizationStencilSolverAssembler(MPI_Comm communicator, bool partitioned,
                                                                                                       std::unique_ptr<DiscretizationStencilStorageInterface<stencil_t>> &&stencils,
                                                                                                       AssemblerKernelArgs&&... assemblerKernelArgs)
    : DiscretizationStencilSolverAssembler(communicator, partitioned, std::forward<AssemblerKernelArgs>(assemblerKernelArgs)...)
#else
template<typename stencil_t, typename solver_kernel_t>
template<typename... AssemblerKernelArgs>
DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::DiscretizationStencilSolverAssembler(std::unique_ptr<DiscretizationStencilStorageInterface<stencil_t>> &&stencils,
                                                                                                       AssemblerKernelArgs&&... assemblerKernelArgs)
    : DiscretizationStencilSolverAssembler(std::forward<AssemblerKernelArgs>(assemblerKernelArgs)...)
#endif
{
    setStencils(std::move(stencils));
    setMatrixSizes();
    setBlockSize();
    setMaximumRowNZ();
}
 
#if BITPIT_ENABLE_MPI==1
template<typename stencil_t, typename solver_kernel_t>
template<typename... AssemblerKernelArgs>
DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::DiscretizationStencilSolverAssembler(AssemblerKernelArgs&&... assemblerKernelArgs)
    : DiscretizationStencilSolverAssembler(MPI_COMM_SELF, false, std::forward<AssemblerKernelArgs>(assemblerKernelArgs)...)
{
}

template<typename stencil_t, typename solver_kernel_t>
template<typename... AssemblerKernelArgs>
DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::DiscretizationStencilSolverAssembler(MPI_Comm communicator, bool partitioned,
                                                                                                       AssemblerKernelArgs&&... assemblerKernelArgs)
    : solver_kernel_type::Assembler(std::forward<AssemblerKernelArgs>(assemblerKernelArgs)...),
      m_partitioned(partitioned), m_communicator(communicator)
{
}
#else
template<typename stencil_t, typename solver_kernel_t>
template<typename... AssemblerKernelArgs>
DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::DiscretizationStencilSolverAssembler(AssemblerKernelArgs&&... assemblerKernelArgs)
    : solver_kernel_type::Assembler(std::forward<AssemblerKernelArgs>(assemblerKernelArgs)...)
{
}
#endif
 
#if BITPIT_ENABLE_MPI==1
template<typename stencil_t, typename solver_kernel_t>
bool DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::isPartitioned() const
{
    return m_partitioned;
}

template<typename stencil_t, typename solver_kernel_t>
const MPI_Comm & DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::getCommunicator() const
{
    return m_communicator;
}
#endif

template<typename stencil_t, typename solver_kernel_t>
typename DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::assembly_options_type DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::getOptions() const
{
    assembly_options_type options;
    options.full   = true;
    options.sorted = false;
 
    return options;
}

template<typename stencil_t, typename solver_kernel_t>
template<typename W, typename V, typename std::enable_if<std::is_fundamental<W>::value>::type *>
void DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::setBlockSize()
{
    setBlockSize(1);
}

template<typename stencil_t, typename solver_kernel_t>
template<typename W, typename V, std::size_t D, typename std::enable_if<std::is_same<std::array<V, D>, W>::value>::type *>
void DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::setBlockSize()
{
    setBlockSize(sizeof(typename StencilVector::weight_type) / sizeof(typename StencilVector::weight_type::value_type));
}

template<typename stencil_t, typename solver_kernel_t>
template<typename W, typename V, typename std::enable_if<std::is_same<std::vector<V>, W>::value>::type *>
void DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::setBlockSize()
{
    // Get block size
    //
    // Block size is evaluated from the constant of the first stencil.
    //
    // The block size is set equal to the square root of the weight size; if the square
    // root of the weight type is not an integer number, an exception is throw.
    if (getRowCount() == 0) {
        throw std::runtime_error("Unable to evaluate the block size.");
    }
 
    const StencilBlock &stencil = getRowStencil(0);
    std::size_t stencilConstantSize = stencil.getConstant().size();
    int blockSize = static_cast<int>(std::round(std::sqrt(stencilConstantSize)));
    if (static_cast<std::size_t>(blockSize * blockSize) != stencilConstantSize) {
        throw std::runtime_error("Weights size should be a square.");
    }
    setBlockSize(blockSize);
 
#ifdef DEBUG
    // Validate block size
    //
    // All weight sizes should match
    for (long i = 0; i < getRowCount(); ++i) {
        const StencilBlock &stencil = getRowStencil(i);
        const StencilBlock::weight_type *weightData = stencil.weightData();
        std::size_t stencilSize = stencil.size();
 
        for (std::size_t k = 0; k < stencilSize; ++k) {
            if (weightData[k].size() != m_blockSize * m_blockSize)) {
                throw std::runtime_error("All stencils weights should have the same size.");
            }
        }
 
        if (stencil.getConstant().size() != m_blockSize * m_blockSize) {
            throw std::runtime_error("The stencils constant should have the same size of the stencil weights.");
        }
    }
#endif
}

template<typename stencil_t, typename solver_kernel_t>
void DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::setBlockSize(int blockSize)
{
    m_blockSize = blockSize;
}

template<typename stencil_t, typename solver_kernel_t>
void DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::setStencils(std::unique_ptr<DiscretizationStencilStorageInterface<stencil_t>> &&stencils)
{
    m_stencils = std::move(stencils);
}

template<typename stencil_t, typename solver_kernel_t>
void DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::setMatrixSizes()
{
    setMatrixSizes(m_stencils->size(), m_stencils->size());
}

template<typename stencil_t, typename solver_kernel_t>
void DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::setMatrixSizes(long nRows, long nCols)
{
    // Set system sizes
    m_nRows = nRows;
    m_nCols = nCols;
 
#if BITPIT_ENABLE_MPI==1
    // Global system sizes
    m_nGlobalRows = nRows;
    m_nGlobalCols = nCols;
    if (m_partitioned) {
        MPI_Allreduce(MPI_IN_PLACE, &m_nGlobalRows, 1, MPI_LONG, MPI_SUM, m_communicator);
        MPI_Allreduce(MPI_IN_PLACE, &m_nGlobalCols, 1, MPI_LONG, MPI_SUM, m_communicator);
    }
 
    // Global offsets
    m_globalRowOffset = 0;
    m_globalColOffset = 0;
    if (m_partitioned) {
        int nProcessors;
        MPI_Comm_size(m_communicator, &nProcessors);
 
        std::vector<long> nRankRows(nProcessors);
        MPI_Allgather(&m_nRows, 1, MPI_LONG, nRankRows.data(), 1, MPI_LONG, m_communicator);
 
        std::vector<long> nRankCols(nProcessors);
        MPI_Allgather(&m_nCols, 1, MPI_LONG, nRankCols.data(), 1, MPI_LONG, m_communicator);
 
        int rank;
        MPI_Comm_rank(m_communicator, &rank);
        for (int i = 0; i < rank; ++i) {
            m_globalRowOffset += nRankRows[i];
            m_globalColOffset += nRankCols[i];
        }
    }
#endif
}

template<typename stencil_t, typename solver_kernel_t>
void DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::setMaximumRowNZ()
{
    long maxRowNZ = 0;
    for (long n = 0; n < getRowCount(); ++n) {
        maxRowNZ = std::max(getRowNZCount(n), maxRowNZ);
    }
 
    setMaximumRowNZ(maxRowNZ);
}

template<typename stencil_t, typename solver_kernel_t>
void DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::setMaximumRowNZ(long maxRowNZ)
{
    m_maxRowNZ = maxRowNZ;
}

template<typename stencil_t, typename solver_kernel_t>
int DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::getBlockSize() const
{
    return m_blockSize;
}

template<typename stencil_t, typename solver_kernel_t>
long DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::getRowCount() const
{
    return m_nRows;
}

template<typename stencil_t, typename solver_kernel_t>
long DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::getColCount() const
{
    return m_nCols;
}

template<typename stencil_t, typename solver_kernel_t>
long DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::getRowElementCount() const
{
    long nRowElements = getBlockSize() * getRowCount();
 
    return nRowElements;
}

template<typename stencil_t, typename solver_kernel_t>
long DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::getColElementCount() const
{
    long nColElements = getBlockSize() * getColCount();
 
    return nColElements;
}
 
#if BITPIT_ENABLE_MPI==1
template<typename stencil_t, typename solver_kernel_t>
long DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::getRowGlobalCount() const
{
    return m_nGlobalRows;
}

template<typename stencil_t, typename solver_kernel_t>
long DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::getColGlobalCount() const
{
    return m_nGlobalCols;
}

template<typename stencil_t, typename solver_kernel_t>
long DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::getRowGlobalElementCount() const
{
    long nElements = getBlockSize() * getRowGlobalCount();
 
    return nElements;
}

template<typename stencil_t, typename solver_kernel_t>
long DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::getColGlobalElementCount() const
{
    long nElements = getBlockSize() * getColGlobalCount();
 
    return nElements;
}

template<typename stencil_t, typename solver_kernel_t>
long DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::getRowGlobalOffset() const
{
    return m_globalRowOffset;
}

template<typename stencil_t, typename solver_kernel_t>
long DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::getColGlobalOffset() const
{
    return m_globalColOffset;
}

template<typename stencil_t, typename solver_kernel_t>
long DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::getRowGlobalElementOffset() const
{
    long offset = getBlockSize() * getRowGlobalOffset();
 
    return offset;
}

template<typename stencil_t, typename solver_kernel_t>
long DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::getColGlobalElementOffset() const
{
    long offset = getBlockSize() * getColGlobalOffset();
 
    return offset;
}
#endif

template<typename stencil_t, typename solver_kernel_t>
long DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::getRowNZCount(long rowIndex) const
{
    const stencil_t &stencil = getRowStencil(rowIndex);
    std::size_t stencilSize = stencil.size();
 
    return stencilSize;
}

template<typename stencil_t, typename solver_kernel_t>
long DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::getMaxRowNZCount() const
{
    return m_maxRowNZ;
}

template<typename stencil_t, typename solver_kernel_t>
void DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::getRowPattern(long rowIndex, ConstProxyVector<long> *pattern) const
{
    // Get stencil information
    const stencil_t &stencil = getRowStencil(rowIndex);
 
    // Get pattern
    getPattern(stencil, pattern);
}

template<typename stencil_t, typename solver_kernel_t>
void DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::getPattern(const stencil_t &stencil, ConstProxyVector<long> *pattern) const
{
    std::size_t stencilSize = stencil.size();
 
    const long *patternData = stencil.patternData();
    pattern->set(patternData, stencilSize);
}

template<typename stencil_t, typename solver_kernel_t>
void DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::getRowValues(long rowIndex, ConstProxyVector<double> *values) const
{
    // Get stencil information
    const stencil_t &stencil = getRowStencil(rowIndex);
 
    // Get values
    getValues(stencil, values);
}

template<typename stencil_t, typename solver_kernel_t>
template<typename U, typename std::enable_if<std::is_fundamental<U>::value>::type *>
void DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::getValues(const stencil_t &stencil, ConstProxyVector<double> *values) const
{
    values->set(stencil.weightData(), m_blockSize * stencil.size());
}

template<typename stencil_t, typename solver_kernel_t>
template<typename U, typename std::enable_if<!std::is_fundamental<U>::value>::type *>
void DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::getValues(const stencil_t &stencil, ConstProxyVector<double> *values) const
{
    std::size_t stencilSize = stencil.size();
    const stencil_weight_type *stencilWeightData = stencil.weightData();
 
    int nBlockElements = m_blockSize * m_blockSize;
    std::size_t nRowValues = m_blockSize * stencilSize;
    std::size_t nValues = nBlockElements * stencilSize;
    values->set(ConstProxyVector<double>::INTERNAL_STORAGE, nValues);
    ConstProxyVector<double>::storage_pointer expandedValuesStorage = values->storedData();
 
 
    for (std::size_t k = 0; k < stencilSize; ++k) {
        const double *weightData = stencilWeightData[k].data();
        for (int i = 0; i < m_blockSize; ++i) {
            int weightOffset = linearalgebra::linearIndexRowMajor(i, 0, m_blockSize, m_blockSize);
            int valuesOffset = linearalgebra::linearIndexRowMajor(i, m_blockSize * k, m_blockSize, nRowValues);
 
            std::copy_n(weightData + weightOffset, m_blockSize, expandedValuesStorage + valuesOffset);
        }
    }
}

template<typename stencil_t, typename solver_kernel_t>
void DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::getRowData(long rowIndex, ConstProxyVector<long> *pattern, ConstProxyVector<double> *values) const
{
    // Get stencil information
    const stencil_t &stencil = getRowStencil(rowIndex);
 
    // Get pattern
    getPattern(stencil, pattern);
 
    // Get values
    getValues(stencil, values);
}

template<typename stencil_t, typename solver_kernel_t>
void DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::getRowConstant(long rowIndex, bitpit::ConstProxyVector<double> *constant) const
{
    // Get stencil information
    const stencil_t &stencil = getRowStencil(rowIndex);
 
    // Get constant
    getConstant(stencil, constant);
}

template<typename stencil_t, typename solver_kernel_t>
template<typename U, typename std::enable_if<std::is_fundamental<U>::value>::type *>
void DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::getConstant(const stencil_t &stencil, bitpit::ConstProxyVector<double> *constant) const
{
    const stencil_weight_type &stencilConstant = stencil.getConstant();
 
    constant->set(ConstProxyVector<double>::INTERNAL_STORAGE, m_blockSize);
    ConstProxyVector<double>::storage_pointer constantStorage = constant->storedData();
    std::copy_n(&stencilConstant, m_blockSize, constantStorage);
}

template<typename stencil_t, typename solver_kernel_t>
template<typename U, typename std::enable_if<!std::is_fundamental<U>::value>::type *>
void DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::getConstant(const stencil_t &stencil, bitpit::ConstProxyVector<double> *constant) const
{
    const stencil_weight_type &stencilConstant = stencil.getConstant();
 
    constant->set(ConstProxyVector<double>::INTERNAL_STORAGE, m_blockSize);
    ConstProxyVector<double>::storage_pointer constantStorage = constant->storedData();
    for (int i = 0; i < m_blockSize; ++i) {
        int offset_i = linearalgebra::linearIndexRowMajor(i, 0, m_blockSize, m_blockSize);
 
        constantStorage[i] = 0;
        for (int j = 0; j < m_blockSize; ++j) {
            constantStorage[i] += stencil.getWeightManager().at(stencilConstant, offset_i + j);
        }
    }
}

template<typename stencil_t, typename solver_kernel_t>
const stencil_t & DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t>::getRowStencil(long rowIndex) const
{
    return m_stencils->at(rowIndex);
}


template<typename stencil_t, typename solver_kernel_t>
void DiscretizationStencilSolver<stencil_t, solver_kernel_t>::clear()
{
    solver_kernel_t::clear();
 
    std::vector<double>().swap(m_constants);
}
 
#if BITPIT_ENABLE_MPI==1
template<typename stencil_t, typename solver_kernel_t>
template<typename stencil_container_t>
void DiscretizationStencilSolver<stencil_t, solver_kernel_t>::assembly(const stencil_container_t &stencils)
{
    assembly(MPI_COMM_SELF, false, stencils);
}

template<typename stencil_t, typename solver_kernel_t>
template<typename stencil_container_t>
void DiscretizationStencilSolver<stencil_t, solver_kernel_t>::assembly(MPI_Comm communicator, bool partitioned, const stencil_container_t &stencils)
#else
template<typename stencil_t, typename solver_kernel_t>
template<typename stencil_container_t>
void DiscretizationStencilSolver<stencil_t, solver_kernel_t>::assembly(const stencil_container_t &stencils)
#endif
{
    // Create the assembler
#if BITPIT_ENABLE_MPI==1
    DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t> assembler(communicator, partitioned, &stencils);
#else
    DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t> assembler(&stencils);
#endif
 
    // Assembly the system
    solver_kernel_t::template assembly<DiscretizationStencilSolver<stencil_t, solver_kernel_t>>(assembler, NaturalSystemMatrixOrdering());
}

template<typename stencil_t, typename solver_kernel_t>
void DiscretizationStencilSolver<stencil_t, solver_kernel_t>::assembly(const Assembler &assembler)
{
    solver_kernel_t::template assembly<DiscretizationStencilSolver<stencil_t, solver_kernel_t>>(assembler, NaturalSystemMatrixOrdering());
}

template<typename stencil_t, typename solver_kernel_t>
void DiscretizationStencilSolver<stencil_t, solver_kernel_t>::matrixAssembly(const Assembler &assembler)
{
    // Assemble matrix
    solver_kernel_t::matrixAssembly(assembler);
 
    // Assemble constants
    assembleConstants(assembler);
}

template<typename stencil_t, typename solver_kernel_t>
void DiscretizationStencilSolver<stencil_t, solver_kernel_t>::matrixUpdate(long nRows, const long *rows, const Assembler &assembler)
{
    // Update the system
    solver_kernel_t::matrixUpdate(nRows, rows, assembler);
 
    // Update the constants
    updateConstants(nRows, rows, assembler);
}

template<typename stencil_t, typename solver_kernel_t>
template<typename stencil_container_t>
void DiscretizationStencilSolver<stencil_t, solver_kernel_t>::update(const stencil_container_t &stencils)
{
    update(this->getRowCount(), nullptr, stencils);
}

template<typename stencil_t, typename solver_kernel_t>
template<typename stencil_container_t>
void DiscretizationStencilSolver<stencil_t, solver_kernel_t>::update(const std::vector<long> &rows, const stencil_container_t &stencils)
{
    update(rows.size(), rows.data(), stencils);
}

template<typename stencil_t, typename solver_kernel_t>
template<typename stencil_container_t>
void DiscretizationStencilSolver<stencil_t, solver_kernel_t>::update(std::size_t nRows, const long *rows, const stencil_container_t &stencils)
{
#if BITPIT_ENABLE_MPI==1
    DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t> assembler(this->getCommunicator(), this->isPartitioned(), &stencils);
#else
    DiscretizationStencilSolverAssembler<stencil_t, solver_kernel_t> assembler(&stencils);
#endif
 
    // Update the system
    solver_kernel_t::template update<DiscretizationStencilSolver<stencil_t, solver_kernel_t>>(nRows, rows, assembler);
}
 

template<typename stencil_t, typename solver_kernel_t>
void DiscretizationStencilSolver<stencil_t, solver_kernel_t>::update(long nRows, const long *rows, const Assembler &assembler)
{
    solver_kernel_t::template update<DiscretizationStencilSolver<stencil_t, solver_kernel_t>>(nRows, rows, assembler);
}

template<typename stencil_t, typename solver_kernel_t>
void DiscretizationStencilSolver<stencil_t, solver_kernel_t>::assembleConstants(const Assembler &assembler)
{
    long nRows = assembler.getRowCount();
    int blockSize = assembler.getBlockSize();
 
    // Create the storage
    m_constants.resize(nRows * blockSize);
 
    // Update the constants
    updateConstants(nRows, nullptr, assembler);
}

template<typename stencil_t, typename solver_kernel_t>
void DiscretizationStencilSolver<stencil_t, solver_kernel_t>::updateConstants(std::size_t nRows, const long *rows,
                                                              const Assembler &assembler)
{
    int blockSize = assembler.getBlockSize();
    ConstProxyVector<double> rowConstant(blockSize);
    for (std::size_t n = 0; n < nRows; ++n) {
        long row;
        if (rows) {
            row = rows[n];
        } else {
            row = n;
        }
 
        assembler.getRowConstant(n, &rowConstant);
        std::copy_n(rowConstant.data(), blockSize, m_constants.data() + row * blockSize);
    }
}

template<typename stencil_t, typename solver_kernel_t>
void DiscretizationStencilSolver<stencil_t, solver_kernel_t>::solve()
{
    // Check if the stencil solver is assembled
    if (!this->isAssembled()) {
        throw std::runtime_error("Unable to solve the system. The stencil solver is not yet assembled.");
    }
 
    // Subtract constant terms to the RHS
    long nUnknowns = this->getBlockSize() * this->getRowCount();
    double *raw_rhs = this->getRHSRawPtr();
    for (long i = 0; i < nUnknowns; ++i) {
        raw_rhs[i] -= m_constants[i];
    }
    this->restoreRHSRawPtr(raw_rhs);
 
    // Solve the system
    solver_kernel_t::solve();
}
 
}