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