98 this->initialized =
true;
101 this->N = s.cells * s.points_dst;
108 n =
new ptrdiff_t[dim];
109 for(
int i = 0; i < dim - 1; i++)
111 n[dim - 1] = N / 2 + 1;
114 ptrdiff_t local_elements = 0;
115 alloc_local = fftw_mpi_local_size(dim, n, comm, &local_elements, &local_start);
116 local_end = local_start + local_elements;
119 int * global_elements =
new int[size];
120 MPI_Allgather(&local_elements, 1, MPI_INTEGER, global_elements, 1, MPI_INTEGER, comm);
123 _indices_proc_rows =
new int[N];
125 for(
int i = 0, c = 0; i < size; i++)
126 for(
int j = 0; j < global_elements[i]; j++, c++)
127 _indices_proc_rows[c] = i;
130 delete[] global_elements;
137 u_real = fftw_alloc_real(2 * alloc_local * dim);
139 this->bsize = 2 * alloc_local;
142 for(
int i = 0; i < 2 * alloc_local * dim; i++)
146 v_real = u_real + 2 * alloc_local;
149 u_comp = fftw_alloc_complex(alloc_local);
150 v_comp = fftw_alloc_complex(alloc_local);
155 w_real = v_real + 2 * alloc_local;
156 w_comp = fftw_alloc_complex(alloc_local);
160 this->e =
new double[N];
161 this->E =
new double[N];
162 this->k =
new double[N];
163 this->K =
new double[N];
164 this->c =
new double[N];
165 this->C =
new double[N];
205 fftw_plan pu = fftw_mpi_plan_dft_r2c(dim, n, u_real, u_comp, comm, FFTW_ESTIMATE);
207 fftw_destroy_plan(pu);
210 fftw_plan pv = fftw_mpi_plan_dft_r2c(dim, n, v_real, v_comp, comm, FFTW_ESTIMATE);
212 fftw_destroy_plan(pv);
217 fftw_plan pw = fftw_mpi_plan_dft_r2c(dim, n, w_real, w_comp, comm, FFTW_ESTIMATE);
219 fftw_destroy_plan(pw);
260 e_spectral += u_comp3(MIN(i, N - i), j, k_)[0] * u_comp3(MIN(i, N - i), j, k_)[0] +
261 u_comp3(MIN(i, N - i), j, k_)[1] * u_comp3(MIN(i, N - i), j, k_)[1] +
262 v_comp3(MIN(i, N - i), j, k_)[0] * v_comp3(MIN(i, N - i), j, k_)[0] +
263 v_comp3(MIN(i, N - i), j, k_)[1] * v_comp3(MIN(i, N - i), j, k_)[1] +
264 w_comp3(MIN(i, N - i), j, k_)[0] * w_comp3(MIN(i, N - i), j, k_)[0] +
265 w_comp3(MIN(i, N - i), j, k_)[1] * w_comp3(MIN(i, N - i), j, k_)[1];
280 MPI_Reduce(&e_physical, &this->e_d, 1, MPI_DOUBLE, MPI_SUM, 0, comm);
281 MPI_Reduce(&e_spectral, &this->e_s, 1, MPI_DOUBLE, MPI_SUM, 0, comm);
290 double scaling = pow(N, dim);
293 for(
int i = 0; i < N; i++)
303 for(
int j = local_start; j < local_end; j++)
305 for(
int i = 0; i < N; i++)
308 double r = sqrt(pow(MIN(i, N - i), 2.0) + pow(MIN(j, N - j), 2.0));
310 int p =
static_cast<int>(std::round(r));
312 e[p] += u_comp2(MIN(i, N - i), j)[0] * u_comp2(MIN(i, N - i), j)[0] +
313 u_comp2(MIN(i, N - i), j)[1] * u_comp2(MIN(i, N - i), j)[1] +
314 v_comp2(MIN(i, N - i), j)[0] * v_comp2(MIN(i, N - i), j)[0] +
315 v_comp2(MIN(i, N - i), j)[1] * v_comp2(MIN(i, N - i), j)[1];
325 for(
int k_ = local_start; k_ < local_end; k_++)
327 for(
int j = 0; j < N; j++)
329 for(
int i = 0; i < N; i++)
333 sqrt(pow(MIN(i, N - i), 2.0) + pow(MIN(j, N - j), 2.0) + pow(MIN(k_, N - k_), 2.0));
335 int p =
static_cast<int>(std::round(r));
337 e[p] += u_comp3(MIN(i, N - i), j, k_)[0] * u_comp3(MIN(i, N - i), j, k_)[0] +
338 u_comp3(MIN(i, N - i), j, k_)[1] * u_comp3(MIN(i, N - i), j, k_)[1] +
339 v_comp3(MIN(i, N - i), j, k_)[0] * v_comp3(MIN(i, N - i), j, k_)[0] +
340 v_comp3(MIN(i, N - i), j, k_)[1] * v_comp3(MIN(i, N - i), j, k_)[1] +
341 w_comp3(MIN(i, N - i), j, k_)[0] * w_comp3(MIN(i, N - i), j, k_)[0] +
342 w_comp3(MIN(i, N - i), j, k_)[1] * w_comp3(MIN(i, N - i), j, k_)[1];
353 MPI_Reduce(e, E, N, MPI_DOUBLE, MPI_SUM, 0, comm);
354 MPI_Reduce(k, K, N, MPI_DOUBLE, MPI_SUM, 0, comm);
355 MPI_Reduce(c, C, N, MPI_DOUBLE, MPI_SUM, 0, comm);
358 for(
int i = 0; i < N; i++)
363 E[i] /= scaling * scaling;
395 int start = local_start;
397 int delta = 2 * (N / 2 + 1) * dealii::Utilities::pow(N, dim - 2);
398 int delta_all = 2 * (N / 2 + 1) * dealii::Utilities::pow(N, dim - 1) *
sizeof(double);
401 int dofs = (end - start) * delta;
402 MPI_Offset disp = 8 *
sizeof(int) + start * delta *
sizeof(
double);
406 MPI_Type_contiguous(dofs, MPI_DOUBLE, &stype);
407 MPI_Type_commit(&stype);
411 MPI_File_open(comm, filename, MPI_MODE_WRONLY, MPI_INFO_NULL, &fh);
413 MPI_File_set_view(fh, disp + delta_all * 0, MPI_DOUBLE, MPI_DOUBLE,
"native", MPI_INFO_NULL);
414 MPI_File_write_all(fh, u_real, dofs, MPI_DOUBLE, MPI_STATUSES_IGNORE);
416 MPI_File_set_view(fh, disp + delta_all * 1, MPI_DOUBLE, MPI_DOUBLE,
"native", MPI_INFO_NULL);
417 MPI_File_write_all(fh, v_real, dofs, MPI_DOUBLE, MPI_STATUSES_IGNORE);
421 MPI_File_set_view(fh, disp + delta_all * 2, MPI_DOUBLE, MPI_DOUBLE,
"native", MPI_INFO_NULL);
422 MPI_File_write_all(fh, w_real, dofs, MPI_DOUBLE, MPI_STATUSES_IGNORE);
426 MPI_Type_free(&stype);
437 int start = local_start;
439 int delta = 2 * (N / 2 + 1) * dealii::Utilities::pow(N, dim - 2);
440 int delta_all = 2 * (N / 2 + 1) * dealii::Utilities::pow(N, dim - 1) *
sizeof(double);
443 int dofs = (end - start) * delta;
444 MPI_Offset disp = 8 *
sizeof(int) + start * delta *
sizeof(
double);
448 MPI_Type_contiguous(dofs, MPI_DOUBLE, &stype);
449 MPI_Type_commit(&stype);
453 MPI_File_open(comm, filename, MPI_MODE_RDONLY, MPI_INFO_NULL, &fh);
455 MPI_File_set_view(fh, disp + delta_all * 0, MPI_DOUBLE, MPI_DOUBLE,
"native", MPI_INFO_NULL);
456 MPI_File_read_all(fh, u_real, dofs, MPI_DOUBLE, MPI_STATUSES_IGNORE);
458 MPI_File_set_view(fh, disp + delta_all * 1, MPI_DOUBLE, MPI_DOUBLE,
"native", MPI_INFO_NULL);
459 MPI_File_read_all(fh, v_real, dofs, MPI_DOUBLE, MPI_STATUSES_IGNORE);
463 MPI_File_set_view(fh, disp + delta_all * 2, MPI_DOUBLE, MPI_DOUBLE,
"native", MPI_INFO_NULL);
464 MPI_File_read_all(fh, w_real, dofs, MPI_DOUBLE, MPI_STATUSES_IGNORE);
468 MPI_Type_free(&stype);