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