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 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++) 66 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, &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, &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) { 128 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Simulation alive at step %D\n",step)); 129 } 130 } else if (step == check_step_alive) { 131 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Simulation dead at step %D\n",step)); 132 } 133 } 134 if (step % viz_interval == 0) { 135 PetscCall(VecView(Xglobal, viewer)); 136 } 137 } 138 } 139 140 PetscCall(PetscViewerDestroy(&viewer)); 141 PetscCall(VecDestroy(&Xglobal)); 142 PetscCall(VecDestroy(&Xlocal)); 143 PetscCall(DMDestroy(&da)); 144 PetscCall(PetscFinalize()); 145 return 0; 146 } 147 148 /*TEST 149 150 test: 151 requires: x 152 nsize: 2 153 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 154 155 TEST*/ 156