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