ASL 0.1.7
Advanced Simulation Library
Loading...
Searching...
No Matches
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>
42#include <utilities/aslTimer.h>
43
44using asl::AVec;
45using 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
159double 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
192typedef float FlT;
193//typedef double FlT;
195
196using asl::AVec;
197using asl::makeAVec;
198
199
200int 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}
void load(int argc, char *argv[])
double dx
Definition aslBlocks.h:66
const DV & getSize() const
Definition aslBlocks.h:208
Numerical method for fluid flow.
Definition aslLBGK.h:78
contains different kernels for preprocessing and posprocessing of data used by LBGK
Definition aslLBGK.h:138
const double realTime() const
Definition aslTimer.h:45
void stop()
Definition aslTimer.h:44
const double processorTime() const
Definition aslTimer.h:46
void start()
Definition aslTimer.h:43
void reset()
Definition aslTimer.h:48
const double processorLoad() const
Definition aslTimer.h:47
Updatable value. This class stores value and its TimeStamp.
Definition aslUValue.h:35
const T & v() const
Definition aslUValue.h:43
void addVector(std::string name, AbstractData &data)
void addScalars(std::string name, AbstractData &data)
float FlT
asl::SPDistanceFunction generateCrystal(asl::Block &bl)
asl::SPDistanceFunction generateBath(asl::Block &bl)
asl::UValue< double > Param
double getWRotation(double t)
asl::SPDistanceFunction generatePlatform(asl::Block &bl)
SPDistanceFunction generateDFConvexPolygonPyramid(std::vector< AVec< double > > points, AVec< double > a)
generates pyramid with convex polygon at its base and apex a
SPDistanceFunction generateDFConvexPolygonPrism(std::vector< AVec< double > > points)
generates infinite prism with convex polygon at its base
std::shared_ptr< DistanceFunction > SPDistanceFunction
Definition aslGeomInc.h:45
SPFDAdvectionDiffusion generateFDAdvectionDiffusion(SPDataWithGhostNodesACLData c, double diffustionCoeff, SPAbstractDataWithGhostNodes v, const VectorTemplate *vt, bool compressibilityCorrection=false)
SPPositionFunction generatePFRotationField(const AVec< double > &axis, const AVec< double > &c)
const VectorTemplate & d3q19()
Vector template.
VectorOfElements generateVEConstant(T a)
Generates VectorOfElements with 1 Element acl::Constant with value a.
AVec< T > makeAVec(T a1)
std::shared_ptr< LBGKUtilities > SPLBGKUtilities
Definition aslLBGK.h:161
std::shared_ptr< LBGK > SPLBGK
Definition aslLBGK.h:133
void initData(SPAbstractData d, double a)
int main()