1 static char help[] = "Demonstrates Conway's Game of Life using a 2d DMDA.\n\n";
2
3 /*
4 At each step in time, the following transitions occur:
5
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.
10
11 https://en.wikipedia.org/wiki/Conway%27s_Game_of_Life
12 */
13
14 #include <petscdm.h>
15 #include <petscdmda.h>
16
17 static const int GLIDER[3][3] = {
18 {0, 1, 0},
19 {0, 1, 1},
20 {1, 0, 1}
21 };
22
main(int argc,char ** argv)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;
31
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();
44
45 PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, NULL, "Life", PETSC_DECIDE, PETSC_DECIDE, 1000, 1000, &viewer));
46
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));
53
54 { /* Initialize */
55 DMDALocalInfo info;
56 PetscScalar **x;
57 PetscInt i, j;
58
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 }
78
79 /* View the initial condition */
80 PetscCall(VecView(Xglobal, viewer));
81
82 { /* Play */
83 PetscInt step;
84
85 for (step = 0; step < steps; step++) {
86 const PetscScalar **x;
87 PetscScalar **y;
88 DMDALocalInfo info;
89 PetscInt i, j;
90
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 }
135
136 PetscCall(PetscViewerDestroy(&viewer));
137 PetscCall(VecDestroy(&Xglobal));
138 PetscCall(VecDestroy(&Xlocal));
139 PetscCall(DMDestroy(&da));
140 PetscCall(PetscFinalize());
141 return 0;
142 }
143
144 /*TEST
145
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
151
152 TEST*/
153