stencilops.h

著名的数学计算类库

  return 2.0 * (*A)[comp] - 5.0 * A.shift(-1,dim)[comp]
    + 4.0 * A.shift(-2,dim)[comp] - A.shift(-3,dim)[comp];
}

BZ_DECLARE_MULTIDIFF(backward32) {
  return 5.0 * (*A)[comp] - 18.0 * A.shift(-1,dim)[comp]
    + 24.0 * A.shift(-2,dim)[comp] - 14.0 * A.shift(-3,dim)[comp]
    + 3.0 * A.shift(-4,dim)[comp];
}

BZ_DECLARE_MULTIDIFF(backward42) {
  return 3.0 * (*A)[comp] - 14.0 * A.shift(-1,dim)[comp]
    + 26.0 * A.shift(-2,dim)[comp] - 24.0 * A.shift(-3,dim)[comp]
    + 11.0 * A.shift(-4,dim)[comp] - 2.0 * A.shift(-5,dim)[comp];
}

/****************************************************************************
 * Backward differences with accuracy O(h^2)  (normalized)
 ****************************************************************************/

BZ_DECLARE_DIFF(backward12n) { return backward12(A,dim) * recip_2; }
BZ_DECLARE_DIFF(backward22n) { return backward22(A,dim); }
BZ_DECLARE_DIFF(backward32n) { return backward32(A,dim) * recip_2; }
BZ_DECLARE_DIFF(backward42n) { return backward42(A,dim); }

/****************************************************************************
 * Backward differences with accuracy O(h^2)  (normalized, multicomponent)
 ****************************************************************************/

BZ_DECLARE_MULTIDIFF(backward12n) { return backward12(A,comp,dim) * recip_2; }
BZ_DECLARE_MULTIDIFF(backward22n) { return backward22(A,comp,dim); }
BZ_DECLARE_MULTIDIFF(backward32n) { return backward32(A,comp,dim) * recip_2; }
BZ_DECLARE_MULTIDIFF(backward42n) { return backward42(A,comp,dim); }

/****************************************************************************
 * Forward differences with accuracy O(h)  
 ****************************************************************************/

BZ_DECLARE_DIFF(forward11) {
  return -(*A) + A.shift(1,dim);
}

BZ_DECLARE_DIFF(forward21) {
  return (*A) - 2.0 * A.shift(1,dim) + A.shift(2,dim);
}

BZ_DECLARE_DIFF(forward31) {
  return -(*A) + 3.0 * A.shift(1,dim) - 3.0 * A.shift(2,dim) + A.shift(3,dim);
}

BZ_DECLARE_DIFF(forward41) {
  return (*A) - 4.0 * A.shift(1,dim) + 6.0 * A.shift(2,dim)
    - 4.0 * A.shift(3,dim) + A.shift(4,dim);
}

/****************************************************************************
 * Forward differences with accuracy O(h)  (multicomponent versions)
 ****************************************************************************/

BZ_DECLARE_MULTIDIFF(forward11) {
  return  -(*A)[comp] + A.shift(1,dim)[comp];
}

BZ_DECLARE_MULTIDIFF(forward21) {
  return (*A)[comp] - 2.0 * A.shift(1,dim)[comp] + A.shift(2,dim)[comp];
}

BZ_DECLARE_MULTIDIFF(forward31) {
  return -(*A)[comp] + 3.0 * A.shift(1,dim)[comp] - 3.0 * A.shift(2,dim)[comp]
    + A.shift(3,dim)[comp];
}

BZ_DECLARE_MULTIDIFF(forward41) {
  return (*A)[comp] - 4.0 * A.shift(1,dim)[comp] + 6.0 * A.shift(2,dim)[comp]
    - 4.0 * A.shift(3,dim)[comp] + A.shift(4,dim)[comp];
}

/****************************************************************************
 * Forward differences with accuracy O(h)     (normalized)
 ****************************************************************************/

BZ_DECLARE_DIFF(forward11n) { return forward11(A,dim); }
BZ_DECLARE_DIFF(forward21n) { return forward21(A,dim); }
BZ_DECLARE_DIFF(forward31n) { return forward31(A,dim); }
BZ_DECLARE_DIFF(forward41n) { return forward41(A,dim); }

/****************************************************************************
 * Forward differences with accuracy O(h)     (multicomponent,normalized)
 ****************************************************************************/

BZ_DECLARE_MULTIDIFF(forward11n) { return forward11(A,comp,dim); }
BZ_DECLARE_MULTIDIFF(forward21n) { return forward21(A,comp,dim); }
BZ_DECLARE_MULTIDIFF(forward31n) { return forward31(A,comp,dim); }
BZ_DECLARE_MULTIDIFF(forward41n) { return forward41(A,comp,dim); }

/****************************************************************************
 * Forward differences with accuracy O(h^2)     
 ****************************************************************************/

BZ_DECLARE_DIFF(forward12) {
  return -3.0 * (*A) + 4.0 * A.shift(1,dim) - A.shift(2,dim);
}

BZ_DECLARE_DIFF(forward22) {
  return 2.0 * (*A) - 5.0 * A.shift(1,dim) + 4.0 * A.shift(2,dim)
    - A.shift(3,dim);
}

BZ_DECLARE_DIFF(forward32) {
  return -5.0 * (*A) + 18.0 * A.shift(1,dim) - 24.0 * A.shift(2,dim) 
    + 14.0 * A.shift(3,dim) - 3.0 * A.shift(4,dim);
}

BZ_DECLARE_DIFF(forward42) {
  return 3.0 * (*A) - 14.0 * A.shift(1,dim) + 26.0 * A.shift(2,dim)
    - 24.0 * A.shift(3,dim) + 11.0 * A.shift(4,dim) - 2.0 * A.shift(5,dim);
}

/****************************************************************************
 * Forward differences with accuracy O(h^2)   (multicomponent versions)
 ****************************************************************************/

BZ_DECLARE_MULTIDIFF(forward12) {
  return -3.0 * (*A)[comp] + 4.0 * A.shift(1,dim)[comp] - A.shift(2,dim)[comp];
}

BZ_DECLARE_MULTIDIFF(forward22) {
  return 2.0 * (*A)[comp] - 5.0 * A.shift(1,dim)[comp]
    + 4.0 * A.shift(2,dim)[comp] - A.shift(3,dim)[comp];
}

BZ_DECLARE_MULTIDIFF(forward32) {
  return -5.0 * (*A)[comp] + 18.0 * A.shift(1,dim)[comp]
    - 24.0 * A.shift(2,dim)[comp] + 14.0 * A.shift(3,dim)[comp]
    - 3.0 * A.shift(4,dim)[comp];
}

BZ_DECLARE_MULTIDIFF(forward42) {
  return 3.0 * (*A)[comp] - 14.0 * A.shift(1,dim)[comp]
    + 26.0 * A.shift(2,dim)[comp] - 24.0 * A.shift(3,dim)[comp]
    + 11.0 * A.shift(4,dim)[comp] - 2.0 * A.shift(5,dim)[comp];
}


/****************************************************************************
 * Forward differences with accuracy O(h^2)     (normalized)
 ****************************************************************************/

BZ_DECLARE_DIFF(forward12n) { return forward12(A,dim) * recip_2; }
BZ_DECLARE_DIFF(forward22n) { return forward22(A,dim); }
BZ_DECLARE_DIFF(forward32n) { return forward32(A,dim) * recip_2; }
BZ_DECLARE_DIFF(forward42n) { return forward42(A,dim); }

/****************************************************************************
 * Forward differences with accuracy O(h^2)     (normalized)
 ****************************************************************************/

BZ_DECLARE_MULTIDIFF(forward12n) { return forward12(A,comp,dim) * recip_2; }
BZ_DECLARE_MULTIDIFF(forward22n) { return forward22(A,comp,dim); }
BZ_DECLARE_MULTIDIFF(forward32n) { return forward32(A,comp,dim) * recip_2; }
BZ_DECLARE_MULTIDIFF(forward42n) { return forward42(A,comp,dim); }

/****************************************************************************
 * Gradient operators
 ****************************************************************************/

template<typename T>
inline TinyVector<_bz_typename T::T_numtype,2> grad2D(T& A) {
  return TinyVector<_bz_typename T::T_numtype,2>(
    central12(A,firstDim),
    central12(A,secondDim));
}

template<typename T>
inline TinyVector<_bz_typename T::T_numtype,2> grad2D4(T& A) {
  return TinyVector<_bz_typename T::T_numtype,2>(
    central14(A,firstDim),
    central14(A,secondDim));
}

template<typename T>
inline TinyVector<_bz_typename T::T_numtype,3> grad3D(T& A) {
  return TinyVector<_bz_typename T::T_numtype,3>(
    central12(A,firstDim),
    central12(A,secondDim),
    central12(A,thirdDim));
}

template<typename T>
inline TinyVector<_bz_typename T::T_numtype,3> grad3D4(T& A) {
  return TinyVector<_bz_typename T::T_numtype,3>(
    central14(A,firstDim),
    central14(A,secondDim),
    central14(A,thirdDim));
}

template<typename T>
inline TinyVector<_bz_typename T::T_numtype,2> grad2Dn(T& A) {
  return TinyVector<_bz_typename T::T_numtype,2>(
    central12n(A,firstDim),
    central12n(A,secondDim));
}

template<typename T>
inline TinyVector<_bz_typename T::T_numtype,2> grad2D4n(T& A) {
  return TinyVector<_bz_typename T::T_numtype,2>(
    central14n(A,firstDim),
    central14n(A,secondDim));
}

template<typename T>
inline TinyVector<_bz_typename T::T_numtype,3> grad3Dn(T& A) {
  return TinyVector<_bz_typename T::T_numtype,3>(
    central12n(A,firstDim),
    central12n(A,secondDim),
    central12n(A,thirdDim));
}

template<typename T>
inline TinyVector<_bz_typename T::T_numtype,3> grad3D4n(T& A) {
  return TinyVector<_bz_typename T::T_numtype,3>(
    central14n(A,firstDim),
    central14n(A,secondDim),
    central14n(A,thirdDim));
}

/****************************************************************************
 * Grad-squared operators
 ****************************************************************************/

template<typename T>
inline TinyVector<_bz_typename T::T_numtype,2> gradSqr2D(T& A) {
  return TinyVector<_bz_typename T::T_numtype,2>(
    central22(A,firstDim),
    central22(A,secondDim));
}

template<typename T>
inline TinyVector<_bz_typename T::T_numtype,2> gradSqr2D4(T& A) {
  return TinyVector<_bz_typename T::T_numtype,2>(
    central24(A,firstDim),
    central24(A,secondDim));
}

template<typename T>
inline TinyVector<_bz_typename T::T_numtype,3> gradSqr3D(T& A) {
  return TinyVector<_bz_typename T::T_numtype,3>(
    central22(A,firstDim),
    central22(A,secondDim),
    central22(A,thirdDim));
}

template<typename T>
inline TinyVector<_bz_typename T::T_numtype,3> gradSqr3D4(T& A) {
  return TinyVector<_bz_typename T::T_numtype,3>(
    central24(A,firstDim),
    central24(A,secondDim),
    central24(A,thirdDim));
}

/****************************************************************************
 * Grad-squared operators (normalized)
 ****************************************************************************/

template<typename T>
inline TinyVector<_bz_typename T::T_numtype,2> gradSqr2Dn(T& A) {
  return gradSqr2D(A);
}

template<typename T>
inline TinyVector<_bz_typename T::T_numtype,2> gradSqr2D4n(T& A) {
  return TinyVector<_bz_typename T::T_numtype,2>(
    central24(A,firstDim) * recip_12,
    central24(A,secondDim) * recip_12);
}

template<typename T>
inline TinyVector<_bz_typename T::T_numtype,3> gradSqr3Dn(T& A) {
  return gradSqr3D(A);
}

template<typename T>
inline TinyVector<_bz_typename T::T_numtype,3> gradSqr3D4n(T& A) {
  return TinyVector<_bz_typename T::T_numtype,3>(
    central24(A,firstDim) * recip_12,
    central24(A,secondDim) * recip_12,
    central24(A,thirdDim) * recip_12);
}

/****************************************************************************
 * Gradient operators on a vector field
 ****************************************************************************/

template<typename T>
inline TinyMatrix<_bz_typename multicomponent_traits<_bz_typename 
    T::T_numtype>::T_element, 3, 3>
Jacobian3D(T& A)
{
    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;

    grad(u,x) = central12(A,u,x);
    grad(u,y) = central12(A,u,y);
    grad(u,z) = central12(A,u,z);
    grad(v,x) = central12(A,v,x);
    grad(v,y) = central12(A,v,y);
    grad(v,z) = central12(A,v,z);
    grad(w,x) = central12(A,w,x);
    grad(w,y) = central12(A,w,y);
    grad(w,z) = central12(A,w,z);

    return grad;
}

template<typename T>
inline TinyMatrix<_bz_typename multicomponent_traits<_bz_typename 
    T::T_numtype>::T_element, 3, 3>
Jacobian3Dn(T& A)
{
    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;
    
    grad(u,x) = central12n(A,u,x);
    grad(u,y) = central12n(A,u,y);
    grad(u,z) = central12n(A,u,z);
    grad(v,x) = central12n(A,v,x);
    grad(v,y) = central12n(A,v,y);
    grad(v,z) = central12n(A,v,z);
    grad(w,x) = central12n(A,w,x);
    grad(w,y) = central12n(A,w,y);
    grad(w,z) = central12n(A,w,z);

    return grad;
}

template<typename T>
inline TinyMatrix<_bz_typename multicomponent_traits<_bz_typename
    T::T_numtype>::T_element, 3, 3>
Jacobian3D4(T& A)
{
    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;
    
    grad(u,x) = central14(A,u,x);
    grad(u,y) = central14(A,u,y);
    grad(u,z) = central14(A,u,z);
    grad(v,x) = central14(A,v,x);
    grad(v,y) = central14(A,v,y);
    grad(v,z) = central14(A,v,z);
    grad(w,x) = central14(A,w,x);
    grad(w,y) = central14(A,w,y);
    grad(w,z) = central14(A,w,z);

    return grad;
}

template<typename T>
inline TinyMatrix<_bz_typename multicomponent_traits<_bz_typename
    T::T_numtype>::T_element, 3, 3>
Jacobian3D4n(T& A)
{
    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;
    
    grad(u,x) = central14n(A,u,x);
    grad(u,y) = central14n(A,u,y);
    grad(u,z) = central14n(A,u,z);
    grad(v,x) = central14n(A,v,x);
    grad(v,y) = central14n(A,v,y);
    grad(v,z) = central14n(A,v,z);
    grad(w,x) = central14n(A,w,x);
    grad(w,y) = central14n(A,w,y);
    grad(w,z) = central14n(A,w,z);

    return grad;
}

/****************************************************************************
 * Curl operators
 ****************************************************************************/

// O(h^2) curl, using central difference

template<typename T>
inline TinyVector<_bz_typename T::T_numtype,3> 
curl(T& vx, T& vy, T& vz) {