Peano 4
Loading...
Searching...
No Matches
DimensionalSplitting.cpp
Go to the documentation of this file.
4
5applications::exahype2::swe::adjoint::DimensionalSplitting::DimensionalSplitting(size_t l_nX, size_t l_nY, double l_dx, double l_dy, double *h, double *hu, double *hv, double *b) :
6l_nX_(l_nX),
7l_nY_(l_nY),
8l_dx_(l_dx),
9l_dy_(l_dy)
10{
11 this->h_ = h;
12 this->hu_ = hu;
13 this->hv_ = hv;
14 this->b_ = b;
15
16 this->h_net_updates_above_ = new double[(l_nX_ + 2)*(l_nY_ * 2)]();
17 this->h_net_updates_below_ = new double[(l_nX_ + 2)*(l_nY_ * 2)]();
18 this->h_net_updates_left_ = new double[(l_nX_ + 2)*(l_nY_ * 2)]();
19 this->h_net_updates_right_ = new double[(l_nX_ + 2)*(l_nY_ * 2)]();
20
21 this->hu_net_updates_left_ = new double[(l_nX_ + 2)*(l_nY_ * 2)]();
22 this->hu_net_updates_right_ = new double[(l_nX_ + 2)*(l_nY_ * 2)]();
23 this->hv_net_updates_above_ = new double[(l_nX_ + 2)*(l_nY_ * 2)]();
24 this->hv_net_updates_below_ = new double[(l_nX_ + 2)*(l_nY_ * 2)]();
25}
26
28{
29 delete this->h_net_updates_above_;
30 delete this->h_net_updates_below_;
31 delete this->h_net_updates_left_;
32 delete this->h_net_updates_right_;
33
34 delete this->hu_net_updates_left_;
35 delete this->hu_net_updates_right_;
36 delete this->hv_net_updates_above_;
37 delete this->hv_net_updates_below_;
38}
39
41{
42 double max_wave_speed = 0.0;
43
44 double max_x_wave_speed = 0.0;
45 // for all columns
46 #ifdef SharedOMP
47 #pragma omp parallel for collapse(2) reduction(max: max_x_wave_speed) schedule(runtime)
48 #endif
49 for (int i = 1; i < l_nX_ + 2; i++)
50 {
51 // for all rows
52 for (int j = 0; j < l_nY_ + 2; j++)
53 {
54 double max_edge_speed = 0.0;
55
57 h_[(i - 1) * (l_nY_ + 2) + j], h_[i * (l_nY_ + 2) + j],
58 hu_[(i - 1) * (l_nY_ + 2) + j], hu_[i * (l_nY_ + 2) + j],
59 b_[(i - 1) * (l_nY_ + 2) + j], b_[i * (l_nY_ + 2) + j],
60 h_net_updates_left_[(i - 1) * (l_nY_ + 2) + j], h_net_updates_right_[(i - 1) * (l_nY_ + 2) + j],
61 hu_net_updates_left_[(i - 1) * (l_nY_ + 2) + j], hu_net_updates_right_[(i - 1) * (l_nY_ + 2) + j],
62 max_edge_speed);
63
64 //update the maximum wave speed of the x-sweep
65 max_x_wave_speed = std::max(max_x_wave_speed, max_edge_speed);
66 }
67 }
68
69 // Steady states in x dimension?
70 if (tarch::la::greater(max_x_wave_speed, 0.0))
71 update_unknowns_x_sweep(0.4 * l_dx_ / max_x_wave_speed);
72
73 double max_y_wave_speed = 0.0;
74 // for all rows
75 #ifdef SharedOMP
76 #pragma omp parallel for collapse(2) reduction(max: max_y_wave_speed) schedule(runtime)
77 #endif
78 for (int j = 1; j < l_nY_ + 2; j++)
79 {
80 // for all columns except most left and right ones (ghost layer)
81 for (int i = 1; i < l_nX_ + 1; i++)
82 {
83 double max_edge_speed = 0.0;
84
86 h_[i * (l_nY_ + 2) + (j - 1)], h_[i * (l_nY_ + 2) + j],
87 hv_[i * (l_nY_ + 2) + (j - 1)], hv_[i * (l_nY_ + 2) + j],
88 b_[i * (l_nY_ + 2) + (j - 1)], b_[i * (l_nY_ + 2) + j],
89 h_net_updates_below_[i * (l_nY_ + 2) + (j - 1)], h_net_updates_above_[i * (l_nY_ + 2) + (j - 1)],
90 hv_net_updates_below_[i * (l_nY_ + 2) + (j - 1)], hv_net_updates_above_[i * (l_nY_ + 2) + (j - 1)],
91 max_edge_speed);
92
93 //update the maximum wave speed of the x-sweep
94 max_y_wave_speed = std::max(max_y_wave_speed, max_edge_speed);
95 }
96 }
97
98 max_wave_speed = std::max(max_x_wave_speed, max_y_wave_speed);
99 if(tarch::la::greater(max_wave_speed, 0.0))
100 {
101 double max_time_step = 0.4 * l_dx_ / max_wave_speed;
102 update_unknowns_y_sweep(max_time_step);
103 return max_time_step;
104 }
105 return max_wave_speed;
106}
107
109{
110 // update cells x-sweep
111 const double factordx = dt / l_dx_;
112 // for all inner columns
113 #ifdef SharedOMP
114 #pragma omp parallel for collapse(2) schedule(static)
115 #endif
116 for (int i = 1; i < l_nX_ + 1; i++)
117 {
118 // for all rows
119 for (int j = 0; j < l_nY_ + 2; j++)
120 {
121 // -dx/dt * (Update of left cell to the right + Update of right cell to the left)
122 h_[i * (l_nY_ + 2) + j] -= factordx * (h_net_updates_right_[(i - 1) * (l_nY_ + 2) + j] + h_net_updates_left_[i * (l_nY_ + 2) + j]);
123 hu_[i * (l_nY_ + 2) + j] -= factordx * (hu_net_updates_right_[(i-1) * (l_nY_ + 2) + j] + hu_net_updates_left_[i * (l_nY_ + 2) + j]);
124 }
125 }
126}
127
129{
130 // update cells x-sweep
131 const double factor_dy = dt / l_dy_;
132 // for all inner columns
133 #ifdef SharedOMP
134 #pragma omp parallel for collapse(2) schedule(static)
135 #endif
136 for (int j = 1; j < l_nY_ + 1; j++)
137 {
138 // for all inner columns
139 for (int i = 1; i < l_nY_ + 1; i++)
140 {
141 // -dy/dt * (Update of cell above to below + Update of cell below to above)
142 h_[i * (l_nY_ + 2) + j] -= factor_dy * (h_net_updates_above_[i * (l_nY_ + 2) + (j - 1)] + h_net_updates_below_[i * (l_nY_ + 2) + j]);
143 hv_[i * (l_nY_ + 2) + j] -= factor_dy * (hv_net_updates_above_[i * (l_nY_ + 2) + (j - 1)] + hv_net_updates_below_[i * (l_nY_ + 2) + j]);
144 }
145 }
146}
And from this we can write down f$ nabla phi_i nabla phi_i dx but since we are constructing matrix let s investigate the f$ j
static void compute_net_updates(double m_h_l, double m_h_r, double m_hu_l, double m_hu_r, double bat_l, double bat_r, double &h_netup_l, double &h_netup_r, double &hu_netup_l, double &hu_netup_r, double &maxEdgeSpeed)
double compute_numerical_fluxes()
computes numerical fluxes for one time step of the simulation.
void update_unknowns_y_sweep(double dt)
updates unknowns after y-sweep.
void update_unknowns_x_sweep(double dt)
updates unknowns after x-sweep.
DimensionalSplitting(size_t l_nX, size_t l_nY, double l_dx, double l_dy, double *h, double *hu, double *hv, double *b)
bool greater(double lhs, double rhs, double tolerance=NUMERICAL_ZERO_DIFFERENCE)