1 #include<Kripke/Kernel/Kernel_3d_GDZ.h>
2 #include<Kripke/Grid.h>
3 #include<Kripke/SubTVec.h>
4
nestingPsi(void) const5 Nesting_Order Kernel_3d_GDZ::nestingPsi(void) const {
6 return NEST_GDZ;
7 }
8
nestingPhi(void) const9 Nesting_Order Kernel_3d_GDZ::nestingPhi(void) const {
10 return NEST_GDZ;
11 }
12
nestingSigt(void) const13 Nesting_Order Kernel_3d_GDZ::nestingSigt(void) const {
14 return NEST_DGZ;
15 }
16
nestingEll(void) const17 Nesting_Order Kernel_3d_GDZ::nestingEll(void) const {
18 return NEST_ZGD;
19 }
20
nestingEllPlus(void) const21 Nesting_Order Kernel_3d_GDZ::nestingEllPlus(void) const {
22 return NEST_ZDG;
23 }
24
nestingSigs(void) const25 Nesting_Order Kernel_3d_GDZ::nestingSigs(void) const {
26 return NEST_GZD;
27 }
28
29
LTimes(Grid_Data * grid_data)30 void Kernel_3d_GDZ::LTimes(Grid_Data *grid_data) {
31 // Outer parameters
32 int nidx = grid_data->total_num_moments;
33
34 // Clear phi
35 for(int ds = 0;ds < grid_data->num_zone_sets;++ ds){
36 grid_data->phi[ds]->clear(0.0);
37 }
38
39 // Loop over Subdomains
40 int num_subdomains = grid_data->subdomains.size();
41 for (int sdom_id = 0; sdom_id < num_subdomains; ++ sdom_id){
42 Subdomain &sdom = grid_data->subdomains[sdom_id];
43
44 // Get dimensioning
45 int num_zones = sdom.num_zones;
46 int num_local_groups = sdom.num_groups;
47 int group0 = sdom.group0;
48 int num_local_directions = sdom.num_directions;
49
50 /* 3D Cartesian Geometry */
51 double * KRESTRICT phi = sdom.phi->ptr(group0, 0, 0);
52 double * KRESTRICT psi_ptr = sdom.psi->ptr();
53
54 for (int group = 0; group < num_local_groups; ++group) {
55 double * KRESTRICT ell_nm = sdom.ell->ptr();
56
57 for(int nm_offset = 0;nm_offset < nidx;++nm_offset){
58 double * KRESTRICT psi = psi_ptr;
59
60 for (int d = 0; d < num_local_directions; d++) {
61 double ell_nm_d = ell_nm[d];
62
63 #ifdef KRIPKE_USE_OPENMP
64 #pragma omp parallel for
65 #endif
66 for(int z = 0;z < num_zones; ++ z){
67 phi[z] += ell_nm_d * psi[z];
68 }
69
70 psi += num_zones;
71 }
72 ell_nm += num_local_directions;
73 phi += num_zones;
74 }
75 psi_ptr += num_zones*num_local_directions;
76 }
77 } // Subdomain
78 }
79
LPlusTimes(Grid_Data * grid_data)80 void Kernel_3d_GDZ::LPlusTimes(Grid_Data *grid_data) {
81 // Outer parameters
82 int nidx = grid_data->total_num_moments;
83
84 // Loop over Subdomains
85 int num_subdomains = grid_data->subdomains.size();
86 for (int sdom_id = 0; sdom_id < num_subdomains; ++ sdom_id){
87 Subdomain &sdom = grid_data->subdomains[sdom_id];
88
89 // Get dimensioning
90 int num_zones = sdom.num_zones;
91 int num_local_groups = sdom.num_groups;
92 int group0 = sdom.group0;
93 int num_local_directions = sdom.num_directions;
94
95 // Get Variables
96 sdom.rhs->clear(0.0);
97
98 /* 3D Cartesian Geometry */
99 double *ell_plus_ptr = sdom.ell_plus->ptr();
100
101 double * KRESTRICT phi_out_ptr = sdom.phi_out->ptr(group0, 0, 0);
102 double * KRESTRICT rhs = sdom.rhs->ptr();
103
104 for (int group = 0; group < num_local_groups; ++group) {
105 double *ell_plus = ell_plus_ptr;
106
107 for (int d = 0; d < num_local_directions; d++) {
108 double * KRESTRICT phi_out = phi_out_ptr;
109
110 for(int nm_offset = 0;nm_offset < nidx;++nm_offset){
111 double ell_plus_d_nm = ell_plus[nm_offset];
112
113 #ifdef KRIPKE_USE_OPENMP
114 #pragma omp parallel for
115 #endif
116 for(int z = 0;z < num_zones; ++ z){
117 rhs[z] += ell_plus_d_nm * phi_out[z];
118 }
119 phi_out += num_zones;
120 }
121 ell_plus += nidx;
122 rhs += num_zones;
123 }
124 phi_out_ptr += num_zones*nidx;
125 }
126 } // Subdomain
127 }
128
129 /**
130 Compute scattering source term phi_out from flux moments in phi.
131 phi_out(gp,z,nm) = sum_g { sigs(g, n, gp) * phi(g,z,nm) }
132
133 we are mapping sigs(g,d,z) to mean:
134 g=source group
135 d=legendre coeff
136 z=destination group
137 */
scattering(Grid_Data * grid_data)138 void Kernel_3d_GDZ::scattering(Grid_Data *grid_data){
139 // Loop over zoneset subdomains
140 for(int zs = 0;zs < grid_data->num_zone_sets;++ zs){
141 // get the phi and phi out references
142 SubTVec &phi = *grid_data->phi[zs];
143 SubTVec &phi_out = *grid_data->phi_out[zs];
144 SubTVec &sigs0 = *grid_data->sigs[0];
145 SubTVec &sigs1 = *grid_data->sigs[1];
146 SubTVec &sigs2 = *grid_data->sigs[2];
147
148 // get material mix information
149 int sdom_id = grid_data->zs_to_sdomid[zs];
150 Subdomain &sdom = grid_data->subdomains[sdom_id];
151 int const * KRESTRICT mixed_to_zones = &sdom.mixed_to_zones[0];
152 int const * KRESTRICT mixed_material = &sdom.mixed_material[0];
153 double const * KRESTRICT mixed_fraction = &sdom.mixed_fraction[0];
154
155 // Zero out source terms
156 phi_out.clear(0.0);
157
158 // grab dimensions
159 int num_mixed = sdom.mixed_to_zones.size();
160 int num_zones = sdom.num_zones;
161 int num_groups = phi.groups;
162 int num_moments = grid_data->total_num_moments;
163 int num_coeff = grid_data->legendre_order+1;
164 int const * KRESTRICT moment_to_coeff = &grid_data->moment_to_coeff[0];
165
166 double *phi_g = phi.ptr();
167 double *sigs0_g_gp = sigs0.ptr();
168 double *sigs1_g_gp = sigs1.ptr();
169 double *sigs2_g_gp = sigs2.ptr();
170 for(int g = 0;g < num_groups;++ g){
171
172 double *phi_out_gp_nm = phi_out.ptr();
173 for(int gp = 0;gp < num_groups;++ gp){
174
175 double *phi_g_nm = phi_g;
176 for(int nm = 0;nm < num_moments;++ nm){
177 // map nm to n
178 int n = moment_to_coeff[nm];
179 double sigs_g_gp_n[3] = {sigs0_g_gp[n], sigs1_g_gp[n], sigs2_g_gp[n]};
180
181 for(int mix = 0;mix < num_mixed;++ mix){
182 int zone = mixed_to_zones[mix];
183 int material = mixed_material[mix];
184 double fraction = mixed_fraction[mix];
185 double sigs_value = sigs_g_gp_n[material];
186
187 phi_out_gp_nm[zone] += sigs_value * phi_g_nm[zone] * fraction;
188 }
189
190 phi_g_nm += num_zones;
191 phi_out_gp_nm += num_zones;
192 }
193 sigs0_g_gp += num_coeff;
194 sigs1_g_gp += num_coeff;
195 sigs2_g_gp += num_coeff;
196 }
197 phi_g += num_moments*num_zones;
198 }
199 }
200 }
201
202
203 /**
204 * Add an isotropic source, with flux of 1, to every zone with Region 1
205 * (or material 0).
206 *
207 * Since it's isotropic, we're just adding this to nm=0.
208 */
source(Grid_Data * grid_data)209 void Kernel_3d_GDZ::source(Grid_Data *grid_data){
210 // Loop over zoneset subdomains
211 for(int zs = 0;zs < grid_data->num_zone_sets;++ zs){
212 // get the phi and phi out references
213 SubTVec &phi_out = *grid_data->phi_out[zs];
214
215 // get material mix information
216 int sdom_id = grid_data->zs_to_sdomid[zs];
217 Subdomain &sdom = grid_data->subdomains[sdom_id];
218 int const * KRESTRICT mixed_to_zones = &sdom.mixed_to_zones[0];
219 int const * KRESTRICT mixed_material = &sdom.mixed_material[0];
220 double const * KRESTRICT mixed_fraction = &sdom.mixed_fraction[0];
221
222 // grab dimensions
223 int num_mixed = sdom.mixed_to_zones.size();
224 int num_zones = sdom.num_zones;
225 int num_groups = phi_out.groups;
226 int num_moments = grid_data->total_num_moments;
227
228 double *phi_out_g_nm0 = phi_out.ptr();
229 for(int g = 0;g < num_groups;++ g){
230 for(int mix = 0;mix < num_mixed;++ mix){
231 int zone = mixed_to_zones[mix];
232 int material = mixed_material[mix];
233 double fraction = mixed_fraction[mix];
234
235 if(material == 0){
236 phi_out_g_nm0[zone] += 1.0 * fraction;
237 }
238 }
239 phi_out_g_nm0 += num_moments * num_zones;
240 }
241 }
242 }
243
244 /* Sweep routine for Diamond-Difference */
245 /* Macros for offsets with fluxes on cell faces */
246 #define I_PLANE_INDEX(j, k) (k)*(local_jmax) + (j)
247 #define J_PLANE_INDEX(i, k) (k)*(local_imax) + (i)
248 #define K_PLANE_INDEX(i, j) (j)*(local_imax) + (i)
249 #define Zonal_INDEX(i, j, k) (i) + (local_imax)*(j) \
250 + (local_imax)*(local_jmax)*(k)
251
sweep(Subdomain * sdom)252 void Kernel_3d_GDZ::sweep(Subdomain *sdom) {
253 int num_directions = sdom->num_directions;
254 int num_groups = sdom->num_groups;
255 int num_zones = sdom->num_zones;
256
257 Directions *direction = sdom->directions;
258
259 int local_imax = sdom->nzones[0];
260 int local_jmax = sdom->nzones[1];
261 int local_kmax = sdom->nzones[2];
262
263 double * dx = &sdom->deltas[0][0];
264 double * dy = &sdom->deltas[1][0];
265 double * dz = &sdom->deltas[2][0];
266
267 // Upwind/Downwind face flux data
268 SubTVec &i_plane = *sdom->plane_data[0];
269 SubTVec &j_plane = *sdom->plane_data[1];
270 SubTVec &k_plane = *sdom->plane_data[2];
271
272 // All directions have same id,jd,kd, since these are all one Direction Set
273 // So pull that information out now
274 Grid_Sweep_Block const &extent = sdom->sweep_block;
275
276 #ifdef KRIPKE_USE_OPENMP
277 #pragma omp parallel for
278 #endif
279 for (int group = 0; group < num_groups; ++group) {
280
281 std::vector<double> xcos_dxi_all(local_imax);
282 std::vector<double> ycos_dyj_all(local_jmax);
283 std::vector<double> zcos_dzk_all(local_kmax);
284 double * KRESTRICT sigt_g = sdom->sigt->ptr(group, 0, 0);
285
286 for (int d = 0; d < num_directions; ++d) {
287 double * KRESTRICT psi_g_d = sdom->psi->ptr(group, d, 0);
288 double * KRESTRICT rhs_g_d = sdom->rhs->ptr(group, d, 0);
289 double * KRESTRICT i_plane_g_d = i_plane.ptr(group, d, 0);
290 double * KRESTRICT j_plane_g_d = j_plane.ptr(group, d, 0);
291 double * KRESTRICT k_plane_g_d = k_plane.ptr(group, d, 0);
292
293 double xcos = direction[d].xcos;
294 double ycos = direction[d].ycos;
295 double zcos = direction[d].zcos;
296 for (int i = 0; i < local_imax; ++i) {
297 double dxi = dx[i + 1];
298 xcos_dxi_all[i] = 2.0 * xcos / dxi;
299 }
300
301 for (int j = 0; j < local_jmax; ++j) {
302 double dyj = dy[j + 1];
303 ycos_dyj_all[j] = 2.0 * ycos / dyj;
304 }
305
306 for (int k = 0; k < local_kmax; ++k) {
307 double dzk = dz[k + 1];
308 zcos_dzk_all[k] = 2.0 * zcos / dzk;
309 }
310
311 /* Perform transport sweep of the grid 1 cell at a time. */
312 for (int k = extent.start_k; k != extent.end_k; k += extent.inc_k) {
313 double zcos_dzk = zcos_dzk_all[k];
314 for (int j = extent.start_j; j != extent.end_j; j += extent.inc_j) {
315 double ycos_dyj = ycos_dyj_all[j];
316 int z_idx = Zonal_INDEX(extent.start_i, j, k);
317 for (int i = extent.start_i; i != extent.end_i; i += extent.inc_i) {
318 double xcos_dxi = xcos_dxi_all[i];
319
320 /* Calculate new zonal flux */
321 double psi_g_d_z = (rhs_g_d[z_idx]
322 + i_plane_g_d[I_PLANE_INDEX(j, k)] * xcos_dxi
323 + j_plane_g_d[J_PLANE_INDEX(i, k)] * ycos_dyj
324 + k_plane_g_d[K_PLANE_INDEX(i, j)] * zcos_dzk)
325 / (xcos_dxi + ycos_dyj + zcos_dzk
326 + sigt_g[z_idx]);
327 psi_g_d[z_idx] = psi_g_d_z;
328
329 /* Apply diamond-difference relationships */
330 i_plane_g_d[I_PLANE_INDEX(j, k)] = 2.0 * psi_g_d_z
331 - i_plane_g_d[I_PLANE_INDEX(j, k)];
332 j_plane_g_d[J_PLANE_INDEX(i, k)] = 2.0 * psi_g_d_z
333 - j_plane_g_d[J_PLANE_INDEX(i, k)];
334 k_plane_g_d[K_PLANE_INDEX(i, j)] = 2.0 * psi_g_d_z
335 - k_plane_g_d[K_PLANE_INDEX(i, j)];
336
337 z_idx += extent.inc_i;
338 }
339 }
340 }
341 }
342 }
343
344 }
345
346