Fix snapshot container system. (idaholab#27101)

modules/stochastic_tools/include/reporters/JacobianContainer.h

@@ -10,7 +10,9 @@
 #pragma once
 
 #include "SnapshotContainerBase.h"
+#include "libmesh/id_types.h"
 #include "libmesh/petsc_vector.h"
+#include "libmesh/vectormap.h"
 
 /**
  * This class is responsible for collecting jacobian row concatenated vectors in one place. The
@@ -26,4 +28,9 @@ class JacobianContainer : public SnapshotContainerBase
 
 protected:
   virtual std::unique_ptr<NumericVector<Number>> collectSnapshot() override;
+
+  std::vector<std::pair<dof_id_type, dof_id_type>> & _sparse_ind;
+
+  NonlinearSystem & _nl_sys;
+  const TagID _tag_id;
 };

modules/stochastic_tools/include/reporters/ResidualContainer.h

@@ -9,7 +9,10 @@
 
 #pragma once
 
+#include "MooseTypes.h"
+#include "NonlinearSystemBase.h"
 #include "SnapshotContainerBase.h"
+#include "SystemBase.h"
 #include "libmesh/petsc_vector.h"
 
 /**
@@ -26,4 +29,7 @@ class ResidualContainer : public SnapshotContainerBase
 
 protected:
   virtual std::unique_ptr<NumericVector<Number>> collectSnapshot() override;
+
+  const NonlinearSystem & _nl_sys;
+  const TagID _tag_id;
 };

modules/stochastic_tools/include/reporters/SnapshotContainerBase.h

@@ -10,6 +10,7 @@
 #pragma once
 
 #include "GeneralReporter.h"
+#include "TaggingInterface.h"
 #include "libmesh/petsc_vector.h"
 #include "libmesh/petsc_matrix.h"
 #include "NonlinearSystemBase.h"

modules/stochastic_tools/include/transfers/SerializedSnapshotTransfer.h

@@ -57,8 +57,6 @@ class SerializedSnapshotTransfer : public StochasticToolsTransfer
   std::vector<SnapshotContainerBase *> _residual_container;
   /// Link to the storage spaces on the subapplications (will only hold one in batch mode)
   std::vector<SnapshotContainerBase *> _jacobian_container;
-  /// Link to the storage spaces on the subapplications (will only hold one in batch mode)
-  std::vector<SnapshotContainerBase *> _jacobian_index_container;
 
 private:
   /// Serialize on the root processor of the subapplication and transfer the result to the main application
@@ -100,6 +98,4 @@ class SerializedSnapshotTransfer : public StochasticToolsTransfer
   /// User-selected switch that determines if we want to serialize on the root of the subapp
   /// only or distribute the solutions between all the ranks of the subapp.
   const bool _serialize_on_root;
-
-  const bool _sparse_representation;
 };

modules/stochastic_tools/include/userobjects/InverseRB.h

@@ -0,0 +1,165 @@
+///* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "DEIMRBMapping.h"
+#include "GeneralUserObject.h"
+#include "MappingInterface.h"
+#include "NonlinearSystemBase.h"
+#include "SystemBase.h"
+#include "libmesh/elem.h"
+#include "libmesh/elem_range.h"
+#include "libmesh/id_types.h"
+#include <memory>
+
+/**
+ * InverseRB is a user object that performs inverse reduced basis mapping.
+ */
+class InverseRB : public GeneralUserObject, public MappingInterface
+{
+public:
+  static InputParameters validParams();
+
+  /**
+   * Constructor for InverseRB.
+   * @param parameters Input parameters
+   */
+  InverseRB(const InputParameters & parameters);
+
+  virtual void initialize() override;
+
+  virtual void execute() override;
+
+  virtual void finalize() override;
+
+  virtual void initialSetup() override;
+
+  // only needed for ElementUserObjects and NodalUseroObjects
+  virtual void threadJoin(const UserObject & /*y*/) override{};
+
+protected:
+  /// Link to the mapping object which provides the inverse mapping function
+  DEIMRBMapping * _mapping;
+
+  std::vector<dof_id_type> _residual_inds;
+
+  std::vector<std::pair<dof_id_type, dof_id_type>> _jacobian_matrix_inds;
+
+  const DofMap & _dof_map;
+
+  const dof_id_type _max_iter;
+
+  const Real _tolerance;
+
+  const Real _relax_factor;
+
+private:
+  /**
+   * Finds the reduced element range for the given DOFs.
+   * @param dofs Vector of DOF indices
+   * @return Vector of pointers to the reduced elements
+   */
+  std::vector<const Elem *> findReducedElemRange(const std::vector<dof_id_type> & dofs);
+
+  /**
+   * Finds the reduced node range for the given DOFs.
+   * @param dofs Vector of DOF indices
+   * @return Vector of pointers to the reduced nodes
+   */
+  std::vector<const Node *> findReducedNodeRange(const std::vector<dof_id_type> & dofs);
+
+  /**
+   * Creates a range object from a vector of items.
+   * @tparam RangeType Type of the range object
+   * @tparam ItemType Type of the items in the vector
+   * @tparam IteratorType Type of the iterator for the range object
+   * @param items Vector of pointers to the items
+   * @return Unique pointer to the created range object
+   */
+  template <typename RangeType, typename ItemType, typename IteratorType>
+  std::unique_ptr<RangeType> createRangeFromVector(const std::vector<const ItemType *> & items);
+
+  /**
+   * Retrieves the reduced Jacobian values from the sparse matrix.
+   * @param jac Sparse Jacobian matrix
+   * @return Vector of reduced Jacobian values
+   */
+  std::vector<Real> getReducedJacValues(const SparseMatrix<Number> & jac);
+
+  /**
+   * Retrieves the reduced residual values from the numeric vector.
+   * @param res Numeric residual vector
+   * @return Reduced residual values
+   */
+  std::vector<Real> getReducedResValues(const NumericVector<Number> & res);
+
+  /**
+   * Computes the reduced Jacobian matrix.
+   * @return Reduced Jacobian
+   */
+  DenseMatrix<Real> computeReducedJacobian();
+
+  /**
+   * Computes the reduced residual vector.
+   * @return Reduced residual
+   */
+  DenseVector<Real> computeReducedResidual();
+
+  /**
+   * Updates the solution vector with the reduced solution.
+   * @param reduced_sol Reduced solution vector
+   */
+  void updateSolution(const DenseVector<Real> & reduced_sol);
+
+  /**
+   * Computes the residual norm.
+   * @return Residual norm value
+   */
+  Real computeResidual();
+
+  /// Vector of reduced Jacobian elements
+  std::vector<const Elem *> _red_jac_elem;
+
+  /// Vector of reduced residual elements
+  std::vector<const Elem *> _red_res_elem;
+
+  /// Vector of reduced Jacobian nodes
+  std::vector<const Node *> _red_jac_node;
+
+  /// Vector of reduced residual nodes
+  std::vector<const Node *> _red_res_node;
+
+  /// Local range of reduced Jacobian elements
+  std::unique_ptr<ConstElemRange> _red_jac_elem_local_range;
+
+  /// Local range of reduced residual elements
+  std::unique_ptr<ConstElemRange> _red_res_elem_local_range;
+
+  /// Local range of reduced Jacobian nodes
+  std::unique_ptr<ConstNodeRange> _red_jac_node_local_range;
+
+  /// Local range of reduced residual nodes
+  std::unique_ptr<ConstNodeRange> _red_res_node_local_range;
+
+  /// Nonlinear system
+  NonlinearSystemBase & _nl_sys;
+
+  /// Jacobian matrix
+  SparseMatrix<Number> * _jac_matrix;
+
+  /// Residual vector
+  NumericVector<Number> * _residual;
+
+  /// Current solution
+  const NumericVector<Number> * _curr_sol;
+
+  /// Linear solver
+  LinearSolver<Number> * _solver;
+};

modules/stochastic_tools/include/utils/QDEIM.h

@@ -0,0 +1,111 @@
+///* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+#include <MooseTypes.h>
+
+// Algorithm is derived from this resource.
+// Z. Drmač and S. Gugercin, "A New Selection Operator for the Discrete Empirical Interpolation
+// Method---Improved A Priori Error Bound and Extensions," SIAM Journal on Scientific Computing,
+// vol. 38, no. 2, pp. A631-A648, 2016, doi: 10.1137/15M1019271. [Online]. Available:
+// https://doi.org/10.1137/15M1019271
+
+class QDEIM
+{
+public:
+  /**
+   * Constructor for QDEIM.
+   * @param ortho_basis_vecs Vector of orthonormal basis vectors
+   */
+  QDEIM(std::vector<DenseVector<double>> & ortho_basis_vecs);
+
+  /**
+   * Computes the selection indices using the QDEIM algorithm.
+   * @param selection_indices Vector to store the computed selection indices
+   */
+  void computeSelection(std::vector<dof_id_type> & selection_indices);
+
+  /**
+   * Computes both the selection indices and a new basis using the QDEIM algorithm.
+   * @param new_basis Pointer to a vector to store the computed new basis vectors
+   * @param selection_indices Vector to store the computed selection indices
+   */
+  void computeSelectionAndBasis(std::vector<DenseVector<double>> * new_basis,
+                                std::vector<dof_id_type> & selection_indices);
+
+private:
+  /**
+   * Builds the orthonormal matrix from the basis vectors.
+   */
+  void buildOrthoMatrix();
+
+  /**
+   * Performs QR decomposition on the orthonormal matrix.
+   */
+  void performQRDecomposition();
+
+  /**
+   * Extracts the selection indices from the QR decomposition.
+   */
+  void extractIndices();
+
+  /**
+   * Constructs the projection matrix using the selected indices.
+   */
+  void constructProjectionMatrix();
+
+  /**
+   * Creates the inverse permutation vector.
+   */
+  void createInversePermutationVector();
+
+  /**
+   * Applies the inverse permutation to the selection vector.
+   */
+  void applyInversePermutation();
+
+  /**
+   * Converts an Eigen vector to a std::vector.
+   * @param sample_indices Vector to store the converted indices
+   */
+  void convertEigenVectorToStdVector(std::vector<dof_id_type> & sample_indices);
+
+  /**
+   * Converts an Eigen matrix to a std::vector of DenseVector.
+   * @param new_basis Vector to store the converted basis vectors
+   */
+  void convertEigenMatrixToStdVector(std::vector<DenseVector<double>> & new_basis);
+
+  /// Vector of orthonormal basis vectors
+  std::vector<DenseVector<double>> & _ortho_basis_vecs;
+
+  /// Number of basis vectors
+  dof_id_type _num_basis;
+
+  /// Length of each basis vector
+  dof_id_type _length_basis;
+
+  /// Orthonormal basis matrix
+  RealEigenMatrix _ortho_basis;
+
+  /// QR decomposition object
+  Eigen::ColPivHouseholderQR<RealEigenMatrix> _qr;
+
+  /// Vector of selected indices
+  Eigen::VectorXi _indices;
+
+  /// Selection vector
+  Eigen::VectorXi _selection_vector;
+
+  /// Projection matrix
+  Eigen::MatrixXd _projection_mat;
+};