Commit 4c1ddb3f authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Speed-up the face slope limiter in Gizmo.

parent 50678501
......@@ -65,7 +65,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
kernel_long_gravity.h vector.h cache.h runner_doiact.h runner_doiact_vec.h runner_doiact_grav.h runner_doiact_fft.h \
runner_doiact_nosort.h units.h intrinsics.h minmax.h kick.h timestep.h drift.h adiabatic_index.h io_properties.h \
dimension.h equation_of_state.h part_type.h periodic.h memswap.h dump.h logger.h \
dimension.h equation_of_state.h part_type.h periodic.h memswap.h dump.h logger.h sign.h \
gravity.h gravity_io.h gravity_cache.h \
gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
gravity/Default/gravity_debug.h gravity/Default/gravity_part.h \
......
......@@ -32,6 +32,10 @@
#ifndef SWIFT_GIZMO_SLOPE_LIMITER_FACE_H
#define SWIFT_GIZMO_SLOPE_LIMITER_FACE_H
#include <float.h>
#include "sign.h"
__attribute__((always_inline)) INLINE static float
hydro_slope_limit_face_quantity(float phi_i, float phi_j, float phi_mid0,
float xij_norm, float r_inv) {
......@@ -39,10 +43,6 @@ hydro_slope_limit_face_quantity(float phi_i, float phi_j, float phi_mid0,
const float psi1 = 0.5f;
const float psi2 = 0.25f;
if (phi_i == phi_j) {
return 0.0f;
}
const float delta1 = psi1 * fabsf(phi_i - phi_j);
const float delta2 = psi2 * fabsf(phi_i - phi_j);
......@@ -53,27 +53,15 @@ hydro_slope_limit_face_quantity(float phi_i, float phi_j, float phi_mid0,
float phiplus, phiminus, phi_mid;
/* if sign(phimax+delta1) == sign(phimax) */
if ((phimax + delta1) * phimax > 0.0f) {
if (same_signf(phimax + delta1, phimax))
phiplus = phimax + delta1;
} else {
if (phimax != 0.f) {
phiplus = phimax / (1.0f + delta1 / fabsf(phimax));
} else {
phiplus = 0.f;
}
}
else
phiplus = phimax / (1.0f + delta1 / (fabsf(phimax) + FLT_MIN));
/* if sign(phimin-delta1) == sign(phimin) */
if ((phimin - delta1) * phimin > 0.0f) {
if (same_signf(phimin - delta1, phimin))
phiminus = phimin - delta1;
} else {
if (phimin != 0.f) {
phiminus = phimin / (1.0f + delta1 / fabsf(phimin));
} else {
phiminus = 0.f;
}
}
else
phiminus = phimin / (1.0f + delta1 / (fabsf(phimin) + FLT_MIN));
if (phi_i < phi_j) {
const float temp = min(phibar + delta2, phi_mid0);
......
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2018 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_SIGN_H
#define SWIFT_SIGN_H
/**
* @brief Return the sign of a floating point number
*
* @param x The number of interest.
* @return >0 if positive, 0 if negative.
*/
__attribute__((always_inline)) INLINE static int signf(float x) {
#ifdef __GNUC__
return !signbit(x);
#else
return (0.f < val) - (val < 0.f);
#endif
}
/**
* @brief Return 1 if two numbers have the same sign, 0 otherwise
*
* @param x The first number
* @param y The second number
*/
__attribute__((always_inline)) INLINE static int same_signf(float x, float y) {
return signf(x) == signf(y);
}
#endif
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment