LoAdSG
RHS.h
1//
2// Created by to35jepo on 3/13/23.
3//
4
5#ifndef SGRUN_RHS_H
6#define SGRUN_RHS_H
7#include "../extemp/vector.h"
8#include "../extemp/multilevelvector.h"
9#include "../stencils/MassStencil.h"
10#include "MatrixVectorInhomogen.h"
11#include "../iterator/RectangularIterator.h"
12#include "../stencils/InterfaceIntegration/IntegratorRHS/integratorRHS.h"
13
14
15
16void CalcHierarchicalBasisForRHS(VectorSparseG &u);
17
18
19class InHomoBoundaryRHS{
20public:
21 InHomoBoundaryRHS(AdaptiveSparseGrid& grid, MultiLevelAdaptiveSparseGrid& mgrid): depthList(grid), nodalMultiLevel(mgrid), nodal(grid){
22 bool prolongations[DimensionSparseGrid]={0};
23 depthList.sortDepthsNeumann(prolongations);
24 };
25
26 void dualTransform(VectorSparseG& neumann, VectorSparseG& dirichlet){
27 dirichlet=0.0;
28 IndexDimension minI;
29 IndexDimension maxI;
30
31 std::list<Depth>* sortedDepths = depthList.getSortierteTiefen();
32 for(Depth T: *sortedDepths) nodalMultiLevel.setMultiLevelValues2(neumann, T);
33
34
35 for(int d=0; d < DimensionSparseGrid; d++) {
36
37 for (Depth T: *sortedDepths) {
38
39
40 Depth T_fine=T;
41 int t = T.at(d) + 1;
42 T_fine.set(t,d);
43
44
45 if(depthList.isIncluded(T_fine)) {
46
47 Matrix M(t);
48
49 IndexDimension I;
50
51 for (int dd = 0; dd < DimensionSparseGrid; dd++) {
52 minI.replace(dd, 0);
53 maxI.replace(dd, 1);
54 }
55 maxI.replace(d,0);
56
57 double PNeumann[M.size_all];
58 double PNodal[M.size_all];
59 unsigned long indices[M.size_all];
60
61
62 for (RectangularIterator iter(minI, maxI, T); iter.goon();++iter) {
63 I = iter.getIndex();
64
65 for (int i = 0; i < M.size_all; i++) {
66
67
68 unsigned long k;
69 nodalMultiLevel.getSparseGrid()->occupied(k, I, T_fine);
70 indices[i]=k;
71 if(i%2==0){
72 nodalMultiLevel.getSparseGrid()->occupied(k, I, T);
73 PNeumann[i]=nodalMultiLevel.getValue(k);
74 } else {
75 PNeumann[i]=nodalMultiLevel.getValue(indices[i]);
76 }
77
78 I = I.nextRight(d, T_fine.at(d));
79 }
80
81 M.solve(PNodal,PNeumann);
82
83 for (int i = 0; i < M.size_all; i++) {
84
85 nodalMultiLevel.setValue(indices[i], PNodal[i]);
86
87 }
88
89 }
90
91 }
92
93 }
94
95 }
96
97 for(Depth T: *sortedDepths){
98 if(T>>0) {
99 nodal.setMultiLevelValues(nodalMultiLevel, T);
100 ConvertDualNodalToPrewavelet(nodal, dirichlet, T);
101
102 }
103 }
104 }
105
106 class Matrix{
107 public:
108 Matrix(int t){
109
110
111 size_all = POW2(t)+1;
112 if (size_all < 6 && size_all >= 1) {
113 size = size_all*size_all;
114 }else {
115 size = size_all * 5;
116 }
117
118 M = new double[size];
119 L = new double[size];
120 R = new double[size];
121
122
123 for(int k=0; k<size; k++){
124 M[k]=0.0;
125 L[k]=0.0;
126 R[k]=0.0;
127 }
128
129
130 if(t>1) {
131 for (int k = 2; k < size_all - 1; k += 2) {
132 M[convert(k, k, size_all)] = 1.0;
133 M[convert(k, k - 1, size_all)] = 0.5;
134 M[convert(k, k + 1, size_all)] = 0.5;
135 }
136 for (int k = 3; k < size_all - 2; k += 2) {
137 M[convert(k, k, size_all)] = 1.0;
138 M[convert(k, k - 1, size_all)] = -0.6;
139 M[convert(k, k - 2, size_all)] = 0.1;
140 M[convert(k, k + 1, size_all)] = -0.6;
141 M[convert(k, k + 2, size_all)] = 0.1;
142 }
143
144 M[convert(0, 0, size_all)] = 1.0;
145 M[convert(0, 1, size_all)] = 0.5;
146 M[convert(1, 0, size_all)] = -1.2;
147 M[convert(1, 1, size_all)] = 1.1;
148 M[convert(1, 2, size_all)] = -0.6;
149 M[convert(1, 3, size_all)] = 0.1;
150
151 M[convert(size_all - 1, size_all - 1, size_all)] = 1.0;
152 M[convert(size_all - 1, size_all - 2, size_all)] = 0.5;
153 M[convert(size_all - 2, size_all - 1, size_all)] = -1.2;
154 M[convert(size_all - 2, size_all - 2, size_all)] = 1.1;
155 M[convert(size_all - 2, size_all - 3, size_all)] = -0.6;
156 M[convert(size_all - 2, size_all - 4, size_all)] = 0.1;
157 }
158 if(t==1){
159 M[convert(0,0,size_all)]=1.0;
160 M[convert(0,1,size_all)]=0.5;
161 M[convert(1,0,size_all)]=-1.0;
162 M[convert(1,1,size_all)]=1.0;
163 M[convert(1,2,size_all)]=-1.0;
164 M[convert(2,1,size_all)]=0.5;
165 M[convert(2,2,size_all)]=1.0;
166 }
167
168 if(t==0){ M[convert(0,0,size_all)]=1.0; M[convert(1,1,size_all)]=1.0;}
169
170 LRcreated = false;
171
172 }
173 ~Matrix(){
174
175 delete[] M;
176 delete[] L;
177 delete[] R;
178
179 }
180 void print() {
181 int I = size_all;
182 for (int i = 0; i < I; i++) {
183 for (int j = 0; j < I; j++) {
184 if (i <= j + 2 && j <= i + 2) {
185 cout << "\t" << M[convert(i, j, I)];
186 } else { cout << "\t"; }
187 }
188 cout << "\n";
189 }
190 cout << endl;
191
192 };
193
194 void solve(double *x, double *b) {
195 if (size_all < 1) {
196 x[0] = b[0];
197 return;
198 }
199 if (!LRcreated) createLR();
200 double y[size_all];
201 for (int i = 0; i < size_all; i++) {
202 y[i] = b[i];
203 }
204
205 forward(y, b, L, size_all);
206 backward(x, y, R, size_all);
207 };
208
209
210 int size_all;
211 private:
212 void createLR() {
213 for (int i = 0; i < size_all; i++) {
214 for (int j = 0; j < size_all; j++) {
215 if (abs(i - j) < 3) {
216 R[convert(i, j, size_all)] = M[convert(i, j, size_all)];
217 }
218 }
219 }
220
221 LR(L, R, size_all);
222 LRcreated = true;
223 };
224
225 int convert(int i,int j,int I){
226
227
228 int savei = i;
229 if (I < 6) { return i * I + j; }
230 if (j > 2 && j < I - 2) {
231 i = i - ((j - 3) + 1);
232 }
233 if (j == I - 2) i = i - j + 3;
234 if (j == I - 1) i = i - j + 4;
235
236 if (I > 5) {
237 if (i > 4 || i < 0) {
238 cout << "error" << endl;
239 cout << " i " << i << " j " << j << endl;
240 cout << "save i " << savei << " I " << I << endl;
241 }
242 }
243
244
245 return i * I + j;
246 }
247
248
249 double *M;
250 double *R;
251 double *L;
252 bool LRcreated;
253
254 int size;
255
256 };
257
258
259
260
261private:
262
263 DepthList depthList;
264 MultiLevelVector nodalMultiLevel;
265 VectorSparseG nodal;
266
267
268
269};
270
271
272
273
274
275
276
277
278class InHomoBoundaryRHSHierarchical{
279public:
280 InHomoBoundaryRHSHierarchical(AdaptiveSparseGrid& grid):nodal(grid){};
281
282 void multiply(VectorSparseG& hierarchical, VectorSparseG& dirichlet);
283
284 void multiply_mass(VectorSparseG& hierarchical, VectorSparseG& dirichlet);
285
286 void multiply_mass_coeff(VectorSparseG &hierarchical, VectorSparseG &dirichlet);
287
288 template<class F>
289 void multiply_mass_coeff_interface(VectorSparseG &hierarchical, VectorSparseG &dirichlet,const F& var_coeff) {
290
291
292 DepthList depthList(*dirichlet.getSparseGrid());
293 DepthList depthList1(*hierarchical.getSparseGrid());
294
295 for (auto it = depthList.begin_all(); it != depthList.end_all(); ++it) {
296 Depth T = *it;
297 auto iter = GetNextSmallerDepthsIterator(T);
298 do {
299 Depth Tlocal = *iter;
300 if (Tlocal >> 0) {
301
302 SingleDepthHashGrid &depthGrid = dirichlet.getSparseGrid()->getMultiDepthHashGrid()->getGridForDepth(
303 Tlocal);
304 const auto &mapping = depthGrid._mapPosToGridPos;
305
306 if (depthGrid.getNumberOfEntries() > 0) {
307
308
309 for (size_t i = 0; i < mapping.size(); i++) {
310
311 IndexDimension INodal = depthGrid._map.getIndexOfTable(i);
312
313 //double sum = 0.0;
314 double sum_omp = 0.0;
315 double sum = 0.0;
316
317
318 for (auto it_inner = depthList1.begin_all();
319 it_inner != depthList1.end_all(); ++it_inner) {
320
321 Depth THier = *it_inner;
322 if (!(THier >> 0)) {
323
324 SingleDepthHashGrid &depthGrid_hier = hierarchical.getSparseGrid()->getMultiDepthHashGrid()->getGridForDepth(
325 THier);
326 const auto &mapping_hier = depthGrid_hier._mapPosToGridPos;
327 const auto &map_hier = depthGrid_hier._map;
328 auto end = mapping_hier.size();
329
330
331 if (depthGrid_hier.getNumberOfEntries() > 0){
332 #pragma omp parallel reduction(+:sum_omp)
333 {
334
335 #pragma omp for schedule(runtime)
336 for (size_t j = 0; j < end; j++) {
337 //IndexDimension IHier = depthGrid_hier._map.getIndexOfTable(j);
338 IndexDimension IHier = map_hier.getIndexOfTable(j);
339
340 array<double, DimensionSparseGrid> p_u{}, p_v{}, l_b{}, r_b{}, h_u{}, h_v{};
341 array<int, DimensionSparseGrid> d_u{}, d_v{};
342
343 for (int d = 0; d < DimensionSparseGrid; d++) {
344 p_u[d] = INodal.coordinate(d);
345 p_v[d] = IHier.coordinate(d);
346 h_u[d] = 1.0 / double(POW2(T.at(d)));
347 h_v[d] = 1.0 / double(POW2(THier.at(d)));
348 l_b[d] = max(p_u[d] - h_u[d], p_v[d] - h_v[d]);
349 l_b[d] = max(l_b[d], 0.0);
350 r_b[d] = min(p_u[d] + h_u[d], p_v[d] + h_v[d]);
351 r_b[d] = min(r_b[d], 1.0);
352 d_u[d] = T.at(d);
353 d_v[d] = THier.at(d);
354 }
355
356 IntegratorRHS<const F, DimensionSparseGrid> integratorRhs(var_coeff,
357 1e-4, 0, 3, 1.0e8);
358 double val = integratorRhs.integration(l_b, r_b,
359 p_u, p_v,
360 h_u,
361 h_v,
362 d_u, d_v);
363
364 //#pragma omp critical
365 sum_omp += val *
366 hierarchical.getValue(mapping_hier[j]);
367
368
369 }
370
371 }
372
373 }
374 }
375 }
376
377 /* if(abs(sum-sum_omp)>1e-10){
378 cout << sum << " " << sum_omp << " " << sum_omp-sum << endl;
379 exit(1);
380
381 }*/
382 nodal.setValue(mapping[i], sum_omp);
383
384 }
385 }
386 }
387
388
389 } while (iter.next());
390
391 ConvertToPrewavelet2(nodal, dirichlet, T);
392 nodal = 0.0;
393 }
394
395 }
396private:
397 inline void ConvertToPrewavelet2(VectorSparseG &Ax_hier, VectorSparseG &Ax, Depth &T) {
398 SingleDepthHashGrid& depthGrid = Ax_hier.getSparseGrid()->getMultiDepthHashGrid()->getGridForDepth(T);
399 const auto& mapping = depthGrid._mapPosToGridPos;
400 // cout << mapping.size() << ", " << depthGrid.getNumberOfEntries() <<endl;
401 for (size_t i = 0; i < mapping.size(); i++)
402 {
403 if(Ax_hier.getSparseGrid()->getActiveTable()[mapping[i]]) {
404 IndexDimension I = depthGrid._map.getIndexOfTable(i);
405 double coeff = 0.0;
406 IndexDimension J;
407
408 for (MultiDimFiveCompass mc; mc.goon(); ++mc) {
409 double basis_coeff = 1.0;
410 J = I.nextFiveP(&mc, T, &basis_coeff);
411 if (mc.goon())
412 coeff = coeff + Ax_hier.getValue(J) * basis_coeff;
413 }
414
415 Ax.setValue(mapping[i], coeff);
416 }
417 }
418 }
419
420
421 VectorSparseG nodal;
422
423
424
425};
426
427#endif //SGRUN_RHS_H
AdaptiveSparseGrid
Definition sparseGrid.h:277
IntegratorRHS
Integrator for the Helmholtz operator with variable coefficients.
Definition integratorRHS.h:26
MultiDimFiveCompass
Definition index.h:436
RectangularIterator
Definition RectangularIterator.h:13