Peano 4
Loading...
Searching...
No Matches
AdjointParser.cpp
Go to the documentation of this file.
1#include "AdjointParser.h"
2#include <iostream>
4#include "../Constants.h"
5
7
8applications::exahype2::swe::parser::AdjointParser::AdjointParser(const char *adjoint_file_path, double start_time){
9 this->parse_adjoint_file(adjoint_file_path, start_time);
10
11}
12
14 delete this->adjoint_h_;
15 delete this->adjoint_hu_;
16 delete this->adjoint_hv_;
17 delete this->adjoint_b_;
18 delete this->adjoint_timestamps_;
19}
20
21void applications::exahype2::swe::parser::AdjointParser::interpolate_values(size_t i_index, double i_x, double i_y, double *interpolated_values){
22 double x_cdf = NetCDFHelper::transformIndexSimulationToCDFRange(i_x, DomainSize(0), min_x_adjoint_, max_x_adjoint_);
23 double y_cdf = NetCDFHelper::transformIndexSimulationToCDFRange(i_y, DomainSize(1), min_y_adjoint_, max_y_adjoint_);
24 if(tarch::la::smaller(x_cdf, min_x_adjoint_) || tarch::la::greater(x_cdf, max_x_adjoint_) || tarch::la::smaller(y_cdf, min_y_adjoint_) || tarch::la::greater(y_cdf, max_y_adjoint_))
25 return;
26 int array_index = NetCDFHelper::transformIndexCDFRangeToArray(x_cdf, y_cdf, min_x_adjoint_, min_y_adjoint_, max_x_adjoint_, max_y_adjoint_, size_x_adjoint_, size_y_adjoint_);
27 int array_index_shift = (number_timestamps_ - i_index - 1) * array_index; // reverese time as adjoint simulation has to be backwards in time but is computed forwards in time
28 interpolated_values[0] = adjoint_h_[array_index_shift]; // get h of adjoint
29 interpolated_values[1] = adjoint_hu_[array_index_shift]; // get hu of adjoint
30 interpolated_values[2] = adjoint_hv_[array_index_shift]; // get hv of adjoint
31 interpolated_values[3] = adjoint_b_[array_index]; // get bathymetry of adjoint
32 return;
33}
34
35double applications::exahype2::swe::parser::AdjointParser::calculate_inner_product(size_t i_index, double i_x, double i_y, const double* Q){
36 if(tarch::la::equals(Q[0], 0.0))
37 {
38 std::cout << "h = 0!" << std::endl;
39 return 0.0;
40 }
41
42 double interpolated_values[4]{0};
43 interpolate_values(i_index, i_x, i_y, interpolated_values);
44 // interpolated_values[0] is h
45 // interpolated_values[1] is hu
46 // interpolated_values[2] is hv
47 // interpolated_values[3] is bathymetry
48
49 if(tarch::la::equals(interpolated_values[0], 0.0))
50 {
51 std::cout << "h_p = 0!" << std::endl;
52 std::cout << interpolated_values[0] << std::endl;
53 return 0.0;
54 }
55
56 return std::abs((Q[0] + Q[3]) * (interpolated_values[0] + interpolated_values[3])
57 + Q[1] * interpolated_values[1]
58 + Q[2] * interpolated_values[2]);
59}
60
61double applications::exahype2::swe::parser::AdjointParser::max_inner_product(double i_t, double i_x, double i_y, const double* Q){
62
63 if(!tarch::la::equals(i_t, t_last_call_))
64 {
65 t_last_call_ = i_t;
66 set_start_index(i_t);
67 }
68
69 double t_nm = adjoint_timestamps_[inner_product_start_index_-1];
70
71 double t_current;
72 double inner_product_current;
73 double t_comp = i_t + (t_f_ - t_s_);
74 double max_inner_product = calculate_inner_product(inner_product_start_index_-1, i_x, i_y, Q);
75
76 for(size_t i = inner_product_start_index_; i < number_timestamps_; i++)
77 {
78 t_current = adjoint_timestamps_[i];
79 if(tarch::la::smaller(i_t, t_current) && tarch::la::smallerEquals(t_nm, t_comp))
80 {
81 inner_product_current = calculate_inner_product(i, i_x, i_y, Q);
82 max_inner_product = std::max(inner_product_current, max_inner_product);
83 t_nm = t_current;
84 }
85 else
86 break;
87 }
88 return max_inner_product;
89}
90
92 double t_nm = adjoint_timestamps_[inner_product_start_index_ - 1];
93 for(size_t i = inner_product_start_index_; i < number_timestamps_; i++)
94 {
95 if(tarch::la::smallerEquals(t_nm, i_t + (t_f_ - t_s_)) && tarch::la::greater(adjoint_timestamps_[i], i_t))
96 {
97 inner_product_start_index_ = i;
98 return;
99 }
100 t_nm = adjoint_timestamps_[i];
101 }
102}
103
104void applications::exahype2::swe::parser::AdjointParser::parse_adjoint_file(const char *adjoint_file_path, double start_time){
105 int ncid;
106 // Open adjoint netCDF file
107 if ((ncid = NetCDFReader::openCdf(adjoint_file_path)) == -1)
108 {
109 std::cout << "Guided AMR could not be initialized: Adjoint file could not be opened at: "<< adjoint_file_path << std::endl;
110 return;
111 }
112
113 // Read number of checkpoints completed by adjoint simulation.
114 // There needs to be at least one record in the netCDF file!
115 if ((number_timestamps_ = NetCDFReader::getDimension(ncid, "time")) == 0)
116 {
117 std::cout << "Guided AMR could not be initialized: No timestamps written so far." << std::endl;
118 return;
119 }
120
121 adjoint_timestamps_ = new double[number_timestamps_];
122
123 // Read time of the last checkpoint written in netCDF.
124 {
125 if (NetCDFReader::readVariable1D(ncid, "time", adjoint_timestamps_) != 0)
126 {
127 std::cout << "Guided AMR could not be initialized: Timestamps not found." << std::endl;
128 return;
129 }
130 }
131
132 // Read x-dimension of adjoint data
133 if ((size_x_adjoint_ = NetCDFReader::getDimension(ncid, "x")) == 0)
134 {
135 std::cout << "Guided AMR could not be initialized: Dimension of the x-Adjoint could not be read." << std::endl;
136 return;
137 }
138
139 // Read y-dimension of adjoint data
140 if ((size_y_adjoint_ = NetCDFReader::getDimension(ncid, "y")) == 0)
141 {
142 std::cout << "Guided AMR could not be initialized: Dimension of the y-Adjoint could not be read." << std::endl;
143 return;
144 }
145
146 adjoint_h_ = new double[number_timestamps_*size_x_adjoint_*size_y_adjoint_];
147 // Read h
148 if (NetCDFReader::readVariable3DHyperslab(ncid, "h", adjoint_h_, size_x_adjoint_, size_y_adjoint_, 0, number_timestamps_) == -1)
149 {
150 std::cout << "Guided AMR could not be initialized: Adjoint h could not be read." << std::endl;
151 return;
152 }
153
154 adjoint_hu_ = new double[number_timestamps_*size_x_adjoint_*size_y_adjoint_];
155 // Read hu
156 if (NetCDFReader::readVariable3DHyperslab(ncid, "hu", adjoint_hu_, size_x_adjoint_, size_y_adjoint_, 0, number_timestamps_) == -1)
157 {
158 std::cout << "Guided AMR could not be initialized: Adjoint hu could not be read." << std::endl;
159 return;
160 }
161
162 adjoint_hv_ = new double[number_timestamps_*size_x_adjoint_*size_y_adjoint_];
163 // Read hv
164 if (NetCDFReader::readVariable3DHyperslab(ncid, "hv", adjoint_hv_, size_y_adjoint_, size_x_adjoint_, 0, number_timestamps_) == -1)
165 {
166 std::cout << "Guided AMR could not be initialized: Adjoint hv could not be read." << std::endl;
167 return;
168 }
169
170 adjoint_b_ = new double[size_x_adjoint_*size_y_adjoint_];
171 // Read b
172 if (NetCDFReader::readVariable2D(ncid, "b", adjoint_b_) == -1)
173 {
174 std::cout << "Guided AMR could not be initialized: Adjoint bathymetry could not be read." << std::endl;
175 return;
176 }
177
178 // Get max and min x
179 {
180 double tmp[size_x_adjoint_]{0};
181 if (NetCDFReader::readVariable1D(ncid, "x", tmp) == -1)
182 {
183 std::cout << "Guided AMR could not be initialized: Adjoint x indices could not be read." << std::endl;
184 return;
185 }
186 min_x_adjoint_ = tmp[0];
187 max_x_adjoint_ = tmp[size_x_adjoint_ - 1];
188 }
189
190 // Get max and min y
191 {
192 double tmp[size_y_adjoint_]{0};
193 if (NetCDFReader::readVariable1D(ncid, "y", tmp) == -1)
194 {
195 std::cout << "Guided AMR could not be initialized: Adjoint y indices could not be read." << std::endl;
196 return;
197 }
198 min_y_adjoint_ = tmp[0];
199 max_y_adjoint_ = tmp[size_y_adjoint_ - 1];
200 }
201
202 // Read end time
203 {
204 double *tmp = nullptr;
205 if((tmp = NetCDFReader::getGlobalAttDouble(ncid, "end_time")) == nullptr)
206 {
207 std::cout << "Guided AMR could not be initialized: Adjoint end_time could not be read." << std::endl;
208 return;
209 }
210 t_f_ = tmp[0];
211 }
212
213 // Read start time
214 if(isinf(start_time))
215 {
216 double *tmp = nullptr;
217 if((tmp = NetCDFReader::getGlobalAttDouble(ncid, "start_time")) == nullptr)
218 {
219 std::cout << "Guided AMR could not be initialized: Adjoint start_time could not be read." << std::endl;
220 return;
221 }
222 t_s_ = tmp[0];
223 }
224 else
225 {
226 t_s_ = start_time;
227 }
228
229 // Close adjoint netCDF file
230 if (NetCDFReader::closeCdf(ncid) == -1)
231 {
232 std::cout << "Guided AMR could not be initialized: Adjoint file could not be closed." << std::endl;
233 return;
234 }
235}
double max_inner_product(double i_t, double i_x, double i_y, const double *Q)
Calculates the maximum inner product of adjoint and current simulation.
void parse_adjoint_file(const char *adjoint_file_path, double start_time=INFINITY)
Parses given netCDF file of adjoint simulation and reads out all necessary values.
AdjointParser(const char *adjoint_file_path, double start_time=INFINITY)
double calculate_inner_product(size_t i_index, double i_x, double i_y, const double *Q)
Calculates the inner product of adjoint and current simulation at given timestamp index.
void interpolate_values(size_t i_index, double i_x, double i_y, double *interpolated_values)
Interpolates spatial indices to fit adjoint solution and returns adjoint solution values at these ind...
void set_start_index(double i_t)
Call this once every time step to set start index for inner product calculation.
int transformIndexCDFRangeToArray(double i_x, double i_y, double i_min_x, double i_min_y, double i_max_x, double i_max_y, size_t i_xDim, size_t i_yDim)
Transforms coordinate pair from netCDF plane to 1D row-wise array index.
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.
int readVariable2D(int ncid, const char *key, double **variable, size_t rows, size_t cols)
Read contents of a 2D double variable from netCDF file into 2D array.
int closeCdf(int ncid)
Closes netCDF file.
int readVariable3DHyperslab(int ncid, const char *key, double *variable, size_t rows, size_t cols, size_t start, size_t end)
Read contents of a 3D double variable from netCDF file into 1D array starting at a specified first-di...
size_t getDimension(int ncid, const char *key)
Get the dimension of a variable.
int readVariable1D(int ncid, const char *key, double *variable)
Read contents of a 1D double variable from netCDF file.
double * getGlobalAttDouble(int ncid, const char *name)
reads the global attribute name of type double and returns an array containing all the associated val...
int openCdf(String filename)
Opens a netCDF file, returns ncid.
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)
bool equals(const Matrix< Rows, Cols, Scalar > &lhs, const Matrix< Rows, Cols, Scalar > &rhs, const Scalar &tolerance=NUMERICAL_ZERO_DIFFERENCE)
Compares to matrices on equality by means of a numerical accuracy.