ASL  0.1.7
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 
30 #include <math/aslTemplates.h>
31 #include <aslGeomInc.h>
32 #include <aslDataInc.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 // typedef to switch to double precision
42 //typedef double FlT;
43 typedef float FlT;
44 
45 using asl::AVec;
46 using asl::makeAVec;
47 
48 class Parameters
49 {
50  private:
51  void init();
52 
53  public:
55 
59 
62 
65 
70 
74 
75 
76  void load(int argc, char * argv[]);
77  Parameters();
78  void updateNumValues();
79 };
80 
81 
83  appParamsManager("multicomponent_flow", "0.1"),
84  size(3),
85  dx(0.0005, "dx", "space step"),
86  dt(1., "dt", "time step"),
87  tSimulation(2e-3, "simulation_time", "simulation time"),
88  tOutput(1e-4, "output_interval", "output interval"),
89  nu(4e-8/1.6, "nu", "viscosity"),
90  tubeL(0.25, "tubeL", "tube's length"),
91  tubeD(0.05, "tubeD", "tube's diameter"),
92  pumpL(0.025, "pumpL", "pump's length"),
93  pumpD(0.03, "pumpD", "pump's diameter"),
94  component1InVel(0.16, "component1_in_velocity", "flow velocity in the component1 input"),
95  component2InVel(0.08, "component2_in_velocity", "flow velocity in the component2 input"),
96  component3InVel(0.1, "component3_in_velocity", "flow velocity in the component3 input")
97 {
98 }
99 
100 
101 void Parameters::load(int argc, char * argv[])
102 {
103  appParamsManager.load(argc, argv);
104  init();
105 }
106 
107 
109 {
110  nuNum = nu.v() * dt.v() / dx.v() / dx.v();
111  size[0] = tubeD.v() / dx.v() + 1;
112  size[1] = (tubeD.v() + 2 * pumpL.v()) / dx.v() + 1;
113  size[2] = tubeL.v() / dx.v() + 1;
114 }
115 
116 
117 void Parameters::init()
118 {
119  if (tubeD.v() < pumpD.v())
120  asl::errorMessage("Tube's diameter is smaller than pump's diameter");
121 
122  updateNumValues();
123 }
124 
125 // Generate geometry of the mixer (cross-coupled pipes)
127 {
128  asl::SPDistanceFunction mixerGeometry;
129  asl::AVec<double> orientation(asl::makeAVec(0., 0., 1.));
130  asl::AVec<double> center(asl::AVec<double>(params.size) * .5 * params.dx.v());
131 
132  mixerGeometry = generateDFCylinderInf(params.tubeD.v() / 2., orientation,
133  center);
134 
135  orientation[1] = 1.0;
136  orientation[2] = 0.0;
137  center[2] = params.pumpD.v() * 1.5;
138  mixerGeometry = mixerGeometry | generateDFCylinderInf(params.pumpD.v() / 2.,
139  orientation, center);
140 
141  return asl::normalize(-(mixerGeometry) | asl::generateDFInBlock(block, 0),
142  params.dx.v());
143 }
144 
145 int main(int argc, char *argv[])
146 {
147  Parameters params;
148  params.load(argc, argv);
149 
150  cout << "Data initialization..." << endl;
151 
152  asl::Block block(params.size, params.dx.v());
153 
154  auto mcfMapMem(asl::generateDataContainerACL_SP<FlT>(block, 1, 1u));
155  asl::initData(mcfMapMem, generateMixer(block, params));
156 
157  auto component1Frac(asl::generateDataContainerACL_SP<FlT>(block, 1, 1u));
158  asl::initData(component1Frac, 0);
159  auto component3Frac(asl::generateDataContainerACL_SP<FlT>(block, 1, 1u));
160  asl::initData(component3Frac, 0);
161 
162 
163  cout << "Finished" << endl;
164 
165  cout << "Numerics initialization..." << endl;
166 
167  auto templ(&asl::d3q15());
168 
169  asl::SPLBGK lbgk(new asl::LBGK(block,
170  acl::generateVEConstant(FlT(params.nuNum.v())),
171  templ));
172 
173  lbgk->init();
174  asl::SPLBGKUtilities lbgkUtil(new asl::LBGKUtilities(lbgk));
175  lbgkUtil->initF(acl::generateVEConstant(.0, .0, .0));
176 
177  auto flowVel(lbgk->getVelocity());
178  auto nmcomponent1(asl::generateFDAdvectionDiffusion(component1Frac, 0.01,
179  flowVel, templ, true));
180  nmcomponent1->init();
181  auto nmcomponent3(asl::generateFDAdvectionDiffusion(component3Frac, 0.01,
182  flowVel, templ));
183  nmcomponent3->init();
184 
185  std::vector<asl::SPNumMethod> bc;
186  std::vector<asl::SPNumMethod> bcV;
187  std::vector<asl::SPNumMethod> bcDif;
188 
189  bc.push_back(generateBCNoSlip(lbgk, mcfMapMem));
190  bc.push_back(generateBCConstantPressure(lbgk, 1., {asl::ZE}));
191  bc.push_back(generateBCConstantPressureVelocity(lbgk, 1.,
192  makeAVec(0., 0., params.component2InVel.v()),
193  {asl::Z0}));
194  bc.push_back(generateBCConstantPressureVelocity(lbgk, 1.,
195  makeAVec(0., -params.component1InVel.v(), 0.),
196  {asl::YE}));
197  bc.push_back(generateBCConstantPressureVelocity(lbgk, 1.,
198  makeAVec(0., params.component3InVel.v(), 0.),
199  {asl::Y0}));
200 
201  bcDif.push_back(generateBCNoSlipVel(lbgk, mcfMapMem));
202  bc.push_back(generateBCConstantGradient(component1Frac, 0., mcfMapMem, templ));
203  bc.push_back(generateBCConstantGradient(component3Frac, 0., mcfMapMem, templ));
204  bc.push_back(generateBCConstantValue(component1Frac, 1., {asl::YE}));
205  bc.push_back(generateBCConstantValue(component3Frac, 0., {asl::YE, asl::Z0, asl::ZE}));
206  bc.push_back(generateBCConstantValue(component1Frac, 0., {asl::Y0, asl::Z0, asl::ZE}));
207  bc.push_back(generateBCConstantValue(component3Frac, 1., {asl::Y0}));
208 // bc.push_back(generateBCConstantGradient(component1Frac, 0.,templ, {asl::ZE}));
209 // bc.push_back(generateBCConstantGradient(component3Frac, 0.,templ, {asl::ZE}));
210 
211  initAll(bc);
212  initAll(bcDif);
213  initAll(bcV);
214 
215  cout << "Finished" << endl;
216  cout << "Computing..." << endl;
217  asl::Timer timer;
218 
219  asl::WriterVTKXML writer("multicomponent_flow");
220  writer.addScalars("map", *mcfMapMem);
221  writer.addScalars("component1", *component1Frac);
222  writer.addScalars("component3", *component3Frac);
223  writer.addScalars("rho", *lbgk->getRho());
224  writer.addVector("v", *flowVel);
225 
226  executeAll(bc);
227  executeAll(bcDif);
228  executeAll(bcV);
229 
230  writer.write();
231 
232  timer.start();
233  for (unsigned int i(1); i < 10001; ++i)
234  {
235  lbgk->execute();
236  executeAll(bcDif);
237  nmcomponent1->execute();
238  nmcomponent3->execute();
239  executeAll(bc);
240 
241  if (!(i%100))
242  {
243  timer.stop();
244  cout << i << "/10000; time left (estimated): " << timer.estimatedRemainder(double(i)/10000.) << endl;
245  executeAll(bcV);
246  writer.write();
247  timer.start();
248  }
249  }
250  timer.stop();
251 
252  cout << "Finished" << endl;
253 
254  cout << "Computation statistic:" << endl;
255  cout << "Real Time = " << timer.realTime() << "; Processor Time = "
256  << timer.processorTime() << "; Processor Load = "
257  << timer.processorLoad() * 100 << "%" << endl;
258 
259  return 0;
260 }
asl::Parameter< double > dt
asl::Block::DV size
const double estimatedRemainder(double completeness)
Returns estimated time till finishing current task based on its current completeness [0....
Definition: aslTimer.h:52
asl::Parameter< double > component1InVel
const double realTime() const
Definition: aslTimer.h:45
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.
void addVector(std::string name, AbstractData &data)
int main(int argc, char *argv[])
asl::Parameter< double > tubeL
asl::Parameter< double > pumpD
void initAll(std::vector< T * > &v)
Definition: aslNumMethod.h:67
AVec< T > makeAVec(T a1)
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
asl::SPDistanceFunction generateMixer(asl::Block &block, Parameters &params)
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
SPFDAdvectionDiffusion generateFDAdvectionDiffusion(SPDataWithGhostNodesACLData c, double diffustionCoeff, SPAbstractDataWithGhostNodes v, const VectorTemplate *vt, bool compressibilityCorrection=false)
asl::Parameter< double > tSimulation
acl::VectorOfElements dx(const TemplateVE &a)
differential operator
float FlT
asl::Parameter< double > tOutput
void load(int argc, char *argv[])
asl::Parameter< double > pumpL
void initData(SPAbstractData d, double a)
asl::UValue< double > nuNum
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)
SPBCond generateBCConstantPressure(SPLBGK nm, double p, const std::vector< SlicesNames > &sl)
void stop()
Definition: aslTimer.h:44
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[])
const double processorLoad() const
Definition: aslTimer.h:47
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)