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