ASL  0.1.7
Advanced Simulation Library
pitot_tube_ice.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>
32 #include <aslGenerators.h>
33 #include <acl/aclGenerators.h>
35 #include <num/aslLBGK.h>
36 #include <num/aslLBGKBC.h>
37 #include <utilities/aslTimer.h>
39 #include <num/aslBasicBC.h>
40 
41 
42 typedef float FlT;
43 //typedef double FlT;
45 
46 using asl::AVec;
47 using asl::makeAVec;
48 
49 class Parameters
50 {
51  private:
52  void init();
53 
54  public:
56 
58 
61 
64 
67 
72 
76 
77  void load(int argc, char * argv[]);
78  Parameters();
79  void updateNumValues();
80 };
81 
82 
84  appParamsManager("pitot_tube_ice", "0.1"),
85  size(3),
86  dx(0.000125, "dx", "space step"),
87  dt(1., "dt", "time step"),
88  tSimulation(2e-3, "simulation_time", "simulation time"),
89  tOutput(1e-4, "output_interval", "output interval"),
90  nu(6.25e-10/4., "nu", "viscosity"),
91  rIn(0.0015, "r_in", "Internal radius, m"),
92  rEx(0.005, "r_ex", "External radius, m"),
93  lCyl(0.002, "l_cyl", "Length of cylindric part, m"),
94  lCone(0.02, "l_cone", "Length of conic part, m"),
95  temperature(253, "temperature", "temperature, K"),
96  humidity(.5, "humidity", "relative humidity, K"),
97  flowVel(0.08, "flow_vel", "flow velocity")
98 {
99 }
100 
101 
102 void Parameters::load(int argc, char * argv[])
103 {
104  appParamsManager.load(argc, argv);
105 
106  init();
107 }
108 
109 
111 {
112  nuNum = nu.v() * dt.v() / dx.v() / dx.v();
113  size[0] = 1.0*(lCyl.v() + lCone.v()) / dx.v() + 1;
114  size[1] = rEx.v() * 2.5 / dx.v() + 1;
115  size[2] = rEx.v() * 2.5 / dx.v() + 1;
116 }
117 
118 
119 void Parameters::init()
120 {
121  if (rEx.v() < rIn.v())
122  asl::errorMessage("External tube's diameter is smaller than internal one");
123 
124  updateNumValues();
125 }
126 
127 // generate geometry
129 {
130  asl::SPDistanceFunction tubeGeometry;
131  asl::AVec<double> orientation(asl::makeAVec(1., 0., 0.));
132  asl::AVec<double> center(asl::AVec<double>(params.size)*.5*params.dx.v());
133  auto centerCyl(center);
134  centerCyl[0] = params.lCyl.v()*.45;
135  auto centerHole(centerCyl+(params.lCone.v()*.6)*orientation);
136  auto lHole(params.lCyl.v()+params.lCone.v());
137 
138  auto apexCone(centerCyl+orientation*(params.lCyl.v()*.49+params.lCone.v()));
139 
140  tubeGeometry = ((generateDFCylinder(params.rEx.v(), orientation*params.lCyl.v(), centerCyl) |
141  generateDFCone(params.rEx.v()*.98, -orientation*params.lCone.v(), apexCone)) &
142  -generateDFCylinder(params.rIn.v(), orientation*lHole, centerHole)) &
143  generateDFPlane(orientation, apexCone-orientation*params.lCone.v()*.5);
144 
145  return asl::normalize(tubeGeometry, params.dx.v());
146 }
147 
148 int main(int argc, char *argv[])
149 {
150  Parameters params;
151  params.load(argc, argv);
152 
153  std::cout << "Data initialization...";
154 
155  asl::Block block(params.size, params.dx.v());
156 
157  auto mcfMapMem(asl::generateDataContainerACL_SP<FlT>(block, 1, 1u));
158  asl::initData(mcfMapMem, generateGeometry(block, params));
159 
160 // auto waterFrac(asl::generateDataContainerACL_SP<FlT>(block, 1, 1u));
161 // asl::initData(waterFrac, 0);
162 
163 
164  std::cout << "Finished" << endl;
165 
166  std::cout << "Flow: Numerics initialization...";
167 
168  auto templ(&asl::d3q15());
169 
170  asl::SPLBGK lbgk(new asl::LBGK(block,
171  acl::generateVEConstant(FlT(params.nuNum.v())),
172  templ));
173 
174  lbgk->init();
175  asl::SPLBGKUtilities lbgkUtil(new asl::LBGKUtilities(lbgk));
176  lbgkUtil->initF(acl::generateVEConstant(-0.9*params.flowVel.v(), params.flowVel.v()*.4, .0));
177 
178  auto flowVel(lbgk->getVelocity());
179 // auto nmWater(asl::generateFDAdvectionDiffusion(waterFrac, 0.01, flowVel, templ, false));
180 // nmWater->init();
181 
182  std::vector<asl::SPNumMethod> bc;
183  std::vector<asl::SPNumMethod> bcV;
184  std::vector<asl::SPNumMethod> bcDif;
185 
186  bc.push_back(generateBCNoSlip(lbgk, mcfMapMem));
187  bc.push_back(generateBCConstantPressure(lbgk,1.,{asl::ZE}));
188  bc.push_back(generateBCConstantPressureVelocity(lbgk, 1.,
189  makeAVec(-params.flowVel.v()*.9,params.flowVel.v()*.3,0.),
191 
192  bcDif.push_back(generateBCNoSlipVel(lbgk, mcfMapMem));
193  bcV.push_back(generateBCNoSlipRho(lbgk, mcfMapMem));
194 // bc.push_back(generateBCConstantGradient(waterFrac, 0., mcfMapMem, templ));
195 // bc.push_back(generateBCConstantValue(waterFrac, 1., {asl::YE}));
196 // bc.push_back(generateBCConstantValue(waterFrac, 0., {asl::Y0, asl::Z0, asl::ZE}));
197 
198  initAll(bc);
199  initAll(bcDif);
200  initAll(bcV);
201 
202  std::cout << "Finished" << endl;
203  std::cout << "Computing..." << endl;
204  asl::Timer timer;
205 
206  asl::WriterVTKXML writer(params.appParamsManager.getDir() + "pitot_tube");
207  writer.addScalars("map", *mcfMapMem);
208 // writer.addScalars("water", *waterFrac);
209  writer.addScalars("rho", *lbgk->getRho());
210  writer.addVector("v", *flowVel);
211 
212  executeAll(bc);
213  executeAll(bcDif);
214  executeAll(bcV);
215 
216  writer.write();
217 
218  timer.start();
219  for(unsigned int i(1); i < 8001; ++i)
220  {
221  lbgk->execute();
222  executeAll(bcDif);
223 // nmWater->execute();
224  executeAll(bc);
225 
226  if(!(i%800))
227  {
228  timer.stop();
229  cout << i << "/8000; time left (estimated): " << timer.estimatedRemainder(double(i)/8000.) << endl;
230  executeAll(bcV);
231  writer.write();
232  timer.start();
233  }
234  }
235  timer.stop();
236 
237  cout << "Finished" << endl;
238 
239  cout << "Computation statistic:" << endl;
240  cout << "Real Time = " << timer.realTime() << "; Processor Time = "
241  << timer.processorTime() << "; Processor Load = "
242  << timer.processorLoad() * 100 << "%" << endl;
243 
244  return 0;
245 }
asl::Block::DV size
SPDistanceFunction generateDFCone(double r, const AVec< double > &l, const AVec< double > &a)
generates cone
const double estimatedRemainder(double completeness)
Returns estimated time till finishing current task based on its current completeness [0....
Definition: aslTimer.h:52
const double realTime() const
Definition: aslTimer.h:45
asl::Parameter< double > humidity
std::shared_ptr< DistanceFunction > SPDistanceFunction
Definition: aslGeomInc.h:44
const double processorTime() const
Definition: aslTimer.h:46
AVec< T > makeAVec(T a1)
asl::ApplicationParametersManager appParamsManager
asl::Parameter< double > nu
Numerical method for fluid flow.
Definition: aslLBGK.h:77
SPDistanceFunction normalize(SPDistanceFunction a, double dx)
const T & v() const
Definition: aslUValue.h:43
void errorMessage(cl_int status, const char *errorMessage)
Prints errorMessage and exits depending on the status.
asl::Parameter< double > rEx
SPNumMethod generateBCNoSlipRho(SPLBGK nmU, SPAbstractDataWithGhostNodes map)
void initAll(std::vector< T * > &v)
Definition: aslNumMethod.h:67
AVec< T > makeAVec(T a1)
std::string getDir()
SPBCond generateBCConstantPressureVelocity(SPLBGK nm, double p, AVec<> v, const std::vector< SlicesNames > &sl)
asl::Parameter< double > lCone
asl::UValue< double > dt
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
const VectorTemplate & d3q15()
Vector template.
const T & v() const
asl::Parameter< double > tSimulation
acl::VectorOfElements dx(const TemplateVE &a)
differential operator
int main(int argc, char *argv[])
asl::UValue< double > Param
float FlT
asl::Parameter< double > tOutput
void load(int argc, char *argv[])
void initData(SPAbstractData d, double a)
asl::Parameter< double > lCyl
asl::SPDistanceFunction generateGeometry(asl::Block &block, Parameters &params)
asl::UValue< double > nuNum
void executeAll(std::vector< T * > &v)
Definition: aslNumMethod.h:55
void addScalars(std::string name, AbstractData &data)
asl::Parameter< double > flowVel
asl::Parameter< double > temperature
SPBCond generateBCConstantPressure(SPLBGK nm, double p, const std::vector< SlicesNames > &sl)
void stop()
Definition: aslTimer.h:44
asl::Parameter< double > dx
asl::Parameter< double > rIn
void start()
Definition: aslTimer.h:43
VectorOfElements generateVEConstant(T a)
Generates VectorOfElements with 1 Element acl::Constant with value a.
void updateNumValues()
SPDistanceFunction generateDFCylinder(double r, const AVec< double > &l, const AVec< double > &c)
generates cylinder
void load(int argc, char *argv[])
const double processorLoad() const
Definition: aslTimer.h:47
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
SPBCond generateBCNoSlip(SPLBGK nm, const std::vector< SlicesNames > &sl)