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