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