source: CIVL/examples/omp/HydroC/hydro_funcs.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.7 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*/
7/*
8
9This software is governed by the CeCILL license under French law and
10abiding by the rules of distribution of free software. You can use,
11modify and/ or redistribute the software under the terms of the CeCILL
12license as circulated by CEA, CNRS and INRIA at the following URL
13"http://www.cecill.info".
14
15As a counterpart to the access to the source code and rights to copy,
16modify and redistribute granted by the license, users are provided only
17with a limited warranty and the software's author, the holder of the
18economic rights, and the successive licensors have only limited
19liability.
20
21In this respect, the user's attention is drawn to the risks associated
22with loading, using, modifying and/or developing or reproducing the
23software by the user in light of its specific status of free software,
24that may mean that it is complicated to manipulate, and that also
25therefore means that it is reserved for developers and experienced
26professionals having in-depth computer knowledge. Users are therefore
27encouraged to load and test the software's suitability as regards their
28requirements in conditions enabling the security of their systems and/or
29data to be ensured and, more generally, to use and operate it in the
30same conditions as regards security.
31
32The fact that you are presently reading this means that you have had
33knowledge of the CeCILL license and that you accept its terms.
34
35*/
36
37#include <stdlib.h>
38#include <unistd.h>
39#include <stdint.h>
40#include <math.h>
41#include <stdio.h>
42#include <string.h>
43#include <assert.h>
44
45#include "utils.h"
46#include "hydro_utils.h"
47#include "hydro_funcs.h"
48
49#include "hydro_numa.h"
50
51void
52hydro_init(hydroparam_t * H, hydrovar_t * Hv) {
53 int i, j;
54 int x, y;
55
56 // *WARNING* : we will use 0 based arrays everywhere since it is C code!
57 H->imin = H->jmin = 0;
58
59 // We add two extra layers left/right/top/bottom
60 H->imax = H->nx + ExtraLayerTot;
61 H->jmax = H->ny + ExtraLayerTot;
62 H->nxt = H->imax - H->imin; // column size in the array
63 H->nyt = H->jmax - H->jmin; // row size in the array
64
65 // maximum direction size
66 H->nxyt = (H->nxt > H->nyt) ? H->nxt : H->nyt;
67 // To make sure that slices are properly aligned, we make the array a
68 // multiple of NDBLE double
69#define NDBLE 16
70 // printf("avant %d ", H->nxyt);
71 // H->nxyt = (H->nxyt + NDBLE - 1) / NDBLE;
72 // H->nxyt *= NDBLE;
73 // printf("apres %d \n", H->nxyt);
74
75 // allocate uold for each conservative variable
76 Hv->uold = (real_t *) DMalloc(H->nvar * H->nxt * H->nyt);
77
78 // wind tunnel with point explosion
79 for (j = H->jmin + ExtraLayer; j < H->jmax - ExtraLayer; j++) {
80 for (i = H->imin + ExtraLayer; i < H->imax - ExtraLayer; i++) {
81 Hv->uold[IHvP(i, j, ID)] = one;
82 Hv->uold[IHvP(i, j, IU)] = zero;
83 Hv->uold[IHvP(i, j, IV)] = zero;
84 Hv->uold[IHvP(i, j, IP)] = 1e-5;
85 }
86 }
87
88 // Initial shock
89 if (H->testCase == 0) {
90 if (H->nproc == 1) {
91 x = (H->imax - H->imin) / 2 + ExtraLayer * 0;
92 y = (H->jmax - H->jmin) / 2 + ExtraLayer * 0;
93 Hv->uold[IHvP(x, y, IP)] = one / H->dx / H->dx;
94 printf("Centered test case : %d %d\n", x, y);
95 } else {
96 x = ((H->globnx) / 2);
97 y = ((H->globny) / 2);
98 if ((x >= H->box[XMIN_BOX]) && (x < H->box[XMAX_BOX]) && (y >= H->box[YMIN_BOX]) && (y < H->box[YMAX_BOX])) {
99 x = x - H->box[XMIN_BOX] + ExtraLayer;
100 y = y - H->box[YMIN_BOX] + ExtraLayer;
101 Hv->uold[IHvP(x, y, IP)] = one / H->dx / H->dx;
102 printf("Centered test case : [%d] %d %d\n", H->mype, x, y);
103 }
104 }
105 }
106 if (H->testCase == 1) {
107 if (H->nproc == 1) {
108 x = ExtraLayer;
109 y = ExtraLayer;
110 Hv->uold[IHvP(x, y, IP)] = one / H->dx / H->dx;
111 printf("Lower corner test case : %d %d\n", x, y);
112 } else {
113 x = ExtraLayer;
114 y = ExtraLayer;
115 if ((x >= H->box[XMIN_BOX]) && (x < H->box[XMAX_BOX]) && (y >= H->box[YMIN_BOX]) && (y < H->box[YMAX_BOX])) {
116 Hv->uold[IHvP(x, y, IP)] = one / H->dx / H->dx;
117 printf("Lower corner test case : [%d] %d %d\n", H->mype, x, y);
118 }
119 }
120 }
121 if (H->testCase == 2) {
122 if (H->nproc == 1) {
123 x = ExtraLayer;
124 y = ExtraLayer;
125 for (j = y; j < H->jmax; j++) {
126 Hv->uold[IHvP(x, j, IP)] = one / H->dx / H->dx;
127 }
128 printf("SOD tube test case\n");
129 } else {
130 x = ExtraLayer;
131 y = ExtraLayer;
132 for (j = 0; j < H->globny; j++) {
133 if ((x >= H->box[XMIN_BOX]) && (x < H->box[XMAX_BOX]) && (j >= H->box[YMIN_BOX]) && (j < H->box[YMAX_BOX])) {
134 y = j - H->box[YMIN_BOX] + ExtraLayer;
135 Hv->uold[IHvP(x, y, IP)] = one / H->dx / H->dx;
136 }
137 }
138 printf("SOD tube test case in //\n");
139 }
140 }
141 if (H->testCase > 2) {
142 printf("Test case not implemented -- aborting !\n");
143 abort();
144 }
145 fflush(stdout);
146} // hydro_init
147
148void
149hydro_finish(const hydroparam_t H, hydrovar_t * Hv) {
150 DFree(&Hv->uold, H.nvar * H.nxt * H.nyt);
151} // hydro_finish
152
153
154static void touchPage(real_t *adr, int lg) {
155 int i;
156#ifndef NOTOUCHPAGE
157#pragma omp parallel for private(i) shared(adr)
158 for(i = 0; i < lg; i++) {
159 adr[i] = 0.0l;
160 }
161#endif
162}
163
164
165void
166allocate_work_space(int ngrid, const hydroparam_t H, hydrowork_t * Hw, hydrovarwork_t * Hvw) {
167 int domain = ngrid * H.nxystep;
168 int domainVar = domain * H.nvar;
169 int domainD = domain * sizeof(real_t);
170 int domainI = domain * sizeof(int);
171 int domainVarD = domainVar * sizeof(real_t);
172 int pageM = 1024*1024;
173
174#define ONEBLOCK 1
175
176#ifndef PAGEOFFSET
177#define PAGEOFFSET sizeof(double)
178#endif
179
180#ifdef ONEBLOCK
181#ifndef TAILLEPAGE
182#define TAILLEPAGE 1024
183#endif
184 int oneBlock = 0;
185 int domainVarM = 0;
186 int domainM = 0;
187 int pageMD = TAILLEPAGE / 8 ;
188 real_t *blockD = 0;
189#endif
190
191 WHERE("allocate_work_space");
192
193#ifdef MOVETHEPAGES
194#ifndef __MIC__
195#define MOVEPAGEVAR(t) force_move_pages(t, domainVar, sizeof(real_t), HYDRO_NUMA_SIZED_BLOCK_RR, pageM)
196#define MOVEPAGE(t) force_move_pages(t, domain, sizeof(real_t), HYDRO_NUMA_SIZED_BLOCK_RR, pageM)
197#else
198#define MOVEPAGEVAR(t)
199#define MOVEPAGE(t)
200#endif
201#else
202#define MOVEPAGEVAR(t)
203#define MOVEPAGE(t)
204#endif
205
206#ifdef ONEBLOCK
207 if (H.mype == 0) fprintf(stdout, "Page offset %d\n", (int) PAGEOFFSET);
208 // determine the right amount of pages to fit all arrays
209 domainVarM = (domainVar + pageMD - 1) / pageMD;
210 domainVarM *= pageMD + PAGEOFFSET;
211 domainM = (domain + pageMD - 1) / pageMD;
212 domainM *= pageMD + PAGEOFFSET;
213
214 oneBlock = 9 * domainVarM + 12 * domainM; // expressed in term of pages of double
215 assert(oneBlock >= (9 * domainVar + 12 * domain));
216#pragma message "ONE BLOCK option"
217 blockD = (real_t *) malloc(oneBlock * sizeof(real_t));
218 assert(blockD != NULL);
219 if (((uint64_t) (&blockD[0]) & 63) == 0) {
220 fprintf(stderr, "ONE block allocated is not aligned \n");
221 }
222 Hvw->u = blockD; touchPage(Hvw->u, domainVar);
223 Hvw->q = Hvw->u + domainVarM; touchPage(Hvw->q, domainVar);
224 Hvw->dq = Hvw->q + domainVarM; touchPage(Hvw->dq, domainVar);
225 Hvw->qxm = Hvw->dq + domainVarM; touchPage(Hvw->qxm, domainVar);
226 Hvw->qxp = Hvw->qxm + domainVarM; touchPage(Hvw->qxp, domainVar);
227 Hvw->qleft = Hvw->qxp + domainVarM; touchPage(Hvw->qleft, domainVar);
228 Hvw->qright = Hvw->qleft + domainVarM; touchPage(Hvw->qright, domainVar);
229 Hvw->qgdnv = Hvw->qright + domainVarM; touchPage(Hvw->qgdnv, domainVar);
230 Hvw->flux = Hvw->qgdnv + domainVarM; touchPage(Hvw->flux, domainVar);
231 Hw->e = Hvw->flux + domainVarM; touchPage(Hw->e, domain);
232 Hw->c = Hw->e + domainM; touchPage(Hw->c, domain);
233 Hw->pstar = Hw->c + domainM; touchPage(Hw->pstar, domain);
234 Hw->rl = Hw->pstar + domainM; touchPage(Hw->rl, domain);
235 Hw->ul = Hw->rl + domainM; touchPage(Hw->ul, domain);
236 Hw->pl = Hw->ul + domainM; touchPage(Hw->pl, domain);
237 Hw->cl = Hw->pl + domainM; touchPage(Hw->cl, domain);
238 Hw->rr = Hw->cl + domainM; touchPage(Hw->rr, domain);
239 Hw->ur = Hw->rr + domainM; touchPage(Hw->ur, domain);
240 Hw->pr = Hw->ur + domainM; touchPage(Hw->pr, domain);
241 Hw->cr = Hw->pr + domainM; touchPage(Hw->cr, domain);
242 Hw->ro = Hw->cr + domainM; touchPage(Hw->ro, domain);
243#else
244 /*
245 force_move_pages(Hvw->u, domainVar, sizeof(double), HYDRO_NUMA_SIZED_BLOCK_RR, pageM);
246 */
247 fprintf(stderr, "Page malloc\n");
248 Hvw->u = DMalloc(domainVar); MOVEPAGEVAR(Hvw->u);
249 Hvw->q = DMalloc(domainVar); MOVEPAGEVAR(Hvw->q);
250 Hvw->dq = DMalloc(domainVar); MOVEPAGEVAR(Hvw->dq);
251 Hvw->qxm = DMalloc(domainVar); MOVEPAGEVAR(Hvw->qxm);
252 Hvw->qxp = DMalloc(domainVar); MOVEPAGEVAR(Hvw->qxp);
253 Hvw->qleft = DMalloc(domainVar); MOVEPAGEVAR(Hvw->qleft);
254 Hvw->qright = DMalloc(domainVar); MOVEPAGEVAR(Hvw->qright);
255 Hvw->qgdnv = DMalloc(domainVar); MOVEPAGEVAR(Hvw->qgdnv);
256 Hvw->flux = DMalloc(domainVar); MOVEPAGEVAR(Hvw->flux);
257 //
258 Hw->e = DMalloc(domain); MOVEPAGE(Hw->e);
259 Hw->c = DMalloc(domain); MOVEPAGE(Hw->c);
260 //
261 Hw->pstar = DMalloc(domain); MOVEPAGE(Hw->pstar);
262 Hw->rl = DMalloc(domain); MOVEPAGE(Hw->rl);
263 Hw->ul = DMalloc(domain); MOVEPAGE(Hw->ul);
264 Hw->pl = DMalloc(domain); MOVEPAGE(Hw->pl);
265 Hw->cl = DMalloc(domain); MOVEPAGE(Hw->cl);
266 Hw->rr = DMalloc(domain); MOVEPAGE(Hw->rr);
267 Hw->ur = DMalloc(domain); MOVEPAGE(Hw->ur);
268 Hw->pr = DMalloc(domain); MOVEPAGE(Hw->pr);
269 Hw->cr = DMalloc(domain); MOVEPAGE(Hw->cr);
270 Hw->ro = DMalloc(domain); MOVEPAGE(Hw->ro);
271#endif
272 Hw->goon = IMalloc(domain);
273 Hw->sgnm = IMalloc(domain);
274
275 // Hw->uo = DMalloc(ngrid);
276 // Hw->po = DMalloc(ngrid);
277 // Hw->co = DMalloc(ngrid);
278 // Hw->rstar = DMalloc(ngrid);
279 // Hw->ustar = DMalloc(ngrid);
280 // Hw->cstar = DMalloc(ngrid);
281 // Hw->wl = DMalloc(ngrid);
282 // Hw->wr = DMalloc(ngrid);
283 // Hw->wo = DMalloc((ngrid));
284 // Hw->spin = DMalloc(ngrid);
285 // Hw->spout = DMalloc(ngrid);
286 // Hw->ushock = DMalloc(ngrid);
287 // Hw->frac = DMalloc(ngrid);
288 // Hw->scr = DMalloc(ngrid);
289 // Hw->delp = DMalloc(ngrid);
290 // Hw->pold = DMalloc(ngrid);
291 // Hw->ind = IMalloc(ngrid);
292 // Hw->ind2 = IMalloc(ngrid);
293} // allocate_work_space
294
295void
296deallocate_work_space(int ngrid, const hydroparam_t H, hydrowork_t * Hw, hydrovarwork_t * Hvw) {
297 int domain = ngrid * H.nxystep;
298 int domainVar = domain * H.nvar;
299 int domainD = domain * sizeof(real_t);
300 int domainI = domain * sizeof(int);
301 int domainVarD = domainVar * sizeof(real_t);
302
303 WHERE("deallocate_work_space");
304#ifdef ONEBLOCK
305 int oneBlock = 0;
306 int domainVarM = 0;
307 int domainM = 0;
308 int pageM = 1024*1024;
309 int pageMD = TAILLEPAGE / 8;
310 real_t *blockD = 0;
311#endif
312
313 //
314#ifdef ONEBLOCK
315 // determine the right amount of pages to fit all arrays
316 domainVarM = (domainVar + pageMD - 1) / pageMD;
317 domainVarM *= pageMD + PAGEOFFSET;
318 domainM = (domain + pageMD - 1) / pageMD;
319 domainM *= pageMD + PAGEOFFSET;
320
321 oneBlock = 9 * domainVarM + 12 * domainM; // expressed in term of pages of double
322 DFree(&Hvw->u, oneBlock);
323 Hvw->q = Hvw->dq = Hvw->qxm = Hvw->qxp = 0;
324 Hvw->qleft = Hvw->qright = Hvw->qgdnv = Hvw->flux = Hw->e = Hw->c =0;
325 Hw->pstar = Hw->rl = Hw->ul = Hw->pl = Hw->cl = Hw->rr = Hw->ur = Hw->pr = Hw->cr = Hw->ro = 0;
326#else
327 DFree(&Hvw->u, domainVar);
328 DFree(&Hvw->q, domainVar);
329 DFree(&Hvw->dq, domainVar);
330 DFree(&Hvw->qxm, domainVar);
331 DFree(&Hvw->qxp, domainVar);
332 DFree(&Hvw->qleft, domainVar);
333 DFree(&Hvw->qright, domainVar);
334 DFree(&Hvw->qgdnv, domainVar);
335 DFree(&Hvw->flux, domainVar);
336 DFree(&Hw->e, domain);
337 DFree(&Hw->c, domain);
338 DFree(&Hw->pstar, domain);
339 DFree(&Hw->rl, domain);
340 DFree(&Hw->ul, domain);
341 DFree(&Hw->pl, domain);
342 DFree(&Hw->cl, domain);
343 DFree(&Hw->rr, domain);
344 DFree(&Hw->ur, domain);
345 DFree(&Hw->pr, domain);
346 DFree(&Hw->cr, domain);
347 DFree(&Hw->ro, domain);
348#endif
349
350 IFree(&Hw->sgnm, domainVar);
351 IFree(&Hw->goon, domain);
352 // Free(Hw->uo);
353 // Free(Hw->po);
354 // Free(Hw->co);
355 // Free(Hw->rstar);
356 // Free(Hw->ustar);
357 // Free(Hw->cstar);
358 // Free(Hw->wl);
359 // Free(Hw->wr);
360 // Free(Hw->wo);
361 // Free(Hw->spin);
362 // Free(Hw->spout);
363 // Free(Hw->ushock);
364 // Free(Hw->frac);
365 // Free(Hw->scr);
366 // Free(Hw->delp);
367 // Free(Hw->pold);
368 // Free(Hw->ind);
369 // Free(Hw->ind2);
370} // deallocate_work_space
371
372
373// EOF
Note: See TracBrowser for help on using the repository browser.