Discontinuous Galerkin Library
#include "dg/algorithm.h"
mpi_derivatives.h
Go to the documentation of this file.
1#pragma once
2
5#include "functions.h"
6#include "derivatives.h"
7#include "mpi_grid.h"
8
9namespace dg{
10
11namespace create{
12
14namespace detail{
15
16
28template<class real_type>
29CooSparseBlockMat<real_type> save_outer_values(EllSparseBlockMat<real_type>& in, const NNCH<real_type>& c)
30{
31 //search outer values in m
32 CooSparseBlockMat<real_type> out( in.num_rows, 6, in.n, in.left_size, in.right_size);
33 int index = in.data.size()/ in.n/in.n;
34 thrust::host_vector<real_type> data_element(in.n*in.n, 0), zero(data_element);
35 bool found=false;
36 for( int i=0; i<in.num_rows; i++)
37 for( int d=0; d<in.blocks_per_line; d++)
38 {
39 if( in.cols_idx[i*in.blocks_per_line+d]==-1)
40 { //change the whole line
41 for( int k=0; k<in.blocks_per_line; k++)
42 {
43 for( int j=0; j<in.n*in.n; j++)
44 data_element[j] = in.data[ in.data_idx[i*in.blocks_per_line+k]*in.n*in.n + j];
45 int col = c.map_index( in.cols_idx[i*in.blocks_per_line+k]);
46 out.add_value( i, col, data_element);
47 in.data_idx[i*in.blocks_per_line+k] = index; //
48 in.cols_idx[i*in.blocks_per_line+k] = 0;
49 }
50 found=true;
51 }
52 if( in.cols_idx[i*in.blocks_per_line+d]==in.num_cols)
53 {
54 for( int k=0; k<in.blocks_per_line; k++)
55 {
56 for( int j=0; j<in.n*in.n; j++)
57 data_element[j] = in.data[ in.data_idx[i*in.blocks_per_line+k]*in.n*in.n + j];
58 //assume col is either 3,4,or 5
59 int col = c.map_index( in.cols_idx[i*in.blocks_per_line+k]);
60 out.add_value( i, col, data_element);
61 in.data_idx[i*in.blocks_per_line+k] = index;
62 in.cols_idx[i*in.blocks_per_line+k] = in.num_cols-1;
63 }
64 found=true;
65 }
66 }
67 if(found)
68 {
69 in.data.insert( in.data.end(), zero.begin(), zero.end());
70 }
71
72 //std::cout << "coo num entries "<<out.num_entries<<"\n";
73 //for( int i=0; i<out.num_entries; i++)
74 // std::cout << "coo entries "<<out.cols_idx[i]<<" "<<out.data_idx[i]<<"\n";
75
76 return out;
77}
78
88template<class real_type>
89EllSparseBlockMat<real_type> distribute_rows( const EllSparseBlockMat<real_type>& src, int coord, const int* howmany)
90{
91 if( howmany[1] == 1)
92 {
93 EllSparseBlockMat<real_type> temp(src);
94 temp.set_left_size( temp.left_size/howmany[0]);
95 temp.set_right_size( temp.right_size/howmany[2]);
96 return temp;
97 }
98 assert( src.num_rows == src.num_cols);
99 int chunk_size = src.num_rows/howmany[1];
100 EllSparseBlockMat<real_type> temp(chunk_size, chunk_size, src.blocks_per_line, src.data.size()/(src.n*src.n), src.n);
101 temp.set_left_size( src.left_size/howmany[0]);
102 temp.set_right_size( src.right_size/howmany[2]);
103 //first copy data elements (even though not all might be needed it doesn't slow down things either)
104 for( unsigned i=0; i<src.data.size(); i++)
105 temp.data[i] = src.data[i];
106 //now grab the right chunk of cols and data indices
107 for( unsigned i=0; i<temp.cols_idx.size(); i++)
108 {
109 temp.data_idx[i] = src.data_idx[ coord*(chunk_size*src.blocks_per_line)+i];
110 temp.cols_idx[i] = src.cols_idx[ coord*(chunk_size*src.blocks_per_line)+i];
111 //data indices are correct but cols are still the global indices (remapping a bit clumsy)
112 //first in the zeroth line the col idx might be (global)num_cols - 1 -> map that to -1
113 if( coord==0 && i/src.blocks_per_line == 0 && temp.cols_idx[i] == src.num_cols-1) temp.cols_idx[i] = -1;
114 //second in the last line the col idx mighty be 0 -> map to (global)num_cols
115 if( coord==(howmany[1]-1)&& (int)i/src.blocks_per_line == temp.num_rows-1 && temp.cols_idx[i] == 0) temp.cols_idx[i] = src.num_cols;
116 //Elements are now in the range -1, 0, 1,..., (global)num_cols
117 //now shift this range to chunk range -1,..,chunk_size
118 temp.cols_idx[i] = (temp.cols_idx[i] - coord*chunk_size );
119 }
120 return temp;
121}
122
123
124} //namespace detail
125
129
139template<class real_type>
141{
142 EllSparseBlockMat<real_type> matrix = dg::create::dx( g.global(), bcx, dir);
143 unsigned vector_dimensions[] = {(unsigned)(g.nx()*g.local().Nx()), (unsigned)(g.ny()*g.local().Ny()), 1}; //x, y, z
144 MPI_Comm comm = g.communicator();
145 int ndims;
146 MPI_Cartdim_get( comm, &ndims);
147 assert( ndims == 2);
148 int dims[ndims], periods[ndims], coords[ndims];
149 MPI_Cart_get( comm, ndims, dims, periods, coords);
150
151 int howmany[] = {dims[1], dims[0], 1}; //left, middle, right
152 //distribute_rows, collective and save_outer_values are aware of howmany[1] == 1
153 EllSparseBlockMat<real_type> inner = detail::distribute_rows(matrix, coords[0], howmany);
154 NNCH<real_type> c( g.nx(), vector_dimensions, comm, 0);
155 CooSparseBlockMat<real_type> outer = detail::save_outer_values(inner,c);
156
158}
159
169template<class real_type>
171{
172 EllSparseBlockMat<real_type> matrix = dg::create::dy( g.global(), bcy, dir);
173 unsigned vector_dimensions[] = {(unsigned)(g.nx()*g.local().Nx()), (unsigned)(g.ny()*g.local().Ny()), 1}; //x, y, z
174 MPI_Comm comm = g.communicator();
175 int ndims;
176 MPI_Cartdim_get( comm, &ndims);
177 assert( ndims == 2);
178 int dims[ndims], periods[ndims], coords[ndims];
179 MPI_Cart_get( comm, ndims, dims, periods, coords);
180
181 int howmany[] = {1, dims[1], dims[0]};
182 EllSparseBlockMat<real_type> inner = detail::distribute_rows(matrix, coords[1], howmany);
183 NNCH<real_type> c( g.ny(), vector_dimensions, comm, 1);
184 CooSparseBlockMat<real_type> outer = detail::save_outer_values(inner,c);
185
187}
188
197template<class real_type>
199{
201 unsigned vector_dimensions[] = {(unsigned)(g.nx()*g.local().Nx()), (unsigned)(g.ny()*g.local().Ny()), 1}; //x, y, z
202 MPI_Comm comm = g.communicator();
203 int ndims;
204 MPI_Cartdim_get( comm, &ndims);
205 assert( ndims == 2);
206 int dims[ndims], periods[ndims], coords[ndims];
207 MPI_Cart_get( comm, ndims, dims, periods, coords);
208
209 int howmany[] = {dims[1], dims[0], 1}; //left, middle, right
210 EllSparseBlockMat<real_type> inner = detail::distribute_rows(matrix, coords[0], howmany);
211 NNCH<real_type> c( g.nx(), vector_dimensions, comm, 0);
212 CooSparseBlockMat<real_type> outer = detail::save_outer_values(inner,c);
213
215}
224template<class real_type>
226{
228 unsigned vector_dimensions[] = {(unsigned)(g.nx()*g.local().Nx()), (unsigned)(g.ny()*g.local().Ny()), 1}; //x, y, z
229 MPI_Comm comm = g.communicator();
230 int ndims;
231 MPI_Cartdim_get( comm, &ndims);
232 assert( ndims == 2);
233 int dims[ndims], periods[ndims], coords[ndims];
234 MPI_Cart_get( comm, ndims, dims, periods, coords);
235
236 int howmany[] = {1, dims[1], dims[0]};
237 EllSparseBlockMat<real_type> inner = detail::distribute_rows(matrix, coords[1], howmany);
238 NNCH<real_type> c( g.ny(), vector_dimensions, comm, 1);
239 CooSparseBlockMat<real_type> outer = detail::save_outer_values(inner,c);
240
242}
243
253template<class real_type>
255{
256 EllSparseBlockMat<real_type> matrix = dg::create::dx( g.global(), bcx, dir);
257 unsigned vector_dimensions[] = {(unsigned)(g.nx()*g.local().Nx()), (unsigned)(g.ny()*g.local().Ny()), (unsigned)(g.nz()*g.local().Nz())}; //x, y, z
258 MPI_Comm comm = g.communicator();
259 int ndims;
260 MPI_Cartdim_get( comm, &ndims);
261 assert( ndims == 3);
262 int dims[ndims], periods[ndims], coords[ndims];
263 MPI_Cart_get( comm, ndims, dims, periods, coords);
264
265 int howmany[] = {dims[2]*dims[1], dims[0], 1};
266 EllSparseBlockMat<real_type> inner = detail::distribute_rows(matrix, coords[0], howmany);
267 NNCH<real_type> c( g.nx(), vector_dimensions, comm, 0);
268 CooSparseBlockMat<real_type> outer = detail::save_outer_values(inner,c);
269
271}
281template<class real_type>
283{
284 EllSparseBlockMat<real_type> matrix = dg::create::dy( g.global(), bcy, dir);
285 unsigned vector_dimensions[] = {(unsigned)(g.nx()*g.local().Nx()), (unsigned)(g.ny()*g.local().Ny()), (unsigned)(g.nz()*g.local().Nz())}; //x, y, z
286 MPI_Comm comm = g.communicator();
287 int ndims;
288 MPI_Cartdim_get( comm, &ndims);
289 assert( ndims == 3);
290 int dims[ndims], periods[ndims], coords[ndims];
291 MPI_Cart_get( comm, ndims, dims, periods, coords);
292
293 int howmany[] = {dims[2], dims[1], dims[0]};
294 EllSparseBlockMat<real_type> inner = detail::distribute_rows(matrix, coords[1], howmany);
295 NNCH<real_type> c( g.ny(), vector_dimensions, comm, 1);
296 CooSparseBlockMat<real_type> outer = detail::save_outer_values(inner,c);
297
299}
309template<class real_type>
311{
312 EllSparseBlockMat<real_type> matrix = dg::create::dz( g.global(), bcz, dir);
313 unsigned vector_dimensions[] = {(unsigned)(g.nx()*g.local().Nx()), (unsigned)(g.ny()*g.local().Ny()), (unsigned)(g.nz()*g.local().Nz())}; //x, y, z
314 MPI_Comm comm = g.communicator();
315 int ndims;
316 MPI_Cartdim_get( comm, &ndims);
317 assert( ndims == 3);
318 int dims[ndims], periods[ndims], coords[ndims];
319 MPI_Cart_get( comm, ndims, dims, periods, coords);
320
321 int howmany[] = {1, dims[2], dims[1]*dims[0]};
322 EllSparseBlockMat<real_type> inner = detail::distribute_rows(matrix, coords[2], howmany);
323 NNCH<real_type> c( g.nz(), vector_dimensions, comm, 2);
324 CooSparseBlockMat<real_type> outer = detail::save_outer_values(inner,c);
325
327}
328
337template<class real_type>
339{
341 unsigned vector_dimensions[] = {(unsigned)(g.nx()*g.local().Nx()), (unsigned)(g.ny()*g.local().Ny()), (unsigned)(g.nz()*g.local().Nz())}; //x, y, z
342 MPI_Comm comm = g.communicator();
343 int ndims;
344 MPI_Cartdim_get( comm, &ndims);
345 assert( ndims == 3);
346 int dims[ndims], periods[ndims], coords[ndims];
347 MPI_Cart_get( comm, ndims, dims, periods, coords);
348
349 int howmany[] = {dims[2]*dims[1], dims[0], 1};
350 EllSparseBlockMat<real_type> inner = detail::distribute_rows(matrix, coords[0], howmany);
351 NNCH<real_type> c( g.nx(), vector_dimensions, comm, 0);
352 CooSparseBlockMat<real_type> outer = detail::save_outer_values(inner,c);
353
355}
356
365template<class real_type>
367{
369 unsigned vector_dimensions[] = {(unsigned)(g.nx()*g.local().Nx()), (unsigned)(g.ny()*g.local().Ny()), (unsigned)(g.nz()*g.local().Nz())}; //x, y, z
370 MPI_Comm comm = g.communicator();
371 int ndims;
372 MPI_Cartdim_get( comm, &ndims);
373 assert( ndims == 3);
374 int dims[ndims], periods[ndims], coords[ndims];
375 MPI_Cart_get( comm, ndims, dims, periods, coords);
376
377 int howmany[] = {dims[2], dims[1], dims[0]};
378 EllSparseBlockMat<real_type> inner = detail::distribute_rows(matrix, coords[1], howmany);
379 NNCH<real_type> c( g.ny(), vector_dimensions, comm, 1);
380 CooSparseBlockMat<real_type> outer = detail::save_outer_values(inner,c);
381
383}
392template<class real_type>
394{
396 unsigned vector_dimensions[] = {(unsigned)(g.nx()*g.local().Nx()), (unsigned)(g.ny()*g.local().Ny()), (unsigned)(g.nz()*g.local().Nz())}; //x, y, z
397 MPI_Comm comm = g.communicator();
398 int ndims;
399 MPI_Cartdim_get( comm, &ndims);
400 assert( ndims == 3);
401 int dims[ndims], periods[ndims], coords[ndims];
402 MPI_Cart_get( comm, ndims, dims, periods, coords);
403
404 int howmany[] = {1, dims[2], dims[1]*dims[0]};
405 EllSparseBlockMat<real_type> inner = detail::distribute_rows(matrix, coords[2], howmany);
406 NNCH<real_type> c( g.nz(), vector_dimensions, comm, 2);
407 CooSparseBlockMat<real_type> outer = detail::save_outer_values(inner,c);
408
410}
411
420template<class real_type>
422{
423 return dx( g, g.bcx(), dir);
424}
425
434template<class real_type>
436{
437 return dx( g, g.bcx(), dir);
438}
446template<class real_type>
448{
449 return jumpX( g, g.bcx());
450}
451
459template<class real_type>
461{
462 return jumpX( g, g.bcx());
463}
464
473template<class real_type>
475{
476 return dy( g, g.bcy(), dir);
477}
478
487template<class real_type>
489{
490 return dy( g, g.bcy(), dir);
491}
492
500template<class real_type>
502{
503 return jumpY( g, g.bcy());
504}
505
513template<class real_type>
515{
516 return jumpY( g, g.bcy());
517}
518
527template<class real_type>
529{
530 return dz( g, g.bcz(), dir);
531}
532
540template<class real_type>
542{
543 return jumpZ( g, g.bcz());
544}
545
546
548
549
550} //namespace create
551} //namespace dg
Convenience functions to create 2D derivatives.
Some utility functions for the dg::evaluate routines.
static DG_DEVICE double zero(double x)
Definition: functions.h:29
EllSparseBlockMat< real_type > dz(const aRealTopology3d< real_type > &g, bc bcz, direction dir=centered)
Create 3d derivative in z-direction.
Definition: derivatives.h:308
EllSparseBlockMat< real_type > dx(const aRealTopology2d< real_type > &g, bc bcx, direction dir=centered)
Create 2d derivative in x-direction.
Definition: derivatives.h:33
EllSparseBlockMat< real_type > jumpZ(const aRealTopology3d< real_type > &g, bc bcz)
Matrix that contains jump terms in Z direction in 3D.
Definition: derivatives.h:190
bc
Switch between boundary conditions.
Definition: enums.h:15
EllSparseBlockMat< real_type > jumpX(const aRealTopology2d< real_type > &g, bc bcx)
Matrix that contains 2d jump terms in X direction.
Definition: derivatives.h:95
EllSparseBlockMat< real_type > dy(const aRealTopology2d< real_type > &g, bc bcy, direction dir=centered)
Create 2d derivative in y-direction.
Definition: derivatives.h:65
direction
Direction of a discrete derivative.
Definition: enums.h:97
EllSparseBlockMat< real_type > jumpY(const aRealTopology2d< real_type > &g, bc bcy)
Matrix that contains 2d jump terms in Y direction.
Definition: derivatives.h:112
@ centered
centered derivative (cell to the left and right and current cell)
Definition: enums.h:100
MPI Grid objects.
MPI matrix classes.
This is the namespace for all functions and classes defined and used by the discontinuous Galerkin li...
Coo Sparse Block Matrix format.
Definition: sparseblockmat.h:179
Ell Sparse Block Matrix format.
Definition: sparseblockmat.h:46
Communicator for asynchronous nearest neighbor communication.
Definition: mpi_vector.h:181
Distributed memory matrix class, asynchronous communication.
Definition: mpi_matrix.h:65
2D MPI abstract grid class
Definition: mpi_grid.h:45
const RealGrid2d< real_type > & local() const
Return a non-MPI grid local for the calling process.
Definition: mpi_grid.h:252
bc bcy() const
global y boundary
Definition: mpi_grid.h:132
unsigned nx() const
number of polynomial coefficients in x
Definition: mpi_grid.h:106
unsigned ny() const
number of polynomial coefficients in y
Definition: mpi_grid.h:108
MPI_Comm communicator() const
Return mpi cartesian communicator that is used in this grid.
Definition: mpi_grid.h:138
const RealGrid2d< real_type > & global() const
Return the global non-MPI grid.
Definition: mpi_grid.h:262
bc bcx() const
global x boundary
Definition: mpi_grid.h:126
3D MPI Grid class
Definition: mpi_grid.h:329
bc bcz() const
global z boundary
Definition: mpi_grid.h:454
unsigned nz() const
number of polynomial coefficients in z
Definition: mpi_grid.h:418
MPI_Comm communicator() const
Return mpi cartesian communicator that is used in this grid.
Definition: mpi_grid.h:459
bc bcy() const
global y boundary
Definition: mpi_grid.h:448
const RealGrid3d< real_type > & global() const
Return the global non-MPI grid.
Definition: mpi_grid.h:562
bc bcx() const
global x boundary
Definition: mpi_grid.h:442
unsigned ny() const
number of polynomial coefficients in y
Definition: mpi_grid.h:416
const RealGrid3d< real_type > & local() const
Return a non-MPI grid local for the calling process.
Definition: mpi_grid.h:560
unsigned nx() const
number of polynomial coefficients in x
Definition: mpi_grid.h:414