support mixed fom/rom: default lspg and default implicit galerkin

include/pressio/ode/concepts/ode_policy_has_call_overload_for_userdefined_action_on_stencil.hpp

@@ -0,0 +1,78 @@
+/*
+//@HEADER
+// ************************************************************************
+//
+// ode_has_const_discrete_time_residual_method_accept_step_time_dt_result_states_return_void.hpp
+//                     		  Pressio
+//                             Copyright 2019
+//    National Technology & Engineering Solutions of Sandia, LLC (NTESS)
+//
+// Under the terms of Contract DE-NA0003525 with NTESS, the
+// U.S. Government retains certain rights in this software.
+//
+// Pressio is licensed under BSD-3-Clause terms of use:
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions
+// are met:
+//
+// 1. Redistributions of source code must retain the above copyright
+// notice, this list of conditions and the following disclaimer.
+//
+// 2. Redistributions in binary form must reproduce the above copyright
+// notice, this list of conditions and the following disclaimer in the
+// documentation and/or other materials provided with the distribution.
+//
+// 3. Neither the name of the copyright holder nor the names of its
+// contributors may be used to endorse or promote products derived
+// from this software without specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+// FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+// COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+// INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+// HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
+// STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
+// IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Questions? Contact Francesco Rizzi (fnrizzi@sandia.gov)
+//
+// ************************************************************************
+//@HEADER
+*/
+
+#ifndef ODE_POLICY_HAS_CALL_OVERLOAD_FOR_USERDEFINED_ACTION_ON_STENCIL_HPP_
+#define ODE_POLICY_HAS_CALL_OVERLOAD_FOR_USERDEFINED_ACTION_ON_STENCIL_HPP_
+
+namespace pressio{ namespace ode{
+
+template<class T>
+using policy_has_call_overload_for_userdefined_action_on_stencil_states_t =
+      decltype
+      (
+       std::declval<T const>()
+       (
+	std::declval<StencilStatesPotentiallyOverwrittenByUser>(),
+	std::declval<StepScheme const &>(),
+	std::declval<typename T::state_type const &>(),
+	std::declval<ImplicitStencilStatesDynamicContainer<typename T::state_type> const & >(),
+	std::declval<ImplicitStencilRightHandSideDynamicContainer<typename T::residual_type> & >(),
+	std::declval< ::pressio::ode::StepEndAt<typename T::independent_variable_type> >(),
+	std::declval< ::pressio::ode::StepCount >(),
+	std::declval< ::pressio::ode::StepSize<typename T::independent_variable_type> >(),
+	std::declval<typename T::residual_type &>(),
+#ifdef PRESSIO_ENABLE_CXX17
+	std::declval< std::optional<typename T::jacobian_type*> >()
+#else
+	std::declval< typename T::jacobian_type* >()
+#endif
+	)
+       );
+
+}}
+#endif

include/pressio/ode/concepts/ode_system.hpp

@@ -4,6 +4,7 @@
 
 #include "ode_predicates_for_system.hpp"
 #include "ode_has_const_discrete_residual_jacobian_method.hpp"
+#include "ode_policy_has_call_overload_for_userdefined_action_on_stencil.hpp"
 
 namespace pressio{ namespace ode{
 
@@ -338,5 +339,18 @@ struct ImplicitResidualJacobianPolicy<
     >
   > : std::true_type{};
 
+template<class T, class = void>
+struct ImplicitResidualJacobianPolicyForUserDefinedStencilStatesAction
+  : std::false_type{};
+
+template<class T>
+struct ImplicitResidualJacobianPolicyForUserDefinedStencilStatesAction<
+  T,
+  ::pressio::mpl::enable_if_t<
+    ImplicitResidualJacobianPolicy<T>::value &&
+    mpl::is_detected< policy_has_call_overload_for_userdefined_action_on_stencil_states_t >::value
+    >
+  > : std::true_type{};
+
 }}
 #endif

include/pressio/ode/concepts/ode_system_cxx20.hpp

@@ -4,6 +4,7 @@
 
 #include <concepts>
 #include "ode_has_const_discrete_residual_jacobian_method.hpp"
+#include "ode_policy_has_call_overload_for_userdefined_action_on_stencil.hpp"
 
 namespace pressio{ namespace ode{
 
@@ -194,6 +195,25 @@ concept ImplicitResidualJacobianPolicy =
       stencilVelocities, rhsEvaluationTime, stepNumber, dt, R, J);
   };
 
+
+template <class T>
+concept ImplicitResidualJacobianPolicyForUserDefinedStencilStatesAction =
+     ImplicitResidualJacobianPolicy<T>
+  && requires(const T & A,
+	      StepScheme schemeName,
+	      const typename T::state_type & predictedState,
+	      const ImplicitStencilStatesDynamicContainer<typename T::state_type> & stencilStatesManager,
+	      ImplicitStencilRightHandSideDynamicContainer<typename T::residual_type> & stencilVelocities,
+	      const ::pressio::ode::StepEndAt<typename T::independent_variable_type> & rhsEvaluationTime,
+	      ::pressio::ode::StepCount stepNumber,
+	      const ::pressio::ode::StepSize<typename T::independent_variable_type> & dt,
+	      typename T::residual_type & R,
+	      std::optional<typename T::jacobian_type *> & J)
+  {
+    A(StencilStatesPotentiallyOverwrittenByUser(), schemeName, predictedState,
+      stencilStatesManager, stencilVelocities, rhsEvaluationTime, stepNumber, dt, R, J);
+  };
+
 //
 // auxiliary stuff
 //

include/pressio/ode/impl/ode_implicit_stepper_standard.hpp

@@ -56,7 +56,7 @@ template<
   class StateType,
   class ResidualType,
   class JacobianType,
-  class ResidualJacobianPolicyType
+  class ResidualJacobianPolicyPossiblyRefType
   >
 class ImplicitStepperStandardImpl
 {
@@ -84,13 +84,15 @@ class ImplicitStepperStandardImpl
   // for cn  : y_n
   ImplicitStencilStatesDynamicContainer<StateType> stencil_states_;
 
-  ::pressio::utils::InstanceOrReferenceWrapper<ResidualJacobianPolicyType> rj_policy_;
+  ::pressio::utils::InstanceOrReferenceWrapper<ResidualJacobianPolicyPossiblyRefType> rj_policy_;
 
   // stencilRightHandSide contains:
   // for bdf1,2: nothing
   // for cn:  f(y_n,t_n) and f(y_np1, t_np1)
   mutable ImplicitStencilRightHandSideDynamicContainer<ResidualType> stencil_rhs_;
 
+  mutable bool userActionCallback_ = false;
+
 public:
   ImplicitStepperStandardImpl() = delete;
   ImplicitStepperStandardImpl(const ImplicitStepperStandardImpl & other)  = default;
@@ -99,35 +101,58 @@ class ImplicitStepperStandardImpl
 
   // *** BDF1 ***//
   ImplicitStepperStandardImpl(::pressio::ode::BDF1,
-			      ResidualJacobianPolicyType && rjPolicyObj)
+			      ResidualJacobianPolicyPossiblyRefType && rjPolicyObj)
     : name_(StepScheme::BDF1),
       recovery_state_{rjPolicyObj.createState()},
       stencil_states_{rjPolicyObj.createState()},
-      rj_policy_(std::forward<ResidualJacobianPolicyType>(rjPolicyObj))
+      rj_policy_(std::forward<ResidualJacobianPolicyPossiblyRefType>(rjPolicyObj))
   {}
 
   // *** BDF2 ***//
   ImplicitStepperStandardImpl(::pressio::ode::BDF2,
-			      ResidualJacobianPolicyType && rjPolicyObj)
+			      ResidualJacobianPolicyPossiblyRefType && rjPolicyObj)
     : name_(StepScheme::BDF2),
       recovery_state_{rjPolicyObj.createState()},
       stencil_states_{rjPolicyObj.createState(),
 		      rjPolicyObj.createState()},
-      rj_policy_(std::forward<ResidualJacobianPolicyType>(rjPolicyObj))
+      rj_policy_(std::forward<ResidualJacobianPolicyPossiblyRefType>(rjPolicyObj))
   {}
 
   // *** CN ***//
   ImplicitStepperStandardImpl(::pressio::ode::CrankNicolson,
-			      ResidualJacobianPolicyType && rjPolicyObj)
+			      ResidualJacobianPolicyPossiblyRefType && rjPolicyObj)
     : name_(StepScheme::CrankNicolson),
       recovery_state_{rjPolicyObj.createState()},
       stencil_states_{rjPolicyObj.createState()},
-      rj_policy_(std::forward<ResidualJacobianPolicyType>(rjPolicyObj)),
+      rj_policy_(std::forward<ResidualJacobianPolicyPossiblyRefType>(rjPolicyObj)),
       stencil_rhs_{rj_policy_.get().createResidual(),
                    rj_policy_.get().createResidual()}
   {}
 
 public:
+
+  template<class SolverType, class UserDefinedActionOnStencilStates, class ...SolverArgs>
+  void operator()(StateType & odeState,
+		  const ::pressio::ode::StepStartAt<independent_variable_type> & stepStartVal,
+		  ::pressio::ode::StepCount stepNumber,
+		  const ::pressio::ode::StepSize<independent_variable_type> & stepSize,
+		  UserDefinedActionOnStencilStates action,
+		  SolverType & solver,
+		  SolverArgs && ...argsForSolver)
+  {
+
+    if (name_ != ::pressio::ode::StepScheme::BDF2){
+      throw std::runtime_error("User-defined action on the stencil states currently only allowed for BDF2");
+    }
+
+    PRESSIOLOG_DEBUG("implicit stepper: do step, user-defined action on stencil states");
+    PRESSIOLOG_WARN("implicit stepper: overload accepting action on stencil states does not yet support strong guarantee");
+    userActionCallback_ = true;
+    doStepImpl(::pressio::ode::BDF2(), action, odeState,
+	       stepStartVal.get(), stepSize.get(), stepNumber.get(),
+	       solver, std::forward<SolverArgs>(argsForSolver)...);
+  }
+
   template<class SolverType, class ...SolverArgs>
   void operator()(StateType & odeState,
 		  const ::pressio::ode::StepStartAt<independent_variable_type> & stepStartVal,
@@ -138,6 +163,8 @@ class ImplicitStepperStandardImpl
   {
     PRESSIOLOG_DEBUG("implicit stepper: do step");
 
+    userActionCallback_ = false;
+
     if (name_==::pressio::ode::StepScheme::BDF1){
       doStepImpl(::pressio::ode::BDF1(),
 		 odeState, stepStartVal.get(), stepSize.get(),
@@ -147,6 +174,7 @@ class ImplicitStepperStandardImpl
 
     else if (name_==::pressio::ode::StepScheme::BDF2){
       doStepImpl(::pressio::ode::BDF2(),
+		 utils::NoOperation<void>(),
 		 odeState, stepStartVal.get(), stepSize.get(),
 		 stepNumber.get(), solver,
 		 std::forward<SolverArgs>(argsForSolver)...);
@@ -172,11 +200,29 @@ class ImplicitStepperStandardImpl
                            jacobian_type* Jo) const
 #endif
   {
-    rj_policy_.get()(name_, odeState, stencil_states_, stencil_rhs_,
-		     ::pressio::ode::StepEndAt<IndVarType>(t_np1_),
-		     ::pressio::ode::StepCount(step_number_),
-		     ::pressio::ode::StepSize<IndVarType>(dt_),
-		     R, Jo);
+    StepEndAt<IndVarType> endAt(t_np1_);
+    StepCount count(step_number_);
+    StepSize<IndVarType> stepsz(dt_);
+
+    if constexpr(mpl::is_detected<
+		 policy_has_call_overload_for_userdefined_action_on_stencil_states_t,
+		 std::remove_reference_t<ResidualJacobianPolicyPossiblyRefType>
+		 >::value)
+    {
+      if (userActionCallback_){
+	rj_policy_.get()(StencilStatesPotentiallyOverwrittenByUser(),
+			 name_, odeState, stencil_states_, stencil_rhs_,
+			 endAt, count, stepsz, R, Jo);
+      }
+      else{
+	rj_policy_.get()(name_, odeState, stencil_states_, stencil_rhs_,
+			 endAt, count, stepsz, R, Jo);
+      }
+    }
+    else{
+      rj_policy_.get()(name_, odeState, stencil_states_, stencil_rhs_,
+		       endAt, count, stepsz, R, Jo);
+    }
   }
 
 private:
@@ -189,6 +235,7 @@ class ImplicitStepperStandardImpl
 		  solver_type & solver,
 		  SolverArgs&& ...argsForSolver)
   {
+    PRESSIOLOG_DEBUG("implicit BDF1 stepper");
 
     /*
       - we are at step = stepNumber
@@ -225,15 +272,17 @@ class ImplicitStepperStandardImpl
     }
   }
 
-  template<class solver_type, class ...SolverArgs>
+  template<class solver_type, class UserDefinedActionOnStencilStates, class ...SolverArgs>
   void doStepImpl(::pressio::ode::BDF2,
+		  UserDefinedActionOnStencilStates action,
 		  state_type & odeState,
 		  const IndVarType & currentTime,
 		  const IndVarType & dt,
 		  const int32_t & stepNumber,
 		  solver_type & solver,
 		  SolverArgs&& ...argsForSolver)
   {
+    PRESSIOLOG_DEBUG("implicit BDF2 stepper");
 
     dt_ = dt;
     t_np1_ = currentTime + dt;
@@ -246,6 +295,7 @@ class ImplicitStepperStandardImpl
      */
 
     if (stepNumber == ::pressio::ode::first_step_value){
+      PRESSIOLOG_DEBUG("implicit BDF2 stepper: beginning/initial step is via BDF1");
 
       /* from t_0 to t_1 and have:
 	 odeState = the initial condition (y0)
@@ -255,6 +305,7 @@ class ImplicitStepperStandardImpl
       ::pressio::ops::deep_copy(odeState_n, odeState);
     }
     else{
+      PRESSIOLOG_DEBUG("implicit BDF2 stepper: actual BDF2 step");
 
       /* for step == 2, we are going from t_1 to t_2 and:
 	 odeState = the state at t1
@@ -268,6 +319,7 @@ class ImplicitStepperStandardImpl
       ::pressio::ops::deep_copy(odeState_n, odeState);
     }
 
+    action(stencil_states_);
     try{
       solver.solve(*this, odeState, std::forward<SolverArgs>(argsForSolver)...);
     }

include/pressio/ode/ode_enum_and_tags.hpp

@@ -96,6 +96,8 @@ class nMinusTwo{};
 class nMinusThree{};
 class nMinusFour{};
 
+struct StencilStatesPotentiallyOverwrittenByUser{};
+
 }}//end namespace pressio::ode
 
 #endif  // ODE_ODE_ENUM_AND_TAGS_HPP_