Skip to content
Snippets Groups Projects
Commit a8725da9 authored by ferrari's avatar ferrari
Browse files

initial commit

parents
No related branches found
No related tags found
No related merge requests found
cmake_minimum_required(VERSION 3.17)
project(Parabolic__)
set(CMAKE_CXX_STANDARD 17)
add_executable(Parabolic__ main.cpp sound_column.cpp sound_column.h)
\ No newline at end of file
main.cpp 0 → 100644
#include <iostream>
#include <vector>
#include <cstring>
#include "sound_column.h"
#include <execution>
#define DEFAULT_DISTANCE 20'000.0f
#define DEFAULT_DEPTH 5'000.0f
#define DEFAULT_MIN_FREQ 25.0f
#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);
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;
int i=1;
while (i < argc){
if (strcmp(argv[i], "--d_max") != 0 || strcmp(argv[i], "-d") != 0)
r_max = atof(argv[++i]);
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)
f_min = atof(argv[++i]);
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){
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;
return 1;
}
else {
std::cout << "Unrecognized argument " << argv[i] << std::endl;
}
i++;
}
return run_simulation(r_max, z_max, f_min, f_max);
}
std::complex<float> gaussian_init(float depth, float k0){
return sqrtf(k0 / 0.6f) * expf(-powf((k0 / 0.6f) * (depth - 1'000.f), 2));
}
float media(float r, float depth){
static const float slope = -0.009f;
if ( depth > slope * r + 2531.f)
return WATER_SPEED / 1650.f;
else if ( depth > slope/2 * r + 3042.f)
return WATER_SPEED / 1700.f;
else if ( depth > slope/3 * r + 4111.f)
return WATER_SPEED / 1800.f;
else
return WATER_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;
column.reserve(nf);
std::ranges::iota_view indexes((size_t) 0, nf);
std::for_each(std::execution::par_unseq, indexes.begin(), indexes.end(),
[&](size_t i) {column[i] = 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);
column[i]->init_source(gaussian_init);});
std::for_each(std::execution::par_unseq, indexes.begin(), indexes.end(),
[&](size_t i) {while (column[i]->step_r(media));});
return 0;
}
\ No newline at end of file
//
// Created by ferrari on 07/12/2020.
//
#include "sound_column.h"
sound_column::sound_column(float dr, size_t nz, size_t np, float dm, float zm, float f, float speed) {
delta_r = dr;
d_max=dm;
z_max=zm;
h2 = zm*zm/(nz*nz);
n_z=nz;
n_p=np;
freq=f;
c0=speed;
pressures.reserve(nz);
_u_.reserve(nz-1);
_l_.reserve(nz-1);
_d_.reserve(nz);
n_next.reserve(nz);
fapbj.reserve(np);
fbmaj.reserve(np);
for(int 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);
fapbj[i] = faj + fbj;
fbmaj[i] = fbj - faj;
}
}
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, z_ind.begin(), z_ind.end(),
[&](size_t i) {n_next[i] = (*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, z_ind.begin(), z_ind.end(),
[&](size_t i) {pressures[i] *= 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, z_ind.begin(), z_ind.end(),
[&](size_t i) {pressures[i] = (*fun)(z_max/n_z*i, 2*M_1_PIf32 * freq / c0);});
}
void sound_column::_step_r_(int 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);
std::for_each(std::execution::par_unseq, z_ind.begin(), z_ind.end(),
[&](size_t i) {_d_[i] = 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));});
// first matmul
// Inv tri matmul
_u_[0] /= _d_[0];
pressures[0] /= _d_[0];
for(size_t i=1; i < n_z-1; i++){
_u_[i] /= _d_[i] - _l_[i-1] * _u_[i-1];
pressures[i] = (pressures[i] - _l_[i-1] * pressures[i-1]) / (_d_[i] - _l_[i-1] * _u_[i-1]);
}
pressures[n_z-1] = (pressures[n_z-1] - _l_[n_z-2] * pressures[n_z-2]) / (_d_[n_z-1] - _l_[n_z-2] * _u_[n_z-2]);
for(ssize_t i=n_z-2; i >= 0; i--) {
pressures[i] -= _u_[i]*pressures[i+1];
}
}
//
// Created by ferrari on 07/12/2020.
//
#ifndef PARABOLIC___SOUND_COLUMN_H
#define PARABOLIC___SOUND_COLUMN_H
#include <vector>
#include <complex>
#include <ranges>
#include <execution>
# define M_1_PIf32 (0.318309886183790671537767526745028724)
class sound_column {
private:
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;
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;
void _step_r_(int ip);
public:
sound_column(float dr, size_t nz, size_t np, float dm, float zm, float f, float speed);
~sound_column();
bool step_r(float (*fun)(float, float));
void init_source(std::complex<float> (*fun)(float, float));
};
#endif //PARABOLIC___SOUND_COLUMN_H
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment