Hi all!
I'm trying to make use of Boost.Units in a simple FDTD simulation
program, but found a strange behavior of some kind. Consider the
following demonstration example where I try to calculate time step
from spatial steps according to the CFL condition [1]:
#include <iostream>
#include
#include
#include
#include
int main() {
namespace bu = boost::units;
namespace si = boost::units::si;
using dimless = bu::quantity;
bu::quantity
dx = 1.0f * si::meter,
dy = 1.0f * si::meter;
bu::quantity
v = 1.0f * si::meter_per_second;
bu::quantity
dt = 0.5f / (v * bu::sqrt(1.0f/(dx*dx) + 1.0f/(dy*dy)));
std::cout << dt << std::endl;
}
The code looks okay to me, but when I try to compile it, I get a scary
error message:
main.cpp:20:19: error: no match for ‘operator/’ (operand types
are ‘float’ and ‘’)
dt = 0.5f / (v * bu::sqrt(1.0f/(dx*dx) + 1.0f/(dy*dy)));
~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
where the second operand type looks like
boost::units::multiply_typeof_helper<
boost::units::quantity<
boost::units::unit<
[{length_base_dimension, static_rational<1l>},
{time_base_dimension, static_rational<-1l, 1l>}],
boost::units::si::system,
float>,
boost::units::quantity<
boost::units::unit<
[{length_base_dimension, static_rational<-1l, 1l>}],
boost::units::si::system,
double>>::type
In the above expression [{}, {}, ...] is a pseudo-code for static
head-tail list of boost::units::dim pairs. Note the
'double' type materialized from barely nowhere: the code explicitly
defines everything as 'float', it's the 'bu::sqrt' to blame for it to
appear. Note that replacing 'bu::sqrt' with 'bu::pow<1, 2>()' does
not help.
As I didn't feel happy with the above error message, I mutated the
code a bit and made it to work with:
bu::quantity dt;
dt = dimless(0.5f) / (v * bu::sqrt(1.0f/(dx*dx) + 1.0f/(dy*dy)));
I had to both (a) add 'dimless', and (b) separate variable declaration
and definition. Doing only (a) didn't help.
Mkay, but what if I need to have a const quantity? (I'm a big fan of
making everything 'const'!). The only one thing I could do was to roll
out my own 'sqrt' and 'pow' functions. Here they are:
#include
#include
#include
#include
#include
template
using static_power_helper = boost::units::unit<
typename boost::units::static_power<
typename boost::units::get_dimension<Unit>::type,
typename boost::units::static_rational::type>::type,
typename boost::units::get_system<Unit>::type>;
template
auto ratpow(boost::units::quantity const &x)
-> boost::units::quantity, Y> {
using std::pow;
return decltype(ratpow<N>(x))::from_value(pow(x.value(), N));
}
template
auto ratpow(boost::units::quantity const &x)
-> boost::units::quantity, Y> {
using std::pow;
return decltype(ratpow(x))::from_value(pow(x.value(), (Y) N / M));
}
template
auto msqrt(boost::units::quantity const &x)
-> decltype(ratpow<1, 2>(x)) {
return ratpow<1, 2>(x);
}
Using them I was able to compile the original form of the expression
(as seen in the very first listing):
bu::quantity const
dt = 0.5f / (v * msqrt(1.0f/(dx*dx) + 1.0f/(dy*dy)));
The questions:
1. Is it possible to work around the problem without custom 'sqrt'?
2. Is the described behavior intentional, -or-
Is there an error in Boost.Units' implementation of 'sqrt'?
3. If the 'boost::unit::sqrt' is fine, is there something wrong with
my implementation?
――― Pavel K.
[1]: https://en.wikipedia.org/wiki/Courant-Friedrichs-Lewy_condition
