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 23 int main(int argc, char **argv) { 24 DM da; 25 PetscViewer viewer; 26 Vec Xlocal, Xglobal; 27 PetscInt glider_loc[2] = {10, 20}, blinker_loc[2] = {20, 10}, two, steps = 100, viz_interval = 1; 28 PetscInt check_step_alive = -1, check_step_dead = -1; 29 PetscBool has_glider, has_blinker; 30 31 PetscFunctionBeginUser; 32 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 33 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Conway's Game of Life", ""); 34 { 35 PetscCall(PetscOptionsIntArray("-glider", "Coordinate at which to center a glider", NULL, glider_loc, (two = 2, &two), &has_glider)); 36 PetscCall(PetscOptionsIntArray("-blinker", "Coordinate at which to center a blinker", NULL, blinker_loc, (two = 2, &two), &has_blinker)); 37 PetscCall(PetscOptionsInt("-steps", "Number of steps to take", NULL, steps, &steps, NULL)); 38 PetscCall(PetscOptionsInt("-viz_interval", "Vizualization interval", NULL, viz_interval, &viz_interval, NULL)); 39 PetscCall(PetscOptionsInt("-check_step_alive", "Step on which to check that the simulation is alive", NULL, check_step_alive, &check_step_alive, NULL)); 40 PetscCall(PetscOptionsInt("-check_step_dead", "Step on which to check that the simulation is dead", NULL, check_step_dead, &check_step_dead, NULL)); 41 } 42 PetscOptionsEnd(); 43 44 PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, NULL, "Life", PETSC_DECIDE, PETSC_DECIDE, 1000, 1000, &viewer)); 45 46 /* Create distributed array and get vectors */ 47 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)); 48 PetscCall(DMSetFromOptions(da)); 49 PetscCall(DMSetUp(da)); 50 PetscCall(DMCreateLocalVector(da, &Xlocal)); 51 PetscCall(DMCreateGlobalVector(da, &Xglobal)); 52 53 { /* Initialize */ 54 DMDALocalInfo info; 55 PetscScalar **x; 56 PetscInt i, j; 57 58 PetscCall(DMDAGetLocalInfo(da, &info)); 59 PetscCall(DMDAVecGetArray(da, Xlocal, &x)); 60 for (j = info.ys; j < info.ys + info.ym; j++) { 61 for (i = info.xs; i < info.xs + info.xm; i++) { 62 if (has_glider && i == glider_loc[0] && j == glider_loc[1]) { 63 PetscInt ii, jj; 64 for (ii = -1; ii <= 1; ii++) 65 for (jj = -1; jj <= 1; jj++) x[j + jj][i + ii] = GLIDER[1 - jj][ii + 1]; 66 } 67 if (has_blinker && i == blinker_loc[0] && j == blinker_loc[1]) { 68 x[j - 1][i] = 1; 69 x[j][i] = 1; 70 x[j + 1][i] = 1; 71 } 72 } 73 } 74 PetscCall(DMDAVecRestoreArray(da, Xlocal, &x)); 75 PetscCall(DMLocalToGlobal(da, Xlocal, ADD_VALUES, Xglobal)); 76 } 77 78 /* View the initial condition */ 79 PetscCall(VecView(Xglobal, viewer)); 80 81 { /* Play */ 82 PetscInt step; 83 84 for (step = 0; step < steps; step++) { 85 const PetscScalar **x; 86 PetscScalar **y; 87 DMDALocalInfo info; 88 PetscInt i, j; 89 90 PetscCall(DMGlobalToLocal(da, Xglobal, INSERT_VALUES, Xlocal)); 91 PetscCall(DMDAGetLocalInfo(da, &info)); 92 PetscCall(DMDAVecGetArrayRead(da, Xlocal, (void *)&x)); 93 PetscCall(DMDAVecGetArrayWrite(da, Xglobal, &y)); 94 for (j = info.ys; j < info.ys + info.ym; j++) { 95 for (i = info.xs; i < info.xs + info.xm; i++) { 96 PetscInt live_neighbors = 0; 97 live_neighbors += PetscRealPart(x[j - 1][i - 1]) > 0; 98 live_neighbors += PetscRealPart(x[j - 1][i]) > 0; 99 live_neighbors += PetscRealPart(x[j - 1][i + 1]) > 0; 100 live_neighbors += PetscRealPart(x[j][i - 1]) > 0; 101 live_neighbors += PetscRealPart(x[j][i + 1]) > 0; 102 live_neighbors += PetscRealPart(x[j + 1][i - 1]) > 0; 103 live_neighbors += PetscRealPart(x[j + 1][i]) > 0; 104 live_neighbors += PetscRealPart(x[j + 1][i + 1]) > 0; 105 if (PetscRealPart(x[j][i]) > 0) { /* Live cell */ 106 switch (live_neighbors) { 107 case 2: 108 case 3: 109 y[j][i] = 1; /* Survive */ 110 break; 111 default: y[j][i] = 0; /* Death */ 112 } 113 } else { /* Dead cell */ 114 if (live_neighbors == 3) y[j][i] = 1; /* Birth */ 115 else y[j][i] = 0; 116 } 117 } 118 } 119 PetscCall(DMDAVecRestoreArrayRead(da, Xlocal, (void *)&x)); 120 PetscCall(DMDAVecRestoreArrayWrite(da, Xglobal, &y)); 121 if (step == check_step_alive || step == check_step_dead) { 122 PetscScalar sum; 123 PetscCall(VecSum(Xglobal, &sum)); 124 if (PetscAbsScalar(sum) > 0.1) { 125 if (step == check_step_dead) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Simulation alive at step %" PetscInt_FMT "\n", step)); 126 } else if (step == check_step_alive) { 127 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Simulation dead at step %" PetscInt_FMT "\n", step)); 128 } 129 } 130 if (step % viz_interval == 0) PetscCall(VecView(Xglobal, viewer)); 131 } 132 } 133 134 PetscCall(PetscViewerDestroy(&viewer)); 135 PetscCall(VecDestroy(&Xglobal)); 136 PetscCall(VecDestroy(&Xlocal)); 137 PetscCall(DMDestroy(&da)); 138 PetscCall(PetscFinalize()); 139 return 0; 140 } 141 142 /*TEST 143 144 test: 145 requires: x 146 nsize: 2 147 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 148 149 TEST*/ 150