/* * Derivative.hpp * * Created on: Oct 5, 2015 * Author: Pietro Incardona */ #ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_DERIVATIVE_HPP_ #define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_DERIVATIVE_HPP_ #define CENTRAL 0 #define CENTRAL_B_ONE_SIDE 1 #define FORWARD 2 #define BACKWARD 3 #define CENTRAL_SYM 4 #include "util/mathutil.hpp" #include "Vector/map_vector.hpp" #include "Grid/comb.hpp" #include "FiniteDifference/util/common.hpp" #include "util/util_num.hpp" #include /*! \brief Derivative second order on h (spacing) * * \tparam d on which dimension derive * \tparam Field which field derive * \tparam impl which implementation * */ template class D { /*! \brief Calculate which colums of the Matrix are non zero * * \param pos position where the derivative is calculated * \param gs Grid info * \param cols non-zero colums calculated by the function * \param coeff coefficent (constant in front of the derivative) * * ### Example * * \snippet FDScheme_unit_tests.hpp Usage of stencil derivative * */ inline static void value(const grid_key_dx & pos, const grid_sm & gs, std::unordered_map & cols, typename Sys_eqs::stype coeff) { std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined"; } /*! \brief Calculate the position where the derivative is calculated * * In case of non staggered case this function just return a null grid_key, in case of staggered, * it calculate how the operator shift the calculation in the cell * * \param pos position where we are calculating the derivative * \param gs Grid info * \param s_pos staggered position of the properties * * \return where (in which cell) the derivative is calculated * */ inline static grid_key_dx position(grid_key_dx & pos, const grid_sm & gs, const comb (& s_pos)[Sys_eqs::nvar]) { std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined"; return pos; } }; /*! \brief Second order central Derivative scheme on direction i * * \verbatim * * -1 +1 * *---+---* * * \endverbatim * */ template class D { public: /*! \brief Calculate which colums of the Matrix are non zero * * \param g_map mapping grid * \param kmap position where the derivative is calculated * \param gs Grid info * \param spacing grid spacing * \param cols non-zero colums calculated by the function * \param coeff coefficent (constant in front of the derivative) * * ### Example * * \snippet FDScheme_unit_tests.hpp Usage of stencil derivative * */ inline static void value(const typename stub_or_real::type & g_map, grid_dist_key_dx & kmap, const grid_sm & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims], std::unordered_map & cols, typename Sys_eqs::stype coeff) { // if the system is staggered the CENTRAL derivative is equivalent to a forward derivative if (is_grid_staggered::value()) { D::value(g_map,kmap,gs,spacing,cols,coeff); return; } long int old_val = kmap.getKeyRef().get(d); kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) + 1); arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[d]/2.0 ); kmap.getKeyRef().set_d(d,old_val); old_val = kmap.getKeyRef().get(d); kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) - 1); arg::value(g_map,kmap,gs,spacing,cols,-coeff/spacing[d]/2.0 ); kmap.getKeyRef().set_d(d,old_val); } /*! \brief Calculate the position where the derivative is calculated * * In case on non staggered case this function just return a null grid_key, in case of staggered, * it calculate how the operator shift in the cell * \verbatim +--$--+ | | # * # | | 0--$--+ # = velocity(y) $= velocity(x) * = pressure \endverbatim * * Consider this 2D staggered cell and a second order central derivative scheme, this lead to * * \f$ \frac{\partial v_y}{\partial x} \f$is calculated on position (*), so the function return the grid_key {0,0} * * \f$ \frac{\partial v_y}{\partial y} \f$is calculated on position (0), so the function return the grid_key {-1,-1} * * \f$ \frac{\partial v_x}{\partial x} \f$is calculated on position (0), so the function return the grid_key {-1,-1} * * \f$ \frac{\partial v_x}{\partial y} \f$is calculated on position (*), so the function return the grid_key {0,0} * * \param pos position where we are calculating the derivative * \param gs Grid info * \param s_pos staggered position of the properties * * \return where (in which cell grid) the derivative is calculated * */ inline static grid_key_dx position(grid_key_dx & pos, const grid_sm & gs, const comb (& s_pos)[Sys_eqs::nvar]) { auto arg_pos = arg::position(pos,gs,s_pos); if (is_grid_staggered::value()) { if (arg_pos.get(d) == -1) { arg_pos.set_d(d,0); return arg_pos; } else { arg_pos.set_d(d,-1); return arg_pos; } } return arg_pos; } }; /*! \brief Second order one sided Derivative scheme on direction i * * \verbatim * * -1.5 2.0 -0.5 * +------*------* * * or * * -0.5 2.0 -1.5 * *------*------+ * * in the bulk * * -1 +1 * *---+---* * * \endverbatim * */ template class D { public: /*! \brief Calculate which colums of the Matrix are non zero * * \param g_map mapping grid points * \param kmap position where the derivative is calculated * \param gs Grid info * \param spacing of the grid * \param cols non-zero colums calculated by the function * \param coeff coefficent (constant in front of the derivative) * * ### Example * * \snippet FDScheme_unit_tests.hpp Usage of stencil derivative * */ static void value(const typename stub_or_real::type & g_map, grid_dist_key_dx & kmap, const grid_sm & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims], std::unordered_map & cols, typename Sys_eqs::stype coeff) { #ifdef SE_CLASS1 if (Sys_eqs::boundary[d] == PERIODIC) std::cerr << __FILE__ << ":" << __LINE__ << " error, it make no sense use one sided derivation with periodic boundary, please use CENTRAL\n"; #endif grid_key_dx pos = g_map.getGKey(kmap); if (pos.get(d) == (long int)gs.size(d)-1 ) { arg::value(g_map,kmap,gs,spacing,cols,1.5*coeff/spacing[d]); long int old_val = kmap.getKeyRef().get(d); kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) - 1); arg::value(g_map,kmap,gs,spacing,cols,-2.0*coeff/spacing[d]); kmap.getKeyRef().set_d(d,old_val); old_val = kmap.getKeyRef().get(d); kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) - 2); arg::value(g_map,kmap,gs,spacing,cols,0.5*coeff/spacing[d]); kmap.getKeyRef().set_d(d,old_val); } else if (pos.get(d) == 0) { arg::value(g_map,kmap,gs,spacing,cols,-1.5*coeff/spacing[d]); long int old_val = kmap.getKeyRef().get(d); kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) + 1); arg::value(g_map,kmap,gs,spacing,cols,2.0*coeff/spacing[d]); kmap.getKeyRef().set_d(d,old_val); old_val = kmap.getKeyRef().get(d); kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) + 2); arg::value(g_map,kmap,gs,spacing,cols,-0.5*coeff/spacing[d]); kmap.getKeyRef().set_d(d,old_val); } else { long int old_val = kmap.getKeyRef().get(d); kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) + 1); arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[d]); kmap.getKeyRef().set_d(d,old_val); old_val = kmap.getKeyRef().get(d); kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) - 1); arg::value(g_map,kmap,gs,spacing,cols,-coeff/spacing[d]); kmap.getKeyRef().set_d(d,old_val); } } /*! \brief Calculate the position where the derivative is calculated * * In case on non staggered case this function just return a null grid_key, in case of staggered, * it calculate how the operator shift in the cell * \verbatim +--$--+ | | # * # | | 0--$--+ # = velocity(y)$ = velocity(x) * = pressure \endverbatim * * In the one side stencil the cell position depend if you are or not at the boundary. * outside the boundary is simply the central scheme, at the boundary it is simply the * staggered position of the property * * \param pos position where we are calculating the derivative * \param gs Grid info * \param s_pos staggered position of the properties * * \return where (in which cell grid) the derivative is calculated * */ inline static grid_key_dx position(grid_key_dx & pos, const grid_sm & gs, const comb (& s_pos)[Sys_eqs::nvar]) { return arg::position(pos,gs,s_pos); } }; /*! \brief First order FORWARD derivative on direction i * * \verbatim * * -1.0 1.0 * +------* * * \endverbatim * */ template class D { public: /*! \brief Calculate which colums of the Matrix are non zero * * \param g_map mapping grid * \param kmap position where the derivative is calculated * \param gs Grid info * \param spacing grid spacing * \param cols non-zero colums calculated by the function * \param coeff coefficent (constant in front of the derivative) * * ### Example * * \snippet FDScheme_unit_tests.hpp Usage of stencil derivative * */ inline static void value(const typename stub_or_real::type & g_map, grid_dist_key_dx & kmap, const grid_sm & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims], std::unordered_map & cols, typename Sys_eqs::stype coeff) { long int old_val = kmap.getKeyRef().get(d); kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) + 1); arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[d]); kmap.getKeyRef().set_d(d,old_val); // backward arg::value(g_map,kmap,gs,spacing,cols,-coeff/spacing[d]); } /*! \brief Calculate the position where the derivative is calculated * * In case of non staggered case this function just return a null grid_key, in case of staggered, * the FORWARD scheme return the position of the staggered property * * \param pos position where we are calculating the derivative * \param gs Grid info * \param s_pos staggered position of the properties * * \return where (in which cell grid) the derivative is calculated * */ inline static grid_key_dx position(grid_key_dx & pos, const grid_sm & gs, const comb (& s_pos)[Sys_eqs::nvar]) { return arg::position(pos,gs,s_pos); } }; /*! \brief First order BACKWARD derivative on direction i * * \verbatim * * -1.0 1.0 * *------+ * * \endverbatim * */ template class D { public: /*! \brief Calculate which colums of the Matrix are non zero * * \param g_map mapping grid * \param kmap position where the derivative is calculated * \param gs Grid info * \param spacing of the grid * \param cols non-zero colums calculated by the function * \param coeff coefficent (constant in front of the derivative) * * ### Example * * \snippet FDScheme_unit_tests.hpp Usage of stencil derivative * */ inline static void value(const typename stub_or_real::type & g_map, grid_dist_key_dx & kmap, const grid_sm & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims], std::unordered_map & cols, typename Sys_eqs::stype coeff) { long int old_val = kmap.getKeyRef().get(d); kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) - 1); arg::value(g_map,kmap,gs,spacing,cols,-coeff/spacing[d]); kmap.getKeyRef().set_d(d,old_val); // forward arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[d]); } /*! \brief Calculate the position where the derivative is calculated * * In case of non staggered case this function just return a null grid_key, in case of staggered, * the BACKWARD scheme return the position of the staggered property * * \param pos position where we are calculating the derivative * \param gs Grid info * \param s_pos staggered position of the properties * * \return where (in which cell grid) the derivative is calculated * */ inline static grid_key_dx position(grid_key_dx & pos, const grid_sm & gs, const comb (& s_pos)[Sys_eqs::nvar]) { return arg::position(pos,gs,s_pos); } }; #endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_DERIVATIVE_HPP_ */