cell.c 15.7 KB
Newer Older
1
/*******************************************************************************
2
 * This file is part of SWIFT.
3
 * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
4
 *
5
6
7
8
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
9
 *
10
11
12
13
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
14
 *
15
16
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
17
 *
18
19
20
21
22
23
24
25
26
 ******************************************************************************/

/* Config parameters. */
#include "../config.h"

/* Some standard headers. */
#include <float.h>
#include <limits.h>
#include <math.h>
27
28
29
30
#include <pthread.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
31

32
33
/* MPI headers. */
#ifdef WITH_MPI
34
#include <mpi.h>
35
36
#endif

37
38
/* Switch off timers. */
#ifdef TIMER
39
#undef TIMER
40
41
#endif

42
43
44
/* This object's header. */
#include "cell.h"

45
/* Local headers. */
46
#include "atomic.h"
47
#include "error.h"
48
#include "hydro.h"
49
50
#include "space.h"
#include "timers.h"
51

52
53
54
/* Global variables. */
int cell_next_tag = 0;

55
56
57
58
59
60
/**
 * @brief Get the size of the cell subtree.
 *
 * @param c The #cell.
 */

61
int cell_getsize(struct cell *c) {
62

Pedro Gonnet's avatar
Pedro Gonnet committed
63
64
  /* Number of cells in this subtree. */
  int count = 1;
65

66
67
  /* Sum up the progeny if split. */
  if (c->split)
Pedro Gonnet's avatar
Pedro Gonnet committed
68
    for (int k = 0; k < 8; k++)
69
70
71
72
73
74
75
      if (c->progeny[k] != NULL) count += cell_getsize(c->progeny[k]);

  /* Return the final count. */
  return count;
}

/**
76
77
78
79
80
81
82
83
84
 * @brief Unpack the data of a given cell and its sub-cells.
 *
 * @param pc An array of packed #pcell.
 * @param c The #cell in which to unpack the #pcell.
 * @param s The #space in which the cells are created.
 *
 * @return The number of cells created.
 */

85
86
87
88
int cell_unpack(struct pcell *pc, struct cell *c, struct space *s) {

  /* Unpack the current pcell. */
  c->h_max = pc->h_max;
89
90
  c->ti_end_min = pc->ti_end_min;
  c->ti_end_max = pc->ti_end_max;
91
92
93
94
  c->count = pc->count;
  c->tag = pc->tag;

  /* Fill the progeny recursively, depth-first. */
Pedro Gonnet's avatar
Pedro Gonnet committed
95
96
  int count = 1;
  for (int k = 0; k < 8; k++)
97
    if (pc->progeny[k] >= 0) {
Pedro Gonnet's avatar
Pedro Gonnet committed
98
      struct cell *temp = space_getcell(s);
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
      temp->count = 0;
      temp->loc[0] = c->loc[0];
      temp->loc[1] = c->loc[1];
      temp->loc[2] = c->loc[2];
      temp->h[0] = c->h[0] / 2;
      temp->h[1] = c->h[1] / 2;
      temp->h[2] = c->h[2] / 2;
      temp->dmin = c->dmin / 2;
      if (k & 4) temp->loc[0] += temp->h[0];
      if (k & 2) temp->loc[1] += temp->h[1];
      if (k & 1) temp->loc[2] += temp->h[2];
      temp->depth = c->depth + 1;
      temp->split = 0;
      temp->dx_max = 0.0;
      temp->nodeID = c->nodeID;
      temp->parent = c;
      c->progeny[k] = temp;
      c->split = 1;
      count += cell_unpack(&pc[pc->progeny[k]], temp, s);
118
119
    }

120
121
122
  /* Return the total number of unpacked cells. */
  return count;
}
123

124
125
126
127
128
129
130
131
132
/**
 * @brief Link the cells recursively to the given part array.
 *
 * @param c The #cell.
 * @param parts The #part array.
 *
 * @return The number of particles linked.
 */

133
int cell_link(struct cell *c, struct part *parts) {
134

135
136
137
  c->parts = parts;

  /* Fill the progeny recursively, depth-first. */
Pedro Gonnet's avatar
Pedro Gonnet committed
138
139
140
141
142
143
144
  if (c->split) {
    int offset = 0;
    for (int k = 0; k < 8; k++) {
      if (c->progeny[k] != NULL)
        offset += cell_link(c->progeny[k], &parts[offset]);
    }
  }
145
146
147
148

  /* Return the total number of unpacked cells. */
  return c->count;
}
149

150
151
152
153
154
155
156
157
158
159
/**
 * @brief Pack the data of the given cell and all it's sub-cells.
 *
 * @param c The #cell.
 * @param pc Pointer to an array of packed cells in which the
 *      cells will be packed.
 *
 * @return The number of packed cells.
 */

160
161
162
163
int cell_pack(struct cell *c, struct pcell *pc) {

  /* Start by packing the data of the current cell. */
  pc->h_max = c->h_max;
164
165
  pc->ti_end_min = c->ti_end_min;
  pc->ti_end_max = c->ti_end_max;
166
167
168
169
  pc->count = c->count;
  c->tag = pc->tag = atomic_inc(&cell_next_tag) % cell_max_tag;

  /* Fill in the progeny, depth-first recursion. */
Pedro Gonnet's avatar
Pedro Gonnet committed
170
171
  int count = 1;
  for (int k = 0; k < 8; k++)
172
173
174
175
176
177
178
179
180
    if (c->progeny[k] != NULL) {
      pc->progeny[k] = count;
      count += cell_pack(c->progeny[k], &pc[count]);
    } else
      pc->progeny[k] = -1;

  /* Return the number of packed cells used. */
  return count;
}
181

182
183
184
185
186
187
/**
 * @brief Lock a cell and hold its parents.
 *
 * @param c The #cell.
 */

188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
int cell_locktree(struct cell *c) {

  TIMER_TIC

  /* First of all, try to lock this cell. */
  if (c->hold || lock_trylock(&c->lock) != 0) {
    TIMER_TOC(timer_locktree);
    return 1;
  }

  /* Did somebody hold this cell in the meantime? */
  if (c->hold) {

    /* Unlock this cell. */
    if (lock_unlock(&c->lock) != 0) error("Failed to unlock cell.");

    /* Admit defeat. */
    TIMER_TOC(timer_locktree);
    return 1;
  }

  /* Climb up the tree and lock/hold/unlock. */
Pedro Gonnet's avatar
Pedro Gonnet committed
210
  struct cell *finger;
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
  for (finger = c->parent; finger != NULL; finger = finger->parent) {

    /* Lock this cell. */
    if (lock_trylock(&finger->lock) != 0) break;

    /* Increment the hold. */
    atomic_inc(&finger->hold);

    /* Unlock the cell. */
    if (lock_unlock(&finger->lock) != 0) error("Failed to unlock cell.");
  }

  /* If we reached the top of the tree, we're done. */
  if (finger == NULL) {
    TIMER_TOC(timer_locktree);
    return 0;
  }

  /* Otherwise, we hit a snag. */
  else {

    /* Undo the holds up to finger. */
Pedro Gonnet's avatar
Pedro Gonnet committed
233
234
    for (struct cell *finger2 = c->parent; finger2 != finger;
         finger2 = finger2->parent)
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
      __sync_fetch_and_sub(&finger2->hold, 1);

    /* Unlock this cell. */
    if (lock_unlock(&c->lock) != 0) error("Failed to unlock cell.");

    /* Admit defeat. */
    TIMER_TOC(timer_locktree);
    return 1;
  }
}

int cell_glocktree(struct cell *c) {

  TIMER_TIC

  /* First of all, try to lock this cell. */
  if (c->ghold || lock_trylock(&c->glock) != 0) {
    TIMER_TOC(timer_locktree);
    return 1;
  }

  /* Did somebody hold this cell in the meantime? */
  if (c->ghold) {

    /* Unlock this cell. */
    if (lock_unlock(&c->glock) != 0) error("Failed to unlock cell.");

    /* Admit defeat. */
    TIMER_TOC(timer_locktree);
    return 1;
  }

  /* Climb up the tree and lock/hold/unlock. */
Pedro Gonnet's avatar
Pedro Gonnet committed
268
  struct cell *finger;
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
  for (finger = c->parent; finger != NULL; finger = finger->parent) {

    /* Lock this cell. */
    if (lock_trylock(&finger->glock) != 0) break;

    /* Increment the hold. */
    __sync_fetch_and_add(&finger->ghold, 1);

    /* Unlock the cell. */
    if (lock_unlock(&finger->glock) != 0) error("Failed to unlock cell.");
  }

  /* If we reached the top of the tree, we're done. */
  if (finger == NULL) {
    TIMER_TOC(timer_locktree);
    return 0;
  }

  /* Otherwise, we hit a snag. */
  else {

    /* Undo the holds up to finger. */
Pedro Gonnet's avatar
Pedro Gonnet committed
291
292
    for (struct cell *finger2 = c->parent; finger2 != finger;
         finger2 = finger2->parent)
293
294
295
296
297
298
299
300
301
302
      __sync_fetch_and_sub(&finger2->ghold, 1);

    /* Unlock this cell. */
    if (lock_unlock(&c->glock) != 0) error("Failed to unlock cell.");

    /* Admit defeat. */
    TIMER_TOC(timer_locktree);
    return 1;
  }
}
303

304
/**
305
 * @brief Unlock a cell's parents.
306
307
308
 *
 * @param c The #cell.
 */
309
310
311
312
313
314
315
316
317

void cell_unlocktree(struct cell *c) {

  TIMER_TIC

  /* First of all, try to unlock this cell. */
  if (lock_unlock(&c->lock) != 0) error("Failed to unlock cell.");

  /* Climb up the tree and unhold the parents. */
Pedro Gonnet's avatar
Pedro Gonnet committed
318
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
319
320
321
322
323
324
325
326
327
328
329
330
331
    __sync_fetch_and_sub(&finger->hold, 1);

  TIMER_TOC(timer_locktree);
}

void cell_gunlocktree(struct cell *c) {

  TIMER_TIC

  /* First of all, try to unlock this cell. */
  if (lock_unlock(&c->glock) != 0) error("Failed to unlock cell.");

  /* Climb up the tree and unhold the parents. */
Pedro Gonnet's avatar
Pedro Gonnet committed
332
  for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
333
334
335
336
337
    __sync_fetch_and_sub(&finger->ghold, 1);

  TIMER_TOC(timer_locktree);
}

338
339
340
341
342
/**
 * @brief Sort the parts into eight bins along the given pivots.
 *
 * @param c The #cell array to be sorted.
 */
343
344
345

void cell_split(struct cell *c) {

Pedro Gonnet's avatar
Pedro Gonnet committed
346
347
348
349
350
  int i, j;
  const int count = c->count, gcount = c->gcount;
  struct part* parts = c->parts;
  struct xpart* xparts = c->xparts;
  struct gpart* gparts = c->gparts;
351
352
353
354
  int left[8], right[8];
  double pivot[3];

  /* Init the pivots. */
Pedro Gonnet's avatar
Pedro Gonnet committed
355
  for (int k = 0; k < 3; k++) pivot[k] = c->loc[k] + c->h[k] / 2;
356
357
358
359
360
361
362
363

  /* Split along the x-axis. */
  i = 0;
  j = count - 1;
  while (i <= j) {
    while (i <= count - 1 && parts[i].x[0] <= pivot[0]) i += 1;
    while (j >= 0 && parts[j].x[0] > pivot[0]) j -= 1;
    if (i < j) {
Pedro Gonnet's avatar
Pedro Gonnet committed
364
      struct part temp = parts[i];
365
366
      parts[i] = parts[j];
      parts[j] = temp;
Pedro Gonnet's avatar
Pedro Gonnet committed
367
      struct xpart xtemp = xparts[i];
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
      xparts[i] = xparts[j];
      xparts[j] = xtemp;
    }
  }
  /* for ( k = 0 ; k <= j ; k++ )
      if ( parts[k].x[0] > pivot[0] )
          error( "cell_split: sorting failed." );
  for ( k = i ; k < count ; k++ )
      if ( parts[k].x[0] < pivot[0] )
          error( "cell_split: sorting failed." ); */
  left[1] = i;
  right[1] = count - 1;
  left[0] = 0;
  right[0] = j;

  /* Split along the y axis, twice. */
Pedro Gonnet's avatar
Pedro Gonnet committed
384
  for (int k = 1; k >= 0; k--) {
385
386
387
388
389
390
    i = left[k];
    j = right[k];
    while (i <= j) {
      while (i <= right[k] && parts[i].x[1] <= pivot[1]) i += 1;
      while (j >= left[k] && parts[j].x[1] > pivot[1]) j -= 1;
      if (i < j) {
Pedro Gonnet's avatar
Pedro Gonnet committed
391
        struct part temp = parts[i];
392
393
        parts[i] = parts[j];
        parts[j] = temp;
Pedro Gonnet's avatar
Pedro Gonnet committed
394
        struct xpart xtemp = xparts[i];
395
396
397
398
399
400
401
402
        xparts[i] = xparts[j];
        xparts[j] = xtemp;
      }
    }
    /* for ( int kk = left[k] ; kk <= j ; kk++ )
        if ( parts[kk].x[1] > pivot[1] ) {
            message( "ival=[%i,%i], i=%i, j=%i." , left[k] , right[k] , i , j );
            error( "sorting failed (left)." );
403
            }
404
405
406
407
408
409
410
411
412
413
    for ( int kk = i ; kk <= right[k] ; kk++ )
        if ( parts[kk].x[1] < pivot[1] )
            error( "sorting failed (right)." ); */
    left[2 * k + 1] = i;
    right[2 * k + 1] = right[k];
    left[2 * k] = left[k];
    right[2 * k] = j;
  }

  /* Split along the z axis, four times. */
Pedro Gonnet's avatar
Pedro Gonnet committed
414
  for (int k = 3; k >= 0; k--) {
415
416
417
418
419
420
    i = left[k];
    j = right[k];
    while (i <= j) {
      while (i <= right[k] && parts[i].x[2] <= pivot[2]) i += 1;
      while (j >= left[k] && parts[j].x[2] > pivot[2]) j -= 1;
      if (i < j) {
Pedro Gonnet's avatar
Pedro Gonnet committed
421
        struct part temp = parts[i];
422
423
        parts[i] = parts[j];
        parts[j] = temp;
Pedro Gonnet's avatar
Pedro Gonnet committed
424
        struct xpart xtemp = xparts[i];
425
426
427
428
429
430
431
432
        xparts[i] = xparts[j];
        xparts[j] = xtemp;
      }
    }
    /* for ( int kk = left[k] ; kk <= j ; kk++ )
        if ( parts[kk].x[2] > pivot[2] ) {
            message( "ival=[%i,%i], i=%i, j=%i." , left[k] , right[k] , i , j );
            error( "sorting failed (left)." );
433
            }
434
435
436
437
438
439
440
441
442
443
444
445
    for ( int kk = i ; kk <= right[k] ; kk++ )
        if ( parts[kk].x[2] < pivot[2] ) {
            message( "ival=[%i,%i], i=%i, j=%i." , left[k] , right[k] , i , j );
            error( "sorting failed (right)." );
            } */
    left[2 * k + 1] = i;
    right[2 * k + 1] = right[k];
    left[2 * k] = left[k];
    right[2 * k] = j;
  }

  /* Store the counts and offsets. */
Pedro Gonnet's avatar
Pedro Gonnet committed
446
  for (int k = 0; k < 8; k++) {
447
448
449
450
451
452
    c->progeny[k]->count = right[k] - left[k] + 1;
    c->progeny[k]->parts = &c->parts[left[k]];
    c->progeny[k]->xparts = &c->xparts[left[k]];
  }

  /* Re-link the gparts. */
Pedro Gonnet's avatar
Pedro Gonnet committed
453
  for (int k = 0; k < count; k++)
454
455
456
457
458
459
460
461
462
463
464
465
466
    if (parts[k].gpart != NULL) parts[k].gpart->part = &parts[k];

  /* Verify that _all_ the parts have been assigned to a cell. */
  /* for ( k = 1 ; k < 8 ; k++ )
      if ( &c->progeny[k-1]->parts[ c->progeny[k-1]->count ] !=
  c->progeny[k]->parts )
          error( "Particle sorting failed (internal consistency)." );
  if ( c->progeny[0]->parts != c->parts )
      error( "Particle sorting failed (left edge)." );
  if ( &c->progeny[7]->parts[ c->progeny[7]->count ] != &c->parts[ count ] )
      error( "Particle sorting failed (right edge)." ); */

  /* Verify a few sub-cells. */
Pedro Gonnet's avatar
Pedro Gonnet committed
467
  /* for (int k = 0 ; k < c->progeny[0]->count ; k++ )
468
469
470
471
      if ( c->progeny[0]->parts[k].x[0] > pivot[0] ||
           c->progeny[0]->parts[k].x[1] > pivot[1] ||
           c->progeny[0]->parts[k].x[2] > pivot[2] )
          error( "Sorting failed (progeny=0)." );
Pedro Gonnet's avatar
Pedro Gonnet committed
472
  for (int k = 0 ; k < c->progeny[1]->count ; k++ )
473
474
475
476
      if ( c->progeny[1]->parts[k].x[0] > pivot[0] ||
           c->progeny[1]->parts[k].x[1] > pivot[1] ||
           c->progeny[1]->parts[k].x[2] <= pivot[2] )
          error( "Sorting failed (progeny=1)." );
Pedro Gonnet's avatar
Pedro Gonnet committed
477
  for (int k = 0 ; k < c->progeny[2]->count ; k++ )
478
479
480
481
482
483
484
485
486
487
488
489
490
491
      if ( c->progeny[2]->parts[k].x[0] > pivot[0] ||
           c->progeny[2]->parts[k].x[1] <= pivot[1] ||
           c->progeny[2]->parts[k].x[2] > pivot[2] )
          error( "Sorting failed (progeny=2)." ); */

  /* Now do the same song and dance for the gparts. */

  /* Split along the x-axis. */
  i = 0;
  j = gcount - 1;
  while (i <= j) {
    while (i <= gcount - 1 && gparts[i].x[0] <= pivot[0]) i += 1;
    while (j >= 0 && gparts[j].x[0] > pivot[0]) j -= 1;
    if (i < j) {
Pedro Gonnet's avatar
Pedro Gonnet committed
492
      struct gpart gtemp = gparts[i];
493
494
      gparts[i] = gparts[j];
      gparts[j] = gtemp;
495
    }
496
497
498
499
500
501
502
  }
  left[1] = i;
  right[1] = gcount - 1;
  left[0] = 0;
  right[0] = j;

  /* Split along the y axis, twice. */
Pedro Gonnet's avatar
Pedro Gonnet committed
503
  for (int k = 1; k >= 0; k--) {
504
505
506
507
508
509
    i = left[k];
    j = right[k];
    while (i <= j) {
      while (i <= right[k] && gparts[i].x[1] <= pivot[1]) i += 1;
      while (j >= left[k] && gparts[j].x[1] > pivot[1]) j -= 1;
      if (i < j) {
Pedro Gonnet's avatar
Pedro Gonnet committed
510
        struct gpart gtemp = gparts[i];
511
512
513
514
515
516
517
518
519
520
521
        gparts[i] = gparts[j];
        gparts[j] = gtemp;
      }
    }
    left[2 * k + 1] = i;
    right[2 * k + 1] = right[k];
    left[2 * k] = left[k];
    right[2 * k] = j;
  }

  /* Split along the z axis, four times. */
Pedro Gonnet's avatar
Pedro Gonnet committed
522
  for (int k = 3; k >= 0; k--) {
523
524
525
526
527
528
    i = left[k];
    j = right[k];
    while (i <= j) {
      while (i <= right[k] && gparts[i].x[2] <= pivot[2]) i += 1;
      while (j >= left[k] && gparts[j].x[2] > pivot[2]) j -= 1;
      if (i < j) {
Pedro Gonnet's avatar
Pedro Gonnet committed
529
        struct gpart gtemp = gparts[i];
530
531
532
533
534
535
536
537
538
539
540
        gparts[i] = gparts[j];
        gparts[j] = gtemp;
      }
    }
    left[2 * k + 1] = i;
    right[2 * k + 1] = right[k];
    left[2 * k] = left[k];
    right[2 * k] = j;
  }

  /* Store the counts and offsets. */
Pedro Gonnet's avatar
Pedro Gonnet committed
541
  for (int k = 0; k < 8; k++) {
542
543
544
545
546
    c->progeny[k]->gcount = right[k] - left[k] + 1;
    c->progeny[k]->gparts = &c->gparts[left[k]];
  }

  /* Re-link the parts. */
Pedro Gonnet's avatar
Pedro Gonnet committed
547
  for (int k = 0; k < gcount; k++)
548
549
    if (gparts[k].id > 0) gparts[k].part->gpart = &gparts[k];
}
550
551
552
553
554
555
556
557
558
559
560

/**
 * @brief Initialises all particles to a valid state even if the ICs were stupid
 *
 * @param c Cell to act upon
 * @param data Unused parameter
 */
void cell_init_parts(struct cell *c, void *data) {

  struct part *p = c->parts;
  struct xpart *xp = c->xparts;
561
  const int count = c->count;
562

563
  for (int i = 0; i < count; ++i) {
564
565
    p[i].ti_begin = 0;
    p[i].ti_end = 0;
566
567
568
    xp[i].v_full[0] = p[i].v[0];
    xp[i].v_full[1] = p[i].v[1];
    xp[i].v_full[2] = p[i].v[2];
569
    hydro_first_init_part(&p[i], &xp[i]);
570
571
    hydro_init_part(&p[i]);
    hydro_reset_acceleration(&p[i]);
572
  }
573
574
  c->ti_end_min = 0;
  c->ti_end_max = 0;
575
576
}

577
/**
578
579
 * @brief Converts hydro quantities to a valid state after the initial density
 *calculation
580
581
582
583
584
585
586
 *
 * @param c Cell to act upon
 * @param data Unused parameter
 */
void cell_convert_hydro(struct cell *c, void *data) {

  struct part *p = c->parts;
587
588

  for (int i = 0; i < c->count; ++i) {
589
590
591
592
    hydro_convert_quantities(&p[i]);
  }
}

Matthieu Schaller's avatar
Matthieu Schaller committed
593
594
595
596
597
598
/**
 * @brief Cleans the links in a given cell.
 *
 * @param c Cell to act upon
 * @param data Unused parameter
 */
599
void cell_clean_links(struct cell *c, void *data) {
Matthieu Schaller's avatar
Matthieu Schaller committed
600
601
  c->density = NULL;
  c->nr_density = 0;
602

Matthieu Schaller's avatar
Matthieu Schaller committed
603
604
605
  c->force = NULL;
  c->nr_force = 0;
}