maia::math Namespace Reference

Typedefs
using	MFloatMatrix2d = Eigen::Matrix< MFloat, 2, 2 >
	all assume dense matrix calculations More...
using	MFloatMatrix3d = Eigen::Matrix< MFloat, 3, 3 >
using	MFloatMatrixXd = Eigen::Matrix< MFloat, Eigen::Dynamic, Eigen::Dynamic >
template<MInt nDim>
using	MFloatMatrix = Eigen::Matrix< MFloat, nDim, nDim >
using	MFloatVector2d = Eigen::Vector< MFloat, 2 >
using	MFloatVector3d = Eigen::Vector< MFloat, 3 >
using	MFloatVectorXd = Eigen::Vector< MFloat, Eigen::Dynamic >
template<MInt nDim>
using	MFloatVector = Eigen::Vector< MFloat, nDim >
using	MTriplet = Eigen::Triplet< MFloat >

Functions
MInt	getGlobalPosFFTW (MInt i0, MInt i1, MInt i2, MInt ny, MInt nz)
MFloat	triadImag (fftw_complex &a, fftw_complex &b, fftw_complex &c)
MFloat	randnormal (MFloat mu, MFloat sigma)
std::complex< MFloat > *	getFourierCoefficient (MFloat *k, MFloat k0, const MFloat Ma)
std::complex< MFloat > *	getFourierCoefficient2 (MFloat *k, MFloat k0)
void	initFftFilter (MFloat *uPhysFieldCoarse, MFloat *vPhysFieldCoarse, MFloat *wPhysFieldCoarse, MInt lx, MInt ly, MInt lz, MInt lxCoarse, MInt lyCoarse, MInt lzCoarse, MInt noPeakModes, const MFloat Ma)
	Generates a velocity field from Fourier-modes using FFTW The disturbances are generated and filtered to a coarser resolution. More...
void	initFft (fftw_complex *uPhysField, fftw_complex *vPhysField, fftw_complex *wPhysField, MInt lx, MInt ly, MInt lz, MInt noPeakModes, const MFloat Ma)
	Generates a velocity field from Fourier-modes using FFTW. More...
void	computeEnergySpectrum (MFloatScratchSpace &velocity, MIntScratchSpace &indices, const ptrdiff_t howmany, MInt locSize, MInt nx, MInt ny, MInt nz, MFloat viscosity, const MPI_Comm comm)
	Compute energy spectrum on unity cube. More...
template<typename T , std::size_t N>
Eigen::MatrixXd	ConvertToEigenMatrix (std::array< std::array< T, N >, N > data)
	convert array to eigen matrix More...
template<typename T , std::size_t N>
void	adjointRow (std::array< std::array< T, N >, N > &m, std::array< T, N > &A, const MInt r)
template void	adjointRow< MFloat, 3 > (std::array< std::array< MFloat, 3 >, 3 > &m, std::array< MFloat, 3 > &A, const MInt)
template<typename T >
void	adjoint1stRow4x4 (std::array< std::array< T, 4 >, 4 > &m, std::array< T, 4 > &A)
template void	adjoint1stRow4x4< MFloat > (std::array< std::array< MFloat, 4 >, 4 > &m, std::array< MFloat, 4 > &A)
template<typename T , std::size_t N>
void	adjoint1stRow (std::array< std::array< T, N >, N > &m, std::array< T, N > &A)
template void	adjoint1stRow< MFloat, 4 > (std::array< std::array< MFloat, 4 >, 4 > &m, std::array< MFloat, 4 > &A)
template void	adjoint1stRow< MFloat, 3 > (std::array< std::array< MFloat, 3 >, 3 > &m, std::array< MFloat, 3 > &A)
template<typename T , std::size_t N>
MFloat	determinant (std::array< T, N > &m)
template MFloat	determinant< MFloat, 2ul > (std::array< MFloat, 2ul > &)
template MFloat	determinant< MFloat, 3ul > (std::array< MFloat, 3ul > &)
template MFloat	determinant< MFloat, 4ul > (std::array< MFloat, 4ul > &)
template<typename T , std::size_t N>
MFloat	determinant (std::array< std::array< T, N >, N > &m)
template MFloat	determinant< MFloat, 2ul > (std::array< std::array< MFloat, 2ul >, 2ul > &)
template MFloat	determinant< MFloat, 3ul > (std::array< std::array< MFloat, 3ul >, 3ul > &)
template MFloat	determinant< MFloat, 4ul > (std::array< std::array< MFloat, 4ul >, 4ul > &)
void	invert (MFloat *A, const MInt m, const MInt n)
template<class T >
void	invert (T &A, T &AInv, const MInt m, const MInt n)
template void	invert< ScratchSpace< MFloat > > (ScratchSpace< MFloat > &, ScratchSpace< MFloat > &, const MInt, const MInt)
template void	invert< maia::tensor::Tensor< MFloat > > (maia::tensor::Tensor< MFloat > &, maia::tensor::Tensor< MFloat > &, const MInt, const MInt)
template<class T >
void	invert (T &A, T &weights, T &AInv, const MInt m, const MInt n)
template void	invert< ScratchSpace< MFloat > > (ScratchSpace< MFloat > &, ScratchSpace< MFloat > &, ScratchSpace< MFloat > &, const MInt, const MInt)
template void	invert< maia::tensor::Tensor< MFloat > > (maia::tensor::Tensor< MFloat > &, maia::tensor::Tensor< MFloat > &, maia::tensor::Tensor< MFloat > &, const MInt, const MInt)
template<class T >
MInt	invertR (T &A, T &weights, T &AInv, const MInt m, const MInt n)
template MInt	invertR< ScratchSpace< MFloat > > (ScratchSpace< MFloat > &, ScratchSpace< MFloat > &, ScratchSpace< MFloat > &, const MInt, const MInt)
void	solveSparseMatrix (MFloat *A_coeff, MInt **pos, const MInt n, const MInt m, MFloat *b, MFloat *x)
void	solveSparseMatrixIterative (MFloat *A_coeff, MInt **pos, const MInt n, const MInt m, MFloat *b, MFloat *x)
void	solveDenseMatrix (MFloat *A_coeff, MInt **pos, const MInt n, const MInt m, MFloat *b, MFloat *x)
void	multiplySparseMatrixVector (MFloat *A_coeff, MInt **pos, const MInt n, const MInt m, MFloat *b, MFloat *x)
template<MInt nDim>
void	solveQR (std::array< std::array< MFloat, nDim >, nDim > &A_, std::array< MFloat, nDim > &b_)
template void	solveQR< 2 > (std::array< std::array< MFloat, 2 >, 2 > &, std::array< MFloat, 2 > &)
template void	solveQR< 3 > (std::array< std::array< MFloat, 3 >, 3 > &, std::array< MFloat, 3 > &)
void	calcEigenValues (MFloat A[3][3], MFloat w[3])
void	calcEigenValues (MFloat **A_in, MFloat *lambda_in, const MInt m)
void	calcEigenVectors (MFloat A[3][3], MFloat Q[3][3], MFloat w[3])
std::vector< MFloat >	svd (MFloat *const A, MFloat *const b, const MInt m, const MInt n, MFloat *const x)
void	quickSortImpl (MInt *a, MInt start, MInt end)
void	quickSort (MInt *a, MInt start, MInt end)
template<MInt nDim, typename T , typename U >
T *	vecAdd (T *const result, const U *const a)
	Linear Algebra. More...
template<MInt nDim, typename T , typename U , typename... Ts>
T *	vecAdd (T *const result, const U *const a, const Ts... b)
template<MInt nDim, typename T , typename... Ts>
void	vecAvg (T *const M, const Ts *const ... args)
	Computes the arithmetic mean of an arbritary number of vectors of dimension nDim and returns it in M. More...
template<typename T >
void	cross (const T *const u, const T *const v, T *const c)
template<typename T >
std::array< T, 3 >	cross (const std::array< T, 3 > &u, const std::array< T, 3 > &v)
template<typename T >
T	cross (const std::array< T, 2 > &u, const std::array< T, 2 > &v)
template<typename T >
T *	cross (const T(&u)[3], const T(&v)[3])
template<typename T >
T	cross (const T(&u)[2], const T(&v)[2])
template<typename T , std::size_t N>
MFloat	norm (const std::array< T, N > &u)
template<typename T >
MFloat	norm (const T *const u, const MInt N)
template<typename T , std::size_t N>
void	normalize (std::array< T, N > &u)
template<typename T >
void	normalize (T *const u, const MInt N)
template<typename T , std::size_t N>
std::array< T, N >	normalized (const std::array< T, N > &u)
template<MInt nDim>
MFloat	distance (const MFloat *a, const MFloat *b)
MFloat	distance (const std::array< MFloat, 2 > a, const std::array< MFloat, 2 > b)
MFloat	distance (const std::array< MFloat, 3 > a, const std::array< MFloat, 3 > b)
MFloat	distance (const MFloat *const a, const std::array< MFloat, 2 > b)
MFloat	distance (const MFloat *const a, const std::array< MFloat, 3 > b)
MFloat	distance (const std::array< MFloat, 2 > a, const MFloat *const b)
MFloat	distance (std::array< MFloat, 3 > a, const MFloat *const b)
void	multiplyMatricesSq (MFloatScratchSpace &m1, MFloatScratchSpace &m2, MFloatScratchSpace &result, MInt dim)
void	multiplyMatrices (MFloatScratchSpace &m1, MFloatScratchSpace &m2, MFloatScratchSpace &result, MInt m1_n, MInt m1_m, MInt m2_n, MInt m2_m)
void	addMatrices (MFloatScratchSpace &m1, MFloatScratchSpace &m2, MFloatScratchSpace &result, MInt dim1, MInt dim2)
MFloat	frobeniusMatrixNormSquared (MFloatScratchSpace &m, MInt dim1, MInt dim2)
MFloat	frobeniusMatrixNorm (MFloatScratchSpace &m, MInt dim1, MInt dim2)
MFloat	besselJ0 (MFloat x)
MFloat	besselJ1 (MFloat x)
MFloat	besselJi (MInt order, MFloat x)
MFloat	lincos (MFloat arg)
void	sortEigenVectors (MFloat A[3][3], MFloat w[3])
	Sorts the eigenvalues and the associated eigenvectors from large to small. More...
void	quatMult (const MFloat *const qA, const MFloat *const qB, MFloat *const qC)
template<typename T >
MInt	sgn (T val)
MFloat	crankAngle (const MFloat time, const MFloat Strouhal, const MFloat offset, const MInt mode)
	help-function for engine calculations which returns the crank-angle for a given time, Strouhals-number and crank-angle offset mode = 0: return CAD in range of (0-720) mode = 1: return accumulated crankAnge in radian More...
void	computeRotationMatrix (MFloatScratchSpace &R, MFloat *q)
	rotation matrix co-rotating(~inertial) frame -> body-fixed frame More...
void	rotation2quaternion (MFloat *rotation, MFloat *quaternion)
void	matrixVectorProduct (MFloat *c, MFloatScratchSpace &A, MFloat *b)
	c=A*b More...
void	matrixVectorProductTranspose (MFloat *c, MFloatScratchSpace &A, MFloat *b)
	c=A^t*b More...
MInt	inverse (MFloat **a, MFloat **ainv, MInt n, const MFloat epsilon)
MInt	quickSortPartition (MInt *a, MInt start, MInt end)
MInt	removeDoubleEntries (MInt *a, MInt size)
void	calculateLegendrePolyAndDeriv (MInt Nmax, MFloat x, MFloat *polynomial, MFloat *derivative)
	Evaluates the Legendre polynomial and its derivative of degree Nmax at point x. More...
void	calculateLegendreGaussNodesAndWeights (MInt Nmax, MFloat *nodes, MFloat *wInt)
	Calculate the Gauss integration nodes and weight for the Legendre polynomials on the interval [-1,1]. More...
MFloat	RBF (const MFloat R, const MFloat R0)
	radial base function More...
MFloat	deltaFun (const MFloat r, const MFloat r0, const MFloat r1)
std::vector< MFloat >	linSpace (const MFloat start, const MFloat end, const MInt num)
MFloat	getSector (MFloat y, MFloat z, MFloat azimuthalAngle)
MFloat	getAngle (MFloat y, MFloat z)

Typedef Documentation

MFloatMatrix

template<MInt nDim>

using maia::math::MFloatMatrix = typedef Eigen::Matrix<MFloat, nDim, nDim>

Definition at line 54 of file maiamath.cpp.

MFloatMatrix2d

using maia::math::MFloatMatrix2d = typedef Eigen::Matrix<MFloat, 2, 2>

Definition at line 50 of file maiamath.cpp.

MFloatMatrix3d

using maia::math::MFloatMatrix3d = typedef Eigen::Matrix<MFloat, 3, 3>

Definition at line 51 of file maiamath.cpp.

MFloatMatrixXd

using maia::math::MFloatMatrixXd = typedef Eigen::Matrix<MFloat, Eigen::Dynamic, Eigen::Dynamic>

Definition at line 52 of file maiamath.cpp.

MFloatVector

template<MInt nDim>

using maia::math::MFloatVector = typedef Eigen::Vector<MFloat, nDim>

Definition at line 60 of file maiamath.cpp.

MFloatVector2d

using maia::math::MFloatVector2d = typedef Eigen::Vector<MFloat, 2>

Definition at line 56 of file maiamath.cpp.

MFloatVector3d

using maia::math::MFloatVector3d = typedef Eigen::Vector<MFloat, 3>

Definition at line 57 of file maiamath.cpp.

MFloatVectorXd

using maia::math::MFloatVectorXd = typedef Eigen::Vector<MFloat, Eigen::Dynamic>

Definition at line 58 of file maiamath.cpp.

MTriplet

using maia::math::MTriplet = typedef Eigen::Triplet<MFloat>

Definition at line 61 of file maiamath.cpp.

Function Documentation

addMatrices()

void maia::math::addMatrices ( MFloatScratchSpace &  m1, MFloatScratchSpace &  m2, MFloatScratchSpace &  result, MInt  dim1, MInt  dim2  )

inline

Definition at line 319 of file maiamath.h.

                                   {
  for(MInt i = 0; i < dim1; i++) {
    for(MInt j = 0; j < dim2; j++) {
      result(i, j) = m1(i, j) + m2(i, j);
    }
  }
}

adjoint1stRow()

template<typename T , std::size_t N>

void maia::math::adjoint1stRow	(	std::array< std::array< T, N >, N > &	m,
		std::array< T, N > &	A
	)

First row of adjoint of a matrix

Template Parameters

T

Type

Parameters

m	Matrix (two dimensional array containing the matrix elements)

Returns

First row of the adjuncte of the matrix m in A

Definition at line 114 of file maiamath.cpp.

                                                                      {
  if constexpr(N == 4) {
    adjoint1stRow4x4(m, A);
  } else {
    adjointRow(m, A, 0);
  }
}

adjoint1stRow4x4()

template<typename T >

void maia::math::adjoint1stRow4x4	(	std::array< std::array< T, 4 >, 4 > &	m,
		std::array< T, 4 > &	A
	)

First row of adjoint of a 4 x 4 matrix

Template Parameters

T

Type

Parameters

m	Matrix (two dimensional array containing the matrix elements)

Returns

First row of the adjuncte of the 4 x 4 matrix m in A

Definition at line 97 of file maiamath.cpp.

                                                                         {
  A[0] = m[1][1] * m[2][2] * m[3][3] + m[1][2] * m[2][3] * m[3][1] + m[1][3] * m[2][1] * m[3][2]
         - m[1][3] * m[2][2] * m[3][1] - m[1][2] * m[2][1] * m[3][3] - m[1][1] * m[2][3] * m[3][2];
  A[1] = -m[0][1] * m[2][2] * m[3][3] - m[0][2] * m[2][3] * m[3][1] - m[0][3] * m[2][1] * m[3][2]
         + m[0][3] * m[2][2] * m[3][1] + m[0][2] * m[2][1] * m[3][3] + m[0][1] * m[2][3] * m[3][2];
  A[2] = m[0][1] * m[1][2] * m[3][3] + m[0][2] * m[1][3] * m[3][1] + m[0][3] * m[1][1] * m[3][2]
         - m[0][3] * m[1][2] * m[3][1] - m[0][2] * m[1][1] * m[3][3] - m[0][1] * m[1][3] * m[3][2];
  A[3] = -m[0][1] * m[1][2] * m[2][3] - m[0][2] * m[1][3] * m[2][1] - m[0][3] * m[1][1] * m[2][2]
         + m[0][3] * m[1][2] * m[2][1] + m[0][2] * m[1][1] * m[2][3] + m[0][1] * m[1][3] * m[2][2];
}

adjoint1stRow4x4< MFloat >()

template void maia::math::adjoint1stRow4x4< MFloat >	(	std::array< std::array< MFloat, 4 >, 4 > &	m,
		std::array< MFloat, 4 > &	A
	)

adjoint1stRow< MFloat, 3 >()

template void maia::math::adjoint1stRow< MFloat, 3 >	(	std::array< std::array< MFloat, 3 >, 3 > &	m,
		std::array< MFloat, 3 > &	A
	)

adjoint1stRow< MFloat, 4 >()

template void maia::math::adjoint1stRow< MFloat, 4 >	(	std::array< std::array< MFloat, 4 >, 4 > &	m,
		std::array< MFloat, 4 > &	A
	)

adjointRow()

template<typename T , std::size_t N>

void maia::math::adjointRow	(	std::array< std::array< T, N >, N > &	m,
		std::array< T, N > &	A,
		const MInt	r
	)

Row of adjoint of a matrix

Template Parameters

T	Type
N	Dimension of the matrix

Parameters

m	Matrix (one dimensional array containing the matrix elements)
r	Row of the adjuncte

Returns

Given row of the adjuncte of the matrix m in A

Definition at line 82 of file maiamath.cpp.

                                                                                 {
  Eigen::Matrix<T, N, N> m2 = ConvertToEigenMatrix(m);
  Eigen::Matrix<T, N, N> adjMatrix(m2.adjoint());
  for(std::size_t i = 0; i < N; i++) {
    A[i] = adjMatrix.row(r)[i];
  }
}

adjointRow< MFloat, 3 >()

template void maia::math::adjointRow< MFloat, 3 >	(	std::array< std::array< MFloat, 3 >, 3 > &	m,
		std::array< MFloat, 3 > &	A,
		const	MInt
	)

besselJ0()

MFloat maia::math::besselJ0 ( MFloat  x )

inline

Definition at line 343 of file maiamath.h.

                                 {
  // This subroutine calculates the First Kind Bessel Function of
  // order 0, for any real number X. The polynomial approximation by
  // series of Chebyshev polynomials is used for 0<X<8 and 0<8/X<1.
  // REFERENCES:
  // M.ABRAMOWITZ,I.A.STEGUN, HANDBOOK OF MATHEMATICAL FUNCTIONS, 1965.
  // C.W.CLENSHAW, NATIONAL PHYSICAL LABORATORY MATHEMATICAL TABLES,
  // VOL.5, 1962.
 
  const MFloat p1 = 1.0, p2 = -0.1098628627e-2, p3 = 0.2734510407e-4, p4 = -0.2073370639e-5, p5 = 0.2093887211e-6,
               q1 = -0.1562499995e-1, q2 = 0.1430488765e-3, q3 = -0.6911147651e-5, q4 = 0.7621095161e-6,
               q5 = -0.9349451520e-7, r1 = 57568490574.0, r2 = -13362590354.0, r3 = 651619640.7, r4 = -11214424.18,
               r5 = 77392.33017, r6 = -184.9052456, s1 = 57568490411.0, s2 = 1029532985.0, s3 = 9494680.718,
               s4 = 59272.64853, s5 = 267.8532712, s6 = 1.0;
  MFloat ax, fr, fs, z, fp, fq1, xx, y, tmp;
 
  if(approx(x, 0.0, std::numeric_limits<MFloat>::epsilon())) {
    return 1.0;
  }
 
  ax = fabs(x);
  if(ax < 8.0) {
    y = x * x;
    fr = r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * r6))));
    fs = s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6))));
    tmp = fr / fs;
  } else {
    z = 8. / ax;
    y = z * z;
    xx = ax - 0.785398164;
    fp = p1 + y * (p2 + y * (p3 + y * (p4 + y * p5)));
    fq1 = q1 + y * (q2 + y * (q3 + y * (q4 + y * q5)));
    tmp = std::sqrt(0.636619772 / ax) * (fp * cos(xx) - z * fq1 * sin(xx));
  }
  return tmp;
}

besselJ1()

MFloat maia::math::besselJ1 ( MFloat  x )

inline

Definition at line 380 of file maiamath.h.

                                 {
  // This subroutine calculates the First Kind Bessel Function of
  // order 1, for any real number X. The polynomial approximation by
  // series of Chebyshev polynomials is used for 0<X<8 and 0<8/X<1.
  // REFERENCES:
  // M.ABRAMOWITZ,I.A.STEGUN, HANDBOOK OF MATHEMATICAL FUNCTIONS, 1965.
  // C.W.CLENSHAW, NATIONAL PHYSICAL LABORATORY MATHEMATICAL TABLES,
  // VOL.5, 1962.
 
  const MFloat p1 = 1.0, p2 = 0.183105e-2, p3 = -0.3516396496e-4, p4 = 0.2457520174e-5, p5 = -0.240337019e-6,
               p6 = 0.636619772, q1 = 0.04687499995, q2 = -0.2002690873e-3, q3 = 0.8449199096e-5, q4 = -0.88228987e-6,
               q5 = 0.105787412e-6, r1 = 72362614232.0, r2 = -7895059235.0, r3 = 242396853.1, r4 = -2972611.439,
               r5 = 15704.48260, r6 = -30.16036606, s1 = 144725228442.0, s2 = 2300535178.0, s3 = 18583304.74,
               s4 = 99447.43394, s5 = 376.9991397, s6 = 1.0;
 
  MFloat ax, fr, fs, y, z, fp, fq1, xx, tmp;
 
  ax = fabs(x);
  if(ax < 8.0) {
    y = x * x;
    fr = r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * r6))));
    fs = s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6))));
    tmp = x * (fr / fs);
  } else {
    z = 8.0 / ax;
    y = z * z;
    xx = ax - 2.35619491;
    fp = p1 + y * (p2 + y * (p3 + y * (p4 + y * p5)));
    fq1 = q1 + y * (q2 + y * (q3 + y * (q4 + y * q5)));
    tmp = sqrt(p6 / ax) * (cos(xx) * fp - z * sin(xx) * fq1) * (x < 0.0 ? -fabs(s6) : fabs(s6));
  }
  return tmp;
}

besselJi()

MFloat maia::math::besselJi ( MInt  order, MFloat  x  )

inline

Definition at line 415 of file maiamath.h.

                                             {
  if(order == 0) {
    return besselJ0(x);
  }
  if(order == 1) {
    return besselJ1(x);
  }
  mTerm(1, AT_, "invalid order");
}

calcEigenValues() [1/2]

void maia::math::calcEigenValues	(	MFloat **	A_in,
		MFloat *	lambda_in,
		const MInt	m
	)

Get the Eigenvalues of matrix A_in (symmetrical)

Parameters

A_in	system matrix
lambda_in	vector of eigen values
m	size of square matrix A

Author

Moritz Waldmann

Definition at line 489 of file maiamath.cpp.

                                                                     {
  Eigen::Map<Eigen::MatrixXd> A(&A_in[0][0], m, m);
  Eigen::Map<Eigen::VectorXd> lambda(lambda_in, m);
 
  // Assume matrix is symmetrical
  Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es(A, Eigen::EigenvaluesOnly);
  lambda = es.eigenvalues().real();
}

calcEigenValues() [2/2]

void maia::math::calcEigenValues	(	MFloat	A[3][3],
		MFloat	w[3]
	)

Get the Eigenvalues of matrix A (symmetrical)

Parameters

A	system matrix
w	Eigenvalues

Author

Sven Berger

Definition at line 475 of file maiamath.cpp.

                                                  {
  Eigen::Map<Eigen::Matrix<MFloat, 3, 3>> mA(&A[0][0]);
  Eigen::Map<Eigen::Matrix<MFloat, 3, 1>> _w(w);
 
  // Assume matrix is symmetrical
  Eigen::SelfAdjointEigenSolver<Eigen::Matrix<MFloat, 3, 3>> es(mA, Eigen::EigenvaluesOnly);
  _w = es.eigenvalues().real();
}

calcEigenVectors()

void maia::math::calcEigenVectors	(	MFloat	A[3][3],
		MFloat	Q[3][3],
		MFloat	w[3]
	)

Get the Eigenvalues and Eigenvectors of matrix A (symmetrical)

Parameters

A	system matrix
w	Eigenvalues
Q	Eigenvectors

Author

Sven Berger

Definition at line 503 of file maiamath.cpp.

                                                                   {
  Eigen::Map<Eigen::Matrix<MFloat, 3, 3>> mA(&A[0][0]);
  Eigen::Map<Eigen::Matrix<MFloat, 3, 3>> mQ(&Q[0][0]);
  Eigen::Map<Eigen::Matrix<MFloat, 3, 1>> _w(w);
 
  // Assume matrix is symmetrical
  Eigen::SelfAdjointEigenSolver<Eigen::Matrix<MFloat, 3, 3>> es(mA);
  mQ = es.eigenvectors().real();
  _w = es.eigenvalues().real();
}

calculateLegendreGaussNodesAndWeights()

void maia::math::calculateLegendreGaussNodesAndWeights ( MInt  Nmax, MFloat *  nodes, MFloat *  wInt  )

inline

Author

Michael Schlottke (mic) mic@a.nosp@m.ia.r.nosp@m.wth-a.nosp@m.ache.nosp@m.n.de

Date

2012-11-07

Parameters

[in]	Nmax	Maximum degree of the polynomials.
[out]	nodes	The resulting integration nodes.
[out]	wInt	The resulting integration weights.

Taken from Kopriva09, p. 64, algorithm 23.

Definition at line 806 of file maiamath.h.

                                                                                          {
  TRACE();
 
  // Reset nodes and weights
  const MInt noNodes = Nmax + 1;
  std::fill_n(&nodes[0], noNodes, F0);
  std::fill_n(&wInt[0], noNodes, F0);
 
  // Set tolerance and number of iterations. According to Kopriva09,
  // 1.0E-15 and 10 should be more than sufficient.
  const MFloat tol = F4 * MFloatEps;
  const MInt noIterations = 10;
 
  // Catch simple cases before going into the full loop
  if(Nmax == 0) {
    nodes[0] = F0;
    wInt[0] = F2;
  } else if(Nmax == 1) {
    nodes[0] = -sqrt(F1B3);
    wInt[0] = 1;
    nodes[1] = -nodes[0];
    wInt[1] = wInt[0];
  } else {
    // Use symmetry property of the roots of the Legendre polynomials
    for(MInt j = 0; j < (Nmax + 1) / 2; j++) {
      // Calculate starting guess for Newton method
      nodes[j] = -cos((F2 * j + F1) / (F2 * Nmax + F2) * PI);
 
      // Use Newton method to find root of Legendre polynomial
      // -> this is also the integration node
      MFloat poly, deriv;
      for(MInt k = 0; k < noIterations; k++) {
        calculateLegendrePolyAndDeriv(Nmax + 1, nodes[j], &poly, &deriv);
        MFloat delta = -poly / deriv;
        nodes[j] += delta;
 
        // Stop iterations if error is small enough
        if(fabs(delta) <= tol * fabs(nodes[j])) {
          break;
        }
      }
 
      // Calculate weight
      calculateLegendrePolyAndDeriv(Nmax + 1, nodes[j], &poly, &deriv);
      wInt[j] = F2 / ((1 - nodes[j] * nodes[j]) * deriv * deriv);
 
      // Set nodes and weights according to symmetry properties
      nodes[Nmax - j] = -nodes[j];
      wInt[Nmax - j] = wInt[j];
    }
  }
 
  // If odd number of nodes (noNodes = Nmax + 1), set center node to
  // origin (0.0) and calculate weight
  if(Nmax % 2 == 0) {
    MFloat poly, deriv;
    calculateLegendrePolyAndDeriv(Nmax + 1, F0, &poly, &deriv);
    nodes[Nmax / 2] = F0;
    wInt[Nmax / 2] = F2 / (deriv * deriv);
  }
}

calculateLegendrePolyAndDeriv()

void maia::math::calculateLegendrePolyAndDeriv ( MInt  Nmax, MFloat  x, MFloat *  polynomial, MFloat *  derivative  )

inline

Author

Michael Schlottke (mic) mic@a.nosp@m.ia.r.nosp@m.wth-a.nosp@m.ache.nosp@m.n.de

Date

2012-11-06

Parameters

[in]	Nmax	The polynomial degree.
[in]	x	The evaluation point.
[out]	polynomial	The resulting value of the Legendre polynomial.
[out]	derivative	The resulting value of the derivative.

Taken from Kopriva09, p. 63, algorithm 22.

Definition at line 754 of file maiamath.h.

                                                                                                       {
  // TRACE();
 
  // Create temporary storage for the polynomial and its derivative
  MFloat poly = F0;
  MFloat deriv = F0;
 
  // Explicitly calculate the first two values of the three term recursion
  if(Nmax == 0) {
    poly = F1;
    deriv = F0;
  } else if(Nmax == 1) {
    poly = x;
    deriv = F1;
  } else {
    MFloat polyLast1, polyLast2, derivLast1, derivLast2;
 
    polyLast2 = F1;
    polyLast1 = x;
    derivLast2 = F0;
    derivLast1 = F1;
 
    // Calculate the polynomial and its derivative for higher degrees
    // (Nmax >= 2)
    for(MInt k = 2; k <= Nmax; k++) {
      poly = (F2 * k - F1) / k * x * polyLast1 - (k - F1) / k * polyLast2;
      deriv = derivLast2 + (F2 * k - F1) * polyLast1;
      polyLast2 = polyLast1;
      polyLast1 = poly;
      derivLast2 = derivLast1;
      derivLast1 = deriv;
    }
  }
 
  // Save results to pointer locations
  *polynomial = poly;
  *derivative = deriv;
}

computeEnergySpectrum()

void maia::math::computeEnergySpectrum ( MFloatScratchSpace &  velocity, MIntScratchSpace &  indices, const ptrdiff_t  howmany, MInt  locSize, MInt  nx, MInt  ny, MInt  nz, MFloat  viscosity, const MPI_Comm  comm  )

inline

Author

Lennart Schneiders

discrete spectrum is computed on unity cube, no spatial scaling required velocity field is computed for u_rms = 1, hence uPhysField subsequently has to be scaled by the magnitude of the fluctuations, e.g. uPhysField *= m_UInfinity

(u,v,w) PhysField: complex velocity in physical space

Parameters

uPhysField	pointer to a fftw_complex for u
vPhysField	pointer to a fftw_complex for v
wPhysField	pointer to a fftw_complex for w

Definition at line 1051 of file maiafftw.h.

                                                                                                                  {
  DEBUG("computeEnergySpectrum entry", MAIA_DEBUG_TRACE_IN);
  TRACE();
 
  const MInt size = nx * ny * nz;
  const MFloat fsize = (MFloat)size;
  const ptrdiff_t rank = 3;
  if(nx % 2 != 0 || ny % 2 != 0 || nz % 2 != 0) {
    std::stringstream errorMessage;
    errorMessage << " FFTInit: domainsize must NOT be an odd number! " << nx << " " << ny << " " << nz << std::endl;
    mTerm(1, AT_, errorMessage.str());
  }
  if(howmany < 3) mTerm(1, AT_, "expecting at least 3D velocity field");
 
  m_log << " --- initializing FFTW --- " << std::endl;
  m_log << " domain size = " << nx << "x" << ny << "x" << nz << " (" << howmany << ")" << std::endl;
 
  // caution: nz must be evenly divisble by number of ranks!
  MInt domainId;
  MInt noDomains;
  MPI_Comm_rank(comm, &domainId);
  MPI_Comm_size(comm, &noDomains);
  MPI_Comm MPI_COMM_FFTW;
  MInt maxRank = 1;
  maxRank = mMin(nz, noDomains);
  while(nz % maxRank != 0) {
    maxRank--;
  }
  MInt color = (domainId < maxRank) ? 0 : MPI_UNDEFINED;
  MPI_Comm_split(comm, color, domainId, &MPI_COMM_FFTW, AT_, "MPI_COMM_FFTW");
  m_log << "no ranks for fft: " << maxRank << std::endl;
 
  MFloat urms0 = F0;
  for(MInt i = 0; i < locSize; i++) {
    urms0 += POW2(velocity(i, 0)) + POW2(velocity(i, 1)) + POW2(velocity(i, 2));
  }
  MPI_Allreduce(MPI_IN_PLACE, &urms0, 1, MPI_DOUBLE, MPI_SUM, comm, AT_, "MPI_IN_PLACE", "urms0");
  urms0 = sqrt(F1B3 * urms0 / fsize);
 
  MInt fftLocSize = 0;
  ptrdiff_t alloc_local = 0, local_n0 = 0, local_0_start = 0;
  fftw_complex* hatField = nullptr;
  fftw_plan plan;
 
  const MFloat time0 = MPI_Wtime();
 
  const ptrdiff_t n[3] = {nx, ny, nz};
#ifndef MAIA_WINDOWS
  alloc_local = (domainId < maxRank) ? fftw_mpi_local_size_many(rank, n, howmany, FFTW_MPI_DEFAULT_BLOCK, MPI_COMM_FFTW,
                                                                &local_n0, &local_0_start)
                                     : 1;
#endif
 
  ScratchSpace<fftw_complex> hatFieldMem(alloc_local, AT_, "hatFieldMem");
  if(domainId < maxRank) {
    hatField = &hatFieldMem[0];
  } else {
    alloc_local = 0;
    local_n0 = 0;
    local_0_start = nx;
  }
 
  ScratchSpace<MLong> offsets(noDomains + 1, AT_, "offsets");
  ScratchSpace<MInt> noRecv(noDomains, AT_, "noRecv");
  ScratchSpace<MInt> noSend(noDomains, AT_, "noSend");
  noRecv.fill(0);
  noSend.fill(0);
  MLong locOffset = (domainId < maxRank) ? ((MInt)local_0_start) * ny * nz : size;
 
  MPI_Allgather(&locOffset, 1, MPI_LONG, &offsets[0], 1, MPI_LONG, comm, AT_, "locOffset", "offsets[0]");
  offsets(noDomains) = size;
  if(domainId < maxRank) {
    if(offsets(domainId) != locOffset || offsets(domainId + 1) != ((MInt)(local_0_start + local_n0)) * ny * nz)
      mTerm(1, AT_, "wrong size 0");
  }
  for(MInt i = 0; i < locSize; i++) {
    MInt j = indices(i);
    MInt nghbrDomain = mMin(maxRank - 1, j / (size / maxRank));
    while(j < offsets(nghbrDomain) || j >= offsets(nghbrDomain + 1)) {
      if(j < offsets(nghbrDomain)) nghbrDomain--;
      if(j >= offsets(nghbrDomain + 1)) nghbrDomain++;
    }
    if(nghbrDomain >= maxRank) mTerm(1, AT_, "wrong domain");
    noSend(nghbrDomain)++;
  }
  MPI_Alltoall(&noSend[0], 1, MPI_INT, &noRecv[0], 1, MPI_INT, comm, AT_, "noSend[0]", "noRecv[0]");
  MInt sendSize = 0;
  MInt recvSize = 0;
  ScratchSpace<MInt> recvOffsets(noDomains + 1, AT_, "recvOffsets");
  ScratchSpace<MInt> sendOffsets(noDomains + 1, AT_, "sendOffsets");
  for(MInt i = 0; i < noDomains; i++) {
    recvOffsets(i) = recvSize;
    sendOffsets(i) = sendSize;
    sendSize += noSend(i);
    recvSize += noRecv(i);
  }
  recvOffsets(noDomains) = recvSize;
  sendOffsets(noDomains) = sendSize;
  ScratchSpace<MFloat> recvData(mMax(1, recvSize), howmany, AT_, "recvData");
  ScratchSpace<MFloat> sendData(mMax(1, sendSize), howmany, AT_, "sendData");
  ScratchSpace<MInt> recvIndices(mMax(1, recvSize), AT_, "recvIndices");
  ScratchSpace<MInt> sendIndices(mMax(1, sendSize), AT_, "sendIndices");
  noSend.fill(0);
  for(MInt i = 0; i < locSize; i++) {
    MInt j = indices(i);
    MInt nghbrDomain = mMin(maxRank - 1, j / (size / maxRank));
    while(j < offsets(nghbrDomain) || j >= offsets(nghbrDomain + 1)) {
      if(j < offsets(nghbrDomain)) nghbrDomain--;
      if(j >= offsets(nghbrDomain + 1)) nghbrDomain++;
    }
    if(nghbrDomain >= maxRank) mTerm(1, AT_, "wrong domain");
    if(sendOffsets(nghbrDomain) + noSend(nghbrDomain) >= sendOffsets(nghbrDomain + 1)) mTerm(1, AT_, "wrong size");
    for(MInt k = 0; k < howmany; k++) {
      sendData(sendOffsets(nghbrDomain) + noSend(nghbrDomain), k) = velocity(i, k);
    }
    sendIndices(sendOffsets(nghbrDomain) + noSend(nghbrDomain)) = j;
    if(j < offsets(nghbrDomain) || j >= offsets(nghbrDomain + 1)) mTerm(1, AT_, "wrong size 01");
    noSend(nghbrDomain)++;
  }
  for(MInt i = 0; i < noDomains; i++) {
    if(noSend(i) != sendOffsets(i + 1) - sendOffsets(i)) mTerm(1, AT_, "wrong size 1");
    if(noRecv(i) != recvOffsets(i + 1) - recvOffsets(i)) mTerm(1, AT_, "wrong size 1");
  }
 
  ScratchSpace<MPI_Request> sendReq(noDomains, AT_, "sendReq");
  sendReq.fill(MPI_REQUEST_NULL);
  MInt cnt = 0;
  for(MInt i = 0; i < noDomains; i++) {
    if(noSend(i) == 0) continue;
    MPI_Issend(&sendData[howmany * sendOffsets(i)], howmany * noSend(i), MPI_DOUBLE, i, 123, comm, &sendReq[cnt], AT_,
               "sendData[howmany * sendOffsets(i)]");
    cnt++;
  }
  for(MInt i = 0; i < noDomains; i++) {
    if(noRecv(i) == 0) continue;
    MPI_Recv(&recvData[howmany * recvOffsets(i)], howmany * noRecv(i), MPI_DOUBLE, i, 123, comm, MPI_STATUS_IGNORE, AT_,
             "recvData[howmany * recvOffsets(i)]");
  }
  if(cnt > 0) MPI_Waitall(cnt, &sendReq[0], MPI_STATUSES_IGNORE, AT_);
  sendReq.fill(MPI_REQUEST_NULL);
  cnt = 0;
  for(MInt i = 0; i < noDomains; i++) {
    if(noSend(i) == 0) continue;
    MPI_Issend(&sendIndices[sendOffsets(i)], noSend(i), MPI_INT, i, 124, comm, &sendReq[cnt], AT_,
               "sendIndices[sendOffsets(i)]");
    cnt++;
  }
  for(MInt i = 0; i < noDomains; i++) {
    if(noRecv(i) == 0) continue;
    MPI_Recv(&recvIndices[recvOffsets(i)], noRecv(i), MPI_INT, i, 124, comm, MPI_STATUS_IGNORE, AT_,
             "recvIndices[recvOffsets(i)]");
  }
  if(cnt > 0) MPI_Waitall(cnt, &sendReq[0], MPI_STATUSES_IGNORE, AT_);
 
  if(domainId < maxRank) {
    for(MInt i = 0; i < recvSize; i++) {
      MInt j = recvIndices(i);
      if(j < ((MInt)local_0_start) * ny * nz || j >= ((MInt)(local_0_start + local_n0)) * ny * nz)
        mTerm(1, AT_, "wrong size 2");
      MInt pos = j - (((MInt)local_0_start) * ny * nz);
      for(MInt k = 0; k < howmany; k++) {
        hatField[howmany * pos + k][0] = recvData(i, k); // Real part of fourier transformed velocity
        hatField[howmany * pos + k][1] = F0;             // Imaginary part of fourier transformed velocity
      }
    }
    fftLocSize = recvSize;
 
    MFloat couplcheck = F0;
    for(MInt i = 0; i < fftLocSize; i++) {
      for(MInt k = 0; k < 3; k++) {
        couplcheck -= hatField[howmany * i + k][0] * hatField[howmany * i + 3 + k][0];
      }
    }
    MPI_Allreduce(MPI_IN_PLACE, &couplcheck, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_FFTW, AT_, "MPI_IN_PLACE", "couplcheck");
    if(domainId == 0)
      std::cerr << "couplecheck fft " << couplcheck << " " << couplcheck * POW3(0.25 / 64.0) << std::endl;
  }
 
  // Perform Fourier transform
  if(domainId < maxRank) {
#ifndef MAIA_WINDOWS
    plan = fftw_mpi_plan_many_dft(rank, n, howmany, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, hatField, hatField,
                                  MPI_COMM_FFTW, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute(plan);
    fftw_destroy_plan(plan);
#endif
    for(MInt i = 0; i < fftLocSize; i++) {
      for(MInt j = 0; j < howmany; j++) {
        hatField[howmany * i + j][0] /= fsize;
        hatField[howmany * i + j][1] /= fsize;
      }
    }
  }
 
  const MFloat time1 = MPI_Wtime();
 
  m_log << "parallel fftw time " << time1 - time0 << std::endl;
 
  // Set coefficients:
  // FFTW stores the coefficients for positive wavenumbers in the first half of the array,
  // and those for negative wavenumbers in reverse order in the second half.
  // [0, 1, ... , N/2-1, N/2, ... , N-1]
  //  - the entry for zero-wavenumber is at position 0
  //  - the k-th entry and the (N-k)th entry correspond to wavenumbers with opposite sign
  //  - the entry at position N/2 corresponds to the Nyquist wavenumber and appears only once
 
  static const MInt referenceCubeSize =
      ((Context::propertyExists("referenceCubeSize")) ? Context::getBasicProperty<MInt>("referenceCubeSize", AT_) : 1);
 
  static const MInt kMinSpec = ((Context::propertyExists("kMin")) ? Context::getBasicProperty<MInt>("kMin", AT_) : 0);
 
  static const MInt kMaxSpec = ((Context::propertyExists("kMax")) ? Context::getBasicProperty<MInt>("kMax", AT_) : 0);
 
  const MFloat k0 = F2 * PI / referenceCubeSize;
  MInt noBins = 0;
  MFloat kMin = 0.0;
  MFloat kMax = 0.0;
  MFloat deltaK = 0.0;
  if(kMinSpec > 0 && kMaxSpec > 0) {
    noBins = kMaxSpec;
    kMin = static_cast<float>(kMinSpec) * k0;
    kMax = static_cast<float>(kMaxSpec) * k0;
    deltaK = (kMax - kMin) / ((MFloat)(kMaxSpec - 1));
  } else {
    noBins = mMax(nx, mMax(ny, nz)) / 2;
    kMin = F1 * k0;
    kMax = noBins * k0;
    deltaK = (kMax - kMin) / ((MFloat)(noBins - 1));
  }
 
  // const MInt noBins = 200;
  // const MFloat kMax = 200.0*k0;
  ScratchSpace<MFloat> spectrumBnd(noBins + 1, AT_, "spectrumBnd");
  ScratchSpace<MFloat> spectrum(noBins, AT_, "spectrum");
  ScratchSpace<MFloat> coupling(noBins, AT_, "coupling");
  ScratchSpace<MFloat> transfer(noBins, AT_, "transfer");
  ScratchSpace<MFloat> transfer2(noBins, AT_, "transfer2");
  ScratchSpace<MFloat> rate(noBins, AT_, "rate");
  ScratchSpace<MFloat> flux(noBins, 2, AT_, "flux");
  ScratchSpace<MFloat> flux2(noBins, 2, AT_, "flux2");
  ScratchSpace<MFloat> spectrumSum(noBins, AT_, "spectrumSum");
  ScratchSpace<MFloat> spectrumSum2(noBins, AT_, "spectrumSum2");
  spectrumBnd[0] = kMin - F1B2 * deltaK;
  const MFloat spectrumBnd00 = kMin - F3B4 * deltaK;
  for(MInt i = 0; i < noBins; i++) {
    spectrumBnd[i + 1] = spectrumBnd[i] + deltaK;
  }
  spectrum.fill(F0);
  coupling.fill(F0);
  transfer.fill(F0);
  transfer2.fill(F0);
  rate.fill(F0);
  flux.fill(F0);
  flux2.fill(F0);
  spectrumSum.fill(F0);
  spectrumSum2.fill(F0);
 
  // beware: very expensive routine
  static const MBool computeTransferRate =
      ((Context::propertyExists("computeTransferRate")) ? Context::getBasicProperty<MBool>("computeTransferRate", AT_)
                                                        : false);
  if(computeTransferRate) {
    MIntScratchSpace noDatPerDomain(noDomains, AT_, "noDatPerDomain");
    MIntScratchSpace dataOffsets(noDomains + 1, AT_, "dataOffsets");
    MPI_Allgather(&fftLocSize, 1, MPI_INT, &noDatPerDomain[0], 1, MPI_INT, comm, AT_, "fftLocSize",
                  "noDatPerDomain[0]");
    dataOffsets.fill(0);
    for(MInt d = 0; d < noDomains; d++) {
      dataOffsets(d + 1) = dataOffsets(d) + noDatPerDomain(d);
    }
 
    MPI_Request sendReqLoc = MPI_REQUEST_NULL;
    fftw_complex* uHatGlobal = fftw_alloc_complex(3 * size);
    fftw_complex* uHatGlobalTmp = nullptr;
    for(MInt i = 0; i < fftLocSize; i++) {
      for(MInt j = 0; j < 3; j++) {
        uHatGlobal[3 * i + j][0] = hatField[howmany * i + j][0];
        uHatGlobal[3 * i + j][1] = hatField[howmany * i + j][1];
      }
    }
    if(fftLocSize > 0) {
      MPI_Issend(&uHatGlobal[0], 6 * fftLocSize, MPI_DOUBLE, 0, 18, comm, &sendReqLoc, AT_, "uHatGlobal[0]");
    }
    if(domainId == 0) {
      uHatGlobalTmp = fftw_alloc_complex(3 * size);
      for(MInt d = 0; d < noDomains; d++) {
        if(noDatPerDomain[d] > 0) {
          MPI_Recv(&(uHatGlobalTmp[3 * dataOffsets[d]][0]), 6 * noDatPerDomain[d], MPI_DOUBLE, d, 18, comm,
                   MPI_STATUSES_IGNORE, AT_, "(uHatGlobalTmp[3 * dataOffsets[d]][0])");
        }
      }
    }
    if(fftLocSize > 0) {
      MPI_Wait(&sendReqLoc, MPI_STATUSES_IGNORE, AT_);
    }
 
    if(domainId == 0) {
      // re-sort by ascending wave numbers
      for(MInt i0 = 0; i0 < nx; i0++) {
        for(MInt i1 = 0; i1 < ny; i1++) {
          for(MInt i2 = 0; i2 < nz; i2++) {
            MInt j0 = (i0 > nx / 2) ? i0 - nx : i0;
            MInt j1 = (i1 > ny / 2) ? i1 - ny : i1;
            MInt j2 = (i2 > nz / 2) ? i2 - nz : i2;
            j0 += nx / 2 - 1;
            j1 += nx / 2 - 1;
            j2 += nx / 2 - 1;
            MInt pos0 = getGlobalPosFFTW(i0, i1, i2, ny, nz);
            MInt pos1 = getGlobalPosFFTW(j0, j1, j2, ny, nz);
            for(MInt j = 0; j < 3; j++) {
              uHatGlobal[3 * pos1 + j][0] = uHatGlobalTmp[3 * pos0 + j][0];
              uHatGlobal[3 * pos1 + j][1] = uHatGlobalTmp[3 * pos0 + j][1];
            }
          }
        }
      }
      fftw_free(uHatGlobalTmp);
    }
 
    MPI_Bcast(&uHatGlobal[0], 6 * size, MPI_DOUBLE, 0, comm, AT_, "uHatGlobal");
 
 
    // MInt domainSize = size/noDomains;
    // MInt posMin = domainSize*domainId;
    // MInt posMax = mMin(size,domainSize*(domainId+1));
    // cerr << domainId << ": pos " << posMin << " " << posMax << std::endl;
 
    MInt noWaveNumbers = 0;
    MLong weightSum = 0;
    for(MInt pos = 0; pos < size; pos++) {
      MInt i2 = pos % nz;
      MInt i1 = (pos / nz) % ny;
      MInt i0 = pos / (ny * nz);
      if(pos != i2 + nz * (i1 + ny * (i0))) mTerm(1, AT_, "index mismatch");
      MFloat k[3];
      MInt j0 = i0 - nx / 2 + 1;
      MInt j1 = i1 - ny / 2 + 1;
      MInt j2 = i2 - nz / 2 + 1;
      k[0] = ((MFloat)j0) * k0;
      k[1] = ((MFloat)j1) * k0;
      k[2] = ((MFloat)j2) * k0;
 
      MFloat kAbs = sqrt(k[0] * k[0] + k[1] * k[1] + k[2] * k[2]);
      if(kAbs > spectrumBnd[noBins] + 1e-12) continue;
      if(kAbs < spectrumBnd[0] - 1e-12) continue;
 
      MInt p0Min = std::max(0, i0 - nx / 2);
      MInt p0Max = std::min(nx, i0 + nx / 2);
      MInt p1Min = std::max(0, i1 - ny / 2);
      MInt p1Max = std::min(ny, i1 + ny / 2);
      MInt p2Min = std::max(0, i2 - nz / 2);
      MInt p2Max = std::min(nz, i2 + nz / 2);
      MLong weight = ((p0Max - p0Min) * (p1Max - p1Min) * (p2Max - p2Min));
 
      weightSum += weight;
      noWaveNumbers++;
    }
    // MInt domainSize = noWaveNumbers/noDomains;
    MLong noDomainsLong = (MLong)noDomains;
    MLong domainIdLong = (MLong)domainId;
    MLong domainSize = weightSum / noDomainsLong;
    MInt posMin = 0;
    MInt posMax = size;
    noWaveNumbers = 0;
    weightSum = 0;
    for(MInt pos = 0; pos < size; pos++) {
      MInt i2 = pos % nz;
      MInt i1 = (pos / nz) % ny;
      MInt i0 = pos / (ny * nz);
      if(pos != i2 + nz * (i1 + ny * (i0))) mTerm(1, AT_, "index mismatch");
      MFloat k[3];
      MInt j0 = i0 - nx / 2 + 1;
      MInt j1 = i1 - ny / 2 + 1;
      MInt j2 = i2 - nz / 2 + 1;
      k[0] = ((MFloat)j0) * k0;
      k[1] = ((MFloat)j1) * k0;
      k[2] = ((MFloat)j2) * k0;
 
      MFloat kAbs = sqrt(k[0] * k[0] + k[1] * k[1] + k[2] * k[2]);
      if(kAbs > spectrumBnd[noBins] + 1e-12) continue;
      if(kAbs < spectrumBnd[0] - 1e-12) continue;
 
      MInt p0Min = std::max(0, i0 - nx / 2);
      MInt p0Max = std::min(nx, i0 + nx / 2);
      MInt p1Min = std::max(0, i1 - ny / 2);
      MInt p1Max = std::min(ny, i1 + ny / 2);
      MInt p2Min = std::max(0, i2 - nz / 2);
      MInt p2Max = std::min(nz, i2 + nz / 2);
      MLong weight = ((p0Max - p0Min) * (p1Max - p1Min) * (p2Max - p2Min));
 
      if(weightSum <= (domainSize * domainIdLong) && (weightSum + weight) > (domainSize * domainIdLong)) posMin = pos;
      if(weightSum <= (domainSize * (domainIdLong + 1)) && (weightSum + weight) > (domainSize * (domainIdLong + 1)))
        posMax = pos;
 
      weightSum += weight;
      noWaveNumbers++;
    }
 
    weightSum = 0;
    for(MInt pos = posMin; pos < posMax; pos++) {
      const MInt i2 = pos % nz;
      const MInt i1 = (pos / nz) % ny;
      const MInt i0 = pos / (ny * nz);
 
      const MInt j0 = i0 - nx / 2 + 1;
      const MInt j1 = i1 - ny / 2 + 1;
      const MInt j2 = i2 - nz / 2 + 1;
      const MFloat k[3] = {((MFloat)j0) * k0, ((MFloat)j1) * k0, ((MFloat)j2) * k0};
      const MFloat kAbs2 = mMax(1e-12, k[0] * k[0] + k[1] * k[1] + k[2] * k[2]);
      const MFloat kAbs = mMax(1e-12, sqrt(k[0] * k[0] + k[1] * k[1] + k[2] * k[2]));
 
      if(kAbs > spectrumBnd[noBins] + 1e-12) continue;
      const MInt bin = ((MInt)((kAbs - spectrumBnd[0] + 1e-12) / deltaK));
      const MInt bin2 = ((MInt)((kAbs - spectrumBnd00 + 1e-12) / deltaK));
 
      const MFloat kjl[3][3] = {
          {1.0 - (k[0] * k[0] / kAbs2), 0.0 - (k[0] * k[1] / kAbs2), 0.0 - (k[0] * k[2] / kAbs2)},
          {0.0 - (k[1] * k[0] / kAbs2), 1.0 - (k[1] * k[1] / kAbs2), 0.0 - (k[1] * k[2] / kAbs2)},
          {0.0 - (k[2] * k[0] / kAbs2), 0.0 - (k[2] * k[1] / kAbs2), 1.0 - (k[2] * k[2] / kAbs2)}};
 
      const MInt p0Min = std::max(0, i0 - nx / 2);
      const MInt p0Max = std::min(nx, i0 + nx / 2);
      const MInt p1Min = std::max(0, i1 - ny / 2);
      const MInt p1Max = std::min(ny, i1 + ny / 2);
      const MInt p2Min = std::max(0, i2 - nz / 2);
      const MInt p2Max = std::min(nz, i2 + nz / 2);
 
      MFloat trans = F0;
 
      for(MInt p0 = p0Min; p0 < p0Max; p0++) {
        for(MInt p1 = p1Min; p1 < p1Max; p1++) {
          for(MInt p2 = p2Min; p2 < p2Max; p2++) {
            const MInt pos1 = p2 + nz * (p1 + ny * (p0));
            const MInt s0 = i0 - p0 + nx / 2 - 1;
            const MInt s1 = i1 - p1 + ny / 2 - 1;
            const MInt s2 = i2 - p2 + nz / 2 - 1;
            const MInt pos2 = s2 + nz * (s1 + ny * (s0));
 
            for(MInt i = 0; i < 3; i++) {
              for(MInt j = 0; j < 3; j++) {
                for(MInt l = 0; l < 3; l++) {
                  trans += k[i] * kjl[j][l]
                           * triadImag(uHatGlobal[3 * pos1 + l], uHatGlobal[3 * pos2 + i], uHatGlobal[3 * pos + j]);
                }
              }
            }
          }
        }
      }
 
      if(bin > -1 && bin < noBins) transfer(bin) += trans;
      if(bin2 > -1 && bin2 < noBins) transfer2(bin2) += trans;
    }
 
    MPI_Allreduce(MPI_IN_PLACE, &transfer[0], noBins, MPI_DOUBLE, MPI_SUM, comm, AT_, "MPI_IN_PLACE", "transfer[0]");
    MPI_Allreduce(MPI_IN_PLACE, &transfer2[0], noBins, MPI_DOUBLE, MPI_SUM, comm, AT_, "MPI_IN_PLACE", "transfer2[0]");
    fftw_free(uHatGlobal);
  }
 
  const MFloat time2 = MPI_Wtime();
 
  if(domainId < maxRank) {
    for(MInt i0 = (MInt)local_0_start; i0 < (MInt)(local_0_start + local_n0); i0++) {
      for(MInt i1 = 0; i1 < ny; i1++) {
        for(MInt i2 = 0; i2 < nz; i2++) {
          MFloat k[3];
          MInt j0 = (i0 > nx / 2) ? i0 - nx : i0;
          MInt j1 = (i1 > ny / 2) ? i1 - ny : i1;
          MInt j2 = (i2 > nz / 2) ? i2 - nz : i2;
          k[0] = ((MFloat)j0) * k0;
          k[1] = ((MFloat)j1) * k0;
          k[2] = ((MFloat)j2) * k0;
 
          MInt pos = i2 + nz * (i1 + ny * (i0 - ((MInt)local_0_start)));
          if(pos < 0 || pos >= fftLocSize) mTerm(1, AT_, "wrong size 3");
 
          MFloat kAbs = sqrt(k[0] * k[0] + k[1] * k[1] + k[2] * k[2]);
          MFloat eng = F0;
          MFloat coupl = F0;
          MFloat dedt = F0;
          for(MInt q = 0; q < 3; q++) {
            eng += F1B2 * (POW2(hatField[howmany * pos + q][0]) + POW2(hatField[howmany * pos + q][1]));
          }
 
          // if not only velocity field is passed as argument, also coupling rate Psi(k) and change of Energy dE/dt is
          // computed
          if(howmany > 3) {
            for(MInt q = 0; q < 3; q++) {
              coupl += (hatField[howmany * pos + q][0] * hatField[howmany * pos + 3 + q][0]
                        + hatField[howmany * pos + q][1] * hatField[howmany * pos + 3 + q][1]);
              dedt += (hatField[howmany * pos + 6 + q][0] * hatField[howmany * pos + q][0]
                       + hatField[howmany * pos + 6 + q][1] * hatField[howmany * pos + q][1]);
            }
          }
          if(kAbs > spectrumBnd[noBins] + 1e-12) continue;
 
          MInt bin = ((MInt)((kAbs - spectrumBnd[0] + 1e-12) / deltaK));
          MInt bin2 = ((MInt)((kAbs - spectrumBnd00 + 1e-12) / deltaK));
          if(bin > -1 && bin < noBins) {
            spectrum(bin) += eng;
            coupling(bin) += coupl;
            rate(bin) += dedt;
            spectrumSum(bin) += F1;
          }
          if(bin2 > -1 && bin2 < noBins) {
            spectrumSum2(bin2) += F1;
          }
        }
      }
    }
 
    MPI_Allreduce(MPI_IN_PLACE, &spectrum[0], noBins, MPI_DOUBLE, MPI_SUM, MPI_COMM_FFTW, AT_, "MPI_IN_PLACE",
                  "spectrum[0]");
    MPI_Allreduce(MPI_IN_PLACE, &coupling[0], noBins, MPI_DOUBLE, MPI_SUM, MPI_COMM_FFTW, AT_, "MPI_IN_PLACE",
                  "coupling[0]");
    MPI_Allreduce(MPI_IN_PLACE, &rate[0], noBins, MPI_DOUBLE, MPI_SUM, MPI_COMM_FFTW, AT_, "MPI_IN_PLACE", "rate[0]");
    MPI_Allreduce(MPI_IN_PLACE, &spectrumSum[0], noBins, MPI_DOUBLE, MPI_SUM, MPI_COMM_FFTW, AT_, "MPI_IN_PLACE",
                  "spectrumSum[0]");
    MPI_Allreduce(MPI_IN_PLACE, &spectrumSum2[0], noBins, MPI_DOUBLE, MPI_SUM, MPI_COMM_FFTW, AT_, "MPI_IN_PLACE",
                  "spectrumSum2[0]");
 
    for(MInt i = 0; i < noBins; i++) {
      MFloat db = spectrumBnd00 - spectrumBnd[0];
      MFloat vol =
          F4B3 * PI * (POW3(spectrumBnd[i + 1]) - POW3(spectrumBnd[i])) / (spectrumBnd[i + 1] - spectrumBnd[i]);
      MFloat vol2 = F4B3 * PI * (POW3(spectrumBnd[i + 1] + db) - POW3(spectrumBnd[i] + db))
                    / (spectrumBnd[i + 1] - spectrumBnd[i]);
      spectrum(i) *= vol / (spectrumSum(i) * POW3(k0));
      coupling(i) *= vol / (spectrumSum(i) * POW3(k0));
      transfer(i) *= vol / (spectrumSum(i) * POW3(k0));
      transfer2(i) *= vol2 / (spectrumSum2(i) * POW3(k0));
      rate(i) *= vol / (spectrumSum(i) * POW3(k0));
    }
 
    for(MInt i = 0; i < noBins; i++) {
      flux(i, 0) = F0;
      flux(i, 1) = F0;
      for(MInt j = 0; j < i; j++) {
        flux(i, 0) += transfer(j) * (spectrumBnd[j + 1] - spectrumBnd[j]) / k0;
      }
      for(MInt j = i; j < noBins; j++) {
        flux(i, 1) += transfer(j) * (spectrumBnd[j + 1] - spectrumBnd[j]) / k0;
      }
      flux2(i, 0) = F0;
      flux2(i, 1) = F0;
      for(MInt j = 0; j < i; j++) {
        flux2(i, 0) += transfer2(j) * (spectrumBnd[j + 1] - spectrumBnd[j]) / k0;
      }
      for(MInt j = i; j < noBins; j++) {
        flux2(i, 1) += transfer2(j) * (spectrumBnd[j + 1] - spectrumBnd[j]) / k0;
      }
    }
 
    if(domainId == 0) {
      std::ofstream ofl;
      std::ofstream ofl2;
      ofl.open("./out/energySpectrum_00" + std::to_string(globalTimeStep), std::ios_base::out | std::ios_base::trunc);
      ofl2.open("./out/energySpectrum_coarse_00" + std::to_string(globalTimeStep),
                std::ios_base::out | std::ios_base::trunc);
      if(ofl.is_open() && ofl.good() && ofl2.is_open() && ofl2.good()) {
        MFloat energy = F0;
        MFloat dissip = F0;
        MFloat length = F0;
        MFloat dedt = F0;
        MFloat force = F0;
        MFloat urms = F0;
        MFloat transf = F0;
        MFloat transf2 = F0;
        for(MInt i = 0; i < noBins; i++) {
          urms += spectrum(i) * (spectrumBnd[i + 1] - spectrumBnd[i]);
        }
        urms = sqrt(F2B3 * urms); // warum hier 2/3 und nicht 1/3
        for(MInt i = 0; i < noBins; i++) {
          MFloat k = F1B2 * (spectrumBnd[i + 1] + spectrumBnd[i]);
          MFloat Ek = spectrum(i);
          if(std::isnan(Ek) || std::isnan(coupling(i)) || std::isnan(rate(i))) continue;
          energy += (spectrumBnd[i + 1] - spectrumBnd[i]) * Ek;
          dissip += F2 * viscosity * POW2(k) * (spectrumBnd[i + 1] - spectrumBnd[i]) * Ek;
          length += (F1B2 * PI / POW2(urms)) * (spectrumBnd[i + 1] - spectrumBnd[i]) * Ek / (k);
          force += (spectrumBnd[i + 1] - spectrumBnd[i]) * coupling(i);
          dedt += (spectrumBnd[i + 1] - spectrumBnd[i]) * rate(i);
          transf += (spectrumBnd[i + 1] - spectrumBnd[i]) * transfer(i);
          transf2 += (spectrumBnd[i + 1] - spectrumBnd[i]) * transfer2(i);
        }
        ofl << "# Energy spectrum at time step " << globalTimeStep << ", total energy Ek=" << energy
            << ", dissipation rate=" << dissip << ", u_rms=" << urms << " (" << urms0
            << "), integral length scale=" << length << ", Kolmogorov time scale=" << sqrt(viscosity / dissip)
            << ", coupling=" << force << ", dedt=" << dedt << ", transfer=" << transf << ", transfer2=" << transf2
            << std::endl;
        ofl << "# 1: wave number k" << std::endl;
        ofl << "# 2: energy E(k)" << std::endl;
        ofl << "# 3: viscous dissipation rate D(k)" << std::endl;
        ofl << "# 4: transfer rate T(k)" << std::endl;
        ofl << "# 5: coupling rate Psi(k)" << std::endl;
        ofl << "# 6: flux -F(k)" << std::endl;
        ofl << "# 7: flux F(k)" << std::endl;
        ofl << "# 8: dE/dt(k)" << std::endl;
        ofl << "# 9: alternative transfer rate T(k)" << std::endl;
        ofl << "# 10: alternative flux -F(k)" << std::endl;
        ofl << "# 11: alternative flux F(k)" << std::endl;
        ofl << "# 12: number of samples in corresponding bin" << std::endl;
        ofl << "# 13: lower boundary of corresponding bin" << std::endl;
        ofl << "# 14: upper boundary of corresponding bin" << std::endl;
        ofl << "# 15: number of samples other spectrum: " << std::endl;
        ofl << "# 16: 4 * PI * k^2" << std::endl;
        ofl << "# 17: 4/3 * PI * (upperBndry^3 - lowerBndry^3)" << std::endl;
        for(MInt i = 0; i < noBins; i++) {
          MFloat k = F1B2 * (spectrumBnd[i + 1] + spectrumBnd[i]);
          MFloat Ek = spectrum(i);
          MFloat coupl = coupling(i);
          MFloat trans = transfer(i);
          MFloat trans2 = transfer2(i);
          ofl << k / k0 << " " << Ek << " " << F2 * POW2(k) * viscosity * Ek << " " << trans << " " << coupl << " "
              << flux(i, 0) << " " << flux(i, 1) << " " << rate(i) << " " << trans2 << " " << flux2(i, 0) << " "
              << flux2(i, 1) << " " << spectrumSum(i) << " " << spectrumBnd[i] << " " << spectrumBnd[i + 1] << " "
              << spectrumSum2(i) << " " << F4 * PI * POW2(k) << " "
              << F4B3 * PI * (POW3(spectrumBnd[i + 1]) - POW3(spectrumBnd[i])) << std::endl;
        }
        for(MInt i = 0; i < noBins / 2; i++) {
          const MInt i0 = 2 * i;
          const MInt i1 = 2 * i + 1;
          const MFloat k = F1B2 * (spectrumBnd[i1 + 1] + spectrumBnd[i0]);
          MFloat vol0 =
              F4B3 * PI * (POW3(spectrumBnd[i0 + 1]) - POW3(spectrumBnd[i0])) / (spectrumBnd[i0 + 1] - spectrumBnd[i0]);
          MFloat vol1 =
              F4B3 * PI * (POW3(spectrumBnd[i1 + 1]) - POW3(spectrumBnd[i1])) / (spectrumBnd[i1 + 1] - spectrumBnd[i1]);
          MFloat vol =
              F4B3 * PI * (POW3(spectrumBnd[i1 + 1]) - POW3(spectrumBnd[i0])) / (spectrumBnd[i1 + 1] - spectrumBnd[i0]);
          MFloat f0 = spectrumSum(i0) * vol / (vol0 * (spectrumSum(i0) + spectrumSum(i1)));
          MFloat f1 = spectrumSum(i1) * vol / (vol1 * (spectrumSum(i0) + spectrumSum(i1)));
 
          ofl2 << k / k0 << " " << f1 * spectrum(i1) + f0 * spectrum(i0) << " "
               << F2 * POW2(k) * viscosity * (f1 * spectrum(i1) + f0 * spectrum(i0)) << " "
               << f1 * transfer(i1) + f0 * transfer(i0) << " " << f1 * coupling(i1) + f0 * coupling(i0) << " "
               << f1 * flux(i1, 0) + f0 * flux(i0, 0) << " " << f1 * flux(i1, 1) + f0 * flux(i0, 1) << " "
               << f1 * rate(i1) + f0 * rate(i0) << " " << f1 * transfer2(i1) + f0 * transfer2(i0) << " "
               << f1 * flux2(i1, 0) + f0 * flux2(i0, 0) << " " << f1 * flux2(i1, 1) + f0 * flux2(i0, 1) << " "
               << spectrumSum(i0) + spectrumSum(i1) << " " << spectrumBnd[i0] << " " << spectrumBnd[i1 + 1] << " "
               << spectrumSum2(i0) + spectrumSum2(i1) << " " << F4 * PI * POW2(k) << " "
               << F4B3 * PI * (POW3(spectrumBnd[i1 + 1]) - POW3(spectrumBnd[i0])) << std::endl;
        }
        ofl.close();
        ofl2.close();
      }
    }
 
    // cleanup
    // fftw_cleanup(); problems occur when calling clean up function!!!
  }
 
  const MFloat time3 = MPI_Wtime();
 
  if(domainId == 0)
    std::cerr << "fft time " << time1 - time0 << " " << time2 - time1 << " " << time3 - time2 << std::endl;
 
  m_log << "fft finished" << std::endl;
 
  DEBUG("computeEnergySpectrum return", MAIA_DEBUG_TRACE_OUT);
}

computeRotationMatrix()

void maia::math::computeRotationMatrix ( MFloatScratchSpace &  R, MFloat *  q  )

inline

Author

Lennart Schneiders

Definition at line 533 of file maiamath.h.

                                                                    {
  R(0, 0) = q[0] * q[0] + q[1] * q[1] - q[2] * q[2] - q[3] * q[3];
  R(1, 1) = q[0] * q[0] - q[1] * q[1] + q[2] * q[2] - q[3] * q[3];
  R(2, 2) = q[0] * q[0] - q[1] * q[1] - q[2] * q[2] + q[3] * q[3];
 
  R(0, 1) = F2 * (q[1] * q[2] + q[0] * q[3]);
  R(0, 2) = F2 * (q[1] * q[3] - q[0] * q[2]);
  R(1, 0) = F2 * (q[1] * q[2] - q[0] * q[3]);
  R(1, 2) = F2 * (q[2] * q[3] + q[0] * q[1]);
  R(2, 0) = F2 * (q[1] * q[3] + q[0] * q[2]);
  R(2, 1) = F2 * (q[2] * q[3] - q[0] * q[1]);
}

ConvertToEigenMatrix()

template<typename T , std::size_t N>

Eigen::MatrixXd maia::math::ConvertToEigenMatrix

(

std::array< std::array< T, N >, N >

data

)

Definition at line 66 of file maiamath.cpp.

                                                                     {
  Eigen::Matrix<T, N, N> m;
  for(std::size_t i = 0; i < N; ++i) {
    m.row(i) = Eigen::VectorXd::Map(&data[i][0], N);
  }
  return m;
}

crankAngle()

MFloat maia::math::crankAngle ( const MFloat  time, const MFloat  Strouhal, const MFloat  offset, const MInt  mode  )

inline

Author

Tim Wegmann

Definition at line 506 of file maiamath.h.

                                                                                                         {
  const MFloat mu2 = Strouhal * F2 * PI;
  MFloat cad = mu2 * time;
 
  if(mode == 0) {
    const MInt maxNoCycles = 20;
    for(MInt cycle = maxNoCycles; cycle > 0; cycle--) {
      if(cad >= 4 * PI * cycle) {
        cad = cad - 4 * PI * cycle;
      }
    }
    cad = cad * 180 / PI;
 
    cad = cad + offset;
 
  } else {
    ASSERT(mode == 1, "Incorrect mode!");
    cad = cad + offset * PI / 180;
  }
 
  return cad;
}

cross() [1/5]

template<typename T >

T maia::math::cross ( const std::array< T, 2 > &  u, const std::array< T, 2 > &  v  )

inline

Definition at line 120 of file maiamath.h.

                                                                 {
  return u[0] * v[1] - u[1] * v[0];
}

cross() [2/5]

template<typename T >

std::array< T, 3 > maia::math::cross ( const std::array< T, 3 > &  u, const std::array< T, 3 > &  v  )

inline

Cross product C = U x V [3D]

Template Parameters

T

Type

Parameters

[in]	u	U
[in]	v	V

Returns

New copy of an array

Definition at line 113 of file maiamath.h.

                                                                              {
  std::array<T, 3> result;
  cross(&u[0], &v[0], &result[0]);
  return result;
}

cross() [3/5]

template<typename T >

void maia::math::cross ( const T *const  u, const T *const  v, T *const  c  )

inline

Cross product C = U x V [3D]

Template Parameters

T

Type

Parameters

[in]	u	U
[in]	v	V
[out]	c	C

Author

Sven Berger

Definition at line 101 of file maiamath.h.

                                                                  {
  c[0] = u[1] * v[2] - v[1] * u[2];
  c[1] = u[2] * v[0] - v[2] * u[0];
  c[2] = u[0] * v[1] - v[0] * u[1];
}

cross() [4/5]

template<typename T >

T maia::math::cross ( const T(&)  u[2], const T(&)  v[2]  )

inline

Definition at line 137 of file maiamath.h.

                                                 {
  return u[0] * v[1] - u[1] * v[0];
}

cross() [5/5]

template<typename T >

T * maia::math::cross ( const T(&)  u[3], const T(&)  v[3]  )

inline

Cross product C = U x V [3D]

Template Parameters

T

Type

Parameters

[in]	u	U
[in]	v	V

Returns

New copy of an array

Definition at line 130 of file maiamath.h.

                                                  {
  T result[3]{};
  cross(&u[0], &v[0], &result[0]);
  return result;
}

deltaFun()

MFloat maia::math::deltaFun ( const MFloat  r, const MFloat  r0, const MFloat  r1  )

inline

Definition at line 885 of file maiamath.h.

                                                                         {
  MFloat R = mMin(F1, mMax(F0, ((r - r0) / (r1 - r0))));
#if MAIA_TRANSITION_FUNCTION == 0
  return (r < r1) ? F0 : F1;
#elif MAIA_TRANSITION_FUNCTION == 1
  return R;
#elif MAIA_TRANSITION_FUNCTION == 2
  return POW2(R) * (3.0 - 2.0 * R);
#elif MAIA_TRANSITION_FUNCTION == 3
  return POW3(R) * (10.0 + R * (6.0 * R - 15.0)); // may cause stability issues
#elif MAIA_TRANSITION_FUNCTION == 4
  return POW3(R) * (6.0 + R * (2.0 * R - 7.0)); // approximates 0.5*(1.0-cos(PI*pow(R,1.25)));
#elif MAIA_TRANSITION_FUNCTION == 5
  return F1B2 * (F1 + sin(PI * (F3B2 + R)));
#endif
}

determinant() [1/2]

template<typename T , std::size_t N>

MFloat maia::math::determinant

(

std::array< std::array< T, N >, N > &

m

)

Determinant of a matrix

Template Parameters

T	Type
N	Dimension of the square matrix

Parameters

m	Matrix (two dimensional array containing the matrix elements)

Returns

Determinant of the matrix

Definition at line 145 of file maiamath.cpp.

                                                   {
  Eigen::Matrix<T, N, N> matrix(&m[0][0]);
  return matrix.determinant();
}

determinant() [2/2]

template<typename T , std::size_t N>

MFloat maia::math::determinant

(

std::array< T, N > &

m

)

Determinant of a matrix

Template Parameters

T	Type
N	Size of the matrix (rows * columns)

Parameters

m	Matrix (one dimensional array containing the matrix elements)

Returns

Determinant of the matrix

Definition at line 130 of file maiamath.cpp.

                                      {
  constexpr MInt dim = N == 2 ? 1 : 2;
  Eigen::Matrix<T, dim, dim> matrix(m.data());
  return matrix.determinant();
}

determinant< MFloat, 2ul >() [1/2]

template MFloat maia::math::determinant< MFloat, 2ul >

(

std::array< MFloat, 2ul > &

)

determinant< MFloat, 2ul >() [2/2]

template MFloat maia::math::determinant< MFloat, 2ul >

(

std::array< std::array< MFloat, 2ul >, 2ul > &

)

determinant< MFloat, 3ul >() [1/2]

template MFloat maia::math::determinant< MFloat, 3ul >

(

std::array< MFloat, 3ul > &

)

determinant< MFloat, 3ul >() [2/2]

template MFloat maia::math::determinant< MFloat, 3ul >

(

std::array< std::array< MFloat, 3ul >, 3ul > &

)

determinant< MFloat, 4ul >() [1/2]

template MFloat maia::math::determinant< MFloat, 4ul >

(

std::array< MFloat, 4ul > &

)

determinant< MFloat, 4ul >() [2/2]

template MFloat maia::math::determinant< MFloat, 4ul >

(

std::array< std::array< MFloat, 4ul >, 4ul > &

)

distance() [1/7]

template<MInt nDim>

MFloat maia::math::distance ( const MFloat *  a, const MFloat *  b  )

inline

Definition at line 249 of file maiamath.h.

                                                         {
  IF_CONSTEXPR(nDim == 2) { return sqrt(POW2(a[0] - b[0]) + POW2(a[1] - b[1])); }
  IF_CONSTEXPR(nDim == 3) { return sqrt(POW2(a[0] - b[0]) + POW2(a[1] - b[1]) + POW2(a[2] - b[2])); }
  return -NAN;
}

distance() [2/7]

MFloat maia::math::distance ( const MFloat *const  a, const std::array< MFloat, 2 >  b  )

inline

Definition at line 263 of file maiamath.h.

                                                                           {
  return sqrt(POW2(a[0] - b[0]) + POW2(a[1] - b[1]));
}

distance() [3/7]

MFloat maia::math::distance ( const MFloat *const  a, const std::array< MFloat, 3 >  b  )

inline

Definition at line 267 of file maiamath.h.

                                                                           {
  return sqrt(POW2(a[0] - b[0]) + POW2(a[1] - b[1]) + POW2(a[2] - b[2]));
}

distance() [4/7]

MFloat maia::math::distance ( const std::array< MFloat, 2 >  a, const MFloat *const  b  )

inline

Definition at line 271 of file maiamath.h.

                                                                           {
  return sqrt(POW2(a[0] - b[0]) + POW2(a[1] - b[1]));
}

distance() [5/7]

MFloat maia::math::distance ( const std::array< MFloat, 2 >  a, const std::array< MFloat, 2 >  b  )

inline

Definition at line 255 of file maiamath.h.

                                                                                 {
  return sqrt(POW2(a[0] - b[0]) + POW2(a[1] - b[1]));
}

distance() [6/7]

MFloat maia::math::distance ( const std::array< MFloat, 3 >  a, const std::array< MFloat, 3 >  b  )

inline

Definition at line 259 of file maiamath.h.

                                                                                 {
  return sqrt(POW2(a[0] - b[0]) + POW2(a[1] - b[1]) + POW2(a[2] - b[2]));
}

distance() [7/7]

MFloat maia::math::distance ( std::array< MFloat, 3 >  a, const MFloat *const  b  )

inline

Definition at line 275 of file maiamath.h.

                                                                     {
  return sqrt(POW2(a[0] - b[0]) + POW2(a[1] - b[1]) + POW2(a[2] - b[2]));
}

frobeniusMatrixNorm()

MFloat maia::math::frobeniusMatrixNorm ( MFloatScratchSpace &  m, MInt  dim1, MInt  dim2  )

inline

Definition at line 339 of file maiamath.h.

                                                                               {
  return sqrt(frobeniusMatrixNormSquared(m, dim1, dim2));
}

frobeniusMatrixNormSquared()

MFloat maia::math::frobeniusMatrixNormSquared ( MFloatScratchSpace &  m, MInt  dim1, MInt  dim2  )

inline

Definition at line 328 of file maiamath.h.

                                                                                      {
  MFloat ret_val = 0.0;
  for(MInt i = 0; i < dim1; i++) {
    for(MInt j = 0; j < dim2; j++) {
      ret_val += m(i, j) * m(i, j);
    }
  }
 
  return ret_val;
}

getAngle()

MFloat maia::math::getAngle ( MFloat  y, MFloat  z  )

inline

Definition at line 925 of file maiamath.h.

                                           {
  MFloat angle = atan(z / y);
  if(y < 0) {
    if(z >= 0) {
      angle += PI;
    } else {
      angle -= PI;
    }
  }
  return angle;
}

getFourierCoefficient()

std::complex< MFloat > * maia::math::getFourierCoefficient ( MFloat *  k, MFloat  k0, const MFloat  Ma  )

inline

brief Generates a single std::complex coefficient of Fourier series

for a given wavenumber k, and a certain energy spectrum (see Appendix of Orszag, 1969)

Definition at line 548 of file maiafftw.h.

                                                                                        {
  TRACE();
 
  MFloat r[6], s[6], kAbs, energy;
  std::complex<MFloat> uHat, vHat, wHat;
  std::complex<MFloat>* fourierCoefficient;
 
  kAbs = sqrt(k[0] * k[0] + k[1] * k[1] + k[2] * k[2]);
 
  // the zero-frequency component is always set to zero, so there is no offset
  if(approx(kAbs, 0.0, MFloatEps)) {
    fourierCoefficient = new std::complex<MFloat>[3];
 
    fourierCoefficient[0] = std::complex<MFloat>(0, 0);
    fourierCoefficient[1] = std::complex<MFloat>(0, 0);
    fourierCoefficient[2] = std::complex<MFloat>(0, 0);
 
 
    return fourierCoefficient;
 
  } else {
    // energy = (kAbs/k0)*(kAbs/k0)*(kAbs/k0)*(kAbs/k0) * exp(-2.0*(kAbs/k0)*(kAbs/k0));
    // energy = pow(kAbs/k0,8.0) * exp(-4.0*(kAbs/k0)*(kAbs/k0)); // set spectral distribution
    energy = pow(kAbs / k0, 4.0) * exp(-2.0 * (kAbs / k0) * (kAbs / k0)); // set spectral distribution
    energy *= exp(2.0) * 0.499 * (Ma * LBCS)
              * (Ma * LBCS); // set maximal fluctuation amplitude to 20% of the freestream velocity (for 128^3: 0.88)
 
    // determine Fourier coefficients:
    // r and s are Independant random vector fields with independant
    // components (zero mean and rms according to energy spectrum).
    // Each vector has three components for k and another three for -k.
 
    for(MInt i = 0; i < 6; i++) {
      r[i] = randnormal(0.0, PI * sqrt(energy) / (SQRT2 * kAbs));
      // r[i] = randNumGen.randNorm(0.0, PI*sqrt(energy)/(SQRT2*kAbs));
      s[i] = randnormal(0.0, PI * sqrt(energy) / (SQRT2 * kAbs));
      // s[i] = randNumGen.randNorm(0.0, PI*sqrt(energy)/(SQRT2*kAbs));
    }
 
    uHat = (1.0 - k[0] * k[0] / (kAbs * kAbs)) * std::complex<MFloat>(r[0] + r[3], s[0] - s[3])
           - k[0] * k[1] / (kAbs * kAbs) * std::complex<MFloat>(r[1] + r[4], s[1] - s[4])
           - k[0] * k[2] / (kAbs * kAbs) * std::complex<MFloat>(r[2] + r[5], s[2] - s[5]);
 
 
    vHat = -k[1] * k[0] / (kAbs * kAbs) * std::complex<MFloat>(r[0] + r[3], s[0] - s[3])
           + (1.0 - k[1] * k[1] / (kAbs * kAbs)) * std::complex<MFloat>(r[1] + r[4], s[1] - s[4])
           - k[1] * k[2] / (kAbs * kAbs) * std::complex<MFloat>(r[2] + r[5], s[2] - s[5]);
 
 
    wHat = -k[2] * k[0] / (kAbs * kAbs) * std::complex<MFloat>(r[0] + r[3], s[0] - s[3])
           - k[2] * k[1] / (kAbs * kAbs) * std::complex<MFloat>(r[1] + r[4], s[1] - s[4])
           + (1.0 - k[2] * k[2] / (kAbs * kAbs)) * std::complex<MFloat>(r[2] + r[5], s[2] - s[5]);
 
    fourierCoefficient = new std::complex<MFloat>[3];
 
    // fourierCoefficient[0] = std::complex<MFloat>(sqrt(2*energy)/SQRT2,sqrt(2*energy)/SQRT2);// uHat;
    // fourierCoefficient[1] = std::complex<MFloat>(sqrt(2*energy)/SQRT2,sqrt(2*energy)/SQRT2);//vHat;
    // fourierCoefficient[2] = std::complex<MFloat>(sqrt(2*energy)/SQRT2,sqrt(2*energy)/SQRT2);//wHat;
 
    fourierCoefficient[0] = uHat;
    fourierCoefficient[1] = vHat;
    fourierCoefficient[2] = wHat;
    return fourierCoefficient;
  }
}

getFourierCoefficient2()

std::complex< MFloat > * maia::math::getFourierCoefficient2 ( MFloat *  k, MFloat  k0  )

inline

brief Generates a single std::complex coefficient of Fourier series

for a given wavenumber k, and a certain energy spectrum

Definition at line 620 of file maiafftw.h.

                                                                        {
  static std::mt19937 randNumGen;
  static std::uniform_real_distribution<> distrib{0.0, 1.0};
 
  MFloat r[3], kAbs, energy;
  std::complex<MFloat> uHat, vHat, wHat;
  std::complex<MFloat>* fourierCoefficient;
 
  kAbs = sqrt(k[0] * k[0] + k[1] * k[1] + k[2] * k[2]);
 
  // the zero-frequency component is always set to zero, so there is no offset
  if(approx(kAbs, 0.0, MFloatEps)) {
    fourierCoefficient = new std::complex<MFloat>[3];
 
    fourierCoefficient[0] = std::complex<MFloat>(0, 0);
    fourierCoefficient[1] = std::complex<MFloat>(0, 0);
    fourierCoefficient[2] = std::complex<MFloat>(0, 0);
 
 
    return fourierCoefficient;
 
  } else {
    // energy = pow(kAbs/k0,4.0) * exp(-2.0*(kAbs/k0)*(kAbs/k0));
    // energy = ( m_Ma*LBCS * m_Ma*LBCS ) * pow(kAbs/k0,4.0) * exp(-2.0*(kAbs/k0)*(kAbs/k0));
    // energy = 0.135335;
    energy = (kAbs / k0) * (kAbs / k0) * (kAbs / k0) * (kAbs / k0) * exp(-2.0 * (kAbs / k0) * (kAbs / k0));
 
    // determine Fourier coefficients:
    for(MInt i = 0; i < 3; i++) {
      r[i] = 2.0 * PI * distrib(randNumGen);
    }
 
    uHat = std::complex<MFloat>(cos(r[0]) * energy, sin(r[0]) * energy);
    vHat = std::complex<MFloat>(cos(r[1]) * energy, sin(r[1]) * energy);
    wHat = std::complex<MFloat>(cos(r[2]) * energy, sin(r[2]) * energy);
 
    fourierCoefficient = new std::complex<MFloat>[3];
 
    fourierCoefficient[0] = uHat;
    fourierCoefficient[1] = vHat;
    fourierCoefficient[2] = wHat;
 
 
    return fourierCoefficient;
  }
}

getGlobalPosFFTW()

MInt maia::math::getGlobalPosFFTW ( MInt  i0, MInt  i1, MInt  i2, MInt  ny, MInt  nz  )

inline

Definition at line 59 of file maiafftw.h.

{ return i2 + nz * (i1 + ny * i0); }

getSector()

MFloat maia::math::getSector ( MFloat  y, MFloat  z, MFloat  azimuthalAngle  )

inline

Definition at line 912 of file maiamath.h.

                                                                   {
  MFloat angle = atan(z / y);
  if(y < 0) {
    if(z >= 0) {
      angle += PI;
    } else {
      angle -= PI;
    }
  }
  angle += PI;
  return (angle - std::fmod(angle, azimuthalAngle / 180.0 * PI)) / (azimuthalAngle / 180.0 * PI);
}

initFft()

void maia::math::initFft ( fftw_complex *  uPhysField, fftw_complex *  vPhysField, fftw_complex *  wPhysField, MInt  lx, MInt  ly, MInt  lz, MInt  noPeakModes, const MFloat  Ma  )

inline

Author

george

Date

09.12.2011

(u,v,w) PhysField: complex velocity in physical space

Parameters

uPhysField	pointer to a fftw_complex for u
vPhysField	pointer to a fftw_complex for v
wPhysField	pointer to a fftw_complex for w

ATTENTION: This version is quite old and should be thoroughly before use!

Definition at line 892 of file maiafftw.h.

                                     {
  TRACE();
 
  MFloat waveVector[3], k0;
 
  std::complex<MFloat>* fourierCoefficient;
 
  fftw_complex *uHatField, *vHatField, *wHatField;
 
  fftw_plan planU, planV, planW;
 
  if(lx % 2 != 0 || ly % 2 != 0 || lz % 2 != 0) {
    std::stringstream errorMessage;
    errorMessage << " FFTInit: domainsize must NOT be an odd number! " << std::endl;
    mTerm(1, AT_, errorMessage.str());
  }
 
  m_log << " --- initializing FFTW --- " << std::endl;
  m_log << " domain size = " << lx << "x" << ly << "x" << lz << std::endl;
 
  // allocate field of Fourier coefficients (is deleted after the transformation)
  uHatField = (fftw_complex*)fftw_malloc(lx * ly * lz * sizeof(fftw_complex));
  vHatField = (fftw_complex*)fftw_malloc(lx * ly * lz * sizeof(fftw_complex));
  wHatField = (fftw_complex*)fftw_malloc(lx * ly * lz * sizeof(fftw_complex));
 
  // create plans for FFTW
  planU = fftw_plan_dft_3d(lx, ly, lz, uHatField, uPhysField, FFTW_BACKWARD, FFTW_MEASURE);
  planV = fftw_plan_dft_3d(lx, ly, lz, vHatField, vPhysField, FFTW_BACKWARD, FFTW_MEASURE);
  planW = fftw_plan_dft_3d(lx, ly, lz, wHatField, wPhysField, FFTW_BACKWARD, FFTW_MEASURE);
 
  // reset coefficients and velocities
  for(MInt p = 0; p < lx; p++) {
    for(MInt q = 0; q < ly; q++) {
      for(MInt r = 0; r < lz; r++) {
        uHatField[r + lz * (q + ly * p)][0] = 0.0;
        uHatField[r + lz * (q + ly * p)][1] = 0.0;
        uPhysField[r + lz * (q + ly * p)][0] = 0.0;
        uPhysField[r + lz * (q + ly * p)][1] = 0.0;
 
        vHatField[r + lz * (q + ly * p)][0] = 0.0;
        vHatField[r + lz * (q + ly * p)][1] = 0.0;
        vPhysField[r + lz * (q + ly * p)][0] = 0.0;
        vPhysField[r + lz * (q + ly * p)][1] = 0.0;
 
        wHatField[r + lz * (q + ly * p)][0] = 0.0;
        wHatField[r + lz * (q + ly * p)][1] = 0.0;
        wPhysField[r + lz * (q + ly * p)][0] = 0.0;
        wPhysField[r + lz * (q + ly * p)][1] = 0.0;
      }
    }
  }
 
  // Set coefficients:
  // FFTW stores the coefficients for positive wavenumbers in the first half of the array,
  // and those for negative wavenumbers in reverse order in the second half.
  // [0, 1, ... , N/2-1, N/2, ... , N-1]
  //  - the entry for zero-wavenumber is at position 0
  //  - the k-th entry and the (N-k)th entry correspond to wavenumbers with opposite sign
  //  - the entry at position N/2 corresponds to the Nyquist wavenumber and appears only once
 
  // peak wave number of energy spectrum
  k0 = 2.0 * PI / (lx / noPeakModes);
 
  for(MInt p = 0; p <= lx / 2; p++) {
    for(MInt q = 0; q <= ly / 2; q++) {
      for(MInt r = 0; r <= lz / 2; r++) {
        // wave-vector: k(p,q,r) = (2 \pi p / lx, 2 \pi q / ly, 2 \pi r / lz)
        waveVector[0] = (p)*2.0 * PI / lx;
        waveVector[1] = (q)*2.0 * PI / ly;
        waveVector[2] = (r)*2.0 * PI / lz;
 
        fourierCoefficient = getFourierCoefficient(waveVector, k0, Ma);
 
        // 1. Positive frequencies:
        uHatField[r + lz * (q + ly * p)][0] = std::real(fourierCoefficient[0]);
        uHatField[r + lz * (q + ly * p)][1] = std::imag(fourierCoefficient[0]);
 
        vHatField[r + lz * (q + ly * p)][0] = std::real(fourierCoefficient[1]);
        vHatField[r + lz * (q + ly * p)][1] = std::imag(fourierCoefficient[1]);
 
        wHatField[r + lz * (q + ly * p)][0] = std::real(fourierCoefficient[2]);
        wHatField[r + lz * (q + ly * p)][1] = std::imag(fourierCoefficient[2]);
 
        // 2. Negative frequencies:
        if(p > 1 && q > 1 && r > 1) {
          if(p < lx / 2 && q < ly / 2 && r < lz / 2) {
            // since the physical velocity field is real, the coefficients for negative frequencies
            // are the std::complex conjugate of those for positive frequencies
            uHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][0] = uHatField[r + lz * (q + ly * p)][0];
            uHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][1] = -uHatField[r + lz * (q + ly * p)][1];
 
            vHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][0] = vHatField[r + lz * (q + ly * p)][0];
            vHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][1] = -vHatField[r + lz * (q + ly * p)][1];
 
            wHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][0] = wHatField[r + lz * (q + ly * p)][0];
            wHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][1] = -wHatField[r + lz * (q + ly * p)][1];
          }
        }
      }
    }
  }
 
  // // test:
  // // cos(x) in a cube with 32x32x32 cells
  // uHatField[(0+lz*(0+ly*2)=2048][0]=90.5;
  // uHatField[(0+lz*(0+ly*30)=30720][0]=90.5;
 
  // Do Fourier transform (backward, see plan definition)
  // Definition in one dimension:
  // u(x) = \sum_{j=0}^{lx-1} \hat{u}_j exp(i 2 \pi j x / lx)
  fftw_execute(planU);
  fftw_execute(planV);
  fftw_execute(planW);
 
  // normalize (this preserves the norm of the basis functions)
  for(MInt p = 0; p < lx; p++) {
    for(MInt q = 0; q < ly; q++) {
      for(MInt r = 0; r < lz; r++) {
        uPhysField[r + lz * (q + ly * p)][0] /= sqrt(MFloat(lx * ly * lz));
        vPhysField[r + lz * (q + ly * p)][0] /= sqrt(MFloat(lx * ly * lz));
        wPhysField[r + lz * (q + ly * p)][0] /= sqrt(MFloat(lx * ly * lz));
 
        uPhysField[r + lz * (q + ly * p)][1] /= sqrt(MFloat(lx * ly * lz));
        vPhysField[r + lz * (q + ly * p)][1] /= sqrt(MFloat(lx * ly * lz));
        wPhysField[r + lz * (q + ly * p)][1] /= sqrt(MFloat(lx * ly * lz));
      }
    }
  }
 
  fftw_destroy_plan(planU);
  fftw_destroy_plan(planV);
  fftw_destroy_plan(planW);
  fftw_free(uHatField);
  fftw_free(vHatField);
  fftw_free(wHatField);
}

initFftFilter()

void maia::math::initFftFilter ( MFloat *  uPhysFieldCoarse, MFloat *  vPhysFieldCoarse, MFloat *  wPhysFieldCoarse, MInt  lx, MInt  ly, MInt  lz, MInt  lxCoarse, MInt  lyCoarse, MInt  lzCoarse, MInt  noPeakModes, const MFloat  Ma  )

inline

Author

george

Date

24.01.2011

(u,v,w)PhysFieldCoarse: real velocity in physical space

Parameters

uPhysFieldCoarse	pointer to a MFloat for u
vPhysFieldCoarse	pointer to a MFloat for v
wPhysFieldCoarse	pointer to a MFloat for w

ATTENTION: This version is quite old and should be thoroughly before use!

Definition at line 681 of file maiafftw.h.

                                           {
  TRACE();
 
 
  MFloat waveVector[3], k0;
 
  std::complex<MFloat>* fourierCoefficient;
 
  fftw_complex *uHatField, *vHatField, *wHatField;
  fftw_complex *uPhysField, *vPhysField, *wPhysField;
 
  fftw_plan planU, planV, planW;
 
  MInt xRatio = lx / lxCoarse;
  MInt yRatio = ly / lyCoarse;
  MInt zRatio = lz / lzCoarse;
 
  if(lx % 2 != 0 || ly % 2 != 0 || lz % 2 != 0) {
    mTerm(1, AT_, " FFTInit: domainsize must NOT be an odd number! ");
  }
 
  m_log << " --- initializing FFTW --- " << std::endl;
  m_log << " domain size = " << lx << "x" << ly << "x" << lz << std::endl;
 
  // allocate field of velocities (is deleted at the end of the method)
  uPhysField = (fftw_complex*)fftw_malloc(lx * ly * lz * sizeof(fftw_complex));
  vPhysField = (fftw_complex*)fftw_malloc(lx * ly * lz * sizeof(fftw_complex));
  wPhysField = (fftw_complex*)fftw_malloc(lx * ly * lz * sizeof(fftw_complex));
 
  // allocate field of Fourier coefficients (is deleted after the transformation)
  uHatField = (fftw_complex*)fftw_malloc(lx * ly * lz * sizeof(fftw_complex));
  vHatField = (fftw_complex*)fftw_malloc(lx * ly * lz * sizeof(fftw_complex));
  wHatField = (fftw_complex*)fftw_malloc(lx * ly * lz * sizeof(fftw_complex));
 
  // create plans for FFTW
  planU = fftw_plan_dft_3d(lx, ly, lz, uHatField, uPhysField, FFTW_BACKWARD, FFTW_MEASURE);
  planV = fftw_plan_dft_3d(lx, ly, lz, vHatField, vPhysField, FFTW_BACKWARD, FFTW_MEASURE);
  planW = fftw_plan_dft_3d(lx, ly, lz, wHatField, wPhysField, FFTW_BACKWARD, FFTW_MEASURE);
 
  // reset coefficients and velocities
  for(MInt p = 0; p < lx; p++) {
    for(MInt q = 0; q < ly; q++) {
      for(MInt r = 0; r < lz; r++) {
        uHatField[r + lz * (q + ly * p)][0] = 0.0;
        uHatField[r + lz * (q + ly * p)][1] = 0.0;
        uPhysField[r + lz * (q + ly * p)][0] = 0.0;
        uPhysField[r + lz * (q + ly * p)][1] = 0.0;
 
        vHatField[r + lz * (q + ly * p)][0] = 0.0;
        vHatField[r + lz * (q + ly * p)][1] = 0.0;
        vPhysField[r + lz * (q + ly * p)][0] = 0.0;
        vPhysField[r + lz * (q + ly * p)][1] = 0.0;
 
        wHatField[r + lz * (q + ly * p)][0] = 0.0;
        wHatField[r + lz * (q + ly * p)][1] = 0.0;
        wPhysField[r + lz * (q + ly * p)][0] = 0.0;
        wPhysField[r + lz * (q + ly * p)][1] = 0.0;
      }
    }
  }
 
  // Set coefficients:
  // FFTW stores the coefficients for positive wavenumbers in the first half of the array,
  // and those for negative wavenumbers in reverse order in the second half.
  // [0, 1, ... , N/2-1, N/2, ... , N-1]
  //  - the entry for zero-wavenumber is at position 0
  //  - the k-th entry and the (N-k)th entry correspond to wavenumbers with opposite sign
  //  - the entry at position N/2 corresponds to the Nyquist wavenumber and appears only once
 
  // peak wave number of energy spectrum
  k0 = 2.0 * PI / (lx / noPeakModes);
 
  for(MInt p = 0; p <= lx / 2; p++) {
    for(MInt q = 0; q <= ly / 2; q++) {
      for(MInt r = 0; r <= lz / 2; r++) {
        // wave-vector: k(p,q,r) = (2 \pi p / lx, 2 \pi q / ly, 2 \pi r / lz)
        waveVector[0] = (p)*2.0 * PI / lx;
        waveVector[1] = (q)*2.0 * PI / ly;
        waveVector[2] = (r)*2.0 * PI / lz;
 
        fourierCoefficient = getFourierCoefficient(waveVector, k0, Ma);
 
        // 1. Positive frequencies:
        uHatField[r + lz * (q + ly * p)][0] = std::real(fourierCoefficient[0]);
        uHatField[r + lz * (q + ly * p)][1] = std::imag(fourierCoefficient[0]);
 
        vHatField[r + lz * (q + ly * p)][0] = std::real(fourierCoefficient[1]);
        vHatField[r + lz * (q + ly * p)][1] = std::imag(fourierCoefficient[1]);
 
        wHatField[r + lz * (q + ly * p)][0] = std::real(fourierCoefficient[2]);
        wHatField[r + lz * (q + ly * p)][1] = std::imag(fourierCoefficient[2]);
 
        // 2. Negative frequencies:
        if(p > 1 && q > 1 && r > 1) {
          if(p < lx / 2 && q < ly / 2 && r < lz / 2) {
            // since the physical velocity field is real, the coefficients for negative frequencies
            // are the std::complex conjugate of those for positive frequencies
            uHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][0] = uHatField[r + lz * (q + ly * p)][0];
            uHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][1] = -uHatField[r + lz * (q + ly * p)][1];
 
            vHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][0] = vHatField[r + lz * (q + ly * p)][0];
            vHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][1] = -vHatField[r + lz * (q + ly * p)][1];
 
            wHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][0] = wHatField[r + lz * (q + ly * p)][0];
            wHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][1] = -wHatField[r + lz * (q + ly * p)][1];
          }
        }
      }
    }
  }
 
  // Do Fourier transform (backward, see plan definition)
  // Definition in one dimension:
  // u(x) = \sum_{j=0}^{lx-1} \hat{u}_j exp(i 2 \pi j x / lx)
  fftw_execute(planU);
  fftw_execute(planV);
  fftw_execute(planW);
 
  // normalize (this preserves the norm of the basis functions)
  for(MInt p = 0; p < lx; p++) {
    for(MInt q = 0; q < ly; q++) {
      for(MInt r = 0; r < lz; r++) {
        uPhysField[r + lz * (q + ly * p)][0] /= sqrt(MFloat(lx * ly * lz));
        vPhysField[r + lz * (q + ly * p)][0] /= sqrt(MFloat(lx * ly * lz));
        wPhysField[r + lz * (q + ly * p)][0] /= sqrt(MFloat(lx * ly * lz));
 
        uPhysField[r + lz * (q + ly * p)][1] /= sqrt(MFloat(lx * ly * lz));
        vPhysField[r + lz * (q + ly * p)][1] /= sqrt(MFloat(lx * ly * lz));
        wPhysField[r + lz * (q + ly * p)][1] /= sqrt(MFloat(lx * ly * lz));
      }
    }
  }
 
  // // test:
  // // cos(x) in a cube with 32x32x32 cells
  // uHatField[(0+lz*(0+ly*2)=2048][0]=90.5;
  // uHatField[(0+lz*(0+ly*30)=30720][0]=90.5;
 
  // reset coarse fields
  for(MInt p = 0; p < lxCoarse; p++) {
    for(MInt q = 0; q < lyCoarse; q++) {
      for(MInt r = 0; r < lzCoarse; r++) {
        uPhysFieldCoarse[r + lzCoarse * (q + lyCoarse * p)] = F0;
        vPhysFieldCoarse[r + lzCoarse * (q + lyCoarse * p)] = F0;
        wPhysFieldCoarse[r + lzCoarse * (q + lyCoarse * p)] = F0;
      }
    }
  }
 
  // transfer physField to physFieldCoarse via arithmetic averaging
  for(MInt p = 0; p < lxCoarse; p++) {
    for(MInt q = 0; q < lyCoarse; q++) {
      for(MInt r = 0; r < lzCoarse; r++) {
        for(MInt i = p * xRatio; i < p * zRatio + xRatio; i++) {
          for(MInt j = q * yRatio; j < q * zRatio + yRatio; j++) {
            for(MInt k = r * zRatio; k < r * zRatio + zRatio; k++) {
              uPhysFieldCoarse[r + lzCoarse * (q + lyCoarse * p)] += uPhysField[k + lz * (j + ly * i)][0];
              vPhysFieldCoarse[r + lzCoarse * (q + lyCoarse * p)] += vPhysField[k + lz * (j + ly * i)][0];
              wPhysFieldCoarse[r + lzCoarse * (q + lyCoarse * p)] += wPhysField[k + lz * (j + ly * i)][0];
            }
          }
        }
      }
    }
  }
  for(MInt p = 0; p < lxCoarse; p++) {
    for(MInt q = 0; q < lyCoarse; q++) {
      for(MInt r = 0; r < lzCoarse; r++) {
        uPhysFieldCoarse[r + lzCoarse * (q + lyCoarse * p)] /= (xRatio * yRatio * zRatio);
        vPhysFieldCoarse[r + lzCoarse * (q + lyCoarse * p)] /= (xRatio * yRatio * zRatio);
        wPhysFieldCoarse[r + lzCoarse * (q + lyCoarse * p)] /= (xRatio * yRatio * zRatio);
      }
    }
  }
 
  fftw_destroy_plan(planU);
  fftw_destroy_plan(planV);
  fftw_destroy_plan(planW);
 
  fftw_free(uHatField);
  fftw_free(vHatField);
  fftw_free(wHatField);
 
  fftw_free(uPhysField);
  fftw_free(vPhysField);
  fftw_free(wPhysField);
}

inverse()

MInt maia::math::inverse ( MFloat **  a, MFloat **  ainv, MInt  n, const MFloat  epsilon  )

inline

Definition at line 587 of file maiamath.h.

                                                                             {
  MInt s;
  MInt pRow = 0; // pivot row
  MBool error = false;
  MFloat f;
  MFloat maximum;
  MInt pivot = 1;
 
  // add the unity matrix
  for(MInt i = 0; i < n; i++) {
    for(MInt j = 0; j < n; j++) {
      a[i][n + j] = F0;
      if(i == j) {
        a[i][n + j] = F1;
      }
    }
  }
 
  // Gauss algorithm
  error = false;
  s = 0;
  while(s < n) {
    maximum = fabs(a[s][s]);
    if(pivot) {
      pRow = s;
      for(MInt i = s + 1; i < n; i++) {
        if(fabs(a[i][s]) > maximum) {
          maximum = fabs(a[i][s]);
          pRow = i;
        }
      }
    }
    if(maximum < epsilon) {
      error = true;
    }
 
    if(error) {
      std::cerr << "Error in matrix inverse computation " << s << " " << a[s][s] << std::endl;
      for(MInt i = 0; i < n; i++) {
        for(MInt j = 0; j < n; j++) {
          std::cerr << a[i][j] << " ";
        }
        std::cerr << std::endl;
      }
      // mTerm(1, AT_, "Error in matrix inverse computation ");
      std::cerr << "Error in matrix inverse computation " << std::endl;
      return 0;
    }
 
    if(pivot) {
      if(pRow != s) // exchange rows if required
      {
        MFloat h;
        for(MInt j = s; j < 2 * n; j++) {
          h = a[s][j];
          a[s][j] = a[pRow][j];
          a[pRow][j] = h;
        }
      }
    }
 
    f = a[s][s];
    for(MInt j = s; j < 2 * n; j++) {
      a[s][j] = a[s][j] / f;
    }
 
    // elimination
    for(MInt i = 0; i < n; i++) {
      if(i != s) {
        f = -a[i][s];
        for(MInt j = s; j < 2 * n; j++) {
          a[i][j] += f * a[s][j];
        }
      }
    }
    s++;
  }
 
  if(error) {
    std::cerr << "Error 2 in inverse matrix computation" << std::endl;
    // mTerm(1,AT_,"Error 2 in inverse matrix computation");
    return 0;
  }
 
  // copy
  for(MInt i = 0; i < n; i++) {
    for(MInt j = 0; j < n; j++) {
      ainv[i][j] = a[i][n + j];
    }
  }
 
  return 1;
}

invert() [1/3]

void maia::math::invert	(	MFloat *	A,
		const MInt	m,
		const MInt	n
	)

Definition at line 171 of file maiamath.cpp.

                                                   {
  Eigen::Map<Eigen::MatrixXd> mA(A, m, n);
  mA = mA.inverse();
}

invert() [2/3]

template<class T >

void maia::math::invert	(	T &	A,
		T &	AInv,
		const MInt	m,
		const MInt	n
	)

Note: don't use this for new code!!!!!!!!!!!!!!!!!!!!!!!!!!! Use instead: A.colPivHouseholderQR.solve(b) (general case) A.LLT.solve(b) (if the matrix is symmetric positive definite) (fastest) A.HouseHolderQR.solve(b) (need for speed) (less accurate)

Get the pseudo-inverse of A by using QR

Parameters

A	system matrix (m x n)
AInv	pseudo-inverse
m	number of columns
n	number of rows

Author

Sven Berger

Definition at line 189 of file maiamath.cpp.

                                                       {
  if(n == 0 || m == 0) {
    std::cerr << globalDomainId() << " Warning: empty eq. sys. (" << m << "x" << n << "), skip." << std::endl;
    return;
  }
  Eigen::Matrix<MFloat, Eigen::Dynamic, Eigen::Dynamic> B(m, n);
  Eigen::Matrix<MFloat, Eigen::Dynamic, Eigen::Dynamic> inv(n, m);
  for(MInt i = 0; i < m; i++) {
    for(MInt j = 0; j < n; j++) {
      B(i, j) = A(i, j);
    }
  }
  //  if (n != m) {
  // system is overdetermined/underdetermined
  Eigen::CompleteOrthogonalDecomposition<Eigen::MatrixXd> orth(B);
  inv = orth.pseudoInverse();
  // todo labels:totest SSE optimization leads to false-positive during sanitizer run
  //  }
  //  else {
  //    inv = B.inverse();
  //  }
  for(MInt i = 0; i < m; i++) {
    for(MInt j = 0; j < n; j++) {
      AInv(i, j) = inv(i, j);
    }
  }
}

invert() [3/3]

template<class T >

void maia::math::invert	(	T &	A,
		T &	weights,
		T &	AInv,
		const MInt	m,
		const MInt	n
	)

Note: don't use this for new code!!!!!!!!!!!!!!!!!!!!!!!!!!! Use instead: A.colPivHouseholderQR.solve(b) (general case) A.LLT.solve(b) (if the matrix is symmetric positive definite) (fastest) A.HouseHolderQR.solve(b) (need for speed) (less accurate)

Get the pseudo-inverse of A by using QR

Parameters

A	system matrix (m x n)
AInv	pseudo-inverse
weights	weights applied to A
m	number of columns
n	number of rows

Author

Sven Berger

Definition at line 235 of file maiamath.cpp.

                                                                   {
  if(n == 0 || m == 0) {
    std::cerr << globalDomainId() << " Warning: empty eq. sys. (" << m << "x" << n << "), skip." << std::endl;
    return;
  }
  Eigen::Matrix<MFloat, Eigen::Dynamic, Eigen::Dynamic> B(m, n);
  Eigen::Matrix<MFloat, Eigen::Dynamic, Eigen::Dynamic> inv(n, m);
  for(MInt i = 0; i < m; i++) {
    for(MInt j = 0; j < n; j++) {
      B(i, j) = weights[i] * A(i, j);
    }
  }
  if(n != m) {
    // system is overdetermined/underdetermined
    Eigen::CompleteOrthogonalDecomposition<Eigen::MatrixXd> orth(B);
    inv = orth.pseudoInverse();
  } else {
    inv = B.inverse();
  }
  for(MInt k = 0; k < n; ++k) {
    for(MInt i = 0; i < m; ++i) {
      AInv(k, i) = inv(k, i) * weights(i);
    }
  }
}

invert< maia::tensor::Tensor< MFloat > >() [1/2]

template void maia::math::invert< maia::tensor::Tensor< MFloat > >	(	maia::tensor::Tensor< MFloat > &	,
		maia::tensor::Tensor< MFloat > &	,
		const	MInt,
		const	MInt
	)

invert< maia::tensor::Tensor< MFloat > >() [2/2]

template void maia::math::invert< maia::tensor::Tensor< MFloat > >	(	maia::tensor::Tensor< MFloat > &	,
		maia::tensor::Tensor< MFloat > &	,
		maia::tensor::Tensor< MFloat > &	,
		const	MInt,
		const	MInt
	)

invert< ScratchSpace< MFloat > >() [1/2]

template void maia::math::invert< ScratchSpace< MFloat > >	(	ScratchSpace< MFloat > &	,
		ScratchSpace< MFloat > &	,
		const	MInt,
		const	MInt
	)

invert< ScratchSpace< MFloat > >() [2/2]

template void maia::math::invert< ScratchSpace< MFloat > >	(	ScratchSpace< MFloat > &	,
		ScratchSpace< MFloat > &	,
		ScratchSpace< MFloat > &	,
		const	MInt,
		const	MInt
	)

invertR()

template<class T >

MInt maia::math::invertR	(	T &	A,
		T &	weights,
		T &	AInv,
		const MInt	m,
		const MInt	n
	)

Note: don't use this for new code!!!!!!!!!!!!!!!!!!!!!!!!!!! Use instead: A.colPivHouseholderQR.solve(b) (general case) A.LLT.solve(b) (if the matrix is symmetric positive definite) (fastest) A.HouseHolderQR.solve(b) (need for speed) (less accurate)

Get the pseudo-inverse of A by using QR

Parameters

A	system matrix (m x n)
AInv	pseudo-inverse
weights	weights applied to A
m	number of columns
n	number of rows

Returns

rank of the matrix A

Author

Sven Berger

Definition at line 280 of file maiamath.cpp.

                                                                    {
  if(n == 0 || m == 0) {
    std::cerr << globalDomainId() << " Warning: empty eq. sys. (" << m << "x" << n << "), skip." << std::endl;
    return -1;
  }
  Eigen::Matrix<MFloat, Eigen::Dynamic, Eigen::Dynamic> B(m, n);
  Eigen::Matrix<MFloat, Eigen::Dynamic, Eigen::Dynamic> inv(n, m);
  for(MInt i = 0; i < m; i++) {
    for(MInt j = 0; j < n; j++) {
      B(i, j) = weights[i] * A(i, j);
    }
  }
 
  Eigen::CompleteOrthogonalDecomposition<Eigen::MatrixXd> orth(B);
  inv = orth.pseudoInverse();
  const MInt rank = orth.rank();
 
  for(MInt k = 0; k < n; ++k) {
    for(MInt i = 0; i < m; ++i) {
      AInv(k, i) = inv(k, i) * weights(i);
    }
  }
  return rank;
}

invertR< ScratchSpace< MFloat > >()

template MInt maia::math::invertR< ScratchSpace< MFloat > >	(	ScratchSpace< MFloat > &	,
		ScratchSpace< MFloat > &	,
		ScratchSpace< MFloat > &	,
		const	MInt,
		const	MInt
	)

lincos()

MFloat maia::math::lincos ( MFloat  arg )

inline

Definition at line 426 of file maiamath.h.

                                 {
  MFloat SIGN = F1;
  // lookuptable has 91 entries from 0 to pi/2
  const MFloat darg = PIB2 / 90;
 
  arg = std::fmod(std::abs(arg), F2 * PI);
 
  if(arg > PI) {
    SIGN *= -F1;
    arg -= PI;
  }
 
  if(arg > PIB2) {
    SIGN *= -F1;
    arg = PI - arg;
  }
 
  const auto lowerInd = static_cast<MInt>(arg / darg);
  const MFloat lowerVal = FTRIG[lowerInd];
  const MFloat higherVal = FTRIG[lowerInd + 1];
  const MFloat deltaArg = arg - lowerInd * darg;
  const MFloat fac = deltaArg / darg;
 
  return (lowerVal * (F1 - fac) + higherVal * fac) * SIGN;
}

linSpace()

std::vector< MFloat > maia::math::linSpace ( const MFloat  start, const MFloat  end, const MInt  num  )

inline

Definition at line 903 of file maiamath.h.

                                                                                        {
  std::vector<MFloat> linspaced(num);
  MFloat delta = (end - start) / (MFloat(num) - F1);
  for(auto i = 0; i < num; i++) {
    linspaced[i] = start + delta * i;
  }
  return linspaced;
}

matrixVectorProduct()

void maia::math::matrixVectorProduct ( MFloat *  c, MFloatScratchSpace &  A, MFloat *  b  )

inline

Author

Lennart Schneiders

Definition at line 561 of file maiamath.h.

                                                                             {
  for(MInt i = 0; i < A.size0(); i++) {
    c[i] = F0;
    for(MInt j = 0; j < A.size1(); j++) {
      c[i] += A(i, j) * b[j];
    }
  }
}

matrixVectorProductTranspose()

void maia::math::matrixVectorProductTranspose ( MFloat *  c, MFloatScratchSpace &  A, MFloat *  b  )

inline

Author

Lennart Schneiders

Definition at line 574 of file maiamath.h.

                                                                                      {
  for(MInt i = 0; i < A.size1(); i++) {
    c[i] = F0;
    for(MInt j = 0; j < A.size0(); j++) {
      c[i] += A(j, i) * b[j];
    }
  }
}

multiplyMatrices()

void maia::math::multiplyMatrices ( MFloatScratchSpace &  m1, MFloatScratchSpace &  m2, MFloatScratchSpace &  result, MInt  m1_n, MInt  m1_m, MInt  m2_n, MInt  m2_m  )

inline

Definition at line 297 of file maiamath.h.

                                                              {
  // m1 has the dimension n x m
  // m2 has the dimension n x m
  if(m1_m != m2_n) mTerm(1, AT_, "Dimension of matrix multiplication does not match!!!");
  // init
  for(MInt i = 0; i < m1_n; i++) {
    for(MInt j = 0; j < m2_m; j++) {
      result(i, j) = F0;
    }
  }
 
  // multiply
  for(MInt i = 0; i < m1_n; i++) {
    for(MInt j = 0; j < m2_m; j++) {
      for(MInt k = 0; k < m1_m; k++) {
        result(i, j) += m1(i, k) * m2(k, j);
      }
    }
  }
}

multiplyMatricesSq()

void maia::math::multiplyMatricesSq ( MFloatScratchSpace &  m1, MFloatScratchSpace &  m2, MFloatScratchSpace &  result, MInt  dim  )

inline

Definition at line 279 of file maiamath.h.

                                                                                                                     {
  // init
  for(MInt i = 0; i < dim; i++) {
    for(MInt j = 0; j < dim; j++) {
      result(i, j) = 0.0;
    }
  }
 
  // multiply
  for(MInt i = 0; i < dim; i++) {
    for(MInt j = 0; j < dim; j++) {
      for(MInt k = 0; k < dim; k++) {
        result(i, j) += m1(i, k) * m2(k, j);
      }
    }
  }
}

multiplySparseMatrixVector()

void maia::math::multiplySparseMatrixVector	(	MFloat *	A_coeff,
		MInt **	pos,
		const MInt	n,
		const MInt	m,
		MFloat *	b,
		MFloat *	x
	)

Solve b = A * x for dense matrix A using Cholesky or LU

Parameters

A_coeff	non-zero coefficients of the system matrix
pos	holds the position of the coefficients in the matrix A
n	number of non-zero coefficients
m	size of system matrix A
b	unknown vector to be calculated by the product A * x
x	known vector

Author

Moritz Waldmann

Definition at line 431 of file maiamath.cpp.

                                                                                                               {
  Eigen::Map<MFloatVectorXd> B(b, m);
  Eigen::Map<MFloatVectorXd> X(x, m);
 
  MFloatMatrixXd A = MFloatMatrixXd::Zero(m, m);
 
  for(MInt i = 0; i < n; i++) {
    A(pos[i][0], pos[i][1]) = A_coeff[i];
  }
  B = A * X;
}

norm() [1/2]

template<typename T , std::size_t N>

MFloat maia::math::norm ( const std::array< T, N > &  u )

inline

Length of a vector (L2-norm)

Template Parameters

T	Type
N	Length of vector

Parameters

u

Vector

Returns

Magnitude of vector

Author

Sven Berger

Definition at line 148 of file maiamath.h.

                                            {
  static_assert(N > 1, "ERROR: Invalid norm call!");
 
  if(N == 2) {
    return std::sqrt(u[0] * u[0] + u[1] * u[1]);
  }
  if(N == 3) {
    return std::sqrt(u[0] * u[0] + u[1] * u[1] + u[2] * u[2]);
  }
  return std::sqrt(std::accumulate(u.begin(), u.end(), 0.0, [](const MFloat& a, const T& b) { return a + b * b; }));
}

norm() [2/2]

template<typename T >

MFloat maia::math::norm ( const T *const  u, const MInt  N  )

inline

Length of a vector (L2-norm)

Template Parameters

T

Type

Parameters

N	Length of vector
u	Vector

Returns

Magnitude of vector

Author

Sven Berger

Definition at line 167 of file maiamath.h.

                                                   {
  ASSERT(N > 1, "ERROR: Invalid norm call!");
 
  if(N == 2) {
    return std::sqrt(u[0] * u[0] + u[1] * u[1]);
  }
  if(N == 3) {
    return std::sqrt(u[0] * u[0] + u[1] * u[1] + u[2] * u[2]);
  }
 
  MFloat tmp = 0;
  for(MInt i = 0; i < N; i++) {
    tmp += u[i] * u[i];
  }
 
  return std::sqrt(tmp);
}

normalize() [1/2]

template<typename T , std::size_t N>

void maia::math::normalize ( std::array< T, N > &  u )

inline

Normalize a given vector (inplace)

Template Parameters

T	Type
N	Length of the vector

Parameters

u	Vector to be normalized

Author

Sven Berger

Definition at line 191 of file maiamath.h.

                                         {
  const MFloat inverse = 1.0 / norm(u);
 
  if(N == 2) {
    u[0] *= inverse;
    u[1] *= inverse;
    return;
  }
  if(N == 3) {
    u[0] *= inverse;
    u[1] *= inverse;
    u[2] *= inverse;
    return;
  }
 
  std::for_each(u.begin(), u.end(), [inverse](T& a) { a *= inverse; });
}

normalize() [2/2]

template<typename T >

void maia::math::normalize ( T *const  u, const MInt  N  )

inline

Normalize a given vector (inplace)

Template Parameters

T

Type

Parameters

N	Length of the vector
u	Vector to be normalized

Author

Sven Berger

Definition at line 215 of file maiamath.h.

                                                {
  const MFloat inverse = 1.0 / norm(u, N);
 
  if(N == 2) {
    u[0] *= inverse;
    u[1] *= inverse;
    return;
  }
  if(N == 3) {
    u[0] *= inverse;
    u[1] *= inverse;
    u[2] *= inverse;
    return;
  }
 
  for(MInt i = 0; i < N; i++) {
    u[i] *= inverse;
  }
}

normalized()

template<typename T , std::size_t N>

std::array< T, N > maia::math::normalized ( const std::array< T, N > &  u )

inline

Normalize a given vector

Template Parameters

T	Type
N	Length of the vector

Parameters

u	Vector to be normalized

Returns

New normalized vector

Author

Sven Berger

Definition at line 242 of file maiamath.h.

                                                          {
  std::array<T, N> n(u);
  normalize(n);
  return n;
}

quatMult()

void maia::math::quatMult ( const MFloat *const  qA, const MFloat *const  qB, MFloat *const  qC  )

inline

Definition at line 482 of file maiamath.h.

                                                                                       {
  MFloat crossProduct[3]{};
  cross(qA, qB, &crossProduct[0]);
 
  const MFloat dot = std::inner_product(qA, &qA[3], qB, 0.0);
 
  for(MInt n = 0; n < 3; n++) {
    qC[n] = qA[3] * qB[n] + qB[3] * qA[n] + crossProduct[n];
  }
  qC[3] = qA[3] * qB[3] - dot;
}

quickSort()

void maia::math::quickSort ( MInt *  a, MInt  start, MInt  end  )

inline

Definition at line 716 of file maiamath.h.

{ maia::math::quickSortImpl(a, start, end); }

quickSortImpl()

void maia::math::quickSortImpl ( MInt *  a, MInt  start, MInt  end  )

inline

Definition at line 704 of file maiamath.h.

                                                         {
  MInt pivot;
 
  if(start >= end) {
    return;
  }
 
  pivot = maia::math::quickSortPartition(a, start, end);
  maia::math::quickSort(a, start, pivot - 1);
  maia::math::quickSort(a, pivot + 1, end);
}

quickSortPartition()

MInt maia::math::quickSortPartition ( MInt *  a, MInt  start, MInt  end  )

inline

Definition at line 681 of file maiamath.h.

                                                              {
  MInt i, temp;
  //---
 
  i = start - 1;
 
  for(MInt j = start; j < end; j++) {
    if(a[j] <= a[end]) {
      i++;
      // swap elements i and j
      temp = a[i];
      a[i] = a[j];
      a[j] = temp;
    }
  }
  // swap the pivot
  temp = a[i + 1];
  a[i + 1] = a[end];
  a[end] = temp;
 
  return i + 1;
}

randnormal()

MFloat maia::math::randnormal ( MFloat  mu, MFloat  sigma  )

inline

Definition at line 502 of file maiafftw.h.

                                                  {
  TRACE();
 
  static MBool deviateAvailable = false; //        flag
  static float storedDeviate;            //        deviate from previous calculation
  MFloat polar, rsquared, var1, var2;
 
  //        If no deviate has been stored, the polar Box-Muller transformation is
  //        performed, producing two independent normally-distributed random
  //        deviates.  One is stored for the next round, and one is returned.
  if(!deviateAvailable) {
    //        choose pairs of uniformly distributed deviates, discarding those
    //        that don't fall within the unit circle
    do {
      var1 = 2.0 * (MFloat(rand()) / MFloat(RAND_MAX)) - 1.0;
      var2 = 2.0 * (MFloat(rand()) / MFloat(RAND_MAX)) - 1.0;
      rsquared = var1 * var1 + var2 * var2;
    } while(rsquared >= 1.0 || approx(rsquared, 0.0, MFloatEps));
 
    //        calculate polar tranformation for each deviate
    polar = sqrt(-2.0 * log(rsquared) / rsquared);
 
    //        store first deviate and set flag
    storedDeviate = var1 * polar;
    deviateAvailable = true;
 
    //        return second deviate
 
    return var2 * polar * sigma + mu;
  }
 
  //        If a deviate is available from a previous call to this function, it is
  //        returned, and the flag is set to false.
  else {
    deviateAvailable = false;
 
    return storedDeviate * sigma + mu;
  }
}

RBF()

MFloat maia::math::RBF ( const MFloat  R, const MFloat  R0  )

inline

Author

Lennart Schneiders

Note

important: provide rapid decrease of return value with R, otherwise stability issues with weighted least squares when inlcuding diagonal neighbor cells!

suggested value for R0 is cellLength

Definition at line 873 of file maiamath.h.

{ return (F1 / sqrt(F1 + POW2(R / R0))); }

removeDoubleEntries()

MInt maia::math::removeDoubleEntries ( MInt *  a, MInt  size  )

inline

Definition at line 722 of file maiamath.h.

                                                    {
  for(MInt i = 1; i < size; i++) {
    // find double entry
    if(a[i] != a[i - 1]) {
      continue;
    }
 
    // shift entries forward
    for(MInt k = i; k < size - 1; k++) {
      a[k] = a[k + 1];
    }
    size--;
    i--;
  }
 
  return size;
}

rotation2quaternion()

void maia::math::rotation2quaternion ( MFloat *  rotation, MFloat *  quaternion  )

inline

Author

Definition at line 550 of file maiamath.h.

                                                                      {
  quaternion[0] = cos(F1B2 * rotation[1]) * cos(F1B2 * (rotation[2] + rotation[0]));
  quaternion[1] = sin(F1B2 * rotation[1]) * sin(F1B2 * (rotation[2] - rotation[0]));
  quaternion[2] = sin(F1B2 * rotation[1]) * cos(F1B2 * (rotation[2] - rotation[0]));
  quaternion[3] = cos(F1B2 * rotation[1]) * sin(F1B2 * (rotation[2] + rotation[0]));
}

sgn()

template<typename T >

MInt maia::math::sgn ( T  val )

inline

Definition at line 495 of file maiamath.h.

                       {
  return (T(0) < val) - (val < T(0));
}

solveDenseMatrix()

void maia::math::solveDenseMatrix	(	MFloat *	A_coeff,
		MInt **	pos,
		const MInt	n,
		const MInt	m,
		MFloat *	b,
		MFloat *	x
	)

Solve x = A^(-1) * b for dense matrix A using LU decomposition

Parameters

A_coeff	non-zero coefficients of the system matrix
pos	holds the position of the coefficients in the matrix A
n	number of non-zero coefficients
m	size of system matrix A
b	known vector to be calculated by the product A * x
x	unknown vector to be calculated by A^(-1) * b

Author

Moritz Waldmann

Definition at line 400 of file maiamath.cpp.

                                                                                                     {
  Eigen::Map<MFloatVectorXd> B(b, m);
  Eigen::Map<MFloatVectorXd> X(x, m);
 
  MFloatMatrixXd A = MFloatMatrixXd::Zero(m, m);
 
  for(MInt i = 0; i < n; i++) {
    A(pos[i][0], pos[i][1]) = A_coeff[i];
  }
 
  Eigen::FullPivLU<MFloatMatrixXd> lu;
  lu.compute(A);
 
  X = lu.solve(B);
  MBool a_solution_exists = (A * X).isApprox(B, 1e-10);
 
  if(a_solution_exists) {
    std::cout << "SOLUTION IS CORRECT" << std::endl;
  } else {
    std::cout << "SOLUTION IS NOT CORRECT" << std::endl;
  }
}

solveQR()

template<MInt nDim>

void maia::math::solveQR	(	std::array< std::array< MFloat, nDim >, nDim > &	A_,
		std::array< MFloat, nDim > &	b_
	)

Solve linear system using QR (square matrix)

Parameters

A	system matrix
b	rhs

Template Parameters

nDim	dimension of the matrix

Author

Sven Berger

Definition at line 449 of file maiamath.cpp.

                                                                                     {
  TRACE();
 
  Eigen::Matrix<MFloat, nDim, nDim> A(&A_[0][0]);
  Eigen::Vector<MFloat, nDim> b(&b_[0]);
  b = A.colPivHouseholderQr().solve(b);
  std::memcpy(&b_[0], b.data(), nDim * sizeof(MFloat));
}

solveQR< 2 >()

template void maia::math::solveQR< 2 >	(	std::array< std::array< MFloat, 2 >, 2 > &	,
		std::array< MFloat, 2 > &
	)

solveQR< 3 >()

template void maia::math::solveQR< 3 >	(	std::array< std::array< MFloat, 3 >, 3 > &	,
		std::array< MFloat, 3 > &
	)

solveSparseMatrix()

void maia::math::solveSparseMatrix	(	MFloat *	A_coeff,
		MInt **	pos,
		const MInt	n,
		const MInt	m,
		MFloat *	b,
		MFloat *	x
	)

Solve b = A * x for sparse matrix A using LU decomposition

Parameters

A_coeff	non-zero coefficients of the system matrix
pos	holds the position of the coefficients in the matrix A
n	number of non-zero coefficients
m	size of system matrix A
b	known vector to be calculated by the product A * x
x	unknown vector to be calculated by A^(-1) * b

Author

Moritz Waldmann

Definition at line 315 of file maiamath.cpp.

                                                                                                      {
  std::vector<MTriplet> coefficients;
  Eigen::Map<MFloatVectorXd> B(b, m);
  Eigen::Map<MFloatVectorXd> X(x, m);
 
  coefficients.reserve(n);
  for(MInt i = 0; i < n; i++) {
    coefficients.emplace_back(MTriplet(pos[i][0], pos[i][1], A_coeff[i]));
  }
 
  // set matrix from coefficients
  Eigen::SparseMatrix<MFloat, Eigen::ColMajor> A(m, m);
  A.setFromTriplets(coefficients.begin(), coefficients.end());
 
  Eigen::SparseLU<Eigen::SparseMatrix<MFloat, Eigen::ColMajor>, Eigen::COLAMDOrdering<MInt>> slu;
  slu.compute(A);
 
  X = slu.solve(B);
  MBool a_solution_exists = (A * X).isApprox(B, 1e-10);
 
  MFloat resultA = F1B2 * X.transpose() * A * X;
  MFloat resultB = X.transpose() * B;
 
  MFloat result = resultA - resultB;
 
  std::cout << "Energie norm is " << resultA << " " << resultB << " " << result << std::endl;
 
  if(a_solution_exists) {
    std::cerr << "SOLUTION IS CORRECT" << std::endl;
  } else {
    std::cerr << "SOLUTION IS NOT CORRECT" << std::endl;
  }
}

solveSparseMatrixIterative()

void maia::math::solveSparseMatrixIterative	(	MFloat *	A_coeff,
		MInt **	pos,
		const MInt	n,
		const MInt	m,
		MFloat *	b,
		MFloat *	x
	)

Solve b = A^(-1) * x for sparse matrix A using BiCGSTAB

Parameters

A_coeff	non-zero coefficients of the system matrix
pos	holds the position of the coefficients in the matrix A
n	number of non-zero coefficients
m	size of system matrix A
b	known vector to be calculated by the product A * x
x	unknown vector to be calculated by A^(-1) * b

Author

Moritz Waldmann

Definition at line 357 of file maiamath.cpp.

                                                                                                               {
  std::vector<MTriplet> coefficients;
  Eigen::Map<MFloatVectorXd> B(b, m);
  Eigen::Map<MFloatVectorXd> X(x, m);
 
  coefficients.reserve(n);
  for(MInt i = 0; i < n; i++) {
    coefficients.emplace_back(MTriplet(pos[i][0], pos[i][1], A_coeff[i]));
  }
 
  // set matrix from coefficients
  Eigen::SparseMatrix<MFloat, Eigen::ColMajor> A(m, m);
  A.setFromTriplets(coefficients.begin(), coefficients.end());
 
  Eigen::BiCGSTAB<Eigen::SparseMatrix<MFloat, Eigen::ColMajor>> BCGST;
  BCGST.compute(A);
 
  X = BCGST.solve(B);
 
  MBool a_solution_exists = (A * X).isApprox(B, 1e-10);
 
  MFloat resultA = F1B2 * X.transpose() * A * X;
  MFloat resultB = X.transpose() * B;
 
  MFloat result = resultA - resultB;
 
  std::cout << "Energie norm is " << resultA << " " << resultB << " " << result << std::endl;
 
  if(a_solution_exists) {
    std::cout << "SOLUTION IS CORRECT" << std::endl;
  } else {
    std::cout << "SOLUTION IS NOT CORRECT" << std::endl;
  }
}

sortEigenVectors()

void maia::math::sortEigenVectors ( MFloat  A[3][3], MFloat  w[3]  )

inline

void sortEigenVectors(MFloat A[3][3], MFloat w[3])

   Parameters:
   A: The symmetric input matrix
   w: Storage buffer for eigenvalues

Definition at line 460 of file maiamath.h.

                                                          {
  const MInt dim = 3;
  MInt k;
  for(MInt i = 0; i < (dim - 1); ++i) {
    MFloat p = w[k = i];
    for(MInt j = i; j < dim; ++j) {
      if(w[j] >= p) {
        p = w[k = j];
      }
    }
    if(k != i) {
      w[k] = w[i];
      w[i] = p;
      for(MInt j = 0; j < dim; ++j) {
        p = A[j][i];
        A[j][i] = A[j][k];
        A[j][k] = p;
      }
    }
  }
}

svd()

std::vector< MFloat > maia::math::svd	(	MFloat *const	A,
		MFloat *const	b,
		const MInt	m,
		const MInt	n,
		MFloat *const	x
	)

Definition at line 514 of file maiamath.cpp.

                                                                                                     {
  Eigen::Map<Eigen::Matrix<MFloat, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> A_eigen(A, m, n);
 
  Eigen::Map<MFloatVectorXd> b_eigen(b, m);
  //      Eigen::Map<maia::math::MFloatVector<Eigen::Dynamic>> S_eigen(&S[0], std::min(m, n));
 
#ifdef useEigenOld
  Eigen::BDCSVD<Eigen::Matrix<MFloat, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> SVD =
      A_eigen.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV);
#else
  Eigen::BDCSVD<Eigen::Matrix<MFloat, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>,
                Eigen::ComputeThinU | Eigen::ComputeThinV>
      SVD = A_eigen.bdcSvd<Eigen::ComputeThinU | Eigen::ComputeThinV>();
#endif
 
  // TODO labels:CONTROLLER this should work as well
  MFloatVectorXd x_eigen = SVD.solve(b_eigen);
  //      Eigen::Matrix<MFloat, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> x_eigen =
  //          SVD.solve(b_eigen);
  MFloatVectorXd S = SVD.singularValues();
 
  for(MInt i = 0; i < n; i++) {
    x[i] = x_eigen(i);
  }
 
  std::vector<MFloat> singularValues(S.size());
  for(MInt i = 0; i < S.size(); i++) {
    singularValues[i] = S[i];
  }
  return singularValues;
}

triadImag()

MFloat maia::math::triadImag ( fftw_complex &  a, fftw_complex &  b, fftw_complex &  c  )

inline

Definition at line 61 of file maiafftw.h.

                                                                           {
  return (a[0] * b[1] * c[0] + a[1] * b[0] * c[0] - a[0] * b[0] * c[1] + a[1] * b[1] * c[1]);
}

vecAdd() [1/2]

template<MInt nDim, typename T , typename U >

T * maia::math::vecAdd ( T *const  result, const U *const  a  )

inline

Adds arbritary number of vectors of dimension nDim and returns result in vector result. Vector result can contain an offset. Usage: Say you want to add vectors a, b, c of dimension 3 and store result in a: vecAdd<3>(a, b, c)

Definition at line 74 of file maiamath.h.

                                                    {
  std::transform(a, a + nDim, result, result, std::plus<T>());
  return result;
}

vecAdd() [2/2]

template<MInt nDim, typename T , typename U , typename... Ts>

T * maia::math::vecAdd ( T *const  result, const U *const  a, const Ts...  b  )

inline

Definition at line 79 of file maiamath.h.

                                                                   {
  std::transform(a, a + nDim, result, vecAdd<nDim>(result, b...), std::plus<T>());
  return result;
}

vecAvg()

template<MInt nDim, typename T , typename... Ts>

void maia::math::vecAvg ( T *const  M, const Ts *const ...  args  )

inline

Definition at line 86 of file maiamath.h.

                                                        {
  std::fill_n(M, nDim, 0.0);
  maia::math::vecAdd<nDim>(M, args...);
  constexpr std::size_t n = sizeof...(Ts);
  std::transform(M, M + nDim, M, std::bind(std::multiplies<T>(), std::placeholders::_1, 1.0 / n));
}