1 static char help[] ="Solves a time-dependent nonlinear PDE.\n"; 2 3 /* ------------------------------------------------------------------------ 4 5 This program solves the two-dimensional time-dependent Bratu problem 6 u_t = u_xx + u_yy + \lambda*exp(u), 7 on the domain 0 <= x,y <= 1, 8 with the boundary conditions 9 u(t,0,y) = 0, u_x(t,1,y) = 0, 10 u(t,x,0) = 0, u_x(t,x,1) = 0, 11 and the initial condition 12 u(0,x,y) = 0. 13 We discretize the right-hand side using finite differences with 14 uniform grid spacings hx,hy: 15 u_xx = (u_{i+1} - 2u_{i} + u_{i-1})/(hx^2) 16 u_yy = (u_{j+1} - 2u_{j} + u_{j-1})/(hy^2) 17 18 ------------------------------------------------------------------------- */ 19 20 #include <petscdmda.h> 21 #include <petscts.h> 22 23 /* 24 User-defined application context - contains data needed by the 25 application-provided call-back routines. 26 */ 27 typedef struct { 28 PetscReal lambda; 29 } AppCtx; 30 31 /* 32 FormIFunctionLocal - Evaluates nonlinear implicit function on local process patch 33 */ 34 static PetscErrorCode FormIFunctionLocal(DMDALocalInfo *info,PetscReal t,PetscScalar **x,PetscScalar **xdot,PetscScalar **f,AppCtx *app) 35 { 36 PetscInt i,j; 37 PetscReal lambda,hx,hy; 38 PetscScalar ut,u,ue,uw,un,us,uxx,uyy; 39 40 PetscFunctionBeginUser; 41 lambda = app->lambda; 42 hx = 1.0/(PetscReal)(info->mx-1); 43 hy = 1.0/(PetscReal)(info->my-1); 44 45 /* 46 Compute RHS function over the locally owned part of the grid 47 */ 48 for (j=info->ys; j<info->ys+info->ym; j++) { 49 for (i=info->xs; i<info->xs+info->xm; i++) { 50 if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) { 51 /* boundary points */ 52 f[j][i] = x[j][i] - (PetscReal)0; 53 } else { 54 /* interior points */ 55 ut = xdot[j][i]; 56 u = x[j][i]; 57 uw = x[j][i-1]; 58 ue = x[j][i+1]; 59 un = x[j+1][i]; 60 us = x[j-1][i]; 61 62 uxx = (uw - 2.0*u + ue)/(hx*hx); 63 uyy = (un - 2.0*u + us)/(hy*hy); 64 f[j][i] = ut - uxx - uyy - lambda*PetscExpScalar(u); 65 } 66 } 67 } 68 PetscFunctionReturn(0); 69 } 70 71 /* 72 FormIJacobianLocal - Evaluates implicit Jacobian matrix on local process patch 73 */ 74 static PetscErrorCode FormIJacobianLocal(DMDALocalInfo *info,PetscReal t,PetscScalar **x,PetscScalar **xdot,PetscScalar shift,Mat jac,Mat jacpre,AppCtx *app) 75 { 76 PetscErrorCode ierr; 77 PetscInt i,j,k; 78 MatStencil col[5],row; 79 PetscScalar v[5],lambda,hx,hy; 80 81 PetscFunctionBeginUser; 82 lambda = app->lambda; 83 hx = 1.0/(PetscReal)(info->mx-1); 84 hy = 1.0/(PetscReal)(info->my-1); 85 86 /* 87 Compute Jacobian entries for the locally owned part of the grid 88 */ 89 for (j=info->ys; j<info->ys+info->ym; j++) { 90 for (i=info->xs; i<info->xs+info->xm; i++) { 91 row.j = j; row.i = i; k = 0; 92 if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) { 93 /* boundary points */ 94 v[0] = 1.0; 95 ierr = MatSetValuesStencil(jacpre,1,&row,1,&row,v,INSERT_VALUES);CHKERRQ(ierr); 96 } else { 97 /* interior points */ 98 v[k] = -1.0/(hy*hy); col[k].j = j-1; col[k].i = i; k++; 99 v[k] = -1.0/(hx*hx); col[k].j = j; col[k].i = i-1; k++; 100 101 v[k] = shift + 2.0/(hx*hx) + 2.0/(hy*hy) - lambda*PetscExpScalar(x[j][i]); 102 col[k].j = j; col[k].i = i; k++; 103 104 v[k] = -1.0/(hx*hx); col[k].j = j; col[k].i = i+1; k++; 105 v[k] = -1.0/(hy*hy); col[k].j = j+1; col[k].i = i; k++; 106 107 ierr = MatSetValuesStencil(jacpre,1,&row,k,col,v,INSERT_VALUES);CHKERRQ(ierr); 108 } 109 } 110 } 111 112 /* 113 Assemble matrix 114 */ 115 ierr = MatAssemblyBegin(jacpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 116 ierr = MatAssemblyEnd(jacpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 117 PetscFunctionReturn(0); 118 } 119 120 121 int main(int argc,char **argv) 122 { 123 TS ts; /* ODE integrator */ 124 DM da; /* DM context */ 125 Vec U; /* solution vector */ 126 DMBoundaryType bt = DM_BOUNDARY_NONE; 127 DMDAStencilType st = DMDA_STENCIL_STAR; 128 PetscInt sw = 1; 129 PetscInt N = 17; 130 PetscInt n = PETSC_DECIDE; 131 AppCtx app; 132 PetscErrorCode ierr; 133 134 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 135 Initialize program 136 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 137 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; 138 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex21 options","");CHKERRQ(ierr); 139 { 140 app.lambda = 6.8; app.lambda = 6.0; 141 ierr = PetscOptionsReal("-lambda","","",app.lambda,&app.lambda,NULL);CHKERRQ(ierr); 142 } 143 ierr = PetscOptionsEnd();CHKERRQ(ierr); 144 145 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 146 Create DM context 147 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 148 ierr = DMDACreate2d(PETSC_COMM_WORLD,bt,bt,st,N,N,n,n,1,sw,NULL,NULL,&da);CHKERRQ(ierr); 149 ierr = DMSetFromOptions(da);CHKERRQ(ierr); 150 ierr = DMSetUp(da);CHKERRQ(ierr); 151 ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0,1.0);CHKERRQ(ierr); 152 153 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 154 Create timestepping solver context 155 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 156 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 157 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 158 ierr = TSSetDM(ts,da);CHKERRQ(ierr); 159 ierr = DMDestroy(&da);CHKERRQ(ierr); 160 161 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 162 ierr = DMDATSSetIFunctionLocal(da,INSERT_VALUES,(DMDATSIFunctionLocal)FormIFunctionLocal,&app);CHKERRQ(ierr); 163 ierr = DMDATSSetIJacobianLocal(da,(DMDATSIJacobianLocal)FormIJacobianLocal,&app);CHKERRQ(ierr); 164 165 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 166 Set solver options 167 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 168 ierr = TSSetType(ts,TSBDF);CHKERRQ(ierr); 169 ierr = TSSetTimeStep(ts,1e-4);CHKERRQ(ierr); 170 ierr = TSSetMaxTime(ts,1.0);CHKERRQ(ierr); 171 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 172 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 173 174 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 175 Set initial conditions 176 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 177 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 178 ierr = DMCreateGlobalVector(da,&U);CHKERRQ(ierr); 179 ierr = VecSet(U,0.0);CHKERRQ(ierr); 180 ierr = TSSetSolution(ts,U);CHKERRQ(ierr); 181 182 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 183 Run timestepping solver 184 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 185 ierr = TSSolve(ts,U);CHKERRQ(ierr); 186 187 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 188 All PETSc objects should be destroyed when they are no longer needed. 189 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 190 ierr = VecDestroy(&U);CHKERRQ(ierr); 191 ierr = TSDestroy(&ts);CHKERRQ(ierr); 192 ierr = PetscFinalize(); 193 return ierr; 194 } 195 196 /*TEST 197 198 testset: 199 requires: !single !complex 200 args: -da_grid_x 5 -da_grid_y 5 -da_refine 2 -dm_view -ts_type bdf -ts_adapt_type none -ts_dt 1e-3 -ts_monitor -ts_max_steps 5 -ts_view -snes_rtol 1e-6 -snes_type ngmres -npc_snes_type fas 201 filter: grep -v "total number of" 202 test: 203 suffix: 1_bdf_ngmres_fas_ms 204 args: -prefix_push npc_fas_levels_ -snes_type ms -snes_max_it 5 -ksp_type preonly -prefix_pop 205 test: 206 suffix: 2_bdf_ngmres_fas_ms 207 args: -prefix_push npc_fas_levels_ -snes_type ms -snes_max_it 5 -ksp_type preonly -prefix_pop 208 nsize: 2 209 test: 210 suffix: 1_bdf_ngmres_fas_ngs 211 args: -prefix_push npc_fas_levels_ -snes_type ngs -snes_max_it 5 -prefix_pop 212 test: 213 suffix: 2_bdf_ngmres_fas_ngs 214 args: -prefix_push npc_fas_levels_ -snes_type ngs -snes_max_it 5 -prefix_pop 215 nsize: 2 216 217 TEST*/ 218