ASL  0.1.6
Advanced Simulation Library
multicomponent_flow.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 <aslDataInc.h>
32 #include <acl/aclGenerators.h>
33 #include <writers/aslVTKFormatWriters.h>
34 #include <num/aslLBGK.h>
35 #include <num/aslLBGKBC.h>
36 #include <utilities/aslTimer.h>
37 #include <num/aslFDAdvectionDiffusion.h>
38 #include <num/aslBasicBC.h>
39 
40 // typedef to switch to double precision
41 //typedef double FlT;
42 typedef float FlT;
43 
44 using asl::AVec;
45 using asl::makeAVec;
46 
47 class Parameters
48 {
49  private:
50  void init();
51 
52  public:
54 
58 
61 
64 
69 
73 
74 
75  void load(int argc, char * argv[]);
76  Parameters();
77  void updateNumValues();
78 };
79 
80 
82  appParamsManager("multicomponent_flow", "0.1"),
83  size(3),
84  dx(0.0005, "dx", "space step"),
85  dt(1., "dt", "time step"),
86  tSimulation(2e-3, "simulation_time", "simulation time"),
87  tOutput(1e-4, "output_interval", "output interval"),
88  nu(4e-8/1.6, "nu", "viscosity"),
89  tubeL(0.25, "tubeL", "tube's length"),
90  tubeD(0.05, "tubeD", "tube's diameter"),
91  pumpL(0.025, "pumpL", "pump's length"),
92  pumpD(0.03, "pumpD", "pump's diameter"),
93  component1InVel(0.16, "component1_in_velocity", "flow velocity in the component1 input"),
94  component2InVel(0.08, "component2_in_velocity", "flow velocity in the component2 input"),
95  component3InVel(0.1, "component3_in_velocity", "flow velocity in the component3 input")
96 {
97 }
98 
99 
100 void Parameters::load(int argc, char * argv[])
101 {
102  appParamsManager.load(argc, argv);
103  init();
104 }
105 
106 
108 {
109  nuNum = nu.v() * dt.v() / dx.v() / dx.v();
110  size[0] = tubeD.v() / dx.v() + 1;
111  size[1] = (tubeD.v() + 2 * pumpL.v()) / dx.v() + 1;
112  size[2] = tubeL.v() / dx.v() + 1;
113 }
114 
115 
116 void Parameters::init()
117 {
118  if (tubeD.v() < pumpD.v())
119  asl::errorMessage("Tube's diameter is smaller than pump's diameter");
120 
121  updateNumValues();
122 }
123 
124 // Generate geometry of the mixer (cross-coupled pipes)
126 {
127  asl::SPDistanceFunction mixerGeometry;
128  asl::AVec<double> orientation(asl::makeAVec(0., 0., 1.));
129  asl::AVec<double> center(asl::AVec<double>(params.size) * .5 * params.dx.v());
130 
131  mixerGeometry = generateDFCylinderInf(params.tubeD.v() / 2., orientation,
132  center);
133 
134  orientation[1] = 1.0;
135  orientation[2] = 0.0;
136  center[2] = params.pumpD.v() * 1.5;
137  mixerGeometry = mixerGeometry | generateDFCylinderInf(params.pumpD.v() / 2.,
138  orientation, center);
139 
140  return asl::normalize(-(mixerGeometry) | asl::generateDFInBlock(block, 0),
141  params.dx.v());
142 }
143 
144 int main(int argc, char *argv[])
145 {
146  Parameters params;
147  params.load(argc, argv);
148 
149  cout << "Data initialization..." << endl;
150 
151  asl::Block block(params.size, params.dx.v());
152 
153  auto mcfMapMem(asl::generateDataContainerACL_SP<FlT>(block, 1, 1u));
154  asl::initData(mcfMapMem, generateMixer(block, params));
155 
156  auto component1Frac(asl::generateDataContainerACL_SP<FlT>(block, 1, 1u));
157  asl::initData(component1Frac, 0);
158  auto component3Frac(asl::generateDataContainerACL_SP<FlT>(block, 1, 1u));
159  asl::initData(component3Frac, 0);
160 
161 
162  cout << "Finished" << endl;
163 
164  cout << "Numerics initialization..." << endl;
165 
166  auto templ(&asl::d3q15());
167 
168  asl::SPLBGK lbgk(new asl::LBGK(block,
169  acl::generateVEConstant(FlT(params.nuNum.v())),
170  templ));
171 
172  lbgk->init();
173  asl::SPLBGKUtilities lbgkUtil(new asl::LBGKUtilities(lbgk));
174  lbgkUtil->initF(acl::generateVEConstant(.0, .0, .0));
175 
176  auto flowVel(lbgk->getVelocity());
177  auto nmcomponent1(asl::generateFDAdvectionDiffusion(component1Frac, 0.01,
178  flowVel, templ, true));
179  nmcomponent1->init();
180  auto nmcomponent3(asl::generateFDAdvectionDiffusion(component3Frac, 0.01,
181  flowVel, templ));
182  nmcomponent3->init();
183 
184  std::vector<asl::SPNumMethod> bc;
185  std::vector<asl::SPNumMethod> bcV;
186  std::vector<asl::SPNumMethod> bcDif;
187 
188  bc.push_back(generateBCNoSlip(lbgk, mcfMapMem));
189  bc.push_back(generateBCConstantPressure(lbgk, 1., {asl::ZE}));
190  bc.push_back(generateBCConstantPressureVelocity(lbgk, 1.,
191  makeAVec(0., 0., params.component2InVel.v()),
192  {asl::Z0}));
193  bc.push_back(generateBCConstantPressureVelocity(lbgk, 1.,
194  makeAVec(0., -params.component1InVel.v(), 0.),
195  {asl::YE}));
196  bc.push_back(generateBCConstantPressureVelocity(lbgk, 1.,
197  makeAVec(0., params.component3InVel.v(), 0.),
198  {asl::Y0}));
199 
200  bcDif.push_back(generateBCNoSlipVel(lbgk, mcfMapMem));
201  bc.push_back(generateBCConstantGradient(component1Frac, 0., mcfMapMem, templ));
202  bc.push_back(generateBCConstantGradient(component3Frac, 0., mcfMapMem, templ));
203  bc.push_back(generateBCConstantValue(component1Frac, 1., {asl::YE}));
204  bc.push_back(generateBCConstantValue(component3Frac, 0., {asl::YE, asl::Z0, asl::ZE}));
205  bc.push_back(generateBCConstantValue(component1Frac, 0., {asl::Y0, asl::Z0, asl::ZE}));
206  bc.push_back(generateBCConstantValue(component3Frac, 1., {asl::Y0}));
207 // bc.push_back(generateBCConstantGradient(component1Frac, 0.,templ, {asl::ZE}));
208 // bc.push_back(generateBCConstantGradient(component3Frac, 0.,templ, {asl::ZE}));
209 
210  initAll(bc);
211  initAll(bcDif);
212  initAll(bcV);
213 
214  cout << "Finished" << endl;
215  cout << "Computing..." << endl;
216  asl::Timer timer;
217 
218  asl::WriterVTKXML writer("multicomponent_flow");
219  writer.addScalars("map", *mcfMapMem);
220  writer.addScalars("component1", *component1Frac);
221  writer.addScalars("component3", *component3Frac);
222  writer.addScalars("rho", *lbgk->getRho());
223  writer.addVector("v", *flowVel);
224 
225  executeAll(bc);
226  executeAll(bcDif);
227  executeAll(bcV);
228 
229  writer.write();
230 
231  timer.start();
232  for (unsigned int i(1); i < 10001; ++i)
233  {
234  lbgk->execute();
235  executeAll(bcDif);
236  nmcomponent1->execute();
237  nmcomponent3->execute();
238  executeAll(bc);
239 
240  if (!(i%100))
241  {
242  timer.stop();
243  cout << i << "/10000; time left (expected): " << timer.getLeftTime(double(i)/10000.) << endl;
244  executeAll(bcV);
245  writer.write();
246  timer.resume();
247  }
248  }
249  timer.stop();
250 
251  cout << "Finished" << endl;
252 
253  cout << "Computation statistic:" << endl;
254  cout << "time = " << timer.getTime() << "; clockTime = "
255  << timer.getClockTime() << "; load = "
256  << timer.getProcessorLoad() * 100 << "%" << endl;
257 
258  return 0;
259 }
asl::Parameter< double > dt
const double getTime() const
Definition: aslTimer.h:46
asl::Block::DV size
asl::Parameter< double > component1InVel
const double getProcessorLoad() const
Definition: aslTimer.h:48
asl::ApplicationParametersManager appParamsManager
asl::Parameter< double > nu
Numerical method for fluid flow.
Definition: aslLBGK.h:77
void resume()
Definition: aslTimer.h:44
const double getClockTime() const
Definition: aslTimer.h:47
SPDistanceFunction normalize(SPDistanceFunction a, double dx)
void errorMessage(cl_int status, const char *errorMessage)
Prints errorMessage and exits depending on the status.
void addVector(std::string name, AbstractData &data)
const T & v() const
Definition: aslUValue.h:43
int main(int argc, char *argv[])
asl::Parameter< double > tubeL
asl::Parameter< double > pumpD
void initAll(std::vector< T * > &v)
Definition: aslNumMethod.h:67
SPBCond generateBCConstantPressureVelocity(SPLBGK nm, double p, AVec<> v, const std::vector< SlicesNames > &sl)
SPDistanceFunction generateDFCylinderInf(double r, const AVec< double > &l, const AVec< double > &c)
generates infinite cylinder
asl::UValue< double > dt
const VectorTemplate & d3q15()
Vector template.
asl::SPDistanceFunction generateMixer(asl::Block &block, Parameters &params)
const T & v() const
SPFDAdvectionDiffusion generateFDAdvectionDiffusion(SPDataWithGhostNodesACLData c, double diffustionCoeff, SPAbstractDataWithGhostNodes v, const VectorTemplate *vt, bool compressibilityCorrection=false)
asl::Parameter< double > tSimulation
float FlT
AVec< T > makeAVec(T a1)
asl::Parameter< double > tOutput
void load(int argc, char *argv[])
asl::Parameter< double > pumpL
void initData(SPAbstractData d, double a)
asl::UValue< double > nuNum
const double getLeftTime(double fractTime)
Definition: aslTimer.h:51
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)
asl::Parameter< double > flowVel
SPBCond generateBCConstantPressure(SPLBGK nm, double p, const std::vector< SlicesNames > &sl)
void stop()
Definition: aslTimer.h:45
asl::Parameter< double > dx
void start()
Definition: aslTimer.h:43
VectorOfElements generateVEConstant(T a)
Generates VectorOfElements with 1 Element acl::Constant with value a.
asl::Parameter< double > tubeD
SPBCond generateBCConstantValue(SPAbstractDataWithGhostNodes d, double v, const std::vector< SlicesNames > &sl)
Bondary condition that puts fixed value in each point.
asl::Parameter< double > component3InVel
void updateNumValues()
void load(int argc, char *argv[])
SPDistanceFunction generateDFInBlock(const Block &b, unsigned int nG)
generates map corresponding to external (ghost) part of the block
asl::Parameter< double > component2InVel
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)