stencilops.h

数值计算工具库,C语言编写的,可以直接调用.

template<class T>
inline TinyVector<_bz_typename T::T_numtype,3>
curln(T& vx, T& vy, T& vz) {
  const int x = firstDim, y = secondDim, z = thirdDim;

  return TinyVector<_bz_typename T::T_numtype,3>(
    (central12(vz,y)-central12(vy,z)) * recip_2,
    (central12(vx,z)-central12(vz,x)) * recip_2,
    (central12(vy,x)-central12(vx,y)) * recip_2);
}

// Multicomponent curl
template<class T>
inline _bz_typename T::T_numtype curl(T& A) {
  const int x = firstDim, y = secondDim, z = thirdDim;

  return T::T_numtype(
    central12(A,z,y)-central12(A,y,z),
    central12(A,x,z)-central12(A,z,x),
    central12(A,y,x)-central12(A,x,y));
}

// Normalized multicomponent curl
template<class T>
inline _bz_typename T::T_numtype curln(T& A) {
  const int x = firstDim, y = secondDim, z = thirdDim;

  return T::T_numtype(
    (central12(A,z,y)-central12(A,y,z)) * recip_2,
    (central12(A,x,z)-central12(A,z,x)) * recip_2,
    (central12(A,y,x)-central12(A,x,y)) * recip_2);
}

// O(h^4) curl, using 4th order central difference
template<class T>
inline TinyVector<_bz_typename T::T_numtype,3>
curl4(T& vx, T& vy, T& vz) {
  const int x = firstDim, y = secondDim, z = thirdDim;

  return TinyVector<_bz_typename T::T_numtype,3>(
    central14(vz,y)-central14(vy,z),
    central14(vx,z)-central14(vz,x),
    central14(vy,x)-central14(vx,y));
}

// O(h^4) curl, using 4th order central difference (multicomponent version)
template<class T>
inline _bz_typename T::T_numtype
curl4(T& A) {
  const int x = firstDim, y = secondDim, z = thirdDim;

  return T::T_numtype(
    central14(A,z,y)-central14(A,y,z),
    central14(A,x,z)-central14(A,z,x),
    central14(A,y,x)-central14(A,x,y));
}

// Normalized O(h^4) curl, using 4th order central difference
template<class T>
inline TinyVector<_bz_typename T::T_numtype,3>
curl4n(T& vx, T& vy, T& vz) {
  const int x = firstDim, y = secondDim, z = thirdDim;

  return TinyVector<_bz_typename T::T_numtype,3>(
    (central14(vz,y)-central14(vy,z)) * recip_2,
    (central14(vx,z)-central14(vz,x)) * recip_2,
    (central14(vy,x)-central14(vx,y)) * recip_2);
}

// O(h^4) curl, using 4th order central difference (normalized multicomponent)
template<class T>
inline _bz_typename T::T_numtype
curl4n(T& A) {
  const int x = firstDim, y = secondDim, z = thirdDim;

  return T::T_numtype(
    (central14(A,z,y)-central14(A,y,z)) * recip_2,
    (central14(A,x,z)-central14(A,z,x)) * recip_2,
    (central14(A,y,x)-central14(A,x,y)) * recip_2);
}



// Two-dimensional curl

template<class T>
inline _bz_typename T::T_numtype
curl(T& vx, T& vy) {
  const int x = firstDim, y = secondDim;

  return central12(vy,x)-central12(vx,y);
}

template<class T>
inline _bz_typename T::T_numtype
curln(T& vx, T& vy) {
  const int x = firstDim, y = secondDim;

  return (central12(vy,x)-central12(vx,y)) * recip_2;
}

// Multicomponent curl
template<class T>
inline _bz_typename T::T_numtype::T_numtype curl2D(T& A) {
  const int x = firstDim, y = secondDim;
  return central12(A,y,x)-central12(A,x,y);
}

template<class T>
inline _bz_typename T::T_numtype::T_numtype curl2Dn(T& A) {
  const int x = firstDim, y = secondDim;
  return (central12(A,y,x)-central12(A,x,y)) * recip_2;
}


// 4th order versions

template<class T>
inline _bz_typename T::T_numtype
curl4(T& vx, T& vy) {
  const int x = firstDim, y = secondDim;

  return central14(vy,x)-central14(vx,y);
}

template<class T>
inline _bz_typename T::T_numtype
curl4n(T& vx, T& vy) {
  const int x = firstDim, y = secondDim;

  return (central14(vy,x)-central14(vx,y)) * recip_12;
}

// Multicomponent curl
template<class T>
inline _bz_typename T::T_numtype::T_numtype curl2D4(T& A) {
  const int x = firstDim, y = secondDim;
  return central14(A,y,x)-central14(A,x,y);
}

template<class T>
inline _bz_typename T::T_numtype::T_numtype curl2D4n(T& A) {
  const int x = firstDim, y = secondDim;
  return (central14(A,y,x)-central14(A,x,y)) * recip_12;
}

/****************************************************************************
 * Divergence
 ****************************************************************************/


BZ_DECLARE_STENCIL_OPERATOR2(div,vx,vy)
  return central12(vx,firstDim) + central12(vy,secondDim);
BZ_END_STENCIL_OPERATOR

BZ_DECLARE_STENCIL_OPERATOR2(divn,vx,vy)
  return (central12(vx,firstDim) + central12(vy,secondDim))
     * recip_2;
BZ_END_STENCIL_OPERATOR

BZ_DECLARE_STENCIL_OPERATOR2(div4,vx,vy)
  return central14(vx,firstDim) + central14(vy,secondDim);
BZ_END_STENCIL_OPERATOR

BZ_DECLARE_STENCIL_OPERATOR2(div4n,vx,vy)
  return (central14(vx,firstDim) + central14(vy,secondDim))
    * recip_12;
BZ_END_STENCIL_OPERATOR

BZ_DECLARE_STENCIL_OPERATOR3(div,vx,vy,vz)
  return central12(vx,firstDim) + central12(vy,secondDim) 
    + central12(vz,thirdDim);
BZ_END_STENCIL_OPERATOR

BZ_DECLARE_STENCIL_OPERATOR3(divn,vx,vy,vz)
  return (central12(vx,firstDim) + central12(vy,secondDim) 
    + central12(vz,thirdDim)) * recip_2;
BZ_END_STENCIL_OPERATOR

BZ_DECLARE_STENCIL_OPERATOR3(div4,vx,vy,vz)
  return central14(vx,firstDim) + central14(vy,secondDim) 
    + central14(vz,thirdDim);
BZ_END_STENCIL_OPERATOR

BZ_DECLARE_STENCIL_OPERATOR3(div4n,vx,vy,vz)
  return (central14(vx,firstDim) + central14(vy,secondDim)
    + central14(vz,thirdDim)) * recip_12;
BZ_END_STENCIL_OPERATOR

template<class T>
inline _bz_typename T::T_numtype::T_numtype
div2D(T& A) {
  const int x = firstDim, y = secondDim;
  return central12(A,x,x) + central12(A,y,y);
}

template<class T>
inline _bz_typename T::T_numtype::T_numtype
div2D4(T& A) {
  const int x = firstDim, y = secondDim;
  return central14(A,x,x) + central14(A,y,y);
}

template<class T>
inline _bz_typename T::T_numtype::T_numtype
div2Dn(T& A) {
  const int x = firstDim, y = secondDim;
  return (central12(A,x,x) + central12(A,y,y)) * recip_2;
}

template<class T>
inline _bz_typename T::T_numtype::T_numtype
div2D4n(T& A) {
  const int x = firstDim, y = secondDim;
  return (central14(A,x,x) + central14(A,y,y)) * recip_12;
}

template<class T>
inline _bz_typename T::T_numtype::T_numtype
div3D(T& A) {
  const int x = firstDim, y = secondDim, z = thirdDim;
  return central12(A,x,x) + central12(A,y,y) + central12(A,z,z);
}

template<class T>
inline _bz_typename T::T_numtype::T_numtype
div3D4(T& A) {
  const int x = firstDim, y = secondDim, z = thirdDim;
  return central14(A,x,x) + central14(A,y,y) + central14(A,z,z);
}

template<class T>
inline _bz_typename T::T_numtype::T_numtype
div3Dn(T& A) {
  const int x = firstDim, y = secondDim, z = thirdDim;
  return (central12(A,x,x) + central12(A,y,y) + central12(A,z,z))
    * recip_2;
}

template<class T>
inline _bz_typename T::T_numtype::T_numtype
div3D4n(T& A) {
  const int x = firstDim, y = secondDim, z = thirdDim;
  return (central14(A,x,x) + central14(A,y,y) + central14(A,z,z))
    * recip_12;
}

/****************************************************************************
 * Mixed Partial derivatives
 ****************************************************************************/

template<class T>
inline _bz_typename T::T_numtype
mixed22(T& A, int x, int y)
{
    return A.shift(-1,x,-1,y) - A.shift(-1,x,1,y)
        -A.shift(1,x,-1,y) + A.shift(1,x,1,y);
}

template<class T>
inline _bz_typename T::T_numtype
mixed22n(T& A, int x, int y)
{
    return mixed22(A, x, y) * recip_4;
}

template<class T>
inline _bz_typename T::T_numtype
mixed24(T& A, int x, int y)
{
    return 64.*(A.shift(-1,x,-1,y) - A.shift(-1,x,1,y)
        -A.shift(1,x,-1,y) + A.shift(1,x,1,y))
        + (A.shift(-2,x,+1,y) - A.shift(-1,x,2,y)
        - A.shift(1,x,2,y)-A.shift(2,x,1,y)
        + A.shift(2,x,-1,y)+A.shift(1,x,-2,y)
        - A.shift(-1,x,-2,y)+A.shift(-2,x,-1,y))
        + 8.*(A.shift(-1,x,1,y)+A.shift(-1,x,2,y)
        -A.shift(2,x,-2,y) + A.shift(2,x,2,y));
}

template<class T>
inline _bz_typename T::T_numtype
mixed24n(T& A, int x, int y)
{
    return mixed24(A,x,y) * recip_144;
}

/****************************************************************************
 * Smoothers
 ****************************************************************************/

// NEEDS_WORK-- put other stencil operators here:
//   Average5pt2D
//   Average7pt3D
// etc.

/****************************************************************************
 * Stencil operators with geometry (experimental)
 ****************************************************************************/

template<class T>
inline _bz_typename multicomponent_traits<_bz_typename
    T::T_numtype>::T_element div3DVec4(T& A, 
    const UniformCubicGeometry<3>& geom)
{
    const int x = 0, y = 1, z = 2;

    return (central14(A, x, firstDim) + central14(A, y, secondDim)
        + central14(A, z, thirdDim)) * recip_12 * geom.recipSpatialStep();
}

template<class T>
inline _bz_typename T::T_numtype Laplacian3D4(T& A, 
    const UniformCubicGeometry<3>& geom)
{
    return Laplacian3D4n(A) * geom.recipSpatialStepPow2();
}

template<class T>
inline _bz_typename T::T_numtype Laplacian3DVec4(T& A,
    const UniformCubicGeometry<3>& geom)
{
    typedef _bz_typename T::T_numtype vector3d;
    typedef _bz_typename multicomponent_traits<vector3d>::T_element 
        T_element;
    const int u = 0, v = 1, w = 2;
    const int x = 0, y = 1, z = 2;

    // central24 is a 5-point stencil
    // This is a 9*5 = 45 point stencil

    T_element t1 = (central24(A,u,x) + central24(A,u,y) + central24(A,u,z))
        * recip_12 * geom.recipSpatialStepPow2();

    T_element t2 = (central24(A,v,x) + central24(A,v,y) + central24(A,v,z))
        * recip_12 * geom.recipSpatialStepPow2();

    T_element t3 = (central24(A,w,x) + central24(A,w,y) + central24(A,w,z))
        * recip_12 * geom.recipSpatialStepPow2();

    return vector3d(t1,t2,t3);
}

template<class T>
inline TinyMatrix<_bz_typename multicomponent_traits<_bz_typename
    T::T_numtype>::T_element, 3, 3>
grad3DVec4(T& A, const UniformCubicGeometry<3>& geom)
{
    const int x=0, y=1, z=2;
    const int u=0, v=1, w=2;

    TinyMatrix<_bz_typename multicomponent_traits<_bz_typename
        T::T_numtype>::T_element, 3, 3> grad;

    // This is a 9*4 = 36 point stencil
    grad(u,x) = central14n(A,u,x) * geom.recipSpatialStep();
    grad(u,y) = central14n(A,u,y) * geom.recipSpatialStep();
    grad(u,z) = central14n(A,u,z) * geom.recipSpatialStep();
    grad(v,x) = central14n(A,v,x) * geom.recipSpatialStep();
    grad(v,y) = central14n(A,v,y) * geom.recipSpatialStep();
    grad(v,z) = central14n(A,v,z) * geom.recipSpatialStep();
    grad(w,x) = central14n(A,w,x) * geom.recipSpatialStep();
    grad(w,y) = central14n(A,w,y) * geom.recipSpatialStep();
    grad(w,z) = central14n(A,w,z) * geom.recipSpatialStep();

    return grad;
}

template<class T>
inline TinyVector<_bz_typename T::T_numtype,3> grad3D4(T& A,
    const UniformCubicGeometry<3>& geom) {
  return TinyVector<_bz_typename T::T_numtype,3>(
    central14(A,firstDim) * recip_12 * geom.recipSpatialStep(),
    central14(A,secondDim) * recip_12 * geom.recipSpatialStep(),
    central14(A,thirdDim) * recip_12 * geom.recipSpatialStep());
}

BZ_NAMESPACE_END

#endif // BZ_ARRAYSTENCILOPS_H