WaveSimPP 1.0
A C library for solving the wave equation and reconstructing the wave field
•All Classes Namespaces Functions Modules Pages
Namespaces | Functions
WaveSimCore

Namespaces

namespace  waveSimCore
 

Functions

boost::multi_array< double, 2 > waveSimCore::get_sigma_1 (boost::multi_array< double, 1 > x, double dx, int nx, int nz, double c_max, int n=15, double R=1e-3, double m=2.0)
 
boost::multi_array< double, 2 > waveSimCore::get_sigma_2 (boost::multi_array< double, 1 > z, double dz, int nx, int nz, double c_max, int n=10, double R=1e-3, double m=2.0)
 
boost::multi_array< double, 2 > waveSimCore::get_profile (double xmin, double xmax, double zmin, double zmax, int nx, int nz, double r)
 Get the velocity profile of the model as a 2D Array.
 
boost::multi_array< double, 2 > waveSimCore::dfdx (boost::multi_array< double, 2 > f, double dx)
 Takes the partial derivative of a 2D matrix f with respect to x.
 
boost::multi_array< double, 2 > waveSimCore::dfdz (boost::multi_array< double, 2 > f, double dz)
 Takes the partial derivative of a 2D matrix f with respect to z.
 
boost::multi_array< double, 2 > waveSimCore::d2fdx2 (boost::multi_array< double, 2 > f, double dx)
 Takes the second partial derivative of a 2D matrix f with respect to x.
 
boost::multi_array< double, 2 > waveSimCore::d2fdz2 (boost::multi_array< double, 2 > f, double dz)
 Takes the second partial derivative of a 2D matrix f with respect to z.
 
boost::multi_array< double, 2 > waveSimCore::divergence (boost::multi_array< double, 2 > f1, boost::multi_array< double, 2 > f2, double dx, double dz)
 
boost::multi_array< double, 3 > waveSimCore::wave_solver (boost::multi_array< double, 2 > c, double dt, double dx, double dz, int nt, int nx, int nz, boost::multi_array< double, 3 > f)
 
boost::multi_array< double, 3 > waveSimCore::wave_solver_complex (boost::multi_array< double, 2 > c, double dt, double dx, double dz, int nt, int nx, int nz, boost::multi_array< double, 3 > f, boost::multi_array< double, 2 > sigma_1, boost::multi_array< double, 2 > sigma_2)
 
boost::multi_array< double, 3 > waveSimCore::ricker (int i_s, int j_s, double f, double amp, double shift, double tmin, double tmax, int nt, int nx, int nz)
 Get the Ricker wavelet as a 3D Array.
 

Detailed Description

Function Documentation

◆ d2fdx2()

boost::multi_array< double, 2 > waveSimCore::d2fdx2 ( boost::multi_array< double, 2 >  f,
double  dx 
)

Takes the second partial derivative of a 2D matrix f with respect to x.

Definition at line 32 of file helper_func.hpp.

33 {
34 boost::multi_array<double, 2> df = dfdx(f, dx);
35 boost::multi_array<double, 2> df2 = dfdx(df, dx);
36 return df2;
37 }
waveSimCore::dfdx
boost::multi_array< double, 2 > dfdx(boost::multi_array< double, 2 > f, double dx)
Takes the partial derivative of a 2D matrix f with respect to x.
Definition: helper_func.hpp:18

◆ d2fdz2()

boost::multi_array< double, 2 > waveSimCore::d2fdz2 ( boost::multi_array< double, 2 >  f,
double  dz 
)

Takes the second partial derivative of a 2D matrix f with respect to z.

Definition at line 40 of file helper_func.hpp.

41 {
42 boost::multi_array<double, 2> df = dfdz(f, dz);
43 boost::multi_array<double, 2> df2 = dfdz(df, dz);
44 return df2;
45 }
waveSimCore::dfdz
boost::multi_array< double, 2 > dfdz(boost::multi_array< double, 2 > f, double dz)
Takes the partial derivative of a 2D matrix f with respect to z.
Definition: helper_func.hpp:25

◆ dfdx()

boost::multi_array< double, 2 > waveSimCore::dfdx ( boost::multi_array< double, 2 >  f,
double  dx 
)

Takes the partial derivative of a 2D matrix f with respect to x.

Definition at line 18 of file helper_func.hpp.

19 {
20 std::vector<boost::multi_array<double, 2>> grad_f = np::gradient(f, {dx, dx});
21 return grad_f[0];
22 }
np::gradient
::value constexpr std::vector< boost::multi_array< T, ND > > gradient(boost::multi_array< T, ND > inArray, std::initializer_list< T > args)
Definition: np.hpp:90

◆ dfdz()

boost::multi_array< double, 2 > waveSimCore::dfdz ( boost::multi_array< double, 2 >  f,
double  dz 
)

Takes the partial derivative of a 2D matrix f with respect to z.

Definition at line 25 of file helper_func.hpp.

26 {
27 std::vector<boost::multi_array<double, 2>> grad_f = np::gradient(f, {dz, dz});
28 return grad_f[1];
29 }

◆ divergence()

boost::multi_array< double, 2 > waveSimCore::divergence ( boost::multi_array< double, 2 >  f1,
boost::multi_array< double, 2 >  f2,
double  dx,
double  dz 
)

Takes the divergence of a 2D matrices fx,fz with respect to x and z
Returns dfx/dx + dfz/dz

Definition at line 49 of file helper_func.hpp.

51 {
52 boost::multi_array<double, 2> f_x = dfdx(f1, dx);
53 boost::multi_array<double, 2> f_z = dfdz(f2, dz);
54 boost::multi_array<double, 2> div = f_x + f_z;
55 return div;
56 }

◆ get_profile()

boost::multi_array< double, 2 > waveSimCore::get_profile ( double  xmin,
double  xmax,
double  zmin,
double  zmax,
int  nx,
int  nz,
double  r 
)

Get the velocity profile of the model as a 2D Array.

Definition at line 16 of file computational.hpp.

17 {
18 boost::multi_array<double, 2> c(boost::extents[nx][nz]);
19
20 boost::multi_array<double, 1> x = np::linspace(xmin, xmax, nx);
21 boost::multi_array<double, 1> z = np::linspace(zmin, zmax, nz);
22
23 const boost::multi_array<double, 1> axis[2] = {x, z};
24 std::vector<boost::multi_array<double, 2>> XZ = np::meshgrid(axis, false, np::ij);
25
26 double x_0 = xmax / 2.0;
27 double z_0 = zmax / 2.0;
28
29 for (int i = 0; i < nx; i++)
30 {
31 for (int j = 0; j < nz; j++)
32 {
33 if (np::pow(XZ[0][i][j] - x_0, 2.0) + np::pow(XZ[1][i][j] - z_0, 2.0) <= np::pow(r, 2.0))
34 c[i][j] = 3.0;
35 else
36 c[i][j] = 3.0;
37 }
38 }
39
40 return c;
41 }
np::pow
::value constexpr boost::multi_array< T, ND > pow(const boost::multi_array< T, ND > &input_array, const T exponent)
Implements the numpy pow function on multi arrays.
Definition: np.hpp:319
np::meshgrid
::value constexpr std::vector< boost::multi_array< T, ND > > meshgrid(const boost::multi_array< T, 1 >(&cinput)[ND], bool sparsing=false, indexing indexing_type=xy)
Definition: np.hpp:184
np::linspace
boost::multi_array< double, 1 > linspace(double start, double stop, long unsigned int num)
Implements the numpy linspace function.
Definition: np.hpp:161

◆ get_sigma_1()

boost::multi_array< double, 2 > waveSimCore::get_sigma_1 ( boost::multi_array< double, 1 >  x,
double  dx,
int  nx,
int  nz,
double  c_max,
int  n = 15,
double  R = 1e-3,
double  m = 2.0 
)

Returns the sigma 1 coefficient for the PML boundary conditions
sigma_1 and sigma_2 are the coefficients for the x and z directions respectively
Check Johnson, Steven G. (2021).

Definition at line 21 of file coeff.hpp.

23 {
24 boost::multi_array<double, 2> sigma_1(boost::extents[nx][nz]);
25 const double PML_width = n * dx;
26
27 const double sigma_max = -c_max * log(R) * (m + 1.0) / np::pow(PML_width, m + 1.0);
28
29 const double x_0 = np::max(x) - PML_width;
30
31 boost::multi_array<double, 1> polynomial(boost::extents[nx]);
32
33 for (int i = 0; i < nx; i++)
34 {
35 if (x[i] > x_0)
36 {
37 polynomial[i] = sigma_max * np::pow(np::abs(x[i] - x_0), m);
38 polynomial[nx - 1 - i] = polynomial[i];
39 }
40 else
41 {
42 polynomial[i] = 0;
43 }
44 }
45
46 for (int i = 0; i < nx; i++)
47 for (int j = 0; j < nz; j++)
48 sigma_1[i][j] = polynomial[i];
49
50 return sigma_1;
51 }
np::abs
::value constexpr T abs(T input)
Implements the numpy abs function for a scalar.
Definition: np.hpp:463
np::max
::value constexpr T max(boost::multi_array< T, ND > const &input_array)
Implements the numpy max function for an n-dimensionl multi array.
Definition: np.hpp:384

◆ get_sigma_2()

boost::multi_array< double, 2 > waveSimCore::get_sigma_2 ( boost::multi_array< double, 1 >  z,
double  dz,
int  nx,
int  nz,
double  c_max,
int  n = 10,
double  R = 1e-3,
double  m = 2.0 
)

Returns the sigma 2 coefficient for the PML boundary conditions
sigma_1 and sigma_2 are the coefficients for the x and z directions respectively
Check Johnson, Steven G. (2021).

Definition at line 56 of file coeff.hpp.

58 {
59 boost::multi_array<double, 2> sigma_2(boost::extents[nx][nz]);
60 const double PML_width = n * dz;
61 const double sigma_max = -c_max * log(R) * (m + 1.0) / np::pow(PML_width, m + 1.0);
62
63 const double z_0 = np::max(z) - PML_width;
64
65 boost::multi_array<double, 1> polynomial(boost::extents[nz]);
66 for (int j = 0; j < nz; j++)
67 {
68 if (z[j] > z_0)
69 {
70 polynomial[j] = sigma_max * np::pow(np::abs(z[j] - z_0), m);
71 polynomial[nz - 1 - j] = polynomial[j];
72 }
73 else
74 {
75 polynomial[j] = 0;
76 }
77 }
78
79 for (int i = 0; i < nx; i++)
80 for (int j = 0; j < nz; j++)
81 sigma_2[i][j] = polynomial[j];
82
83 return sigma_2;
84 }

◆ ricker()

boost::multi_array< double, 3 > waveSimCore::ricker ( int  i_s,
int  j_s,
double  f,
double  amp,
double  shift,
double  tmin,
double  tmax,
int  nt,
int  nx,
int  nz 
)

Get the Ricker wavelet as a 3D Array.

Definition at line 14 of file source.hpp.

16 {
17 const double pi = 3.141592654;
18
19 boost::multi_array<double, 1> t = np::linspace(tmin, tmax, nt);
20 boost::multi_array<double, 1> pft2 = np::pow(pi * f * (t - shift), 2.0);
21 boost::multi_array<double, 1> r = amp * (1.0 - 2.0 * pft2) * np::exp(-1.0 * pft2);
22
23 int dimensions_x[] = {nx};
24 boost::multi_array<double, 1> x = np::zeros<double>(dimensions_x);
25
26 int dimensions_z[] = {nz};
27 boost::multi_array<double, 1> z = np::zeros<double>(dimensions_z);
28
29 x[i_s] = 1.0;
30 z[j_s] = 1.0;
31
32 const boost::multi_array<double, 1> axis[3] = {r, x, z};
33 std::vector<boost::multi_array<double, 3>> RXZ = np::meshgrid(axis, false, np::ij);
34 boost::multi_array<double, 3> source = RXZ[0] * RXZ[1] * RXZ[2];
35
36 return source;
37 }
np::exp
::value constexpr boost::multi_array< T, ND > exp(const boost::multi_array< T, ND > &input_array)
Implements the numpy exp function on multi arrays.
Definition: np.hpp:289

◆ wave_solver()

boost::multi_array< double, 3 > waveSimCore::wave_solver ( boost::multi_array< double, 2 >  c,
double  dt,
double  dx,
double  dz,
int  nt,
int  nx,
int  nz,
boost::multi_array< double, 3 >  f 
)

This function solves the wave equation using the methods explained in the readme
The boundary conditions are so that waves bounce at the boundary

Definition at line 27 of file solver.hpp.

30 {
31
32 const boost::multi_array<double, 2> CX = np::pow(c * dt / dx, 2.0);
33 const boost::multi_array<double, 2> CZ = np::pow(c * dt / dz, 2.0);
34 const double Cf = np::pow(dt, 2.0);
35
36 int dimensions_1[] = {nt, nx, nz};
37 boost::multi_array<double, 3> u = np::zeros<double>(dimensions_1);
38 int dimensions_4[] = {nx, nz};
39 boost::multi_array<double, 2> u_xx = np::zeros<double>(dimensions_4);
40 int dimensions_5[] = {nx, nz};
41 boost::multi_array<double, 2> u_zz = np::zeros<double>(dimensions_5);
42
43 for (int n = 1; n < nt - 1; n++)
44 {
45 for (int i = 1; i < nx - 1; i++)
46 {
47 for (int j = 1; j < nz - 1; j++)
48 {
49 u_xx[i][j] = u[n][i + 1][j] + u[n][i - 1][j] - 2.0 * u[n][i][j];
50 u_zz[i][j] = u[n][i][j + 1] + u[n][i][j - 1] - 2.0 * u[n][i][j];
51 }
52 }
53
54 // ! Update u
55 for (int i = 0; i < nx; i++)
56 {
57 for (int j = 0; j < nz; j++)
58 {
59 u[n + 1][i][j] = 2.0 * u[n][i][j] + CX[i][j] * u_xx[i][j] + CZ[i][j] * u_zz[i][j] + Cf * f[n][i][j] - u[n - 1][i][j];
60 }
61 }
62
63 // Dirichlet boundary condition
64 for (int i = 0; i < nx; i++)
65 {
66 u[n + 1][i][0] = 0.0;
67 u[n + 1][i][nz - 1] = 0.0;
68 }
69 for (int j = 0; j < nz; j++)
70 {
71 u[n + 1][0][j] = 0.0;
72 u[n + 1][nx - 1][j] = 0.0;
73 }
74 }
75 return u;
76 }

◆ wave_solver_complex()

boost::multi_array< double, 3 > waveSimCore::wave_solver_complex ( boost::multi_array< double, 2 >  c,
double  dt,
double  dx,
double  dz,
int  nt,
int  nx,
int  nz,
boost::multi_array< double, 3 >  f,
boost::multi_array< double, 2 >  sigma_1,
boost::multi_array< double, 2 >  sigma_2 
)

Solves the wave equation using a more complex solver compared to wave_solver
Also has different boundary conditions. In this one the boundary extends to infinity

Definition at line 22 of file solver_complex.hpp.

26 {
27
28 const boost::multi_array<double, 2> C1 = 1.0 + (dt * (sigma_1 + sigma_2) * 0.5);
29 const boost::multi_array<double, 2> C2 = (sigma_1 * sigma_2 * np::pow(dt, 2.0)) - 2.0;
30 const boost::multi_array<double, 2> C3 = 1.0 - (dt * (sigma_1 + sigma_2) * 0.5);
31 const boost::multi_array<double, 2> C4 = np::pow(dt * c, 2.0);
32 const boost::multi_array<double, 2> C5 = 1.0 + (dt * sigma_1 * 0.5);
33 const boost::multi_array<double, 2> C6 = 1.0 + (dt * sigma_2 * 0.5);
34 const boost::multi_array<double, 2> C7 = 1.0 - (dt * sigma_1 * 0.5);
35 const boost::multi_array<double, 2> C8 = 1.0 - (dt * sigma_2 * 0.5);
36
37 int dimensions_1[] = {nt, nx, nz};
38 boost::multi_array<double, 3> u = np::zeros<double>(dimensions_1);
39
40 int dimensions_2[] = {nx, nz};
41 boost::multi_array<double, 2> q_1 = np::zeros<double>(dimensions_2);
42 int dimensions_3[] = {nx, nz};
43 boost::multi_array<double, 2> q_2 = np::zeros<double>(dimensions_3);
44
45 int dimensions_4[] = {nx, nz};
46 boost::multi_array<double, 2> u_xx = np::zeros<double>(dimensions_4);
47 int dimensions_5[] = {nx, nz};
48 boost::multi_array<double, 2> u_zz = np::zeros<double>(dimensions_5);
49
50 boost::multi_array<double, 2> f_n(boost::extents[nx][nz]);
51 boost::multi_array<double, 2> u_n(boost::extents[nx][nz]);
52 boost::multi_array<double, 2> u_n_1(boost::extents[nx][nz]);
53
54 for (int n = 1; n < nt - 1; n++)
55 {
56
57 for (int i = 0; i < nx; i++)
58 {
59 for (int j = 0; j < nz; j++)
60 {
61 f_n[i][j] = f[n][i][j];
62 u_n[i][j] = u[n][i][j];
63 u_n_1[i][j] = u[n - 1][i][j];
64 }
65 }
66
67 boost::multi_array<double, 2> div = divergence(q_1 * sigma_1, q_2 * sigma_2, dx, dz);
68 boost::multi_array<double, 2> dq_1dx = dfdx(q_1, dx);
69 boost::multi_array<double, 2> dq_2dz = dfdz(q_2, dz);
70 u_xx = d2fdx2(u_n, dx); // (nx, nz)
71 u_zz = d2fdz2(u_n, dz); // (nx, nz)
72 boost::multi_array<double, 2> dudx = dfdx(u_n, dx);
73 boost::multi_array<double, 2> dudz = dfdz(u_n, dz);
74
75 // for (int i = 1; i < nx-1; i++)
76 // {
77 // for (int j = 1; j < nz-1; j++)
78 // {
79 // u_xx[i][j] = u_n[i+1][j] + u_n[i-1][j] - 2.0 * u_n[i][j];
80 // u_zz[i][j] = u_n[i][j+1] + u_n[i][j-1] - 2.0 * u_n[i][j];
81 // }
82 // }
83
84 // ! Update u
85 for (int i = 0; i < nx; i++)
86 {
87 for (int j = 0; j < nz; j++)
88 {
89 u[n + 1][i][j] = (C4[i][j] * ((u_xx[i][j] / np::pow(dx, 2.0)) + (u_zz[i][j] / np::pow(dz, 2.0)) - div[i][j] + sigma_2[i][j] * dq_1dx[i][j] + sigma_1[i][j] * dq_2dz[i][j] + f_n[i][j]) - (C2[i][j] * u_n[i][j]) - (C3[i][j] * u_n_1[i][j])) / C1[i][j];
90 }
91 }
92
93 // ! Update q_1, q_2
94 for (int i = 0; i < nx; i++)
95 {
96 for (int j = 0; j < nz; j++)
97 {
98 q_1[i][j] = (dt * dudx[i][j] + C7[i][j] * q_1[i][j]) / C5[i][j];
99 q_2[i][j] = (dt * dudz[i][j] + C8[i][j] * q_2[i][j]) / C6[i][j];
100 }
101 }
102 //
103 // Dirichlet boundary condition
104 for (int i = 0; i < nx; i++)
105 {
106 u[n + 1][i][0] = 0.0;
107 u[n + 1][i][nz - 1] = 0.0;
108 }
109 for (int j = 0; j < nz; j++)
110 {
111 u[n + 1][0][j] = 0.0;
112 u[n + 1][nx - 1][j] = 0.0;
113 }
114 }
115 return u;
116 }
waveSimCore::d2fdx2
boost::multi_array< double, 2 > d2fdx2(boost::multi_array< double, 2 > f, double dx)
Takes the second partial derivative of a 2D matrix f with respect to x.
Definition: helper_func.hpp:32
waveSimCore::d2fdz2
boost::multi_array< double, 2 > d2fdz2(boost::multi_array< double, 2 > f, double dz)
Takes the second partial derivative of a 2D matrix f with respect to z.
Definition: helper_func.hpp:40
waveSimCore::divergence
boost::multi_array< double, 2 > divergence(boost::multi_array< double, 2 > f1, boost::multi_array< double, 2 > f2, double dx, double dz)
Definition: helper_func.hpp:49
Generated by doxygen 1.9.6