flowKDPGrowth.cc

Example: flowKDPGrowth

/*
* Advanced Simulation Library <http://asl.org.il>
*
* Copyright 2015 Avtech Scientific <http://avtechscientific.com>
*
*
* This file is part of Advanced Simulation Library (ASL).
*
* ASL is free software: you can redistribute it and/or modify it
* under the terms of the GNU Affero General Public License as
* published by the Free Software Foundation, version 3 of the License.
*
* ASL is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with ASL. If not, see <http://www.gnu.org/licenses/>.
*
*/
#include <utilities/aslParametersManager.h>
#include <math/aslTemplates.h>
#include <aslGeomInc.h>
#include <math/aslPositionFunction.h>
#include <aslDataInc.h>
#include <acl/aclGenerators.h>
#include <acl/aclMath/aclVectorOfElements.h>
#include <writers/aslVTKFormatWriters.h>
#include <num/aslLBGK.h>
#include <num/aslLBGKBC.h>
#include <num/aslBasicBC.h>
#include <num/aslCrystalGrowthBC.h>
#include <num/aslFDAdvectionDiffusion.h>
#include <utilities/aslTimer.h>
using asl::AVec;
using asl::makeAVec;
asl::SPDistanceFunction generateBath(asl::Block & bl)
{
// double hBath(2.);
double rBath(1.);
double dx(bl.dx);
asl::AVec<int> size(bl.getSize());
asl::AVec<>center(.5*dx*AVec<>(size));
auto bath(-(generateDFCylinderInf(rBath, makeAVec(0.,0.,1.), dx*AVec<>(size)*.5) &
generateDFPlane(makeAVec(0.,0.,1.), center*1.99) &
generateDFPlane(makeAVec(0.,0.,-1.), center*0.)));
return normalize(bath, dx);
}
asl::SPDistanceFunction generatePlatform(asl::Block & bl)
{
double rDisk(.9);
double hDisk(0.1);
double rAxis(0.05);
double hAxis(.5);
double wPillar(.2);
double dPillar(.1);
double dx(bl.dx);
asl::AVec<int> size(bl.getSize());
asl::AVec<>center(.5*dx*AVec<>(size));
vector<asl::AVec<>> pillar1{makeAVec(wPillar*.5, dPillar*.5,0.),
makeAVec(-wPillar*.5, dPillar*.5,0.),
makeAVec(-wPillar*.5, -dPillar*.5,0.),
makeAVec(wPillar*.5, -dPillar*.5,0.)};
vector<asl::AVec<>> pillar2{makeAVec(dPillar*.5, wPillar*.5,0.),
makeAVec(-dPillar*.5, wPillar*.5,0.),
makeAVec(-dPillar*.5, -wPillar*.5,0.),
makeAVec(dPillar*.5, -wPillar*.5,0.)};
vector<asl::AVec<>> pillarC{makeAVec(center[0]+rDisk-dPillar*.5, center[1], 0.),
makeAVec(center[0]-rDisk+dPillar*.5, center[1], 0.),
makeAVec(center[0], center[1]+rDisk-dPillar*.5,0.),
makeAVec(center[0], center[1]-rDisk+dPillar*.5,0.)};
vector<vector<asl::AVec<>>> pillarsPoints(4);
for(unsigned int i(0); i<4; ++i)
pillarsPoints[i].resize(4);
for(unsigned int i(0); i<4; ++i)
{
pillarsPoints[0][i] = pillar2[i] + pillarC[0];
pillarsPoints[1][i] = pillar2[i] + pillarC[1];
pillarsPoints[2][i] = pillar1[i] + pillarC[2];
pillarsPoints[3][i] = pillar1[i] + pillarC[3];
}
auto diskBottom(generateDFCylinder(rDisk,
makeAVec(0., 0., hDisk),
makeAVec(center[0], center[1], .5*hDisk)));
auto diskTop(generateDFCylinder(rDisk,
makeAVec(0., 0., hDisk),
makeAVec(center[0], center[1], -.5*hDisk - hAxis + dx*size[2])));
auto axis(generateDFCylinder(rAxis,
makeAVec(0., 0., hAxis+hDisk*.5),
makeAVec(center[0], center[1], - .5*hAxis - hDisk*.25 + dx*size[2])));
auto dfPillar1(generateDFConvexPolygonPrism(pillarsPoints[0]));
auto dfPillar2(generateDFConvexPolygonPrism(pillarsPoints[1]));
auto dfPillar3(generateDFConvexPolygonPrism(pillarsPoints[2]));
auto dfPillar4(generateDFConvexPolygonPrism(pillarsPoints[3]));
auto dfPillars((dfPillar1 | dfPillar2 | dfPillar3 | dfPillar4) &
generateDFPlane(makeAVec(0.,0.,-1.), makeAVec(0.,0.,.5*hDisk)) &
generateDFPlane(makeAVec(0.,0.,1.), makeAVec(0.,0.,-.5*hDisk - hAxis + dx*size[2])));
return normalize(diskBottom | diskTop | axis | dfPillars, dx);
}
asl::SPDistanceFunction generateCrystal(asl::Block & bl)
{
double aCrystal(.5);
double hCrystalBase(.5);
double hCrystalPyramid(.5);
double hDisk(0.1);
double dx(bl.dx);
asl::AVec<int> size(bl.getSize());
asl::AVec<>center(.5*dx*AVec<>(size));
auto crystalB(asl::generateDFConvexPolygonPrism({center+makeAVec( aCrystal, aCrystal,0.),
center+makeAVec(-aCrystal, aCrystal,0.),
center+makeAVec(-aCrystal, -aCrystal,0.),
center+makeAVec( aCrystal, -aCrystal,0.)}) &
generateDFPlane(makeAVec(0.,0.,-1.), makeAVec(0.,0., hDisk-.001)) &
generateDFPlane(makeAVec(0.,0., 1.), makeAVec(0.,0., hDisk + hCrystalBase)));
auto cCrPyrBase(makeAVec(center[0],center[1],hDisk+hCrystalBase-.01));
auto crystalT(asl::generateDFConvexPolygonPyramid({cCrPyrBase+makeAVec( aCrystal, aCrystal,0.),
cCrPyrBase+makeAVec(-aCrystal, aCrystal,0.),
cCrPyrBase+makeAVec(-aCrystal, -aCrystal,0.),
cCrPyrBase+makeAVec( aCrystal, -aCrystal,0.)},
cCrPyrBase+makeAVec(0.,0.,hCrystalPyramid)));
return normalize(crystalB | crystalT, dx);
// return crystalB | crystalT;
}
double getWRotation(double t)
{
double tPeriod(128);
double wMax(6.*3.14*2./60.);
double tPlato(tPeriod * .25);
double tAcceleration(tPeriod * .1);
double tStop(tPeriod * .05);
double intPart;
double tRel(modf(t/tPeriod, &intPart));
double x(0);
if(tRel<=tAcceleration)
x = tRel / tAcceleration;
if(tRel>tAcceleration && tRel<=tAcceleration+tPlato)
x = 1.;
if(tRel>tAcceleration+tPlato && tRel<=2.*tAcceleration+tPlato)
x = (2.*tAcceleration + tPlato - tRel) / tAcceleration;
if(tRel>2.*tAcceleration+tPlato && tRel<=2.*tAcceleration+tPlato+tStop)
x = 0;
if(tRel>2.*tAcceleration+tPlato+tStop && tRel<=3.*tAcceleration+tPlato+tStop)
x = -(tRel-2.*tAcceleration-tPlato-tStop) / tAcceleration;
if(tRel>3.*tAcceleration+tPlato+tStop && tRel<=3.*tAcceleration+2.*tPlato+tStop)
x = -1.;
if(tRel>3.*tAcceleration+2.*tPlato+tStop && tRel<=4.*tAcceleration+2.*tPlato+tStop)
x = -(4.*tAcceleration+2.*tPlato+tStop-tRel)/tAcceleration;
if(tRel>4.*tAcceleration+2.*tPlato+tStop)
x = 0;
return wMax*x;
// flux = -9.32e-5*(1.170 - c); c_0=0.326 ceq=0.267
}
typedef float FlT;
//typedef double FlT;
typedef asl::UValue<double> Param;
using asl::AVec;
using asl::makeAVec;
int main(int argc, char* argv[])
{
// Optionally add appParamsManager to be able to manipulate at least
// hardware parameters(platform/device) through command line/parameters file
asl::ApplicationParametersManager appParamsManager("flowKDPGrowth",
"1.0");
appParamsManager.load(argc, argv);
Param dx(.02);
Param dt(0.8e-2);
Param nu(1e-2);
Param difC(1e-2/300.);
// Param w(48.*3.14*2./60.);
// Angular velocity
Param w(6.*3.14*2./60.);
Param nuNum(nu.v()*dt.v()/dx.v()/dx.v());
Param difCNum(difC.v()*dt.v()/dx.v()/dx.v());
// Angular velocity in one iteration
Param wNum(w.v()*dt.v());
Param c0(0.326);
AVec<int> size(asl::makeAVec(105.,105.,100.));
AVec<> gSize(dx.v()*AVec<>(size));
std::cout << "Data initialization...";
auto templ(&asl::d3q19());
asl::Block block(size,dx.v());
auto bathMap(asl::generateDataContainerACL_SP<FlT>(block, 1, 1u));
asl::initData(bathMap, generateBath(block));
auto platformCrysMap(asl::generateDataContainerACL_SP<FlT>(block, 1, 1u));
asl::initData(platformCrysMap, generatePlatform(block) | generateCrystal(block));
auto bathPlatformMap(asl::generateDataContainerACL_SP<FlT>(block, 1, 1u));
asl::initData(bathPlatformMap, generateBath(block) | generatePlatform(block));
auto bathPlatformCrystalMap(asl::generateDataContainerACL_SP<FlT>(block, 1, 1u));
asl::initData(bathPlatformCrystalMap, generateBath(block) | generatePlatform(block) | generateCrystal(block));
auto crystalMap(asl::generateDataContainerACL_SP<FlT>(block, 1, 1u));
asl::initData(crystalMap, generateCrystal(block));
auto cField(asl::generateDataContainerACL_SP<FlT>(block, 1, 1u));
asl::initData(cField, c0.v());
std::cout << "Finished" << endl;
std::cout << "Numerics initialization...";
asl::SPLBGK lbgk(new asl::LBGK(block,
acl::generateVEConstant(FlT(nuNum.v())),
templ));
// Set angular velocity in lbgk
lbgk->setOmega(acl::generateVEConstant(makeAVec(0.,0.,wNum.v())));
lbgk->init();
asl::SPLBGKUtilities lbgkUtil(new asl::LBGKUtilities(lbgk));
lbgkUtil->initF(acl::generateVEConstant(.0,.0,.0));
auto nmDif(asl::generateFDAdvectionDiffusion(cField,
difCNum.v(),
lbgk->getVelocity(),
templ,
true));
nmDif->init();
std::vector<asl::SPNumMethod> bc;
std::vector<asl::SPNumMethod> bcV;
std::vector<asl::SPNumMethod> bcDif;
// Position Function Angular Velocity Field
auto vfBath(asl::generatePFRotationField(makeAVec(0.,0., wNum.v()/dx.v()), .5*gSize));
// Boundary condition
bc.push_back(generateBCVelocity(lbgk, vfBath, bathMap));
// Boundary condition for visualization
bcV.push_back(generateBCVelocityVel(lbgk, vfBath, bathMap));
bc.push_back(generateBCNoSlip(lbgk, platformCrysMap));
bcV.push_back(generateBCNoSlipVel(lbgk, platformCrysMap));
bcDif.push_back(generateBCConstantGradient(cField,
0.,
bathPlatformMap,
bathPlatformCrystalMap,
templ));
bcDif.push_back(generateBCLinearGrowth2(cField, 1.17,
-9.32e-6/difC.v()*dx.v(),
crystalMap,
bathPlatformCrystalMap,
templ));
// bcDif.push_back(generateBCConstantGradient2(cField, .1, crystalMap, bathPlatformCrystalMap, templ));
initAll(bc);
initAll(bcV);
initAll(bcDif);
std::cout << "Finished" << endl;
std::cout << "Computing...";
asl::Timer timer;
asl::Timer timerBC;
asl::WriterVTKXML writer("flowKDPGrowthRes");
writer.addScalars("mapBath", *bathMap);
writer.addScalars("mapPlatformCrys", *platformCrysMap);
writer.addScalars("mapBathPlatformCrystal", *bathPlatformCrystalMap);
writer.addScalars("mapCrys", *crystalMap);
writer.addScalars("rho", *lbgk->getRho());
writer.addScalars("c", *cField);
writer.addVector("v", *lbgk->getVelocity());
executeAll(bcV);
executeAll(bc);
executeAll(bcDif);
writer.write();
timer.start();
timerBC.reset();
for (unsigned int i(0); i <= 8001 ; ++i)
{
lbgk->execute();
timerBC.start();
executeAll(bcV);
executeAll(bc);
timerBC.stop();
nmDif->execute();
timerBC.start();
executeAll(bcDif);
timerBC.stop();
if (!(i%2000))
{
cout << i << endl;
writer.write();
}
}
timer.stop();
cout << "Finished" << endl;
cout << "Computation statistic:" << endl;
cout << "Real Time = " << timer.realTime() << "; Processor Time = "
<< timer.processorTime() << "; Processor Load = "
<< timer.processorLoad() * 100 << "%" << endl;
return 0;
}
asl::AVec
Definition: aslVectorsDynamicLength.h:40
asl::AVec::normalize
const AVec normalize(const AVec< T > &a)
Definition: aslVectorsDynamicLengthOperations.h:152
asl::AVec::makeAVec
AVec< T > makeAVec(T a1)
Definition: aslVectorsDynamicLength.h:176
asl::ApplicationParametersManager
Definition: aslParametersManager.h:159
asl::Block
Definition: aslBlocks.h:57
asl::Block::dx
double dx
Definition: aslBlocks.h:66
asl::Block::getSize
const DV & getSize() const
Definition: aslBlocks.h:208
asl::LBGK
Numerical method for fluid flow.
Definition: aslLBGK.h:78
asl::LBGKUtilities
contains different kernels for preprocessing and posprocessing of data used by LBGK
Definition: aslLBGK.h:138
asl::Timer
Definition: aslTimer.h:34
asl::Timer::realTime
const double realTime() const
Definition: aslTimer.h:45
asl::Timer::stop
void stop()
Definition: aslTimer.h:44
asl::Timer::processorTime
const double processorTime() const
Definition: aslTimer.h:46
asl::Timer::start
void start()
Definition: aslTimer.h:43
asl::Timer::reset
void reset()
Definition: aslTimer.h:48
asl::Timer::processorLoad
const double processorLoad() const
Definition: aslTimer.h:47
asl::UValue< double >
asl::WriterVTKXML
Definition: aslVTKFormatWriters.h:42
main
int main(int argc, char *argv[])
Definition: flowKDPGrowth.cc:200
FlT
float FlT
Definition: flowKDPGrowth.cc:192
generateCrystal
asl::SPDistanceFunction generateCrystal(asl::Block &bl)
Definition: flowKDPGrowth.cc:130
generateBath
asl::SPDistanceFunction generateBath(asl::Block &bl)
Definition: flowKDPGrowth.cc:47
Param
asl::UValue< double > Param
Definition: flowKDPGrowth.cc:194
getWRotation
double getWRotation(double t)
Definition: flowKDPGrowth.cc:159
generatePlatform
asl::SPDistanceFunction generatePlatform(asl::Block &bl)
Definition: flowKDPGrowth.cc:68
asl::dx
acl::VectorOfElements dx(const TemplateVE &a)
differential operator
asl::generateBCConstantGradient
SPBCond generateBCConstantGradient(SPAbstractDataWithGhostNodes d, double v, const VectorTemplate *const t, const std::vector< SlicesNames > &sl)
Bondary condition that makes fixed gradient <>
asl::generateDFConvexPolygonPrism
SPDistanceFunction generateDFConvexPolygonPrism(std::vector< AVec< double >> points)
generates infinite prism with convex polygon at its base
asl::generateDFConvexPolygonPyramid
SPDistanceFunction generateDFConvexPolygonPyramid(std::vector< AVec< double >> points, AVec< double > a)
generates pyramid with convex polygon at its base and apex a
asl::generateDFCylinder
SPDistanceFunction generateDFCylinder(double r, const AVec< double > &l, const AVec< double > &c)
generates cylinder
asl::generateDFPlane
SPDistanceFunction generateDFPlane(const AVec< double > &n, const AVec< double > &p0)
asl::generateDFCylinderInf
SPDistanceFunction generateDFCylinderInf(double r, const AVec< double > &l, const AVec< double > &c)
generates infinite cylinder
asl::SPDistanceFunction
std::shared_ptr< DistanceFunction > SPDistanceFunction
Definition: aslGeomInc.h:44
asl::generateFDAdvectionDiffusion
SPFDAdvectionDiffusion generateFDAdvectionDiffusion(SPDataWithGhostNodesACLData c, double diffustionCoeff, SPAbstractDataWithGhostNodes v, const VectorTemplate *vt, bool compressibilityCorrection=false)
asl::generatePFRotationField
SPPositionFunction generatePFRotationField(const AVec< double > &axis, const AVec< double > &c)
asl::d3q19
const VectorTemplate & d3q19()
Vector template.
asl::generateBCVelocityVel
SPNumMethod generateBCVelocityVel(SPLBGK nm, SPPositionFunction v, SPAbstractDataWithGhostNodes map)
asl::generateBCNoSlipVel
SPNumMethod generateBCNoSlipVel(SPLBGK nmU, SPAbstractDataWithGhostNodes map)
asl::generateBCVelocity
SPNumMethod generateBCVelocity(SPLBGK nm, SPPositionFunction v, SPAbstractDataWithGhostNodes map)
asl::generateBCNoSlip
SPBCond generateBCNoSlip(SPLBGK nm, const std::vector< SlicesNames > &sl)
acl::generateVEConstant
VectorOfElements generateVEConstant(T a)
Generates VectorOfElements with 1 Element acl::Constant with value a.
asl::makeAVec
AVec< T > makeAVec(T a1)
Definition: aslVectorsDynamicLength.h:176
asl::generateBCLinearGrowth2
SPNumMethod generateBCLinearGrowth2(SPAbstractDataWithGhostNodes d, double cEq, double beta, SPAbstractDataWithGhostNodes map, const VectorTemplate *const t)
asl::initAll
void initAll(std::vector< T * > &v)
Definition: aslNumMethod.h:67
asl::SPLBGK
std::shared_ptr< LBGK > SPLBGK
Definition: aslLBGK.h:133
asl::initData
void initData(SPAbstractData d, double a)
asl::SPLBGKUtilities
std::shared_ptr< LBGKUtilities > SPLBGKUtilities
Definition: aslLBGK.h:161
asl::executeAll
void executeAll(std::vector< T * > &v)
Definition: aslNumMethod.h:55