Discontinuous Galerkin Library
#include "dg/algorithm.h"
Loading...
Searching...
No Matches
dx.h
Go to the documentation of this file.
1#pragma once
2
3#include <cassert>
4
5#include "dg/backend/sparseblockmat.h"
6#include "grid.h"
7#include "functions.h"
8#include "operator.h"
9#include "weights.h"
10
14namespace dg
15{
16namespace create
17{
19
20
32template<class real_type>
33EllSparseBlockMat<real_type, thrust::host_vector> dx_symm(int n, int N, real_type h, bc bcx)
34{
35 constexpr int invalid_idx = EllSparseBlockMat<real_type, thrust::host_vector>::invalid_index;
36
37 SquareMatrix<real_type> l = create::lilj<real_type>(n);
38 SquareMatrix<real_type> r = create::rirj<real_type>(n);
39 SquareMatrix<real_type> lr = create::lirj<real_type>(n);
40 SquareMatrix<real_type> rl = create::rilj<real_type>(n);
41 SquareMatrix<real_type> d = create::pidxpj<real_type>(n);
42 SquareMatrix<real_type> t = create::pipj_inv<real_type>(n);
43 t *= 2./h;
44
45 SquareMatrix< real_type> a = 1./2.*t*(d-d.transpose());
46 //bcx = PER
47 SquareMatrix<real_type> a_bound_right(a), a_bound_left(a);
48 //left boundary
49 if( bcx == DIR || bcx == DIR_NEU )
50 a_bound_left += 0.5*t*l;
51 else if (bcx == NEU || bcx == NEU_DIR)
52 a_bound_left -= 0.5*t*l;
53 //right boundary
54 if( bcx == DIR || bcx == NEU_DIR)
55 a_bound_right -= 0.5*t*r;
56 else if( bcx == NEU || bcx == DIR_NEU)
57 a_bound_right += 0.5*t*r;
58 if( bcx == PER ) //periodic bc
59 a_bound_left = a_bound_right = a;
60 SquareMatrix<real_type> b = t*(1./2.*rl);
61 SquareMatrix<real_type> bp = t*(-1./2.*lr); //pitfall: T*-m^T is NOT -(T*m)^T
62 //transform to XSPACE
63 SquareMatrix<real_type> backward=dg::DLT<real_type>::backward(n);
64 SquareMatrix<real_type> forward=dg::DLT<real_type>::forward(n);
65 a = backward*a*forward, a_bound_left = backward*a_bound_left*forward;
66 b = backward*b*forward, a_bound_right = backward*a_bound_right*forward;
67 bp = backward*bp*forward;
68 //assemble the matrix
69 if( bcx != PER)
70 {
71 EllSparseBlockMat<real_type, thrust::host_vector> A(N, N, 3, 5, n);
72 for( int i=0; i<n; i++)
73 for( int j=0; j<n; j++)
74 {
75 A.data[(0*n+i)*n+j] = bp(i,j);
76 A.data[(1*n+i)*n+j] = a(i,j);
77 A.data[(2*n+i)*n+j] = b(i,j);
78 A.data[(3*n+i)*n+j] = a_bound_left(i,j);
79 A.data[(4*n+i)*n+j] = a_bound_right(i,j);
80 }
81 A.data_idx[0*3+0] = 3; //a_bound_left
82 A.cols_idx[0*3+0] = 0;
83 A.data_idx[0*3+1] = 2; //b
84 A.cols_idx[0*3+1] = 1;
85 A.data_idx[0*3+2] = 2; //
86 A.cols_idx[0*3+2] = invalid_idx; //prevent unnecessary data fetch
87 for( int i=1; i<N-1; i++)
88 for( int d=0; d<3; d++)
89 {
90 A.data_idx[i*3+d] = d; //bp, a, b
91 A.cols_idx[i*3+d] = i+d-1;
92 }
93 A.data_idx[(N-1)*3+0] = 0; //bp
94 A.cols_idx[(N-1)*3+0] = N-2;
95 A.data_idx[(N-1)*3+1] = 4; //a_bound_right
96 A.cols_idx[(N-1)*3+1] = N-1;
97 A.data_idx[(N-1)*3+2] = 4; //0
98 A.cols_idx[(N-1)*3+2] = invalid_idx; //prevent unnecessary data fetch
99 return A;
100
101 }
102 else //periodic
103 {
104 EllSparseBlockMat<real_type, thrust::host_vector> A(N, N, 3, 3, n);
105 for( int i=0; i<n; i++)
106 for( int j=0; j<n; j++)
107 {
108 A.data[(0*n+i)*n+j] = bp(i,j);
109 A.data[(1*n+i)*n+j] = a(i,j);
110 A.data[(2*n+i)*n+j] = b(i,j);
111 }
112 for( int i=0; i<N; i++)
113 for( int d=0; d<3; d++)
114 {
115 A.data_idx[i*3+d] = d; //bp, a, b
116 A.cols_idx[i*3+d] = (i+d-1+N)%N;
117 }
118 return A;
119 }
120};
121
132template<class real_type>
133EllSparseBlockMat<real_type, thrust::host_vector> dx_plus( int n, int N, real_type h, bc bcx )
134{
135 constexpr int invalid_idx = EllSparseBlockMat<real_type, thrust::host_vector>::invalid_index;
136
137 SquareMatrix<real_type> l = create::lilj<real_type>(n);
138 SquareMatrix<real_type> r = create::rirj<real_type>(n);
139 SquareMatrix<real_type> lr = create::lirj<real_type>(n);
140 SquareMatrix<real_type> rl = create::rilj<real_type>(n);
141 SquareMatrix<real_type> d = create::pidxpj<real_type>(n);
142 SquareMatrix<real_type> t = create::pipj_inv<real_type>(n);
143 t *= 2./h;
144 SquareMatrix<real_type> a = t*(-l-d.transpose());
145 //if( dir == backward) a = -a.transpose();
146 SquareMatrix<real_type> a_bound_left = a; //PER, NEU and NEU_DIR
147 SquareMatrix<real_type> a_bound_right = a; //PER, DIR and NEU_DIR
148 if( bcx == dg::DIR || bcx == dg::DIR_NEU)
149 a_bound_left = t*(-d.transpose());
150 if( bcx == dg::NEU || bcx == dg::DIR_NEU)
151 a_bound_right = t*(d);
152 SquareMatrix<real_type> b = t*rl;
153 //transform to XSPACE
154 SquareMatrix<real_type> backward=dg::DLT<real_type>::backward(n);
155 SquareMatrix<real_type> forward=dg::DLT<real_type>::forward(n);
156 a = backward*a*forward, a_bound_left = backward*a_bound_left*forward;
157 b = backward*b*forward, a_bound_right = backward*a_bound_right*forward;
158 //assemble the matrix
159 if( bcx != PER)
160 {
161 EllSparseBlockMat<real_type, thrust::host_vector> A(N, N, 2, 4, n);
162 for( int i=0; i<n; i++)
163 for( int j=0; j<n; j++)
164 {
165 A.data[(0*n+i)*n+j] = a(i,j);
166 A.data[(1*n+i)*n+j] = b(i,j);
167 A.data[(2*n+i)*n+j] = a_bound_left(i,j);
168 A.data[(3*n+i)*n+j] = a_bound_right(i,j);
169 }
170 A.data_idx[0*2+0] = 2; //a_bound_left
171 A.cols_idx[0*2+0] = 0;
172 A.data_idx[0*2+1] = 1; //b
173 A.cols_idx[0*2+1] = 1;
174 for( int i=1; i<N-1; i++) //a
175 for( int d=0; d<2; d++)
176 {
177 A.data_idx[i*2+d] = d; //a, b
178 A.cols_idx[i*2+d] = i+d; //0,1
179 }
180 A.data_idx[(N-1)*2+0] = 3; //a_bound_right
181 A.cols_idx[(N-1)*2+0] = N-1;
182 A.data_idx[(N-1)*2+1] = 3; //0
183 A.cols_idx[(N-1)*2+1] = invalid_idx; //prevent unnecessary data fetch
184 return A;
185
186 }
187 else //periodic
188 {
189 EllSparseBlockMat<real_type, thrust::host_vector> A(N, N, 2, 2, n);
190 for( int i=0; i<n; i++)
191 for( int j=0; j<n; j++)
192 {
193 A.data[(0*n+i)*n+j] = a(i,j);
194 A.data[(1*n+i)*n+j] = b(i,j);
195 }
196 for( int i=0; i<N; i++)
197 for( int d=0; d<2; d++)
198 {
199 A.data_idx[i*2+d] = d; //a, b
200 A.cols_idx[i*2+d] = (i+d+N)%N;
201 }
202 return A;
203 }
204};
205
216template<class real_type>
217EllSparseBlockMat<real_type, thrust::host_vector> dx_minus( int n, int N, real_type h, bc bcx )
218{
219 constexpr int invalid_idx = EllSparseBlockMat<real_type, thrust::host_vector>::invalid_index;
220 SquareMatrix<real_type> l = create::lilj<real_type>(n);
221 SquareMatrix<real_type> r = create::rirj<real_type>(n);
222 SquareMatrix<real_type> lr = create::lirj<real_type>(n);
223 SquareMatrix<real_type> rl = create::rilj<real_type>(n);
224 SquareMatrix<real_type> d = create::pidxpj<real_type>(n);
225 SquareMatrix<real_type> t = create::pipj_inv<real_type>(n);
226 t *= 2./h;
227 SquareMatrix<real_type> a = t*(l+d);
228 //if( dir == backward) a = -a.transpose();
229 SquareMatrix<real_type> a_bound_right = a; //PER, NEU and DIR_NEU
230 SquareMatrix<real_type> a_bound_left = a; //PER, DIR and DIR_NEU
231 if( bcx == dg::DIR || bcx == dg::NEU_DIR)
232 a_bound_right = t*(-d.transpose());
233 if( bcx == dg::NEU || bcx == dg::NEU_DIR)
234 a_bound_left = t*d;
235 SquareMatrix<real_type> bp = -t*lr;
236 //transform to XSPACE
237 SquareMatrix<real_type> backward=dg::DLT<real_type>::backward(n);
238 SquareMatrix<real_type> forward=dg::DLT<real_type>::forward(n);
239 a = backward*a*forward, a_bound_left = backward*a_bound_left*forward;
240 bp = backward*bp*forward, a_bound_right = backward*a_bound_right*forward;
241
242 //assemble the matrix
243 if(bcx != dg::PER)
244 {
245 EllSparseBlockMat<real_type, thrust::host_vector> A(N, N, 2, 4, n);
246 for( int i=0; i<n; i++)
247 for( int j=0; j<n; j++)
248 {
249 A.data[(0*n+i)*n+j] = bp(i,j);
250 A.data[(1*n+i)*n+j] = a(i,j);
251 A.data[(2*n+i)*n+j] = a_bound_left(i,j);
252 A.data[(3*n+i)*n+j] = a_bound_right(i,j);
253 }
254 A.data_idx[0*2+0] = 2; //a_bound_left
255 A.cols_idx[0*2+0] = 0;
256 A.data_idx[0*2+1] = 2; //0
257 A.cols_idx[0*2+1] = invalid_idx; //prevent data fetch
258 for( int i=1; i<N-1; i++) //a
259 for( int d=0; d<2; d++)
260 {
261 A.data_idx[i*2+d] = d; //bp, a
262 A.cols_idx[i*2+d] = i+d-1;
263 }
264 A.data_idx[(N-1)*2+0] = 0; //bp
265 A.cols_idx[(N-1)*2+0] = N-2;
266 A.data_idx[(N-1)*2+1] = 3; //a_bound_right
267 A.cols_idx[(N-1)*2+1] = N-1;
268 return A;
269
270 }
271 else //periodic
272 {
273 EllSparseBlockMat<real_type, thrust::host_vector> A(N, N, 2, 2, n);
274 for( int i=0; i<n; i++)
275 for( int j=0; j<n; j++)
276 {
277 A.data[(0*n+i)*n+j] = bp(i,j);
278 A.data[(1*n+i)*n+j] = a(i,j);
279 }
280 for( int i=0; i<N; i++)
281 for( int d=0; d<2; d++)
282 {
283 A.data_idx[i*2+d] = d; //bp, a
284 A.cols_idx[i*2+d] = (i+d-1+N)%N; //-1, 0
285 }
286 return A;
287 }
288};
289
300template<class real_type>
301EllSparseBlockMat<real_type, thrust::host_vector> jump( int n, int N, real_type h, bc bcx)
302{
303 constexpr int invalid_idx = EllSparseBlockMat<real_type, thrust::host_vector>::invalid_index;
304 SquareMatrix<real_type> l = create::lilj<real_type>(n);
305 SquareMatrix<real_type> r = create::rirj<real_type>(n);
306 SquareMatrix<real_type> lr = create::lirj<real_type>(n);
307 SquareMatrix<real_type> rl = create::rilj<real_type>(n);
308 SquareMatrix<real_type> a = l+r;
309 SquareMatrix<real_type> a_bound_left = a;//DIR and PER
310 if( bcx == NEU || bcx == NEU_DIR)
311 a_bound_left = r;
312 SquareMatrix<real_type> a_bound_right = a; //DIR and PER
313 if( bcx == NEU || bcx == DIR_NEU)
314 a_bound_right = l;
315 SquareMatrix<real_type> b = -rl;
316 SquareMatrix<real_type> bp = -lr;
317 //transform to XSPACE
318 SquareMatrix<real_type> t = create::pipj_inv<real_type>(n);
319 t *= 2./h;
320 SquareMatrix<real_type> backward=dg::DLT<real_type>::backward(n);
321 SquareMatrix<real_type> forward=dg::DLT<real_type>::forward(n);
322 a = backward*t*a*forward, a_bound_left = backward*t*a_bound_left*forward;
323 b = backward*t*b*forward, a_bound_right = backward*t*a_bound_right*forward;
324 bp = backward*t*bp*forward;
325 //assemble the matrix
326 if(bcx != dg::PER)
327 {
328 EllSparseBlockMat<real_type, thrust::host_vector> A(N, N, 3, 5, n);
329 for( int i=0; i<n; i++)
330 for( int j=0; j<n; j++)
331 {
332 A.data[(0*n+i)*n+j] = bp(i,j);
333 A.data[(1*n+i)*n+j] = a(i,j);
334 A.data[(2*n+i)*n+j] = b(i,j);
335 A.data[(3*n+i)*n+j] = a_bound_left(i,j);
336 A.data[(4*n+i)*n+j] = a_bound_right(i,j);
337 }
338 A.data_idx[0*3+0] = 3; //a_bound_left
339 A.cols_idx[0*3+0] = 0;
340 A.data_idx[0*3+1] = 2; //b
341 A.cols_idx[0*3+1] = 1;
342 A.data_idx[0*3+2] = 2; //0
343 A.cols_idx[0*3+2] = invalid_idx; //prevent unnecessary data fetch
344 for( int i=1; i<N-1; i++) //a
345 for( int d=0; d<3; d++)
346 {
347 A.data_idx[i*3+d] = d; //bp, a, b
348 A.cols_idx[i*3+d] = i+d-1;
349 }
350 A.data_idx[(N-1)*3+0] = 0; //bp
351 A.cols_idx[(N-1)*3+0] = N-2;
352 A.data_idx[(N-1)*3+1] = 4; //a_bound_right
353 A.cols_idx[(N-1)*3+1] = N-1;
354 A.data_idx[(N-1)*3+2] = 4; //0
355 A.cols_idx[(N-1)*3+2] = invalid_idx; //prevent unnecessary data fetch
356 return A;
357
358 }
359 else //periodic
360 {
361 EllSparseBlockMat<real_type, thrust::host_vector> A(N, N, 3, 3, n);
362 for( int i=0; i<n; i++)
363 for( int j=0; j<n; j++)
364 {
365 A.data[(0*n+i)*n+j] = bp(i,j);
366 A.data[(1*n+i)*n+j] = a(i,j);
367 A.data[(2*n+i)*n+j] = b(i,j);
368 }
369 for( int i=0; i<N; i++)
370 for( int d=0; d<3; d++)
371 {
372 A.data_idx[i*3+d] = d; //bp, a, b
373 A.cols_idx[i*3+d] = (i+d-1+N)%N;
374 }
375 return A;
376 }
377};
378
390template<class real_type>
391EllSparseBlockMat<real_type, thrust::host_vector> dx_normed( int n, int N, real_type h, bc bcx, direction dir )
392{
393 if( dir == centered)
394 return create::dx_symm(n, N, h, bcx);
395 else if (dir == forward)
396 return create::dx_plus(n, N, h, bcx);
397 else if (dir == backward)
398 return create::dx_minus(n, N, h, bcx);
399 return EllSparseBlockMat<real_type, thrust::host_vector>();
400}
402
403} //namespace create
404} //namespace dg
405
functions.h
Some utility functions for the dg::evaluate routines.
grid.h
base topology classes
dg::create::jump
EllSparseBlockMat< real_type, thrust::host_vector > jump(unsigned coord, const aRealTopology< real_type, Nd > &g, dg::bc bc)
Create a jump matrix along given coordinate.
Definition derivatives.h:70
dg::bc
bc
Switch between boundary conditions.
Definition enums.h:15
dg::direction
direction
Direction of a discrete derivative.
Definition enums.h:97
dg::NEU_DIR
@ NEU_DIR
Neumann on left, Dirichlet on right boundary.
Definition enums.h:19
dg::PER
@ PER
periodic boundaries
Definition enums.h:16
dg::NEU
@ NEU
Neumann on both boundaries.
Definition enums.h:20
dg::DIR
@ DIR
homogeneous dirichlet boundaries
Definition enums.h:17
dg::DIR_NEU
@ DIR_NEU
Dirichlet on left, Neumann on right boundary.
Definition enums.h:18
dg::backward
@ backward
backward derivative (cell to the left and current cell)
Definition enums.h:99
dg::forward
@ forward
forward derivative (cell to the right and current cell)
Definition enums.h:98
dg::centered
@ centered
centered derivative (cell to the left and right and current cell)
Definition enums.h:100
dg
This is the namespace for all functions and classes defined and used by the discontinuous Galerkin li...
dg::DLT::backward
static std::vector< real_type > backward(unsigned n)
Return backward DLT trafo matrix.
Definition dlt.h:139
dg::DLT::forward
static std::vector< real_type > forward(unsigned n)
Return forward DLT trafo matrix.
Definition dlt.h:195
dg::EllSparseBlockMat::invalid_index
static constexpr int invalid_index
Definition sparseblockmat.h:49
weights.h
Creation functions for integration weights and their inverse.