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