ASL  0.1.6
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 
28 #include <utilities/aslParametersManager.h>
29 #include <math/aslTemplates.h>
30 #include <aslGeomInc.h>
31 #include <math/aslPositionFunction.h>
32 #include <aslDataInc.h>
33 #include <acl/aclGenerators.h>
34 #include <acl/aclMath/aclVectorOfElements.h>
35 #include <writers/aslVTKFormatWriters.h>
36 #include <num/aslLBGK.h>
37 #include <num/aslLBGKBC.h>
38 #include <num/aslBasicBC.h>
39 #include <num/aslCrystalGrowthBC.h>
40 #include <num/aslFDAdvectionDiffusion.h>
41 #include <utilities/aslTimer.h>
42 
43 using asl::AVec;
44 using asl::makeAVec;
45 
47 {
48 
49 // double hBath(2.);
50  double rBath(1.);
51 
52  double dx(bl.dx);
53 
54  asl::AVec<int> size(bl.getSize());
55 
56  asl::AVec<>center(.5*dx*AVec<>(size));
57 
58 
59 
60  auto bath(-(generateDFCylinderInf(rBath, makeAVec(0.,0.,1.), dx*AVec<>(size)*.5) &
61  generateDFPlane(makeAVec(0.,0.,1.), center*1.99) &
62  generateDFPlane(makeAVec(0.,0.,-1.), center*0.)));
63 
64  return normalize(bath, dx);
65 }
66 
68 {
69  double rDisk(.9);
70  double hDisk(0.1);
71 
72  double rAxis(0.05);
73  double hAxis(.5);
74 
75  double wPillar(.2);
76  double dPillar(.1);
77 
78  double dx(bl.dx);
79  asl::AVec<int> size(bl.getSize());
80  asl::AVec<>center(.5*dx*AVec<>(size));
81 
82  vector<asl::AVec<>> pillar1{makeAVec(wPillar*.5, dPillar*.5,0.),
83  makeAVec(-wPillar*.5, dPillar*.5,0.),
84  makeAVec(-wPillar*.5, -dPillar*.5,0.),
85  makeAVec(wPillar*.5, -dPillar*.5,0.)};
86 
87  vector<asl::AVec<>> pillar2{makeAVec(dPillar*.5, wPillar*.5,0.),
88  makeAVec(-dPillar*.5, wPillar*.5,0.),
89  makeAVec(-dPillar*.5, -wPillar*.5,0.),
90  makeAVec(dPillar*.5, -wPillar*.5,0.)};
91 
92  vector<asl::AVec<>> pillarC{makeAVec(center[0]+rDisk-dPillar*.5, center[1], 0.),
93  makeAVec(center[0]-rDisk+dPillar*.5, center[1], 0.),
94  makeAVec(center[0], center[1]+rDisk-dPillar*.5,0.),
95  makeAVec(center[0], center[1]-rDisk+dPillar*.5,0.)};
96  vector<vector<asl::AVec<>>> pillarsPoints(4);
97  for(unsigned int i(0); i<4; ++i)
98  pillarsPoints[i].resize(4);
99 
100  for(unsigned int i(0); i<4; ++i)
101  {
102  pillarsPoints[0][i] = pillar2[i] + pillarC[0];
103  pillarsPoints[1][i] = pillar2[i] + pillarC[1];
104  pillarsPoints[2][i] = pillar1[i] + pillarC[2];
105  pillarsPoints[3][i] = pillar1[i] + pillarC[3];
106  }
107 
108 
109  auto diskBottom(generateDFCylinder(rDisk,
110  makeAVec(0., 0., hDisk),
111  makeAVec(center[0], center[1], .5*hDisk)));
112  auto diskTop(generateDFCylinder(rDisk,
113  makeAVec(0., 0., hDisk),
114  makeAVec(center[0], center[1], -.5*hDisk - hAxis + dx*size[2])));
115  auto axis(generateDFCylinder(rAxis,
116  makeAVec(0., 0., hAxis+hDisk*.5),
117  makeAVec(center[0], center[1], - .5*hAxis - hDisk*.25 + dx*size[2])));
118  auto dfPillar1(generateDFConvexPolygonPrism(pillarsPoints[0]));
119  auto dfPillar2(generateDFConvexPolygonPrism(pillarsPoints[1]));
120  auto dfPillar3(generateDFConvexPolygonPrism(pillarsPoints[2]));
121  auto dfPillar4(generateDFConvexPolygonPrism(pillarsPoints[3]));
122  auto dfPillars((dfPillar1 | dfPillar2 | dfPillar3 | dfPillar4) &
123  generateDFPlane(makeAVec(0.,0.,-1.), makeAVec(0.,0.,.5*hDisk)) &
124  generateDFPlane(makeAVec(0.,0.,1.), makeAVec(0.,0.,-.5*hDisk - hAxis + dx*size[2])));
125 
126  return normalize(diskBottom | diskTop | axis | dfPillars, dx);
127 }
128 
130 {
131 
132  double aCrystal(.5);
133  double hCrystalBase(.5);
134  double hCrystalPyramid(.5);
135 
136  double hDisk(0.1);
137 
138  double dx(bl.dx);
139  asl::AVec<int> size(bl.getSize());
140  asl::AVec<>center(.5*dx*AVec<>(size));
141 
142  auto crystalB(asl::generateDFConvexPolygonPrism({center+makeAVec( aCrystal, aCrystal,0.),
143  center+makeAVec(-aCrystal, aCrystal,0.),
144  center+makeAVec(-aCrystal, -aCrystal,0.),
145  center+makeAVec( aCrystal, -aCrystal,0.)}) &
146  generateDFPlane(makeAVec(0.,0.,-1.), makeAVec(0.,0., hDisk-.001)) &
147  generateDFPlane(makeAVec(0.,0., 1.), makeAVec(0.,0., hDisk + hCrystalBase)));
148  auto cCrPyrBase(makeAVec(center[0],center[1],hDisk+hCrystalBase-.01));
149  auto crystalT(asl::generateDFConvexPolygonPyramid({cCrPyrBase+makeAVec( aCrystal, aCrystal,0.),
150  cCrPyrBase+makeAVec(-aCrystal, aCrystal,0.),
151  cCrPyrBase+makeAVec(-aCrystal, -aCrystal,0.),
152  cCrPyrBase+makeAVec( aCrystal, -aCrystal,0.)},
153  cCrPyrBase+makeAVec(0.,0.,hCrystalPyramid)));
154  return normalize(crystalB | crystalT, dx);
155 // return crystalB | crystalT;
156 }
157 
158 double getWRotation(double t)
159 {
160  double tPeriod(128);
161  double wMax(6.*3.14*2./60.);
162  double tPlato(tPeriod * .25);
163  double tAcceleration(tPeriod * .1);
164  double tStop(tPeriod * .05);
165 
166  double intPart;
167  double tRel(modf(t/tPeriod, &intPart));
168  double x(0);
169  if(tRel<=tAcceleration)
170  x = tRel / tAcceleration;
171  if(tRel>tAcceleration && tRel<=tAcceleration+tPlato)
172  x = 1.;
173  if(tRel>tAcceleration+tPlato && tRel<=2.*tAcceleration+tPlato)
174  x = (2.*tAcceleration + tPlato - tRel) / tAcceleration;
175  if(tRel>2.*tAcceleration+tPlato && tRel<=2.*tAcceleration+tPlato+tStop)
176  x = 0;
177  if(tRel>2.*tAcceleration+tPlato+tStop && tRel<=3.*tAcceleration+tPlato+tStop)
178  x = -(tRel-2.*tAcceleration-tPlato-tStop) / tAcceleration;
179  if(tRel>3.*tAcceleration+tPlato+tStop && tRel<=3.*tAcceleration+2.*tPlato+tStop)
180  x = -1.;
181  if(tRel>3.*tAcceleration+2.*tPlato+tStop && tRel<=4.*tAcceleration+2.*tPlato+tStop)
182  x = -(4.*tAcceleration+2.*tPlato+tStop-tRel)/tAcceleration;
183  if(tRel>4.*tAcceleration+2.*tPlato+tStop)
184  x = 0;
185  return wMax*x;
186 
187 // flux = -9.32e-5*(1.170 - c); c_0=0.326 ceq=0.267
188 }
189 
190 
191 typedef float FlT;
192 //typedef double FlT;
194 
195 using asl::AVec;
196 using asl::makeAVec;
197 
198 
199 int main(int argc, char* argv[])
200 {
201  // Optionally add appParamsManager to be able to manipulate at least
202  // hardware parameters(platform/device) through command line/parameters file
203  asl::ApplicationParametersManager appParamsManager("flowKDPGrowth",
204  "1.0");
205  appParamsManager.load(argc, argv);
206 
207  Param dx(.02);
208  Param dt(0.8e-2);
209  Param nu(1e-2);
210  Param difC(1e-2/300.);
211 // Param w(48.*3.14*2./60.);
212  // Angular velocity
213  Param w(6.*3.14*2./60.);
214 
215  Param nuNum(nu.v()*dt.v()/dx.v()/dx.v());
216  Param difCNum(difC.v()*dt.v()/dx.v()/dx.v());
217 
218  // Angular velocity in one iteration
219  Param wNum(w.v()*dt.v());
220 
221  Param c0(0.326);
222 
223  AVec<int> size(asl::makeAVec(105.,105.,100.));
224 
225  AVec<> gSize(dx.v()*AVec<>(size));
226 
227  std::cout << "Data initialization...";
228 
229  auto templ(&asl::d3q19());
230  asl::Block block(size,dx.v());
231 
232  auto bathMap(asl::generateDataContainerACL_SP<FlT>(block, 1, 1u));
233  asl::initData(bathMap, generateBath(block));
234  auto platformCrysMap(asl::generateDataContainerACL_SP<FlT>(block, 1, 1u));
235  asl::initData(platformCrysMap, generatePlatform(block) | generateCrystal(block));
236  auto bathPlatformMap(asl::generateDataContainerACL_SP<FlT>(block, 1, 1u));
237  asl::initData(bathPlatformMap, generateBath(block) | generatePlatform(block));
238  auto bathPlatformCrystalMap(asl::generateDataContainerACL_SP<FlT>(block, 1, 1u));
239  asl::initData(bathPlatformCrystalMap, generateBath(block) | generatePlatform(block) | generateCrystal(block));
240  auto crystalMap(asl::generateDataContainerACL_SP<FlT>(block, 1, 1u));
241  asl::initData(crystalMap, generateCrystal(block));
242 
243  auto cField(asl::generateDataContainerACL_SP<FlT>(block, 1, 1u));
244  asl::initData(cField, c0.v());
245 
246  std::cout << "Finished" << endl;
247 
248  std::cout << "Numerics initialization...";
249 
250  asl::SPLBGK lbgk(new asl::LBGK(block,
251  acl::generateVEConstant(FlT(nuNum.v())),
252  templ));
253  // Set angular velocity in lbgk
254  lbgk->setOmega(acl::generateVEConstant(makeAVec(0.,0.,wNum.v())));
255  lbgk->init();
256  asl::SPLBGKUtilities lbgkUtil(new asl::LBGKUtilities(lbgk));
257  lbgkUtil->initF(acl::generateVEConstant(.0,.0,.0));
258 
259  auto nmDif(asl::generateFDAdvectionDiffusion(cField,
260  difCNum.v(),
261  lbgk->getVelocity(),
262  templ,
263  true));
264 
265  nmDif->init();
266  std::vector<asl::SPNumMethod> bc;
267  std::vector<asl::SPNumMethod> bcV;
268  std::vector<asl::SPNumMethod> bcDif;
269 
270  // Position Function Angular Velocity Field
271  auto vfBath(asl::generatePFRotationField(makeAVec(0.,0., wNum.v()/dx.v()), .5*gSize));
272  // Boundary condition
273  bc.push_back(generateBCVelocity(lbgk, vfBath, bathMap));
274  // Boundary condition for visualization
275  bcV.push_back(generateBCVelocityVel(lbgk, vfBath, bathMap));
276  bc.push_back(generateBCNoSlip(lbgk, platformCrysMap));
277  bcV.push_back(generateBCNoSlipVel(lbgk, platformCrysMap));
278  bcDif.push_back(generateBCConstantGradient(cField,
279  0.,
280  bathPlatformMap,
281  bathPlatformCrystalMap,
282  templ));
283  bcDif.push_back(generateBCLinearGrowth2(cField, 1.17,
284  -9.32e-6/difC.v()*dx.v(),
285  crystalMap,
286  bathPlatformCrystalMap,
287  templ));
288 // bcDif.push_back(generateBCConstantGradient2(cField, .1, crystalMap, bathPlatformCrystalMap, templ));
289  initAll(bc);
290  initAll(bcV);
291  initAll(bcDif);
292 
293  std::cout << "Finished" << endl;
294  std::cout << "Computing...";
295  asl::Timer timer;
296  asl::Timer timerBC;
297 
298  asl::WriterVTKXML writer("flowKDPGrowthRes");
299  writer.addScalars("mapBath", *bathMap);
300  writer.addScalars("mapPlatformCrys", *platformCrysMap);
301  writer.addScalars("mapBathPlatformCrystal", *bathPlatformCrystalMap);
302  writer.addScalars("mapCrys", *crystalMap);
303  writer.addScalars("rho", *lbgk->getRho());
304  writer.addScalars("c", *cField);
305  writer.addVector("v", *lbgk->getVelocity());
306 
307  executeAll(bcV);
308  executeAll(bc);
309  executeAll(bcDif);
310 
311  writer.write();
312 
313  timer.start();
314  timerBC.reset();
315  for (unsigned int i(0); i <= 8001 ; ++i)
316  {
317  lbgk->execute();
318  timerBC.resume();
319  executeAll(bcV);
320  executeAll(bc);
321  timerBC.stop();
322  nmDif->execute();
323  timerBC.resume();
324  executeAll(bcDif);
325  timerBC.stop();
326 
327  if (!(i%2000))
328  {
329  cout << i << endl;
330  writer.write();
331  }
332  }
333  timer.stop();
334 
335  std::cout << "Finished" << endl;
336 
337  cout << "time=" << timer.getTime() << "; clockTime="
338  << timer.getClockTime() << "; load="
339  << timer.getProcessorLoad() * 100 << "%; timeBC = "
340  << timerBC.getTime() << endl;
341 
342  std::cout << "Output...";
343  std::cout << "Finished" << endl;
344  std::cout << "Ok" << endl;
345 
346  return 0;
347 }
const double getTime() const
Definition: aslTimer.h:46
asl::UValue< double > Param
const double getProcessorLoad() const
Definition: aslTimer.h:48
Numerical method for fluid flow.
Definition: aslLBGK.h:77
void resume()
Definition: aslTimer.h:44
const double getClockTime() const
Definition: aslTimer.h:47
const AVec normalize(const AVec< T > &a)
void addVector(std::string name, AbstractData &data)
const T & v() const
Definition: aslUValue.h:43
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)
float FlT
SPFDAdvectionDiffusion generateFDAdvectionDiffusion(SPDataWithGhostNodesACLData c, double diffustionCoeff, SPAbstractDataWithGhostNodes v, const VectorTemplate *vt, bool compressibilityCorrection=false)
acl::VectorOfElements dx(const TemplateVE &a)
differential operator
AVec< T > makeAVec(T a1)
void reset()
Definition: aslTimer.h:49
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 DV & getSize() const
Definition: aslBlocks.h:208
void stop()
Definition: aslTimer.h:45
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 VectorTemplate & d3q19()
Vector template.
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)