v0.13.1
acoustic.cpp
Go to the documentation of this file.
1#include <cmath>
2#include <cstdlib>
3#include <iostream>
4#include <vector>
5using namespace std;
6
7#include "../../FTensor.hpp"
8
9void initial(double P1[], double P2[], double P3[], double c[], const int N)
10{
11 int cr = N / 2 - 1;
12 int cc = 7.0 * N / 8.0 - 1;
13 float s2 = 64.0 * 9.0 / pow(N / 2.0, 2.0);
14
15 for(int i = 0; i < N; ++i)
16 for(int j = 0; j < N; ++j)
17 {
18 c[i + N * j] = 0.2;
19 P1[i + N * j] = 0.0;
20 P2[i + N * j] = exp(-(pow(i - cr, 2.0) + pow(j - cc, 2.0)) * s2);
21 P3[i + N * j] = 0.0;
22 }
23
24 const int blockLeft = 0;
25 const int blockRight = 2 * N / 5.0;
26 const int blockTop = N / 3.0;
27 const int blockBottom = 2 * N / 3.0;
28
29 for(int i = blockTop; i < blockBottom; ++i)
30 for(int j = blockLeft; j < blockRight; ++j)
31 c[i + N * j] = 0.5;
32
33 int channelLeft = 4 * N / 5.0;
34 int channelRight = N;
35 int channel1Height = 3 * N / 8.0;
36 int channel2Height = 5 * N / 8.0;
37 for(int i = channelLeft; i < channelRight; ++i)
38 for(int j = channel1Height; j < channel2Height; ++j)
39 c[i + N * j] = 0;
40}
41
42void FTen(double P1_[], double P2_[], double P3_[], double c_[], const int N)
43{
44 FTensor::Tensor0<double *> P1(P1_), P2(P2_), P3(P3_), c(c_);
45
50 for(int j = 0; j < N; ++j)
51 for(int i = 0; i < N; ++i)
52 {
53 if(i != 0 && i != N - 1 && j != 0 && j != N - 1)
54 {
55 // P3 = 2*P2 - P1;
56 P3 = 2 * P2 - P1
57 + c
58 * (dd(P2, l, m, d_ijk, d_xyz)(0, 0)
59 + dd(P2, l, m, d_ijk, d_xyz)(1, 1));
60 }
61 ++P1;
62 ++P2;
63 ++P3;
64 ++c;
65 }
66 // cout << P3_[(N/2-1)+N*((7*N)/8-1)] << endl;
67}
68
69void simple(double P1[], double P2[], double P3[], double c[], const int N)
70{
71 for(int j = 1; j < N - 1; ++j)
72 for(int i = 1; i < N - 1; ++i)
73 {
74 // P3[i+N*j] = 2*P2[i+N*j] - P1[i+N*j];
75 P3[i + N * j] = (2 - 4 * c[i + N * j]) * P2[i + N * j] - P1[i + N * j]
76 + c[i + N * j]
77 * (P2[i - 1 + N * j] + P2[i + 1 + N * j]
78 + P2[i + N * (j - 1)] + P2[i + N * (j + 1)]);
79 }
80 // cout << P3[(N/2-1)+N*((7*N)/8-1)] << endl;
81}
82
83int main(int argc, char *argv[])
84{
85 const int N = atoi(argv[1]);
86
87 vector<double> P1(N * N), P2(N * N), P3(N * N), c(N * N);
88
89 initial(&c[0], &P1[0], &P2[0], &P3[0], N);
90
91 const int niters = atoi(argv[2]);
92
93 for(int i = 0; i < niters; ++i)
94 {
95 // simple(&P1[0],&P2[0],&P3[0],&c[0],N);
96 // simple(&P2[0],&P3[0],&P1[0],&c[0],N);
97 // simple(&P3[0],&P1[0],&P2[0],&c[0],N);
98 FTen(&P1[0], &P2[0], &P3[0], &c[0], N);
99 FTen(&P2[0], &P3[0], &P1[0], &c[0], N);
100 FTen(&P3[0], &P1[0], &P2[0], &c[0], N);
101 }
102}
int main(int argc, char *argv[])
Definition: acoustic.cpp:83
void initial(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:9
void FTen(double P1_[], double P2_[], double P3_[], double c_[], const int N)
Definition: acoustic.cpp:42
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
STL namespace.
const int N
Definition: speed_test.cpp:3