#include <UT_MultigridArray.h>

Public Member Functions
	UT_MultigridArrayT ()

	~UT_MultigridArrayT ()=default

	UT_MultigridArrayT (const UT_MultigridArrayT &a)
	Construct as copy. More...

void	init (const UT_Vector3I &res, const UT_Vector3T< T > &spacing, const UT_Vector3I &boundariesNeg, const UT_Vector3I &boundariesPos, int level=0, UT_Vector3I oddCoarsenings=UT_Vector3I(0, 0, 0))

void	initFromVoxels (const UT_VoxelArray< T > &f, const UT_Vector3T< T > &spacing, const UT_Vector3I &boundariesNeg, const UT_Vector3I &boundariesPos)

void	copyToVoxels (UT_VoxelArray< T > &f) const

void	match (const UT_MultigridArrayT &src)

bool	isInit ()
	Whether this array has been initialized - default constructor does not. More...

T	operator() (int i, int j, int k) const

T &	operator() (int i, int j, int k)

T	operator() (const UT_Vector3I &idx) const

T &	operator() (const UT_Vector3I &idx)

T	ghostValueNeg (int axis, const UT_Vector3I &idx, const UT_Vector3I &idxadj) const

T	ghostValuePos (int axis, const UT_Vector3I &idx, const UT_Vector3I &idxadj) const

T	ghostValueNeg (int axis, T U, T Uadj) const

T	ghostValuePos (int axis, T U, T Uadj) const

UT_MultigridArrayT< T > &	operator= (const UT_MultigridArrayT &a)

UT_MultigridArrayT< T > &	operator+= (const UT_MultigridArrayT &a)

void	zero ()

T	norm (int n) const

exint	getRes (int axis) const

UT_Vector3I	getRes () const

const UT_Vector3T< T > &	getSpacing () const

const UT_Vector3I &	getBoundariesNeg () const

const UT_Vector3I &	getBoundariesPos () const

int	getLevel () const

UT_Vector3I	getParity () const

exint	numElements () const
	Returns the total number of elements in the array. More...

exint	getMaxCoarse (exint minPerAxis) const

bool	is1D () const
	Returns whether this is a 1-D array. More...

bool	shouldMultiThread () const

const T *	data () const

T *	data ()

void	computeNorms (fpreal64 &norminf, fpreal64 &norm2) const
	Computes the 0 and 2-norm of this in one pass. More...

void	computeResidualNorms (const UT_MultigridArrayT &x, fpreal64 &norminf, fpreal64 &norm2) const

	THREADED_METHOD2_CONST (UT_MultigridArrayT, shouldMultiThread(), subtractApplyLaplacian, const UT_MultigridArrayT &, x, UT_MultigridArrayT &, r) void subtractApplyLaplacianPartial(const UT_MultigridArrayT &x
	Computes rhe residual r = b - A * x, where b == *this. More...

	THREADED_METHOD2_CONST (UT_MultigridArrayT, shouldMultiThread(), smoothLaplacianGaussSeidel, int, parity, UT_MultigridArrayT &, x) void smoothLaplacianGaussSeidelPartial(int parity

	THREADED_METHOD2_CONST (UT_MultigridArrayT, shouldMultiThread(), coarsenAlongAxis, UT_MultigridArrayT &, uh, int, axis) void coarsenAlongAxisPartial(UT_MultigridArrayT &uh

	THREADED_METHOD2_CONST (UT_MultigridArrayT, shouldMultiThread(), interpolateAlongAxis, UT_MultigridArrayT &, u, int, axis) void interpolateAlongAxisPartial(UT_MultigridArrayT &u

template<typename S >
void	buildMatrix (UT_MatrixT< S > &A) const

void	directSolve (UT_MultigridArrayT &x) const

UT_Vector3I	coarsen (exint minPerAxis, UT_MultigridArrayT &uh) const

void	interpolate (const UT_Vector3I &coarsenAxis, const UT_Vector3I &parity, UT_MultigridArrayT &u) const

void	vcycle (exint minPerAxis, int nSmoothDown, int nSmoothUp, UT_MultigridArrayT &x, bool smoothTopLevelDown=true) const

void	fullMultigrid (exint minPerAxis, int nSmoothDown, int nSmoothUp, UT_MultigridArrayT &x) const

int	solvePoisson (fpreal64 abstol, fpreal64 reltol, int miniter, int maxiter, UT_MultigridArrayT &x, UT_ValArray< fpreal64 > resnorminf, UT_ValArray< fpreal64 > resnorm2, bool startFMG=true) const

Public Attributes
UT_MultigridArrayT &	r

UT_MultigridArrayT const UT_JobInfo &info	const

UT_MultigridArrayT &	x

int	axis

int const UT_JobInfo &info	const

Protected Member Functions
bool	solvePoissonCholesky (UT_MultigridArrayT &x) const

bool	solvePoissonLU (UT_MultigridArrayT &x) const

bool	solvePoissonSVD (UT_MultigridArrayT &x) const

void	initStorage ()

void	initLaplacian ()

	THREADED_METHOD2_CONST (UT_MultigridArrayT, shouldMultiThread(), computeNormsInternal, fpreal64 , norminf, fpreal64 , norm2) void computeNormsInternalPartial(fpreal64 *norminf

	THREADED_METHOD3_CONST (UT_MultigridArrayT, shouldMultiThread(), computeResidualNormsInternal, const UT_MultigridArrayT &, x, fpreal64 , norminf, fpreal64 , norm2) void computeResidualNormsInternalPartial(const UT_MultigridArrayT &x

Protected Attributes
fpreal64 *	norm2

fpreal64 const UT_JobInfo &info	const

fpreal64 *	norminf

fpreal64 fpreal64 *	norm2

fpreal64 fpreal64 const UT_JobInfo &info	const

UT_VectorT< T >	myVec
	UT_VectorT that holds the array data. More...

UT_Vector3I	myRes
	Number of elements in each dimension. More...

UT_Vector3T< T >	mySpacing
	Grid spacing in each dimension. More...

exint	myYStride
	Array stride in y and z axes. More...

exint	myZStride

T	myInvdx2
	Laplacian operator values calculated based on grid spacing. More...

T	myInvdy2

T	myInvdz2

T	myDiag

T	myOmega

UT_Vector3I	myBoundariesNeg

UT_Vector3I	myBoundariesPos

UT_Vector3T< T >	myBoundModNeg

UT_Vector3T< T >	myBoundModPos

UT_Vector3T< T >	myBoundAdjModNeg

UT_Vector3T< T >	myBoundAdjModPos

UT_Vector3T< T >	myDiagModNeg

UT_Vector3T< T >	myDiagModPos

int	myLevel

bool	myAllOpen

bool	myAllClosed

UT_Vector3I	myOddCoarsenings

Static Protected Attributes
static const exint	PARALLEL_BLOCK_SIZE_VEC = 512

static const exint	PARALLEL_BLOCK_SIZE_Z = 8

Detailed Description

template<typename T>
class UT_MultigridArrayT< T >

UT_MultigridArrayT This provides an array that can solve the 2D or 3D Poisson equation using the multigrid algorithm. There are currently only explicit template instantations for fpreal32 and fpreal64.

The created array has elements indexed from 0, ie: [0..xdiv-1].

Definition at line 44 of file UT_MultigridArray.h.

Constructor & Destructor Documentation

template<typename T>

UT_MultigridArrayT< T >::UT_MultigridArrayT ( )

Default constructor, array will not be initialized after this, i.e. isInit() == false

template<typename T>

UT_MultigridArrayT< T >::~UT_MultigridArrayT ( )

default

template<typename T>

UT_MultigridArrayT< T >::UT_MultigridArrayT ( const UT_MultigridArrayT< T > & a )

inline

Construct as copy.

Definition at line 54 of file UT_MultigridArray.h.

Member Function Documentation

template<typename T>

template<typename S >

void UT_MultigridArrayT< T >::buildMatrix ( UT_MatrixT< S > & A ) const

Build the dense matrix in A that corresponds to this grid's Laplacian operator.

template<typename T>

UT_Vector3I UT_MultigridArrayT< T >::coarsen	(	exint	minPerAxis,
		UT_MultigridArrayT< T > &	uh
	)		const

Coarsen this array into a lower resolution version of the grid. coarsenAxis is a vector of booleans specifying which axes to coarsen. Where coarsenAxis[axis] is true, ensure that uh[axis] = getRes(axis) / 2

template<typename T>

void UT_MultigridArrayT< T >::computeNorms	(	fpreal64 &	norminf,
		fpreal64 &	norm2
	)		const

Computes the 0 and 2-norm of this in one pass.

template<typename T>

void UT_MultigridArrayT< T >::computeResidualNorms	(	const UT_MultigridArrayT< T > &	x,
		fpreal64 &	norminf,
		fpreal64 &	norm2
	)		const

Computes the inf and 2-norm of the residual r = b - A * x, where b == *this

template<typename T>

void UT_MultigridArrayT< T >::copyToVoxels ( UT_VoxelArray< T > & f ) const

Copy the array data to the provided UT_VoxelArray. The UT_VoxelArray will be resized if the number of voxels doesn't match the number of elements in this array.

template<typename T>

const T* UT_MultigridArrayT< T >::data ( ) const

inline

Definition at line 236 of file UT_MultigridArray.h.

template<typename T>

T* UT_MultigridArrayT< T >::data ( )

inline

Definition at line 241 of file UT_MultigridArray.h.

template<typename T>

void UT_MultigridArrayT< T >::directSolve ( UT_MultigridArrayT< T > & x ) const

Directly solve the Poisson equation that this array represents, using either Cholesky, LU, or SVD decomposition, depending on the nature of the Poisson problem. Should only be called for very coarse grids.

template<typename T>

void UT_MultigridArrayT< T >::fullMultigrid	(	exint	minPerAxis,
		int	nSmoothDown,
		int	nSmoothUp,
		UT_MultigridArrayT< T > &	x
	)		const

Performs recursive Full Multigrid to generate an initial guess for the solution of the Poisson equation represented by this array. At each level, nSmoothDown and nSmoothUp iterations of smoothing are done.

template<typename T>

const UT_Vector3I& UT_MultigridArrayT< T >::getBoundariesNeg ( ) const

inline

Definition at line 188 of file UT_MultigridArray.h.

template<typename T>

const UT_Vector3I& UT_MultigridArrayT< T >::getBoundariesPos ( ) const

inline

Definition at line 193 of file UT_MultigridArray.h.

template<typename T>

int UT_MultigridArrayT< T >::getLevel ( ) const

inline

Definition at line 198 of file UT_MultigridArray.h.

template<typename T>

exint UT_MultigridArrayT< T >::getMaxCoarse ( exint minPerAxis ) const

Calculates the maximum coarse grid elements in the coarsest grid, based on the minPerAxis parameter.

template<typename T>

UT_Vector3I UT_MultigridArrayT< T >::getParity ( ) const

inline

Definition at line 203 of file UT_MultigridArray.h.

template<typename T>

exint UT_MultigridArrayT< T >::getRes ( int axis ) const

inline

Definition at line 173 of file UT_MultigridArray.h.

template<typename T>

UT_Vector3I UT_MultigridArrayT< T >::getRes ( ) const

inline

Definition at line 178 of file UT_MultigridArray.h.

template<typename T>

const UT_Vector3T<T>& UT_MultigridArrayT< T >::getSpacing ( ) const

inline

Definition at line 183 of file UT_MultigridArray.h.

template<typename T>

T UT_MultigridArrayT< T >::ghostValueNeg	(	int	axis,
		const UT_Vector3I &	idx,
		const UT_Vector3I &	idxadj
	)		const

inline

Returns the ghost value for the negative side boundaries, given the values at idx and idxadj, which should represent the index for a border cell and the one adjacent to it along the provided axis.

Definition at line 126 of file UT_MultigridArray.h.

template<typename T>

T UT_MultigridArrayT< T >::ghostValueNeg	(	int	axis,
		T	U,
		T	Uadj
	)		const

inline

Returns the ghost value for the negative side boundaries, given the provided values U and Uadj, which should be the values at the border cell and the one adjacent to it along the provided axis.

Definition at line 146 of file UT_MultigridArray.h.

template<typename T>

T UT_MultigridArrayT< T >::ghostValuePos	(	int	axis,
		const UT_Vector3I &	idx,
		const UT_Vector3I &	idxadj
	)		const

inline

Returns the ghost value for the positive side boundaries, given the values at idx and idxadj, which should represent the index for a border cell and the one adjacent to it along the provided axis.

Definition at line 136 of file UT_MultigridArray.h.

template<typename T>

T UT_MultigridArrayT< T >::ghostValuePos	(	int	axis,
		T	U,
		T	Uadj
	)		const

inline

Returns the ghost value for the positive side boundaries, given the provided values U and Uadj, which should be the values at the border cell and the one adjacent to it along the provided axis.

Definition at line 154 of file UT_MultigridArray.h.

template<typename T>

void UT_MultigridArrayT< T >::init	(	const UT_Vector3I &	res,
		const UT_Vector3T< T > &	spacing,
		const UT_Vector3I &	boundariesNeg,
		const UT_Vector3I &	boundariesPos,
		int	level = `0`,
		UT_Vector3I	oddCoarsenings = `UT_Vector3I(0, 0, 0)`
	)

Initialize array to given size and spacing NB: For the 2d case, UT_MultigridArrays MUST have the shape (nx, 1, ny).

template<typename T>

void UT_MultigridArrayT< T >::initFromVoxels	(	const UT_VoxelArray< T > &	f,
		const UT_Vector3T< T > &	spacing,
		const UT_Vector3I &	boundariesNeg,
		const UT_Vector3I &	boundariesPos
	)

Initializes this array from the provided UT_VoxelArray, assuming the provided grid spacing. After this function returns, isInit() will return true. If this array has already been initialized, it will be resized as needed to match the UT_VoxelArray. NB: In the 2d case, UT_MultigridArrays MUST have the shape (nx, 1, ny). This function ensures that will be the case, but if a provided 2D UT_VoxelArray has a different shape, e.g. (nx, ny, 1), then the shapes of the two arrays will not match after this function, although the number of voxels will. copyToVoxels will only resize if the number of voxels differs, so initFromVoxels and copyToVoxels with the same UT_VoxelArray will work fine in the 2D and 3D case.

template<typename T>

void UT_MultigridArrayT< T >::initLaplacian ( )

protected

template<typename T>

void UT_MultigridArrayT< T >::initStorage ( )

protected

template<typename T>

void UT_MultigridArrayT< T >::interpolate	(	const UT_Vector3I &	coarsenAxis,
		const UT_Vector3I &	parity,
		UT_MultigridArrayT< T > &	u
	)		const

Interpolate this array into a fine version of the grid. interpolateAxis is a vector of booleans specifying which axes to interpolate Where interpolateAxis[axis] is true, ensure that u[axis] = getRes(axis) * 2 + parity[axis].

template<typename T>

bool UT_MultigridArrayT< T >::is1D ( ) const

inline

Returns whether this is a 1-D array.

Definition at line 222 of file UT_MultigridArray.h.

template<typename T>

bool UT_MultigridArrayT< T >::isInit ( )

inline

Whether this array has been initialized - default constructor does not.

Definition at line 97 of file UT_MultigridArray.h.

template<typename T>

void UT_MultigridArrayT< T >::match ( const UT_MultigridArrayT< T > & src )

Initializes this to be the same dimensions, boundary conditions, etc, of the given array. The values of this may be reset to zero.

template<typename T>

T UT_MultigridArrayT< T >::norm ( int n ) const

inline

Definition at line 168 of file UT_MultigridArray.h.

template<typename T>

exint UT_MultigridArrayT< T >::numElements ( ) const

inline

Returns the total number of elements in the array.

Definition at line 212 of file UT_MultigridArray.h.

template<typename T>

T UT_MultigridArrayT< T >::operator()	(	int	i,
		int	j,
		int	k
	)		const

inline

Definition at line 103 of file UT_MultigridArray.h.

template<typename T>

T& UT_MultigridArrayT< T >::operator()	(	int	i,
		int	j,
		int	k
	)

inline

Definition at line 108 of file UT_MultigridArray.h.

template<typename T>

T UT_MultigridArrayT< T >::operator() ( const UT_Vector3I & idx ) const

inline

Definition at line 113 of file UT_MultigridArray.h.

template<typename T>

T& UT_MultigridArrayT< T >::operator() ( const UT_Vector3I & idx )

inline

Definition at line 118 of file UT_MultigridArray.h.

template<typename T>

UT_MultigridArrayT<T>& UT_MultigridArrayT< T >::operator+= ( const UT_MultigridArrayT< T > & a )

template<typename T>

UT_MultigridArrayT<T>& UT_MultigridArrayT< T >::operator= ( const UT_MultigridArrayT< T > & a )

template<typename T>

bool UT_MultigridArrayT< T >::shouldMultiThread ( ) const

inline

Definition at line 231 of file UT_MultigridArray.h.

template<typename T>

int UT_MultigridArrayT< T >::solvePoisson	(	fpreal64	abstol,
		fpreal64	reltol,
		int	miniter,
		int	maxiter,
		UT_MultigridArrayT< T > &	x,
		UT_ValArray< fpreal64 > *	resnorminf,
		UT_ValArray< fpreal64 > *	resnorm2,
		bool	startFMG = `true`
	)		const

Iteratively solve the Poisson equation, assuming *this is the right hand side. This function will do at least miniter iterations, then continue to iterate until the 2-norm of the residual has been reduced by reltol, the inf-norm of the residual is below abstol, or maxiters is reached. The optional resnorm0 and resnorm2 parameters will contain the 0- and 2-norms of the residuals for each iteration. If startFMG is true, the first iteration will be full-multigrid; otherwise, only v-cycles are used.

template<typename T>

bool UT_MultigridArrayT< T >::solvePoissonCholesky ( UT_MultigridArrayT< T > & x ) const

protected

Directly solves the Poisson equation for this right-hand side, storing the result in x, using UT_MatrixSolver::choleskyDecomp and UT_MatrixSolver::choleskyBackSub. This should only be called for relatively small N when the coarsest grid level has been reached, and only when the result of buildMatrix is positive (semi-definite)

template<typename T>

bool UT_MultigridArrayT< T >::solvePoissonLU ( UT_MultigridArrayT< T > & x ) const

protected

Directly solves the Poisson equation for this right-hand side, using LU decomposition. This should only be called for relatively small N when the coarsest grid level has been reached, and only when the result of buildMatrix is non-singular.

template<typename T>

bool UT_MultigridArrayT< T >::solvePoissonSVD ( UT_MultigridArrayT< T > & x ) const

protected

Directly solves the Poisson equation for this right-hand side, using SVD decomposition. This can be used when the result of buildMatrix is singular. This should only be called for relatively small N when the coarsest grid level has been reached.

template<typename T>

UT_MultigridArrayT< T >::THREADED_METHOD2_CONST	(	UT_MultigridArrayT< T >	,
		shouldMultiThread()	,
		subtractApplyLaplacian	,
		const UT_MultigridArrayT< T > &	,
		x	,
		UT_MultigridArrayT< T > &	,
		r
	)		const

Computes rhe residual r = b - A * x, where b == *this.

template<typename T>

UT_MultigridArrayT< T >::THREADED_METHOD2_CONST	(	UT_MultigridArrayT< T >	,
		shouldMultiThread()	,
		smoothLaplacianGaussSeidel	,
		int	,
		parity	,
		UT_MultigridArrayT< T > &	,
		x
	)

Performs red-black Gauss-Seidel smoothing according to parity, where b == *this, and the implicit matrix is the Laplacian operator. This should be called twice, with parity=0 and 1, for a full smoothing cycle.

template<typename T>

UT_MultigridArrayT< T >::THREADED_METHOD2_CONST	(	UT_MultigridArrayT< T >	,
		shouldMultiThread()	,
		coarsenAlongAxis	,
		UT_MultigridArrayT< T > &	,
		uh	,
		int	,
		axis
	)

Performs coarsening using full-weighting along one axis. On output uh.shape[axis] == this->shape[axis] / 2

template<typename T>

UT_MultigridArrayT< T >::THREADED_METHOD2_CONST	(	UT_MultigridArrayT< T >	,
		shouldMultiThread()	,
		interpolateAlongAxis	,
		UT_MultigridArrayT< T > &	,
		u	,
		int	,
		axis
	)

Interpolate this array into u using inverse of full-weighting along one axis. On output u.shape[axis] == this->shape[axis] * 2 + parity.

template<typename T>

UT_MultigridArrayT< T >::THREADED_METHOD2_CONST	(	UT_MultigridArrayT< T >	,
		shouldMultiThread()	,
		computeNormsInternal	,
		fpreal64 *	,
		norminf	,
		fpreal64 *	,
		norm2
	)

protected

template<typename T>

UT_MultigridArrayT< T >::THREADED_METHOD3_CONST	(	UT_MultigridArrayT< T >	,
		shouldMultiThread()	,
		computeResidualNormsInternal	,
		const UT_MultigridArrayT< T > &	,
		x	,
		fpreal64 *	,
		norminf	,
		fpreal64 *	,
		norm2
	)		const

protected

template<typename T>

void UT_MultigridArrayT< T >::vcycle	(	exint	minPerAxis,
		int	nSmoothDown,
		int	nSmoothUp,
		UT_MultigridArrayT< T > &	x,
		bool	smoothTopLevelDown = `true`
	)		const

Perform recursive V-cycle for multigrid algorithm, until reaching numElements() < maxcoarse, at which point a direct solve is done using directSolve.. At each level, nSmoothDown and nSmoothUp iterations of Gauss-Seidel smoothing are done, unless smoothTopLevelDown is false, in which case no smoothing is done for the top-level grid on the down-cycle. This can be used as an optimization when calling vcycle iteratively; the previous vcycle call will have smoothed the top-level grid as its last step, so you can often skip smoothing that same grid on the next vcycle without loss of convergence.

template<typename T>

void UT_MultigridArrayT< T >::zero ( )

inline

Definition at line 163 of file UT_MultigridArray.h.