@@ -57,11 +57,11 @@ namespace libmpdataxx
5757
5858 void fill_halos_vctr_alng (arrvec_t <arr_t > &av, const rng_t &j, const bool ad = false )
5959 {
60- using namespace idxperm ;
60+ using namespace idxperm ;
6161 // zero velocity condition
6262 for (int i = this ->left_halo_vctr .first (), n = halo; i <= this ->left_halo_vctr .last () - (ad ? 1 : 0 ); ++i, --n)
6363 {
64- av[d](pi<d>(i, j)) = -av[d](pi<d>(this ->left_edge_sclr + n - h, j));
64+ av[d](pi<d>(i, j)) = -av[d](pi<d>(this ->left_edge_sclr + n - h, j));
6565 }
6666 }
6767
@@ -74,8 +74,32 @@ namespace libmpdataxx
7474 void fill_halos_flux (arrvec_t <arr_t > &av, const rng_t &j)
7575 {
7676 using namespace idxperm ;
77- av[d](pi<d>(this ->left_halo_vctr .last (), j)) = -av[d](pi<d>(this ->left_edge_sclr + h, j));
77+ av[d](pi<d>(this ->left_halo_vctr .last (), j)) = -av[d](pi<d>(this ->left_edge_sclr + h, j));
78+ }
79+
80+ void fill_halos_sgs_div (arr_t &a, const rng_t &j)
81+ {
82+ using namespace idxperm ;
83+ a (pi<d>(this ->left_edge_sclr - h, j)) = 2 * a (pi<d>(this ->left_edge_sclr + h, j)) - a (pi<d>(this ->left_edge_sclr + 1 + h, j));
84+ }
85+
86+ void fill_halos_sgs_vctr (arrvec_t <arr_t > &av, const arr_t &, const rng_t &j, const int offset = 0 )
87+ {
88+ // fill halos for a staggered field so that it has zero value on tke edge
89+ // that is 0.5 * (a(edge-h) + a(edge+h)) = 0
90+ using namespace idxperm ;
91+ const auto &a = av[offset + d];
92+ a (pi<d>(this ->left_edge_sclr - h, j)) = - a (pi<d>(this ->left_edge_sclr + h, j));
7893 }
94+
95+ void fill_halos_sgs_tnsr (arrvec_t <arr_t > &av, const arr_t &, const arr_t &, const rng_t &j, const real_t di)
96+ {
97+ using namespace idxperm ;
98+ const auto &a = av[d];
99+ a (pi<d>(this ->left_edge_sclr - h, j)) = 2 * a (pi<d>(this ->left_edge_sclr + h, j))
100+ - a (pi<d>(this ->left_edge_sclr + 1 + h, j));
101+ }
102+
79103 };
80104
81105 template <typename real_t , int halo, bcond_e knd, drctn_e dir, int n_dims, int d>
@@ -124,11 +148,11 @@ namespace libmpdataxx
124148
125149 void fill_halos_vctr_alng (arrvec_t <arr_t > &av, const rng_t &j, const bool ad = false )
126150 {
127- using namespace idxperm ;
151+ using namespace idxperm ;
128152 // zero velocity condition
129153 for (int i = this ->rght_halo_vctr .first () + (ad ? 1 : 0 ), n = 1 ; i <= this ->rght_halo_vctr .last (); ++i, ++n)
130154 {
131- av[d](pi<d>(i, j)) = -av[d](pi<d>(this ->rght_edge_sclr - n + h, j));
155+ av[d](pi<d>(i, j)) = -av[d](pi<d>(this ->rght_edge_sclr - n + h, j));
132156 }
133157 }
134158
@@ -140,9 +164,32 @@ namespace libmpdataxx
140164
141165 void fill_halos_flux (arrvec_t <arr_t > &av, const rng_t &j)
142166 {
143- using namespace idxperm ;
167+ using namespace idxperm ;
144168 // zero flux condition
145- av[d](pi<d>(this ->rght_halo_vctr .first (), j)) = -av[d](pi<d>(this ->rght_edge_sclr - h, j));
169+ av[d](pi<d>(this ->rght_halo_vctr .first (), j)) = -av[d](pi<d>(this ->rght_edge_sclr - h, j));
170+ }
171+
172+ void fill_halos_sgs_div (arr_t &a, const rng_t &j)
173+ {
174+ using namespace idxperm ;
175+ a (pi<d>(this ->rght_edge_sclr + h, j)) = 2 * a (pi<d>(this ->rght_edge_sclr - h, j)) - a (pi<d>(this ->rght_edge_sclr - 1 - h, j));
176+ }
177+
178+ void fill_halos_sgs_vctr (arrvec_t <arr_t > &av, const arr_t &, const rng_t &j, const int offset = 0 )
179+ {
180+ // fill halos for a staggered field so that it has zero value on tke edge
181+ // that is 0.5 * (a(edge-h) + a(edge+h)) = 0
182+ using namespace idxperm ;
183+ const auto &a = av[offset + d];
184+ a (pi<d>(this ->rght_edge_sclr + h, j)) = - a (pi<d>(this ->rght_edge_sclr - h, j));
185+ }
186+
187+ void fill_halos_sgs_tnsr (arrvec_t <arr_t > &av, const arr_t &, const arr_t &, const rng_t &j, const real_t di)
188+ {
189+ using namespace idxperm ;
190+ const auto &a = av[d];
191+ a (pi<d>(this ->rght_edge_sclr + h, j)) = 2 * a (pi<d>(this ->rght_edge_sclr - h, j))
192+ - a (pi<d>(this ->rght_edge_sclr - 1 - h, j));
146193 }
147194 };
148195 } // namespace bcond
0 commit comments