ASL  0.1.7
Advanced Simulation Library
flowKDPGrowth.cc
Go to the documentation of this file.
1 /*
2  * Advanced Simulation Library <http://asl.org.il>
3  *
4  * Copyright 2015 Avtech Scientific <http://avtechscientific.com>
5  *
6  *
7  * This file is part of Advanced Simulation Library (ASL).
8  *
9  * ASL is free software: you can redistribute it and/or modify it
10  * under the terms of the GNU Affero General Public License as
11  * published by the Free Software Foundation, version 3 of the License.
12  *
13  * ASL is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU Affero General Public License for more details.
17  *
18  * You should have received a copy of the GNU Affero General Public License
19  * along with ASL. If not, see <http://www.gnu.org/licenses/>.
20  *
21  */
22 
23 
30 #include <math/aslTemplates.h>
31 #include <aslGeomInc.h>
33 #include <aslDataInc.h>
34 #include <acl/aclGenerators.h>
37 #include <num/aslLBGK.h>
38 #include <num/aslLBGKBC.h>
39 #include <num/aslBasicBC.h>
40 #include <num/aslCrystalGrowthBC.h>
42 #include <utilities/aslTimer.h>
43 
44 using asl::AVec;
45 using asl::makeAVec;
46 
48 {
49 
50 // double hBath(2.);
51  double rBath(1.);
52 
53  double dx(bl.dx);
54 
55  asl::AVec<int> size(bl.getSize());
56 
57  asl::AVec<>center(.5*dx*AVec<>(size));
58 
59 
60 
61  auto bath(-(generateDFCylinderInf(rBath, makeAVec(0.,0.,1.), dx*AVec<>(size)*.5) &
62  generateDFPlane(makeAVec(0.,0.,1.), center*1.99) &
63  generateDFPlane(makeAVec(0.,0.,-1.), center*0.)));
64 
65  return normalize(bath, dx);
66 }
67 
69 {
70  double rDisk(.9);
71  double hDisk(0.1);
72 
73  double rAxis(0.05);
74  double hAxis(.5);
75 
76  double wPillar(.2);
77  double dPillar(.1);
78 
79  double dx(bl.dx);
80  asl::AVec<int> size(bl.getSize());
81  asl::AVec<>center(.5*dx*AVec<>(size));
82 
83  vector<asl::AVec<>> pillar1{makeAVec(wPillar*.5, dPillar*.5,0.),
84  makeAVec(-wPillar*.5, dPillar*.5,0.),
85  makeAVec(-wPillar*.5, -dPillar*.5,0.),
86  makeAVec(wPillar*.5, -dPillar*.5,0.)};
87 
88  vector<asl::AVec<>> pillar2{makeAVec(dPillar*.5, wPillar*.5,0.),
89  makeAVec(-dPillar*.5, wPillar*.5,0.),
90  makeAVec(-dPillar*.5, -wPillar*.5,0.),
91  makeAVec(dPillar*.5, -wPillar*.5,0.)};
92 
93  vector<asl::AVec<>> pillarC{makeAVec(center[0]+rDisk-dPillar*.5, center[1], 0.),
94  makeAVec(center[0]-rDisk+dPillar*.5, center[1], 0.),
95  makeAVec(center[0], center[1]+rDisk-dPillar*.5,0.),
96  makeAVec(center[0], center[1]-rDisk+dPillar*.5,0.)};
97  vector<vector<asl::AVec<>>> pillarsPoints(4);
98  for(unsigned int i(0); i<4; ++i)
99  pillarsPoints[i].resize(4);
100 
101  for(unsigned int i(0); i<4; ++i)
102  {
103  pillarsPoints[0][i] = pillar2[i] + pillarC[0];
104  pillarsPoints[1][i] = pillar2[i] + pillarC[1];
105  pillarsPoints[2][i] = pillar1[i] + pillarC[2];
106  pillarsPoints[3][i] = pillar1[i] + pillarC[3];
107  }
108 
109 
110  auto diskBottom(generateDFCylinder(rDisk,
111  makeAVec(0., 0., hDisk),
112  makeAVec(center[0], center[1], .5*hDisk)));
113  auto diskTop(generateDFCylinder(rDisk,
114  makeAVec(0., 0., hDisk),
115  makeAVec(center[0], center[1], -.5*hDisk - hAxis + dx*size[2])));
116  auto axis(generateDFCylinder(rAxis,
117  makeAVec(0., 0., hAxis+hDisk*.5),
118  makeAVec(center[0], center[1], - .5*hAxis - hDisk*.25 + dx*size[2])));
119  auto dfPillar1(generateDFConvexPolygonPrism(pillarsPoints[0]));
120  auto dfPillar2(generateDFConvexPolygonPrism(pillarsPoints[1]));
121  auto dfPillar3(generateDFConvexPolygonPrism(pillarsPoints[2]));
122  auto dfPillar4(generateDFConvexPolygonPrism(pillarsPoints[3]));
123  auto dfPillars((dfPillar1 | dfPillar2 | dfPillar3 | dfPillar4) &
124  generateDFPlane(makeAVec(0.,0.,-1.), makeAVec(0.,0.,.5*hDisk)) &
125  generateDFPlane(makeAVec(0.,0.,1.), makeAVec(0.,0.,-.5*hDisk - hAxis + dx*size[2])));
126 
127  return normalize(diskBottom | diskTop | axis | dfPillars, dx);
128 }
129 
131 {
132 
133  double aCrystal(.5);
134  double hCrystalBase(.5);
135  double hCrystalPyramid(.5);
136 
137  double hDisk(0.1);
138 
139  double dx(bl.dx);
140  asl::AVec<int> size(bl.getSize());
141  asl::AVec<>center(.5*dx*AVec<>(size));
142 
143  auto crystalB(asl::generateDFConvexPolygonPrism({center+makeAVec( aCrystal, aCrystal,0.),
144  center+makeAVec(-aCrystal, aCrystal,0.),
145  center+makeAVec(-aCrystal, -aCrystal,0.),
146  center+makeAVec( aCrystal, -aCrystal,0.)}) &
147  generateDFPlane(makeAVec(0.,0.,-1.), makeAVec(0.,0., hDisk-.001)) &
148  generateDFPlane(makeAVec(0.,0., 1.), makeAVec(0.,0., hDisk + hCrystalBase)));
149  auto cCrPyrBase(makeAVec(center[0],center[1],hDisk+hCrystalBase-.01));
150  auto crystalT(asl::generateDFConvexPolygonPyramid({cCrPyrBase+makeAVec( aCrystal, aCrystal,0.),
151  cCrPyrBase+makeAVec(-aCrystal, aCrystal,0.),
152  cCrPyrBase+makeAVec(-aCrystal, -aCrystal,0.),
153  cCrPyrBase+makeAVec( aCrystal, -aCrystal,0.)},
154  cCrPyrBase+makeAVec(0.,0.,hCrystalPyramid)));
155  return normalize(crystalB | crystalT, dx);
156 // return crystalB | crystalT;
157 }
158 
159 double getWRotation(double t)
160 {
161  double tPeriod(128);
162  double wMax(6.*3.14*2./60.);
163  double tPlato(tPeriod * .25);
164  double tAcceleration(tPeriod * .1);
165  double tStop(tPeriod * .05);
166 
167  double intPart;
168  double tRel(modf(t/tPeriod, &intPart));
169  double x(0);
170  if(tRel<=tAcceleration)
171  x = tRel / tAcceleration;
172  if(tRel>tAcceleration && tRel<=tAcceleration+tPlato)
173  x = 1.;
174  if(tRel>tAcceleration+tPlato && tRel<=2.*tAcceleration+tPlato)
175  x = (2.*tAcceleration + tPlato - tRel) / tAcceleration;
176  if(tRel>2.*tAcceleration+tPlato && tRel<=2.*tAcceleration+tPlato+tStop)
177  x = 0;
178  if(tRel>2.*tAcceleration+tPlato+tStop && tRel<=3.*tAcceleration+tPlato+tStop)
179  x = -(tRel-2.*tAcceleration-tPlato-tStop) / tAcceleration;
180  if(tRel>3.*tAcceleration+tPlato+tStop && tRel<=3.*tAcceleration+2.*tPlato+tStop)
181  x = -1.;
182  if(tRel>3.*tAcceleration+2.*tPlato+tStop && tRel<=4.*tAcceleration+2.*tPlato+tStop)
183  x = -(4.*tAcceleration+2.*tPlato+tStop-tRel)/tAcceleration;
184  if(tRel>4.*tAcceleration+2.*tPlato+tStop)
185  x = 0;
186  return wMax*x;
187 
188 // flux = -9.32e-5*(1.170 - c); c_0=0.326 ceq=0.267
189 }
190 
191 
192 typedef float FlT;
193 //typedef double FlT;
195 
196 using asl::AVec;
197 using asl::makeAVec;
198 
199 
200 int main(int argc, char* argv[])
201 {
202  // Optionally add appParamsManager to be able to manipulate at least
203  // hardware parameters(platform/device) through command line/parameters file
204  asl::ApplicationParametersManager appParamsManager("flowKDPGrowth",
205  "1.0");
206  appParamsManager.load(argc, argv);
207 
208  Param dx(.02);
209  Param dt(0.8e-2);
210  Param nu(1e-2);
211  Param difC(1e-2/300.);
212 // Param w(48.*3.14*2./60.);
213  // Angular velocity
214  Param w(6.*3.14*2./60.);
215 
216  Param nuNum(nu.v()*dt.v()/dx.v()/dx.v());
217  Param difCNum(difC.v()*dt.v()/dx.v()/dx.v());
218 
219  // Angular velocity in one iteration
220  Param wNum(w.v()*dt.v());
221 
222  Param c0(0.326);
223 
224  AVec<int> size(asl::makeAVec(105.,105.,100.));
225 
226  AVec<> gSize(dx.v()*AVec<>(size));
227 
228  std::cout << "Data initialization...";
229 
230  auto templ(&asl::d3q19());
231  asl::Block block(size,dx.v());
232 
233  auto bathMap(asl::generateDataContainerACL_SP<FlT>(block, 1, 1u));
234  asl::initData(bathMap, generateBath(block));
235  auto platformCrysMap(asl::generateDataContainerACL_SP<FlT>(block, 1, 1u));
236  asl::initData(platformCrysMap, generatePlatform(block) | generateCrystal(block));
237  auto bathPlatformMap(asl::generateDataContainerACL_SP<FlT>(block, 1, 1u));
238  asl::initData(bathPlatformMap, generateBath(block) | generatePlatform(block));
239  auto bathPlatformCrystalMap(asl::generateDataContainerACL_SP<FlT>(block, 1, 1u));
240  asl::initData(bathPlatformCrystalMap, generateBath(block) | generatePlatform(block) | generateCrystal(block));
241  auto crystalMap(asl::generateDataContainerACL_SP<FlT>(block, 1, 1u));
242  asl::initData(crystalMap, generateCrystal(block));
243 
244  auto cField(asl::generateDataContainerACL_SP<FlT>(block, 1, 1u));
245  asl::initData(cField, c0.v());
246 
247  std::cout << "Finished" << endl;
248 
249  std::cout << "Numerics initialization...";
250 
251  asl::SPLBGK lbgk(new asl::LBGK(block,
252  acl::generateVEConstant(FlT(nuNum.v())),
253  templ));
254  // Set angular velocity in lbgk
255  lbgk->setOmega(acl::generateVEConstant(makeAVec(0.,0.,wNum.v())));
256  lbgk->init();
257  asl::SPLBGKUtilities lbgkUtil(new asl::LBGKUtilities(lbgk));
258  lbgkUtil->initF(acl::generateVEConstant(.0,.0,.0));
259 
260  auto nmDif(asl::generateFDAdvectionDiffusion(cField,
261  difCNum.v(),
262  lbgk->getVelocity(),
263  templ,
264  true));
265 
266  nmDif->init();
267  std::vector<asl::SPNumMethod> bc;
268  std::vector<asl::SPNumMethod> bcV;
269  std::vector<asl::SPNumMethod> bcDif;
270 
271  // Position Function Angular Velocity Field
272  auto vfBath(asl::generatePFRotationField(makeAVec(0.,0., wNum.v()/dx.v()), .5*gSize));
273  // Boundary condition
274  bc.push_back(generateBCVelocity(lbgk, vfBath, bathMap));
275  // Boundary condition for visualization
276  bcV.push_back(generateBCVelocityVel(lbgk, vfBath, bathMap));
277  bc.push_back(generateBCNoSlip(lbgk, platformCrysMap));
278  bcV.push_back(generateBCNoSlipVel(lbgk, platformCrysMap));
279  bcDif.push_back(generateBCConstantGradient(cField,
280  0.,
281  bathPlatformMap,
282  bathPlatformCrystalMap,
283  templ));
284  bcDif.push_back(generateBCLinearGrowth2(cField, 1.17,
285  -9.32e-6/difC.v()*dx.v(),
286  crystalMap,
287  bathPlatformCrystalMap,
288  templ));
289 // bcDif.push_back(generateBCConstantGradient2(cField, .1, crystalMap, bathPlatformCrystalMap, templ));
290  initAll(bc);
291  initAll(bcV);
292  initAll(bcDif);
293 
294  std::cout << "Finished" << endl;
295  std::cout << "Computing...";
296  asl::Timer timer;
297  asl::Timer timerBC;
298 
299  asl::WriterVTKXML writer("flowKDPGrowthRes");
300  writer.addScalars("mapBath", *bathMap);
301  writer.addScalars("mapPlatformCrys", *platformCrysMap);
302  writer.addScalars("mapBathPlatformCrystal", *bathPlatformCrystalMap);
303  writer.addScalars("mapCrys", *crystalMap);
304  writer.addScalars("rho", *lbgk->getRho());
305  writer.addScalars("c", *cField);
306  writer.addVector("v", *lbgk->getVelocity());
307 
308  executeAll(bcV);
309  executeAll(bc);
310  executeAll(bcDif);
311 
312  writer.write();
313 
314  timer.start();
315  timerBC.reset();
316  for (unsigned int i(0); i <= 8001 ; ++i)
317  {
318  lbgk->execute();
319  timerBC.start();
320  executeAll(bcV);
321  executeAll(bc);
322  timerBC.stop();
323  nmDif->execute();
324  timerBC.start();
325  executeAll(bcDif);
326  timerBC.stop();
327 
328  if (!(i%2000))
329  {
330  cout << i << endl;
331  writer.write();
332  }
333  }
334  timer.stop();
335 
336 
337  cout << "Finished" << endl;
338 
339  cout << "Computation statistic:" << endl;
340  cout << "Real Time = " << timer.realTime() << "; Processor Time = "
341  << timer.processorTime() << "; Processor Load = "
342  << timer.processorLoad() * 100 << "%" << endl;
343 
344  return 0;
345 }
asl::UValue< double > Param
const double realTime() const
Definition: aslTimer.h:45
std::shared_ptr< DistanceFunction > SPDistanceFunction
Definition: aslGeomInc.h:44
const DV & getSize() const
Definition: aslBlocks.h:208
const double processorTime() const
Definition: aslTimer.h:46
AVec< T > makeAVec(T a1)
Numerical method for fluid flow.
Definition: aslLBGK.h:77
const T & v() const
Definition: aslUValue.h:43
const AVec normalize(const AVec< T > &a)
void addVector(std::string name, AbstractData &data)
void initAll(std::vector< T * > &v)
Definition: aslNumMethod.h:67
AVec< T > makeAVec(T a1)
SPDistanceFunction generateDFConvexPolygonPrism(std::vector< AVec< double >> points)
generates infinite prism with convex polygon at its base
double getWRotation(double t)
double dx
Definition: aslBlocks.h:66
SPDistanceFunction generateDFCylinderInf(double r, const AVec< double > &l, const AVec< double > &c)
generates infinite cylinder
SPDistanceFunction generateDFPlane(const AVec< double > &n, const AVec< double > &p0)
std::shared_ptr< LBGKUtilities > SPLBGKUtilities
Definition: aslLBGK.h:161
std::shared_ptr< LBGK > SPLBGK
Definition: aslLBGK.h:133
float FlT
SPFDAdvectionDiffusion generateFDAdvectionDiffusion(SPDataWithGhostNodesACLData c, double diffustionCoeff, SPAbstractDataWithGhostNodes v, const VectorTemplate *vt, bool compressibilityCorrection=false)
acl::VectorOfElements dx(const TemplateVE &a)
differential operator
void reset()
Definition: aslTimer.h:48
SPDistanceFunction generateDFConvexPolygonPyramid(std::vector< AVec< double >> points, AVec< double > a)
generates pyramid with convex polygon at its base and apex a
asl::SPDistanceFunction generateBath(asl::Block &bl)
SPNumMethod generateBCVelocityVel(SPLBGK nm, SPPositionFunction v, SPAbstractDataWithGhostNodes map)
void initData(SPAbstractData d, double a)
SPNumMethod generateBCLinearGrowth2(SPAbstractDataWithGhostNodes d, double cEq, double beta, SPAbstractDataWithGhostNodes map, const VectorTemplate *const t)
SPBCond generateBCConstantGradient(SPAbstractDataWithGhostNodes d, double v, const VectorTemplate *const t, const std::vector< SlicesNames > &sl)
Bondary condition that makes fixed gradient.
void executeAll(std::vector< T * > &v)
Definition: aslNumMethod.h:55
void addScalars(std::string name, AbstractData &data)
const VectorTemplate & d3q19()
Vector template.
void stop()
Definition: aslTimer.h:44
void start()
Definition: aslTimer.h:43
VectorOfElements generateVEConstant(T a)
Generates VectorOfElements with 1 Element acl::Constant with value a.
SPPositionFunction generatePFRotationField(const AVec< double > &axis, const AVec< double > &c)
SPDistanceFunction generateDFCylinder(double r, const AVec< double > &l, const AVec< double > &c)
generates cylinder
int main(int argc, char *argv[])
asl::SPDistanceFunction generatePlatform(asl::Block &bl)
void load(int argc, char *argv[])
const double processorLoad() const
Definition: aslTimer.h:47
asl::SPDistanceFunction generateCrystal(asl::Block &bl)
SPNumMethod generateBCNoSlipVel(SPLBGK nmU, SPAbstractDataWithGhostNodes map)
for velocity field
contains different kernels for preprocessing and posprocessing of data used by LBGK
Definition: aslLBGK.h:137
SPNumMethod generateBCVelocity(SPLBGK nm, SPPositionFunction v, SPAbstractDataWithGhostNodes map)
SPBCond generateBCNoSlip(SPLBGK nm, const std::vector< SlicesNames > &sl)