1 2 3 static char help[] = "Time-dependent PDE in 2d. Simplified from ex7.c for illustrating how to use TS on a structured domain. \n"; 4 /* 5 u_t = uxx + uyy 6 0 < x < 1, 0 < y < 1; 7 At t=0: u(x,y) = exp(c*r*r*r), if r=PetscSqrtReal((x-.5)*(x-.5) + (y-.5)*(y-.5)) < .125 8 u(x,y) = 0.0 if r >= .125 9 10 mpiexec -n 2 ./ex13 -da_grid_x 40 -da_grid_y 40 -ts_max_steps 2 -snes_monitor -ksp_monitor 11 mpiexec -n 1 ./ex13 -snes_fd_color -ts_monitor_draw_solution 12 mpiexec -n 2 ./ex13 -ts_type sundials -ts_monitor 13 */ 14 15 #include <petscdm.h> 16 #include <petscdmda.h> 17 #include <petscts.h> 18 19 /* 20 User-defined data structures and routines 21 */ 22 typedef struct { 23 PetscReal c; 24 } AppCtx; 25 26 extern PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*); 27 extern PetscErrorCode RHSJacobian(TS,PetscReal,Vec,Mat,Mat,void*); 28 extern PetscErrorCode FormInitialSolution(DM,Vec,void*); 29 30 int main(int argc,char **argv) 31 { 32 TS ts; /* nonlinear solver */ 33 Vec u,r; /* solution, residual vector */ 34 Mat J; /* Jacobian matrix */ 35 PetscInt steps; /* iterations for convergence */ 36 PetscErrorCode ierr; 37 DM da; 38 PetscReal ftime,dt; 39 AppCtx user; /* user-defined work context */ 40 41 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 42 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 43 Create distributed array (DMDA) to manage parallel grid and vectors 44 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 45 ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,8,8,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);CHKERRQ(ierr); 46 ierr = DMSetFromOptions(da);CHKERRQ(ierr); 47 ierr = DMSetUp(da);CHKERRQ(ierr); 48 49 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 50 Extract global vectors from DMDA; 51 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 52 ierr = DMCreateGlobalVector(da,&u);CHKERRQ(ierr); 53 ierr = VecDuplicate(u,&r);CHKERRQ(ierr); 54 55 /* Initialize user application context */ 56 user.c = -30.0; 57 58 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 59 Create timestepping solver context 60 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 61 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 62 ierr = TSSetDM(ts,da);CHKERRQ(ierr); 63 ierr = TSSetType(ts,TSBEULER);CHKERRQ(ierr); 64 ierr = TSSetRHSFunction(ts,r,RHSFunction,&user);CHKERRQ(ierr); 65 66 /* Set Jacobian */ 67 ierr = DMSetMatType(da,MATAIJ);CHKERRQ(ierr); 68 ierr = DMCreateMatrix(da,&J);CHKERRQ(ierr); 69 ierr = TSSetRHSJacobian(ts,J,J,RHSJacobian,NULL);CHKERRQ(ierr); 70 71 ftime = 1.0; 72 ierr = TSSetMaxTime(ts,ftime);CHKERRQ(ierr); 73 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 74 75 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 76 Set initial conditions 77 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 78 ierr = FormInitialSolution(da,u,&user);CHKERRQ(ierr); 79 dt = .01; 80 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); 81 82 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 83 Set runtime options 84 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 85 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 86 87 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 88 Solve nonlinear system 89 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 90 ierr = TSSolve(ts,u);CHKERRQ(ierr); 91 ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr); 92 ierr = TSGetStepNumber(ts,&steps);CHKERRQ(ierr); 93 94 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 95 Free work space. 96 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 97 ierr = MatDestroy(&J);CHKERRQ(ierr); 98 ierr = VecDestroy(&u);CHKERRQ(ierr); 99 ierr = VecDestroy(&r);CHKERRQ(ierr); 100 ierr = TSDestroy(&ts);CHKERRQ(ierr); 101 ierr = DMDestroy(&da);CHKERRQ(ierr); 102 103 ierr = PetscFinalize(); 104 return ierr; 105 } 106 /* ------------------------------------------------------------------- */ 107 /* 108 RHSFunction - Evaluates nonlinear function, F(u). 109 110 Input Parameters: 111 . ts - the TS context 112 . U - input vector 113 . ptr - optional user-defined context, as set by TSSetFunction() 114 115 Output Parameter: 116 . F - function vector 117 */ 118 PetscErrorCode RHSFunction(TS ts,PetscReal ftime,Vec U,Vec F,void *ptr) 119 { 120 /* PETSC_UNUSED AppCtx *user=(AppCtx*)ptr; */ 121 DM da; 122 PetscErrorCode ierr; 123 PetscInt i,j,Mx,My,xs,ys,xm,ym; 124 PetscReal two = 2.0,hx,hy,sx,sy; 125 PetscScalar u,uxx,uyy,**uarray,**f; 126 Vec localU; 127 128 PetscFunctionBeginUser; 129 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 130 ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr); 131 ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr); 132 133 hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx); 134 hy = 1.0/(PetscReal)(My-1); sy = 1.0/(hy*hy); 135 136 /* 137 Scatter ghost points to local vector,using the 2-step process 138 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 139 By placing code between these two statements, computations can be 140 done while messages are in transition. 141 */ 142 ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 143 ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 144 145 /* Get pointers to vector data */ 146 ierr = DMDAVecGetArrayRead(da,localU,&uarray);CHKERRQ(ierr); 147 ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr); 148 149 /* Get local grid boundaries */ 150 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 151 152 /* Compute function over the locally owned part of the grid */ 153 for (j=ys; j<ys+ym; j++) { 154 for (i=xs; i<xs+xm; i++) { 155 if (i == 0 || j == 0 || i == Mx-1 || j == My-1) { 156 f[j][i] = uarray[j][i]; 157 continue; 158 } 159 u = uarray[j][i]; 160 uxx = (-two*u + uarray[j][i-1] + uarray[j][i+1])*sx; 161 uyy = (-two*u + uarray[j-1][i] + uarray[j+1][i])*sy; 162 f[j][i] = uxx + uyy; 163 } 164 } 165 166 /* Restore vectors */ 167 ierr = DMDAVecRestoreArrayRead(da,localU,&uarray);CHKERRQ(ierr); 168 ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr); 169 ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr); 170 ierr = PetscLogFlops(11.0*ym*xm);CHKERRQ(ierr); 171 PetscFunctionReturn(0); 172 } 173 174 /* --------------------------------------------------------------------- */ 175 /* 176 RHSJacobian - User-provided routine to compute the Jacobian of 177 the nonlinear right-hand-side function of the ODE. 178 179 Input Parameters: 180 ts - the TS context 181 t - current time 182 U - global input vector 183 dummy - optional user-defined context, as set by TSetRHSJacobian() 184 185 Output Parameters: 186 J - Jacobian matrix 187 Jpre - optionally different preconditioning matrix 188 str - flag indicating matrix structure 189 */ 190 PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat J,Mat Jpre,void *ctx) 191 { 192 PetscErrorCode ierr; 193 DM da; 194 DMDALocalInfo info; 195 PetscInt i,j; 196 PetscReal hx,hy,sx,sy; 197 198 PetscFunctionBeginUser; 199 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 200 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 201 hx = 1.0/(PetscReal)(info.mx-1); sx = 1.0/(hx*hx); 202 hy = 1.0/(PetscReal)(info.my-1); sy = 1.0/(hy*hy); 203 for (j=info.ys; j<info.ys+info.ym; j++) { 204 for (i=info.xs; i<info.xs+info.xm; i++) { 205 PetscInt nc = 0; 206 MatStencil row,col[5]; 207 PetscScalar val[5]; 208 row.i = i; row.j = j; 209 if (i == 0 || j == 0 || i == info.mx-1 || j == info.my-1) { 210 col[nc].i = i; col[nc].j = j; val[nc++] = 1.0; 211 } else { 212 col[nc].i = i-1; col[nc].j = j; val[nc++] = sx; 213 col[nc].i = i+1; col[nc].j = j; val[nc++] = sx; 214 col[nc].i = i; col[nc].j = j-1; val[nc++] = sy; 215 col[nc].i = i; col[nc].j = j+1; val[nc++] = sy; 216 col[nc].i = i; col[nc].j = j; val[nc++] = -2*sx - 2*sy; 217 } 218 ierr = MatSetValuesStencil(Jpre,1,&row,nc,col,val,INSERT_VALUES);CHKERRQ(ierr); 219 } 220 } 221 ierr = MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 222 ierr = MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 223 if (J != Jpre) { 224 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 225 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 226 } 227 PetscFunctionReturn(0); 228 } 229 230 /* ------------------------------------------------------------------- */ 231 PetscErrorCode FormInitialSolution(DM da,Vec U,void* ptr) 232 { 233 AppCtx *user=(AppCtx*)ptr; 234 PetscReal c=user->c; 235 PetscErrorCode ierr; 236 PetscInt i,j,xs,ys,xm,ym,Mx,My; 237 PetscScalar **u; 238 PetscReal hx,hy,x,y,r; 239 240 PetscFunctionBeginUser; 241 ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr); 242 243 hx = 1.0/(PetscReal)(Mx-1); 244 hy = 1.0/(PetscReal)(My-1); 245 246 /* Get pointers to vector data */ 247 ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr); 248 249 /* Get local grid boundaries */ 250 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 251 252 /* Compute function over the locally owned part of the grid */ 253 for (j=ys; j<ys+ym; j++) { 254 y = j*hy; 255 for (i=xs; i<xs+xm; i++) { 256 x = i*hx; 257 r = PetscSqrtReal((x-.5)*(x-.5) + (y-.5)*(y-.5)); 258 if (r < .125) u[j][i] = PetscExpReal(c*r*r*r); 259 else u[j][i] = 0.0; 260 } 261 } 262 263 /* Restore vectors */ 264 ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr); 265 PetscFunctionReturn(0); 266 } 267 268 /*TEST 269 270 test: 271 args: -ts_max_steps 5 -ts_monitor 272 273 test: 274 suffix: 2 275 args: -ts_max_steps 5 -ts_monitor 276 277 test: 278 suffix: 3 279 args: -ts_max_steps 5 -snes_fd_color -ts_monitor 280 281 TEST*/ 282