v0.10.0
Functions
acoustic.cpp File Reference
#include <cmath>
#include <cstdlib>
#include <iostream>
#include <vector>
#include "../../FTensor.hpp"

Go to the source code of this file.

Functions

void initial (double P1[], double P2[], double P3[], double c[], const int N)
 
void FTen (double P1_[], double P2_[], double P3_[], double c_[], const int N)
 
void simple (double P1[], double P2[], double P3[], double c[], const int N)
 
int main (int argc, char *argv[])
 

Function Documentation

◆ FTen()

void FTen ( double  P1_[],
double  P2_[],
double  P3_[],
double  c_[],
const int  N 
)

Definition at line 42 of file acoustic.cpp.

43 {
44  FTensor::Tensor0<double *> P1(P1_), P2(P2_), P3(P3_), c(c_);
45 
46  FTensor::Tensor1<int, 2> d_ijk(1, N);
47  FTensor::Tensor1<double, 2> d_xyz(1, 1);
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 }
static Index< 'l', 3 > l
static Index< 'm', 3 > m
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
static Index< 'i', 3 > i
static Index< 'j', 3 > j
const int N
Definition: speed_test.cpp:3

◆ initial()

void initial ( double  P1[],
double  P2[],
double  P3[],
double  c[],
const int  N 
)

Definition at line 9 of file acoustic.cpp.

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 }
static Index< 'i', 3 > i
static Index< 'j', 3 > j
const int N
Definition: speed_test.cpp:3

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 83 of file acoustic.cpp.

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 }
void FTen(double P1_[], double P2_[], double P3_[], double c_[], const int N)
Definition: acoustic.cpp:42
void initial(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:9
static Index< 'i', 3 > i
const int N
Definition: speed_test.cpp:3

◆ simple()

void simple ( double  P1[],
double  P2[],
double  P3[],
double  c[],
const int  N 
)
Examples
hello_world.cpp, lesson3_poisson.cpp, lesson6_radiation.cpp, lesson7_plastic.cpp, and lesson8_contact.cpp.

Definition at line 69 of file acoustic.cpp.

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 }
static Index< 'i', 3 > i
static Index< 'j', 3 > j
const int N
Definition: speed_test.cpp:3