ASL
0.1.7
Advanced Simulation Library
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;
}
Generated by
1.8.15