search

GEOS-DEV / GEOS

GEOS Simulation Framework
GNU Lesser General Public License v2.1
206 stars 83 forks source link

Redefine LAI set/insert/add interface functions...or delete them entirely #1675

Open rrsettgast opened 2 years ago

rrsettgast commented 2 years ago

Describe the issue The LAI objects have set/insert/add interface functions that take values, pointers, and arrayslice. This either requires the interface function to copy the data and move it to the correct memory space, or it requires the caller to manage the memory motion rather than allowing the interface functions to manage the memory motion as part of their responsibility. For example, MatrixBase contains the following interface functions:

  virtual void add( globalIndex const rowIndex,
                    globalIndex const colIndex,
                    real64 const value ) = 0;

  virtual void set( globalIndex const rowIndex,
                    globalIndex const colIndex,
                    real64 const value ) = 0;

  virtual void insert( globalIndex const rowIndex,
                       globalIndex const colIndex,
                       real64 const value ) = 0;

These are single point interface functions. They are very inefficient, and require that the interface function copy the data, then move it to the correct memory space. These should be removed.

virtual void add( globalIndex const rowIndex,
                    globalIndex const * colIndices,
                    real64 const * values,
                    localIndex const size ) = 0;

  virtual void set( globalIndex const rowIndex,
                    globalIndex const * colIndices,
                    real64 const * values,
                    localIndex const size ) = 0;

  virtual void insert( globalIndex const rowIndex,
                       globalIndex const * colIndices,
                       real64 const * values,
                       localIndex const size ) = 0;

  virtual void add( globalIndex const rowIndex,
                    arraySlice1d< globalIndex const > const & colIndices,
                    arraySlice1d< real64 const > const & values ) = 0;

  virtual void set( globalIndex const rowIndex,
                    arraySlice1d< globalIndex const > const & colIndices,
                    arraySlice1d< real64 const > const & values ) = 0;

  virtual void insert( globalIndex const rowIndex,
                       arraySlice1d< globalIndex const > const & colIndices,
                       arraySlice1d< real64 const > const & values ) = 0;

These are data insertions for a single row. Pointers/arraySlice make this interface hard to manage memory motion. 

  virtual void add( arraySlice1d< globalIndex const > const & rowIndices,
                    arraySlice1d< globalIndex const > const & colIndices,
                    arraySlice2d< real64 const > const & values ) override;

  virtual void set( arraySlice1d< globalIndex const > const & rowIndices,
                    arraySlice1d< globalIndex const > const & colIndices,
                    arraySlice2d< real64 const > const & values ) override;

  virtual void insert( arraySlice1d< globalIndex const > const & rowIndices,
                       arraySlice1d< globalIndex const > const & colIndices,
                       arraySlice2d< real64 const > const & values ) override;

  virtual void add( globalIndex const * rowIndices,
                    globalIndex const * colIndices,
                    real64 const * values,
                    localIndex const numRows,
                    localIndex const numCols ) override;

  virtual void set( globalIndex const * rowIndices,
                    globalIndex const * colIndices,
                    real64 const * values,
                    localIndex const numRows,
                    localIndex const numCols ) override;

  virtual void insert( globalIndex const * rowIndices,
                       globalIndex const * colIndices,
                       real64 const * values,
                       localIndex const numRows,
                       localIndex const numCols ) override;

Dense data insertions. These are not as bad since they should be launched inside a kernel...but they are not __device__ callable.

Proposed cleanup Remove all of these functions. We should just be using CRSMatrix objects and creating these objects from CRSMatrix. We can then focus on getting a shallow copy from CRSMatrix to the derived classes of MatrixBase.

joshua-white commented 2 years ago

It would be nice to better encapsulate our current assembly pattern, which is pretty raw, e.g.

https://github.com/GEOSX/GEOSX/blob/e1f643dc63297bde35e2f0eb3e2f3e394b4629ce/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsSmallStrainQuasiStaticKernel.hpp#L290-L305

https://github.com/GEOSX/GEOSX/blob/e1f643dc63297bde35e2f0eb3e2f3e394b4629ce/src/coreComponents/physicsSolvers/multiphysics/CompositionalMultiphaseReservoir.cpp#L279-L289

castelletto1 commented 2 years ago

Proposed cleanup Remove all of these functions. We should just be using CRSMatrix objects and creating these objects from CRSMatrix. We can then focus on getting a shallow copy from CRSMatrix to the derived classes of MatrixBase.

For the hypre interface I would push the cleanup even further. The proposed changes are removing the need for the IJ layer. I think we can redesign the interface manipulating directly the low level hypre_ParCSRMatrix object - which is what we already do in plenty of methods.

This method in mfem (link) is a good example for what has to change in the current create() method:

https://github.com/GEOSX/GEOSX/blob/2baf4927d912262c60bb03c36d73a13f2564473e/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMatrix.hpp#L131-L132

This effort may require a dedicated PR.