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

151: TEST*/