Line data Source code
1 : #ifndef _Error_ESTIMATION_H_
2 : #define _Error_ESTIMATION_H_
3 :
4 : namespace coupling {
5 : namespace error {
6 :
7 : class ErrorEstimation;
8 :
9 : }
10 : } // namespace coupling
11 :
12 : /** This class provides the basic eqations for analytical error estimation based
13 : *in the statistical mechanic. There are two typs of error estimation in this
14 : *class bases on velocity ans density fluctuation. Also it can be used to
15 : *estimate the required number of MD-simulations as a function of required
16 : *maximum error.
17 : * @brief This class is used to analyse the error and predict the required
18 : *number of instances.
19 : * @author Vahid Jafarittmer
20 : */
21 : class coupling::error::ErrorEstimation {
22 : public:
23 : /** Constructor: Initializes the private member of the class.
24 : * @param velocity
25 : * @param temperature
26 : * @param numberOfParticle
27 : * @param particleMass
28 : * @param soundSpeed
29 : * @param numberOfSamples
30 : * @param cellVolume
31 : */
32 16 : ErrorEstimation(double velocity, double temperature, double numberOfParticle, double particleMass, double soundSpeed, double numberOfSamples,
33 : double cellVolume)
34 16 : : _velocity(velocity), _temperature(temperature), _numberOfParticle(numberOfParticle), _particleMass(particleMass), _soundSpeed(soundSpeed),
35 16 : _numberOfSamples(numberOfSamples), _cellVolume(cellVolume), _desiredAbsErrorVelocity(0.08), _desiredRelErrorVelocity(0.1),
36 16 : _desiredAbsErrorDensity(0.05), _desiredRelErrorDensity(0.05), _desiredAbsErrorTemperature(0.05), _desiredRelErrorTemperature(0.05),
37 16 : _desiredAbsErrorPressure(0.05), _desiredRelErrorPressure(0.05), _gamma(1.667), _C_v(1.4973), _k(1) {}
38 :
39 : /** Destructor:
40 : */
41 16 : ~ErrorEstimation() {
42 : #ifdef DEBUG_Error
43 : std::cout << " Error estimation deconstructed" << std::endl;
44 : #endif
45 16 : }
46 :
47 : /** The quantiy, its error we want to analyse
48 : * @enum errorBaseQuantity
49 : */
50 : enum errorBaseQuantity {
51 : Velocity = 0 /**< error in Velocity*/,
52 : Density = 1 /**< error in Density*/,
53 : Temperature = 2 /**< error in Temperature*/,
54 : Pressure = 3 /**< error in Pressure*/
55 : };
56 : /** relative/absolute error
57 : * @enum errorType
58 : */
59 : enum errorType {
60 : Relative = 0 /**< Relative error */,
61 : Absolute = 1 /**< Absolute error */
62 : };
63 :
64 : /** This function predict the relative or absolute error of the quantites
65 : *veocty, density, tempetature or pressure based on its arguments.
66 : * @param errorBaseQuantity
67 : * @param errorType
68 : */
69 28 : double getError(errorBaseQuantity BaseQuantity, errorType Type) {
70 :
71 28 : double error;
72 :
73 28 : if (BaseQuantity == Velocity) {
74 :
75 8 : error = getErrorVelocity(_numberOfSamples, Type == Relative ? _velocity : 1, _temperature, _numberOfParticle, _particleMass);
76 :
77 20 : } else if (BaseQuantity == Density) {
78 :
79 8 : error = getErrorDensity(_numberOfSamples, _soundSpeed, _temperature, Type == Relative ? _numberOfParticle : 1 / _numberOfParticle, _particleMass);
80 :
81 12 : } else if (BaseQuantity == Temperature) {
82 :
83 8 : error = getErrorTemperature(_numberOfSamples, _numberOfParticle, Type == Relative ? 1 : _temperature);
84 :
85 4 : } else if (BaseQuantity == Pressure) {
86 :
87 4 : error = getErrorPressure(_numberOfSamples, _numberOfParticle, _temperature, _soundSpeed, _particleMass, _cellVolume, 1);
88 :
89 : } else {
90 :
91 0 : std::cout << "ERROR coupling::ErrorEstimation(): Base quantity invalid! " << std::endl;
92 0 : exit(EXIT_FAILURE);
93 : }
94 :
95 28 : return error;
96 : }
97 :
98 : /** This function predict the required number of MD instances to keep the
99 : *relative or absolute error of the quantites veocty, density, tempetature or
100 : *pressure based on its arguments under a certain value. The desired error
101 : *value has to be set before, this function is called.
102 : * @param errorBaseQuantity
103 : * @param errorType
104 : */
105 12 : double getCorrectorNumberOfSamples(errorBaseQuantity BaseQuantity, errorType Type) {
106 :
107 12 : if (BaseQuantity == Velocity) {
108 :
109 8 : return requiredSamplesV(Type == Relative ? _desiredRelErrorVelocity : _desiredAbsErrorVelocity, _temperature, _soundSpeed,
110 8 : Type == Relative ? _velocity : 1, _numberOfParticle, _particleMass);
111 :
112 4 : } else if (BaseQuantity == Density) {
113 :
114 4 : return requiredSamplesD(_desiredRelErrorDensity, _temperature, _soundSpeed, _velocity, _particleMass, _numberOfParticle);
115 :
116 0 : } else if (BaseQuantity == Temperature) {
117 :
118 : return 3.0;
119 :
120 0 : } else if (BaseQuantity == Pressure) {
121 :
122 : return 4.0;
123 : }
124 :
125 0 : std::cout << "Error Estimation::getCorrectorNumberOfSamples not a Valid quantity" << std::endl;
126 0 : std::exit(EXIT_FAILURE);
127 : return 0.0;
128 : }
129 :
130 : /** This function estimates the velocity error
131 : * @param numberOfSamples
132 : * @param velocity
133 : * @param temperature
134 : * @param numberOfParticle
135 : * @param particleMass
136 : * @return error veocty error
137 : */
138 8 : double getErrorVelocity(double numberOfSamples, double velocity, double temperature, double numberOfParticle, double particleMass) {
139 8 : double error = 1 / (velocitySNR(velocity, temperature, numberOfParticle, particleMass) * std::sqrt(numberOfSamples));
140 8 : return error;
141 : }
142 :
143 : /** This function estimates the density error
144 : * @param numberOfSamples
145 : * @param soundSpeed
146 : * @param temperature
147 : * @param numberOfParticle
148 : * @param particleMass
149 : * @return error density error
150 : */
151 8 : double getErrorDensity(double numberOfSamples, double soundSpeed, double temperature, double numberOfParticle, double particleMass) {
152 8 : double refSP = referenceSoundSpeed(temperature, particleMass);
153 8 : double Ac = acousticNumber(soundSpeed, refSP);
154 8 : double error = 1 / std::sqrt(numberOfParticle * numberOfSamples) / Ac;
155 : // std::cout << "numberOfParticle " << numberOfParticle <<
156 : //"numberOfSamples " << numberOfSamples << "Ac " << Ac << std::endl;
157 8 : return error;
158 : }
159 :
160 : /** This function estimates the temperature error
161 : * @param numberOfSamples
162 : * @param temperature
163 : * @param numberOfParticle
164 : * @return error temperature error
165 : */
166 8 : double getErrorTemperature(double numberOfSamples, double numberOfParticle, double temperature) {
167 :
168 : // double deltaT = k*temperature*temperature/_C_v/numberOfParticle;
169 :
170 8 : double error = std::sqrt(_k / (_C_v * numberOfParticle * numberOfSamples)) * temperature;
171 8 : return error;
172 : }
173 :
174 : /** This function estimates the pressure error
175 : * @param numberOfSamples
176 : * @param numberOfParticle
177 : * @param temperature
178 : * @param soundSpeed
179 : * @param particleMass
180 : * @param cellVolume
181 : * @param pressure
182 : * @return error pressure error
183 : */
184 4 : double getErrorPressure(double numberOfSamples, double numberOfParticle, double temperature, double soundSpeed, double particleMass, double cellVolume,
185 : double pressure) {
186 :
187 4 : double Ac = acousticNumber(soundSpeed, referenceSoundSpeed(temperature, particleMass));
188 4 : double referenceP = referencePressure(temperature, numberOfParticle, cellVolume);
189 :
190 4 : double error = referenceP / pressure * Ac * std::sqrt(_gamma / (numberOfParticle * numberOfSamples));
191 4 : return error;
192 : }
193 :
194 : /** This function estimates the number of MD instances required to keep the
195 : *error of the veloctiy under _desiredAbsErrorVelocity or
196 : *_desiredrelErrorVelocity
197 : * @param desiredError
198 : * @param temperature
199 : * @param soundSpeed
200 : * @param velocity
201 : * @param numberOfParticle
202 : * @param particleMass
203 : * @return desiredNumber required number of MD instances
204 : */
205 8 : double requiredSamplesV(double desiredError, double temperature, double soundSpeed, double velocity, double numberOfParticle, double particleMass) {
206 : // double refSP= referenceSoundSpeed( temperature, particleMass);
207 : // double Ac = acousticNumber(soundSpeed, refSP);
208 8 : double SNR = velocitySNR(velocity, temperature, numberOfParticle, particleMass);
209 8 : double desiredNumber = 1 / ((SNR * desiredError) * (SNR * desiredError));
210 :
211 8 : return desiredNumber;
212 : }
213 :
214 : /** This function estimates the number of MD instances required to keep the
215 : *error of the density under _desiredAbsErrorDensity or
216 : *_desiredrelErrorDensity
217 : * @param desiredError
218 : * @param temperature
219 : * @param soundSpeed
220 : * @param velocity
221 : * @param numberOfParticle
222 : * @param particleMass
223 : * @return desiredNumber required number of MD instances
224 : */
225 4 : double requiredSamplesD(double desiredError, double temperature, double soundSpeed, double velocity, double particleMass, double numberOfParticle) {
226 4 : double refSP = referenceSoundSpeed(temperature, particleMass);
227 4 : double Ac = acousticNumber(soundSpeed, refSP);
228 4 : double desiredNumber = 1 / numberOfParticle / (desiredError * desiredError) / (Ac * Ac);
229 :
230 4 : return desiredNumber;
231 : }
232 :
233 : /** This function estimates the number of MD instances required to keep the
234 : *error of the temperature under _desiredAbsErrorTemperature or
235 : *_desiredrelErrorTemperature
236 : * @param desiredError
237 : * @param temperature
238 : * @param soundSpeed
239 : * @param velocity
240 : * @param numberOfParticle
241 : * @param particleMass
242 : * @return desiredNumber required number of MD instances
243 : */
244 : double reqiredSamplesT(double desiredError) {
245 :
246 : double desiredNumber = (_k / (_C_v * _numberOfParticle * desiredError)) * (_k / (_C_v * _numberOfParticle * desiredError));
247 :
248 : return desiredNumber;
249 : }
250 :
251 : /** signal-to-noise ratio as the average fluid velocity over the standard
252 : *deviation (square root of the veloyity fluctuation=
253 : * @param velocity
254 : * @param temperature
255 : * @param numberOfParticle
256 : * @param particleMass
257 : * @return desiredNumber required number of MD instances
258 : * @sa veloyityFluctuation
259 : */
260 16 : double velocitySNR(double velocity, double temperature, double numberOfParticle, double particleMass) {
261 :
262 16 : return velocity / veloyityFluctuation(temperature, numberOfParticle, particleMass);
263 : }
264 :
265 : /** veloyity fluctuation from equilibrium statistical mechanics
266 : * @param temperature
267 : * @param numberOfParticle
268 : * @param particleMass
269 : * @return veloyity fluctuation
270 : * @sa velocitySNR
271 : */
272 16 : double veloyityFluctuation(double temperature, double numberOfParticle, double particleMass) {
273 :
274 8 : return std::sqrt(_k * temperature / particleMass / numberOfParticle);
275 : }
276 :
277 : /** sets the required absloute error of the velocity @param error */
278 4 : void setAbsVelocityError(double error) { _desiredAbsErrorVelocity = error; }
279 :
280 : /** sets the required relative error of the velocity @param error */
281 4 : void setRelVelocityError(double error) { _desiredRelErrorVelocity = error; }
282 :
283 : /** sets the required absloute error of the density @param error */
284 : void setAbsDensityError(double error) { _desiredAbsErrorDensity = error; }
285 :
286 : /** sets the required relative error of the density @param error */
287 12 : void setRelDensityError(double error) { _desiredRelErrorDensity = error; }
288 :
289 : /** sets the required absloute error of the temperature @param error */
290 : void setAbsTemperatureError(double error) { _desiredAbsErrorTemperature = error; }
291 :
292 : /** sets the required relative error of the temperature @param error */
293 : void setRelTemperatureError(double error) { _desiredRelErrorTemperature = error; }
294 :
295 : /** sets the required absloute error of the pressure @param error */
296 : void setAbsPressureError(double error) { _desiredAbsErrorPressure = error; }
297 :
298 : /** sets the required relative error of the pressure @param error */
299 : void setRelPressureError(double error) { _desiredRelErrorPressure = error; }
300 :
301 : private:
302 : /** This function estimates the velocity error
303 : * @return velocity error
304 : */
305 : double getErrorVelocity() {
306 : double er = getErrorVelocity(_numberOfSamples, _velocity, _temperature, _numberOfParticle, _particleMass);
307 : return er;
308 : }
309 :
310 : /** This function estimates the density error
311 : * @return density error
312 : */
313 : double getErrorDensity() {
314 : double er = getErrorDensity(_numberOfSamples, _soundSpeed, _temperature, _numberOfParticle, _particleMass);
315 : return er;
316 : }
317 :
318 : /** This function estimates the temperature error
319 : * @return temperature error
320 : */
321 : double getErrorTemperature() {
322 : double er = getErrorTemperature(_numberOfSamples, _numberOfParticle, _temperature);
323 : return er;
324 : }
325 :
326 16 : double acousticNumber(double soundSpeed, double soundSpeed_ref) { return (soundSpeed / soundSpeed_ref); }
327 :
328 : //------------------------- Reference Properties
329 : //---------------------------------------------------
330 : /** calculates the sound speed of a reference ideal gas at the same
331 : * temperature
332 : * @param temperature
333 : * @param particleMass
334 : * @return referenc sound speed
335 : */
336 16 : double referenceSoundSpeed(double temperature, double particleMass) { return std::sqrt(_gamma * _k * temperature / particleMass); }
337 :
338 : /** calculates thepressure of an ideal gas under the same conditions
339 : * @param temperature
340 : * @param numberOfParticle
341 : * @param cellVolume
342 : * @return reference pressure
343 : */
344 4 : double referencePressure(double temperature, double numberOfParticle, double cellVolume) {
345 4 : return std::sqrt(numberOfParticle * _k * temperature / cellVolume);
346 : }
347 :
348 : double _velocity;
349 : double _temperature;
350 : double _numberOfParticle;
351 : double _particleMass;
352 : double _soundSpeed;
353 : double _numberOfSamples;
354 : double _cellVolume;
355 :
356 : double _desiredAbsErrorVelocity; //=0.08;
357 : double _desiredRelErrorVelocity; // =0.1;
358 :
359 : double _desiredAbsErrorDensity; //=0.05;
360 : double _desiredRelErrorDensity; //=0.05;
361 :
362 : double _desiredAbsErrorTemperature; //=0.05;
363 : double _desiredRelErrorTemperature; //=0.05;
364 :
365 : double _desiredAbsErrorPressure; //=0.05;
366 : double _desiredRelErrorPressure; //=0.05;
367 :
368 : double _gamma;
369 : double _C_v;
370 : /** Boltzman Constant unit: [J/K] */
371 : double _k;
372 : };
373 :
374 : #endif
|
