1 static const char help[] = "Time-dependent PDE in 1d. Simplified from ex15.c for illustrating how to solve DAEs. \n"; 2 /* 3 u_t = uxx 4 0 < x < 1; 5 At t=0: u(x) = exp(c*r*r*r), if r=PetscSqrtReal((x-.5)*(x-.5)) < .125 6 u(x) = 0.0 if r >= .125 7 8 9 Boundary conditions: 10 Dirichlet BC: 11 At x=0, x=1, u = 0.0 12 13 Neumann BC: 14 At x=0, x=1: du(x,t)/dx = 0 15 16 mpiexec -n 2 ./ex17 -da_grid_x 40 -ts_max_steps 2 -snes_monitor -ksp_monitor 17 ./ex17 -da_grid_x 40 -monitor_solution 18 ./ex17 -da_grid_x 100 -ts_type theta -ts_theta_theta 0.5 # Midpoint is not L-stable 19 ./ex17 -jac_type fd_coloring -da_grid_x 500 -boundary 1 20 ./ex17 -da_grid_x 100 -ts_type gl -ts_adapt_type none -ts_max_steps 2 21 */ 22 23 #include <petscdm.h> 24 #include <petscdmda.h> 25 #include <petscts.h> 26 27 typedef enum {JACOBIAN_ANALYTIC,JACOBIAN_FD_COLORING,JACOBIAN_FD_FULL} JacobianType; 28 static const char *const JacobianTypes[] = {"analytic","fd_coloring","fd_full","JacobianType","fd_",0}; 29 30 /* 31 User-defined data structures and routines 32 */ 33 typedef struct { 34 PetscReal c; 35 PetscInt boundary; /* Type of boundary condition */ 36 PetscBool viewJacobian; 37 } AppCtx; 38 39 static PetscErrorCode FormIFunction(TS,PetscReal,Vec,Vec,Vec,void*); 40 static PetscErrorCode FormIJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*); 41 static PetscErrorCode FormInitialSolution(TS,Vec,void*); 42 43 int main(int argc,char **argv) 44 { 45 TS ts; /* nonlinear solver */ 46 Vec u; /* solution, residual vectors */ 47 Mat J; /* Jacobian matrix */ 48 PetscInt nsteps; 49 PetscReal vmin,vmax,norm; 50 PetscErrorCode ierr; 51 DM da; 52 PetscReal ftime,dt; 53 AppCtx user; /* user-defined work context */ 54 JacobianType jacType; 55 56 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 57 58 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 59 Create distributed array (DMDA) to manage parallel grid and vectors 60 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 61 ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,11,1,1,NULL,&da);CHKERRQ(ierr); 62 ierr = DMSetFromOptions(da);CHKERRQ(ierr); 63 ierr = DMSetUp(da);CHKERRQ(ierr); 64 65 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 66 Extract global vectors from DMDA; 67 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 68 ierr = DMCreateGlobalVector(da,&u);CHKERRQ(ierr); 69 70 /* Initialize user application context */ 71 user.c = -30.0; 72 user.boundary = 0; /* 0: Dirichlet BC; 1: Neumann BC */ 73 user.viewJacobian = PETSC_FALSE; 74 75 ierr = PetscOptionsGetInt(NULL,NULL,"-boundary",&user.boundary,NULL);CHKERRQ(ierr); 76 ierr = PetscOptionsHasName(NULL,NULL,"-viewJacobian",&user.viewJacobian);CHKERRQ(ierr); 77 78 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 79 Create timestepping solver context 80 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 81 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 82 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 83 ierr = TSSetType(ts,TSTHETA);CHKERRQ(ierr); 84 ierr = TSThetaSetTheta(ts,1.0);CHKERRQ(ierr); /* Make the Theta method behave like backward Euler */ 85 ierr = TSSetIFunction(ts,NULL,FormIFunction,&user);CHKERRQ(ierr); 86 87 ierr = DMSetMatType(da,MATAIJ);CHKERRQ(ierr); 88 ierr = DMCreateMatrix(da,&J);CHKERRQ(ierr); 89 jacType = JACOBIAN_ANALYTIC; /* use user-provide Jacobian */ 90 91 ierr = TSSetDM(ts,da);CHKERRQ(ierr); /* Use TSGetDM() to access. Setting here allows easy use of geometric multigrid. */ 92 93 ftime = 1.0; 94 ierr = TSSetMaxTime(ts,ftime);CHKERRQ(ierr); 95 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 96 97 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 98 Set initial conditions 99 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 100 ierr = FormInitialSolution(ts,u,&user);CHKERRQ(ierr); 101 ierr = TSSetSolution(ts,u);CHKERRQ(ierr); 102 dt = .01; 103 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); 104 105 106 /* Use slow fd Jacobian or fast fd Jacobian with colorings. 107 Note: this requirs snes which is not created until TSSetUp()/TSSetFromOptions() is called */ 108 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Options for Jacobian evaluation",NULL);CHKERRQ(ierr); 109 ierr = PetscOptionsEnum("-jac_type","Type of Jacobian","",JacobianTypes,(PetscEnum)jacType,(PetscEnum*)&jacType,0);CHKERRQ(ierr); 110 ierr = PetscOptionsEnd();CHKERRQ(ierr); 111 if (jacType == JACOBIAN_ANALYTIC) { 112 ierr = TSSetIJacobian(ts,J,J,FormIJacobian,&user);CHKERRQ(ierr); 113 } else if (jacType == JACOBIAN_FD_COLORING) { 114 SNES snes; 115 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 116 ierr = SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,0);CHKERRQ(ierr); 117 } else if (jacType == JACOBIAN_FD_FULL) { 118 SNES snes; 119 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 120 ierr = SNESSetJacobian(snes,J,J,SNESComputeJacobianDefault,&user);CHKERRQ(ierr); 121 } 122 123 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 124 Set runtime options 125 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 126 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 127 128 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 129 Integrate ODE system 130 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 131 ierr = TSSolve(ts,u);CHKERRQ(ierr); 132 133 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 134 Compute diagnostics of the solution 135 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 136 ierr = VecNorm(u,NORM_1,&norm);CHKERRQ(ierr); 137 ierr = VecMax(u,NULL,&vmax);CHKERRQ(ierr); 138 ierr = VecMin(u,NULL,&vmin);CHKERRQ(ierr); 139 ierr = TSGetStepNumber(ts,&nsteps);CHKERRQ(ierr); 140 ierr = TSGetTime(ts,&ftime);CHKERRQ(ierr); 141 ierr = PetscPrintf(PETSC_COMM_WORLD,"timestep %D: time %g, solution norm %g, max %g, min %g\n",nsteps,(double)ftime,(double)norm,(double)vmax,(double)vmin);CHKERRQ(ierr); 142 143 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 144 Free work space. 145 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 146 ierr = MatDestroy(&J);CHKERRQ(ierr); 147 ierr = VecDestroy(&u);CHKERRQ(ierr); 148 ierr = TSDestroy(&ts);CHKERRQ(ierr); 149 ierr = DMDestroy(&da);CHKERRQ(ierr); 150 ierr = PetscFinalize(); 151 return ierr; 152 } 153 /* ------------------------------------------------------------------- */ 154 static PetscErrorCode FormIFunction(TS ts,PetscReal ftime,Vec U,Vec Udot,Vec F,void *ptr) 155 { 156 AppCtx *user=(AppCtx*)ptr; 157 DM da; 158 PetscErrorCode ierr; 159 PetscInt i,Mx,xs,xm; 160 PetscReal hx,sx; 161 PetscScalar *u,*udot,*f; 162 Vec localU; 163 164 PetscFunctionBeginUser; 165 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 166 ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr); 167 ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr); 168 169 hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx); 170 171 /* 172 Scatter ghost points to local vector,using the 2-step process 173 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 174 By placing code between these two statements, computations can be 175 done while messages are in transition. 176 */ 177 ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 178 ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 179 180 /* Get pointers to vector data */ 181 ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr); 182 ierr = DMDAVecGetArrayRead(da,Udot,&udot);CHKERRQ(ierr); 183 ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr); 184 185 /* Get local grid boundaries */ 186 ierr = DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr); 187 188 /* Compute function over the locally owned part of the grid */ 189 for (i=xs; i<xs+xm; i++) { 190 if (user->boundary == 0) { /* Dirichlet BC */ 191 if (i == 0 || i == Mx-1) f[i] = u[i]; /* F = U */ 192 else f[i] = udot[i] + (2.*u[i] - u[i-1] - u[i+1])*sx; 193 } else { /* Neumann BC */ 194 if (i == 0) f[i] = u[0] - u[1]; 195 else if (i == Mx-1) f[i] = u[i] - u[i-1]; 196 else f[i] = udot[i] + (2.*u[i] - u[i-1] - u[i+1])*sx; 197 } 198 } 199 200 /* Restore vectors */ 201 ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr); 202 ierr = DMDAVecRestoreArrayRead(da,Udot,&udot);CHKERRQ(ierr); 203 ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr); 204 ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr); 205 PetscFunctionReturn(0); 206 } 207 208 /* --------------------------------------------------------------------- */ 209 /* 210 IJacobian - Compute IJacobian = dF/dU + a dF/dUdot 211 */ 212 PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat J,Mat Jpre,void *ctx) 213 { 214 PetscErrorCode ierr; 215 PetscInt i,rstart,rend,Mx; 216 PetscReal hx,sx; 217 AppCtx *user = (AppCtx*)ctx; 218 DM da; 219 MatStencil col[3],row; 220 PetscInt nc; 221 PetscScalar vals[3]; 222 223 PetscFunctionBeginUser; 224 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 225 ierr = MatGetOwnershipRange(Jpre,&rstart,&rend);CHKERRQ(ierr); 226 ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr); 227 hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx); 228 for (i=rstart; i<rend; i++) { 229 nc = 0; 230 row.i = i; 231 if (user->boundary == 0 && (i == 0 || i == Mx-1)) { 232 col[nc].i = i; vals[nc++] = 1.0; 233 } else if (user->boundary > 0 && i == 0) { /* Left Neumann */ 234 col[nc].i = i; vals[nc++] = 1.0; 235 col[nc].i = i+1; vals[nc++] = -1.0; 236 } else if (user->boundary > 0 && i == Mx-1) { /* Right Neumann */ 237 col[nc].i = i-1; vals[nc++] = -1.0; 238 col[nc].i = i; vals[nc++] = 1.0; 239 } else { /* Interior */ 240 col[nc].i = i-1; vals[nc++] = -1.0*sx; 241 col[nc].i = i; vals[nc++] = 2.0*sx + a; 242 col[nc].i = i+1; vals[nc++] = -1.0*sx; 243 } 244 ierr = MatSetValuesStencil(Jpre,1,&row,nc,col,vals,INSERT_VALUES);CHKERRQ(ierr); 245 } 246 247 ierr = MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 248 ierr = MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 249 if (J != Jpre) { 250 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 251 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 252 } 253 if (user->viewJacobian) { 254 ierr = PetscPrintf(PETSC_COMM_WORLD,"Jpre:\n");CHKERRQ(ierr); 255 ierr = MatView(Jpre,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 256 } 257 PetscFunctionReturn(0); 258 } 259 260 /* ------------------------------------------------------------------- */ 261 PetscErrorCode FormInitialSolution(TS ts,Vec U,void *ptr) 262 { 263 AppCtx *user=(AppCtx*)ptr; 264 PetscReal c =user->c; 265 DM da; 266 PetscErrorCode ierr; 267 PetscInt i,xs,xm,Mx; 268 PetscScalar *u; 269 PetscReal hx,x,r; 270 271 PetscFunctionBeginUser; 272 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 273 ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr); 274 275 hx = 1.0/(PetscReal)(Mx-1); 276 277 /* Get pointers to vector data */ 278 ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr); 279 280 /* Get local grid boundaries */ 281 ierr = DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr); 282 283 /* Compute function over the locally owned part of the grid */ 284 for (i=xs; i<xs+xm; i++) { 285 x = i*hx; 286 r = PetscSqrtReal((x-.5)*(x-.5)); 287 if (r < .125) u[i] = PetscExpReal(c*r*r*r); 288 else u[i] = 0.0; 289 } 290 291 /* Restore vectors */ 292 ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr); 293 PetscFunctionReturn(0); 294 } 295 296 /*TEST 297 298 test: 299 requires: !single 300 args: -da_grid_x 40 -ts_max_steps 2 -snes_monitor_short -ksp_monitor_short -ts_monitor 301 302 test: 303 suffix: 2 304 requires: !single 305 args: -da_grid_x 100 -ts_type theta -ts_theta_theta 0.5 306 307 308 TEST*/ 309 310