fix

Pressio · Nov 13, 2024 · fc02ae9 · fc02ae9
1 parent fbc8dc2
commit fc02ae9
Show file tree

Hide file tree

Showing 3 changed files with 89 additions and 11 deletions.
diff --git a/include/pressio/rom/impl/galerkin_steady_system_default.hpp b/include/pressio/rom/impl/galerkin_steady_system_default.hpp
@@ -72,7 +72,7 @@ class GalerkinSteadyDefaultSystem
       decltype(
 	       std::declval<_FomSystemType const>().residual(
 							     std::declval<state_type const>(),
-							     std::declval<residual_type>()
+							     std::declval<residual_type &>()
 							     )
 	       )
     >
@@ -122,9 +122,10 @@ class GalerkinSteadyDefaultSystem
     std::is_same_v<
       void,
       decltype(
-	       std::declval<_FomSystemType const>().applyJacobian(
-							     std::declval<state_type const>(),
-							     std::declval<residual_type>()
+	       std::declval<_FomSystemType const>().jacobianAction(
+							     std::declval<state_type const &>(),
+							     std::declval<state_type const &>(),
+							     std::declval<state_type &>()
 							     )
 	       )
       >
@@ -133,20 +134,22 @@ class GalerkinSteadyDefaultSystem
 		OperandT const & reducedOperand,
 		ResultT & out) const
   {
-#if 0
-    const auto & phi = trialSubspace_.get().basisOfTranslatedSpace();
     trialSubspace_.get().mapFromReducedState(reducedState, fomState_);
 
     auto auxVec = pressio::ops::clone(fomState_);
     auto auxVec2 = pressio::ops::clone(fomState_);
-    pressio::ops::set_sero(auxVec);
-    trialSubspace_.get().mapFromReducedState(reducedOperand, auxVec2_);
+    pressio::ops::set_zero(auxVec);
+    trialSubspace_.get().mapFromReducedState(reducedOperand, auxVec2);
 
-    fomSystem_.applyJacobian(fomState_, reducedOperand, auxVec);
+    fomSystem_.get().jacobianAction(fomState_, reducedOperand, auxVec);
 
+    using scalar_t = typename ::pressio::Traits<basis_matrix_type>::scalar_type;
+    scalar_t alpha{1};
+    scalar_t beta{0};
+
+    const auto & phi = trialSubspace_.get().basisOfTranslatedSpace();
     ::pressio::ops::product(::pressio::transpose(),
-			    1, phi, fomJacAction_, beta, outVec);
-#endif
+			    alpha, phi, auxVec2, beta, out);
   }
 
 private:

diff --git a/tests/functional_small/rom/CMakeLists.txt b/tests/functional_small/rom/CMakeLists.txt
@@ -29,6 +29,7 @@ endif()
 if(PRESSIO_ENABLE_TPL_EIGEN)
   set(SOURCES_GALERKIN_STEADY
     ${CMAKE_CURRENT_SOURCE_DIR}/galerkin_steady/main1.cc
+    ${CMAKE_CURRENT_SOURCE_DIR}/galerkin_steady/main2.cc
     ${CMAKE_CURRENT_SOURCE_DIR}/galerkin_steady/main3.cc
     ${CMAKE_CURRENT_SOURCE_DIR}/galerkin_steady/main4.cc)
   add_serial_utest(${TESTING_LEVEL}_rom_galerkin_steady ${SOURCES_GALERKIN_STEADY})

diff --git a/tests/functional_small/rom/galerkin_steady/main2.cc b/tests/functional_small/rom/galerkin_steady/main2.cc
@@ -0,0 +1,74 @@
+
+#include <gtest/gtest.h>
+#include "pressio/rom_subspaces.hpp"
+#include "pressio/rom_galerkin_steady.hpp"
+
+namespace{
+
+struct MyFom
+{
+  using state_type        = Eigen::VectorXd;
+  using residual_type     = state_type;
+  int N_ = {};
+
+  MyFom(int N): N_(N){}
+  residual_type createResidual() const{ return residual_type(N_); }
+
+  Eigen::MatrixXd createResultOfJacobianActionOn(const Eigen::MatrixXd & B) const{
+    Eigen::MatrixXd A(N_, B.cols());
+    return A;
+  }
+
+  void residual(const state_type & u,
+		residual_type & r) const
+  {}
+
+  template<class OperandType, class ResultType>
+  void jacobianAction(const state_type &,
+		      OperandType const &,
+		      ResultType &) const
+  {}
+
+  void residualAndJacobianAction(const state_type & u,
+				 residual_type & r,
+				 const Eigen::MatrixXd & B,
+				 std::optional<Eigen::MatrixXd *> Ain) const
+  {
+  }
+
+};
+
+}
+
+TEST(rom_galerkin_steady, default_matrix_free)
+{
+  using namespace pressio;
+
+  /*
+    just supposed to compile for now
+  */
+
+  log::initialize(logto::terminal);
+  log::setVerbosity({log::level::debug});
+
+  constexpr int N = 8;
+  using fom_t = MyFom;
+  fom_t fomSystem(N);
+
+  using phi_t = Eigen::Matrix<double, -1,-1>;
+  phi_t phi(N, 3);
+
+  using reduced_state_type = Eigen::VectorXd;
+  typename fom_t::state_type shift(N);
+  auto space = rom::create_trial_column_subspace<reduced_state_type>(phi, shift, false);
+
+  auto problem = rom::galerkin::create_steady_problem(space, fomSystem);
+
+  using tag = linearsolvers::iterative::GMRES;
+  auto nonLinSolver = experimental::create_matrixfree_newtonkrylov_solver<tag>(problem);
+
+  auto romState = space.createReducedState();
+  nonLinSolver.solve(romState);
+
+  log::finalize();
+}