static char help[] ="Solves a time-dependent nonlinear PDE.\n"; /* ------------------------------------------------------------------------ This program solves the two-dimensional time-dependent Bratu problem u_t = u_xx + u_yy + \lambda*exp(u), on the domain 0 <= x,y <= 1, with the boundary conditions u(t,0,y) = 0, u_x(t,1,y) = 0, u(t,x,0) = 0, u_x(t,x,1) = 0, and the initial condition u(0,x,y) = 0. We discretize the right-hand side using finite differences with uniform grid spacings hx,hy: u_xx = (u_{i+1} - 2u_{i} + u_{i-1})/(hx^2) u_yy = (u_{j+1} - 2u_{j} + u_{j-1})/(hy^2) ------------------------------------------------------------------------- */ #include #include /* User-defined application context - contains data needed by the application-provided call-back routines. */ typedef struct { PetscReal lambda; } AppCtx; /* FormIFunctionLocal - Evaluates nonlinear implicit function on local process patch */ static PetscErrorCode FormIFunctionLocal(DMDALocalInfo *info,PetscReal t,PetscScalar **x,PetscScalar **xdot,PetscScalar **f,AppCtx *app) { PetscInt i,j; PetscReal lambda,hx,hy; PetscScalar ut,u,ue,uw,un,us,uxx,uyy; PetscFunctionBeginUser; lambda = app->lambda; hx = 1.0/(PetscReal)(info->mx-1); hy = 1.0/(PetscReal)(info->my-1); /* Compute RHS function over the locally owned part of the grid */ for (j=info->ys; jys+info->ym; j++) { for (i=info->xs; ixs+info->xm; i++) { if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) { /* boundary points */ f[j][i] = x[j][i] - (PetscReal)0; } else { /* interior points */ ut = xdot[j][i]; u = x[j][i]; uw = x[j][i-1]; ue = x[j][i+1]; un = x[j+1][i]; us = x[j-1][i]; uxx = (uw - 2.0*u + ue)/(hx*hx); uyy = (un - 2.0*u + us)/(hy*hy); f[j][i] = ut - uxx - uyy - lambda*PetscExpScalar(u); } } } PetscFunctionReturn(0); } /* FormIJacobianLocal - Evaluates implicit Jacobian matrix on local process patch */ static PetscErrorCode FormIJacobianLocal(DMDALocalInfo *info,PetscReal t,PetscScalar **x,PetscScalar **xdot,PetscScalar shift,Mat jac,Mat jacpre,AppCtx *app) { PetscErrorCode ierr; PetscInt i,j,k; MatStencil col[5],row; PetscScalar v[5],lambda,hx,hy; PetscFunctionBeginUser; lambda = app->lambda; hx = 1.0/(PetscReal)(info->mx-1); hy = 1.0/(PetscReal)(info->my-1); /* Compute Jacobian entries for the locally owned part of the grid */ for (j=info->ys; jys+info->ym; j++) { for (i=info->xs; ixs+info->xm; i++) { row.j = j; row.i = i; k = 0; if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) { /* boundary points */ v[0] = 1.0; ierr = MatSetValuesStencil(jacpre,1,&row,1,&row,v,INSERT_VALUES);CHKERRQ(ierr); } else { /* interior points */ v[k] = -1.0/(hy*hy); col[k].j = j-1; col[k].i = i; k++; v[k] = -1.0/(hx*hx); col[k].j = j; col[k].i = i-1; k++; v[k] = shift + 2.0/(hx*hx) + 2.0/(hy*hy) - lambda*PetscExpScalar(x[j][i]); col[k].j = j; col[k].i = i; k++; v[k] = -1.0/(hx*hx); col[k].j = j; col[k].i = i+1; k++; v[k] = -1.0/(hy*hy); col[k].j = j+1; col[k].i = i; k++; ierr = MatSetValuesStencil(jacpre,1,&row,k,col,v,INSERT_VALUES);CHKERRQ(ierr); } } } /* Assemble matrix */ ierr = MatAssemblyBegin(jacpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(jacpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); PetscFunctionReturn(0); } int main(int argc,char **argv) { TS ts; /* ODE integrator */ DM da; /* DM context */ Vec U; /* solution vector */ DMBoundaryType bt = DM_BOUNDARY_NONE; DMDAStencilType st = DMDA_STENCIL_STAR; PetscInt sw = 1; PetscInt N = 17; PetscInt n = PETSC_DECIDE; AppCtx app; PetscErrorCode ierr; /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Initialize program - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex21 options","");CHKERRQ(ierr); { app.lambda = 6.8; app.lambda = 6.0; ierr = PetscOptionsReal("-lambda","","",app.lambda,&app.lambda,NULL);CHKERRQ(ierr); } ierr = PetscOptionsEnd();CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Create DM context - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = DMDACreate2d(PETSC_COMM_WORLD,bt,bt,st,N,N,n,n,1,sw,NULL,NULL,&da);CHKERRQ(ierr); ierr = DMSetFromOptions(da);CHKERRQ(ierr); ierr = DMSetUp(da);CHKERRQ(ierr); ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0,1.0);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Create timestepping solver context - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); ierr = TSSetDM(ts,da);CHKERRQ(ierr); ierr = DMDestroy(&da);CHKERRQ(ierr); ierr = TSGetDM(ts,&da);CHKERRQ(ierr); ierr = DMDATSSetIFunctionLocal(da,INSERT_VALUES,(DMDATSIFunctionLocal)FormIFunctionLocal,&app);CHKERRQ(ierr); ierr = DMDATSSetIJacobianLocal(da,(DMDATSIJacobianLocal)FormIJacobianLocal,&app);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Set solver options - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = TSSetType(ts,TSBDF);CHKERRQ(ierr); ierr = TSSetTimeStep(ts,1e-4);CHKERRQ(ierr); ierr = TSSetMaxTime(ts,1.0);CHKERRQ(ierr); ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); ierr = TSSetFromOptions(ts);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Set initial conditions - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = TSGetDM(ts,&da);CHKERRQ(ierr); ierr = DMCreateGlobalVector(da,&U);CHKERRQ(ierr); ierr = VecSet(U,0.0);CHKERRQ(ierr); ierr = TSSetSolution(ts,U);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Run timestepping solver - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = TSSolve(ts,U);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - All PETSc objects should be destroyed when they are no longer needed. - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = VecDestroy(&U);CHKERRQ(ierr); ierr = TSDestroy(&ts);CHKERRQ(ierr); ierr = PetscFinalize(); return ierr; } /*TEST testset: requires: !single !complex 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 filter: grep -v "total number of" test: suffix: 1_bdf_ngmres_fas_ms args: -prefix_push npc_fas_levels_ -snes_type ms -snes_max_it 5 -ksp_type preonly -prefix_pop test: suffix: 2_bdf_ngmres_fas_ms args: -prefix_push npc_fas_levels_ -snes_type ms -snes_max_it 5 -ksp_type preonly -prefix_pop nsize: 2 test: suffix: 1_bdf_ngmres_fas_ngs args: -prefix_push npc_fas_levels_ -snes_type ngs -snes_max_it 5 -prefix_pop test: suffix: 2_bdf_ngmres_fas_ngs args: -prefix_push npc_fas_levels_ -snes_type ngs -snes_max_it 5 -prefix_pop nsize: 2 TEST*/