diff --git a/main.cpp b/main.cpp
index c66983d48ee620d879e46f7b06fa462798404f10..f5be00a6630acaaac22b66b6a24f8f5ccb46d7f5 100644
--- a/main.cpp
+++ b/main.cpp
@@ -1,5 +1,7 @@
#include <iostream>
+#include <fstream>
#include <vector>
+#include <string>
#include <cstring>
#include "sound_column.h"
#include <execution>
@@ -11,30 +13,41 @@
#define DEFAULT_MAX_FREQ 12'500.0f
#define WATER_SPEED 1'500.f
-std::complex<float> gaussian_init(float, float);
-float media(float, float);
-int run_simulation(float, float, float, float);
+void gaussian_init(sound_column&);
+bool step_r(sound_column&);
+void _step_r_(sound_column&, size_t);
+float media(float, float, float);
+std::vector<sound_column*>* run_simulation(float, float, float, float);
+int save_results(std::string&, std::vector<sound_column*>&);
int main(int argc, char **argv) {
float r_max = DEFAULT_DISTANCE, z_max = DEFAULT_DEPTH, f_min = DEFAULT_MIN_FREQ, f_max = DEFAULT_MAX_FREQ;
+ bool dry_run = false;
+ std::string save_path = "";
int i=1;
while (i < argc){
- if (strcmp(argv[i], "--d_max") != 0 || strcmp(argv[i], "-d") != 0)
+ if (strcmp(argv[i], "--r_max") == 0 || strcmp(argv[i], "-r") == 0)
r_max = atof(argv[++i]);
- else if (strcmp(argv[i], "--f_max") != 0 || strcmp(argv[i], "-F") != 0)
+ else if (strcmp(argv[i], "--f_max") == 0 || strcmp(argv[i], "-F") == 0)
f_max = atof(argv[++i]);
- else if (strcmp(argv[i], "--f_min") != 0 || strcmp(argv[i], "-f") != 0)
+ else if (strcmp(argv[i], "--f_min") == 0 || strcmp(argv[i], "-f") == 0)
f_min = atof(argv[++i]);
- else if (strcmp(argv[i], "--z_max") != 0 || strcmp(argv[i], "-z") != 0)
+ else if (strcmp(argv[i], "--z_max") == 0 || strcmp(argv[i], "-z") == 0)
z_max = atof(argv[++i]);
- else if (strcmp(argv[i], "--help") != 0 || strcmp(argv[i], "-h") != 0){
+ else if (strcmp(argv[i], "--save_path") == 0 || strcmp(argv[i], "-s") == 0)
+ save_path = argv[++i];
+ else if (strcmp(argv[i], "--dry_run") == 0)
+ dry_run = true;
+ else if (strcmp(argv[i], "--help") == 0 || strcmp(argv[i], "-h") == 0){
std::cout << "Simulation of wave propagation using the parabolic equation." << std::endl
<< "--r_max, -r : Set the maximum range of the simulation. Default :" << DEFAULT_DISTANCE << std::endl
<< "--f_max, -F : Set the maximum frequency of the simulation. Default :" << DEFAULT_MAX_FREQ << std::endl
- << "--f_min, -f : Set the minimum frequency of the simulation. Default :" << DEFAULT_MIN_FREQ << std::endl
- << "--z_max, -z : Set the maximum depth of the simulation. Default :" << DEFAULT_DEPTH << std::endl;
+ << "--f_min, -f : Set the minimum frequency and the frequency step of the simulation. Default :" << DEFAULT_MIN_FREQ << std::endl
+ << "--z_max, -z : Set the maximum depth of the simulation. Default :" << DEFAULT_DEPTH << std::endl
+ << "--save_path, -s : Set the path to the output file. Default will not save anything" << std::endl
+ << "--dry_run : Run the simulation but does not save anything, even if save_path is set" << std::endl;
return 1;
}
else {
@@ -43,40 +56,109 @@ int main(int argc, char **argv) {
i++;
}
+ std::vector<sound_column*>* results(run_simulation(r_max, z_max, f_min, f_max));
- return run_simulation(r_max, z_max, f_min, f_max);
+ if (not dry_run && save_path.length()){
+ std::cout << "Saving results to " << save_path << std::endl;
+ return save_results(save_path, *results);
+ }
+ else
+ return 0;
+}
+
+
+std::vector<sound_column*>* run_simulation(float r_max, float z_max, float f_min, float f_max){
+ auto nf = size_t(f_max / f_min);
+ float delta_r = r_max/10'000;
+
+ static std::vector<sound_column*> column(nf, nullptr);
+
+ //std::ranges::iota_view indexes((size_t) 0, nf);
+ std::cout << "Initializing " << nf << " pressure columns" << std::endl;
+ std::for_each(std::execution::par, column.begin(), column.end(), //par_unseq does not compile
+ [&](sound_column*& c) {int i = &c - &column[0];
+ c = new sound_column(delta_r, (size_t) 1.10 * z_max*(4*(i+1)*f_min)/WATER_SPEED, 7, r_max, z_max, (i+1)*f_min, WATER_SPEED);
+ gaussian_init(*c);});
+ std::cout << "Running the simulation" << std::endl;
+ std::for_each(std::execution::par, column.begin(), column.end(),
+ [&](sound_column*& c) {while (not step_r(*c));std::cout << '.' << std::flush;});
+ std::cout << std::endl << "Done" << std::endl;
+ return &column;
}
-std::complex<float> gaussian_init(float depth, float k0){
- return sqrtf(k0 / 0.6f) * expf(-powf((k0 / 0.6f) * (depth - 1'000.f), 2));
+void gaussian_init(sound_column& column){
+ float k0 = 2*M_1_PIf32 * column.freq / column.c0;
+ std::for_each(std::execution::par_unseq, column.pressures.begin(), column.pressures.end(),
+ [&, k0](std::complex<float>& p) {int i = &p - &column.pressures[0];
+ p = sqrtf(k0 / 0.6f) * expf(-powf((k0 / 0.6f) * (column.z_max/column.n_z*i - 1'000.f), 2));});
}
-float media(float r, float depth){
- static const float slope = -0.009f;
+bool step_r(sound_column& column) {
+ if (column.r >= column.d_max)
+ return true;
+ //std::ranges::iota_view z_ind((size_t) 0, n_z);
+ column.r += column.delta_r;
+ std::for_each(std::execution::par_unseq, column.n_next.begin(), column.n_next.end(), //par_unseq does not compile
+ [&](float& n) {int i = &n - &column.n_next[0]; n = media(column.r, column.z_max/column.n_z*i, column.c0);});
+ for(size_t i=0; i < column.n_p;i++){
+ _step_r_(column, i);
+ }
+ // Hank
+ std::for_each(std::execution::par_unseq, column.pressures.begin(), column.pressures.end(),
+ [&](std::complex<float>& p) {p *= std::complex<float>(2.f / (M_1_PIf32 * M_1_PIf32 * column.freq / column.c0 * column.r), 0.0f) *
+ std::exp(std::complex<float>(0.0f, 2.f* M_1_PIf32 * column.freq / column.c0 * column.r - M_1_PIf32 / 4));});
+ return false;
+}
+
+void _step_r_(sound_column& column, size_t ip) {
+// std::ranges::iota_view z_ind((size_t) 0, n_z);
+ //tri diag init
+ std::fill(std::execution::par_unseq, column._u_.begin(), column._u_.end(), column.fbmaj[ip]/column.h2);
+ std::fill(std::execution::par_unseq, column._l_.begin(), column._l_.end(), column.fbmaj[ip]/column.h2);
+
+ float k02 = powf(2*M_1_PIf32 * column.freq / column.c0, 2);
+ std::for_each(std::execution::par_unseq, column._d_.begin(), column._d_.end(),
+ [&, k02, ip](std::complex<float>& d) {int i = &d - &column._d_[0];
+ d = std::complex<float>(2 * k02, 0.0f) + column.fbmaj[ip] /
+ std::complex<float>(-2.f/column.h2 + k02 * (powf(column.n_next[i], 2) - 1.f), 0.0f);});
+
+ // first matmul
+
+ // Inv tri matmul
+ column._u_[0] /= column._d_[0];
+ column.pressures[0] /= column._d_[0];
+ for(size_t i=1; i < column.n_z-1; i++){
+ column._u_[i] /= column._d_[i] - column._l_[i-1] * column._u_[i-1];
+ column.pressures[i] = (column.pressures[i] - column._l_[i-1] * column.pressures[i-1]) / (column._d_[i] - column._l_[i-1] * column._u_[i-1]);
+ }
+ column.pressures[column.n_z-1] = (column.pressures[column.n_z-1] - column._l_[column.n_z-2] * column.pressures[column.n_z-2]) / (column._d_[column.n_z-1] - column._l_[column.n_z-2] * column._u_[column.n_z-2]);
+ for(ssize_t i=column.n_z-2; i >= 0; i--) {
+ column.pressures[i] -= column._u_[i]*column.pressures[i+1];
+ }
+}
+
+float media(float r, float depth, float speed){
+ const float slope = -0.009f;
if ( depth > slope * r + 2531.f)
- return WATER_SPEED / 1650.f;
+ return speed / 1650.f;
else if ( depth > slope/2 * r + 3042.f)
- return WATER_SPEED / 1700.f;
+ return speed / 1700.f;
else if ( depth > slope/3 * r + 4111.f)
- return WATER_SPEED / 1800.f;
+ return speed / 1800.f;
else
- return WATER_SPEED / (1500.f*(1.f+0.00737f*(depth-1.f+expf(-depth))));
+ return speed / (1500.f*(1.f+0.00737f*(depth-1.f+expf(-depth))));
}
-int run_simulation(float r_max, float z_max, float f_min, float f_max){
- auto nf = size_t(f_max / f_min);
- float delta_r = r_max/10'000;
-
- std::vector<sound_column*> column(nf, nullptr);
-
- //std::ranges::iota_view indexes((size_t) 0, nf);
+int save_results(std::string& path, std::vector<sound_column*>& results){
+ std::fstream file;
+ file.open(path, std::fstream::out | std::fstream::trunc);
+ for(auto c:results){
+ for(auto val:c->pressures){
+ file << val << ',';
+ }
+ file << std::endl;
+ }
- std::for_each(std::execution::par_unseq, column.begin(), column.end(),
- [&](sound_column*& c) {int i = &c - &column[0];
- c = new sound_column(delta_r, (size_t) 1.10 * z_max*(4*(i+1)*f_min)/WATER_SPEED, 7, r_max, z_max, (i+1)*f_min, WATER_SPEED);
- c->init_source(gaussian_init);});
- std::for_each(std::execution::par_unseq, column.begin(), column.end(),
- [&](sound_column*& c) {while (c->step_r(media));});
return 0;
}
diff --git a/sound_column.cpp b/sound_column.cpp
index 03724ad53080fb73db9de52c5dedbaabe84f91ec..ded23b1b256f3370adf3383877ccafa85015d880 100644
--- a/sound_column.cpp
+++ b/sound_column.cpp
@@ -21,7 +21,7 @@ sound_column::sound_column(float dr, size_t nz, size_t np, float dm, float zm, f
n_next.resize(nz);
fapbj.resize(np);
fbmaj.resize(np);
- for(int i=0; i < np; i++) {
+ for(size_t i=0; i < np; i++) {
faj = std::complex<float>(0, 4* M_1_PIf32 * freq / c0) / float((2.f * n_p + 1.f) *
powf(sinf(float(i+1) * M_1_PIf32 / (2.f * n_p + 1.f)), 2) * delta_r);
fbj = 2 * powf(cosf(float(i+1) * M_1_PIf32 / (2.f * n_p + 1.f)),2);
@@ -31,41 +31,28 @@ sound_column::sound_column(float dr, size_t nz, size_t np, float dm, float zm, f
}
}
+sound_column::sound_column() = default;
+
sound_column::~sound_column() = default;
-bool sound_column::step_r(float (*fun)(float, float)) {
- if (r >= d_max)
- return true;
- //std::ranges::iota_view z_ind((size_t) 0, n_z);
- r += delta_r;
- std::for_each(std::execution::par_unseq, n_next.begin(), n_next.end(),
- [&](float& n) {int i = &n - &n_next[0]; n = (*fun)(r, z_max/n_z*i);});
- for(int i=0; i < n_p;i++){
- _step_r_(i);
- }
- // Hank
- std::for_each(std::execution::par_unseq, pressures.begin(), pressures.end(),
- [&](std::complex<float>& p) {p *= static_cast<std::complex<float>>(2.f / (M_1_PIf32 * M_1_PIf32 * freq / c0 * r)) *
- std::exp(std::complex<float>(0, 2.f* M_1_PIf32 * freq / c0 * r - M_1_PIf32 / 4));});
- return false;
-}
-void sound_column::init_source(std::complex<float> (*fun)(float, float)) {
-// std::ranges::iota_view z_ind((size_t) 0, n_z);
- std::for_each(std::execution::par_unseq, pressures.begin(), pressures.end(),
- [&](std::complex<float>& p) {int i = &p - &pressures[0];
- p = (*fun)(z_max/n_z*i, 2*M_1_PIf32 * freq / c0);});
-}
+//void sound_column::init_source(std::complex<float> (*fun)(float, float)) {
+//// std::ranges::iota_view z_ind((size_t) 0, n_z);
+// std::for_each(std::execution::par_unseq, pressures.begin(), pressures.end(),
+// [&](std::complex<float>& p) {int i = &p - &pressures[0];
+// p = (*fun)(z_max/n_z*i, 2*M_1_PIf32 * freq / c0);});
+//}
-void sound_column::_step_r_(int ip) {
+/*void sound_column::_step_r_(size_t ip) {
// std::ranges::iota_view z_ind((size_t) 0, n_z);
//tri diag init
std::fill(std::execution::par_unseq, _u_.begin(), _u_.end(), fbmaj[ip]/h2);
std::fill(std::execution::par_unseq, _l_.begin(), _l_.end(), fbmaj[ip]/h2);
+ float double_k0 = 4*M_1_PIf32 * freq / c0;
std::for_each(std::execution::par_unseq, _d_.begin(), _d_.end(),
- [&](std::complex<float>& d) {int i = &d - &_d_[0];d = static_cast<std::complex<float>>(4* M_1_PIf32 * freq / c0) + fbmaj[ip] /
- static_cast<std::complex<float>>(-2.f/h2 + powf(4* M_1_PIf32 * freq / c0, 2) * (powf(n_next[i], 2) - 1.f));});
+ [&, double_k0, ip](std::complex<float>& d) {int i = &d - &_d_[0];d = std::complex<float>(double_k0,0.0f) + fbmaj[ip] /
+ std::complex<float>(-2.f/h2 + powf(double_k0, 2) * (powf(n_next[i], 2) - 1.f), 0.0f);});
// first matmul
@@ -80,4 +67,4 @@ void sound_column::_step_r_(int ip) {
for(ssize_t i=n_z-2; i >= 0; i--) {
pressures[i] -= _u_[i]*pressures[i+1];
}
-}
+}*/
diff --git a/sound_column.h b/sound_column.h
index 23049d87a0b584eb9849c1306543cb1274c6a78a..76243a9ef00b82e3f0486f37e976ab04c26aa647 100644
--- a/sound_column.h
+++ b/sound_column.h
@@ -10,38 +10,36 @@
#include <execution>
#include <algorithm>
-# define M_1_PIf32 (0.318309886183790671537767526745028724)
+#ifndef M_1_PIf32
+#define M_1_PIf32 (0.318309886183790671537767526745028724)
+#endif
class sound_column {
-private:
- float delta_r;
- size_t n_z;
- size_t n_p;
+public:
+ float delta_r{};
+ size_t n_z{};
+ size_t n_p{};
float r=0;
- float d_max;
- float z_max;
- float h2;
- float freq;
- float c0;
- std::vector<std::complex<float>> pressures;
- std::complex<float> faj, fbj;
+ float d_max{};
+ float z_max{};
+ float h2{};
+ std::complex<float> faj{}, fbj{};
std::vector<std::complex<float>> fapbj; // faj + fbj
std::vector<std::complex<float>> fbmaj; // fbj - faj
std::vector<std::complex<float>> _u_;
std::vector<std::complex<float>> _l_;
std::vector<std::complex<float>> _d_;
- std::vector<float> n_next;
- std::complex<float> (*fun)(float, float) source_fun;
- void _step_r_(int ip);
-
+ std::vector<float> n_next{};
+// void _step_r_(size_t ip);
-public:
sound_column(float dr, size_t nz, size_t np, float dm, float zm, float f, float speed);
+ sound_column();
~sound_column();
- bool step_r(float (*fun)(float, float));
- void init_source(std::complex<float> (*fun)(float, float));
- void init_source();
+// void init_source(std::complex<float> (*fun)(float, float));
+ float freq{};
+ float c0{};
+ std::vector<std::complex<float>> pressures;
};