forked from bartvdbraak/blender
810 lines
24 KiB
C++
810 lines
24 KiB
C++
|
|
|
|
// DO NOT EDIT !
|
|
// This file is generated using the MantaFlow preprocessor (prep generate).
|
|
|
|
/******************************************************************************
|
|
*
|
|
* MantaFlow fluid solver framework
|
|
* Copyright 2011 Tobias Pfaff, Nils Thuerey
|
|
*
|
|
* This program is free software, distributed under the terms of the
|
|
* Apache License, Version 2.0
|
|
* http://www.apache.org/licenses/LICENSE-2.0
|
|
*
|
|
* Conjugate gradient solver
|
|
*
|
|
******************************************************************************/
|
|
|
|
#ifndef _CONJUGATEGRADIENT_H
|
|
#define _CONJUGATEGRADIENT_H
|
|
|
|
#include "vectorbase.h"
|
|
#include "grid.h"
|
|
#include "kernel.h"
|
|
#include "multigrid.h"
|
|
|
|
namespace Manta {
|
|
|
|
static const bool CG_DEBUG = false;
|
|
|
|
//! Basic CG interface
|
|
class GridCgInterface {
|
|
public:
|
|
enum PreconditionType { PC_None = 0, PC_ICP, PC_mICP, PC_MGP };
|
|
|
|
GridCgInterface() : mUseL2Norm(true){};
|
|
virtual ~GridCgInterface(){};
|
|
|
|
// solving functions
|
|
virtual bool iterate() = 0;
|
|
virtual void solve(int maxIter) = 0;
|
|
|
|
// precond
|
|
virtual void setICPreconditioner(
|
|
PreconditionType method, Grid<Real> *A0, Grid<Real> *Ai, Grid<Real> *Aj, Grid<Real> *Ak) = 0;
|
|
virtual void setMGPreconditioner(PreconditionType method, GridMg *MG) = 0;
|
|
|
|
// access
|
|
virtual Real getSigma() const = 0;
|
|
virtual Real getIterations() const = 0;
|
|
virtual Real getResNorm() const = 0;
|
|
virtual void setAccuracy(Real set) = 0;
|
|
virtual Real getAccuracy() const = 0;
|
|
|
|
//! force reinit upon next iterate() call, can be used for doing multiple solves
|
|
virtual void forceReinit() = 0;
|
|
|
|
void setUseL2Norm(bool set)
|
|
{
|
|
mUseL2Norm = set;
|
|
}
|
|
|
|
protected:
|
|
// use l2 norm of residualfor threshold? (otherwise uses max norm)
|
|
bool mUseL2Norm;
|
|
};
|
|
|
|
//! Run single iteration of the cg solver
|
|
/*! the template argument determines the type of matrix multiplication,
|
|
typically a ApplyMatrix kernel, another one is needed e.g. for the
|
|
mesh-based wave equation solver */
|
|
template<class APPLYMAT> class GridCg : public GridCgInterface {
|
|
public:
|
|
//! constructor
|
|
GridCg(Grid<Real> &dst,
|
|
Grid<Real> &rhs,
|
|
Grid<Real> &residual,
|
|
Grid<Real> &search,
|
|
const FlagGrid &flags,
|
|
Grid<Real> &tmp,
|
|
std::vector<Grid<Real> *> matrixAVec,
|
|
std::vector<Grid<Real> *> rhsVec = {});
|
|
~GridCg(){};
|
|
|
|
void doInit();
|
|
bool iterate();
|
|
void solve(int maxIter);
|
|
//! init pointers, and copy values from "normal" matrix
|
|
void setICPreconditioner(
|
|
PreconditionType method, Grid<Real> *A0, Grid<Real> *Ai, Grid<Real> *Aj, Grid<Real> *Ak);
|
|
void setMGPreconditioner(PreconditionType method, GridMg *MG);
|
|
void forceReinit()
|
|
{
|
|
mInited = false;
|
|
}
|
|
|
|
// Accessors
|
|
Real getSigma() const
|
|
{
|
|
return mSigma;
|
|
}
|
|
Real getIterations() const
|
|
{
|
|
return mIterations;
|
|
}
|
|
|
|
Real getResNorm() const
|
|
{
|
|
return mResNorm;
|
|
}
|
|
|
|
void setAccuracy(Real set)
|
|
{
|
|
mAccuracy = set;
|
|
}
|
|
Real getAccuracy() const
|
|
{
|
|
return mAccuracy;
|
|
}
|
|
|
|
protected:
|
|
bool mInited;
|
|
int mIterations;
|
|
// grids
|
|
Grid<Real> &mDst;
|
|
Grid<Real> &mRhs;
|
|
Grid<Real> &mResidual;
|
|
Grid<Real> &mSearch;
|
|
const FlagGrid &mFlags;
|
|
Grid<Real> &mTmp;
|
|
|
|
//! shape of A matrix defined here (e.g. diagonal, positive neighbor cells, etc)
|
|
std::vector<Grid<Real> *> mMatrixA;
|
|
//! shape of rhs vector defined here (e.g. 1 rhs for regular fluids solve, 3 rhs for viscosity)
|
|
std::vector<Grid<Real> *> mVecRhs;
|
|
|
|
PreconditionType mPcMethod;
|
|
//! preconditioning grids
|
|
Grid<Real> *mpPCA0, *mpPCAi, *mpPCAj, *mpPCAk;
|
|
GridMg *mMG;
|
|
|
|
//! sigma / residual
|
|
Real mSigma;
|
|
//! accuracy of solver (max. residuum)
|
|
Real mAccuracy;
|
|
//! norm of the residual
|
|
Real mResNorm;
|
|
}; // GridCg
|
|
|
|
//! Kernel: Apply symmetric stored Matrix
|
|
|
|
struct ApplyMatrix : public KernelBase {
|
|
ApplyMatrix(const FlagGrid &flags,
|
|
Grid<Real> &dst,
|
|
const Grid<Real> &src,
|
|
const std::vector<Grid<Real> *> matrixA,
|
|
const std::vector<Grid<Real> *> vecRhs)
|
|
: KernelBase(&flags, 0), flags(flags), dst(dst), src(src), matrixA(matrixA), vecRhs(vecRhs)
|
|
{
|
|
runMessage();
|
|
run();
|
|
}
|
|
inline void op(IndexInt idx,
|
|
const FlagGrid &flags,
|
|
Grid<Real> &dst,
|
|
const Grid<Real> &src,
|
|
const std::vector<Grid<Real> *> matrixA,
|
|
const std::vector<Grid<Real> *> vecRhs) const
|
|
{
|
|
unusedParameter(vecRhs); // Not needed in this matrix application
|
|
|
|
if (matrixA.size() != 4)
|
|
errMsg("ConjugateGrad: Invalid A matrix in apply matrix step");
|
|
Grid<Real> &A0 = *matrixA[0];
|
|
Grid<Real> &Ai = *matrixA[1];
|
|
Grid<Real> &Aj = *matrixA[2];
|
|
Grid<Real> &Ak = *matrixA[3];
|
|
|
|
if (!flags.isFluid(idx)) {
|
|
dst[idx] = src[idx];
|
|
return;
|
|
}
|
|
|
|
dst[idx] = src[idx] * A0[idx] + src[idx - X] * Ai[idx - X] + src[idx + X] * Ai[idx] +
|
|
src[idx - Y] * Aj[idx - Y] + src[idx + Y] * Aj[idx] + src[idx - Z] * Ak[idx - Z] +
|
|
src[idx + Z] * Ak[idx];
|
|
}
|
|
inline const FlagGrid &getArg0()
|
|
{
|
|
return flags;
|
|
}
|
|
typedef FlagGrid type0;
|
|
inline Grid<Real> &getArg1()
|
|
{
|
|
return dst;
|
|
}
|
|
typedef Grid<Real> type1;
|
|
inline const Grid<Real> &getArg2()
|
|
{
|
|
return src;
|
|
}
|
|
typedef Grid<Real> type2;
|
|
inline const std::vector<Grid<Real> *> &getArg3()
|
|
{
|
|
return matrixA;
|
|
}
|
|
typedef std::vector<Grid<Real> *> type3;
|
|
inline const std::vector<Grid<Real> *> &getArg4()
|
|
{
|
|
return vecRhs;
|
|
}
|
|
typedef std::vector<Grid<Real> *> type4;
|
|
void runMessage()
|
|
{
|
|
debMsg("Executing kernel ApplyMatrix ", 3);
|
|
debMsg("Kernel range"
|
|
<< " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
|
|
4);
|
|
};
|
|
void operator()(const tbb::blocked_range<IndexInt> &__r) const
|
|
{
|
|
for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
|
|
op(idx, flags, dst, src, matrixA, vecRhs);
|
|
}
|
|
void run()
|
|
{
|
|
tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
|
|
}
|
|
const FlagGrid &flags;
|
|
Grid<Real> &dst;
|
|
const Grid<Real> &src;
|
|
const std::vector<Grid<Real> *> matrixA;
|
|
const std::vector<Grid<Real> *> vecRhs;
|
|
};
|
|
|
|
//! Kernel: Apply symmetric stored Matrix. 2D version
|
|
|
|
struct ApplyMatrix2D : public KernelBase {
|
|
ApplyMatrix2D(const FlagGrid &flags,
|
|
Grid<Real> &dst,
|
|
const Grid<Real> &src,
|
|
const std::vector<Grid<Real> *> matrixA,
|
|
const std::vector<Grid<Real> *> vecRhs)
|
|
: KernelBase(&flags, 0), flags(flags), dst(dst), src(src), matrixA(matrixA), vecRhs(vecRhs)
|
|
{
|
|
runMessage();
|
|
run();
|
|
}
|
|
inline void op(IndexInt idx,
|
|
const FlagGrid &flags,
|
|
Grid<Real> &dst,
|
|
const Grid<Real> &src,
|
|
const std::vector<Grid<Real> *> matrixA,
|
|
const std::vector<Grid<Real> *> vecRhs) const
|
|
{
|
|
unusedParameter(vecRhs); // Not needed in this matrix application
|
|
|
|
if (matrixA.size() != 3)
|
|
errMsg("ConjugateGrad: Invalid A matrix in apply matrix step");
|
|
Grid<Real> &A0 = *matrixA[0];
|
|
Grid<Real> &Ai = *matrixA[1];
|
|
Grid<Real> &Aj = *matrixA[2];
|
|
|
|
if (!flags.isFluid(idx)) {
|
|
dst[idx] = src[idx];
|
|
return;
|
|
}
|
|
|
|
dst[idx] = src[idx] * A0[idx] + src[idx - X] * Ai[idx - X] + src[idx + X] * Ai[idx] +
|
|
src[idx - Y] * Aj[idx - Y] + src[idx + Y] * Aj[idx];
|
|
}
|
|
inline const FlagGrid &getArg0()
|
|
{
|
|
return flags;
|
|
}
|
|
typedef FlagGrid type0;
|
|
inline Grid<Real> &getArg1()
|
|
{
|
|
return dst;
|
|
}
|
|
typedef Grid<Real> type1;
|
|
inline const Grid<Real> &getArg2()
|
|
{
|
|
return src;
|
|
}
|
|
typedef Grid<Real> type2;
|
|
inline const std::vector<Grid<Real> *> &getArg3()
|
|
{
|
|
return matrixA;
|
|
}
|
|
typedef std::vector<Grid<Real> *> type3;
|
|
inline const std::vector<Grid<Real> *> &getArg4()
|
|
{
|
|
return vecRhs;
|
|
}
|
|
typedef std::vector<Grid<Real> *> type4;
|
|
void runMessage()
|
|
{
|
|
debMsg("Executing kernel ApplyMatrix2D ", 3);
|
|
debMsg("Kernel range"
|
|
<< " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
|
|
4);
|
|
};
|
|
void operator()(const tbb::blocked_range<IndexInt> &__r) const
|
|
{
|
|
for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
|
|
op(idx, flags, dst, src, matrixA, vecRhs);
|
|
}
|
|
void run()
|
|
{
|
|
tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
|
|
}
|
|
const FlagGrid &flags;
|
|
Grid<Real> &dst;
|
|
const Grid<Real> &src;
|
|
const std::vector<Grid<Real> *> matrixA;
|
|
const std::vector<Grid<Real> *> vecRhs;
|
|
};
|
|
|
|
struct ApplyMatrixViscosityU : public KernelBase {
|
|
ApplyMatrixViscosityU(const FlagGrid &flags,
|
|
Grid<Real> &dst,
|
|
const Grid<Real> &src,
|
|
const std::vector<Grid<Real> *> matrixA,
|
|
const std::vector<Grid<Real> *> vecRhs)
|
|
: KernelBase(&flags, 1), flags(flags), dst(dst), src(src), matrixA(matrixA), vecRhs(vecRhs)
|
|
{
|
|
runMessage();
|
|
run();
|
|
}
|
|
inline void op(int i,
|
|
int j,
|
|
int k,
|
|
const FlagGrid &flags,
|
|
Grid<Real> &dst,
|
|
const Grid<Real> &src,
|
|
const std::vector<Grid<Real> *> matrixA,
|
|
const std::vector<Grid<Real> *> vecRhs) const
|
|
{
|
|
if (matrixA.size() != 15)
|
|
errMsg("ConjugateGrad: Invalid A matrix in apply matrix step");
|
|
Grid<Real> &A0 = *matrixA[0];
|
|
Grid<Real> &Aplusi = *matrixA[1];
|
|
Grid<Real> &Aplusj = *matrixA[2];
|
|
Grid<Real> &Aplusk = *matrixA[3];
|
|
Grid<Real> &Aminusi = *matrixA[4];
|
|
Grid<Real> &Aminusj = *matrixA[5];
|
|
Grid<Real> &Aminusk = *matrixA[6];
|
|
|
|
if (vecRhs.size() != 2)
|
|
errMsg("ConjugateGrad: Invalid rhs vector in apply matrix step");
|
|
Grid<Real> &srcV = *vecRhs[0];
|
|
Grid<Real> &srcW = *vecRhs[1];
|
|
|
|
dst(i, j, k) = src(i, j, k) * A0(i, j, k) + src(i + 1, j, k) * Aplusi(i, j, k) +
|
|
src(i, j + 1, k) * Aplusj(i, j, k) + src(i, j, k + 1) * Aplusk(i, j, k) +
|
|
src(i - 1, j, k) * Aminusi(i, j, k) + src(i, j - 1, k) * Aminusj(i, j, k) +
|
|
src(i, j, k - 1) * Aminusk(i, j, k);
|
|
|
|
dst(i, j, k) += srcV(i, j + 1, k) * (*matrixA[7])(i, j, k) +
|
|
srcV(i - 1, j + 1, k) * (*matrixA[8])(i, j, k) +
|
|
srcV(i, j, k) * (*matrixA[9])(i, j, k) +
|
|
srcV(i - 1, j, k) * (*matrixA[10])(i, j, k) +
|
|
srcW(i, j, k + 1) * (*matrixA[11])(i, j, k) +
|
|
srcW(i - 1, j, k + 1) * (*matrixA[12])(i, j, k) +
|
|
srcW(i, j, k) * (*matrixA[13])(i, j, k) +
|
|
srcW(i - 1, j, k) * (*matrixA[14])(i, j, k);
|
|
}
|
|
inline const FlagGrid &getArg0()
|
|
{
|
|
return flags;
|
|
}
|
|
typedef FlagGrid type0;
|
|
inline Grid<Real> &getArg1()
|
|
{
|
|
return dst;
|
|
}
|
|
typedef Grid<Real> type1;
|
|
inline const Grid<Real> &getArg2()
|
|
{
|
|
return src;
|
|
}
|
|
typedef Grid<Real> type2;
|
|
inline const std::vector<Grid<Real> *> &getArg3()
|
|
{
|
|
return matrixA;
|
|
}
|
|
typedef std::vector<Grid<Real> *> type3;
|
|
inline const std::vector<Grid<Real> *> &getArg4()
|
|
{
|
|
return vecRhs;
|
|
}
|
|
typedef std::vector<Grid<Real> *> type4;
|
|
void runMessage()
|
|
{
|
|
debMsg("Executing kernel ApplyMatrixViscosityU ", 3);
|
|
debMsg("Kernel range"
|
|
<< " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
|
|
4);
|
|
};
|
|
void operator()(const tbb::blocked_range<IndexInt> &__r) const
|
|
{
|
|
const int _maxX = maxX;
|
|
const int _maxY = maxY;
|
|
if (maxZ > 1) {
|
|
for (int k = __r.begin(); k != (int)__r.end(); k++)
|
|
for (int j = 1; j < _maxY; j++)
|
|
for (int i = 1; i < _maxX; i++)
|
|
op(i, j, k, flags, dst, src, matrixA, vecRhs);
|
|
}
|
|
else {
|
|
const int k = 0;
|
|
for (int j = __r.begin(); j != (int)__r.end(); j++)
|
|
for (int i = 1; i < _maxX; i++)
|
|
op(i, j, k, flags, dst, src, matrixA, vecRhs);
|
|
}
|
|
}
|
|
void run()
|
|
{
|
|
if (maxZ > 1)
|
|
tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
|
|
else
|
|
tbb::parallel_for(tbb::blocked_range<IndexInt>(1, maxY), *this);
|
|
}
|
|
const FlagGrid &flags;
|
|
Grid<Real> &dst;
|
|
const Grid<Real> &src;
|
|
const std::vector<Grid<Real> *> matrixA;
|
|
const std::vector<Grid<Real> *> vecRhs;
|
|
};
|
|
|
|
struct ApplyMatrixViscosityV : public KernelBase {
|
|
ApplyMatrixViscosityV(const FlagGrid &flags,
|
|
Grid<Real> &dst,
|
|
const Grid<Real> &src,
|
|
const std::vector<Grid<Real> *> matrixA,
|
|
const std::vector<Grid<Real> *> vecRhs)
|
|
: KernelBase(&flags, 1), flags(flags), dst(dst), src(src), matrixA(matrixA), vecRhs(vecRhs)
|
|
{
|
|
runMessage();
|
|
run();
|
|
}
|
|
inline void op(int i,
|
|
int j,
|
|
int k,
|
|
const FlagGrid &flags,
|
|
Grid<Real> &dst,
|
|
const Grid<Real> &src,
|
|
const std::vector<Grid<Real> *> matrixA,
|
|
const std::vector<Grid<Real> *> vecRhs) const
|
|
{
|
|
if (matrixA.size() != 15)
|
|
errMsg("ConjugateGrad: Invalid A matrix in apply matrix step");
|
|
Grid<Real> &A0 = *matrixA[0];
|
|
Grid<Real> &Aplusi = *matrixA[1];
|
|
Grid<Real> &Aplusj = *matrixA[2];
|
|
Grid<Real> &Aplusk = *matrixA[3];
|
|
Grid<Real> &Aminusi = *matrixA[4];
|
|
Grid<Real> &Aminusj = *matrixA[5];
|
|
Grid<Real> &Aminusk = *matrixA[6];
|
|
|
|
if (vecRhs.size() != 2)
|
|
errMsg("ConjugateGrad: Invalid rhs vector in apply matrix step");
|
|
Grid<Real> &srcU = *vecRhs[0];
|
|
Grid<Real> &srcW = *vecRhs[1];
|
|
|
|
dst(i, j, k) = src(i, j, k) * A0(i, j, k) + src(i + 1, j, k) * Aplusi(i, j, k) +
|
|
src(i, j + 1, k) * Aplusj(i, j, k) + src(i, j, k + 1) * Aplusk(i, j, k) +
|
|
src(i - 1, j, k) * Aminusi(i, j, k) + src(i, j - 1, k) * Aminusj(i, j, k) +
|
|
src(i, j, k - 1) * Aminusk(i, j, k);
|
|
|
|
dst(i, j, k) += srcU(i + 1, j, k) * (*matrixA[7])(i, j, k) +
|
|
srcU(i + 1, j - 1, k) * (*matrixA[8])(i, j, k) +
|
|
srcU(i, j, k) * (*matrixA[9])(i, j, k) +
|
|
srcU(i, j - 1, k) * (*matrixA[10])(i, j, k) +
|
|
srcW(i, j, k + 1) * (*matrixA[11])(i, j, k) +
|
|
srcW(i, j - 1, k + 1) * (*matrixA[12])(i, j, k) +
|
|
srcW(i, j, k) * (*matrixA[13])(i, j, k) +
|
|
srcW(i, j - 1, k) * (*matrixA[14])(i, j, k);
|
|
}
|
|
inline const FlagGrid &getArg0()
|
|
{
|
|
return flags;
|
|
}
|
|
typedef FlagGrid type0;
|
|
inline Grid<Real> &getArg1()
|
|
{
|
|
return dst;
|
|
}
|
|
typedef Grid<Real> type1;
|
|
inline const Grid<Real> &getArg2()
|
|
{
|
|
return src;
|
|
}
|
|
typedef Grid<Real> type2;
|
|
inline const std::vector<Grid<Real> *> &getArg3()
|
|
{
|
|
return matrixA;
|
|
}
|
|
typedef std::vector<Grid<Real> *> type3;
|
|
inline const std::vector<Grid<Real> *> &getArg4()
|
|
{
|
|
return vecRhs;
|
|
}
|
|
typedef std::vector<Grid<Real> *> type4;
|
|
void runMessage()
|
|
{
|
|
debMsg("Executing kernel ApplyMatrixViscosityV ", 3);
|
|
debMsg("Kernel range"
|
|
<< " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
|
|
4);
|
|
};
|
|
void operator()(const tbb::blocked_range<IndexInt> &__r) const
|
|
{
|
|
const int _maxX = maxX;
|
|
const int _maxY = maxY;
|
|
if (maxZ > 1) {
|
|
for (int k = __r.begin(); k != (int)__r.end(); k++)
|
|
for (int j = 1; j < _maxY; j++)
|
|
for (int i = 1; i < _maxX; i++)
|
|
op(i, j, k, flags, dst, src, matrixA, vecRhs);
|
|
}
|
|
else {
|
|
const int k = 0;
|
|
for (int j = __r.begin(); j != (int)__r.end(); j++)
|
|
for (int i = 1; i < _maxX; i++)
|
|
op(i, j, k, flags, dst, src, matrixA, vecRhs);
|
|
}
|
|
}
|
|
void run()
|
|
{
|
|
if (maxZ > 1)
|
|
tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
|
|
else
|
|
tbb::parallel_for(tbb::blocked_range<IndexInt>(1, maxY), *this);
|
|
}
|
|
const FlagGrid &flags;
|
|
Grid<Real> &dst;
|
|
const Grid<Real> &src;
|
|
const std::vector<Grid<Real> *> matrixA;
|
|
const std::vector<Grid<Real> *> vecRhs;
|
|
};
|
|
|
|
struct ApplyMatrixViscosityW : public KernelBase {
|
|
ApplyMatrixViscosityW(const FlagGrid &flags,
|
|
Grid<Real> &dst,
|
|
const Grid<Real> &src,
|
|
const std::vector<Grid<Real> *> matrixA,
|
|
const std::vector<Grid<Real> *> vecRhs)
|
|
: KernelBase(&flags, 1), flags(flags), dst(dst), src(src), matrixA(matrixA), vecRhs(vecRhs)
|
|
{
|
|
runMessage();
|
|
run();
|
|
}
|
|
inline void op(int i,
|
|
int j,
|
|
int k,
|
|
const FlagGrid &flags,
|
|
Grid<Real> &dst,
|
|
const Grid<Real> &src,
|
|
const std::vector<Grid<Real> *> matrixA,
|
|
const std::vector<Grid<Real> *> vecRhs) const
|
|
{
|
|
if (matrixA.size() != 15)
|
|
errMsg("ConjugateGrad: Invalid A matrix in apply matrix step");
|
|
Grid<Real> &A0 = *matrixA[0];
|
|
Grid<Real> &Aplusi = *matrixA[1];
|
|
Grid<Real> &Aplusj = *matrixA[2];
|
|
Grid<Real> &Aplusk = *matrixA[3];
|
|
Grid<Real> &Aminusi = *matrixA[4];
|
|
Grid<Real> &Aminusj = *matrixA[5];
|
|
Grid<Real> &Aminusk = *matrixA[6];
|
|
|
|
if (vecRhs.size() != 2)
|
|
errMsg("ConjugateGrad: Invalid rhs vector in apply matrix step");
|
|
Grid<Real> &srcU = *vecRhs[0];
|
|
Grid<Real> &srcV = *vecRhs[1];
|
|
|
|
dst(i, j, k) = src(i, j, k) * A0(i, j, k) + src(i + 1, j, k) * Aplusi(i, j, k) +
|
|
src(i, j + 1, k) * Aplusj(i, j, k) + src(i, j, k + 1) * Aplusk(i, j, k) +
|
|
src(i - 1, j, k) * Aminusi(i, j, k) + src(i, j - 1, k) * Aminusj(i, j, k) +
|
|
src(i, j, k - 1) * Aminusk(i, j, k);
|
|
|
|
dst(i, j, k) += srcU(i + 1, j, k) * (*matrixA[7])(i, j, k) +
|
|
srcU(i + 1, j, k - 1) * (*matrixA[8])(i, j, k) +
|
|
srcU(i, j, k) * (*matrixA[9])(i, j, k) +
|
|
srcU(i, j, k - 1) * (*matrixA[10])(i, j, k) +
|
|
srcV(i, j + 1, k) * (*matrixA[11])(i, j, k) +
|
|
srcV(i, j + 1, k - 1) * (*matrixA[12])(i, j, k) +
|
|
srcV(i, j, k) * (*matrixA[13])(i, j, k) +
|
|
srcV(i, j, k - 1) * (*matrixA[14])(i, j, k);
|
|
}
|
|
inline const FlagGrid &getArg0()
|
|
{
|
|
return flags;
|
|
}
|
|
typedef FlagGrid type0;
|
|
inline Grid<Real> &getArg1()
|
|
{
|
|
return dst;
|
|
}
|
|
typedef Grid<Real> type1;
|
|
inline const Grid<Real> &getArg2()
|
|
{
|
|
return src;
|
|
}
|
|
typedef Grid<Real> type2;
|
|
inline const std::vector<Grid<Real> *> &getArg3()
|
|
{
|
|
return matrixA;
|
|
}
|
|
typedef std::vector<Grid<Real> *> type3;
|
|
inline const std::vector<Grid<Real> *> &getArg4()
|
|
{
|
|
return vecRhs;
|
|
}
|
|
typedef std::vector<Grid<Real> *> type4;
|
|
void runMessage()
|
|
{
|
|
debMsg("Executing kernel ApplyMatrixViscosityW ", 3);
|
|
debMsg("Kernel range"
|
|
<< " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
|
|
4);
|
|
};
|
|
void operator()(const tbb::blocked_range<IndexInt> &__r) const
|
|
{
|
|
const int _maxX = maxX;
|
|
const int _maxY = maxY;
|
|
if (maxZ > 1) {
|
|
for (int k = __r.begin(); k != (int)__r.end(); k++)
|
|
for (int j = 1; j < _maxY; j++)
|
|
for (int i = 1; i < _maxX; i++)
|
|
op(i, j, k, flags, dst, src, matrixA, vecRhs);
|
|
}
|
|
else {
|
|
const int k = 0;
|
|
for (int j = __r.begin(); j != (int)__r.end(); j++)
|
|
for (int i = 1; i < _maxX; i++)
|
|
op(i, j, k, flags, dst, src, matrixA, vecRhs);
|
|
}
|
|
}
|
|
void run()
|
|
{
|
|
if (maxZ > 1)
|
|
tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
|
|
else
|
|
tbb::parallel_for(tbb::blocked_range<IndexInt>(1, maxY), *this);
|
|
}
|
|
const FlagGrid &flags;
|
|
Grid<Real> &dst;
|
|
const Grid<Real> &src;
|
|
const std::vector<Grid<Real> *> matrixA;
|
|
const std::vector<Grid<Real> *> vecRhs;
|
|
};
|
|
|
|
/* NOTE: Use this template for new matrix application kernels
|
|
|
|
//! Template for matrix application kernels
|
|
KERNEL()
|
|
void ApplyMatrixTemplate (const FlagGrid& flags, Grid<Real>& dst, const Grid<Real>& src,
|
|
const std::vector<Grid<Real> *> matrixA, const std::vector<Grid<Real> *> vecRhs)
|
|
{
|
|
// The kernel must define how to use the grids from the matrixA and vecRhs lists
|
|
}
|
|
|
|
*/
|
|
|
|
//! Kernel: Construct the matrix for the poisson equation
|
|
|
|
struct MakeLaplaceMatrix : public KernelBase {
|
|
MakeLaplaceMatrix(const FlagGrid &flags,
|
|
Grid<Real> &A0,
|
|
Grid<Real> &Ai,
|
|
Grid<Real> &Aj,
|
|
Grid<Real> &Ak,
|
|
const MACGrid *fractions = 0)
|
|
: KernelBase(&flags, 1), flags(flags), A0(A0), Ai(Ai), Aj(Aj), Ak(Ak), fractions(fractions)
|
|
{
|
|
runMessage();
|
|
run();
|
|
}
|
|
inline void op(int i,
|
|
int j,
|
|
int k,
|
|
const FlagGrid &flags,
|
|
Grid<Real> &A0,
|
|
Grid<Real> &Ai,
|
|
Grid<Real> &Aj,
|
|
Grid<Real> &Ak,
|
|
const MACGrid *fractions = 0) const
|
|
{
|
|
if (!flags.isFluid(i, j, k))
|
|
return;
|
|
|
|
if (!fractions) {
|
|
// diagonal, A0
|
|
if (!flags.isObstacle(i - 1, j, k))
|
|
A0(i, j, k) += 1.;
|
|
if (!flags.isObstacle(i + 1, j, k))
|
|
A0(i, j, k) += 1.;
|
|
if (!flags.isObstacle(i, j - 1, k))
|
|
A0(i, j, k) += 1.;
|
|
if (!flags.isObstacle(i, j + 1, k))
|
|
A0(i, j, k) += 1.;
|
|
if (flags.is3D() && !flags.isObstacle(i, j, k - 1))
|
|
A0(i, j, k) += 1.;
|
|
if (flags.is3D() && !flags.isObstacle(i, j, k + 1))
|
|
A0(i, j, k) += 1.;
|
|
|
|
// off-diagonal entries
|
|
if (flags.isFluid(i + 1, j, k))
|
|
Ai(i, j, k) = -1.;
|
|
if (flags.isFluid(i, j + 1, k))
|
|
Aj(i, j, k) = -1.;
|
|
if (flags.is3D() && flags.isFluid(i, j, k + 1))
|
|
Ak(i, j, k) = -1.;
|
|
}
|
|
else {
|
|
// diagonal
|
|
A0(i, j, k) += fractions->get(i, j, k).x;
|
|
A0(i, j, k) += fractions->get(i + 1, j, k).x;
|
|
A0(i, j, k) += fractions->get(i, j, k).y;
|
|
A0(i, j, k) += fractions->get(i, j + 1, k).y;
|
|
if (flags.is3D())
|
|
A0(i, j, k) += fractions->get(i, j, k).z;
|
|
if (flags.is3D())
|
|
A0(i, j, k) += fractions->get(i, j, k + 1).z;
|
|
|
|
// off-diagonal entries
|
|
if (flags.isFluid(i + 1, j, k))
|
|
Ai(i, j, k) = -fractions->get(i + 1, j, k).x;
|
|
if (flags.isFluid(i, j + 1, k))
|
|
Aj(i, j, k) = -fractions->get(i, j + 1, k).y;
|
|
if (flags.is3D() && flags.isFluid(i, j, k + 1))
|
|
Ak(i, j, k) = -fractions->get(i, j, k + 1).z;
|
|
}
|
|
}
|
|
inline const FlagGrid &getArg0()
|
|
{
|
|
return flags;
|
|
}
|
|
typedef FlagGrid type0;
|
|
inline Grid<Real> &getArg1()
|
|
{
|
|
return A0;
|
|
}
|
|
typedef Grid<Real> type1;
|
|
inline Grid<Real> &getArg2()
|
|
{
|
|
return Ai;
|
|
}
|
|
typedef Grid<Real> type2;
|
|
inline Grid<Real> &getArg3()
|
|
{
|
|
return Aj;
|
|
}
|
|
typedef Grid<Real> type3;
|
|
inline Grid<Real> &getArg4()
|
|
{
|
|
return Ak;
|
|
}
|
|
typedef Grid<Real> type4;
|
|
inline const MACGrid *getArg5()
|
|
{
|
|
return fractions;
|
|
}
|
|
typedef MACGrid type5;
|
|
void runMessage()
|
|
{
|
|
debMsg("Executing kernel MakeLaplaceMatrix ", 3);
|
|
debMsg("Kernel range"
|
|
<< " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
|
|
4);
|
|
};
|
|
void operator()(const tbb::blocked_range<IndexInt> &__r) const
|
|
{
|
|
const int _maxX = maxX;
|
|
const int _maxY = maxY;
|
|
if (maxZ > 1) {
|
|
for (int k = __r.begin(); k != (int)__r.end(); k++)
|
|
for (int j = 1; j < _maxY; j++)
|
|
for (int i = 1; i < _maxX; i++)
|
|
op(i, j, k, flags, A0, Ai, Aj, Ak, fractions);
|
|
}
|
|
else {
|
|
const int k = 0;
|
|
for (int j = __r.begin(); j != (int)__r.end(); j++)
|
|
for (int i = 1; i < _maxX; i++)
|
|
op(i, j, k, flags, A0, Ai, Aj, Ak, fractions);
|
|
}
|
|
}
|
|
void run()
|
|
{
|
|
if (maxZ > 1)
|
|
tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
|
|
else
|
|
tbb::parallel_for(tbb::blocked_range<IndexInt>(1, maxY), *this);
|
|
}
|
|
const FlagGrid &flags;
|
|
Grid<Real> &A0;
|
|
Grid<Real> &Ai;
|
|
Grid<Real> &Aj;
|
|
Grid<Real> &Ak;
|
|
const MACGrid *fractions;
|
|
};
|
|
|
|
} // namespace Manta
|
|
|
|
#endif
|