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 { 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, (char *)0, 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 151 TEST*/ 152