My Project
simulation.cpp
Go to the documentation of this file.
1 #include"simulation.hpp"
2 #include<omp.h>
3 
13 Simulation::Simulation(string parameters_path):input_parameters_path(parameters_path)
14 {
15  input_parameters = Parameters(parameters_path);
16  //arma::arma_rng::set_seed_random();
17  arma::arma_rng::set_seed(input_parameters.seed);
20 
21 }
22 
23 Simulation::Simulation(string parameters_path, bool friction_parallelized):input_parameters_path(parameters_path)
24 {
25  if(!friction_parallelized)
26  {
27  cout << "Called with friction_parallized parameter, but it is set to false." << endl;
28  cout << "Please run without this parameter if the friction parallization version is not wanted." << endl;
29  }
30  else
31  {
32  cout << "Running the friction measurement in parallel for different slider speeds." << endl;
33  input_parameters = Parameters(parameters_path);
34  //arma::arma_rng::set_seed_random();
35  arma::arma_rng::set_seed(input_parameters.seed);
38  }
39 
40 
41 }
42 
43 // START: SET UP VECTORS AND MATRICES //
45 {
46 
49 
52 
53 
55 }
56 
58 {
59  block_matrix.set_size(input_parameters.N, 2);
61  block_matrix = 2.0*block_matrix.randu() - 1.0;
62  }
63 
65  {
66  block_matrix = block_matrix.randu() * (-1.0);
67  }
68 
69  else
70  {
71  block_matrix.randu();
72  }
73 }
74 
76 {
78  block_matrix_derivative.fill(1.0);
79 }
80 
82 {
83  pad_vector = arma::zeros<arma::rowvec>(2);
84 }
85 
87 {
88  pad_vector_derivative = arma::zeros<arma::rowvec>(2);
89 }
90 
92 {
93  // SAVE PURPOSES //
98  {
99  loggingSpecificBlocks = true;
100  }
102  {
105  } else
106  {
109  }
110  printf("Number of time steps: %i\n",number_of_save_steps);
111  printf("Logging specific blocks: %d\n", loggingSpecificBlocks);
112  block_position_to_file.print();
113 
117  {
119  }else
120  {
122 
123  }
124 
125 }
126 
128 {
130  {
131  if(input_parameters.interval == 0)
132  {
136  }
137  else
138  {
139  updateSliderSpeedDtInterval = (int)(input_parameters.interval/input_parameters.dt); // used for updating slider speed
140  }
141 
142  }else
143  {
144  updateSliderSpeedDtInterval = (int)(input_parameters.interval/input_parameters.dt); // used for updating slider speed
145  }
146  printf("Slider speed increment: %.10f\n",input_parameters.increment);
147 }
148 
149 // END: SET UP VECTORS AND MATRICES //
150 
151 // START: CONSTANTS //
152 
154 {
156  k_p = input_parameters.m_scale_P*input_parameters.m_k_P0/input_parameters.N; // Stationary spring between pad and blocks
157  m_c_crit = 2*sqrt((input_parameters.m_k_P0+input_parameters.N*k_p)*input_parameters.m_mass_x); // Critical damping coefficient
158  c_p = input_parameters.m_zeta*m_c_crit; // Damper between pad and surface
160 
162 }
163 
164 // END: CONSTANTS //
165 
166 // START: HELP FUNCTIONS CALCULATION //
167 
178 double Simulation::phi(double y)
179 {
180 
181  if (y == 0){
183  printf("[Friction zero counter: %i]", friction_zero_counter);
184  return 0; // Maximum static friction force is handled elsewhere
185  }
186 
187  return (y < 0 ? -1 : 1)*(1 - input_parameters.m_sigma) / (1 + abs(y) / (1 - input_parameters.m_sigma));
188 
189 }
190 
192 {
194  {
196  }
197 }
198 
199 
201  int time_step
202 )
203 {
204  if(time_step % updateSliderSpeedDtInterval == 0)
205  {
206  printf("[Slider speed: %.6f]",input_parameters.slider_speed);
209  {
210  if(time_step % input_parameters.save_interval_dt == 0)
211  {
213  }
215  cout << input_parameters.slider_speed << endl;
216  }
218  }
219 }
220 
221 // END: HELP FUNCTIONS CALCULATION //
222 
223 // START: F FUNCTIONS NUMERICAL SCHEME //
241  const arma::mat& current_block_matrix,
242  const arma::rowvec& current_pad_vector,
243  arma::mat& current_block_matrix_derivative)
244 {
245 
246  arma::colvec spring_forces = arma::zeros<arma::colvec>(input_parameters.N);
247 
248  for (int i = 0; i < input_parameters.N; i++) {
249 
251  {
252  // Add contribution from neighboring springs, if not deactivated
253  spring_forces(i)
254  = k_c*(i == 0 ? 0 : (current_block_matrix(i - 1, 0) - current_block_matrix(i, 0)))
255  + k_c*(i == input_parameters.N - 1 ? 0 : (current_block_matrix(i + 1, 0) - current_block_matrix(i, 0)));
256 
257  }
258 
260  {
261  // Add contribution from stationary springs, if not deactivated
262  spring_forces(i) -= k_p*(current_block_matrix(i, 0) - current_pad_vector(0));
263  }
264 
266  {
267  // Add friction force
268  spring_forces(i) -= m_u*phi(2*input_parameters.m_alpha*input_parameters.slider_speed + current_block_matrix(i, 1));
269  }
270  }
271 
272  current_block_matrix_derivative.col(0) = current_block_matrix.col(1);
273  current_block_matrix_derivative.col(1) = spring_forces / m_u;
274 }
275 
290  const arma::mat& current_block_matrix,
291  const arma::rowvec& current_pad_vector,
292  arma::rowvec& current_pad_vector_derivative)
293 {
294 
295  double pad_second_derivative = 0;
296 
298  {
299  pad_second_derivative = calculatePadSupportDamper(current_pad_vector);
300  }
301 
303  {
304  pad_second_derivative += calculatePadFriction(current_block_matrix, current_pad_vector);
305  }
306 
307 
308  pad_second_derivative += -input_parameters.m_k_P0*current_pad_vector(0);
309 
310  current_pad_vector_derivative(0) = current_pad_vector(1);
311  current_pad_vector_derivative(1) = pad_second_derivative / input_parameters.m_mass_x;
312 
313  /* cout << pad_second_derivative << endl;
314  current_pad_vector.print("current_pad_vector");
315  pad_vector_derivative.print("pad_vector_derivative"); */
316 
317 }
318 
332 // ONE DEGREE-OF-FREEDOM //
334  const arma::rowvec& current_pad_vector,
335  arma::rowvec& current_pad_vector_derivative
336 ){
337  double pad_second_derivative = 0;
338 
340  {
341  pad_second_derivative = calculatePadSupportDamper(current_pad_vector);
342  }
343 
345  {
346  pad_second_derivative += -phi(current_pad_vector(1) + input_parameters.slider_speed);
347  }
348 
349 
350  pad_second_derivative += -input_parameters.m_k_P0*current_pad_vector(0);
351 
352  current_pad_vector_derivative(0) = current_pad_vector(1);
353  current_pad_vector_derivative(1) = pad_second_derivative / input_parameters.m_mass_x;
354 }
355 
356 double Simulation::calculatePadSupportDamper(arma::rowvec current_pad_vector)
357 {
358  return -c_p*current_pad_vector(1);
359 }
360 
361 double Simulation::calculatePadFriction(arma::mat current_block_matrix,
362  arma::rowvec current_pad_vector)
363 {
364  return k_p*arma::sum(current_block_matrix.col(0) - current_pad_vector(0));
365 }
366 
367 
386 void Simulation::midpointMethod(int time_step)
387 {
388  arma::mat block_matrix_half_step = block_matrix;
389  arma::rowvec pad_vector_half_step = pad_vector;
390 
393 
394  block_matrix_half_step = block_matrix + 0.5 * input_parameters.dt * block_matrix_derivative;
395  pad_vector_half_step = pad_vector + 0.5 * input_parameters.dt * pad_vector_derivative;
396 
397  function_u(block_matrix_half_step, pad_vector_half_step, block_matrix_derivative);
398  function_x(block_matrix_half_step, pad_vector_half_step, pad_vector_derivative);
399 
402 
403  if(time_step % input_parameters.save_interval_dt == 0)
404  {
405  writeValuesToLoggers(time_step);
406  }
407 }
408 
410  int time_step
411 ){
412  arma::rowvec pad_vector_half_step = pad_vector;
413 
415 
416  pad_vector_half_step = pad_vector + 0.5 * input_parameters.dt * pad_vector_derivative;
417 
418  function_x_odof(pad_vector_half_step, pad_vector_derivative);
419 
421 
422  if(time_step % input_parameters.save_interval_dt == 0)
423  {
424  writeValuesToLoggers(time_step);
425  }
426 
427 }
428 
429 void Simulation::writeValuesToLoggers(const int& time_step){
431  {
433  {
436  }else
437  {
440  }
443  }
445  index_to_file +=1;
447 }
448 
449 void Simulation::printMidPointMethod(arma::mat block_matrix, arma::rowvec pad_vector, string explanation_text)
450 {
451  block_matrix.print("Block matrix" + explanation_text + ": ");
452  pad_vector.print("Pad vector" + explanation_text + ": ");
453 }
454 
456 {
457  double percentageTransformer = 100.0/(double)number_of_time_steps;
458  cout << "\n\n" << endl;
459  for(int i = 0; i < number_of_time_steps; i++)
460  {
461  printf("\r[Percentage complete: %.0f%%]",(double)i*percentageTransformer);
463  {
464  midpointMethod(i);
465  }else
466  {
468  }
469  updateSliderSpeed(i+1);
470  }
471  cout << "\n\n" << endl;
472 
473  saveToCsv(input_parameters.file_name, "results/");
474 }
475 
476 void Simulation::saveToCsv(string filename, string result_path = "results/")
477 {
479  {
480  block_position_to_file.save(result_path + filename + "_midpoint_block_position.csv", arma::csv_ascii);
481  block_velocity_to_file.save(result_path + filename + "_midpoint_block_velocity.csv", arma::csv_ascii);
482  pad_position_to_file.save(result_path + filename + "_midpoint_pad_position.csv", arma::csv_ascii);
483  pad_velocity_to_file.save(result_path + filename + "_midpoint_pad_velocity.csv", arma::csv_ascii);
484  }
486  {
487  pad_friction_to_file.save(result_path + filename + "_midpoint_pad_friction.csv", arma::csv_ascii);
488  }
489 }
490 
491 // START: Parallel friction calculation //
492 
494  int timestep,
495  arma::mat& block_matrix_parallel,
496  arma::rowvec& pad_vector_parallel,
497  arma::mat block_matrix_derivative_parallel,
498  arma::rowvec pad_vector_derivative_parallel
499 )
500 {
501  arma::mat block_matrix_half_step = block_matrix_parallel;
502  arma::rowvec pad_vector_half_step = pad_vector_parallel;
503 
504  // HERE NOW
505 
506 
507  function_u(block_matrix_parallel, pad_vector_parallel, block_matrix_derivative_parallel);
508  function_x(block_matrix_parallel, pad_vector_parallel, pad_vector_derivative_parallel);
509 
510  block_matrix_half_step = block_matrix_parallel + 0.5 * input_parameters.dt * block_matrix_derivative_parallel;
511  pad_vector_half_step = pad_vector_parallel + 0.5 * input_parameters.dt * pad_vector_derivative_parallel;
512 
513  function_u(block_matrix_half_step, pad_vector_half_step, block_matrix_derivative_parallel);
514  function_x(block_matrix_half_step, pad_vector_half_step, pad_vector_derivative_parallel);
515 
516  block_matrix_parallel = block_matrix_parallel + input_parameters.dt * block_matrix_derivative_parallel;
517  pad_vector_parallel = pad_vector_parallel + input_parameters.dt * pad_vector_derivative_parallel;
518 }
519 
521  int time_step,
522  double slider_speed)
523 {
524  arma::mat block_matrix_parallel;
525  arma::mat block_matrix_derivative_parallel;
526 
527  arma::rowvec pad_vector_parallel;
528  arma::rowvec pad_vector_derivative_parallel;
529  block_matrix_parallel.set_size(input_parameters.N, 2);
531  block_matrix_parallel = 2.0*block_matrix_parallel.randu() - 1.0;
532  }
533 
535  {
536  block_matrix_parallel = block_matrix_parallel.randu() * (-1.0);
537  }
538 
539  else
540  {
541  block_matrix_parallel.randu();
542  }
543  block_matrix_derivative_parallel.set_size(input_parameters.N, 2);
544  block_matrix_derivative_parallel.fill(1.0);
545 
546  pad_vector_parallel = arma::zeros<arma::rowvec>(2);
547  pad_vector_derivative_parallel = arma::zeros<arma::rowvec>(2);
548 
549  arma::colvec pad_friction_parallel_to_file;
550  pad_friction_parallel_to_file.set_size(time_step);
551 
552  for(int i = 0; i < time_step; i++)
553  {
554  midpointMethodParallel(time_step, block_matrix_parallel, pad_vector_parallel, block_matrix_derivative_parallel, pad_vector_derivative_parallel);
555  pad_friction_parallel_to_file(i) = calculatePadFriction(block_matrix_parallel, pad_vector_parallel);
556  }
557 
558  saveToCsvParallel(pad_friction_parallel_to_file, "results/", slider_speed);
559 }
560 
562  const arma::colvec& armadillo_vector,
563  const string& result_path,
564  const double& slider_speed
565  )
566 {
567  armadillo_vector.save(result_path + input_parameters.file_name + "_midpoint_pad_friction" + to_string(slider_speed) + ".csv", arma::csv_ascii);
568 }
569 
571 {
572  int number_of_slider_speed_intervals = floor(input_parameters.max_time / input_parameters.interval);
573  cout << number_of_slider_speed_intervals;
574  double slider_speed_parallel = input_parameters.slider_speed;
575 
576  double time_steps_per_slider_speed = input_parameters.interval / input_parameters.dt;
577 
578  omp_set_num_threads(4);
579  #pragma omp parallel for
580  for(int i = 1; i < number_of_slider_speed_intervals + 1; i++)
581  {
582  slider_speed_parallel = input_parameters.slider_speed * i;
583  setupAndRunFrictionInParallel(time_steps_per_slider_speed, slider_speed_parallel);
584  }
585 
586 }
587 
588 // END: Parallel friction calculation //
589 
590 // END: F FUNCTIONS NUMERICAL SCHEME //
591 
593 {
595 }
596 
597 // DEBUGGING FUNCTIONS //
598 
599 void DebugSimulation::printMatrix(arma::mat matrix, bool array_format = false)
600 {
601  cout << endl;
602  if(!array_format)
603  {
604  matrix.print("Matrix: ");
605  }else
606  {
607  int count = 0;
608  arma::mat::iterator it_end = matrix.end();
609  for(arma::mat::iterator it = matrix.begin(); it != it_end; ++it)
610  {
611  cout << it[0] << endl;
612  count += 1;
613  }
614  cout << count << endl;
615  }
616  cout << endl;
617 }
618 
619 void DebugSimulation::printRowVector(arma::rowvec vector){
620  cout << endl;
621  vector.print("Row vector: ");
622  cout << endl;
623 }
624 
626 {
628  /* printMatrix(block_matrix);
629  printRowVector(pad_vector); */
630  // function_u(block_matrix, pad_vector);
631  // function_x(block_matrix, pad_vector);
632  /* printMatrix(block_matrix_derivative);
633  printRowVector(pad_vector_derivative); */
635 }
636 
637 // START: EXAMPLE //
638 // END: EXAMPLE //
double increment
void midpointMethodParallel(int timestep, arma::mat &block_matrix_parallel, arma::rowvec &pad_vector_parallel, arma::mat block_matrix_derivative_parallel, arma::rowvec pad_vector_derivative_parallel)
Definition: simulation.cpp:493
void function_u(const arma::mat &current_block_matrix, const arma::rowvec &current_pad_vector, arma::mat &current_block_matrix_derivative)
Definition: simulation.cpp:240
void setupBlockMatrixDerivative()
Definition: simulation.cpp:75
arma::mat block_matrix_derivative
Definition: simulation.hpp:99
void function_x_odof(const arma::rowvec &current_pad_vector, arma::rowvec &current_pad_vector_derivative)
Definition: simulation.cpp:333
double calculatePadFriction(arma::mat current_block_matrix, arma::rowvec current_pad_vector)
Definition: simulation.cpp:361
double k_p
Definition: simulation.hpp:109
bool debug_no_neighbor_springs
arma::uvec block_velocity_column_index
Definition: simulation.hpp:118
double m_scale_P
double m_scale_C
double threshold_speed
int friction_zero_counter
Definition: simulation.hpp:129
void calculateFrictionInParallel()
Definition: simulation.cpp:570
arma::colvec pad_friction_to_file
Definition: simulation.hpp:124
arma::mat block_velocity_to_file
Definition: simulation.hpp:121
arma::mat block_position_to_file
Definition: simulation.hpp:120
arma::uvec loggedBlocks
double phi(double y)
Definition: simulation.cpp:178
void saveToCsvParallel(const arma::colvec &armadillo_vector, const string &result_path, const double &slider_speed)
Definition: simulation.cpp:561
arma::rowvec pad_vector_derivative
Definition: simulation.hpp:102
double calculatePadSupportDamper(arma::rowvec current_pad_vector)
Definition: simulation.cpp:356
string input_parameters_path
Definition: simulation.hpp:94
void function_x(const arma::mat &current_block_matrix, const arma::rowvec &current_pad_vector, arma::rowvec &current_pad_vector_derivative)
Definition: simulation.cpp:289
void setupConstants()
Definition: simulation.cpp:153
double m_u
Definition: simulation.hpp:108
void setupBlockMatrix()
Definition: simulation.cpp:57
Parameters input_parameters
Definition: simulation.hpp:96
void setupPadVectorDerivative()
Definition: simulation.cpp:86
void printMatrix(arma::mat matrix, bool arra_format)
Definition: simulation.cpp:599
int index_pad_friction_to_file
Definition: simulation.hpp:126
bool debug_only_negative_initial
double slider_speed
arma::colvec pad_position_to_file
Definition: simulation.hpp:122
void main()
Definition: simulation.cpp:592
int updateSliderSpeedDtInterval
Definition: simulation.hpp:114
void midpointMethod(int time_step)
Definition: simulation.cpp:386
int number_of_time_steps
Definition: simulation.hpp:104
void setupPadVector()
Definition: simulation.cpp:81
void updateSliderSpeed(int time_step)
Definition: simulation.cpp:200
double end_speed_continuous
bool debug_continuous_slider_speed
void printMidPointMethod(arma::mat block_matrix, arma::rowvec pad_vector, string explanation_text="")
Definition: simulation.cpp:449
int number_of_save_steps
Definition: simulation.hpp:105
void setupAndRunFrictionInParallel(int time_step, double slider_speed)
Definition: simulation.cpp:520
arma::rowvec pad_vector
Definition: simulation.hpp:101
double k_c
Definition: simulation.hpp:112
void oneDegreeOfFreedomMidpointMethod(int time_step)
Definition: simulation.cpp:409
arma::colvec pad_velocity_to_file
Definition: simulation.hpp:123
bool debug_one_degree_freedom_mode
bool debug_only_write_friction
void saveToCsv(string filename, string result_path)
Definition: simulation.cpp:476
bool debug_no_stationary_springs
void runMidpointMethod()
Definition: simulation.cpp:455
int number_of_interval_save_step
Definition: simulation.hpp:106
arma::uvec block_position_column_index
Definition: simulation.hpp:117
double start_speed_continuous
Simulation(string parameters_path)
Definition: simulation.cpp:13
void writeValuesToLoggers(const int &time_step)
Definition: simulation.cpp:429
void printRowVector(arma::rowvec vector)
Definition: simulation.cpp:619
bool debug_negative_intial_values
double m_c_crit
Definition: simulation.hpp:110
void setSpeedAfterThreshold()
Definition: simulation.cpp:191
string file_name
bool debug_no_friction
bool debug_no_damper
double m_scale_mass
double c_p
Definition: simulation.hpp:111
void setupSaveVectorsMatrices()
Definition: simulation.cpp:91
int index_to_file
Definition: simulation.hpp:127
void setupSliderSpeedMode()
Definition: simulation.cpp:127
void setupModelVectorsMatrices()
Definition: simulation.cpp:44
arma::mat block_matrix
Definition: simulation.hpp:98
bool loggingSpecificBlocks
Definition: simulation.hpp:116