ASL  0.1.6
Advanced Simulation Library
locomotive.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>
36 #include <num/aslLBGK.h>
37 #include <num/aslLBGKBC.h>
38 #include <utilities/aslTimer.h>
40 
41 
42 // typedef to switch to double precision
43 //typedef double FlT;
44 typedef float FlT;
45 
46 using asl::AVec;
47 using asl::makeAVec;
48 
49 // Generate geometry of the tunnel
51 {
52 
53  // Set length of the tunnel to the length (X size) of the block
54  double l(bl.getBPosition()[0] - bl.position[0] + bl.dx);
55  // Set radius of the tunnel to the ca. half of the block's height (Z size)
56  double rTunnel((bl.getBPosition()[2] - bl.position[2]) / 2.1);
57 
58  // Center of the tunnel (described as cylinder cut by a plane)
59  asl::AVec<> center(.5 * (bl.getBPosition() + bl.position));
60  center[1] = bl.position[1] + .25 * rTunnel;
61 
62  // Center of the ground plane (that cuts the cylinder)
63  asl::AVec<> centerG(center);
64  centerG[1] = bl.position[1];
65 
66  /* DF = DistanceFunction (part of the geometrical module of ASL)
67  1. Genarate cylinder
68  2. Generate ground plane
69  3. Conjunction of the cylinder and the plane ('&' - operator)
70  4. Space inversion ('-' - operator) */
71  auto tunnel(-(generateDFCylinder(rTunnel, makeAVec(l, 0., 0.), center) &
72  generateDFPlane(makeAVec(0., -1., 0.), centerG)));
73 
74  // Normalize DistanceFunction to the range [-1; 1]
75  return normalize(tunnel, bl.dx);
76 }
77 
78 
79 int main(int argc, char* argv[])
80 {
81  /* Convenience facility to manage simulation parameters (and also
82  hardware parameters defining platform and device for computations)
83  through command line and/or parameters file.
84  See `locomotive --help` for more information */
85  asl::ApplicationParametersManager appParamsManager("locomotive",
86  "1.0");
87 
88  /* Important: declare Parameters only after declaring
89  ApplicationParametersManager instance because each Parameter adds itself
90  to it automatically!
91  0.08 - default value; will be used if nothing else is provided during
92  runtime through command line or parameters file.
93  "dx" - option key; is used to specify this parameter through command line
94  and/or parameters file, like `locomotive --dx 0.05`
95  "space step" - option description; is used in the help output:
96  `locomotive -h` and as comment on parameters file generation:
97  `locomotive -g ./defaultParameters.ini`
98  "m" - parameter units; is used to complement the option description mentioned
99  above. Might be used for automatic unit conversion in future (to this end
100  it is recommended to use the notation of the Boost::Units library). */
101  asl::Parameter<FlT> dx(0.08, "dx", "space step", "m");
102  asl::Parameter<FlT> dt(1., "dt", "time step", "s");
103  asl::Parameter<FlT> nu(.001, "nu", "kinematic viscosity", "m^2/s");
104  asl::Parameter<unsigned int> iterations(10001, "iterations", "iterations number");
105  asl::Parameter<string> input("input", "path to the geometry input file");
106 
107  /* Load previously declared Parameters from command line and/or
108  parameters file. Use default values if neither is provided. */
109  appParamsManager.load(argc, argv);
110 
111  /* Set the size of the block to 40x10x15 m (in accordance with the
112  locomotive size read later on from the input file) */
113  AVec<int> size(makeAVec(40., 10., 15.) * (1. / dx.v()));
114 
115  /* Create block and shift it in accordance with the
116  position of the locomotive in the input file */
117  asl::Block bl(size, dx.v(), makeAVec(-30., 8.58, 1.53));
118 
119  // Define dimensionless viscosity value
120  FlT nuNum(nu.v() * dt.v() / dx.v() / dx.v());
121 
122  cout << "Data initialization... " << flush;
123 
124  // Read geometry of the locomotive from the file `locomotive.stl`
125  auto locomotive(asl::readSurface(input.v(), bl));
126 
127  // Create block for further use
128  asl::Block block(locomotive->getInternalBlock());
129 
130  // Generate memory data container for the tunnel
131  auto tunnelMap(asl::generateDataContainerACL_SP<FlT>(block, 1, 1u));
132  // Place generated geometry of the tunnel into the tunnel data container
133  asl::initData(tunnelMap, generateTunnel(block));
134 
135  // Data container for air friction field
136  auto forceField(asl::generateDataContainerACL_SP<FlT>(block, 3, 1u));
137  // Initialization
138  asl::initData(forceField, makeAVec(0., 0., 0.));
139 
140  cout << "Finished" << endl;
141 
142  cout << "Numerics initialization... " << flush;
143 
144  // NOTE: the problem is considered in the reference frame related to the locomotive
145 
146  // Generate numerical method for air flow - LBGK (lattice Bhatnagar–Gross–Krook)
147  asl::SPLBGK lbgk(new asl::LBGKTurbulence(block,
148  acl::generateVEConstant(FlT(nu.v())),
149  &asl::d3q15()));
150  lbgk->init();
151  // Generate an instance for LBGK data initialization
152  asl::SPLBGKUtilities lbgkUtil(new asl::LBGKUtilities(lbgk));
153  // Initialize the LBGK internal data with the flow velocity of (0.1, 0, 0) in [lattice units]
154  lbgkUtil->initF(acl::generateVEConstant(.1, .0, .0));
155 
156  auto vfTunnel(asl::generatePFConstant(makeAVec(0.1, 0., 0.)));
157 
158  std::vector<asl::SPNumMethod> bc;
159  std::vector<asl::SPNumMethod> bcV;
160 
161  // Generate boundary conditions for the tunnel geometry. Constant velocity BC
162  bc.push_back(generateBCVelocity(lbgk, vfTunnel, tunnelMap));
163  // Generate boundary conditions for the tunnel geometry. Constant velocity BC
164  // This BC is used for visualization.
165  bcV.push_back(generateBCVelocityVel(lbgk, vfTunnel, tunnelMap));
166  bcV.push_back(generateBCNoSlipRho(lbgk, tunnelMap));
167 
168  // Generate boundary conditions for the locomotive geometry. Non-slip BC
169  bc.push_back(generateBCNoSlip(lbgk, locomotive));
170  bcV.push_back(generateBCNoSlipVel(lbgk, locomotive));
171  bcV.push_back(generateBCNoSlipRho(lbgk, locomotive));
172 
173  // Generate constant presure BC for in and out planes of the tunnel
174  bc.push_back(generateBCConstantPressureVelocity(lbgk, 1.,
175  makeAVec(0.1, 0., 0.),
176  {asl::X0, asl::XE}));
177 
178  // Initialization and building of all BC
179  initAll(bc);
180  initAll(bcV);
181 
182  // Generate a numerical method for computation of the air force field that acts on the locomotive
183  auto computeForce(generateComputeSurfaceForce(lbgk, forceField, locomotive));
184  computeForce->init();
185 
186  cout << "Finished" << endl;
187  cout << "Computing..." << endl;
188  asl::Timer timer;
189 
190  // Initialization of the output system
191  // Write the output to the directory containing the input parameters file (default "./")
192  asl::WriterVTKXML writer(appParamsManager.getDir() + "locomotive");
193  writer.addScalars("map", *locomotive);
194  writer.addScalars("tunnel", *tunnelMap);
195  writer.addScalars("rho", *lbgk->getRho());
196  writer.addVector("v", *lbgk->getVelocity());
197  writer.addVector("force", *forceField);
198 
199  // Execute all BC
200  executeAll(bc);
201  executeAll(bcV);
202  computeForce->execute();
203 
204  // First data output
205  writer.write();
206 
207  timer.start();
208  // Iteration loop
209  for (unsigned int i(1); i < iterations.v(); ++i)
210  {
211  // One iteration (timestep) of bulk numerical procedure
212  lbgk->execute();
213  // Execution of the BC procedures
214  executeAll(bc);
215  // Output and analysis scope
216  if (!(i%1000))
217  {
218  cout << i << endl;
219  // Execution of the visualization BC procedures
220  executeAll(bcV);
221  // Computation of the force field
222  computeForce->execute();
223  // Data writing
224  writer.write();
225  }
226  }
227  timer.stop();
228 
229  cout << "Finished" << endl;
230 
231  cout << "Computation statistic:" << endl;
232  cout << "time = " << timer.getTime() << "; clockTime = "
233  << timer.getClockTime() << "; load = "
234  << timer.getProcessorLoad() * 100 << "%" << endl;
235 
236  return 0;
237 }
int main(int argc, char *argv[])
Definition: locomotive.cc:79
float FlT
AVec< T > makeAVec(T a1)
const AVec normalize(const AVec< T > &a)
void initAll(std::vector< T *> &v)
Definition: aslNumMethod.h:67
SPNumMethod generateBCNoSlipRho(SPLBGK nmU, SPAbstractDataWithGhostNodes map)
AVec< T > makeAVec(T a1)
std::string getDir()
SPDataWithGhostNodesACLData readSurface(const string &fileName, double dx, acl::CommandQueue queue=acl::hardware.defaultQueue)
SPBCond generateBCConstantPressureVelocity(SPLBGK nm, double p, AVec<> v, const std::vector< SlicesNames > &sl)
double dx
Definition: aslBlocks.h:66
SPDistanceFunction generateDFPlane(const AVec< double > &n, const AVec< double > &p0)
const double getProcessorLoad() const
Definition: aslTimer.h:48
const VectorTemplate & d3q15()
Vector template.
const T & v() const
cl_int flush(void)
Definition: cl.hpp:7042
acl::VectorOfElements dx(const TemplateVE &a)
differential operator
const double getClockTime() const
Definition: aslTimer.h:47
void executeAll(std::vector< T *> &v)
Definition: aslNumMethod.h:55
asl::SPDistanceFunction generateTunnel(asl::Block &bl)
Definition: locomotive.cc:50
SPNumMethod generateBCVelocityVel(SPLBGK nm, SPPositionFunction v, SPAbstractDataWithGhostNodes map)
void initData(SPAbstractData d, double a)
const double getTime() const
Definition: aslTimer.h:46
void addScalars(std::string name, AbstractData &data)
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.
SPNumMethod generateComputeSurfaceForce(SPLBGK nm, SPDataWithGhostNodesACLData fF, SPAbstractDataWithGhostNodes map)
float FlT
Definition: locomotive.cc:44
const V getBPosition() const
Definition: aslBlocks.h:214
SPDistanceFunction generateDFCylinder(double r, const AVec< double > &l, const AVec< double > &c)
generates cylinder
void load(int argc, char *argv[])
SPNumMethod generateBCNoSlipVel(SPLBGK nmU, SPAbstractDataWithGhostNodes map)
for velocity field
SPPositionFunction generatePFConstant(const AVec< double > &a)
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)