ASL 0.1.7
Advanced Simulation Library
Loading...
Searching...
No Matches
acousticWaves.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 <aslDataInc.h>
29#include <acl/aclGenerators.h>
31#include <num/aslFDElasticity.h>
33#include <num/aslBasicBC.h>
34#include <utilities/aslTimer.h>
36#include <math/aslTemplates.h>
37#include <aslGeomInc.h>
38
39
40typedef float FlT;
41//typedef asl::UValue<FlT> Param;
42
44{
45 private:
46 void init();
47
48 public:
50
52
64
67
71
72 void load(int argc, char * argv[]);
73 Parameters();
74 void updateNumValues();
75};
76
77
79 appParamsManager("acousticWaves", "0.1"),
80 size(3),
81 dx(1e-3,"dx", "dx"),
82 bulkModulus(160e9,"bulk_modulus", "bulk modulus"),
83 shearModulus(79e9,"shear_modulus", "shear modulus"),
84 rho(7800,"rho", "density"),
85 tubeL(.2,"tube_length", "pipe length" "m"),
86 tubeDEx(0.021, "tube_diameter_external", "external pipe diameter" "m"),
87// tubeDIn(0.0157,"tube_diameter_internal", "internal pipe diameter" "m"),
88 tubeDIn(0.0107,"tube_diameter_internal", "internal pipe diameter" "m"),
89 hole1Pos(0.1,"hole_1_position", "position of first hole" "m"),
90 hole2Pos(0.15,"hole_2_position", "position of second hole" "m"),
91 hole1D(15e-3,"hole_1_diameter", "diameter of first hole" "m"),
92 hole2D(15e-3,"hole_2_diameter", "diameter of second hole" "m"),
93 tSimulation(8e-5, "simulation_time", "simulation time"),
94 tOutput(1e-6, "output_interval", "output interval")
95{
96}
97
98
99void Parameters::load(int argc, char * argv[])
100{
101 appParamsManager.load(argc, argv);
102
103 init();
104}
105
106
108{
109 double vs(sqrt((bulkModulus.v()+2.*shearModulus.v())/rho.v()));
110 dt=dx.v()/vs*.1;
111 cout << vs << "; " << dx.v() << "; " << dt.v() << endl;
112 bulkMNum = bulkModulus.v()/rho.v()/dx.v()/dx.v();
113 shearMNum = shearModulus.v()/rho.v()/dx.v()/dx.v();
114 size = asl::makeAVec(tubeL.v() / dx.v() + 1,
115 tubeDEx.v() / dx.v() + 1,
116 tubeDEx.v() / dx.v() + 1);
117}
118
119
120void Parameters::init()
121{
122// if (tubeD.v() < pumpD.v())
123// asl::errorMessage("Tube's diameter is smaller than pump's diameter");
125}
126
128{
129 asl::SPDistanceFunction pipeGeometry;
130 asl::AVec<double> orientation(asl::makeAVec(1., 0., 0.));
131 asl::AVec<double> lVec(asl::makeAVec(params.tubeL.v()+2.*params.dx.v(), 0., 0.));
132 asl::AVec<double> h1Orientation(asl::makeAVec(0., 1., 0.));
133 asl::AVec<double> h2Orientation(asl::makeAVec(0., 0., 1.));
134 asl::AVec<double> center(asl::AVec<double>(params.size)*.5*params.dx.v());
135 double wallMid((params.tubeDEx.v()+params.tubeDIn.v())*.25);
136 double wallTh((params.tubeDEx.v())*.5);
137 asl::AVec<double> h1Center(center - (center*orientation)*orientation +
138 params.hole1Pos.v()*orientation +
139 h1Orientation*wallMid);
140 asl::AVec<double> h2Center(center - (center*orientation)*orientation +
141 params.hole2Pos.v()*orientation +
142 h2Orientation*wallMid);
143
144 pipeGeometry = asl::generateDFCylinder(params.tubeDEx.v() / 2., lVec, center) &
145 (-asl::generateDFCylinderInf(params.tubeDIn.v() / 2., orientation, center));
146 pipeGeometry = pipeGeometry &
147 (-asl::generateDFCylinder(params.hole1D.v() / 2., h1Orientation * wallTh, h1Center));
148 pipeGeometry = pipeGeometry &
149 (-asl::generateDFCylinder(params.hole2D.v() / 2., h2Orientation * wallTh, h2Center));
150 return asl::normalize(-pipeGeometry, params.dx.v());
151}
152
154{
155 float a(it<200. ? 1.-cos(it*6.28/200.) : 0);
156 return asl::makeAVec(a,0.f,0.f);
157}
158
159
160int main(int argc, char* argv[])
161{
162 Parameters params;
163 params.load(argc, argv);
164
165 std::cout << "Data initialization... " << flush;
166
167 asl::Block block(params.size, params.dx.v());
168 auto displacement(asl::generateDataContainerACL_SP<FlT>(block, 3, 1u));
169 asl::initData(displacement, asl::makeAVec(0.,0.,0.));
170
171 auto mapMem(asl::generateDataContainerACL_SP<FlT>(block, 1, 1u));
172 asl::initData(mapMem, generatePipe(block, params));
173
174
175 asl::WriterVTKXML writer(params.appParamsManager.getDir() + "acousticWaves");
176 writer.addScalars("map", *mapMem);
177 writer.addVector("displacement", *displacement);
178 writer.write();
179
180 std::cout << "Finished" << endl;
181
182 std::cout << "Numerics initialization... " << flush;
183
184 auto elasticity(generateFDElasticityRelax(displacement,
185 params.bulkMNum.v(),
186 params.shearMNum.v(),
187 params.dt.v(),
188 &asl::d3q19()));
189/* auto elasticity(generateFDElasticity(displacement,
190 params.bulkMNum.v(),
191 params.shearMNum.v(),
192 params.dt.v(),
193 &asl::d3q19()));*/
194 elasticity->setDumpingFactor(acl::generateVEConstant(.9999));
195 elasticity->init();
196
197
198 std::vector<asl::SPNumMethod> bc;
199
200 bc.push_back(generateBCZeroStress(elasticity, mapMem));
201// bc.push_back(generateBCConstantGradient(displacement,asl::makeAVec(0.,0.,0.),mapMem,&asl::d3q19()));
202// bc.push_back(generateBCConstantGradient(elasticity->getPressureData(),0.,mapMem,&asl::d3q19()));
203// bc.push_back(generateBCZeroStressP(elasticity, mapMem));
205 bc.push_back(asl::generateBCConstantValue(displacement, pres, {asl::X0}));
206
207 initAll(bc);
208
209 std::cout << "Finished" << endl;
210 std::cout << "Computing..." << endl;
211 asl::Timer timer;
212
213 executeAll(bc);
214
215 timer.start();
216 double tOutPrev(0);
217 cout << params.dt.v() << endl;
218 for (double t(0); t < params.tSimulation.v(); t+=params.dt.v())
219 {
220 elasticity->execute();
221 pres=getAmplitude(t/params.dt.v());
222 executeAll(bc);
223 if(t - params.tOutput.v()>=tOutPrev)
224 {
225 timer.stop();
226 tOutPrev=t;
227 cout << t << "/" << params.tSimulation.v() << "; time left (estimated): " << timer.estimatedRemainder(t/params.tSimulation.v()) << endl;
228 writer.write();
229 timer.start();
230
231 }
232 }
233 timer.stop();
234
235 cout << "Finished" << endl;
236
237 cout << "Computation statistic:" << endl;
238 cout << "Real Time = " << timer.realTime() << "; Processor Time = "
239 << timer.processorTime() << "; Processor Load = "
240 << timer.processorLoad() * 100 << "%" << endl;
241
242 return 0;
243}
float FlT
asl::SPDistanceFunction generatePipe(asl::Block &block, Parameters &params)
asl::AVec< float > getAmplitude(double it)
asl::Parameter< double > tubeL
asl::Parameter< double > bulkModulus
asl::Parameter< double > tubeDIn
asl::Parameter< double > tOutput
asl::Parameter< double > tSimulation
asl::Parameter< double > hole2Pos
asl::Parameter< double > shearModulus
asl::UValue< double > bulkMNum
asl::Block::DV size
asl::UValue< double > dt
asl::Parameter< double > tubeDEx
asl::Parameter< double > dx
asl::Parameter< double > hole2D
asl::ApplicationParametersManager appParamsManager
void updateNumValues()
asl::Parameter< double > hole1D
asl::Parameter< double > rho
asl::Parameter< double > hole1Pos
asl::UValue< double > shearMNum
void load(int argc, char *argv[])
void load(int argc, char *argv[])
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
const double estimatedRemainder(double completeness)
Returns estimated time till finishing current task based on its current completeness [0....
Definition aslTimer.h:52
Updatable value. This class stores value and its TimeStamp.
Definition aslUValue.h:35
const T & v() const
Definition aslUValue.h:43
void addVector(std::string name, AbstractData &data)
void addScalars(std::string name, AbstractData &data)
@ X0
Definition aslBCond.h:309
SPBCond generateBCConstantValue(SPAbstractDataWithGhostNodes d, double v, const std::vector< SlicesNames > &sl)
Bondary condition that puts fixed value in each point.
SPDistanceFunction generateDFCylinder(double r, const AVec< double > &l, const AVec< double > &c)
generates cylinder
SPDistanceFunction normalize(SPDistanceFunction a, double dx)
SPDistanceFunction generateDFCylinderInf(double r, const AVec< double > &l, const AVec< double > &c)
generates infinite cylinder
std::shared_ptr< DistanceFunction > SPDistanceFunction
Definition aslGeomInc.h:45
const VectorTemplate & d3q19()
Vector template.
VectorOfElements generateVEConstant(T a)
Generates VectorOfElements with 1 Element acl::Constant with value a.
AVec< T > makeAVec(T a1)
void initData(SPAbstractData d, double a)
int main()