version working for nvc++ -stdpar -fast

main.cpp

 #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;
+}
+
+
+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));});
 }
 
+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);});
 
-std::complex<float> gaussian_init(float depth, float k0){
-    return sqrtf(k0 / 0.6f) * expf(-powf((k0 / 0.6f) * (depth - 1'000.f), 2));
+    //  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){
-    static const float slope = -0.009f;
+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;
 }

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];
     }
-}
+}*/

sound_column.h

@@ -10,38 +10,36 @@
 #include <execution>
 #include <algorithm>
 
+#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;
 };