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