Skip to content

Commit

Permalink
Merge pull request #613 from Pressio/fixes_for_internal_use
Browse files Browse the repository at this point in the history 
Fixes for internal use
  • Loading branch information
@fnrizzi
fnrizzi authored Jun 27, 2023
2 parents e046442 + be07540 commit a6a86b6
Show file tree
Hide file tree
Showing 13 changed files with 141 additions and 161 deletions.
2 changes: 1 addition & 1 deletion include/pressio/ops/tpetra/ops_level3.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -214,7 +214,7 @@ product(::pressio::transpose /*unused*/,
const sc_t alpha_(alpha);
const sc_t beta_(beta);
const auto has_beta = beta_ != zero;
beta_t tmp = zero;
sc_t tmp = zero;

// A dot A = A^T*A, which yields a symmetric matrix
// only need to compute half and fill remaining entries accordingly
Expand Down
55 changes: 43 additions & 12 deletions include/pressio/rom/impl/galerkin_unsteady_system_fully_discrete_fom.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,14 +84,10 @@ class GalerkinDefaultFullyDiscreteSystem

const auto & ynp1 = fomStatesManager_(::pressio::ode::nPlusOne());
const auto & yn = fomStatesManager_(::pressio::ode::n());
const auto phi = trialSubspace_.get().basisOfTranslatedSpace();
const bool computeJacobian = (bool) galerkinJacobian;
try
{
fomSystem_.get().discreteTimeResidualAndJacobianAction(currentStepNumber, time_np1,
dt, fomResidual_, phi,
computeJacobian, fomJacAction_,
ynp1, yn);

try{
queryFomOperators(currentStepNumber, time_np1, dt, computeJacobian, ynp1, yn);
}
catch (::pressio::eh::DiscreteTimeResidualFailureUnrecoverable const & e){
throw ::pressio::eh::ResidualEvaluationFailureUnrecoverable();
Expand Down Expand Up @@ -121,13 +117,10 @@ class GalerkinDefaultFullyDiscreteSystem
const auto & ynp1 = fomStatesManager_(::pressio::ode::nPlusOne());
const auto & yn = fomStatesManager_(::pressio::ode::n());
const auto & ynm1 = fomStatesManager_(::pressio::ode::nMinusOne());
const auto phi = trialSubspace_.get().basisOfTranslatedSpace();
const bool computeJacobian = bool(galerkinJacobian);

try{
fomSystem_.get().discreteTimeResidualAndJacobianAction(currentStepNumber, time_np1,
dt, fomResidual_, phi,
computeJacobian, fomJacAction_,
ynp1, yn, ynm1);
queryFomOperators(currentStepNumber, time_np1, dt, computeJacobian, ynp1, yn, ynm1);
}
catch (::pressio::eh::DiscreteTimeResidualFailureUnrecoverable const & e){
throw ::pressio::eh::ResidualEvaluationFailureUnrecoverable();
Expand All @@ -136,7 +129,45 @@ class GalerkinDefaultFullyDiscreteSystem
computeReducedOperators(galerkinResidual, galerkinJacobian);
}


private:
template<typename step_t, class ...States>
void queryFomOperators(const step_t & currentStepNumber,
const independent_variable_type & time_np1,
const independent_variable_type & dt,
bool computeJacobian,
States && ... states) const
{
const auto phi = trialSubspace_.get().basisOfTranslatedSpace();

#ifdef PRESSIO_ENABLE_CXX17
if (computeJacobian){
using op_ja_t = std::optional<fom_jac_action_result_type*>;
fomSystem_.get().discreteTimeResidualAndJacobianAction(currentStepNumber, time_np1,
dt, fomResidual_, phi,
op_ja_t{&fomJacAction_},
std::forward<States>(states)...);
}
else{
fomSystem_.get().discreteTimeResidualAndJacobianAction(currentStepNumber, time_np1,
dt, fomResidual_, phi, {},
std::forward<States>(states)...);
}
#else
if (computeJacobian){
fomSystem_.get().discreteTimeResidualAndJacobianAction(currentStepNumber, time_np1,
dt, fomResidual_, phi, &fomJacAction_,
std::forward<States>(states)...);
}
else{
fomSystem_.get().discreteTimeResidualAndJacobianAction(currentStepNumber, time_np1,
dt, fomResidual_, phi, nullptr,
std::forward<States>(states)...);
}
#endif

}


void computeReducedOperators(discrete_residual_type & galerkinResidual,
#ifdef PRESSIO_ENABLE_CXX17
Expand Down
18 changes: 2 additions & 16 deletions include/pressio/rom/impl/lspg_unsteady_fully_discrete_system.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -118,18 +118,11 @@ class LspgFullyDiscreteSystem
const auto & ynp1 = fomStatesManager_(::pressio::ode::nPlusOne());
const auto & yn = fomStatesManager_(::pressio::ode::n());
const auto phi = trialSubspace_.get().basisOfTranslatedSpace();
const bool computeJacobian = bool(Jo);

try
{
fomSystem_.get().discreteTimeResidualAndJacobianAction(currentStepNumber, time_np1, dt,
R, phi, computeJacobian,
#ifdef PRESSIO_ENABLE_CXX17
*Jo.value(),
#else
*Jo,
#endif
ynp1, yn);
R, phi, Jo, ynp1, yn);
}
catch (::pressio::eh::DiscreteTimeResidualFailureUnrecoverable const & e){
throw ::pressio::eh::ResidualEvaluationFailureUnrecoverable();
Expand Down Expand Up @@ -157,17 +150,10 @@ class LspgFullyDiscreteSystem
const auto & yn = fomStatesManager_(::pressio::ode::n());
const auto & ynm1 = fomStatesManager_(::pressio::ode::nMinusOne());
const auto phi = trialSubspace_.get().basisOfTranslatedSpace();
const bool computeJacobian = bool(Jo);

try{
fomSystem_.get().discreteTimeResidualAndJacobianAction(currentStepNumber, time_np1, dt,
R, phi, computeJacobian,
#ifdef PRESSIO_ENABLE_CXX17
*Jo.value(),
#else
*Jo,
#endif
ynp1, yn, ynm1);
R, phi, Jo, ynp1, yn, ynm1);
}
catch (::pressio::eh::DiscreteTimeResidualFailureUnrecoverable const & e){
throw ::pressio::eh::ResidualEvaluationFailureUnrecoverable();
Expand Down
34 changes: 21 additions & 13 deletions include/pressio/rom/impl/ode_has_const_discrete_residual_jacobian_action_method.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,14 +28,17 @@ struct has_const_discrete_residual_jacobian_action_method<
(
std::declval<T const>().discreteTimeResidualAndJacobianAction
(
std::declval<StepType const &>(),
std::declval<IndVarType const &>(),
std::declval<IndVarType const &>(),
std::declval<ResidualType &>(),
std::declval<ManifoldJacobian const &>(),
std::declval<bool>(),
std::declval<JacobianType &>(),
std::declval<StateType const&>()
std::declval<StepType const &>(),
std::declval<IndVarType const &>(),
std::declval<IndVarType const &>(),
std::declval<ResidualType &>(),
std::declval<ManifoldJacobian const &>(),
#ifdef PRESSIO_ENABLE_CXX17
std::declval< std::optional<JacobianType*> >(),
#else
std::declval<JacobianType*>(),
#endif
std::declval<StateType const&>()
)
)
>::value
Expand All @@ -58,8 +61,11 @@ struct has_const_discrete_residual_jacobian_action_method<
std::declval<IndVarType const &>(),
std::declval<ResidualType &>(),
std::declval<ManifoldJacobian const &>(),
std::declval<bool>(),
std::declval<JacobianType &>(),
#ifdef PRESSIO_ENABLE_CXX17
std::declval< std::optional<JacobianType*> >(),
#else
std::declval<JacobianType*>(),
#endif
std::declval<StateType const&>(),
std::declval<StateType const&>()
)
Expand All @@ -84,8 +90,11 @@ struct has_const_discrete_residual_jacobian_action_method<
std::declval<IndVarType const &>(),
std::declval<ResidualType &>(),
std::declval<ManifoldJacobian const &>(),
std::declval<bool>(),
std::declval<JacobianType &>(),
#ifdef PRESSIO_ENABLE_CXX17
std::declval< std::optional<JacobianType*> >(),
#else
std::declval<JacobianType*>(),
#endif
std::declval<StateType const&>(),
std::declval<StateType const&>(),
std::declval<StateType const&>()
Expand All @@ -95,6 +104,5 @@ struct has_const_discrete_residual_jacobian_action_method<
>
> : std::true_type{};


}}
#endif // ROM_PREDICATES_HPP_
95 changes: 2 additions & 93 deletions include/pressio/rom/predicates.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,8 @@
#ifndef ROM_PREDICATES_HPP_
#define ROM_PREDICATES_HPP_

#include "./impl/ode_has_const_discrete_residual_jacobian_action_method.hpp"

namespace pressio{ namespace rom{

#define PRESSIO_ROM_IMPL_HAS_CONST_CREATE_RETURN_RESULT(NAMEA, NAMEB) \
Expand Down Expand Up @@ -247,99 +249,6 @@ struct has_const_create_apply_mass_matrix_result_method_accept_operand_return_re
> : std::true_type{};
// ---------------------------------------------------------------


template <
class T,
int n,
class StepType,
class IndVarType,
class StateType,
class ResidualType,
class ManifoldJacobian,
class JacobianType,
class = void
>
struct has_const_discrete_residual_jacobian_action_method : std::false_type{};

template <
class T, class StepType, class IndVarType, class StateType,
class ResidualType, class ManifoldJacobian, class JacobianType>
struct has_const_discrete_residual_jacobian_action_method<
T, 1, StepType, IndVarType, StateType, ResidualType, ManifoldJacobian, JacobianType,
::pressio::mpl::enable_if_t<
std::is_void<
decltype
(
std::declval<T const>().discreteTimeResidualAndJacobianAction
(
std::declval<StepType const &>(),
std::declval<IndVarType const &>(),
std::declval<IndVarType const &>(),
std::declval<ResidualType &>(),
std::declval<ManifoldJacobian const &>(),
std::declval<bool>(),
std::declval<JacobianType &>(),
std::declval<StateType const&>()
)
)
>::value
>
> : std::true_type{};

template <
class T, class StepType, class IndVarType, class StateType,
class ResidualType, class ManifoldJacobian, class JacobianType>
struct has_const_discrete_residual_jacobian_action_method<
T, 2, StepType, IndVarType, StateType, ResidualType, ManifoldJacobian, JacobianType,
::pressio::mpl::enable_if_t<
std::is_void<
decltype
(
std::declval<T const>().discreteTimeResidualAndJacobianAction
(
std::declval<StepType const &>(),
std::declval<IndVarType const &>(),
std::declval<IndVarType const &>(),
std::declval<ResidualType &>(),
std::declval<ManifoldJacobian const &>(),
std::declval<bool>(),
std::declval<JacobianType &>(),
std::declval<StateType const&>(),
std::declval<StateType const&>()
)
)
>::value
>
> : std::true_type{};

template <
class T, class StepType, class IndVarType, class StateType,
class ResidualType, class ManifoldJacobian, class JacobianType>
struct has_const_discrete_residual_jacobian_action_method<
T, 3, StepType, IndVarType, StateType, ResidualType, ManifoldJacobian, JacobianType,
::pressio::mpl::enable_if_t<
std::is_void<
decltype
(
std::declval<T const>().discreteTimeResidualAndJacobianAction
(
std::declval<StepType const &>(),
std::declval<IndVarType const &>(),
std::declval<IndVarType const &>(),
std::declval<ResidualType &>(),
std::declval<ManifoldJacobian const &>(),
std::declval<bool>(),
std::declval<JacobianType &>(),
std::declval<StateType const&>(),
std::declval<StateType const&>(),
std::declval<StateType const&>()
)
)
>::value
>
> : std::true_type{};
// ---------------------------------------------------------------

template <class T, class StateType, class IndVarType, class RhsType, class = void>
struct has_const_rhs_method_accept_state_indvar_result_return_void
: std::false_type{};
Expand Down
14 changes: 7 additions & 7 deletions include/pressio/solvers_nonlinear/impl/nonlinear_least_squares.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -135,13 +135,13 @@ class NonLinLeastSquares : public RegistryType
// currently we don't have the diagonostics stuff all flushed out
// so we limit it to work for a specific case
const auto & publicDiags = diagnostics_.publicNames();
assert(publicDiags.size() == 6);
assert(publicDiags[0] == Diagnostic::objectiveAbsolute);
assert(publicDiags[1] == Diagnostic::objectiveRelative);
assert(publicDiags[2] == Diagnostic::correctionAbsolutel2Norm);
assert(publicDiags[3] == Diagnostic::correctionRelativel2Norm);
assert(publicDiags[4] == Diagnostic::gradientAbsolutel2Norm);
assert(publicDiags[5] == Diagnostic::gradientRelativel2Norm);
// assert(publicDiags.size() == 6);
// assert(publicDiags[0] == Diagnostic::objectiveAbsolute);
// assert(publicDiags[1] == Diagnostic::objectiveRelative);
// assert(publicDiags[2] == Diagnostic::correctionAbsolutel2Norm);
// assert(publicDiags[3] == Diagnostic::correctionRelativel2Norm);
// assert(publicDiags[4] == Diagnostic::gradientAbsolutel2Norm);
// assert(publicDiags[5] == Diagnostic::gradientRelativel2Norm);

// check that the stopping criterion uses a metric already
// supported in the diagonstics, otherwise add it
Expand Down
6 changes: 6 additions & 0 deletions include/pressio/solvers_nonlinear/solvers_create_gauss_newton.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,6 +98,8 @@ auto create_gauss_newton_solver(const SystemType & system, /*(1)*/
const std::vector<Diagnostic> defaultDiagnostics =
{Diagnostic::objectiveAbsolute,
Diagnostic::objectiveRelative,
Diagnostic::residualAbsolutel2Norm,
Diagnostic::residualRelativel2Norm,
Diagnostic::correctionAbsolutel2Norm,
Diagnostic::correctionRelativel2Norm,
Diagnostic::gradientAbsolutel2Norm,
Expand Down Expand Up @@ -159,6 +161,8 @@ auto create_gauss_newton_solver(const SystemType & system, /*(2)*/
const std::vector<Diagnostic> defaultDiagnostics =
{Diagnostic::objectiveAbsolute,
Diagnostic::objectiveRelative,
Diagnostic::residualAbsolutel2Norm,
Diagnostic::residualRelativel2Norm,
Diagnostic::correctionAbsolutel2Norm,
Diagnostic::correctionRelativel2Norm,
Diagnostic::gradientAbsolutel2Norm,
Expand Down Expand Up @@ -217,6 +221,8 @@ auto create_gauss_newton_qr_solver(const SystemType & system,
const std::vector<Diagnostic> defaultDiagnostics =
{Diagnostic::objectiveAbsolute,
Diagnostic::objectiveRelative,
Diagnostic::residualAbsolutel2Norm,
Diagnostic::residualRelativel2Norm,
Diagnostic::correctionAbsolutel2Norm,
Diagnostic::correctionRelativel2Norm,
Diagnostic::gradientAbsolutel2Norm,
Expand Down
3 changes: 2 additions & 1 deletion include/pressio/solvers_nonlinear/solvers_create_levenberg_marquardt.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -100,12 +100,13 @@ auto create_levenberg_marquardt_solver(const SystemType & system,
const std::vector<Diagnostic> defaultDiagnostics =
{Diagnostic::objectiveAbsolute,
Diagnostic::objectiveRelative,
Diagnostic::residualAbsolutel2Norm,
Diagnostic::residualRelativel2Norm,
Diagnostic::correctionAbsolutel2Norm,
Diagnostic::correctionRelativel2Norm,
Diagnostic::gradientAbsolutel2Norm,
Diagnostic::gradientRelativel2Norm};


using tag = nonlinearsolvers::impl::LevenbergMarquardtNormalEqTag;
using state_t = typename SystemType::state_type;
using reg_t = nonlinearsolvers::impl::RegistryLevMarNormalEqs<SystemType, LinearSolverType>;
Expand Down
19 changes: 15 additions & 4 deletions tests/functional_small/rom/galerkin_unsteady_implicit/main4.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -218,8 +218,11 @@ class MyFom
double dt,
discrete_residual_type & R,
const phi_type & B,
bool computeJacobian,
phi_type & JA,
#ifdef PRESSIO_ENABLE_CXX17
std::optional<phi_type*> JA,
#else
phi_type* JA,
#endif
const state_type & y_np1,
const state_type & y_n ) const
{
Expand All @@ -231,11 +234,19 @@ class MyFom

R = y_np1 - y_n - dt*f;

if (computeJacobian){
#ifdef PRESSIO_ENABLE_CXX17
if (bool(JA)){
auto appJac = B;
appJac.array() += time;
JA = (B - dt*appJac);
*JA.value() = (B - dt*appJac);
}
#else
if (JA){
auto appJac = B;
appJac.array() += time;
*JA = (B - dt*appJac);
}
#endif
}
};
}
Expand Down
Loading

0 comments on commit a6a86b6

Please sign in to comment.