wip matrix free galerkin

Pressio · Nov 13, 2024 · c41f310 · c41f310
1 parent 8c95f00
commit c41f310
Show file tree

Hide file tree

Showing 2 changed files with 58 additions and 0 deletions.
diff --git a/include/pressio/rom/galerkin_steady.hpp b/include/pressio/rom/galerkin_steady.hpp
@@ -29,6 +29,7 @@ auto create_steady_problem(const TrialSubspaceType & trialSpace,   /*(1)*/
   return return_type(trialSpace, fomSystem);
 }
 
+
 // ------------------------------------------------------------------------
 // hyper-reduced
 // ------------------------------------------------------------------------

diff --git a/include/pressio/rom/impl/galerkin_steady_system_default.hpp b/include/pressio/rom/impl/galerkin_steady_system_default.hpp
@@ -65,6 +65,30 @@ class GalerkinSteadyDefaultSystem
     return impl::CreateGalerkinJacobian<jacobian_type>()(trialSubspace_.get().dimension());
   }
 
+  template<class _FomSystemType = FomSystemType>
+  std::enable_if_t<
+    std::is_same_v<
+      void,
+      decltype(
+	       std::declval<_FomSystemType const>().residual(
+							     std::declval<state_type const>(),
+							     std::declval<residual_type>()
+							     )
+	       )
+      >
+    >
+  >
+  residual(const state_type & reducedState,
+	   residual_type & reducedResidual) const
+  {
+    const auto & phi = trialSubspace_.get().basisOfTranslatedSpace();
+    trialSubspace_.get().mapFromReducedState(reducedState, fomState_);
+
+    fomSystem_.get().residual(fomState_, fomResidual_);
+    ::pressio::ops::product(::pressio::transpose(),
+			    1, phi, fomResidual_, 0, reducedResidual);
+  }
+
   void residualAndJacobian(const state_type & reducedState,
 			   residual_type & reducedResidual,
 			   std::optional<jacobian_type*> reducedJacobian) const
@@ -94,6 +118,39 @@ class GalerkinSteadyDefaultSystem
     }
   }
 
+  template<class OperandT, class ResultT, class _FomSystemType = FomSystemType>
+  std::enable_if_t<
+    std::is_same_v<
+      void,
+      decltype(
+	       std::declval<_FomSystemType const>().applyJacobian(
+							     std::declval<state_type const>(),
+							     std::declval<residual_type>()
+							     )
+	       )
+      >
+    >
+  >
+  applyJacobian(const state_type & reducedState,
+		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_);
+
+    fomSystem_.applyJacobian(fomState_, reducedOperand, auxVec);
+
+    ::pressio::ops::product(::pressio::transpose(),
+			    1, phi, fomJacAction_, beta, outVec);
+
+  }
+
 private:
   std::reference_wrapper<const TrialSubspaceType> trialSubspace_;
   std::reference_wrapper<const FomSystemType> fomSystem_;