roystgnr / MetaPhysicL

Metaprogramming and operator-overloaded classes for numerical simulations
Other
22 stars 12 forks source link

Failing std::pow instantiation #37

Closed dschwen closed 5 years ago

dschwen commented 5 years ago

Instantiation of std::pow with a DualNumber base and an immediate double constant exponent fails with no matching function for call to 'if_else'

See long error below. (@lindsayad have you encountered this?

/Users/schwd/Programs/libmesh/installed/include/metaphysicl/dualnumber.h:415:795: error: no matching function for call to 'if_else'
template <typename T, typename D, typename T2, typename D2> inline typename CompareTypes<DualNumber<T,D>,DualNumber<T2,D2> >::supertype pow (const DualNumber<T,D>& a, const DualNumber<T2,D2>& b) { typedef typename CompareTypes<T,T2>::supertype TS; typedef typename CompareTypes<DualNumber<T,D>,DualNumber<T2,D2> >::supertype type; TS funcval = std::pow(a.value(), b.value()); return type(funcval, funcval * (b.value() * a.derivatives() / a.value() + MetaPhysicL::if_else(b.derivatives(), b.derivatives() * std::log(a.value()), b.derivatives()))); } template <typename T, typename D> inline DualNumber<T,D> pow (const DualNumber<T,D>& a, const DualNumber<T,D>& b) { T funcval = std::pow(a.value(), b.value()); return DualNumber<T,D>(funcval, funcval * (b.value() * a.derivatives() / a.value() + MetaPhysicL::if_else(b.derivatives(), b.derivatives() * std::log(a.value()), b.derivatives()))); } template <typename T, typename T2, typename D> inline typename CompareTypes<DualNumber<T2,D>,T,true>::supertype pow (const T& a, const DualNumber<T2,D>& b) { typedef typename CompareTypes<DualNumber<T2,D>,T,true>::supertype type; type newa(a); return std::pow(newa, b); } template <typename T, typename T2, typename D> inline typename CompareTypes<DualNumber<T,D>,T2>::supertype pow (const DualNumber<T,D>& a, const T2& b) { typedef typename CompareTypes<DualNumber<T,D>,T2>::supertype type; type newb(b); return std::pow(a, newb); }
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          ^~~~~~~~~~~~~~~~~~~~
/Users/schwd/Programs/libmesh/installed/include/metaphysicl/dualnumber.h:415:1412: note: in instantiation of function template specialization 'std::pow<double, MetaPhysicL::NumberArray<50, double> >' requested here
template <typename T, typename D, typename T2, typename D2> inline typename CompareTypes<DualNumber<T,D>,DualNumber<T2,D2> >::supertype pow (const DualNumber<T,D>& a, const DualNumber<T2,D2>& b) { typedef typename CompareTypes<T,T2>::supertype TS; typedef typename CompareTypes<DualNumber<T,D>,DualNumber<T2,D2> >::supertype type; TS funcval = std::pow(a.value(), b.value()); return type(funcval, funcval * (b.value() * a.derivatives() / a.value() + MetaPhysicL::if_else(b.derivatives(), b.derivatives() * std::log(a.value()), b.derivatives()))); } template <typename T, typename D> inline DualNumber<T,D> pow (const DualNumber<T,D>& a, const DualNumber<T,D>& b) { T funcval = std::pow(a.value(), b.value()); return DualNumber<T,D>(funcval, funcval * (b.value() * a.derivatives() / a.value() + MetaPhysicL::if_else(b.derivatives(), b.derivatives() * std::log(a.value()), b.derivatives()))); } template <typename T, typename T2, typename D> inline typename CompareTypes<DualNumber<T2,D>,T,true>::supertype pow (const T& a, const DualNumber<T2,D>& b) { typedef typename CompareTypes<DualNumber<T2,D>,T,true>::supertype type; type newa(a); return std::pow(newa, b); } template <typename T, typename T2, typename D> inline typename CompareTypes<DualNumber<T,D>,T2>::supertype pow (const DualNumber<T,D>& a, const T2& b) { typedef typename CompareTypes<DualNumber<T,D>,T2>::supertype type; type newb(b); return std::pow(a, newb); }
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   ^
/Users/schwd/Programs/moose/modules/tensor_mechanics/src/materials/ADComputeFiniteStrain.C:80:26: note: in instantiation of function template specialization 'std::pow<double, double, MetaPhysicL::NumberArray<50, double> >' requested here
      _Fhat[_qp] *= std::pow(ave_Fhat.det() / _Fhat[_qp].det(), 1.0 / 3.0);
dschwen commented 5 years ago

Simple code example (using the MOOSE typedefs)

  ADReal a;
  ADReal b;
  auto r1 = std::pow(a, b);

also fails with

ed/include/metaphysicl/dualnumber.h:415:795: error: no matching function for call to 'if_else'
template <typename T, typename D, typename T2, typename D2> inline typename CompareTypes<DualNumber<T,D>,DualNumber<T2,D2> >::supertype pow (const DualNumber<T,D>& a, const DualNumber<T2,D2>& b) { typedef typename CompareTypes<T,T2>::supertype TS; typedef typename CompareTypes<DualNumber<T,D>,DualNumber<T2,D2> >::supertype type; TS funcval = std::pow(a.value(), b.value()); return type(funcval, funcval * (b.value() * a.derivatives() / a.value() + MetaPhysicL::if_else(b.derivatives(), b.derivatives() * std::log(a.value()), b.derivatives()))); } template <typename T, typename D> inline DualNumber<T,D> pow (const DualNumber<T,D>& a, const DualNumber<T,D>& b) { T funcval = std::pow(a.value(), b.value()); return DualNumber<T,D>(funcval, funcval * (b.value() * a.derivatives() / a.value() + MetaPhysicL::if_else(b.derivatives(), b.derivatives() * std::log(a.value()), b.derivatives()))); } template <typename T, typename T2, typename D> inline typename CompareTypes<DualNumber<T2,D>,T,true>::supertype pow (const T& a, const DualNumber<T2,D>& b) { typedef typename CompareTypes<DualNumber<T2,D>,T,true>::supertype type; type newa(a); return std::pow(newa, b); } template <typename T, typename T2, typename D> inline typename CompareTypes<DualNumber<T,D>,T2>::supertype pow (const DualNumber<T,D>& a, const T2& b) { typedef typename CompareTypes<DualNumber<T,D>,T2>::supertype type; type newb(b); return std::pow(a, newb); }
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          ^~~~~~~~~~~~~~~~~~~~
/Users/schwd/Programs/moose/modules/tensor_mechanics/src/materials/ComputeFiniteStrain.C:44:18: note: in instantiation of function template specialization 'std::pow<double, MetaPhysicL::NumberArray<50, double> >' requested here
  auto r1 = std::pow(a, b);

What ever happened to SFINAE?!!!

dschwen commented 5 years ago

With

  ADReal a;
  ADReal b;
  auto r1 = std::pow<ADReal, ADReal>(a, b);

I get

/Users/schwd/Programs/libmesh/installed/include/metaphysicl/dualnumber.h:32:
/Users/schwd/Programs/libmesh/installed/include/metaphysicl/dualnumber_decl.h:143:40: error: no viable conversion from returned value of type 'const MetaPhysicL::NumberArray<50, double>' to function return type 'double'
  static T value(const T2& v) { return v; }
                                       ^
/Users/schwd/Programs/libmesh/installed/include/metaphysicl/dualnumber.h:98:37: note: in instantiation of function template specialization 'MetaPhysicL::DualNumberConstructor<double, MetaPhysicL::NumberArray<50, double> >::value<MetaPhysicL::NumberArray<50, double> >' requested here
  _val (DualNumberConstructor<T,D>::value(val)),

It looks like std::pow is just hosed (it is implemented in dualnumber.h though)

roystgnr commented 5 years ago

Huh. It looks like we don't actually have any unit tests taking pow(DualNumber,double), but we do have a swath of pow(NumberArray<DualNumber>, double, which should use the former under the hood. Can you add something in MetaPhysicL/test/ that replicates the problem there? I'll try to do the same.

lindsayad commented 5 years ago

This will currently fail if dualnumber is included before numberarray...

dschwen commented 5 years ago

Found two cases of "wrong" include order in MOOSE and swapping it fixes it.

lindsayad commented 5 years ago

@roystgnr Perhaps we should add numberarray_decl and numbervector_decl and put these at the top of dualnumber. Or I guess manually put some of these declarations at the top...

roystgnr commented 5 years ago

"Make every header depend on every single decl header" is what that would boil down to, one added include at a time, which isn't a solution that scales very well.

Perhaps we could avoid this bug in the future by removing an include? If we no longer make foo.h include foo_decl.h, that would force users to directly include foo_decl.h rather than rely on indirect inclusion, which would make it harder to accidentally leave a decl header until too late.

I feel like the right solution would be some sort of static_assert that detects this sort of failure and gives a more sane error message, but I can't figure out how to do that.

roystgnr commented 5 years ago

Was this resolved by #38?

lindsayad commented 5 years ago

Yep!