1429 lines
43 KiB
C++
1429 lines
43 KiB
C++
|
|
|
|
// DO NOT EDIT !
|
|
// This file is generated using the MantaFlow preprocessor (prep generate).
|
|
|
|
/******************************************************************************
|
|
*
|
|
* MantaFlow fluid solver framework
|
|
* Copyright 2020 Sebastian Barschkis, 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
|
|
*
|
|
* Accurate Viscous Free Surfaces for Buckling, Coiling, and Rotating Liquids
|
|
* Batty et al., SCA 2008
|
|
*
|
|
******************************************************************************/
|
|
|
|
#include "conjugategrad.h"
|
|
#include "general.h"
|
|
#include "grid.h"
|
|
#include "vectorbase.h"
|
|
|
|
#include <chrono>
|
|
|
|
#if OPENMP == 1 || TBB == 1
|
|
# define ENABLE_PARALLEL 0
|
|
#endif
|
|
|
|
#if ENABLE_PARALLEL == 1
|
|
# include <thread>
|
|
# include <algorithm>
|
|
|
|
static const int manta_num_threads = std::thread::hardware_concurrency();
|
|
|
|
# define parallel_block \
|
|
{ \
|
|
std::vector<std::thread> threads; \
|
|
{
|
|
|
|
# define do_parallel threads.push_back( std::thread([&]() {
|
|
# define do_end \
|
|
} ) );
|
|
|
|
# define block_end \
|
|
} \
|
|
for (auto &thread : threads) { \
|
|
thread.join(); \
|
|
} \
|
|
}
|
|
|
|
#endif
|
|
|
|
#define FOR_INT_IJK(num) \
|
|
for (int k_off = 0; k_off < num; ++k_off) \
|
|
for (int j_off = 0; j_off < num; ++j_off) \
|
|
for (int i_off = 0; i_off < num; ++i_off)
|
|
|
|
using namespace std;
|
|
|
|
namespace Manta {
|
|
|
|
//! Assumes phi0<0 and phi1>=0, phi2>=0, and phi3>=0 or vice versa.
|
|
//! In particular, phi0 must not equal any of phi1, phi2 or phi3.
|
|
static Real sortedTetFraction(Real phi0, Real phi1, Real phi2, Real phi3)
|
|
{
|
|
return phi0 * phi0 * phi0 / ((phi0 - phi1) * (phi0 - phi2) * (phi0 - phi3));
|
|
}
|
|
|
|
//! Assumes phi0<0, phi1<0, and phi2>=0, and phi3>=0 or vice versa.
|
|
//! In particular, phi0 and phi1 must not equal any of phi2 and phi3.
|
|
static Real sortedPrismFraction(Real phi0, Real phi1, Real phi2, Real phi3)
|
|
{
|
|
Real a = phi0 / (phi0 - phi2);
|
|
Real b = phi0 / (phi0 - phi3);
|
|
Real c = phi1 / (phi1 - phi3);
|
|
Real d = phi1 / (phi1 - phi2);
|
|
return a * b * (1 - d) + b * (1 - c) * d + c * d;
|
|
}
|
|
|
|
Real volumeFraction(Real phi0, Real phi1, Real phi2, Real phi3)
|
|
{
|
|
sort(phi0, phi1, phi2, phi3);
|
|
if (phi3 <= 0)
|
|
return 1;
|
|
else if (phi2 <= 0)
|
|
return 1 - sortedTetFraction(phi3, phi2, phi1, phi0);
|
|
else if (phi1 <= 0)
|
|
return sortedPrismFraction(phi0, phi1, phi2, phi3);
|
|
else if (phi0 <= 0)
|
|
return sortedTetFraction(phi0, phi1, phi2, phi3);
|
|
else
|
|
return 0;
|
|
}
|
|
|
|
//! The average of the two possible decompositions of the cube into five tetrahedra.
|
|
Real volumeFraction(Real phi000,
|
|
Real phi100,
|
|
Real phi010,
|
|
Real phi110,
|
|
Real phi001,
|
|
Real phi101,
|
|
Real phi011,
|
|
Real phi111)
|
|
{
|
|
return (volumeFraction(phi000, phi001, phi101, phi011) +
|
|
volumeFraction(phi000, phi101, phi100, phi110) +
|
|
volumeFraction(phi000, phi010, phi011, phi110) +
|
|
volumeFraction(phi101, phi011, phi111, phi110) +
|
|
2 * volumeFraction(phi000, phi011, phi101, phi110) +
|
|
volumeFraction(phi100, phi101, phi001, phi111) +
|
|
volumeFraction(phi100, phi001, phi000, phi010) +
|
|
volumeFraction(phi100, phi110, phi111, phi010) +
|
|
volumeFraction(phi001, phi111, phi011, phi010) +
|
|
2 * volumeFraction(phi100, phi111, phi001, phi010)) /
|
|
12;
|
|
}
|
|
|
|
//! Kernel loop over grid with 2x base resolution!
|
|
|
|
struct KnEstimateVolumeFraction : public KernelBase {
|
|
KnEstimateVolumeFraction(Grid<Real> &volumes,
|
|
const Grid<Real> &phi,
|
|
const Vec3 &startCentre,
|
|
const Real dx)
|
|
: KernelBase(&volumes, 0), volumes(volumes), phi(phi), startCentre(startCentre), dx(dx)
|
|
{
|
|
runMessage();
|
|
run();
|
|
}
|
|
inline void op(int i,
|
|
int j,
|
|
int k,
|
|
Grid<Real> &volumes,
|
|
const Grid<Real> &phi,
|
|
const Vec3 &startCentre,
|
|
const Real dx) const
|
|
{
|
|
const Vec3 centre = startCentre + Vec3(i, j, k) * 0.5;
|
|
const Real offset = 0.5 * dx;
|
|
const int order = 1; // is sufficient
|
|
|
|
Real phi000 = phi.getInterpolatedHi(centre + Vec3(-offset, -offset, -offset), order);
|
|
Real phi001 = phi.getInterpolatedHi(centre + Vec3(-offset, -offset, +offset), order);
|
|
Real phi010 = phi.getInterpolatedHi(centre + Vec3(-offset, +offset, -offset), order);
|
|
Real phi011 = phi.getInterpolatedHi(centre + Vec3(-offset, +offset, +offset), order);
|
|
Real phi100 = phi.getInterpolatedHi(centre + Vec3(+offset, -offset, -offset), order);
|
|
Real phi101 = phi.getInterpolatedHi(centre + Vec3(+offset, -offset, +offset), order);
|
|
Real phi110 = phi.getInterpolatedHi(centre + Vec3(+offset, +offset, -offset), order);
|
|
Real phi111 = phi.getInterpolatedHi(centre + Vec3(+offset, +offset, +offset), order);
|
|
|
|
volumes(i, j, k) = volumeFraction(
|
|
phi000, phi100, phi010, phi110, phi001, phi101, phi011, phi111);
|
|
}
|
|
inline Grid<Real> &getArg0()
|
|
{
|
|
return volumes;
|
|
}
|
|
typedef Grid<Real> type0;
|
|
inline const Grid<Real> &getArg1()
|
|
{
|
|
return phi;
|
|
}
|
|
typedef Grid<Real> type1;
|
|
inline const Vec3 &getArg2()
|
|
{
|
|
return startCentre;
|
|
}
|
|
typedef Vec3 type2;
|
|
inline const Real &getArg3()
|
|
{
|
|
return dx;
|
|
}
|
|
typedef Real type3;
|
|
void runMessage()
|
|
{
|
|
debMsg("Executing kernel KnEstimateVolumeFraction ", 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 = 0; j < _maxY; j++)
|
|
for (int i = 0; i < _maxX; i++)
|
|
op(i, j, k, volumes, phi, startCentre, dx);
|
|
}
|
|
else {
|
|
const int k = 0;
|
|
for (int j = __r.begin(); j != (int)__r.end(); j++)
|
|
for (int i = 0; i < _maxX; i++)
|
|
op(i, j, k, volumes, phi, startCentre, dx);
|
|
}
|
|
}
|
|
void run()
|
|
{
|
|
if (maxZ > 1)
|
|
tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
|
|
else
|
|
tbb::parallel_for(tbb::blocked_range<IndexInt>(0, maxY), *this);
|
|
}
|
|
Grid<Real> &volumes;
|
|
const Grid<Real> φ
|
|
const Vec3 &startCentre;
|
|
const Real dx;
|
|
};
|
|
|
|
struct KnUpdateVolumeGrid : public KernelBase {
|
|
KnUpdateVolumeGrid(Grid<Real> &cVolLiquid,
|
|
Grid<Real> &uVolLiquid,
|
|
Grid<Real> &vVolLiquid,
|
|
Grid<Real> &wVolLiquid,
|
|
Grid<Real> &exVolLiquid,
|
|
Grid<Real> &eyVolLiquid,
|
|
Grid<Real> &ezVolLiquid,
|
|
const Grid<Real> &src)
|
|
: KernelBase(&cVolLiquid, 0),
|
|
cVolLiquid(cVolLiquid),
|
|
uVolLiquid(uVolLiquid),
|
|
vVolLiquid(vVolLiquid),
|
|
wVolLiquid(wVolLiquid),
|
|
exVolLiquid(exVolLiquid),
|
|
eyVolLiquid(eyVolLiquid),
|
|
ezVolLiquid(ezVolLiquid),
|
|
src(src)
|
|
{
|
|
runMessage();
|
|
run();
|
|
}
|
|
inline void op(int i,
|
|
int j,
|
|
int k,
|
|
Grid<Real> &cVolLiquid,
|
|
Grid<Real> &uVolLiquid,
|
|
Grid<Real> &vVolLiquid,
|
|
Grid<Real> &wVolLiquid,
|
|
Grid<Real> &exVolLiquid,
|
|
Grid<Real> &eyVolLiquid,
|
|
Grid<Real> &ezVolLiquid,
|
|
const Grid<Real> &src) const
|
|
{
|
|
// Work out c
|
|
cVolLiquid(i, j, k) = 0;
|
|
FOR_INT_IJK(2)
|
|
{
|
|
cVolLiquid(i, j, k) += src(2 * i + i_off, 2 * j + j_off, 2 * k + k_off);
|
|
}
|
|
cVolLiquid(i, j, k) /= 8;
|
|
|
|
// Work out u
|
|
if (i >= 1) {
|
|
uVolLiquid(i, j, k) = 0;
|
|
int base_i = 2 * i - 1;
|
|
int base_j = 2 * j;
|
|
int base_k = 2 * k;
|
|
FOR_INT_IJK(2)
|
|
{
|
|
uVolLiquid(i, j, k) += src(base_i + i_off, base_j + j_off, base_k + k_off);
|
|
}
|
|
uVolLiquid(i, j, k) /= 8;
|
|
}
|
|
|
|
// v
|
|
if (j >= 1) {
|
|
vVolLiquid(i, j, k) = 0;
|
|
int base_i = 2 * i;
|
|
int base_j = 2 * j - 1;
|
|
int base_k = 2 * k;
|
|
FOR_INT_IJK(2)
|
|
{
|
|
vVolLiquid(i, j, k) += src(base_i + i_off, base_j + j_off, base_k + k_off);
|
|
}
|
|
vVolLiquid(i, j, k) /= 8;
|
|
}
|
|
|
|
// w
|
|
if (k >= 1) {
|
|
wVolLiquid(i, j, k) = 0;
|
|
int base_i = 2 * i;
|
|
int base_j = 2 * j;
|
|
int base_k = 2 * k - 1;
|
|
FOR_INT_IJK(2)
|
|
{
|
|
wVolLiquid(i, j, k) += src(base_i + i_off, base_j + j_off, base_k + k_off);
|
|
}
|
|
wVolLiquid(i, j, k) /= 8;
|
|
}
|
|
|
|
// e-x
|
|
if (j >= 1 && k >= 1) {
|
|
exVolLiquid(i, j, k) = 0;
|
|
int base_i = 2 * i;
|
|
int base_j = 2 * j - 1;
|
|
int base_k = 2 * k - 1;
|
|
FOR_INT_IJK(2)
|
|
{
|
|
exVolLiquid(i, j, k) += src(base_i + i_off, base_j + j_off, base_k + k_off);
|
|
}
|
|
exVolLiquid(i, j, k) /= 8;
|
|
}
|
|
|
|
// e-y
|
|
if (i >= 1 && k >= 1) {
|
|
eyVolLiquid(i, j, k) = 0;
|
|
int base_i = 2 * i - 1;
|
|
int base_j = 2 * j;
|
|
int base_k = 2 * k - 1;
|
|
FOR_INT_IJK(2)
|
|
{
|
|
eyVolLiquid(i, j, k) += src(base_i + i_off, base_j + j_off, base_k + k_off);
|
|
}
|
|
eyVolLiquid(i, j, k) /= 8;
|
|
}
|
|
|
|
// e-z
|
|
if (i >= 1 && j >= 1) {
|
|
ezVolLiquid(i, j, k) = 0;
|
|
int base_i = 2 * i - 1;
|
|
int base_j = 2 * j - 1;
|
|
int base_k = 2 * k;
|
|
FOR_INT_IJK(2)
|
|
{
|
|
ezVolLiquid(i, j, k) += src(base_i + i_off, base_j + j_off, base_k + k_off);
|
|
}
|
|
ezVolLiquid(i, j, k) /= 8;
|
|
}
|
|
}
|
|
inline Grid<Real> &getArg0()
|
|
{
|
|
return cVolLiquid;
|
|
}
|
|
typedef Grid<Real> type0;
|
|
inline Grid<Real> &getArg1()
|
|
{
|
|
return uVolLiquid;
|
|
}
|
|
typedef Grid<Real> type1;
|
|
inline Grid<Real> &getArg2()
|
|
{
|
|
return vVolLiquid;
|
|
}
|
|
typedef Grid<Real> type2;
|
|
inline Grid<Real> &getArg3()
|
|
{
|
|
return wVolLiquid;
|
|
}
|
|
typedef Grid<Real> type3;
|
|
inline Grid<Real> &getArg4()
|
|
{
|
|
return exVolLiquid;
|
|
}
|
|
typedef Grid<Real> type4;
|
|
inline Grid<Real> &getArg5()
|
|
{
|
|
return eyVolLiquid;
|
|
}
|
|
typedef Grid<Real> type5;
|
|
inline Grid<Real> &getArg6()
|
|
{
|
|
return ezVolLiquid;
|
|
}
|
|
typedef Grid<Real> type6;
|
|
inline const Grid<Real> &getArg7()
|
|
{
|
|
return src;
|
|
}
|
|
typedef Grid<Real> type7;
|
|
void runMessage()
|
|
{
|
|
debMsg("Executing kernel KnUpdateVolumeGrid ", 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 = 0; j < _maxY; j++)
|
|
for (int i = 0; i < _maxX; i++)
|
|
op(i,
|
|
j,
|
|
k,
|
|
cVolLiquid,
|
|
uVolLiquid,
|
|
vVolLiquid,
|
|
wVolLiquid,
|
|
exVolLiquid,
|
|
eyVolLiquid,
|
|
ezVolLiquid,
|
|
src);
|
|
}
|
|
else {
|
|
const int k = 0;
|
|
for (int j = __r.begin(); j != (int)__r.end(); j++)
|
|
for (int i = 0; i < _maxX; i++)
|
|
op(i,
|
|
j,
|
|
k,
|
|
cVolLiquid,
|
|
uVolLiquid,
|
|
vVolLiquid,
|
|
wVolLiquid,
|
|
exVolLiquid,
|
|
eyVolLiquid,
|
|
ezVolLiquid,
|
|
src);
|
|
}
|
|
}
|
|
void run()
|
|
{
|
|
if (maxZ > 1)
|
|
tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
|
|
else
|
|
tbb::parallel_for(tbb::blocked_range<IndexInt>(0, maxY), *this);
|
|
}
|
|
Grid<Real> &cVolLiquid;
|
|
Grid<Real> &uVolLiquid;
|
|
Grid<Real> &vVolLiquid;
|
|
Grid<Real> &wVolLiquid;
|
|
Grid<Real> &exVolLiquid;
|
|
Grid<Real> &eyVolLiquid;
|
|
Grid<Real> &ezVolLiquid;
|
|
const Grid<Real> &src;
|
|
};
|
|
|
|
void computeWeights(const Grid<Real> &phi,
|
|
Grid<Real> &doubleSized,
|
|
Grid<Real> &cVolLiquid,
|
|
Grid<Real> &uVolLiquid,
|
|
Grid<Real> &vVolLiquid,
|
|
Grid<Real> &wVolLiquid,
|
|
Grid<Real> &exVolLiquid,
|
|
Grid<Real> &eyVolLiquid,
|
|
Grid<Real> &ezVolLiquid,
|
|
Real dx)
|
|
{
|
|
KnEstimateVolumeFraction(doubleSized, phi, Vec3(0.25 * dx, 0.25 * dx, 0.25 * dx), 0.5 * dx);
|
|
KnUpdateVolumeGrid(cVolLiquid,
|
|
uVolLiquid,
|
|
vVolLiquid,
|
|
wVolLiquid,
|
|
exVolLiquid,
|
|
eyVolLiquid,
|
|
ezVolLiquid,
|
|
doubleSized);
|
|
}
|
|
|
|
struct KnUpdateFaceStates : public KernelBase {
|
|
KnUpdateFaceStates(const FlagGrid &flags,
|
|
Grid<int> &uState,
|
|
Grid<int> &vState,
|
|
Grid<int> &wState)
|
|
: KernelBase(&flags, 0), flags(flags), uState(uState), vState(vState), wState(wState)
|
|
{
|
|
runMessage();
|
|
run();
|
|
}
|
|
inline void op(int i,
|
|
int j,
|
|
int k,
|
|
const FlagGrid &flags,
|
|
Grid<int> &uState,
|
|
Grid<int> &vState,
|
|
Grid<int> &wState) const
|
|
{
|
|
bool curObs = flags.isObstacle(i, j, k);
|
|
uState(i, j, k) = (i > 0 && !flags.isObstacle(i - 1, j, k) && !curObs) ?
|
|
FlagGrid::TypeFluid :
|
|
FlagGrid::TypeObstacle;
|
|
vState(i, j, k) = (j > 0 && !flags.isObstacle(i, j - 1, k) && !curObs) ?
|
|
FlagGrid::TypeFluid :
|
|
FlagGrid::TypeObstacle;
|
|
wState(i, j, k) = (k > 0 && !flags.isObstacle(i, j, k - 1) && !curObs) ?
|
|
FlagGrid::TypeFluid :
|
|
FlagGrid::TypeObstacle;
|
|
}
|
|
inline const FlagGrid &getArg0()
|
|
{
|
|
return flags;
|
|
}
|
|
typedef FlagGrid type0;
|
|
inline Grid<int> &getArg1()
|
|
{
|
|
return uState;
|
|
}
|
|
typedef Grid<int> type1;
|
|
inline Grid<int> &getArg2()
|
|
{
|
|
return vState;
|
|
}
|
|
typedef Grid<int> type2;
|
|
inline Grid<int> &getArg3()
|
|
{
|
|
return wState;
|
|
}
|
|
typedef Grid<int> type3;
|
|
void runMessage()
|
|
{
|
|
debMsg("Executing kernel KnUpdateFaceStates ", 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 = 0; j < _maxY; j++)
|
|
for (int i = 0; i < _maxX; i++)
|
|
op(i, j, k, flags, uState, vState, wState);
|
|
}
|
|
else {
|
|
const int k = 0;
|
|
for (int j = __r.begin(); j != (int)__r.end(); j++)
|
|
for (int i = 0; i < _maxX; i++)
|
|
op(i, j, k, flags, uState, vState, wState);
|
|
}
|
|
}
|
|
void run()
|
|
{
|
|
if (maxZ > 1)
|
|
tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
|
|
else
|
|
tbb::parallel_for(tbb::blocked_range<IndexInt>(0, maxY), *this);
|
|
}
|
|
const FlagGrid &flags;
|
|
Grid<int> &uState;
|
|
Grid<int> &vState;
|
|
Grid<int> &wState;
|
|
};
|
|
|
|
struct KnApplyVelocities : public KernelBase {
|
|
KnApplyVelocities(MACGrid &dst,
|
|
const Grid<int> &uState,
|
|
const Grid<int> &vState,
|
|
const Grid<int> &wState,
|
|
Grid<Real> &srcU,
|
|
Grid<Real> &srcV,
|
|
Grid<Real> &srcW)
|
|
: KernelBase(&dst, 0),
|
|
dst(dst),
|
|
uState(uState),
|
|
vState(vState),
|
|
wState(wState),
|
|
srcU(srcU),
|
|
srcV(srcV),
|
|
srcW(srcW)
|
|
{
|
|
runMessage();
|
|
run();
|
|
}
|
|
inline void op(int i,
|
|
int j,
|
|
int k,
|
|
MACGrid &dst,
|
|
const Grid<int> &uState,
|
|
const Grid<int> &vState,
|
|
const Grid<int> &wState,
|
|
Grid<Real> &srcU,
|
|
Grid<Real> &srcV,
|
|
Grid<Real> &srcW) const
|
|
{
|
|
dst(i, j, k).x = (uState(i, j, k) == FlagGrid::TypeFluid) ? srcU(i, j, k) : 0;
|
|
dst(i, j, k).y = (vState(i, j, k) == FlagGrid::TypeFluid) ? srcV(i, j, k) : 0;
|
|
if (dst.is3D())
|
|
dst(i, j, k).z = (wState(i, j, k) == FlagGrid::TypeFluid) ? srcW(i, j, k) : 0;
|
|
}
|
|
inline MACGrid &getArg0()
|
|
{
|
|
return dst;
|
|
}
|
|
typedef MACGrid type0;
|
|
inline const Grid<int> &getArg1()
|
|
{
|
|
return uState;
|
|
}
|
|
typedef Grid<int> type1;
|
|
inline const Grid<int> &getArg2()
|
|
{
|
|
return vState;
|
|
}
|
|
typedef Grid<int> type2;
|
|
inline const Grid<int> &getArg3()
|
|
{
|
|
return wState;
|
|
}
|
|
typedef Grid<int> type3;
|
|
inline Grid<Real> &getArg4()
|
|
{
|
|
return srcU;
|
|
}
|
|
typedef Grid<Real> type4;
|
|
inline Grid<Real> &getArg5()
|
|
{
|
|
return srcV;
|
|
}
|
|
typedef Grid<Real> type5;
|
|
inline Grid<Real> &getArg6()
|
|
{
|
|
return srcW;
|
|
}
|
|
typedef Grid<Real> type6;
|
|
void runMessage()
|
|
{
|
|
debMsg("Executing kernel KnApplyVelocities ", 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 = 0; j < _maxY; j++)
|
|
for (int i = 0; i < _maxX; i++)
|
|
op(i, j, k, dst, uState, vState, wState, srcU, srcV, srcW);
|
|
}
|
|
else {
|
|
const int k = 0;
|
|
for (int j = __r.begin(); j != (int)__r.end(); j++)
|
|
for (int i = 0; i < _maxX; i++)
|
|
op(i, j, k, dst, uState, vState, wState, srcU, srcV, srcW);
|
|
}
|
|
}
|
|
void run()
|
|
{
|
|
if (maxZ > 1)
|
|
tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
|
|
else
|
|
tbb::parallel_for(tbb::blocked_range<IndexInt>(0, maxY), *this);
|
|
}
|
|
MACGrid &dst;
|
|
const Grid<int> &uState;
|
|
const Grid<int> &vState;
|
|
const Grid<int> &wState;
|
|
Grid<Real> &srcU;
|
|
Grid<Real> &srcV;
|
|
Grid<Real> &srcW;
|
|
};
|
|
|
|
void solveViscosity(const FlagGrid &flags,
|
|
MACGrid &vel,
|
|
Grid<Real> &cVolLiquid,
|
|
Grid<Real> &uVolLiquid,
|
|
Grid<Real> &vVolLiquid,
|
|
Grid<Real> &wVolLiquid,
|
|
Grid<Real> &exVolLiquid,
|
|
Grid<Real> &eyVolLiquid,
|
|
Grid<Real> &ezVolLiquid,
|
|
Grid<Real> &viscosity,
|
|
const Real dt,
|
|
const Real dx,
|
|
const Real cgAccuracy,
|
|
const Real cgMaxIterFac)
|
|
{
|
|
const Real factor = dt * square(1.0 / dx);
|
|
const int maxIter = (int)(cgMaxIterFac * flags.getSize().max()) * (flags.is3D() ? 1 : 4);
|
|
GridCg<ApplyMatrixViscosityU> *uGcg;
|
|
GridCg<ApplyMatrixViscosityV> *vGcg;
|
|
GridCg<ApplyMatrixViscosityW> *wGcg;
|
|
|
|
// Tmp grids for CG solve in U, V, W dimensions
|
|
FluidSolver *parent = flags.getParent();
|
|
Grid<Real> uResidual(parent);
|
|
Grid<Real> vResidual(parent);
|
|
Grid<Real> wResidual(parent);
|
|
Grid<Real> uSearch(parent);
|
|
Grid<Real> vSearch(parent);
|
|
Grid<Real> wSearch(parent);
|
|
Grid<Real> uTmp(parent);
|
|
Grid<Real> vTmp(parent);
|
|
Grid<Real> wTmp(parent);
|
|
Grid<Real> uRhs(parent);
|
|
Grid<Real> vRhs(parent);
|
|
Grid<Real> wRhs(parent);
|
|
|
|
// A matrix U grids
|
|
Grid<Real> uA0(parent); // diagonal elements in A matrix
|
|
Grid<Real> uAplusi(parent); // neighbor at i+1
|
|
Grid<Real> uAplusj(parent); // neighbor at j+1
|
|
Grid<Real> uAplusk(parent); // neighbor at k+1
|
|
Grid<Real> uAminusi(parent); // neighbor at i-1
|
|
Grid<Real> uAminusj(parent); // neighbor at j-1
|
|
Grid<Real> uAminusk(parent); // neighbor at k-1
|
|
Grid<Real> uAhelper1(parent); // additional helper grids for off diagonal elements
|
|
Grid<Real> uAhelper2(parent);
|
|
Grid<Real> uAhelper3(parent);
|
|
Grid<Real> uAhelper4(parent);
|
|
Grid<Real> uAhelper5(parent);
|
|
Grid<Real> uAhelper6(parent);
|
|
Grid<Real> uAhelper7(parent);
|
|
Grid<Real> uAhelper8(parent);
|
|
|
|
// A matrix V grids
|
|
Grid<Real> vA0(parent);
|
|
Grid<Real> vAplusi(parent);
|
|
Grid<Real> vAplusj(parent);
|
|
Grid<Real> vAplusk(parent);
|
|
Grid<Real> vAminusi(parent);
|
|
Grid<Real> vAminusj(parent);
|
|
Grid<Real> vAminusk(parent);
|
|
Grid<Real> vAhelper1(parent);
|
|
Grid<Real> vAhelper2(parent);
|
|
Grid<Real> vAhelper3(parent);
|
|
Grid<Real> vAhelper4(parent);
|
|
Grid<Real> vAhelper5(parent);
|
|
Grid<Real> vAhelper6(parent);
|
|
Grid<Real> vAhelper7(parent);
|
|
Grid<Real> vAhelper8(parent);
|
|
|
|
// A matrix W grids
|
|
Grid<Real> wA0(parent);
|
|
Grid<Real> wAplusi(parent);
|
|
Grid<Real> wAplusj(parent);
|
|
Grid<Real> wAplusk(parent);
|
|
Grid<Real> wAminusi(parent);
|
|
Grid<Real> wAminusj(parent);
|
|
Grid<Real> wAminusk(parent);
|
|
Grid<Real> wAhelper1(parent);
|
|
Grid<Real> wAhelper2(parent);
|
|
Grid<Real> wAhelper3(parent);
|
|
Grid<Real> wAhelper4(parent);
|
|
Grid<Real> wAhelper5(parent);
|
|
Grid<Real> wAhelper6(parent);
|
|
Grid<Real> wAhelper7(parent);
|
|
Grid<Real> wAhelper8(parent);
|
|
|
|
// Solution grids for CG solvers
|
|
Grid<Real> uSolution(parent);
|
|
Grid<Real> vSolution(parent);
|
|
Grid<Real> wSolution(parent);
|
|
|
|
// Save state of voxel face (fluid or obstacle)
|
|
Grid<int> uState(parent);
|
|
Grid<int> vState(parent);
|
|
Grid<int> wState(parent);
|
|
|
|
// Save state of voxel face (fluid or obstacle)
|
|
KnUpdateFaceStates(flags, uState, vState, wState);
|
|
|
|
// Shorter names for flags, we will use them often
|
|
int isFluid = FlagGrid::TypeFluid;
|
|
int isObstacle = FlagGrid::TypeObstacle;
|
|
|
|
// Main viscosity loop: construct A matrices and rhs's in all dimensions
|
|
FOR_IJK_BND(flags, 1)
|
|
{
|
|
|
|
// U-terms: 2u_xx+ v_xy +uyy + u_zz + w_xz
|
|
if (uState(i, j, k) == isFluid) {
|
|
|
|
uRhs(i, j, k) = uVolLiquid(i, j, k) * vel(i, j, k).x;
|
|
uA0(i, j, k) = uVolLiquid(i, j, k);
|
|
|
|
Real viscRight = viscosity(i, j, k);
|
|
Real viscLeft = viscosity(i - 1, j, k);
|
|
Real volRight = cVolLiquid(i, j, k);
|
|
Real volLeft = cVolLiquid(i - 1, j, k);
|
|
|
|
Real viscTop = 0.25 * (viscosity(i - 1, j + 1, k) + viscosity(i - 1, j, k) +
|
|
viscosity(i, j + 1, k) + viscosity(i, j, k));
|
|
Real viscBottom = 0.25 * (viscosity(i - 1, j, k) + viscosity(i - 1, j - 1, k) +
|
|
viscosity(i, j, k) + viscosity(i, j - 1, k));
|
|
Real volTop = ezVolLiquid(i, j + 1, k);
|
|
Real volBottom = ezVolLiquid(i, j, k);
|
|
|
|
Real viscFront = 0.25 * (viscosity(i - 1, j, k + 1) + viscosity(i - 1, j, k) +
|
|
viscosity(i, j, k + 1) + viscosity(i, j, k));
|
|
Real viscBack = 0.25 * (viscosity(i - 1, j, k) + viscosity(i - 1, j, k - 1) +
|
|
viscosity(i, j, k) + viscosity(i, j, k - 1));
|
|
Real volFront = eyVolLiquid(i, j, k + 1);
|
|
Real volBack = eyVolLiquid(i, j, k);
|
|
|
|
Real factorRight = 2 * factor * viscRight * volRight;
|
|
Real factorLeft = 2 * factor * viscLeft * volLeft;
|
|
Real factorTop = factor * viscTop * volTop;
|
|
Real factorBottom = factor * viscBottom * volBottom;
|
|
Real factorFront = factor * viscFront * volFront;
|
|
Real factorBack = factor * viscBack * volBack;
|
|
|
|
// u_x_right
|
|
uA0(i, j, k) += factorRight;
|
|
if (uState(i + 1, j, k) == isFluid) {
|
|
uAplusi(i, j, k) += -factorRight;
|
|
}
|
|
else if (uState(i + 1, j, k) == isObstacle) {
|
|
uRhs(i, j, k) -= -vel(i + 1, j, k).x * factorRight;
|
|
}
|
|
|
|
// u_x_left
|
|
uA0(i, j, k) += factorLeft;
|
|
if (uState(i - 1, j, k) == isFluid) {
|
|
uAminusi(i, j, k) += -factorLeft;
|
|
}
|
|
else if (uState(i - 1, j, k) == isObstacle) {
|
|
uRhs(i, j, k) -= -vel(i - 1, j, k).x * factorLeft;
|
|
}
|
|
|
|
// u_y_top
|
|
uA0(i, j, k) += factorTop;
|
|
if (uState(i, j + 1, k) == isFluid) {
|
|
uAplusj(i, j, k) += -factorTop;
|
|
}
|
|
else if (uState(i, j + 1, k) == isObstacle) {
|
|
uRhs(i, j, k) -= -vel(i, j + 1, k).x * factorTop;
|
|
}
|
|
|
|
// u_y_bottom
|
|
uA0(i, j, k) += factorBottom;
|
|
if (uState(i, j - 1, k) == isFluid) {
|
|
uAminusj(i, j, k) += -factorBottom;
|
|
}
|
|
else if (uState(i, j - 1, k) == isObstacle) {
|
|
uRhs(i, j, k) -= -vel(i, j - 1, k).x * factorBottom;
|
|
}
|
|
|
|
// u_z_front
|
|
uA0(i, j, k) += factorFront;
|
|
if (uState(i, j, k + 1) == isFluid) {
|
|
uAplusk(i, j, k) += -factorFront;
|
|
}
|
|
else if (uState(i, j, k + 1) == isObstacle) {
|
|
uRhs(i, j, k) -= -vel(i, j, k + 1).x * factorFront;
|
|
}
|
|
|
|
// u_z_back
|
|
uA0(i, j, k) += factorBack;
|
|
if (uState(i, j, k - 1) == isFluid) {
|
|
uAminusk(i, j, k) += -factorBack;
|
|
}
|
|
else if (uState(i, j, k - 1) == isObstacle) {
|
|
uRhs(i, j, k) -= -vel(i, j, k - 1).x * factorBack;
|
|
}
|
|
|
|
// v_x_top
|
|
if (vState(i, j + 1, k) == isFluid) {
|
|
uAhelper1(i, j, k) += -factorTop;
|
|
}
|
|
else if (vState(i, j + 1, k) == isObstacle) {
|
|
uRhs(i, j, k) -= -vel(i, j + 1, k).y * factorTop;
|
|
}
|
|
|
|
if (vState(i - 1, j + 1, k) == isFluid) {
|
|
uAhelper2(i, j, k) += factorTop;
|
|
}
|
|
else if (vState(i - 1, j + 1, k) == isObstacle) {
|
|
uRhs(i, j, k) -= vel(i - 1, j + 1, k).y * factorTop;
|
|
}
|
|
|
|
// v_x_bottom
|
|
if (vState(i, j, k) == isFluid) {
|
|
uAhelper3(i, j, k) += factorBottom;
|
|
}
|
|
else if (vState(i, j, k) == isObstacle) {
|
|
uRhs(i, j, k) -= vel(i, j, k).y * factorBottom;
|
|
}
|
|
|
|
if (vState(i - 1, j, k) == isFluid) {
|
|
uAhelper4(i, j, k) += -factorBottom;
|
|
}
|
|
else if (vState(i - 1, j, k) == isObstacle) {
|
|
uRhs(i, j, k) -= -vel(i - 1, j, k).y * factorBottom;
|
|
}
|
|
|
|
// w_x_front
|
|
if (wState(i, j, k + 1) == isFluid) {
|
|
uAhelper5(i, j, k) += -factorFront;
|
|
}
|
|
else if (wState(i, j, k + 1) == isObstacle) {
|
|
uRhs(i, j, k) -= -vel(i, j, k + 1).z * factorFront;
|
|
}
|
|
|
|
if (wState(i - 1, j, k + 1) == isFluid) {
|
|
uAhelper6(i, j, k) += factorFront;
|
|
}
|
|
else if (wState(i - 1, j, k + 1) == isObstacle) {
|
|
uRhs(i, j, k) -= vel(i - 1, j, k + 1).z * factorFront;
|
|
}
|
|
|
|
// w_x_back
|
|
if (wState(i, j, k) == isFluid) {
|
|
uAhelper7(i, j, k) += factorBack;
|
|
}
|
|
else if (wState(i, j, k) == isObstacle) {
|
|
uRhs(i, j, k) -= vel(i, j, k).z * factorBack;
|
|
}
|
|
|
|
if (wState(i - 1, j, k) == isFluid) {
|
|
uAhelper8(i, j, k) += -factorBack;
|
|
}
|
|
else if (wState(i - 1, j, k) == isObstacle) {
|
|
uRhs(i, j, k) -= -vel(i - 1, j, k).z * factorBack;
|
|
}
|
|
}
|
|
|
|
// V-terms: vxx + 2vyy + vzz + u_yx + w_yz
|
|
if (vState(i, j, k) == isFluid) {
|
|
|
|
vRhs(i, j, k) = vVolLiquid(i, j, k) * vel(i, j, k).y;
|
|
vA0(i, j, k) = vVolLiquid(i, j, k);
|
|
|
|
Real viscRight = 0.25 * (viscosity(i, j - 1, k) + viscosity(i + 1, j - 1, k) +
|
|
viscosity(i, j, k) + viscosity(i + 1, j, k));
|
|
Real viscLeft = 0.25 * (viscosity(i, j - 1, k) + viscosity(i - 1, j - 1, k) +
|
|
viscosity(i, j, k) + viscosity(i - 1, j, k));
|
|
Real volRight = ezVolLiquid(i + 1, j, k);
|
|
Real volLeft = ezVolLiquid(i, j, k);
|
|
|
|
Real viscTop = viscosity(i, j, k);
|
|
Real viscBottom = viscosity(i, j - 1, k);
|
|
Real volTop = cVolLiquid(i, j, k);
|
|
Real volBottom = cVolLiquid(i, j - 1, k);
|
|
|
|
Real viscFront = 0.25 * (viscosity(i, j - 1, k) + viscosity(i, j - 1, k + 1) +
|
|
viscosity(i, j, k) + viscosity(i, j, k + 1));
|
|
Real viscBack = 0.25 * (viscosity(i, j - 1, k) + viscosity(i, j - 1, k - 1) +
|
|
viscosity(i, j, k) + viscosity(i, j, k - 1));
|
|
Real volFront = exVolLiquid(i, j, k + 1);
|
|
Real volBack = exVolLiquid(i, j, k);
|
|
|
|
Real factorRight = factor * viscRight * volRight;
|
|
Real factorLeft = factor * viscLeft * volLeft;
|
|
Real factorTop = 2 * factor * viscTop * volTop;
|
|
Real factorBottom = 2 * factor * viscBottom * volBottom;
|
|
Real factorFront = factor * viscFront * volFront;
|
|
Real factorBack = factor * viscBack * volBack;
|
|
|
|
// v_x_right
|
|
vA0(i, j, k) += factorRight;
|
|
if (vState(i + 1, j, k) == isFluid) {
|
|
vAplusi(i, j, k) += -factorRight;
|
|
}
|
|
else if (vState(i + 1, j, k) == isObstacle) {
|
|
vRhs(i, j, k) -= -vel(i + 1, j, k).y * factorRight;
|
|
}
|
|
|
|
// v_x_left
|
|
vA0(i, j, k) += factorLeft;
|
|
if (vState(i - 1, j, k) == isFluid) {
|
|
vAminusi(i, j, k) += -factorLeft;
|
|
}
|
|
else if (vState(i - 1, j, k) == isObstacle) {
|
|
vRhs(i, j, k) -= -vel(i - 1, j, k).y * factorLeft;
|
|
}
|
|
|
|
// vy_top
|
|
vA0(i, j, k) += factorTop;
|
|
if (vState(i, j + 1, k) == isFluid) {
|
|
vAplusj(i, j, k) += -factorTop;
|
|
}
|
|
else if (vState(i, j + 1, k) == isObstacle) {
|
|
vRhs(i, j, k) -= -vel(i, j + 1, k).y * factorTop;
|
|
}
|
|
|
|
// vy_bottom
|
|
vA0(i, j, k) += factorBottom;
|
|
if (vState(i, j - 1, k) == isFluid) {
|
|
vAminusj(i, j, k) += -factorBottom;
|
|
}
|
|
else if (vState(i, j - 1, k) == isObstacle) {
|
|
vRhs(i, j, k) -= -vel(i, j - 1, k).y * factorBottom;
|
|
}
|
|
|
|
// v_z_front
|
|
vA0(i, j, k) += factorFront;
|
|
if (vState(i, j, k + 1) == isFluid) {
|
|
vAplusk(i, j, k) += -factorFront;
|
|
}
|
|
else if (vState(i, j, k + 1) == isObstacle) {
|
|
vRhs(i, j, k) -= -vel(i, j, k + 1).y * factorFront;
|
|
}
|
|
|
|
// v_z_back
|
|
vA0(i, j, k) += factorBack;
|
|
if (vState(i, j, k - 1) == isFluid) {
|
|
vAminusk(i, j, k) += -factorBack;
|
|
}
|
|
else if (vState(i, j, k - 1) == isObstacle) {
|
|
vRhs(i, j, k) -= -vel(i, j, k - 1).y * factorBack;
|
|
}
|
|
|
|
// u_y_right
|
|
if (uState(i + 1, j, k) == isFluid) {
|
|
vAhelper1(i, j, k) += -factorRight;
|
|
}
|
|
else if (uState(i + 1, j, k) == isObstacle) {
|
|
vRhs(i, j, k) -= -vel(i + 1, j, k).x * factorRight;
|
|
}
|
|
|
|
if (uState(i + 1, j - 1, k) == isFluid) {
|
|
vAhelper2(i, j, k) += factorRight;
|
|
}
|
|
else if (uState(i + 1, j - 1, k) == isObstacle) {
|
|
vRhs(i, j, k) -= vel(i + 1, j - 1, k).x * factorRight;
|
|
}
|
|
|
|
// u_y_left
|
|
if (uState(i, j, k) == isFluid) {
|
|
vAhelper3(i, j, k) += factorLeft;
|
|
}
|
|
else if (uState(i, j, k) == isObstacle) {
|
|
vRhs(i, j, k) -= vel(i, j, k).x * factorLeft;
|
|
}
|
|
|
|
if (uState(i, j - 1, k) == isFluid) {
|
|
vAhelper4(i, j, k) += -factorLeft;
|
|
}
|
|
else if (uState(i, j - 1, k) == isObstacle) {
|
|
vRhs(i, j, k) -= -vel(i, j - 1, k).x * factorLeft;
|
|
}
|
|
|
|
// w_y_front
|
|
if (wState(i, j, k + 1) == isFluid) {
|
|
vAhelper5(i, j, k) += -factorFront;
|
|
}
|
|
else if (wState(i, j, k + 1) == isObstacle) {
|
|
vRhs(i, j, k) -= -vel(i, j, k + 1).z * factorFront;
|
|
}
|
|
|
|
if (wState(i, j - 1, k + 1) == isFluid) {
|
|
vAhelper6(i, j, k) += factorFront;
|
|
}
|
|
else if (wState(i, j - 1, k + 1) == isObstacle) {
|
|
vRhs(i, j, k) -= vel(i, j - 1, k + 1).z * factorFront;
|
|
}
|
|
|
|
// w_y_back
|
|
if (wState(i, j, k) == isFluid) {
|
|
vAhelper7(i, j, k) += factorBack;
|
|
}
|
|
else if (wState(i, j, k) == isObstacle) {
|
|
vRhs(i, j, k) -= vel(i, j, k).z * factorBack;
|
|
}
|
|
|
|
if (wState(i, j - 1, k) == isFluid) {
|
|
vAhelper8(i, j, k) += -factorBack;
|
|
}
|
|
else if (wState(i, j - 1, k) == isObstacle) {
|
|
vRhs(i, j, k) -= -vel(i, j - 1, k).z * factorBack;
|
|
}
|
|
}
|
|
|
|
// W-terms: wxx+ wyy+ 2wzz + u_zx + v_zy
|
|
if (wState(i, j, k) == isFluid) {
|
|
|
|
wRhs(i, j, k) = wVolLiquid(i, j, k) * vel(i, j, k).z;
|
|
wA0(i, j, k) = wVolLiquid(i, j, k);
|
|
|
|
Real viscRight = 0.25 * (viscosity(i, j, k) + viscosity(i, j, k - 1) +
|
|
viscosity(i + 1, j, k) + viscosity(i + 1, j, k - 1));
|
|
Real viscLeft = 0.25 * (viscosity(i, j, k) + viscosity(i, j, k - 1) +
|
|
viscosity(i - 1, j, k) + viscosity(i - 1, j, k - 1));
|
|
Real volRight = eyVolLiquid(i + 1, j, k);
|
|
Real volLeft = eyVolLiquid(i, j, k);
|
|
|
|
Real viscTop = 0.25 * (viscosity(i, j, k) + viscosity(i, j, k - 1) + viscosity(i, j + 1, k) +
|
|
viscosity(i, j + 1, k - 1));
|
|
Real viscBottom = 0.25 * (viscosity(i, j, k) + viscosity(i, j, k - 1) +
|
|
viscosity(i, j - 1, k) + viscosity(i, j - 1, k - 1));
|
|
Real volTop = exVolLiquid(i, j + 1, k);
|
|
Real volBottom = exVolLiquid(i, j, k);
|
|
|
|
Real viscFront = viscosity(i, j, k);
|
|
Real viscBack = viscosity(i, j, k - 1);
|
|
Real volFront = cVolLiquid(i, j, k);
|
|
Real volBack = cVolLiquid(i, j, k - 1);
|
|
|
|
Real factorRight = factor * viscRight * volRight;
|
|
Real factorLeft = factor * viscLeft * volLeft;
|
|
Real factorTop = factor * viscTop * volTop;
|
|
Real factorBottom = factor * viscBottom * volBottom;
|
|
Real factorFront = 2 * factor * viscFront * volFront;
|
|
Real factorBack = 2 * factor * viscBack * volBack;
|
|
|
|
// w_x_right
|
|
wA0(i, j, k) += factorRight;
|
|
if (wState(i + 1, j, k) == isFluid) {
|
|
wAplusi(i, j, k) += -factorRight;
|
|
}
|
|
else if (wState(i + 1, j, k) == isObstacle) {
|
|
wRhs(i, j, k) -= -vel(i + 1, j, k).z * factorRight;
|
|
}
|
|
|
|
// w_x_left
|
|
wA0(i, j, k) += factorLeft;
|
|
if (wState(i - 1, j, k) == isFluid) {
|
|
wAminusi(i, j, k) += -factorLeft;
|
|
}
|
|
else if (wState(i - 1, j, k) == isObstacle) {
|
|
wRhs(i, j, k) -= -vel(i - 1, j, k).z * factorLeft;
|
|
}
|
|
|
|
// w_y_top
|
|
wA0(i, j, k) += factorTop;
|
|
if (wState(i, j + 1, k) == isFluid) {
|
|
wAplusj(i, j, k) += -factorTop;
|
|
}
|
|
else if (wState(i, j + 1, k) == isObstacle) {
|
|
wRhs(i, j, k) -= -vel(i, j + 1, k).z * factorTop;
|
|
}
|
|
|
|
// w_y_bottom
|
|
wA0(i, j, k) += factorBottom;
|
|
if (wState(i, j - 1, k) == isFluid) {
|
|
wAminusj(i, j, k) += -factorBottom;
|
|
}
|
|
else if (wState(i, j - 1, k) == isObstacle) {
|
|
wRhs(i, j, k) -= -vel(i, j - 1, k).z * factorBottom;
|
|
}
|
|
|
|
// w_z_front
|
|
wA0(i, j, k) += factorFront;
|
|
if (wState(i, j, k + 1) == isFluid) {
|
|
wAplusk(i, j, k) += -factorFront;
|
|
}
|
|
else if (wState(i, j, k + 1) == isObstacle) {
|
|
wRhs(i, j, k) -= -vel(i, j, k + 1).z * factorFront;
|
|
}
|
|
|
|
// w_z_back
|
|
wA0(i, j, k) += factorBack;
|
|
if (wState(i, j, k - 1) == isFluid) {
|
|
wAminusk(i, j, k) += -factorBack;
|
|
}
|
|
else if (wState(i, j, k - 1) == isObstacle) {
|
|
wRhs(i, j, k) -= -vel(i, j, k - 1).z * factorBack;
|
|
}
|
|
|
|
// u_z_right
|
|
if (uState(i + 1, j, k) == isFluid) {
|
|
wAhelper1(i, j, k) += -factorRight;
|
|
}
|
|
else if (uState(i + 1, j, k) == isObstacle) {
|
|
wRhs(i, j, k) -= -vel(i + 1, j, k).x * factorRight;
|
|
}
|
|
|
|
if (uState(i + 1, j, k - 1) == isFluid) {
|
|
wAhelper2(i, j, k) += factorRight;
|
|
}
|
|
else if (uState(i + 1, j, k - 1) == isObstacle) {
|
|
wRhs(i, j, k) -= vel(i + 1, j, k - 1).x * factorRight;
|
|
}
|
|
|
|
// u_z_left
|
|
if (uState(i, j, k) == isFluid) {
|
|
wAhelper3(i, j, k) += factorLeft;
|
|
}
|
|
else if (uState(i, j, k) == isObstacle) {
|
|
wRhs(i, j, k) -= vel(i, j, k).x * factorLeft;
|
|
}
|
|
|
|
if (uState(i, j, k - 1) == isFluid) {
|
|
wAhelper4(i, j, k) += -factorLeft;
|
|
}
|
|
else if (uState(i, j, k - 1) == isObstacle) {
|
|
wRhs(i, j, k) -= -vel(i, j, k - 1).x * factorLeft;
|
|
}
|
|
|
|
// v_z_top
|
|
if (vState(i, j + 1, k) == isFluid) {
|
|
wAhelper5(i, j, k) += -factorTop;
|
|
}
|
|
else if (vState(i, j + 1, k) == isObstacle) {
|
|
wRhs(i, j, k) -= -vel(i, j + 1, k).y * factorTop;
|
|
}
|
|
|
|
if (vState(i, j + 1, k - 1) == isFluid) {
|
|
wAhelper6(i, j, k) += factorTop;
|
|
}
|
|
else if (vState(i, j + 1, k - 1) == isObstacle) {
|
|
wRhs(i, j, k) -= vel(i, j + 1, k - 1).y * factorTop;
|
|
}
|
|
|
|
// v_z_bottom
|
|
if (vState(i, j, k) == isFluid) {
|
|
wAhelper7(i, j, k) += factorBottom;
|
|
}
|
|
else if (vState(i, j, k) == isObstacle) {
|
|
wRhs(i, j, k) -= vel(i, j, k).y * factorBottom;
|
|
}
|
|
|
|
if (vState(i, j, k - 1) == isFluid) {
|
|
wAhelper8(i, j, k) += -factorBottom;
|
|
}
|
|
else if (vState(i, j, k - 1) == isObstacle) {
|
|
wRhs(i, j, k) -= -vel(i, j, k - 1).y * factorBottom;
|
|
}
|
|
}
|
|
}
|
|
|
|
// CG solver for U
|
|
if (flags.is3D()) {
|
|
vector<Grid<Real> *> uMatA{&uA0,
|
|
&uAplusi,
|
|
&uAplusj,
|
|
&uAplusk,
|
|
&uAminusi,
|
|
&uAminusj,
|
|
&uAminusk,
|
|
&uAhelper1,
|
|
&uAhelper2,
|
|
&uAhelper3,
|
|
&uAhelper4,
|
|
&uAhelper5,
|
|
&uAhelper6,
|
|
&uAhelper7,
|
|
&uAhelper8};
|
|
vector<Grid<Real> *> uVecRhs{&vRhs, &wRhs};
|
|
uGcg = new GridCg<ApplyMatrixViscosityU>(
|
|
uSolution, uRhs, uResidual, uSearch, flags, uTmp, uMatA, uVecRhs);
|
|
}
|
|
else {
|
|
errMsg("Viscosity: 2D Matrix application not yet supported in viscosity solver");
|
|
}
|
|
|
|
// CG solver for V
|
|
if (flags.is3D()) {
|
|
vector<Grid<Real> *> vMatA{&vA0,
|
|
&vAplusi,
|
|
&vAplusj,
|
|
&vAplusk,
|
|
&vAminusi,
|
|
&vAminusj,
|
|
&vAminusk,
|
|
&vAhelper1,
|
|
&vAhelper2,
|
|
&vAhelper3,
|
|
&vAhelper4,
|
|
&vAhelper5,
|
|
&vAhelper6,
|
|
&vAhelper7,
|
|
&vAhelper8};
|
|
vector<Grid<Real> *> vVecRhs{&uRhs, &wRhs};
|
|
vGcg = new GridCg<ApplyMatrixViscosityV>(
|
|
vSolution, vRhs, vResidual, vSearch, flags, vTmp, vMatA, vVecRhs);
|
|
}
|
|
else {
|
|
errMsg("Viscosity: 2D Matrix application not yet supported in viscosity solver");
|
|
}
|
|
|
|
// CG solver for W
|
|
if (flags.is3D()) {
|
|
vector<Grid<Real> *> wMatA{&wA0,
|
|
&wAplusi,
|
|
&wAplusj,
|
|
&wAplusk,
|
|
&wAminusi,
|
|
&wAminusj,
|
|
&wAminusk,
|
|
&wAhelper1,
|
|
&wAhelper2,
|
|
&wAhelper3,
|
|
&wAhelper4,
|
|
&wAhelper5,
|
|
&wAhelper6,
|
|
&wAhelper7,
|
|
&wAhelper8};
|
|
vector<Grid<Real> *> wVecRhs{&uRhs, &vRhs};
|
|
wGcg = new GridCg<ApplyMatrixViscosityW>(
|
|
wSolution, wRhs, wResidual, wSearch, flags, wTmp, wMatA, wVecRhs);
|
|
}
|
|
else {
|
|
errMsg("Viscosity: 2D Matrix application not yet supported in viscosity solver");
|
|
}
|
|
|
|
// Same accuracy for all dimensions
|
|
uGcg->setAccuracy(cgAccuracy);
|
|
vGcg->setAccuracy(cgAccuracy);
|
|
wGcg->setAccuracy(cgAccuracy);
|
|
|
|
// CG solve. Preconditioning not supported yet. Instead, U, V, W can optionally be solved in
|
|
// parallel.
|
|
for (int uIter = 0, vIter = 0, wIter = 0; uIter < maxIter || vIter < maxIter || wIter < maxIter;
|
|
uIter++, vIter++, wIter++) {
|
|
#if ENABLE_PARALLEL == 1
|
|
parallel_block do_parallel
|
|
#endif
|
|
if (uIter < maxIter && !uGcg->iterate()) uIter = maxIter;
|
|
#if ENABLE_PARALLEL == 1
|
|
do_end do_parallel
|
|
#endif
|
|
if (vIter < maxIter && !vGcg->iterate()) vIter = maxIter;
|
|
#if ENABLE_PARALLEL == 1
|
|
do_end do_parallel
|
|
#endif
|
|
if (wIter < maxIter && !wGcg->iterate()) wIter = maxIter;
|
|
#if ENABLE_PARALLEL == 1
|
|
do_end block_end
|
|
#endif
|
|
|
|
// Make sure that next CG iteration has updated rhs grids
|
|
uRhs.copyFrom(uSearch);
|
|
vRhs.copyFrom(vSearch);
|
|
wRhs.copyFrom(wSearch);
|
|
}
|
|
debMsg(
|
|
"Viscosity: solveViscosity() done. "
|
|
"Iterations (u,v,w): ("
|
|
<< uGcg->getIterations() << "," << vGcg->getIterations() << "," << wGcg->getIterations()
|
|
<< "), "
|
|
"Residual (u,v,w): ("
|
|
<< uGcg->getResNorm() << "," << vGcg->getResNorm() << "," << wGcg->getResNorm() << ")",
|
|
2);
|
|
|
|
delete uGcg;
|
|
delete vGcg;
|
|
delete wGcg;
|
|
|
|
// Apply solutions to global velocity grid
|
|
KnApplyVelocities(vel, uState, vState, wState, uSolution, vSolution, wSolution);
|
|
}
|
|
|
|
//! To use the viscosity plugin, scenes must call this function before solving pressure.
|
|
//! Note that the 'volumes' grid uses 2x the base resolution
|
|
|
|
void applyViscosity(const FlagGrid &flags,
|
|
const Grid<Real> &phi,
|
|
MACGrid &vel,
|
|
Grid<Real> &volumes,
|
|
Grid<Real> &viscosity,
|
|
const Real cgAccuracy = 1e-9,
|
|
const Real cgMaxIterFac = 1.5)
|
|
{
|
|
const Real dx = flags.getParent()->getDx();
|
|
const Real dt = flags.getParent()->getDt();
|
|
|
|
// Reserve temp grids for volume weight calculation
|
|
FluidSolver *parent = flags.getParent();
|
|
Grid<Real> cVolLiquid(parent);
|
|
Grid<Real> uVolLiquid(parent);
|
|
Grid<Real> vVolLiquid(parent);
|
|
Grid<Real> wVolLiquid(parent);
|
|
Grid<Real> exVolLiquid(parent);
|
|
Grid<Real> eyVolLiquid(parent);
|
|
Grid<Real> ezVolLiquid(parent);
|
|
|
|
// Ensure final weight grid gets cleared at every step
|
|
volumes.clear();
|
|
|
|
// Save viscous fluid volume in double-sized volumes grid
|
|
computeWeights(phi,
|
|
volumes,
|
|
cVolLiquid,
|
|
uVolLiquid,
|
|
vVolLiquid,
|
|
wVolLiquid,
|
|
exVolLiquid,
|
|
eyVolLiquid,
|
|
ezVolLiquid,
|
|
dx);
|
|
|
|
// Set up A matrix and rhs. Solve with CG. Update velocity grid.
|
|
solveViscosity(flags,
|
|
vel,
|
|
cVolLiquid,
|
|
uVolLiquid,
|
|
vVolLiquid,
|
|
wVolLiquid,
|
|
exVolLiquid,
|
|
eyVolLiquid,
|
|
ezVolLiquid,
|
|
viscosity,
|
|
dt,
|
|
dx,
|
|
cgAccuracy,
|
|
cgMaxIterFac);
|
|
}
|
|
static PyObject *_W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
|
|
{
|
|
try {
|
|
PbArgs _args(_linargs, _kwds);
|
|
FluidSolver *parent = _args.obtainParent();
|
|
bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
|
|
pbPreparePlugin(parent, "applyViscosity", !noTiming);
|
|
PyObject *_retval = nullptr;
|
|
{
|
|
ArgLocker _lock;
|
|
const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
|
|
const Grid<Real> &phi = *_args.getPtr<Grid<Real>>("phi", 1, &_lock);
|
|
MACGrid &vel = *_args.getPtr<MACGrid>("vel", 2, &_lock);
|
|
Grid<Real> &volumes = *_args.getPtr<Grid<Real>>("volumes", 3, &_lock);
|
|
Grid<Real> &viscosity = *_args.getPtr<Grid<Real>>("viscosity", 4, &_lock);
|
|
const Real cgAccuracy = _args.getOpt<Real>("cgAccuracy", 5, 1e-9, &_lock);
|
|
const Real cgMaxIterFac = _args.getOpt<Real>("cgMaxIterFac", 6, 1.5, &_lock);
|
|
_retval = getPyNone();
|
|
applyViscosity(flags, phi, vel, volumes, viscosity, cgAccuracy, cgMaxIterFac);
|
|
_args.check();
|
|
}
|
|
pbFinalizePlugin(parent, "applyViscosity", !noTiming);
|
|
return _retval;
|
|
}
|
|
catch (std::exception &e) {
|
|
pbSetError("applyViscosity", e.what());
|
|
return 0;
|
|
}
|
|
}
|
|
static const Pb::Register _RP_applyViscosity("", "applyViscosity", _W_0);
|
|
extern "C" {
|
|
void PbRegister_applyViscosity()
|
|
{
|
|
KEEP_UNUSED(_RP_applyViscosity);
|
|
}
|
|
}
|
|
|
|
} // namespace Manta
|
|
|
|
#if ENABLE_PARALLEL == 1
|
|
|
|
# undef parallel_block
|
|
# undef do_parallel
|
|
# undef do_end
|
|
# undef block_end
|
|
# undef parallel_for
|
|
# undef parallel_end
|
|
|
|
#endif
|