Skip to content
Snippets Groups Projects
Commit c9f3ad26 authored by 数学の武士's avatar 数学の武士
Browse files

float built in functions instead of double

parent eb5c7056
Branches
Tags
No related merge requests found
......@@ -9,22 +9,22 @@ FUCNTION_DECLARATION_SPECIFIER void get_convection_lim(T &zeta_lim, T &Rib_lim,
{
T psi_m, psi_h, f_m, f_h, c;
c = pow(param.Pr_t_inf_inv / param.Pr_t_0_inv, 4);
c = powf(param.Pr_t_inf_inv / param.Pr_t_0_inv, 4);
zeta_lim = (2.0 * param.alpha_h - c * param.alpha_m -
sqrt( (c * param.alpha_m)*(c * param.alpha_m) + 4.0 * c * param.alpha_h * (param.alpha_h - param.alpha_m))) / (2.0 * param.alpha_h*param.alpha_h);
sqrtf( (c * param.alpha_m)*(c * param.alpha_m) + 4.0 * c * param.alpha_h * (param.alpha_h - param.alpha_m))) / (2.0 * param.alpha_h*param.alpha_h);
f_m_lim = pow(1.0 - param.alpha_m * zeta_lim, 0.25);
f_h_lim = sqrt(1.0 - param.alpha_h * zeta_lim);
f_m_lim = powf(1.0 - param.alpha_m * zeta_lim, 0.25);
f_h_lim = sqrtf(1.0 - param.alpha_h * zeta_lim);
f_m = zeta_lim / h0_m;
f_h = zeta_lim / h0_t;
if (fabs(B) < 1.0e-10) f_h = f_m;
f_m = pow(1.0 - param.alpha_m * f_m, 0.25);
f_h = sqrt(1.0 - param.alpha_h_fix * f_h);
f_m = powf(1.0 - param.alpha_m * f_m, 0.25);
f_h = sqrtf(1.0 - param.alpha_h_fix * f_h);
psi_m = 2.0 * (atan(f_m_lim) - atan(f_m)) + log((f_m_lim - 1.0) * (f_m + 1.0)/((f_m_lim + 1.0) * (f_m - 1.0)));
psi_h = log((f_h_lim - 1.0) * (f_h + 1.0)/((f_h_lim + 1.0) * (f_h - 1.0))) / param.Pr_t_0_inv;
psi_m = 2.0 * (atanf(f_m_lim) - atanf(f_m)) + logf((f_m_lim - 1.0) * (f_m + 1.0)/((f_m_lim + 1.0) * (f_m - 1.0)));
psi_h = logf((f_h_lim - 1.0) * (f_h + 1.0)/((f_h_lim + 1.0) * (f_h - 1.0))) / param.Pr_t_0_inv;
Rib_lim = zeta_lim * psi_h / (psi_m * psi_m);
}
......@@ -36,12 +36,12 @@ FUCNTION_DECLARATION_SPECIFIER void get_psi_stable(T &psi_m, T &psi_h, T &zeta,
{
T Rib_coeff, psi0_m, psi0_h, phi, c;
psi0_m = log(h0_m);
psi0_m = logf(h0_m);
psi0_h = B / psi0_m;
Rib_coeff = param.beta_m * Rib;
c = (psi0_h + 1.0) / param.Pr_t_0_inv - 2.0 * Rib_coeff;
zeta = psi0_m * (sqrt(c*c + 4.0 * Rib_coeff * (1.0 - Rib_coeff)) - c) / (2.0 * param.beta_m * (1.0 - Rib_coeff));
zeta = psi0_m * (sqrtf(c*c + 4.0 * Rib_coeff * (1.0 - Rib_coeff)) - c) / (2.0 * param.beta_m * (1.0 - Rib_coeff));
phi = param.beta_m * zeta;
......@@ -58,8 +58,8 @@ FUCNTION_DECLARATION_SPECIFIER void get_psi_convection(T &psi_m, T &psi_h, T &ze
{
T zeta0_m, zeta0_h, f0_m, f0_h, p_m, p_h, a_m, a_h, c_lim, f;
p_m = 2.0 * atan(f_m_conv_lim) + log((f_m_conv_lim - 1.0) / (f_m_conv_lim + 1.0));
p_h = log((f_h_conv_lim - 1.0) / (f_h_conv_lim + 1.0));
p_m = 2.0 * atanf(f_m_conv_lim) + logf((f_m_conv_lim - 1.0) / (f_m_conv_lim + 1.0));
p_h = logf((f_h_conv_lim - 1.0) / (f_h_conv_lim + 1.0));
zeta = zeta_conv_lim;
......@@ -70,13 +70,13 @@ FUCNTION_DECLARATION_SPECIFIER void get_psi_convection(T &psi_m, T &psi_h, T &ze
if (fabs(B) < 1.0e-10)
zeta0_h = zeta0_m;
f0_m = pow(1.0 - param.alpha_m * zeta0_m, 0.25);
f0_h = sqrt(1.0 - param.alpha_h_fix * zeta0_h);
f0_m = powf(1.0 - param.alpha_m * zeta0_m, 0.25);
f0_h = sqrtf(1.0 - param.alpha_h_fix * zeta0_h);
a_m = -2.0*atan(f0_m) + log((f0_m + 1.0)/(f0_m - 1.0));
a_h = log((f0_h + 1.0)/(f0_h - 1.0));
a_m = -2.0*atanf(f0_m) + logf((f0_m + 1.0)/(f0_m - 1.0));
a_h = logf((f0_h + 1.0)/(f0_h - 1.0));
c_lim = pow(zeta_conv_lim / zeta, 1.0 / 3.0);
c_lim = powf(zeta_conv_lim / zeta, 1.0 / 3.0);
f = 3.0 * (1.0 - c_lim);
psi_m = f / f_m_conv_lim + p_m + a_m;
......@@ -95,8 +95,8 @@ FUCNTION_DECLARATION_SPECIFIER void get_psi_neutral(T &psi_m, T &psi_h, T &zeta,
const sfx_param& param)
{
zeta = 0.0;
psi_m = log(h0_m);
psi_h = log(h0_t) / param.Pr_t_0_inv;
psi_m = logf(h0_m);
psi_h = logf(h0_t) / param.Pr_t_0_inv;
if (fabs(B) < 1.0e-10)
psi_h = psi_m / param.Pr_t_0_inv;
}
......@@ -109,8 +109,8 @@ FUCNTION_DECLARATION_SPECIFIER void get_psi_semi_convection(T &psi_m, T &psi_h,
{
T zeta0_m, zeta0_h, f0_m, f0_h, f_m, f_h;
psi_m = log(h0_m);
psi_h = log(h0_t);
psi_m = logf(h0_m);
psi_h = logf(h0_t);
if (fabs(B) < 1.0e-10)
psi_h = psi_m;
......@@ -124,17 +124,17 @@ FUCNTION_DECLARATION_SPECIFIER void get_psi_semi_convection(T &psi_m, T &psi_h,
if (fabs(B) < 1.0e-10)
zeta0_h = zeta0_m;
f_m = pow(1.0 - param.alpha_m * zeta, 0.25e0);
f_h = sqrt(1.0 - param.alpha_h_fix * zeta);
f_m = powf(1.0 - param.alpha_m * zeta, 0.25e0);
f_h = sqrtf(1.0 - param.alpha_h_fix * zeta);
f0_m = pow(1.0 - param.alpha_m * zeta0_m, 0.25e0);
f0_h = sqrt(1.0 - param.alpha_h_fix * zeta0_h);
f0_m = powf(1.0 - param.alpha_m * zeta0_m, 0.25e0);
f0_h = sqrtf(1.0 - param.alpha_h_fix * zeta0_h);
f0_m = sfx_math::max(f0_m, T(1.000001e0));
f0_h = sfx_math::max(f0_h, T(1.000001e0));
psi_m = log((f_m - 1.0e0)*(f0_m + 1.0e0)/((f_m + 1.0e0)*(f0_m - 1.0e0))) + 2.0e0*(atan(f_m) - atan(f0_m));
psi_h = log((f_h - 1.0e0)*(f0_h + 1.0e0)/((f_h + 1.0e0)*(f0_h - 1.0e0))) / param.Pr_t_0_inv;
psi_m = logf((f_m - 1.0e0)*(f0_m + 1.0e0)/((f_m + 1.0e0)*(f0_m - 1.0e0))) + 2.0e0*(atanf(f_m) - atanf(f0_m));
psi_h = logf((f_h - 1.0e0)*(f0_h + 1.0e0)/((f_h + 1.0e0)*(f0_h - 1.0e0))) / param.Pr_t_0_inv;
if (i == maxiters_convection + 1)
break;
......@@ -153,23 +153,23 @@ FUCNTION_DECLARATION_SPECIFIER void get_psi(T &psi_m, T &psi_h,
if (zeta >= 0.0)
{
q_m = pow((1.0 - param.b_m) / param.b_m, 1.0 / 3.0);
q_h = sqrt(param.c_h * param.c_h - 4.0);
q_m = powf((1.0 - param.b_m) / param.b_m, 1.0 / 3.0);
q_h = sqrtf(param.c_h * param.c_h - 4.0);
x_m = pow(1.0 + zeta, 1.0 / 3.0);
x_m = powf(1.0 + zeta, 1.0 / 3.0);
x_h = zeta;
psi_m = -3.0 * (param.a_m / param.b_m) * (x_m - 1.0) + 0.5 * (param.a_m / param.b_m) * q_m * (2.0 * log((x_m + q_m) / (1.0 + q_m)) - log((x_m * x_m - x_m * q_m + q_m * q_m) / (1.0 - q_m + q_m * q_m)) + 2.0 * sqrt(3.0) * (atan((2.0 * x_m - q_m) / (sqrt(3.0) * q_m)) - atan((2.0 - q_m) / (sqrt(3.0) * q_m))));
psi_m = -3.0 * (param.a_m / param.b_m) * (x_m - 1.0) + 0.5 * (param.a_m / param.b_m) * q_m * (2.0 * logf((x_m + q_m) / (1.0 + q_m)) - logf((x_m * x_m - x_m * q_m + q_m * q_m) / (1.0 - q_m + q_m * q_m)) + 2.0 * sqrtf(3.0) * (atanf((2.0 * x_m - q_m) / (sqrtf(3.0) * q_m)) - atanf((2.0 - q_m) / (sqrtf(3.0) * q_m))));
psi_h = -0.5 * param.b_h * log(1.0 + param.c_h * x_h + x_h * x_h) + ((-param.a_h / q_h) + ((param.b_h * param.c_h) / (2.0 * q_h))) * (log((2.0 * x_h + param.c_h - q_h) / (2.0 * x_h + param.c_h + q_h)) - log((param.c_h - q_h) / (param.c_h + q_h)));
psi_h = -0.5 * param.b_h * logf(1.0 + param.c_h * x_h + x_h * x_h) + ((-param.a_h / q_h) + ((param.b_h * param.c_h) / (2.0 * q_h))) * (logf((2.0 * x_h + param.c_h - q_h) / (2.0 * x_h + param.c_h + q_h)) - logf((param.c_h - q_h) / (param.c_h + q_h)));
}
else
{
x_m = pow(1.0 - param.alpha_m * zeta, 0.25);
x_h = pow(1.0 - param.alpha_h * zeta, 0.25);
x_m = powf(1.0 - param.alpha_m * zeta, 0.25);
x_h = powf(1.0 - param.alpha_h * zeta, 0.25);
psi_m = (4.0 * atan(1.0) / 2.0) + 2.0 * log(0.5 * (1.0 + x_m)) + log(0.5 * (1.0 + x_m * x_m)) - 2.0 * atan(x_m);
psi_h = 2.0 * log(0.5 * (1.0 + x_h * x_h));
psi_m = (4.0 * atanf(1.0) / 2.0) + 2.0 * logf(0.5 * (1.0 + x_m)) + logf(0.5 * (1.0 + x_m * x_m)) - 2.0 * atanf(x_m);
psi_h = 2.0 * logf(0.5 * (1.0 + x_h * x_h));
}
}
......@@ -183,28 +183,28 @@ FUCNTION_DECLARATION_SPECIFIER void get_psi_mh(T &psi_m, T &psi_h,
if (zeta_m >= 0.0)
{
q_m = pow((1.0 - param.b_m) / param.b_m, 1.0 / 3.0);
x_m = pow(1.0 + zeta_m, 1.0 / 3.0);
q_m = powf((1.0 - param.b_m) / param.b_m, 1.0 / 3.0);
x_m = powf(1.0 + zeta_m, 1.0 / 3.0);
psi_m = -3.0 * (param.a_m / param.b_m) * (x_m - 1.0) + 0.5 * (param.a_m / param.b_m) * q_m * (2.0 * log((x_m + q_m) / (1.0 + q_m)) - log((x_m * x_m - x_m * q_m + q_m * q_m) / (1.0 - q_m + q_m * q_m)) + 2.0 * sqrt(3.0) * (atan((2.0 * x_m - q_m) / (sqrt(3.0) * q_m)) - atan((2.0 - q_m) / (sqrt(3.0) * q_m))));
psi_m = -3.0 * (param.a_m / param.b_m) * (x_m - 1.0) + 0.5 * (param.a_m / param.b_m) * q_m * (2.0 * logf((x_m + q_m) / (1.0 + q_m)) - logf((x_m * x_m - x_m * q_m + q_m * q_m) / (1.0 - q_m + q_m * q_m)) + 2.0 * sqrtf(3.0) * (atanf((2.0 * x_m - q_m) / (sqrtf(3.0) * q_m)) - atanf((2.0 - q_m) / (sqrtf(3.0) * q_m))));
}
else
{
x_m = pow(1.0 - param.alpha_m * zeta_m, 0.25);
psi_m = (4.0 * atan(1.0) / 2.0) + 2.0 * log(0.5 * (1.0 + x_m)) + log(0.5 * (1.0 + x_m * x_m)) - 2.0 * atan(x_m);
x_m = powf(1.0 - param.alpha_m * zeta_m, 0.25);
psi_m = (4.0 * atanf(1.0) / 2.0) + 2.0 * logf(0.5 * (1.0 + x_m)) + logf(0.5 * (1.0 + x_m * x_m)) - 2.0 * atanf(x_m);
}
if (zeta_h >= 0.0)
{
q_h = sqrt(param.c_h * param.c_h - 4.0);
q_h = sqrtf(param.c_h * param.c_h - 4.0);
x_h = zeta_h;
psi_h = -0.5 * param.b_h * log(1.0 + param.c_h * x_h + x_h * x_h) + ((-param.a_h / q_h) + ((param.b_h * param.c_h) / (2.0 * q_h))) * (log((2.0 * x_h + param.c_h - q_h) / (2.0 * x_h + param.c_h + q_h)) - log((param.c_h - q_h) / (param.c_h + q_h)));
psi_h = -0.5 * param.b_h * logf(1.0 + param.c_h * x_h + x_h * x_h) + ((-param.a_h / q_h) + ((param.b_h * param.c_h) / (2.0 * q_h))) * (logf((2.0 * x_h + param.c_h - q_h) / (2.0 * x_h + param.c_h + q_h)) - logf((param.c_h - q_h) / (param.c_h + q_h)));
}
else
{
x_h = pow(1.0 - param.alpha_h * zeta_h, 0.25);
psi_h = 2.0 * log(0.5 * (1.0 + x_h * x_h));
x_h = powf(1.0 - param.alpha_h * zeta_h, 0.25);
psi_h = 2.0 * logf(0.5 * (1.0 + x_h * x_h));
}
}
......@@ -220,9 +220,9 @@ FUCNTION_DECLARATION_SPECIFIER void get_dynamic_scales(T &Udyn, T &Tdyn, T &Qdyn
const T gamma = T(0.61);
T psi_m, psi_h, psi0_m, psi0_h, Linv;
Udyn = param.kappa * U / log(z / z0_m);
Tdyn = param.kappa * dT * param.Pr_t_0_inv / log(z / z0_t);
Qdyn = param.kappa * dQ * param.Pr_t_0_inv / log(z / z0_t);
Udyn = param.kappa * U / logf(z / z0_m);
Tdyn = param.kappa * dT * param.Pr_t_0_inv / logf(z / z0_t);
Qdyn = param.kappa * dQ * param.Pr_t_0_inv / logf(z / z0_t);
zeta = 0.0;
// --- no wind
......@@ -242,9 +242,9 @@ FUCNTION_DECLARATION_SPECIFIER void get_dynamic_scales(T &Udyn, T &Tdyn, T &Qdyn
get_psi_mh(psi0_m, psi0_h, z0_m * Linv, z0_t * Linv,
param);
Udyn = param.kappa * U / (log(z / z0_m) - (psi_m - psi0_m));
Tdyn = param.kappa * dT * param.Pr_t_0_inv / (log(z / z0_t) - (psi_h - psi0_h));
Qdyn = param.kappa * dQ * param.Pr_t_0_inv / (log(z / z0_t) - (psi_h - psi0_h));
Udyn = param.kappa * U / (logf(z / z0_m) - (psi_m - psi0_m));
Tdyn = param.kappa * dT * param.Pr_t_0_inv / (logf(z / z0_t) - (psi_h - psi0_h));
Qdyn = param.kappa * dQ * param.Pr_t_0_inv / (logf(z / z0_t) - (psi_h - psi0_h));
if (Udyn < 1e-5)
break;
......@@ -261,12 +261,12 @@ FUCNTION_DECLARATION_SPECIFIER void get_phi(T &phi_m, T &phi_h,
{
if (zeta >= 0.0)
{
phi_m = 1.0 + (param.a_m * zeta * pow(1.0 + zeta, 1.0 / 3.0) ) / (1.0 + param.b_m * zeta);
phi_m = 1.0 + (param.a_m * zeta * powf(1.0 + zeta, 1.0 / 3.0) ) / (1.0 + param.b_m * zeta);
phi_h = 1.0 + (param.a_h * zeta + param.b_h * zeta * zeta) / (1.0 + param.c_h * zeta + zeta * zeta);
}
else
{
phi_m = pow(1.0 - param.alpha_m * zeta, -0.25);
phi_h = pow(1.0 - param.alpha_h * zeta, -0.5);
phi_m = powf(1.0 - param.alpha_m * zeta, -0.25);
phi_h = powf(1.0 - param.alpha_h * zeta, -0.5);
}
}
......@@ -8,9 +8,9 @@ FUCNTION_DECLARATION_SPECIFIER void get_thermal_roughness(T &z0_t, T &B,
const struct sfx_surface_param& param,
const int surface_type)
{
// --- define B = log(z0_m / z0_t)
B = (Re <= param.Re_rough_min) ? param.B1_rough * log(param.B3_rough * Re) + param.B2_rough :
param.B4_rough * (pow(Re, param.B2_rough));
// --- define B = logf(z0_m / z0_t)
B = (Re <= param.Re_rough_min) ? param.B1_rough * logf(param.B3_rough * Re) + param.B2_rough :
param.B4_rough * (powf(Re, param.B2_rough));
// --- apply max restriction based on surface type
if (surface_type == param.surface_ocean) B = sfx_math::min(B, param.B_max_ocean);
......@@ -32,22 +32,22 @@ FUCNTION_DECLARATION_SPECIFIER void get_charnock_roughness(T &z0_m, T &u_dyn0,
Uc = U;
a = 0.0;
b = 25.0;
c_min = log(surface_param.h_charnock) / surface_param.kappa;
c_min = logf(surface_param.h_charnock) / surface_param.kappa;
for (int i = 0; i < maxiters; i++)
{
f = surface_param.c1_charnock - 2.0 * log(Uc);
f = surface_param.c1_charnock - 2.0 * logf(Uc);
for (int j = 0; j < maxiters; j++)
{
c = (f + 2.0 * log(b)) / surface_param.kappa;
c = (f + 2.0 * logf(b)) / surface_param.kappa;
if (U <= 8.0e0)
a = log(1.0 + surface_param.c2_charnock * ( pow(b / Uc, 3) ) ) / surface_param.kappa;
a = logf(1.0 + surface_param.c2_charnock * ( powf(b / Uc, 3) ) ) / surface_param.kappa;
c = sfx_math::max(c - a, c_min);
b = c;
}
z0_m = surface_param.h_charnock * exp(-c * surface_param.kappa);
z0_m = sfx_math::max(z0_m, T(0.000015e0));
Uc = U * log(surface_param.h_charnock / z0_m) / log(h / z0_m);
Uc = U * logf(surface_param.h_charnock / z0_m) / logf(h / z0_m);
}
u_dyn0 = Uc / c;
......
......@@ -52,7 +52,7 @@ __global__ void sfx_kernel::compute_flux(sfxDataVecTypeC sfx,
if (surface_type == surface.surface_land)
{
h0_m = h / z0_m;
u_dyn0 = U * model.kappa / log(h0_m);
u_dyn0 = U * model.kappa / logf(h0_m);
}
Re = u_dyn0 * z0_m / phys.nu_air;
......@@ -79,7 +79,7 @@ __global__ void sfx_kernel::compute_flux(sfxDataVecTypeC sfx,
{
get_psi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, zeta_conv_lim, f_m_conv_lim, f_h_conv_lim, model, numerics.maxiters_convection);
fval = pow(zeta_conv_lim / zeta, 1.0/3.0);
fval = powf(zeta_conv_lim / zeta, 1.0/3.0);
phi_m = fval / f_m_conv_lim;
phi_h = fval / (model.Pr_t_0_inv * f_h_conv_lim);
}
......@@ -94,8 +94,8 @@ __global__ void sfx_kernel::compute_flux(sfxDataVecTypeC sfx,
{
get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, model, numerics.maxiters_convection);
phi_m = pow(1.0 - model.alpha_m * zeta, -0.25);
phi_h = 1.0 / (model.Pr_t_0_inv * sqrt(1.0 - model.alpha_h_fix * zeta));
phi_m = powf(1.0 - model.alpha_m * zeta, -0.25);
phi_h = 1.0 / (model.Pr_t_0_inv * sqrtf(1.0 - model.alpha_h_fix * zeta));
}
Cm = model.kappa / psi_m;
......
......@@ -53,7 +53,7 @@ __global__ void sfx_kernel::compute_flux(sfxDataVecTypeC sfx,
if (surface_type == surface.surface_land)
{
h0_m = h / z0_m;
u_dyn0 = U * model.kappa / log(h0_m);
u_dyn0 = U * model.kappa / logf(h0_m);
}
Re = u_dyn0 * z0_m / phys.nu_air;
......
......@@ -55,7 +55,7 @@ void FluxEsm<T, memIn, memOut, MemType::CPU>::compute_flux()
if (surface_type == surface.surface_land)
{
h0_m = h / z0_m;
u_dyn0 = U * model.kappa / log(h0_m);
u_dyn0 = U * model.kappa / logf(h0_m);
}
Re = u_dyn0 * z0_m / phys.nu_air;
......@@ -82,7 +82,7 @@ void FluxEsm<T, memIn, memOut, MemType::CPU>::compute_flux()
{
get_psi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, zeta_conv_lim, f_m_conv_lim, f_h_conv_lim, model, numerics.maxiters_convection);
fval = pow(zeta_conv_lim / zeta, 1.0/3.0);
fval = powf(zeta_conv_lim / zeta, 1.0/3.0);
phi_m = fval / f_m_conv_lim;
phi_h = fval / (model.Pr_t_0_inv * f_h_conv_lim);
}
......@@ -97,8 +97,8 @@ void FluxEsm<T, memIn, memOut, MemType::CPU>::compute_flux()
{
get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, model, numerics.maxiters_convection);
phi_m = pow(1.0 - model.alpha_m * zeta, -0.25);
phi_h = 1.0 / (model.Pr_t_0_inv * sqrt(1.0 - model.alpha_h_fix * zeta));
phi_m = powf(1.0 - model.alpha_m * zeta, -0.25);
phi_h = 1.0 / (model.Pr_t_0_inv * sqrtf(1.0 - model.alpha_h_fix * zeta));
}
Cm = model.kappa / psi_m;
......
......@@ -56,7 +56,7 @@ void FluxSheba<T, memIn, memOut, MemType::CPU>::compute_flux()
if (surface_type == surface.surface_land)
{
h0_m = h / z0_m;
u_dyn0 = U * model.kappa / log(h0_m);
u_dyn0 = U * model.kappa / logf(h0_m);
}
Re = u_dyn0 * z0_m / phys.nu_air;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment