Page MenuHome
Differential D3850 Diff 17619 intern/mantaflow/intern/manta_develop/preprocessed/omp/plugin/kepsilon.cpp

Changeset View

Standalone View

View Options

intern/mantaflow/intern/manta_develop/preprocessed/omp/plugin/kepsilon.cpp

  • This file was added.
// 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
*
* Turbulence modeling plugins
*
******************************************************************************/
#include "grid.h"
#include "commonkernels.h"
#include "vortexsheet.h"
#include "conjugategrad.h"
using namespace std;
namespace Manta {
// k-epsilon model constants
const Real keCmu = 0.09;
const Real keC1 = 1.44;
const Real keC2 = 1.92;
const Real keS1 = 1.0;
const Real keS2 = 1.3;
// k-epsilon limiters
const Real keU0 = 1.0;
const Real keImin = 2e-3;
const Real keImax = 1.0;
const Real keNuMin = 1e-3;
const Real keNuMax = 5.0;
//! clamp k and epsilon to limits
struct KnTurbulenceClamp : public KernelBase {
KnTurbulenceClamp(Grid<Real>& kgrid, Grid<Real>& egrid, Real minK, Real maxK, Real minNu, Real maxNu) : KernelBase(&kgrid,0) ,kgrid(kgrid),egrid(egrid),minK(minK),maxK(maxK),minNu(minNu),maxNu(maxNu) {
runMessage(); run(); }
inline void op(IndexInt idx, Grid<Real>& kgrid, Grid<Real>& egrid, Real minK, Real maxK, Real minNu, Real maxNu ) {
Real eps = egrid[idx];
Real ke = clamp(kgrid[idx],minK,maxK);
Real nu = keCmu*square(ke)/eps;
if (nu > maxNu)
eps = keCmu*square(ke)/maxNu;
if (nu < minNu)
eps = keCmu*square(ke)/minNu;
kgrid[idx] = ke;
egrid[idx] = eps;
} inline Grid<Real>& getArg0() {
return kgrid; }
typedef Grid<Real> type0;inline Grid<Real>& getArg1() {
return egrid; }
typedef Grid<Real> type1;inline Real& getArg2() {
return minK; }
typedef Real type2;inline Real& getArg3() {
return maxK; }
typedef Real type3;inline Real& getArg4() {
return minNu; }
typedef Real type4;inline Real& getArg5() {
return maxNu; }
typedef Real type5; void runMessage() { debMsg("Executing kernel KnTurbulenceClamp ", 3); debMsg("Kernel range" << " x "<< maxX << " y "<< maxY << " z "<< minZ<<" - "<< maxZ << " " , 4); }; void run() {
const IndexInt _sz = size;
#pragma omp parallel
{
#pragma omp for
for (IndexInt i = 0; i < _sz; i++) op(i,kgrid,egrid,minK,maxK,minNu,maxNu); }
}
Grid<Real>& kgrid; Grid<Real>& egrid; Real minK; Real maxK; Real minNu; Real maxNu; }
;
//! Compute k-epsilon production term P = 2*nu_T*sum_ij(Sij^2) and the turbulent viscosity nu_T=C_mu*k^2/eps
struct KnComputeProduction : public KernelBase {
KnComputeProduction(const MACGrid& vel, const Grid<Vec3>& velCenter, const Grid<Real>& ke, const Grid<Real>& eps, Grid<Real>& prod, Grid<Real>& nuT, Grid<Real>* strain, Real pscale = 1.0f) : KernelBase(&vel,1) ,vel(vel),velCenter(velCenter),ke(ke),eps(eps),prod(prod),nuT(nuT),strain(strain),pscale(pscale) {
runMessage(); run(); }
inline void op(int i, int j, int k, const MACGrid& vel, const Grid<Vec3>& velCenter, const Grid<Real>& ke, const Grid<Real>& eps, Grid<Real>& prod, Grid<Real>& nuT, Grid<Real>* strain, Real pscale = 1.0f ) {
Real curEps = eps(i,j,k);
if (curEps > 0) {
// turbulent viscosity: nu_T = C_mu * k^2/eps
Real curNu = keCmu * square(ke(i,j,k)) / curEps;
// compute Sij = 1/2 * (dU_i/dx_j + dU_j/dx_i)
Vec3 diag = Vec3(vel(i+1,j,k).x, vel(i,j+1,k).y, vel(i,j,k+1).z) - vel(i,j,k);
Vec3 ux = 0.5*(velCenter(i+1,j,k)-velCenter(i-1,j,k));
Vec3 uy = 0.5*(velCenter(i,j+1,k)-velCenter(i,j-1,k));
Vec3 uz = 0.5*(velCenter(i,j,k+1)-velCenter(i,j,k-1));
Real S12 = 0.5*(ux.y+uy.x);
Real S13 = 0.5*(ux.z+uz.x);
Real S23 = 0.5*(uy.z+uz.y);
Real S2 = square(diag.x) + square(diag.y) + square(diag.z) +
2.0*square(S12) + 2.0*square(S13) + 2.0*square(S23);
// P = 2*nu_T*sum_ij(Sij^2)
prod(i,j,k) = 2.0 * curNu * S2 * pscale;
nuT(i,j,k) = curNu;
if (strain) (*strain)(i,j,k) = sqrt(S2);
}
else {
prod(i,j,k) = 0;
nuT(i,j,k) = 0;
if (strain) (*strain)(i,j,k) = 0;
}
} inline const MACGrid& getArg0() {
return vel; }
typedef MACGrid type0;inline const Grid<Vec3>& getArg1() {
return velCenter; }
typedef Grid<Vec3> type1;inline const Grid<Real>& getArg2() {
return ke; }
typedef Grid<Real> type2;inline const Grid<Real>& getArg3() {
return eps; }
typedef Grid<Real> type3;inline Grid<Real>& getArg4() {
return prod; }
typedef Grid<Real> type4;inline Grid<Real>& getArg5() {
return nuT; }
typedef Grid<Real> type5;inline Grid<Real>* getArg6() {
return strain; }
typedef Grid<Real> type6;inline Real& getArg7() {
return pscale; }
typedef Real type7; void runMessage() { debMsg("Executing kernel KnComputeProduction ", 3); debMsg("Kernel range" << " x "<< maxX << " y "<< maxY << " z "<< minZ<<" - "<< maxZ << " " , 4); }; void run() {
const int _maxX = maxX; const int _maxY = maxY; if (maxZ > 1) {
#pragma omp parallel
{
#pragma omp for
for (int k=minZ; k < maxZ; k++) for (int j=1; j < _maxY; j++) for (int i=1; i < _maxX; i++) op(i,j,k,vel,velCenter,ke,eps,prod,nuT,strain,pscale); }
}
else {
const int k=0;
#pragma omp parallel
{
#pragma omp for
for (int j=1; j < _maxY; j++) for (int i=1; i < _maxX; i++) op(i,j,k,vel,velCenter,ke,eps,prod,nuT,strain,pscale); }
}
}
const MACGrid& vel; const Grid<Vec3>& velCenter; const Grid<Real>& ke; const Grid<Real>& eps; Grid<Real>& prod; Grid<Real>& nuT; Grid<Real>* strain; Real pscale; }
;
//! Compute k-epsilon production term P = 2*nu_T*sum_ij(Sij^2) and the turbulent viscosity nu_T=C_mu*k^2/eps
void KEpsilonComputeProduction(const MACGrid& vel, Grid<Real>& k, Grid<Real>& eps, Grid<Real>& prod, Grid<Real>& nuT, Grid<Real>* strain=0, Real pscale = 1.0f) {
// get centered velocity grid
Grid<Vec3> vcenter(k.getParent());
GetCentered(vcenter, vel);
FillInBoundary(vcenter,1);
// compute limits
const Real minK = 1.5*square(keU0)*square(keImin);
const Real maxK = 1.5*square(keU0)*square(keImax);
KnTurbulenceClamp(k, eps, minK, maxK, keNuMin, keNuMax);
KnComputeProduction(vel, vcenter, k, eps, prod, nuT, strain, pscale);
} 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, "KEpsilonComputeProduction" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; const MACGrid& vel = *_args.getPtr<MACGrid >("vel",0,&_lock); Grid<Real>& k = *_args.getPtr<Grid<Real> >("k",1,&_lock); Grid<Real>& eps = *_args.getPtr<Grid<Real> >("eps",2,&_lock); Grid<Real>& prod = *_args.getPtr<Grid<Real> >("prod",3,&_lock); Grid<Real>& nuT = *_args.getPtr<Grid<Real> >("nuT",4,&_lock); Grid<Real>* strain = _args.getPtrOpt<Grid<Real> >("strain",5,0,&_lock); Real pscale = _args.getOpt<Real >("pscale",6,1.0f,&_lock); _retval = getPyNone(); KEpsilonComputeProduction(vel,k,eps,prod,nuT,strain,pscale); _args.check(); }
pbFinalizePlugin(parent,"KEpsilonComputeProduction", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("KEpsilonComputeProduction",e.what()); return 0; }
}
static const Pb::Register _RP_KEpsilonComputeProduction ("","KEpsilonComputeProduction",_W_0);
extern "C" {
void PbRegister_KEpsilonComputeProduction() {
KEEP_UNUSED(_RP_KEpsilonComputeProduction); }
}
//! Integrate source terms of k-epsilon equation
struct KnAddTurbulenceSource : public KernelBase {
KnAddTurbulenceSource(Grid<Real>& kgrid, Grid<Real>& egrid, const Grid<Real>& pgrid, Real dt) : KernelBase(&kgrid,0) ,kgrid(kgrid),egrid(egrid),pgrid(pgrid),dt(dt) {
runMessage(); run(); }
inline void op(IndexInt idx, Grid<Real>& kgrid, Grid<Real>& egrid, const Grid<Real>& pgrid, Real dt ) {
Real eps = egrid[idx], prod = pgrid[idx], ke = kgrid[idx];
if (ke <= 0) ke = 1e-3; // pre-clamp to avoid nan
Real newK = ke + dt*(prod - eps);
Real newEps = eps + dt*(prod * keC1 - eps * keC2) * (eps / ke);
if (newEps <= 0) newEps = 1e-4; // pre-clamp to avoid nan
kgrid[idx] = newK;
egrid[idx] = newEps;
} inline Grid<Real>& getArg0() {
return kgrid; }
typedef Grid<Real> type0;inline Grid<Real>& getArg1() {
return egrid; }
typedef Grid<Real> type1;inline const Grid<Real>& getArg2() {
return pgrid; }
typedef Grid<Real> type2;inline Real& getArg3() {
return dt; }
typedef Real type3; void runMessage() { debMsg("Executing kernel KnAddTurbulenceSource ", 3); debMsg("Kernel range" << " x "<< maxX << " y "<< maxY << " z "<< minZ<<" - "<< maxZ << " " , 4); }; void run() {
const IndexInt _sz = size;
#pragma omp parallel
{
#pragma omp for
for (IndexInt i = 0; i < _sz; i++) op(i,kgrid,egrid,pgrid,dt); }
}
Grid<Real>& kgrid; Grid<Real>& egrid; const Grid<Real>& pgrid; Real dt; }
;
//! Integrate source terms of k-epsilon equation
void KEpsilonSources(Grid<Real>& k, Grid<Real>& eps, Grid<Real>& prod) {
Real dt = k.getParent()->getDt();
KnAddTurbulenceSource(k, eps, prod, dt);
// compute limits
const Real minK = 1.5*square(keU0)*square(keImin);
const Real maxK = 1.5*square(keU0)*square(keImax);
KnTurbulenceClamp(k, eps, minK, maxK, keNuMin, keNuMax);
} static PyObject* _W_1 (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, "KEpsilonSources" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; Grid<Real>& k = *_args.getPtr<Grid<Real> >("k",0,&_lock); Grid<Real>& eps = *_args.getPtr<Grid<Real> >("eps",1,&_lock); Grid<Real>& prod = *_args.getPtr<Grid<Real> >("prod",2,&_lock); _retval = getPyNone(); KEpsilonSources(k,eps,prod); _args.check(); }
pbFinalizePlugin(parent,"KEpsilonSources", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("KEpsilonSources",e.what()); return 0; }
}
static const Pb::Register _RP_KEpsilonSources ("","KEpsilonSources",_W_1);
extern "C" {
void PbRegister_KEpsilonSources() {
KEEP_UNUSED(_RP_KEpsilonSources); }
}
//! Initialize the domain or boundary conditions
void KEpsilonBcs(const FlagGrid& flags, Grid<Real>& k, Grid<Real>& eps, Real intensity, Real nu, bool fillArea) {
// compute limits
const Real vk = 1.5*square(keU0)*square(intensity);
const Real ve = keCmu*square(vk) / nu;
FOR_IDX(k) {
if (fillArea || flags.isObstacle(idx)) {
k[idx] = vk;
eps[idx] = ve;
}
}
} static PyObject* _W_2 (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, "KEpsilonBcs" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; const FlagGrid& flags = *_args.getPtr<FlagGrid >("flags",0,&_lock); Grid<Real>& k = *_args.getPtr<Grid<Real> >("k",1,&_lock); Grid<Real>& eps = *_args.getPtr<Grid<Real> >("eps",2,&_lock); Real intensity = _args.get<Real >("intensity",3,&_lock); Real nu = _args.get<Real >("nu",4,&_lock); bool fillArea = _args.get<bool >("fillArea",5,&_lock); _retval = getPyNone(); KEpsilonBcs(flags,k,eps,intensity,nu,fillArea); _args.check(); }
pbFinalizePlugin(parent,"KEpsilonBcs", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("KEpsilonBcs",e.what()); return 0; }
}
static const Pb::Register _RP_KEpsilonBcs ("","KEpsilonBcs",_W_2);
extern "C" {
void PbRegister_KEpsilonBcs() {
KEEP_UNUSED(_RP_KEpsilonBcs); }
}
//! Gradient diffusion smoothing. Not unconditionally stable -- should probably do substepping etc.
void ApplyGradDiff(const Grid<Real>& grid, Grid<Real>& res, const Grid<Real>& nu, Real dt, Real sigma) {
// should do this (but requires better boundary handling)
/*MACGrid grad(grid.getParent());
GradientOpMAC(grad, grid);
grad *= nu;
DivergenceOpMAC(res, grad);
res *= dt/sigma; */
LaplaceOp(res, grid);
res *= nu;
res *= dt/sigma;
}
//! Compute k-epsilon turbulent viscosity
void KEpsilonGradientDiffusion(Grid<Real>& k, Grid<Real>& eps, Grid<Real>& nuT, Real sigmaU=4.0, MACGrid* vel=0) {
Real dt = k.getParent()->getDt();
Grid<Real> res(k.getParent());
// gradient diffusion of k
ApplyGradDiff(k, res, nuT, dt, keS1);
k += res;
// gradient diffusion of epsilon
ApplyGradDiff(eps, res, nuT, dt, keS2);
eps += res;
// gradient diffusion of velocity
if (vel) {
Grid<Real> vc(k.getParent());
for (int c=0; c<3; c++) {
GetComponent(*vel, vc, c);
ApplyGradDiff(vc, res, nuT, dt, sigmaU);
vc += res;
SetComponent(*vel, vc, c);
}
}
} static PyObject* _W_3 (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, "KEpsilonGradientDiffusion" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; Grid<Real>& k = *_args.getPtr<Grid<Real> >("k",0,&_lock); Grid<Real>& eps = *_args.getPtr<Grid<Real> >("eps",1,&_lock); Grid<Real>& nuT = *_args.getPtr<Grid<Real> >("nuT",2,&_lock); Real sigmaU = _args.getOpt<Real >("sigmaU",3,4.0,&_lock); MACGrid* vel = _args.getPtrOpt<MACGrid >("vel",4,0,&_lock); _retval = getPyNone(); KEpsilonGradientDiffusion(k,eps,nuT,sigmaU,vel); _args.check(); }
pbFinalizePlugin(parent,"KEpsilonGradientDiffusion", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("KEpsilonGradientDiffusion",e.what()); return 0; }
}
static const Pb::Register _RP_KEpsilonGradientDiffusion ("","KEpsilonGradientDiffusion",_W_3);
extern "C" {
void PbRegister_KEpsilonGradientDiffusion() {
KEEP_UNUSED(_RP_KEpsilonGradientDiffusion); }
}
} // namespace