Peano 4
Loading...
Searching...
No Matches
SWEAdjoint.cpp
Go to the documentation of this file.
1#include "SWEAdjoint.hpp"
4
6{
7 // set ghost layer in 0th row and (l_nY_ + 1)th row
8 #ifdef SharedOMP
9 #pragma omp parallel for schedule(static)
10 #endif
11 for(int i = 0; i < l_nX_ + 2; i++)
12 {
13 h_[i * (l_nY_ + 2)] = h_[i * (l_nY_ + 2) + 1];
14 hu_[i * (l_nY_ + 2)] = hu_[i * (l_nY_ + 2) + 1];
15 hv_[i * (l_nY_ + 2)] = hv_[i * (l_nY_ + 2) + 1];
16
17 h_[i * (l_nY_+ 2) + (l_nY_ + 1)] = h_[i * (l_nY_+ 2) + l_nY_];
18 hu_[i * (l_nY_+ 2) + (l_nY_ + 1)] = hu_[i * (l_nY_+ 2) + l_nY_];
19 hv_[i * (l_nY_+ 2) + (l_nY_ + 1)] = hv_[i * (l_nY_+ 2) + l_nY_];
20 }
21
22 // set ghost layer in 0th column and (l_nX_ + 1)th column
23 #ifdef SharedOMP
24 #pragma omp parallel for schedule(static)
25 #endif
26 for(int j = 0; j < l_nY_ + 2; j++)
27 {
28 h_[j] = h_[l_nY_ + 2 + j];
29 hu_[j] = hu_[l_nY_ + 2 + j];
30 hv_[j] = hv_[l_nY_ + 2 + j];
31
32 h_[(l_nX_ + 1) * (l_nY_ + 2) + j] = h_[l_nX_ * (l_nY_ + 2) + j];
33 hu_[(l_nX_ + 1) * (l_nY_ + 2) + j] = hu_[l_nX_ * (l_nY_ + 2) + j];
34 hv_[(l_nX_ + 1) * (l_nY_ + 2) + j] = hv_[l_nX_ * (l_nY_ + 2) + j];
35 }
36}
37
39{
40 double x_cdf = applications::exahype2::swe::parser::NetCDFHelper::transformIndexSimulationToCDFRange(x, l_nX_, min_x_bathymetry_, max_x_bathymetry_);
41 double y_cdf = applications::exahype2::swe::parser::NetCDFHelper::transformIndexSimulationToCDFRange(y, l_nY_, min_y_bathymetry_, max_y_bathymetry_);
42 double x_min = std::min(area_of_interest_coord_0_, area_of_interest_coord_1_);
43 double x_max = std::max(area_of_interest_coord_0_, area_of_interest_coord_1_);
44 double y_min = std::min(area_of_interest_coord_2_, area_of_interest_coord_3_);
45 double y_max = std::max(area_of_interest_coord_2_, area_of_interest_coord_3_);
46 return tarch::la::smaller(x_min, x_cdf) && tarch::la::greater(x_max, x_cdf) && tarch::la::smaller(y_min, y_cdf) && tarch::la::greater(y_max, y_cdf);
47}
48
50{
51 double x_cdf = applications::exahype2::swe::parser::NetCDFHelper::transformIndexSimulationToCDFRange(x, l_nX_, min_x_bathymetry_, max_x_bathymetry_);
52 double y_cdf = applications::exahype2::swe::parser::NetCDFHelper::transformIndexSimulationToCDFRange(y, l_nY_, min_y_bathymetry_, max_y_bathymetry_);
53 double x_center = area_of_interest_coord_0_;
54 double y_center = area_of_interest_coord_1_;
55 double radius = area_of_interest_coord_2_;
56 double distance_x = x_cdf - x_center;
57 double distance_y = y_cdf - y_center;
58 return tarch::la::smaller(sqrt((distance_x * distance_x) + (distance_y * distance_y)), radius);
59}
60
62{
63 auto dimensional_splitting = DimensionalSplitting(l_nX_, l_nY_, cell_size_, cell_size_, h_, hu_, hv_, b_);
64 double t = 0.0;
65 size_t checkpoint = 0;
66 netcdf_writer_->writeTimeStep(h_, hu_, hv_, t, checkpoint);
67 checkpoint++;
68 while(t < end_time_)
69 {
70 std::cout << "Simulating time: " << t << std::endl;
71 t += dimensional_splitting.compute_numerical_fluxes();
72 set_ghost_layer();
73 if(checkpoints_ == nullptr || tarch::la::smallerEquals(checkpoints_[checkpoint], t))
74 {
75 std::cout << "Writing time step to file" << std::endl;
76 netcdf_writer_->writeTimeStep(h_, hu_, hv_, t, checkpoint);
77 checkpoint++;
78 }
79 }
80}
81
83 const char *output_file_name,
84 const char *bathymetry_file_path,
85 double width_x, double width_y,
86 double cell_size,
87 double end_time,
88 double area_of_interest_coord_0, double area_of_interest_coord_1, double area_of_interest_coord_2, double area_of_interest_coord_3,
89 double start_time,
90 size_t checkpoints,
91 size_t flush,
92 const char* nc_key_x, const char* nc_key_y, const char* nc_key_bathymetry
93 ):
94 width_x_(width_x),
95 width_y_(width_y),
96 cell_size_(cell_size),
97 end_time_(end_time),
98 area_of_interest_coord_0_(area_of_interest_coord_0),
99 area_of_interest_coord_1_(area_of_interest_coord_1),
100 area_of_interest_coord_2_(area_of_interest_coord_2),
101 area_of_interest_coord_3_(area_of_interest_coord_3),
102 start_time_(start_time),
103 number_of_checkpoints_(checkpoints)
104{
105 // calculate number of cells for each dimension
106 l_nX_ = floor(width_x_ / cell_size_);
107 l_nY_ = floor(width_y_ / cell_size_);
108
109 topology_parser_ = new applications::exahype2::swe::parser::TopologyParser(bathymetry_file_path, l_nX_, l_nY_, nc_key_x, nc_key_y, nc_key_bathymetry);
110
115
116 // create arrays to hold values of simulation
117 h_ = new double[(l_nX_ + 2) * (l_nY_ + 2)];
118 hu_ = new double[(l_nX_ + 2) * (l_nY_ + 2)](); // initialize to 0
119 hv_ = new double[(l_nX_ + 2) * (l_nY_ + 2)](); // initialize to 0
120 b_ = new double[(l_nX_ + 2) * (l_nY_ + 2)];
121
122 // pick the correct hump criterion
124
125 // initialize bathymetry column wise
126 // initialize water height using bathymetry (and hump criterion in area of interest)
127 #ifdef SharedOMP
128 #pragma omp parallel for collapse(2) schedule(static)
129 #endif
130 for(int i = 1; i < l_nX_ + 1; i++)
131 {
132 for(int j = 1; j < l_nY_ + 1; j++)
133 {
134 double x = i;
135 double y = j;
136 b_[i * (l_nY_ + 2) + j] = topology_parser_->sample_bathymetry(x, y);
137 h_[i * (l_nY_ + 2) + j] = std::max((this->*hump_criterion)(x,y) * 1.0 - b_[i * (l_nY_ + 2) + j], 0.0);
138 }
139 }
140
141 // delete topology parser, no longer needed
142 delete topology_parser_;
143
144 // set ghost layers for bathymetry
145 // for 0th column and (l_nX_ + 1)th column
146 #ifdef SharedOMP
147 #pragma omp parallel for schedule(static)
148 #endif
149 for(int j = 0; j < l_nY_ + 2; j++)
150 {
151 b_[j] = b_[l_nY_ + 2 + j];
152 b_[(l_nX_ + 1) * (l_nY_ + 2) + j] = b_[l_nX_ * (l_nY_ + 2) + j];
153 }
154
155 // set ghost layer for 0th row and (l_nY_ + 1)th row
156 #ifdef SharedOMP
157 #pragma omp parallel for schedule(static)
158 #endif
159 for(int i = 0; i < l_nX_ + 2; i++)
160 {
161 b_[i * (l_nY_ + 2)] = b_[i * (l_nY_ + 2) + 1];
162 b_[i * (l_nY_ + 2) + (l_nY_ + 1)] = b_[i * (l_nY_ + 2) + l_nY_];
163 }
164
165 // set ghost layer for all other variables
167
168 // calculate checkpoints, if not every timestep results in output
170 {
171 checkpoints_ = new double[number_of_checkpoints_ + 1];
172 double time_width = end_time_ / number_of_checkpoints_;
173 for(int i = 0; i <= number_of_checkpoints_; i++)
174 checkpoints_[i] = time_width * i;
175 }
176
178 output_file_name,
179 bathymetry_file_path,
180 l_nX_, l_nY_,
181 cell_size_,
182 b_,
184 end_time, start_time_,
185 flush
186 );
187
188 std::cout << "Starting simulation with the following parameters:" << std::endl
189 << "width: [" << width_x_ << ", " << width_y_ << "]" << std::endl
190 << "dimension: [" << l_nX_ << ", " << l_nY_ << "]" << std::endl
191 << "cell_size: " << cell_size_ << std::endl
192 << "end_time: " << end_time_ << std::endl
193 << "hump_criterion: " << (hump_criterion == &applications::exahype2::swe::adjoint::SWEAdjoint::hump_criterion_circle ? "circle" : "square") << std::endl
194 << "output_file_name: " << output_file_name << std::endl
195 << "checkpoints: " << (number_of_checkpoints_ > 0 ? "" + number_of_checkpoints_ : "inf") << std::endl;
196 this->run_simulation();
198 #ifdef ADJOINT_ONLY
199 exit(0);
200 #endif
201}
202
204{
205 delete h_;
206 delete hu_;
207 delete hv_;
208 delete b_;
209 delete netcdf_writer_;
210 delete checkpoints_;
211}
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
void set_ghost_layer()
Set the ghost layer of all simulation arrays.
Definition SWEAdjoint.cpp:5
SWEAdjoint(const char *output_file_name, const char *bathymetry_file_path, double width_x, double width_y, double cell_size, double end_time, double area_of_interest_coord_0, double area_of_interest_coord_1, double area_of_interest_coord_2, double area_of_interest_coord_3=INFINITY, double start_time=0, size_t checkpoints=0, size_t flush=0, const char *nc_key_x="x", const char *nc_key_y="y", const char *nc_key_bathymetry="z")
Construct a new swe adjoint object.
bool hump_criterion_circle(double x, double y)
Define circular area of interest.
applications::exahype2::swe::parser::TopologyParser * topology_parser_
bool hump_criterion_square(double x, double y)
Define square area of interest.
applications::exahype2::swe::adjoint::NetCDFWriter * netcdf_writer_
double transformIndexSimulationToCDFRange(double i_z, double i_nZ, double i_min_z_cdf, double i_max_z_cdf)
Transforms coordinate from simulation plane to netCDF plane.
static bool smaller(double lhs, double rhs, double tolerance=NUMERICAL_ZERO_DIFFERENCE) InlineMethod
Smaller operator for floating point values.
bool greater(double lhs, double rhs, double tolerance=NUMERICAL_ZERO_DIFFERENCE)
bool smallerEquals(double lhs, double rhs, double tolerance=NUMERICAL_ZERO_DIFFERENCE)