source: CIVL/examples/omp/HydroC/make_boundary.c@ 1aaefd4

main test-branch
Last change on this file since 1aaefd4 was ea777aa, checked in by Alex Wilton <awilton@…>, 3 years ago

Moved examples, include, build_default.properties, common.xml, and README out from dev.civl.com into the root of the repo.

git-svn-id: ​svn://vsl.cis.udel.edu/civl/trunk@5704 fb995dde-84ed-4084-dfe6-e5aef3e2452c
  • Property mode set to 100644
File size: 12.9 KB
Line 
1/*
2 A simple 2D hydro code
3 (C) Romain Teyssier : CEA/IRFU -- original F90 code
4 (C) Pierre-Francois Lavallee : IDRIS -- original F90 code
5 (C) Guillaume Colin de Verdiere : CEA/DAM -- for the C version
6 (C) Adèle Villiermet : CINES -- for FTI integration
7*/
8/*
9
10This software is governed by the CeCILL license under French law and
11abiding by the rules of distribution of free software. You can use,
12modify and/ or redistribute the software under the terms of the CeCILL
13license as circulated by CEA, CNRS and INRIA at the following URL
14"http://www.cecill.info".
15
16As a counterpart to the access to the source code and rights to copy,
17modify and redistribute granted by the license, users are provided only
18with a limited warranty and the software's author, the holder of the
19economic rights, and the successive licensors have only limited
20liability.
21
22In this respect, the user's attention is drawn to the risks associated
23with loading, using, modifying and/or developing or reproducing the
24software by the user in light of its specific status of free software,
25that may mean that it is complicated to manipulate, and that also
26therefore means that it is reserved for developers and experienced
27professionals having in-depth computer knowledge. Users are therefore
28encouraged to load and test the software's suitability as regards their
29requirements in conditions enabling the security of their systems and/or
30data to be ensured and, more generally, to use and operate it in the
31same conditions as regards security.
32
33The fact that you are presently reading this means that you have had
34knowledge of the CeCILL license and that you accept its terms.
35
36*/
37
38#ifdef MPI
39#include <mpi.h>
40#if FTI>0
41#include <fti.h>
42#endif
43#endif
44#include <stdlib.h>
45#include <unistd.h>
46#include <math.h>
47#include <stdio.h>
48#include <string.h>
49#include <strings.h>
50#include <assert.h>
51
52
53#include "parametres.h"
54#include "make_boundary.h"
55#include "perfcnt.h"
56#include "utils.h"
57
58static int
59pack_arrayv(const int xmin, const hydroparam_t H, hydrovar_t * Hv, real_t *buffer);
60static int
61unpack_arrayv(const int xmin, const hydroparam_t H, hydrovar_t * Hv, real_t *buffer);
62static int
63pack_arrayh(const int xmin, const hydroparam_t H, hydrovar_t * Hv, real_t *buffer);
64static int
65unpack_arrayh(const int xmin, const hydroparam_t H, hydrovar_t * Hv, real_t *buffer);
66
67int
68pack_arrayv(const int xmin, const hydroparam_t H, hydrovar_t * Hv, real_t *buffer) {
69 int ivar, i, j, p = 0;
70 for (ivar = 0; ivar < H.nvar; ivar++) {
71 for (j = 0; j < H.nyt; j++) {
72 // #warning "GATHER to vectorize ?"
73 for (i = xmin; i < xmin + ExtraLayer; i++) {
74 buffer[p++] = Hv->uold[IHv(i, j, ivar)];
75 }
76 }
77 }
78 return p;
79}
80
81int
82unpack_arrayv(const int xmin, const hydroparam_t H, hydrovar_t * Hv, real_t *buffer) {
83 int ivar, i, j, p = 0;
84 for (ivar = 0; ivar < H.nvar; ivar++) {
85 for (j = 0; j < H.nyt; j++) {
86 // #warning "SCATTER to vectorize ?"
87 for (i = xmin; i < xmin + ExtraLayer; i++) {
88 Hv->uold[IHv(i, j, ivar)] = buffer[p++];
89 }
90 }
91 }
92 return p;
93}
94
95int
96pack_arrayh(const int ymin, const hydroparam_t H, hydrovar_t * Hv, real_t *buffer) {
97 int ivar, i, j, p = 0;
98 for (ivar = 0; ivar < H.nvar; ivar++) {
99 for (j = ymin; j < ymin + ExtraLayer; j++) {
100 // #warning "GATHER to vectorize ?"
101 // #pragma simd
102 for (i = 0; i < H.nxt; i++) {
103 buffer[p++] = Hv->uold[IHv(i, j, ivar)];
104 }
105 }
106 }
107 return p;
108}
109
110int
111unpack_arrayh(const int ymin, const hydroparam_t H, hydrovar_t * Hv, real_t *buffer) {
112 int ivar, i, j, p = 0;
113 for (ivar = 0; ivar < H.nvar; ivar++) {
114 for (j = ymin; j < ymin + ExtraLayer; j++) {
115 // #warning "SCATTER to vectorize ?"
116 for (i = 0; i < H.nxt; i++) {
117 Hv->uold[IHv(i, j, ivar)] = buffer[p++];
118 }
119 }
120 }
121 return p;
122}
123
124#define VALPERLINE 11
125int
126print_bufferh(FILE * fic, const int ymin, const hydroparam_t H, hydrovar_t * Hv, real_t *buffer) {
127 int ivar, i, j, p = 0, nbr = 1;
128 for (ivar = 3; ivar < H.nvar; ivar++) {
129 fprintf(fic, "BufferH v=%d\n", ivar);
130 for (j = ymin; j < ymin + ExtraLayer; j++) {
131 for (i = 0; i < H.nxt; i++) {
132 fprintf(fic, "%13.6le ", buffer[p++]);
133 nbr++;
134 if (nbr == VALPERLINE) {
135 fprintf(fic, "\n");
136 nbr = 1;
137 }
138 }
139 }
140 if (nbr != 1)
141 fprintf(fic, "\n");
142 }
143 return p;
144}
145
146void
147make_boundary(int idim, const hydroparam_t H, hydrovar_t * Hv) {
148
149 // - - - - - - - - - - - - - - - - - - -
150 // Cette portion de code est à vérifier
151 // détail. J'ai des doutes sur la conversion
152 // des index depuis fortran.
153 // - - - - - - - - - - - - - - - - - - -
154 int i, ivar, i0, j, j0, err, size;
155 real_t sign;
156 real_t sendbufld[ExtraLayerTot * H.nxyt * H.nvar];
157 real_t sendbufru[ExtraLayerTot * H.nxyt * H.nvar];
158 // real_t *sendbufru, *sendbufld;
159 real_t recvbufru[ExtraLayerTot * H.nxyt * H.nvar];
160 real_t recvbufld[ExtraLayerTot * H.nxyt * H.nvar];
161 // real_t *recvbufru, *recvbufld;
162#ifdef MPI
163 MPI_Request requests[4];
164 MPI_Status status[4];
165 MPI_Datatype mpiFormat = MPI_DOUBLE;
166#endif
167 int reqcnt = 0;
168
169 static FILE *fic = NULL;
170
171 WHERE("make_boundary");
172
173#ifdef MPI
174 if (sizeof(real_t) == sizeof(float)) mpiFormat = MPI_FLOAT;
175#endif
176
177 if (idim == 1) {
178#ifdef MPI
179 i = ExtraLayer;
180 size = pack_arrayv(i, H, Hv, sendbufld);
181 i = H.nx;
182 size = pack_arrayv(i, H, Hv, sendbufru);
183#if FTI==0
184 if (H.box[RIGHT_BOX] != -1) {
185 MPI_Isend(sendbufru, size, mpiFormat, H.box[RIGHT_BOX], 123, MPI_COMM_WORLD, &requests[reqcnt]);
186 reqcnt++;
187 }
188 if (H.box[LEFT_BOX] != -1) {
189 MPI_Isend(sendbufld, size, mpiFormat, H.box[LEFT_BOX], 246, MPI_COMM_WORLD, &requests[reqcnt]);
190 reqcnt++;
191 }
192 if (H.box[RIGHT_BOX] != -1) {
193 MPI_Irecv(recvbufru, size, mpiFormat, H.box[RIGHT_BOX], 246, MPI_COMM_WORLD, &requests[reqcnt]);
194 reqcnt++;
195 }
196 if (H.box[LEFT_BOX] != -1) {
197 MPI_Irecv(recvbufld, size, mpiFormat, H.box[LEFT_BOX], 123, MPI_COMM_WORLD, &requests[reqcnt]);
198 reqcnt++;
199 }
200#endif
201#if FTI>0
202 if (H.box[RIGHT_BOX] != -1) {
203 MPI_Isend(sendbufru, size, mpiFormat, H.box[RIGHT_BOX], 123, FTI_COMM_WORLD, &requests[reqcnt]);
204 reqcnt++;
205 }
206 if (H.box[LEFT_BOX] != -1) {
207 MPI_Isend(sendbufld, size, mpiFormat, H.box[LEFT_BOX], 246, FTI_COMM_WORLD, &requests[reqcnt]);
208 reqcnt++;
209 }
210 if (H.box[RIGHT_BOX] != -1) {
211 MPI_Irecv(recvbufru, size, mpiFormat, H.box[RIGHT_BOX], 246, FTI_COMM_WORLD, &requests[reqcnt]);
212 reqcnt++;
213 }
214 if (H.box[LEFT_BOX] != -1) {
215 MPI_Irecv(recvbufld, size, mpiFormat, H.box[LEFT_BOX], 123, FTI_COMM_WORLD, &requests[reqcnt]);
216 reqcnt++;
217 }
218#endif
219 err = MPI_Waitall(reqcnt, requests, status);
220 assert(err == MPI_SUCCESS);
221
222 if (H.box[RIGHT_BOX] != -1) {
223 {
224 i = H.nx + ExtraLayer;
225 size = unpack_arrayv(i, H, Hv, recvbufru);
226 }
227 }
228
229 if (H.box[LEFT_BOX] != -1) {
230 {
231 i = 0;
232 size = unpack_arrayv(i, H, Hv, recvbufld);
233 }
234 }
235#endif
236
237 if (H.boundary_left > 0) {
238 // Left boundary
239 for (ivar = 0; ivar < H.nvar; ivar++) {
240 for (i = 0; i < ExtraLayer; i++) {
241 sign = 1.0;
242 if (H.boundary_left == 1) {
243 i0 = ExtraLayerTot - i - 1;
244 if (ivar == IU) {
245 sign = -1.0;
246 }
247 } else if (H.boundary_left == 2) {
248 i0 = 2;
249 } else {
250 i0 = H.nx + i;
251 }
252// #pragma simd
253 for (j = H.jmin + ExtraLayer; j < H.jmax - ExtraLayer; j++) {
254 Hv->uold[IHv(i, j, ivar)] = Hv->uold[IHv(i0, j, ivar)] * sign;
255 }
256 }
257 }
258 {
259 int nops = H.nvar * ExtraLayer * ((H.jmax - ExtraLayer) - (H.jmin + ExtraLayer));
260 FLOPS(1 * nops, 0 * nops, 0 * nops, 0 * nops);
261 }
262 }
263
264 if (H.boundary_right > 0) {
265 // Right boundary
266 for (ivar = 0; ivar < H.nvar; ivar++) {
267 for (i = H.nx + ExtraLayer; i < H.nx + ExtraLayerTot; i++) {
268 sign = 1.0;
269 if (H.boundary_right == 1) {
270 i0 = 2 * H.nx + ExtraLayerTot - i - 1;
271 if (ivar == IU) {
272 sign = -1.0;
273 }
274 } else if (H.boundary_right == 2) {
275 i0 = H.nx + ExtraLayer;
276 } else {
277 i0 = i - H.nx;
278 }
279// #pragma simd
280 for (j = H.jmin + ExtraLayer; j < H.jmax - ExtraLayer; j++) {
281 Hv->uold[IHv(i, j, ivar)] = Hv->uold[IHv(i0, j, ivar)] * sign;
282 } // for j
283 } // for i
284 }
285 {
286 int nops = H.nvar * ((H.jmax - ExtraLayer) - (H.jmin + ExtraLayer)) * ((H.nx + ExtraLayerTot) - (H.nx + ExtraLayer));
287 FLOPS(1 * nops, 0 * nops, 0 * nops, 0 * nops);
288 }
289 }
290 } else {
291#ifdef MPI
292 {
293 if (fic) {
294 fprintf(fic, "- = - = - = - Avant\n");
295 printuoldf(fic, H, Hv);
296 }
297 }
298 j = ExtraLayer;
299 size = pack_arrayh(j, H, Hv, sendbufld);
300 // fprintf(stderr, "%d prep %d\n", H.mype, j);
301 if (fic) {
302 fprintf(fic, "%d prep %d\n", H.mype, j);
303 print_bufferh(fic, j, H, Hv, sendbufld);
304 }
305 j = H.ny;
306 size = pack_arrayh(j, H, Hv, sendbufru);
307 // fprintf(stderr, "%d prep %d (s=%d)\n", H.mype, j, size);
308 if (fic) {
309 fprintf(fic, "%d prep %d\n", H.mype, j);
310 print_bufferh(fic, j, H, Hv, sendbufru);
311 }
312#if FTI==0
313 if (H.box[DOWN_BOX] != -1) {
314 MPI_Isend(sendbufld, size, mpiFormat, H.box[DOWN_BOX], 123, MPI_COMM_WORLD, &requests[reqcnt]);
315 reqcnt++;
316 }
317 if (H.box[UP_BOX] != -1) {
318 MPI_Isend(sendbufru, size, mpiFormat, H.box[UP_BOX], 246, MPI_COMM_WORLD, &requests[reqcnt]);
319 reqcnt++;
320 }
321 if (H.box[DOWN_BOX] != -1) {
322 MPI_Irecv(recvbufld, size, mpiFormat, H.box[DOWN_BOX], 246, MPI_COMM_WORLD, &requests[reqcnt]);
323 reqcnt++;
324 }
325 if (H.box[UP_BOX] != -1) {
326 MPI_Irecv(recvbufru, size, mpiFormat, H.box[UP_BOX], 123, MPI_COMM_WORLD, &requests[reqcnt]);
327 reqcnt++;
328 }
329#endif
330#if FTI>0
331 if (H.box[DOWN_BOX] != -1) {
332 MPI_Isend(sendbufld, size, mpiFormat, H.box[DOWN_BOX], 123, FTI_COMM_WORLD, &requests[reqcnt]);
333 reqcnt++;
334 }
335 if (H.box[UP_BOX] != -1) {
336 MPI_Isend(sendbufru, size, mpiFormat, H.box[UP_BOX], 246, FTI_COMM_WORLD, &requests[reqcnt]);
337 reqcnt++;
338 }
339 if (H.box[DOWN_BOX] != -1) {
340 MPI_Irecv(recvbufld, size, mpiFormat, H.box[DOWN_BOX], 246, FTI_COMM_WORLD, &requests[reqcnt]);
341 reqcnt++;
342 }
343 if (H.box[UP_BOX] != -1) {
344 MPI_Irecv(recvbufru, size, mpiFormat, H.box[UP_BOX], 123, FTI_COMM_WORLD, &requests[reqcnt]);
345 reqcnt++;
346 }
347#endif
348 err = MPI_Waitall(reqcnt, requests, status);
349 assert(err == MPI_SUCCESS);
350
351 if (H.box[DOWN_BOX] != -1) {
352 {
353 j = 0;
354 unpack_arrayh(j, H, Hv, recvbufld);
355 if (fic) {
356 fprintf(fic, "%d down %d\n", H.mype, j);
357 print_bufferh(fic, j, H, Hv, recvbufld);
358 }
359 // fprintf(stderr, "%d down %d\n", H.mype, j);
360 }
361 }
362 if (H.box[UP_BOX] != -1) {
363 {
364 j = H.ny + ExtraLayer;
365 unpack_arrayh(j, H, Hv, recvbufru);
366 if (fic) {
367 fprintf(fic, "%d up %d\n", H.mype, j);
368 print_bufferh(fic, j, H, Hv, recvbufru);
369 }
370 // fprintf(stderr, "%d up %d\n", H.mype, j);
371 }
372 }
373 // if (H.mype == 0)
374 {
375 if (fic) {
376 fprintf(fic, "- = - = - = - Apres\n");
377 printuoldf(fic, H, Hv);
378 }
379 }
380#endif
381 // Lower boundary
382 if (H.boundary_down > 0) {
383 j0 = 0;
384 for (ivar = 0; ivar < H.nvar; ivar++) {
385 for (j = 0; j < ExtraLayer; j++) {
386 sign = 1.0;
387 if (H.boundary_down == 1) {
388 j0 = ExtraLayerTot - j - 1;
389 if (ivar == IV) {
390 sign = -1.0;
391 }
392 } else if (H.boundary_down == 2) {
393 j0 = ExtraLayerTot;
394 } else {
395 j0 = H.ny + j;
396 }
397// #pragma simd
398 for (i = H.imin + ExtraLayer; i < H.imax - ExtraLayer; i++) {
399 Hv->uold[IHv(i, j, ivar)] = Hv->uold[IHv(i, j0, ivar)] * sign;
400 }
401 }
402 }
403 {
404 int nops = H.nvar * ((ExtraLayer) - (0)) * ((H.imax - ExtraLayer) - (H.imin + ExtraLayer));
405 FLOPS(1 * nops, 0 * nops, 0 * nops, 0 * nops);
406 }
407 }
408 // Upper boundary
409 if (H.boundary_up > 0) {
410 for (ivar = 0; ivar < H.nvar; ivar++) {
411 for (j = H.ny + ExtraLayer; j < H.ny + ExtraLayerTot; j++) {
412 sign = 1.0;
413 if (H.boundary_up == 1) {
414 j0 = 2 * H.ny + ExtraLayerTot - j - 1;
415 if (ivar == IV) {
416 sign = -1.0;
417 }
418 } else if (H.boundary_up == 2) {
419 j0 = H.ny + 1;
420 } else {
421 j0 = j - H.ny;
422 }
423// #pragma simd
424 for (i = H.imin + ExtraLayer; i < H.imax - ExtraLayer; i++) {
425 Hv->uold[IHv(i, j, ivar)] = Hv->uold[IHv(i, j0, ivar)] * sign;
426 }
427 }
428 }
429 {
430 int nops = H.nvar * ((H.ny + ExtraLayerTot) - (H.ny + ExtraLayer)) * ((H.imax - ExtraLayer) - (H.imin + ExtraLayer));
431 FLOPS(1 * nops, 0 * nops, 0 * nops, 0 * nops);
432 }
433 }
434 }
435}
436
437// make_boundary
438//EOF
Note: See TracBrowser for help on using the repository browser.