Skip to content

Commit a6a86b6

Browse files
author
Francesco Rizzi
authored
Merge pull request #613 from Pressio/fixes_for_internal_use
Fixes for internal use
2 parents e046442 + be07540 commit a6a86b6

File tree

13 files changed

+141
-161
lines changed

13 files changed

+141
-161
lines changed

‎include/pressio/ops/tpetra/ops_level3.hpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -214,7 +214,7 @@ product(::pressio::transpose /*unused*/,
214214
const sc_t alpha_(alpha);
215215
const sc_t beta_(beta);
216216
const auto has_beta = beta_ != zero;
217-
beta_t tmp = zero;
217+
sc_t tmp = zero;
218218

219219
// A dot A = A^T*A, which yields a symmetric matrix
220220
// only need to compute half and fill remaining entries accordingly

‎include/pressio/rom/impl/galerkin_unsteady_system_fully_discrete_fom.hpp‎

Lines changed: 43 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -84,14 +84,10 @@ class GalerkinDefaultFullyDiscreteSystem
8484

8585
const auto & ynp1 = fomStatesManager_(::pressio::ode::nPlusOne());
8686
const auto & yn = fomStatesManager_(::pressio::ode::n());
87-
const auto phi = trialSubspace_.get().basisOfTranslatedSpace();
8887
const bool computeJacobian = (bool) galerkinJacobian;
89-
try
90-
{
91-
fomSystem_.get().discreteTimeResidualAndJacobianAction(currentStepNumber, time_np1,
92-
dt, fomResidual_, phi,
93-
computeJacobian, fomJacAction_,
94-
ynp1, yn);
88+
89+
try{
90+
queryFomOperators(currentStepNumber, time_np1, dt, computeJacobian, ynp1, yn);
9591
}
9692
catch (::pressio::eh::DiscreteTimeResidualFailureUnrecoverable const & e){
9793
throw ::pressio::eh::ResidualEvaluationFailureUnrecoverable();
@@ -121,13 +117,10 @@ class GalerkinDefaultFullyDiscreteSystem
121117
const auto & ynp1 = fomStatesManager_(::pressio::ode::nPlusOne());
122118
const auto & yn = fomStatesManager_(::pressio::ode::n());
123119
const auto & ynm1 = fomStatesManager_(::pressio::ode::nMinusOne());
124-
const auto phi = trialSubspace_.get().basisOfTranslatedSpace();
125120
const bool computeJacobian = bool(galerkinJacobian);
121+
126122
try{
127-
fomSystem_.get().discreteTimeResidualAndJacobianAction(currentStepNumber, time_np1,
128-
dt, fomResidual_, phi,
129-
computeJacobian, fomJacAction_,
130-
ynp1, yn, ynm1);
123+
queryFomOperators(currentStepNumber, time_np1, dt, computeJacobian, ynp1, yn, ynm1);
131124
}
132125
catch (::pressio::eh::DiscreteTimeResidualFailureUnrecoverable const & e){
133126
throw ::pressio::eh::ResidualEvaluationFailureUnrecoverable();
@@ -136,7 +129,45 @@ class GalerkinDefaultFullyDiscreteSystem
136129
computeReducedOperators(galerkinResidual, galerkinJacobian);
137130
}
138131

132+
139133
private:
134+
template<typename step_t, class ...States>
135+
void queryFomOperators(const step_t & currentStepNumber,
136+
const independent_variable_type & time_np1,
137+
const independent_variable_type & dt,
138+
bool computeJacobian,
139+
States && ... states) const
140+
{
141+
const auto phi = trialSubspace_.get().basisOfTranslatedSpace();
142+
143+
#ifdef PRESSIO_ENABLE_CXX17
144+
if (computeJacobian){
145+
using op_ja_t = std::optional<fom_jac_action_result_type*>;
146+
fomSystem_.get().discreteTimeResidualAndJacobianAction(currentStepNumber, time_np1,
147+
dt, fomResidual_, phi,
148+
op_ja_t{&fomJacAction_},
149+
std::forward<States>(states)...);
150+
}
151+
else{
152+
fomSystem_.get().discreteTimeResidualAndJacobianAction(currentStepNumber, time_np1,
153+
dt, fomResidual_, phi, {},
154+
std::forward<States>(states)...);
155+
}
156+
#else
157+
if (computeJacobian){
158+
fomSystem_.get().discreteTimeResidualAndJacobianAction(currentStepNumber, time_np1,
159+
dt, fomResidual_, phi, &fomJacAction_,
160+
std::forward<States>(states)...);
161+
}
162+
else{
163+
fomSystem_.get().discreteTimeResidualAndJacobianAction(currentStepNumber, time_np1,
164+
dt, fomResidual_, phi, nullptr,
165+
std::forward<States>(states)...);
166+
}
167+
#endif
168+
169+
}
170+
140171

141172
void computeReducedOperators(discrete_residual_type & galerkinResidual,
142173
#ifdef PRESSIO_ENABLE_CXX17

‎include/pressio/rom/impl/lspg_unsteady_fully_discrete_system.hpp‎

Lines changed: 2 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -118,18 +118,11 @@ class LspgFullyDiscreteSystem
118118
const auto & ynp1 = fomStatesManager_(::pressio::ode::nPlusOne());
119119
const auto & yn = fomStatesManager_(::pressio::ode::n());
120120
const auto phi = trialSubspace_.get().basisOfTranslatedSpace();
121-
const bool computeJacobian = bool(Jo);
122121

123122
try
124123
{
125124
fomSystem_.get().discreteTimeResidualAndJacobianAction(currentStepNumber, time_np1, dt,
126-
R, phi, computeJacobian,
127-
#ifdef PRESSIO_ENABLE_CXX17
128-
*Jo.value(),
129-
#else
130-
*Jo,
131-
#endif
132-
ynp1, yn);
125+
R, phi, Jo, ynp1, yn);
133126
}
134127
catch (::pressio::eh::DiscreteTimeResidualFailureUnrecoverable const & e){
135128
throw ::pressio::eh::ResidualEvaluationFailureUnrecoverable();
@@ -157,17 +150,10 @@ class LspgFullyDiscreteSystem
157150
const auto & yn = fomStatesManager_(::pressio::ode::n());
158151
const auto & ynm1 = fomStatesManager_(::pressio::ode::nMinusOne());
159152
const auto phi = trialSubspace_.get().basisOfTranslatedSpace();
160-
const bool computeJacobian = bool(Jo);
161153

162154
try{
163155
fomSystem_.get().discreteTimeResidualAndJacobianAction(currentStepNumber, time_np1, dt,
164-
R, phi, computeJacobian,
165-
#ifdef PRESSIO_ENABLE_CXX17
166-
*Jo.value(),
167-
#else
168-
*Jo,
169-
#endif
170-
ynp1, yn, ynm1);
156+
R, phi, Jo, ynp1, yn, ynm1);
171157
}
172158
catch (::pressio::eh::DiscreteTimeResidualFailureUnrecoverable const & e){
173159
throw ::pressio::eh::ResidualEvaluationFailureUnrecoverable();

‎include/pressio/rom/impl/ode_has_const_discrete_residual_jacobian_action_method.hpp‎

Lines changed: 21 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -28,14 +28,17 @@ struct has_const_discrete_residual_jacobian_action_method<
2828
(
2929
std::declval<T const>().discreteTimeResidualAndJacobianAction
3030
(
31-
std::declval<StepType const &>(),
32-
std::declval<IndVarType const &>(),
33-
std::declval<IndVarType const &>(),
34-
std::declval<ResidualType &>(),
35-
std::declval<ManifoldJacobian const &>(),
36-
std::declval<bool>(),
37-
std::declval<JacobianType &>(),
38-
std::declval<StateType const&>()
31+
std::declval<StepType const &>(),
32+
std::declval<IndVarType const &>(),
33+
std::declval<IndVarType const &>(),
34+
std::declval<ResidualType &>(),
35+
std::declval<ManifoldJacobian const &>(),
36+
#ifdef PRESSIO_ENABLE_CXX17
37+
std::declval< std::optional<JacobianType*> >(),
38+
#else
39+
std::declval<JacobianType*>(),
40+
#endif
41+
std::declval<StateType const&>()
3942
)
4043
)
4144
>::value
@@ -58,8 +61,11 @@ struct has_const_discrete_residual_jacobian_action_method<
5861
std::declval<IndVarType const &>(),
5962
std::declval<ResidualType &>(),
6063
std::declval<ManifoldJacobian const &>(),
61-
std::declval<bool>(),
62-
std::declval<JacobianType &>(),
64+
#ifdef PRESSIO_ENABLE_CXX17
65+
std::declval< std::optional<JacobianType*> >(),
66+
#else
67+
std::declval<JacobianType*>(),
68+
#endif
6369
std::declval<StateType const&>(),
6470
std::declval<StateType const&>()
6571
)
@@ -84,8 +90,11 @@ struct has_const_discrete_residual_jacobian_action_method<
8490
std::declval<IndVarType const &>(),
8591
std::declval<ResidualType &>(),
8692
std::declval<ManifoldJacobian const &>(),
87-
std::declval<bool>(),
88-
std::declval<JacobianType &>(),
93+
#ifdef PRESSIO_ENABLE_CXX17
94+
std::declval< std::optional<JacobianType*> >(),
95+
#else
96+
std::declval<JacobianType*>(),
97+
#endif
8998
std::declval<StateType const&>(),
9099
std::declval<StateType const&>(),
91100
std::declval<StateType const&>()
@@ -95,6 +104,5 @@ struct has_const_discrete_residual_jacobian_action_method<
95104
>
96105
> : std::true_type{};
97106

98-
99107
}}
100108
#endif // ROM_PREDICATES_HPP_

‎include/pressio/rom/predicates.hpp‎

Lines changed: 2 additions & 93 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,8 @@
22
#ifndef ROM_PREDICATES_HPP_
33
#define ROM_PREDICATES_HPP_
44

5+
#include "./impl/ode_has_const_discrete_residual_jacobian_action_method.hpp"
6+
57
namespace pressio{ namespace rom{
68

79
#define PRESSIO_ROM_IMPL_HAS_CONST_CREATE_RETURN_RESULT(NAMEA, NAMEB) \
@@ -247,99 +249,6 @@ struct has_const_create_apply_mass_matrix_result_method_accept_operand_return_re
247249
> : std::true_type{};
248250
// ---------------------------------------------------------------
249251

250-
251-
template <
252-
class T,
253-
int n,
254-
class StepType,
255-
class IndVarType,
256-
class StateType,
257-
class ResidualType,
258-
class ManifoldJacobian,
259-
class JacobianType,
260-
class = void
261-
>
262-
struct has_const_discrete_residual_jacobian_action_method : std::false_type{};
263-
264-
template <
265-
class T, class StepType, class IndVarType, class StateType,
266-
class ResidualType, class ManifoldJacobian, class JacobianType>
267-
struct has_const_discrete_residual_jacobian_action_method<
268-
T, 1, StepType, IndVarType, StateType, ResidualType, ManifoldJacobian, JacobianType,
269-
::pressio::mpl::enable_if_t<
270-
std::is_void<
271-
decltype
272-
(
273-
std::declval<T const>().discreteTimeResidualAndJacobianAction
274-
(
275-
std::declval<StepType const &>(),
276-
std::declval<IndVarType const &>(),
277-
std::declval<IndVarType const &>(),
278-
std::declval<ResidualType &>(),
279-
std::declval<ManifoldJacobian const &>(),
280-
std::declval<bool>(),
281-
std::declval<JacobianType &>(),
282-
std::declval<StateType const&>()
283-
)
284-
)
285-
>::value
286-
>
287-
> : std::true_type{};
288-
289-
template <
290-
class T, class StepType, class IndVarType, class StateType,
291-
class ResidualType, class ManifoldJacobian, class JacobianType>
292-
struct has_const_discrete_residual_jacobian_action_method<
293-
T, 2, StepType, IndVarType, StateType, ResidualType, ManifoldJacobian, JacobianType,
294-
::pressio::mpl::enable_if_t<
295-
std::is_void<
296-
decltype
297-
(
298-
std::declval<T const>().discreteTimeResidualAndJacobianAction
299-
(
300-
std::declval<StepType const &>(),
301-
std::declval<IndVarType const &>(),
302-
std::declval<IndVarType const &>(),
303-
std::declval<ResidualType &>(),
304-
std::declval<ManifoldJacobian const &>(),
305-
std::declval<bool>(),
306-
std::declval<JacobianType &>(),
307-
std::declval<StateType const&>(),
308-
std::declval<StateType const&>()
309-
)
310-
)
311-
>::value
312-
>
313-
> : std::true_type{};
314-
315-
template <
316-
class T, class StepType, class IndVarType, class StateType,
317-
class ResidualType, class ManifoldJacobian, class JacobianType>
318-
struct has_const_discrete_residual_jacobian_action_method<
319-
T, 3, StepType, IndVarType, StateType, ResidualType, ManifoldJacobian, JacobianType,
320-
::pressio::mpl::enable_if_t<
321-
std::is_void<
322-
decltype
323-
(
324-
std::declval<T const>().discreteTimeResidualAndJacobianAction
325-
(
326-
std::declval<StepType const &>(),
327-
std::declval<IndVarType const &>(),
328-
std::declval<IndVarType const &>(),
329-
std::declval<ResidualType &>(),
330-
std::declval<ManifoldJacobian const &>(),
331-
std::declval<bool>(),
332-
std::declval<JacobianType &>(),
333-
std::declval<StateType const&>(),
334-
std::declval<StateType const&>(),
335-
std::declval<StateType const&>()
336-
)
337-
)
338-
>::value
339-
>
340-
> : std::true_type{};
341-
// ---------------------------------------------------------------
342-
343252
template <class T, class StateType, class IndVarType, class RhsType, class = void>
344253
struct has_const_rhs_method_accept_state_indvar_result_return_void
345254
: std::false_type{};

‎include/pressio/solvers_nonlinear/impl/nonlinear_least_squares.hpp‎

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -135,13 +135,13 @@ class NonLinLeastSquares : public RegistryType
135135
// currently we don't have the diagonostics stuff all flushed out
136136
// so we limit it to work for a specific case
137137
const auto & publicDiags = diagnostics_.publicNames();
138-
assert(publicDiags.size() == 6);
139-
assert(publicDiags[0] == Diagnostic::objectiveAbsolute);
140-
assert(publicDiags[1] == Diagnostic::objectiveRelative);
141-
assert(publicDiags[2] == Diagnostic::correctionAbsolutel2Norm);
142-
assert(publicDiags[3] == Diagnostic::correctionRelativel2Norm);
143-
assert(publicDiags[4] == Diagnostic::gradientAbsolutel2Norm);
144-
assert(publicDiags[5] == Diagnostic::gradientRelativel2Norm);
138+
// assert(publicDiags.size() == 6);
139+
// assert(publicDiags[0] == Diagnostic::objectiveAbsolute);
140+
// assert(publicDiags[1] == Diagnostic::objectiveRelative);
141+
// assert(publicDiags[2] == Diagnostic::correctionAbsolutel2Norm);
142+
// assert(publicDiags[3] == Diagnostic::correctionRelativel2Norm);
143+
// assert(publicDiags[4] == Diagnostic::gradientAbsolutel2Norm);
144+
// assert(publicDiags[5] == Diagnostic::gradientRelativel2Norm);
145145

146146
// check that the stopping criterion uses a metric already
147147
// supported in the diagonstics, otherwise add it

‎include/pressio/solvers_nonlinear/solvers_create_gauss_newton.hpp‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -98,6 +98,8 @@ auto create_gauss_newton_solver(const SystemType & system, /*(1)*/
9898
const std::vector<Diagnostic> defaultDiagnostics =
9999
{Diagnostic::objectiveAbsolute,
100100
Diagnostic::objectiveRelative,
101+
Diagnostic::residualAbsolutel2Norm,
102+
Diagnostic::residualRelativel2Norm,
101103
Diagnostic::correctionAbsolutel2Norm,
102104
Diagnostic::correctionRelativel2Norm,
103105
Diagnostic::gradientAbsolutel2Norm,
@@ -159,6 +161,8 @@ auto create_gauss_newton_solver(const SystemType & system, /*(2)*/
159161
const std::vector<Diagnostic> defaultDiagnostics =
160162
{Diagnostic::objectiveAbsolute,
161163
Diagnostic::objectiveRelative,
164+
Diagnostic::residualAbsolutel2Norm,
165+
Diagnostic::residualRelativel2Norm,
162166
Diagnostic::correctionAbsolutel2Norm,
163167
Diagnostic::correctionRelativel2Norm,
164168
Diagnostic::gradientAbsolutel2Norm,
@@ -217,6 +221,8 @@ auto create_gauss_newton_qr_solver(const SystemType & system,
217221
const std::vector<Diagnostic> defaultDiagnostics =
218222
{Diagnostic::objectiveAbsolute,
219223
Diagnostic::objectiveRelative,
224+
Diagnostic::residualAbsolutel2Norm,
225+
Diagnostic::residualRelativel2Norm,
220226
Diagnostic::correctionAbsolutel2Norm,
221227
Diagnostic::correctionRelativel2Norm,
222228
Diagnostic::gradientAbsolutel2Norm,

‎include/pressio/solvers_nonlinear/solvers_create_levenberg_marquardt.hpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -100,12 +100,13 @@ auto create_levenberg_marquardt_solver(const SystemType & system,
100100
const std::vector<Diagnostic> defaultDiagnostics =
101101
{Diagnostic::objectiveAbsolute,
102102
Diagnostic::objectiveRelative,
103+
Diagnostic::residualAbsolutel2Norm,
104+
Diagnostic::residualRelativel2Norm,
103105
Diagnostic::correctionAbsolutel2Norm,
104106
Diagnostic::correctionRelativel2Norm,
105107
Diagnostic::gradientAbsolutel2Norm,
106108
Diagnostic::gradientRelativel2Norm};
107109

108-
109110
using tag = nonlinearsolvers::impl::LevenbergMarquardtNormalEqTag;
110111
using state_t = typename SystemType::state_type;
111112
using reg_t = nonlinearsolvers::impl::RegistryLevMarNormalEqs<SystemType, LinearSolverType>;

‎tests/functional_small/rom/galerkin_unsteady_implicit/main4.cc‎

Lines changed: 15 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -218,8 +218,11 @@ class MyFom
218218
double dt,
219219
discrete_residual_type & R,
220220
const phi_type & B,
221-
bool computeJacobian,
222-
phi_type & JA,
221+
#ifdef PRESSIO_ENABLE_CXX17
222+
std::optional<phi_type*> JA,
223+
#else
224+
phi_type* JA,
225+
#endif
223226
const state_type & y_np1,
224227
const state_type & y_n ) const
225228
{
@@ -231,11 +234,19 @@ class MyFom
231234

232235
R = y_np1 - y_n - dt*f;
233236

234-
if (computeJacobian){
237+
#ifdef PRESSIO_ENABLE_CXX17
238+
if (bool(JA)){
235239
auto appJac = B;
236240
appJac.array() += time;
237-
JA = (B - dt*appJac);
241+
*JA.value() = (B - dt*appJac);
238242
}
243+
#else
244+
if (JA){
245+
auto appJac = B;
246+
appJac.array() += time;
247+
*JA = (B - dt*appJac);
248+
}
249+
#endif
239250
}
240251
};
241252
}

0 commit comments

Comments
 (0)