Actual source code: ex2.c

      1: static char help[] = "Demonstrates Conway's Game of Life using a 2d DMDA.\n\n";

      3: /*
      4:  At each step in time, the following transitions occur:

      6:     Any live cell with fewer than two live neighbours dies, as if by underpopulation.
      7:     Any live cell with two or three live neighbours lives on to the next generation.
      8:     Any live cell with more than three live neighbours dies, as if by overpopulation.
      9:     Any dead cell with exactly three live neighbours becomes a live cell, as if by reproduction.

     11:  https://en.wikipedia.org/wiki/Conway%27s_Game_of_Life
     12: */

     14: #include <petscdm.h>
     15: #include <petscdmda.h>

     17: static const int GLIDER[3][3] = {
     18:   {0, 1, 0},
     19:   {0, 1, 1},
     20:   {1, 0, 1}
     21: };

     23: int main(int argc, char **argv)
     24: {
     25:   DM          da;
     26:   PetscViewer viewer;
     27:   Vec         Xlocal, Xglobal;
     28:   PetscInt    glider_loc[2] = {10, 20}, blinker_loc[2] = {20, 10}, two, steps = 100, viz_interval = 1;
     29:   PetscInt    check_step_alive = -1, check_step_dead = -1;
     30:   PetscBool   has_glider, has_blinker;

     32:   PetscFunctionBeginUser;
     33:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
     34:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Conway's Game of Life", "");
     35:   {
     36:     PetscCall(PetscOptionsIntArray("-glider", "Coordinate at which to center a glider", NULL, glider_loc, (two = 2, &two), &has_glider));
     37:     PetscCall(PetscOptionsIntArray("-blinker", "Coordinate at which to center a blinker", NULL, blinker_loc, (two = 2, &two), &has_blinker));
     38:     PetscCall(PetscOptionsInt("-steps", "Number of steps to take", NULL, steps, &steps, NULL));
     39:     PetscCall(PetscOptionsInt("-viz_interval", "Visualization interval", NULL, viz_interval, &viz_interval, NULL));
     40:     PetscCall(PetscOptionsInt("-check_step_alive", "Step on which to check that the simulation is alive", NULL, check_step_alive, &check_step_alive, NULL));
     41:     PetscCall(PetscOptionsInt("-check_step_dead", "Step on which to check that the simulation is dead", NULL, check_step_dead, &check_step_dead, NULL));
     42:   }
     43:   PetscOptionsEnd();

     45:   PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, NULL, "Life", PETSC_DECIDE, PETSC_DECIDE, 1000, 1000, &viewer));

     47:   /* Create distributed array and get vectors */
     48:   PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, DMDA_STENCIL_BOX, 30, 30, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da));
     49:   PetscCall(DMSetFromOptions(da));
     50:   PetscCall(DMSetUp(da));
     51:   PetscCall(DMCreateLocalVector(da, &Xlocal));
     52:   PetscCall(DMCreateGlobalVector(da, &Xglobal));

     54:   { /* Initialize */
     55:     DMDALocalInfo info;
     56:     PetscScalar **x;
     57:     PetscInt      i, j;

     59:     PetscCall(DMDAGetLocalInfo(da, &info));
     60:     PetscCall(DMDAVecGetArray(da, Xlocal, &x));
     61:     for (j = info.ys; j < info.ys + info.ym; j++) {
     62:       for (i = info.xs; i < info.xs + info.xm; i++) {
     63:         if (has_glider && i == glider_loc[0] && j == glider_loc[1]) {
     64:           PetscInt ii, jj;
     65:           for (ii = -1; ii <= 1; ii++)
     66:             for (jj = -1; jj <= 1; jj++) x[j + jj][i + ii] = GLIDER[1 - jj][ii + 1];
     67:         }
     68:         if (has_blinker && i == blinker_loc[0] && j == blinker_loc[1]) {
     69:           x[j - 1][i] = 1;
     70:           x[j][i]     = 1;
     71:           x[j + 1][i] = 1;
     72:         }
     73:       }
     74:     }
     75:     PetscCall(DMDAVecRestoreArray(da, Xlocal, &x));
     76:     PetscCall(DMLocalToGlobal(da, Xlocal, ADD_VALUES, Xglobal));
     77:   }

     79:   /* View the initial condition */
     80:   PetscCall(VecView(Xglobal, viewer));

     82:   { /* Play */
     83:     PetscInt step;

     85:     for (step = 0; step < steps; step++) {
     86:       const PetscScalar **x;
     87:       PetscScalar       **y;
     88:       DMDALocalInfo       info;
     89:       PetscInt            i, j;

     91:       PetscCall(DMGlobalToLocal(da, Xglobal, INSERT_VALUES, Xlocal));
     92:       PetscCall(DMDAGetLocalInfo(da, &info));
     93:       PetscCall(DMDAVecGetArrayRead(da, Xlocal, (void *)&x));
     94:       PetscCall(DMDAVecGetArrayWrite(da, Xglobal, &y));
     95:       for (j = info.ys; j < info.ys + info.ym; j++) {
     96:         for (i = info.xs; i < info.xs + info.xm; i++) {
     97:           PetscInt live_neighbors = 0;
     98:           live_neighbors += PetscRealPart(x[j - 1][i - 1]) > 0;
     99:           live_neighbors += PetscRealPart(x[j - 1][i]) > 0;
    100:           live_neighbors += PetscRealPart(x[j - 1][i + 1]) > 0;
    101:           live_neighbors += PetscRealPart(x[j][i - 1]) > 0;
    102:           live_neighbors += PetscRealPart(x[j][i + 1]) > 0;
    103:           live_neighbors += PetscRealPart(x[j + 1][i - 1]) > 0;
    104:           live_neighbors += PetscRealPart(x[j + 1][i]) > 0;
    105:           live_neighbors += PetscRealPart(x[j + 1][i + 1]) > 0;
    106:           if (PetscRealPart(x[j][i]) > 0) { /* Live cell */
    107:             switch (live_neighbors) {
    108:             case 2:
    109:             case 3:
    110:               y[j][i] = 1; /* Survive */
    111:               break;
    112:             default:
    113:               y[j][i] = 0; /* Death */
    114:             }
    115:           } else {                                /* Dead cell */
    116:             if (live_neighbors == 3) y[j][i] = 1; /* Birth */
    117:             else y[j][i] = 0;
    118:           }
    119:         }
    120:       }
    121:       PetscCall(DMDAVecRestoreArrayRead(da, Xlocal, (void *)&x));
    122:       PetscCall(DMDAVecRestoreArrayWrite(da, Xglobal, &y));
    123:       if (step == check_step_alive || step == check_step_dead) {
    124:         PetscScalar sum;
    125:         PetscCall(VecSum(Xglobal, &sum));
    126:         if (PetscAbsScalar(sum) > 0.1) {
    127:           if (step == check_step_dead) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Simulation alive at step %" PetscInt_FMT "\n", step));
    128:         } else if (step == check_step_alive) {
    129:           PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Simulation dead at step %" PetscInt_FMT "\n", step));
    130:         }
    131:       }
    132:       if (step % viz_interval == 0) PetscCall(VecView(Xglobal, viewer));
    133:     }
    134:   }

    136:   PetscCall(PetscViewerDestroy(&viewer));
    137:   PetscCall(VecDestroy(&Xglobal));
    138:   PetscCall(VecDestroy(&Xlocal));
    139:   PetscCall(DMDestroy(&da));
    140:   PetscCall(PetscFinalize());
    141:   return 0;
    142: }

    144: /*TEST

    146:    test:
    147:       requires: x
    148:       nsize: 2
    149:       args: -glider 5,6 -blinker 12,12 -steps 35 -check_step_alive 31 -check_step_dead 32 -da_grid_x 20 -da_grid_y 20 -nox
    150:       output_file: output/empty.out

    152: TEST*/