1 #include<Kripke/Kernel/Kernel_3d_GZD.h>
2 #include<Kripke/Grid.h>
3 #include<Kripke/SubTVec.h>
4 
nestingPsi(void) const5 Nesting_Order Kernel_3d_GZD::nestingPsi(void) const {
6   return NEST_GZD;
7 }
8 
nestingPhi(void) const9 Nesting_Order Kernel_3d_GZD::nestingPhi(void) const {
10   return NEST_GZD;
11 }
12 
nestingSigt(void) const13 Nesting_Order Kernel_3d_GZD::nestingSigt(void) const {
14   return NEST_DGZ;
15 }
16 
nestingEll(void) const17 Nesting_Order Kernel_3d_GZD::nestingEll(void) const {
18   return NEST_ZDG;
19 }
20 
nestingEllPlus(void) const21 Nesting_Order Kernel_3d_GZD::nestingEllPlus(void) const {
22   return NEST_ZDG;
23 }
24 
nestingSigs(void) const25 Nesting_Order Kernel_3d_GZD::nestingSigs(void) const {
26   return NEST_GZD;
27 }
28 
29 
LTimes(Grid_Data * grid_data)30 void Kernel_3d_GZD::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     int num_groups_zones = num_local_groups*num_zones;
50 
51     /* 3D Cartesian Geometry */
52     double *ell_ptr = sdom.ell->ptr();
53 
54 #ifdef KRIPKE_USE_OPENMP
55 #pragma omp parallel for
56 #endif
57     for(int gz = 0;gz < num_groups_zones; ++ gz){
58       double * KRESTRICT psi = sdom.psi->ptr() + gz*num_local_directions;
59       double * KRESTRICT phi = sdom.phi->ptr(group0, 0, 0) + gz*nidx;
60       double * KRESTRICT ell_d = ell_ptr;
61 
62       for (int d = 0; d < num_local_directions; d++) {
63         double psi_d = psi[d];
64 
65         for(int nm_offset = 0;nm_offset < nidx;++nm_offset){
66           phi[nm_offset] += ell_d[nm_offset] * psi_d;
67         }
68         ell_d += nidx;
69       }
70 
71     }
72   } // Subdomain
73 }
74 
LPlusTimes(Grid_Data * grid_data)75 void Kernel_3d_GZD::LPlusTimes(Grid_Data *grid_data) {
76   // Outer parameters
77   int nidx = grid_data->total_num_moments;
78 
79   // Loop over Subdomains
80   int num_subdomains = grid_data->subdomains.size();
81   for (int sdom_id = 0; sdom_id < num_subdomains; ++ sdom_id){
82     Subdomain &sdom = grid_data->subdomains[sdom_id];
83 
84     // Get dimensioning
85     int num_zones = sdom.num_zones;
86     int num_local_groups = sdom.num_groups;
87     int group0 = sdom.group0;
88     int num_local_directions = sdom.num_directions;
89     int num_groups_zones = num_local_groups*num_zones;
90 
91     sdom.rhs->clear(0.0);
92 
93     /* 3D Cartesian Geometry */
94     double * KRESTRICT ell_plus_ptr = sdom.ell_plus->ptr();
95 
96 #ifdef KRIPKE_USE_OPENMP
97 #pragma omp parallel for
98 #endif
99     for(int gz = 0;gz < num_groups_zones; ++ gz){
100       double * KRESTRICT rhs = sdom.rhs->ptr(0, 0, 0) + gz*num_local_directions;
101       double * KRESTRICT phi_out = sdom.phi_out->ptr(group0, 0, 0) + gz*nidx;
102       double * KRESTRICT ell_plus_d = ell_plus_ptr;
103 
104       for (int d = 0; d < num_local_directions; d++) {
105 
106         for(int nm_offset = 0;nm_offset < nidx;++nm_offset){
107           rhs[d] += ell_plus_d[nm_offset] * phi_out[nm_offset];
108         }
109         ell_plus_d += nidx;
110       }
111     }
112   } // Subdomain
113 }
114 
115 
116 /**
117   Compute scattering source term phi_out from flux moments in phi.
118   phi_out(gp,z,nm) = sum_g { sigs(g, n, gp) * phi(g,z,nm) }
119 
120   we are mapping sigs(g,d,z) to mean:
121     g=source group
122     d=legendre coeff
123     z=destination group
124 */
scattering(Grid_Data * grid_data)125 void Kernel_3d_GZD::scattering(Grid_Data *grid_data){
126   // Loop over zoneset subdomains
127   for(int zs = 0;zs < grid_data->num_zone_sets;++ zs){
128     // get the phi and phi out references
129     SubTVec &phi = *grid_data->phi[zs];
130     SubTVec &phi_out = *grid_data->phi_out[zs];
131     SubTVec &sigs0 = *grid_data->sigs[0];
132     SubTVec &sigs1 = *grid_data->sigs[1];
133     SubTVec &sigs2 = *grid_data->sigs[2];
134 
135     // get material mix information
136     int sdom_id = grid_data->zs_to_sdomid[zs];
137     Subdomain &sdom = grid_data->subdomains[sdom_id];
138     int const * KRESTRICT mixed_to_zones = &sdom.mixed_to_zones[0];
139     int const * KRESTRICT mixed_material = &sdom.mixed_material[0];
140     double const * KRESTRICT mixed_fraction = &sdom.mixed_fraction[0];
141 
142     // Zero out source terms
143     phi_out.clear(0.0);
144 
145     // grab dimensions
146     int num_mixed = sdom.mixed_to_zones.size();
147     int num_zones = sdom.num_zones;
148     int num_groups = phi.groups;
149     int num_coeff = grid_data->legendre_order+1;
150     int num_moments = grid_data->total_num_moments;
151     int const * KRESTRICT moment_to_coeff = &grid_data->moment_to_coeff[0];
152 
153     double *phi_g = phi.ptr();
154     double *sigs0_g_gp = sigs0.ptr();
155     double *sigs1_g_gp = sigs1.ptr();
156     double *sigs2_g_gp = sigs2.ptr();
157     for(int g = 0;g < num_groups;++ g){
158 
159       double *phi_out_gp = phi_out.ptr();
160       for(int gp = 0;gp < num_groups;++ gp){
161 
162         double *sigs_g_gp[3] = {
163           sigs0_g_gp,
164           sigs1_g_gp,
165           sigs2_g_gp
166         };
167 
168         for(int mix = 0;mix < num_mixed;++ mix){
169           int zone = mixed_to_zones[mix];
170           int material = mixed_material[mix];
171           double fraction = mixed_fraction[mix];
172           double *sigs_g_gp_mat = sigs_g_gp[material];
173           double *phi_g_z = phi_g + zone*num_moments;
174           double *phi_out_gp_z = phi_out_gp + zone*num_moments;
175 
176           for(int nm = 0;nm < num_moments;++ nm){
177             // map nm to n
178             int n = moment_to_coeff[nm];
179 
180             phi_out_gp_z[nm] += sigs_g_gp_mat[n] * phi_g_z[nm] * fraction;
181           }
182         }
183         sigs0_g_gp += num_coeff;
184         sigs1_g_gp += num_coeff;
185         sigs2_g_gp += num_coeff;
186         phi_out_gp += num_zones*num_moments;
187       }
188       phi_g += num_zones*num_moments;
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_GZD::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_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_g[zone*num_moments] += 1.0 * fraction;
228         }
229       }
230       phi_out_g += num_moments * num_zones;
231     }
232   }
233 }
234 
235 /* Sweep routine for Diamond-Difference */
236 /* Macros for offsets with fluxes on cell faces */
237 #define I_PLANE_INDEX(j, k) (k)*(local_jmax) + (j)
238 #define J_PLANE_INDEX(i, k) (k)*(local_imax) + (i)
239 #define K_PLANE_INDEX(i, j) (j)*(local_imax) + (i)
240 #define Zonal_INDEX(i, j, k) (i) + (local_imax)*(j) \
241   + (local_imax)*(local_jmax)*(k)
242 
sweep(Subdomain * sdom)243 void Kernel_3d_GZD::sweep(Subdomain *sdom) {
244   int num_directions = sdom->num_directions;
245   int num_groups = sdom->num_groups;
246   int num_zones = sdom->num_zones;
247 
248   Directions *direction = sdom->directions;
249 
250   int local_imax = sdom->nzones[0];
251   int local_jmax = sdom->nzones[1];
252   int local_kmax = sdom->nzones[2];
253 
254   double *dx = &sdom->deltas[0][0];
255   double *dy = &sdom->deltas[1][0];
256   double *dz = &sdom->deltas[2][0];
257 
258   // Upwind/Downwind face flux data
259   SubTVec &i_plane = *sdom->plane_data[0];
260   SubTVec &j_plane = *sdom->plane_data[1];
261   SubTVec &k_plane = *sdom->plane_data[2];
262 
263   // All directions have same id,jd,kd, since these are all one Direction Set
264   // So pull that information out now
265   Grid_Sweep_Block const &extent = sdom->sweep_block;
266 
267 #ifdef KRIPKE_USE_OPENMP
268 #pragma omp parallel for
269 #endif
270   for (int group = 0; group < num_groups; ++group) {
271     double *sigt_g = sdom->sigt->ptr(group, 0, 0);
272 
273     /*  Perform transport sweep of the grid 1 cell at a time.   */
274     for (int k = extent.start_k; k != extent.end_k; k += extent.inc_k) {
275       double dzk = dz[k + 1];
276       double two_dz = 2.0 / dzk;
277       for (int j = extent.start_j; j != extent.end_j; j += extent.inc_j) {
278         double dyj = dy[j + 1];
279         double two_dy = 2.0 / dyj;
280         for (int i = extent.start_i; i != extent.end_i; i += extent.inc_i) {
281           double dxi = dx[i + 1];
282           double two_dx = 2.0 / dxi;
283 
284           int z = Zonal_INDEX(i, j, k);
285 
286           double * KRESTRICT psi_g_z = sdom->psi->ptr(group, 0, z);
287           double * KRESTRICT rhs_g_z = sdom->rhs->ptr(group, 0, z);
288 
289           double * KRESTRICT psi_lf_g_z = i_plane.ptr(group, 0, I_PLANE_INDEX(j, k));
290           double * KRESTRICT psi_fr_g_z = j_plane.ptr(group, 0, J_PLANE_INDEX(i, k));
291           double * KRESTRICT psi_bo_g_z = k_plane.ptr(group, 0, K_PLANE_INDEX(i, j));
292 
293           for (int d = 0; d < num_directions; ++d) {
294             double xcos = direction[d].xcos;
295             double ycos = direction[d].ycos;
296             double zcos = direction[d].zcos;
297 
298             double zcos_dzk = zcos * two_dz;
299             double ycos_dyj = ycos * two_dy;
300             double xcos_dxi = xcos * two_dx;
301 
302             /* Calculate new zonal flux */
303             double psi_g_z_d = (rhs_g_z[d] + psi_lf_g_z[d] * xcos_dxi
304                 + psi_fr_g_z[d] * ycos_dyj + psi_bo_g_z[d] * zcos_dzk)
305                 / (xcos_dxi + ycos_dyj + zcos_dzk
306                     + sigt_g[Zonal_INDEX(i, j, k)]);
307 
308             psi_g_z[d] = psi_g_z_d;
309 
310             /* Apply diamond-difference relationships */
311             psi_lf_g_z[d] = 2.0 * psi_g_z_d - psi_lf_g_z[d];
312             psi_fr_g_z[d] = 2.0 * psi_g_z_d - psi_fr_g_z[d];
313             psi_bo_g_z[d] = 2.0 * psi_g_z_d - psi_bo_g_z[d];
314           }
315         }
316       }
317     }
318   } // group
319 }
320 
321 
322