OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
Eno_Weno.hpp
1/*
2 *
3 * Based on the description in Osher and Fedkiw, Level Set Methods and Dynamic Implicit Surfaces
4 * For HJ ENO method - Sec. 3.3
5 * For HJ WENO method - Sec. 3.4
6*/
7
8// Created by Sachin, modified by Justina
9
10#ifndef __ENO_WENO_HPP__
11#define __ENO_WENO_HPP__
12
13#include "Grid/grid_dist_id.hpp"
14
15static double adjustWeights(double v1, double v2, double v3, double v4, double v5)
16{
17 double phix1, phix2, phix3;
18 double s1, s2, s3;
19 double a1, a2, a3;
20 double eps = 1e-6;
21 double w1, w2, w3;
22
23 // Eqs. (3.25), (3.26), (3.27)
24 phix1 = (v1 / 3.0) - (7.0 * v2 / 6.0) + (11.0 * v3 / 6.0);
25 phix2 = -(v2 / 6.0) + (5.0 * v3 / 6.0) + (v4 / 3.0);
26 phix3 = (v3 / 3.0) + (5.0 * v4 / 6.0) - (v5 / 6.0);
27
28 // Eqs. (3.32), (3.33), (3.34)
29 s1 = (13.0 / 12.0) * (v1 - 2*v2 + v3) * (v1 - 2*v2 + v3)
30 + 0.25 * (v1 - 4*v2 + 3*v3) * (v1 - 4*v2 + 3*v3);
31 s2 = (13.0 / 12.0) * (v2 - 2*v3 + v4) * (v2 - 2*v3 + v4)
32 + 0.25 * (v2 - v4) * (v2 - v4);
33 s3 = (13.0 / 12.0) * (v3 - 2*v4 + v5) * (v3 - 2*v4 + v5)
34 + 0.25 * (3*v3 - 4*v4 + v5) * (3*v3 - 4*v4 + v5);
35
36 // Eqs. (3.35), (3.36), (3.37)
37 a1 = 0.1 / ((s1 + eps)*(s1 + eps));
38 a2 = 0.6 / ((s2 + eps)*(s2 + eps));
39 a3 = 0.3 / ((s3 + eps)*(s3 + eps));
40
41 // Eqs. (3.39), (3.40), (3.41)
42 w1 = a1 / (a1 + a2 + a3);
43 w2 = a2 / (a1 + a2 + a3);
44 w3 = a3 / (a1 + a2 + a3);
45
46 // Eq. (3.28)
47 return (w1 * phix1 + w2 * phix2 + w3 * phix3);
48}
49
50template<size_t Field, typename grid_type, typename key_type>
51double WENO_5_Plus(grid_type & grid, key_type key, size_t x)
52{
53 double c2, c3, c4, c5, c6;
54 double coeff = 1.0 / grid.spacing(x);
55 double df;
56
57 c2 = (grid.template get<Field>(key.move(x, -1)) - grid.template get<Field>(key.move(x, -2))) * coeff;
58 c3 = (grid.template get<Field>(key) - grid.template get<Field>(key.move(x, -1))) * coeff;
59 c4 = (grid.template get<Field>(key.move(x,1)) - grid.template get<Field>(key)) * coeff;
60 c5 = (grid.template get<Field>(key.move(x, 2)) - grid.template get<Field>(key.move(x, 1))) * coeff;
61 c6 = (grid.template get<Field>(key.move(x, 3)) - grid.template get<Field>(key.move(x, 2))) * coeff;
62 df = adjustWeights(c6, c5, c4, c3, c2);
63
64 return df;
65}
66
67template<size_t Field, typename grid_type, typename key_type>
68double WENO_5_Minus(grid_type & grid, key_type key, size_t x)
69{
70 double c1, c2, c3, c4, c5;
71 double coeff = 1.0 / grid.spacing(x);
72 double df;
73
74 c1 = (grid.template get<Field>(key.move(x, -2)) - grid.template get<Field>(key.move(x, -3))) * coeff;
75 c2 = (grid.template get<Field>(key.move(x, -1)) - grid.template get<Field>(key.move(x, -2))) * coeff;
76 c3 = (grid.template get<Field>(key) - grid.template get<Field>(key.move(x, -1))) * coeff;
77 c4 = (grid.template get<Field>(key.move(x,1)) - grid.template get<Field>(key)) * coeff;
78 c5 = (grid.template get<Field>(key.move(x, 2)) - grid.template get<Field>(key.move(x, 1))) * coeff;
79
80 df = adjustWeights(c1, c2, c3, c4, c5);
81
82 return df;
83}
84
85template<size_t Field, typename grid_type, typename key_type>
86double ENO_3_Plus(grid_type & grid, key_type key, size_t x)
87{
88 double q1x, q2x, q3x;
89 double gridsize = grid.spacing(x);
90 double coeff = 1.0 / gridsize;
91
92 q1x = (grid.template get<Field>(key.move(x,1)) - grid.template get<Field>(key)) * coeff;
93 double d2ix = (grid.template get<Field>(key.move(x,1)) - 2*grid.template get<Field>(key) + grid.template get<Field>(key.move(x,-1))) * 0.5 * coeff * coeff;
94 double d2ip1x = (grid.template get<Field>(key.move(x,2)) - 2*grid.template get<Field>(key.move(x,1)) + grid.template get<Field>(key)) * 0.5 * coeff * coeff;
95 if(abs(d2ix) <= abs(d2ip1x))
96 {
97 q2x = -d2ix * gridsize; // (3.22)
98 double d3p1hx = (grid.template get<Field>(key.move(x,1)) - 3*grid.template get<Field>(key) + 3*grid.template get<Field>(key.move(x,-1)) - grid.template get<Field>(key.move(x,-2))) * coeff * coeff * coeff / 6.0;
99 double d3p3hx = (grid.template get<Field>(key.move(x,2)) - 3*grid.template get<Field>(key.move(x,1)) + 3*grid.template get<Field>(key) - grid.template get<Field>(key.move(x,-1))) * coeff * coeff * coeff / 6.0;
100 if(abs(d3p1hx) <= abs(d3p3hx))
101 q3x = -d3p1hx * gridsize * gridsize;
102 else
103 q3x = -d3p3hx * gridsize * gridsize;
104 }
105 else
106 {
107 q2x = -d2ip1x * gridsize;
108 double d3p1hx = (grid.template get<Field>(key.move(x,2)) - 3*grid.template get<Field>(key.move(x,1)) + 3*grid.template get<Field>(key) - grid.template get<Field>(key.move(x,-1))) * coeff * coeff * coeff / 6.0;
109 double d3p3hx = (grid.template get<Field>(key.move(x,3)) - 3*grid.template get<Field>(key.move(x,2)) + 3*grid.template get<Field>(key.move(x,1)) - grid.template get<Field>(key)) * coeff * coeff * coeff / 6.0;
110 if(abs(d3p1hx) <= abs(d3p3hx))
111 q3x = d3p1hx * 2 * gridsize * gridsize;
112 else
113 q3x = d3p3hx * 2 * gridsize * gridsize;
114
115 }
116 return (q1x + q2x + q3x);
117}
118
119template<size_t Field, typename grid_type, typename key_type>
120double ENO_3_Minus(grid_type & grid, key_type key, size_t x)
121{
122 double q1x, q2x, q3x;
123 double gridsize = grid.spacing(x);
124 double coeff = 1.0 / gridsize;
125
126 q1x = (grid.template get<Field>(key) - grid.template get<Field>(key.move(x, -1))) * coeff;
127 double d2im1x = (grid.template get<Field>(key) - 2*grid.template get<Field>(key.move(x, -1)) + grid.template get<Field>(key.move(x, -2))) * 0.5 * coeff * coeff;
128 double d2ix = (grid.template get<Field>(key.move(x,1)) - 2*grid.template get<Field>(key) + grid.template get<Field>(key.move(x,-1))) * 0.5 * coeff * coeff;
129
130 if(abs(d2im1x) <= abs(d2ix))
131 {
132 q2x = d2im1x * gridsize;
133 double d3p1hx = (grid.template get<Field>(key) - 3*grid.template get<Field>(key.move(x,-1)) + 3*grid.template get<Field>(key.move(x,-2)) - grid.template get<Field>(key.move(x,-3))) * coeff * coeff * coeff / 6.0;
134 double d3p3hx = (grid.template get<Field>(key.move(x,1)) - 3*grid.template get<Field>(key) + 3*grid.template get<Field>(key.move(x,-1)) - grid.template get<Field>(key.move(x,-2))) * coeff * coeff * coeff / 6.0;
135 if(abs(d3p1hx) <= abs(d3p3hx))
136 q3x = d3p1hx * 2 * gridsize * gridsize;
137 else
138 q3x = d3p3hx * 2 * gridsize * gridsize;
139 }
140 else
141 {
142 q2x = d2ix * gridsize;
143 double d3p1hx = (grid.template get<Field>(key.move(x,1)) - 3*grid.template get<Field>(key) + 3*grid.template get<Field>(key.move(x,-1)) - grid.template get<Field>(key.move(x,-2))) * coeff * coeff * coeff / 6.0;
144 double d3p3hx = (grid.template get<Field>(key.move(x,2)) - 3*grid.template get<Field>(key.move(x,1)) + 3*grid.template get<Field>(key) - grid.template get<Field>(key.move(x,-1))) * coeff * coeff * coeff / 6.0;
145 if(abs(d3p1hx) <= abs(d3p3hx))
146 q3x = -d3p1hx * gridsize * gridsize;
147 else
148 q3x = -d3p3hx * gridsize * gridsize;
149 }
150
151 return (q1x + q2x + q3x);
152}
153
154#endif
155
grid_dist_id
This is a distributed grid.
Definition grid_dist_id.hpp:278