1c4762a1bSJed Brown 2c4762a1bSJed Brown static char help[] ="Solves a time-dependent nonlinear PDE. Uses implicit\n\ 3c4762a1bSJed Brown timestepping. Runtime options include:\n\ 4c4762a1bSJed Brown -M <xg>, where <xg> = number of grid points\n\ 5c4762a1bSJed Brown -debug : Activate debugging printouts\n\ 6c4762a1bSJed Brown -nox : Deactivate x-window graphics\n\n"; 7c4762a1bSJed Brown 8c4762a1bSJed Brown /* 9c4762a1bSJed Brown Concepts: TS^time-dependent nonlinear problems 10c4762a1bSJed Brown Processors: n 11c4762a1bSJed Brown */ 12c4762a1bSJed Brown 13c4762a1bSJed Brown /* ------------------------------------------------------------------------ 14c4762a1bSJed Brown 15c4762a1bSJed Brown This program solves the PDE 16c4762a1bSJed Brown 17c4762a1bSJed Brown u * u_xx 18c4762a1bSJed Brown u_t = --------- 19c4762a1bSJed Brown 2*(t+1)^2 20c4762a1bSJed Brown 21c4762a1bSJed Brown on the domain 0 <= x <= 1, with boundary conditions 22c4762a1bSJed Brown u(t,0) = t + 1, u(t,1) = 2*t + 2, 23c4762a1bSJed Brown and initial condition 24c4762a1bSJed Brown u(0,x) = 1 + x*x. 25c4762a1bSJed Brown 26c4762a1bSJed Brown The exact solution is: 27c4762a1bSJed Brown u(t,x) = (1 + x*x) * (1 + t) 28c4762a1bSJed Brown 29c4762a1bSJed Brown Note that since the solution is linear in time and quadratic in x, 30c4762a1bSJed Brown the finite difference scheme actually computes the "exact" solution. 31c4762a1bSJed Brown 32c4762a1bSJed Brown We use by default the backward Euler method. 33c4762a1bSJed Brown 34c4762a1bSJed Brown ------------------------------------------------------------------------- */ 35c4762a1bSJed Brown 36c4762a1bSJed Brown /* 37c4762a1bSJed Brown Include "petscts.h" to use the PETSc timestepping routines. Note that 38c4762a1bSJed Brown this file automatically includes "petscsys.h" and other lower-level 39c4762a1bSJed Brown PETSc include files. 40c4762a1bSJed Brown 41c4762a1bSJed Brown Include the "petscdmda.h" to allow us to use the distributed array data 42c4762a1bSJed Brown structures to manage the parallel grid. 43c4762a1bSJed Brown */ 44c4762a1bSJed Brown #include <petscts.h> 45c4762a1bSJed Brown #include <petscdm.h> 46c4762a1bSJed Brown #include <petscdmda.h> 47c4762a1bSJed Brown #include <petscdraw.h> 48c4762a1bSJed Brown 49c4762a1bSJed Brown /* 50c4762a1bSJed Brown User-defined application context - contains data needed by the 51c4762a1bSJed Brown application-provided callback routines. 52c4762a1bSJed Brown */ 53c4762a1bSJed Brown typedef struct { 54c4762a1bSJed Brown MPI_Comm comm; /* communicator */ 55c4762a1bSJed Brown DM da; /* distributed array data structure */ 56c4762a1bSJed Brown Vec localwork; /* local ghosted work vector */ 57c4762a1bSJed Brown Vec u_local; /* local ghosted approximate solution vector */ 58c4762a1bSJed Brown Vec solution; /* global exact solution vector */ 59c4762a1bSJed Brown PetscInt m; /* total number of grid points */ 60c4762a1bSJed Brown PetscReal h; /* mesh width: h = 1/(m-1) */ 61c4762a1bSJed Brown PetscBool debug; /* flag (1 indicates activation of debugging printouts) */ 62c4762a1bSJed Brown } AppCtx; 63c4762a1bSJed Brown 64c4762a1bSJed Brown /* 65c4762a1bSJed Brown User-defined routines, provided below. 66c4762a1bSJed Brown */ 67c4762a1bSJed Brown extern PetscErrorCode InitialConditions(Vec,AppCtx*); 68c4762a1bSJed Brown extern PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*); 69c4762a1bSJed Brown extern PetscErrorCode RHSJacobian(TS,PetscReal,Vec,Mat,Mat,void*); 70c4762a1bSJed Brown extern PetscErrorCode Monitor(TS,PetscInt,PetscReal,Vec,void*); 71c4762a1bSJed Brown extern PetscErrorCode ExactSolution(PetscReal,Vec,AppCtx*); 72c4762a1bSJed Brown 73c4762a1bSJed Brown int main(int argc,char **argv) 74c4762a1bSJed Brown { 75c4762a1bSJed Brown AppCtx appctx; /* user-defined application context */ 76c4762a1bSJed Brown TS ts; /* timestepping context */ 77c4762a1bSJed Brown Mat A; /* Jacobian matrix data structure */ 78c4762a1bSJed Brown Vec u; /* approximate solution vector */ 79c4762a1bSJed Brown PetscInt time_steps_max = 100; /* default max timesteps */ 80c4762a1bSJed Brown PetscReal dt; 81c4762a1bSJed Brown PetscReal time_total_max = 100.0; /* default max total time */ 82c4762a1bSJed Brown PetscBool mymonitor = PETSC_FALSE; 83c4762a1bSJed Brown PetscReal bounds[] = {1.0, 3.3}; 84c4762a1bSJed Brown 85c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 86c4762a1bSJed Brown Initialize program and set problem parameters 87c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 88c4762a1bSJed Brown 89*b122ec5aSJacob Faibussowitsch CHKERRQ(PetscInitialize(&argc,&argv,(char*)0,help)); 905f80ce2aSJacob Faibussowitsch CHKERRQ(PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),1,bounds)); 91c4762a1bSJed Brown 92c4762a1bSJed Brown appctx.comm = PETSC_COMM_WORLD; 93c4762a1bSJed Brown appctx.m = 60; 94c4762a1bSJed Brown 955f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-M",&appctx.m,NULL)); 965f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsHasName(NULL,NULL,"-debug",&appctx.debug)); 975f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsHasName(NULL,NULL,"-mymonitor",&mymonitor)); 98c4762a1bSJed Brown 99c4762a1bSJed Brown appctx.h = 1.0/(appctx.m-1.0); 100c4762a1bSJed Brown 101c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 102c4762a1bSJed Brown Create vector data structures 103c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 104c4762a1bSJed Brown 105c4762a1bSJed Brown /* 106c4762a1bSJed Brown Create distributed array (DMDA) to manage parallel grid and vectors 107c4762a1bSJed Brown and to set up the ghost point communication pattern. There are M 108c4762a1bSJed Brown total grid values spread equally among all the processors. 109c4762a1bSJed Brown */ 1105f80ce2aSJacob Faibussowitsch CHKERRQ(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,appctx.m,1,1,NULL,&appctx.da)); 1115f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetFromOptions(appctx.da)); 1125f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetUp(appctx.da)); 113c4762a1bSJed Brown 114c4762a1bSJed Brown /* 115c4762a1bSJed Brown Extract global and local vectors from DMDA; we use these to store the 116c4762a1bSJed Brown approximate solution. Then duplicate these for remaining vectors that 117c4762a1bSJed Brown have the same types. 118c4762a1bSJed Brown */ 1195f80ce2aSJacob Faibussowitsch CHKERRQ(DMCreateGlobalVector(appctx.da,&u)); 1205f80ce2aSJacob Faibussowitsch CHKERRQ(DMCreateLocalVector(appctx.da,&appctx.u_local)); 121c4762a1bSJed Brown 122c4762a1bSJed Brown /* 123c4762a1bSJed Brown Create local work vector for use in evaluating right-hand-side function; 124c4762a1bSJed Brown create global work vector for storing exact solution. 125c4762a1bSJed Brown */ 1265f80ce2aSJacob Faibussowitsch CHKERRQ(VecDuplicate(appctx.u_local,&appctx.localwork)); 1275f80ce2aSJacob Faibussowitsch CHKERRQ(VecDuplicate(u,&appctx.solution)); 128c4762a1bSJed Brown 129c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 130c4762a1bSJed Brown Create timestepping solver context; set callback routine for 131c4762a1bSJed Brown right-hand-side function evaluation. 132c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 133c4762a1bSJed Brown 1345f80ce2aSJacob Faibussowitsch CHKERRQ(TSCreate(PETSC_COMM_WORLD,&ts)); 1355f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetProblemType(ts,TS_NONLINEAR)); 1365f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetRHSFunction(ts,NULL,RHSFunction,&appctx)); 137c4762a1bSJed Brown 138c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 139c4762a1bSJed Brown Set optional user-defined monitoring routine 140c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 141c4762a1bSJed Brown 142c4762a1bSJed Brown if (mymonitor) { 1435f80ce2aSJacob Faibussowitsch CHKERRQ(TSMonitorSet(ts,Monitor,&appctx,NULL)); 144c4762a1bSJed Brown } 145c4762a1bSJed Brown 146c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 147c4762a1bSJed Brown For nonlinear problems, the user can provide a Jacobian evaluation 148c4762a1bSJed Brown routine (or use a finite differencing approximation). 149c4762a1bSJed Brown 150c4762a1bSJed Brown Create matrix data structure; set Jacobian evaluation routine. 151c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 152c4762a1bSJed Brown 1535f80ce2aSJacob Faibussowitsch CHKERRQ(MatCreate(PETSC_COMM_WORLD,&A)); 1545f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,appctx.m,appctx.m)); 1555f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetFromOptions(A)); 1565f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetUp(A)); 1575f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetRHSJacobian(ts,A,A,RHSJacobian,&appctx)); 158c4762a1bSJed Brown 159c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 160c4762a1bSJed Brown Set solution vector and initial timestep 161c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 162c4762a1bSJed Brown 163c4762a1bSJed Brown dt = appctx.h/2.0; 1645f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetTimeStep(ts,dt)); 165c4762a1bSJed Brown 166c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 167c4762a1bSJed Brown Customize timestepping solver: 168c4762a1bSJed Brown - Set the solution method to be the Backward Euler method. 169c4762a1bSJed Brown - Set timestepping duration info 170c4762a1bSJed Brown Then set runtime options, which can override these defaults. 171c4762a1bSJed Brown For example, 172c4762a1bSJed Brown -ts_max_steps <maxsteps> -ts_max_time <maxtime> 173c4762a1bSJed Brown to override the defaults set by TSSetMaxSteps()/TSSetMaxTime(). 174c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 175c4762a1bSJed Brown 1765f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetType(ts,TSBEULER)); 1775f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetMaxSteps(ts,time_steps_max)); 1785f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetMaxTime(ts,time_total_max)); 1795f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 1805f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetFromOptions(ts)); 181c4762a1bSJed Brown 182c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 183c4762a1bSJed Brown Solve the problem 184c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 185c4762a1bSJed Brown 186c4762a1bSJed Brown /* 187c4762a1bSJed Brown Evaluate initial conditions 188c4762a1bSJed Brown */ 1895f80ce2aSJacob Faibussowitsch CHKERRQ(InitialConditions(u,&appctx)); 190c4762a1bSJed Brown 191c4762a1bSJed Brown /* 192c4762a1bSJed Brown Run the timestepping solver 193c4762a1bSJed Brown */ 1945f80ce2aSJacob Faibussowitsch CHKERRQ(TSSolve(ts,u)); 195c4762a1bSJed Brown 196c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 197c4762a1bSJed Brown Free work space. All PETSc objects should be destroyed when they 198c4762a1bSJed Brown are no longer needed. 199c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 200c4762a1bSJed Brown 2015f80ce2aSJacob Faibussowitsch CHKERRQ(TSDestroy(&ts)); 2025f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&u)); 2035f80ce2aSJacob Faibussowitsch CHKERRQ(MatDestroy(&A)); 2045f80ce2aSJacob Faibussowitsch CHKERRQ(DMDestroy(&appctx.da)); 2055f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&appctx.localwork)); 2065f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&appctx.solution)); 2075f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&appctx.u_local)); 208c4762a1bSJed Brown 209c4762a1bSJed Brown /* 210c4762a1bSJed Brown Always call PetscFinalize() before exiting a program. This routine 211c4762a1bSJed Brown - finalizes the PETSc libraries as well as MPI 212c4762a1bSJed Brown - provides summary and diagnostic information if certain runtime 213c4762a1bSJed Brown options are chosen (e.g., -log_view). 214c4762a1bSJed Brown */ 215*b122ec5aSJacob Faibussowitsch CHKERRQ(PetscFinalize()); 216*b122ec5aSJacob Faibussowitsch return 0; 217c4762a1bSJed Brown } 218c4762a1bSJed Brown /* --------------------------------------------------------------------- */ 219c4762a1bSJed Brown /* 220c4762a1bSJed Brown InitialConditions - Computes the solution at the initial time. 221c4762a1bSJed Brown 222c4762a1bSJed Brown Input Parameters: 223c4762a1bSJed Brown u - uninitialized solution vector (global) 224c4762a1bSJed Brown appctx - user-defined application context 225c4762a1bSJed Brown 226c4762a1bSJed Brown Output Parameter: 227c4762a1bSJed Brown u - vector with solution at initial time (global) 228c4762a1bSJed Brown */ 229c4762a1bSJed Brown PetscErrorCode InitialConditions(Vec u,AppCtx *appctx) 230c4762a1bSJed Brown { 231c4762a1bSJed Brown PetscScalar *u_localptr,h = appctx->h,x; 232c4762a1bSJed Brown PetscInt i,mybase,myend; 233c4762a1bSJed Brown 234c4762a1bSJed Brown /* 235c4762a1bSJed Brown Determine starting point of each processor's range of 236c4762a1bSJed Brown grid values. 237c4762a1bSJed Brown */ 2385f80ce2aSJacob Faibussowitsch CHKERRQ(VecGetOwnershipRange(u,&mybase,&myend)); 239c4762a1bSJed Brown 240c4762a1bSJed Brown /* 241c4762a1bSJed Brown Get a pointer to vector data. 242c4762a1bSJed Brown - For default PETSc vectors, VecGetArray() returns a pointer to 243c4762a1bSJed Brown the data array. Otherwise, the routine is implementation dependent. 244c4762a1bSJed Brown - You MUST call VecRestoreArray() when you no longer need access to 245c4762a1bSJed Brown the array. 246c4762a1bSJed Brown - Note that the Fortran interface to VecGetArray() differs from the 247c4762a1bSJed Brown C version. See the users manual for details. 248c4762a1bSJed Brown */ 2495f80ce2aSJacob Faibussowitsch CHKERRQ(VecGetArray(u,&u_localptr)); 250c4762a1bSJed Brown 251c4762a1bSJed Brown /* 252c4762a1bSJed Brown We initialize the solution array by simply writing the solution 253c4762a1bSJed Brown directly into the array locations. Alternatively, we could use 254c4762a1bSJed Brown VecSetValues() or VecSetValuesLocal(). 255c4762a1bSJed Brown */ 256c4762a1bSJed Brown for (i=mybase; i<myend; i++) { 257c4762a1bSJed Brown x = h*(PetscReal)i; /* current location in global grid */ 258c4762a1bSJed Brown u_localptr[i-mybase] = 1.0 + x*x; 259c4762a1bSJed Brown } 260c4762a1bSJed Brown 261c4762a1bSJed Brown /* 262c4762a1bSJed Brown Restore vector 263c4762a1bSJed Brown */ 2645f80ce2aSJacob Faibussowitsch CHKERRQ(VecRestoreArray(u,&u_localptr)); 265c4762a1bSJed Brown 266c4762a1bSJed Brown /* 267c4762a1bSJed Brown Print debugging information if desired 268c4762a1bSJed Brown */ 269c4762a1bSJed Brown if (appctx->debug) { 2705f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(appctx->comm,"initial guess vector\n")); 2715f80ce2aSJacob Faibussowitsch CHKERRQ(VecView(u,PETSC_VIEWER_STDOUT_WORLD)); 272c4762a1bSJed Brown } 273c4762a1bSJed Brown 274c4762a1bSJed Brown return 0; 275c4762a1bSJed Brown } 276c4762a1bSJed Brown /* --------------------------------------------------------------------- */ 277c4762a1bSJed Brown /* 278c4762a1bSJed Brown ExactSolution - Computes the exact solution at a given time. 279c4762a1bSJed Brown 280c4762a1bSJed Brown Input Parameters: 281c4762a1bSJed Brown t - current time 282c4762a1bSJed Brown solution - vector in which exact solution will be computed 283c4762a1bSJed Brown appctx - user-defined application context 284c4762a1bSJed Brown 285c4762a1bSJed Brown Output Parameter: 286c4762a1bSJed Brown solution - vector with the newly computed exact solution 287c4762a1bSJed Brown */ 288c4762a1bSJed Brown PetscErrorCode ExactSolution(PetscReal t,Vec solution,AppCtx *appctx) 289c4762a1bSJed Brown { 290c4762a1bSJed Brown PetscScalar *s_localptr,h = appctx->h,x; 291c4762a1bSJed Brown PetscInt i,mybase,myend; 292c4762a1bSJed Brown 293c4762a1bSJed Brown /* 294c4762a1bSJed Brown Determine starting and ending points of each processor's 295c4762a1bSJed Brown range of grid values 296c4762a1bSJed Brown */ 2975f80ce2aSJacob Faibussowitsch CHKERRQ(VecGetOwnershipRange(solution,&mybase,&myend)); 298c4762a1bSJed Brown 299c4762a1bSJed Brown /* 300c4762a1bSJed Brown Get a pointer to vector data. 301c4762a1bSJed Brown */ 3025f80ce2aSJacob Faibussowitsch CHKERRQ(VecGetArray(solution,&s_localptr)); 303c4762a1bSJed Brown 304c4762a1bSJed Brown /* 305c4762a1bSJed Brown Simply write the solution directly into the array locations. 306c4762a1bSJed Brown Alternatively, we could use VecSetValues() or VecSetValuesLocal(). 307c4762a1bSJed Brown */ 308c4762a1bSJed Brown for (i=mybase; i<myend; i++) { 309c4762a1bSJed Brown x = h*(PetscReal)i; 310c4762a1bSJed Brown s_localptr[i-mybase] = (t + 1.0)*(1.0 + x*x); 311c4762a1bSJed Brown } 312c4762a1bSJed Brown 313c4762a1bSJed Brown /* 314c4762a1bSJed Brown Restore vector 315c4762a1bSJed Brown */ 3165f80ce2aSJacob Faibussowitsch CHKERRQ(VecRestoreArray(solution,&s_localptr)); 317c4762a1bSJed Brown return 0; 318c4762a1bSJed Brown } 319c4762a1bSJed Brown /* --------------------------------------------------------------------- */ 320c4762a1bSJed Brown /* 321c4762a1bSJed Brown Monitor - User-provided routine to monitor the solution computed at 322c4762a1bSJed Brown each timestep. This example plots the solution and computes the 323c4762a1bSJed Brown error in two different norms. 324c4762a1bSJed Brown 325c4762a1bSJed Brown Input Parameters: 326c4762a1bSJed Brown ts - the timestep context 327c4762a1bSJed Brown step - the count of the current step (with 0 meaning the 328c4762a1bSJed Brown initial condition) 329c4762a1bSJed Brown time - the current time 330c4762a1bSJed Brown u - the solution at this timestep 331c4762a1bSJed Brown ctx - the user-provided context for this monitoring routine. 332c4762a1bSJed Brown In this case we use the application context which contains 333c4762a1bSJed Brown information about the problem size, workspace and the exact 334c4762a1bSJed Brown solution. 335c4762a1bSJed Brown */ 336c4762a1bSJed Brown PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec u,void *ctx) 337c4762a1bSJed Brown { 338c4762a1bSJed Brown AppCtx *appctx = (AppCtx*) ctx; /* user-defined application context */ 339c4762a1bSJed Brown PetscReal en2,en2s,enmax; 340c4762a1bSJed Brown PetscDraw draw; 341c4762a1bSJed Brown 342c4762a1bSJed Brown /* 343e1dfdf8eSBarry Smith We use the default X Windows viewer 344c4762a1bSJed Brown PETSC_VIEWER_DRAW_(appctx->comm) 345c4762a1bSJed Brown that is associated with the current communicator. This saves 346c4762a1bSJed Brown the effort of calling PetscViewerDrawOpen() to create the window. 347c4762a1bSJed Brown Note that if we wished to plot several items in separate windows we 348c4762a1bSJed Brown would create each viewer with PetscViewerDrawOpen() and store them in 349c4762a1bSJed Brown the application context, appctx. 350c4762a1bSJed Brown 351c4762a1bSJed Brown PetscReal buffering makes graphics look better. 352c4762a1bSJed Brown */ 3535f80ce2aSJacob Faibussowitsch CHKERRQ(PetscViewerDrawGetDraw(PETSC_VIEWER_DRAW_(appctx->comm),0,&draw)); 3545f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDrawSetDoubleBuffer(draw)); 3555f80ce2aSJacob Faibussowitsch CHKERRQ(VecView(u,PETSC_VIEWER_DRAW_(appctx->comm))); 356c4762a1bSJed Brown 357c4762a1bSJed Brown /* 358c4762a1bSJed Brown Compute the exact solution at this timestep 359c4762a1bSJed Brown */ 3605f80ce2aSJacob Faibussowitsch CHKERRQ(ExactSolution(time,appctx->solution,appctx)); 361c4762a1bSJed Brown 362c4762a1bSJed Brown /* 363c4762a1bSJed Brown Print debugging information if desired 364c4762a1bSJed Brown */ 365c4762a1bSJed Brown if (appctx->debug) { 3665f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(appctx->comm,"Computed solution vector\n")); 3675f80ce2aSJacob Faibussowitsch CHKERRQ(VecView(u,PETSC_VIEWER_STDOUT_WORLD)); 3685f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(appctx->comm,"Exact solution vector\n")); 3695f80ce2aSJacob Faibussowitsch CHKERRQ(VecView(appctx->solution,PETSC_VIEWER_STDOUT_WORLD)); 370c4762a1bSJed Brown } 371c4762a1bSJed Brown 372c4762a1bSJed Brown /* 373c4762a1bSJed Brown Compute the 2-norm and max-norm of the error 374c4762a1bSJed Brown */ 3755f80ce2aSJacob Faibussowitsch CHKERRQ(VecAXPY(appctx->solution,-1.0,u)); 3765f80ce2aSJacob Faibussowitsch CHKERRQ(VecNorm(appctx->solution,NORM_2,&en2)); 377c4762a1bSJed Brown en2s = PetscSqrtReal(appctx->h)*en2; /* scale the 2-norm by the grid spacing */ 3785f80ce2aSJacob Faibussowitsch CHKERRQ(VecNorm(appctx->solution,NORM_MAX,&enmax)); 379c4762a1bSJed Brown 380c4762a1bSJed Brown /* 381c4762a1bSJed Brown PetscPrintf() causes only the first processor in this 382c4762a1bSJed Brown communicator to print the timestep information. 383c4762a1bSJed Brown */ 3845f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(appctx->comm,"Timestep %D: time = %g 2-norm error = %g max norm error = %g\n",step,(double)time,(double)en2s,(double)enmax)); 385c4762a1bSJed Brown 386c4762a1bSJed Brown /* 387c4762a1bSJed Brown Print debugging information if desired 388c4762a1bSJed Brown */ 389c4762a1bSJed Brown if (appctx->debug) { 3905f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(appctx->comm,"Error vector\n")); 3915f80ce2aSJacob Faibussowitsch CHKERRQ(VecView(appctx->solution,PETSC_VIEWER_STDOUT_WORLD)); 392c4762a1bSJed Brown } 393c4762a1bSJed Brown return 0; 394c4762a1bSJed Brown } 395c4762a1bSJed Brown /* --------------------------------------------------------------------- */ 396c4762a1bSJed Brown /* 397c4762a1bSJed Brown RHSFunction - User-provided routine that evalues the right-hand-side 398c4762a1bSJed Brown function of the ODE. This routine is set in the main program by 399c4762a1bSJed Brown calling TSSetRHSFunction(). We compute: 400c4762a1bSJed Brown global_out = F(global_in) 401c4762a1bSJed Brown 402c4762a1bSJed Brown Input Parameters: 403c4762a1bSJed Brown ts - timesteping context 404c4762a1bSJed Brown t - current time 405c4762a1bSJed Brown global_in - vector containing the current iterate 406c4762a1bSJed Brown ctx - (optional) user-provided context for function evaluation. 407c4762a1bSJed Brown In this case we use the appctx defined above. 408c4762a1bSJed Brown 409c4762a1bSJed Brown Output Parameter: 410c4762a1bSJed Brown global_out - vector containing the newly evaluated function 411c4762a1bSJed Brown */ 412c4762a1bSJed Brown PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec global_in,Vec global_out,void *ctx) 413c4762a1bSJed Brown { 414c4762a1bSJed Brown AppCtx *appctx = (AppCtx*) ctx; /* user-defined application context */ 415c4762a1bSJed Brown DM da = appctx->da; /* distributed array */ 416c4762a1bSJed Brown Vec local_in = appctx->u_local; /* local ghosted input vector */ 417c4762a1bSJed Brown Vec localwork = appctx->localwork; /* local ghosted work vector */ 418c4762a1bSJed Brown PetscInt i,localsize; 419c4762a1bSJed Brown PetscMPIInt rank,size; 420c4762a1bSJed Brown PetscScalar *copyptr,sc; 421c4762a1bSJed Brown const PetscScalar *localptr; 422c4762a1bSJed Brown 423c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 424c4762a1bSJed Brown Get ready for local function computations 425c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 426c4762a1bSJed Brown /* 427c4762a1bSJed Brown Scatter ghost points to local vector, using the 2-step process 428c4762a1bSJed Brown DMGlobalToLocalBegin(), DMGlobalToLocalEnd(). 429c4762a1bSJed Brown By placing code between these two statements, computations can be 430c4762a1bSJed Brown done while messages are in transition. 431c4762a1bSJed Brown */ 4325f80ce2aSJacob Faibussowitsch CHKERRQ(DMGlobalToLocalBegin(da,global_in,INSERT_VALUES,local_in)); 4335f80ce2aSJacob Faibussowitsch CHKERRQ(DMGlobalToLocalEnd(da,global_in,INSERT_VALUES,local_in)); 434c4762a1bSJed Brown 435c4762a1bSJed Brown /* 436c4762a1bSJed Brown Access directly the values in our local INPUT work array 437c4762a1bSJed Brown */ 4385f80ce2aSJacob Faibussowitsch CHKERRQ(VecGetArrayRead(local_in,&localptr)); 439c4762a1bSJed Brown 440c4762a1bSJed Brown /* 441c4762a1bSJed Brown Access directly the values in our local OUTPUT work array 442c4762a1bSJed Brown */ 4435f80ce2aSJacob Faibussowitsch CHKERRQ(VecGetArray(localwork,©ptr)); 444c4762a1bSJed Brown 445c4762a1bSJed Brown sc = 1.0/(appctx->h*appctx->h*2.0*(1.0+t)*(1.0+t)); 446c4762a1bSJed Brown 447c4762a1bSJed Brown /* 448c4762a1bSJed Brown Evaluate our function on the nodes owned by this processor 449c4762a1bSJed Brown */ 4505f80ce2aSJacob Faibussowitsch CHKERRQ(VecGetLocalSize(local_in,&localsize)); 451c4762a1bSJed Brown 452c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 453c4762a1bSJed Brown Compute entries for the locally owned part 454c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 455c4762a1bSJed Brown 456c4762a1bSJed Brown /* 457c4762a1bSJed Brown Handle boundary conditions: This is done by using the boundary condition 458c4762a1bSJed Brown u(t,boundary) = g(t,boundary) 459c4762a1bSJed Brown for some function g. Now take the derivative with respect to t to obtain 460c4762a1bSJed Brown u_{t}(t,boundary) = g_{t}(t,boundary) 461c4762a1bSJed Brown 462c4762a1bSJed Brown In our case, u(t,0) = t + 1, so that u_{t}(t,0) = 1 463c4762a1bSJed Brown and u(t,1) = 2t+ 2, so that u_{t}(t,1) = 2 464c4762a1bSJed Brown */ 4655f80ce2aSJacob Faibussowitsch CHKERRMPI(MPI_Comm_rank(appctx->comm,&rank)); 4665f80ce2aSJacob Faibussowitsch CHKERRMPI(MPI_Comm_size(appctx->comm,&size)); 467dd400576SPatrick Sanan if (rank == 0) copyptr[0] = 1.0; 468c4762a1bSJed Brown if (rank == size-1) copyptr[localsize-1] = 2.0; 469c4762a1bSJed Brown 470c4762a1bSJed Brown /* 471c4762a1bSJed Brown Handle the interior nodes where the PDE is replace by finite 472c4762a1bSJed Brown difference operators. 473c4762a1bSJed Brown */ 474c4762a1bSJed Brown for (i=1; i<localsize-1; i++) copyptr[i] = localptr[i] * sc * (localptr[i+1] + localptr[i-1] - 2.0*localptr[i]); 475c4762a1bSJed Brown 476c4762a1bSJed Brown /* 477c4762a1bSJed Brown Restore vectors 478c4762a1bSJed Brown */ 4795f80ce2aSJacob Faibussowitsch CHKERRQ(VecRestoreArrayRead(local_in,&localptr)); 4805f80ce2aSJacob Faibussowitsch CHKERRQ(VecRestoreArray(localwork,©ptr)); 481c4762a1bSJed Brown 482c4762a1bSJed Brown /* 483c4762a1bSJed Brown Insert values from the local OUTPUT vector into the global 484c4762a1bSJed Brown output vector 485c4762a1bSJed Brown */ 4865f80ce2aSJacob Faibussowitsch CHKERRQ(DMLocalToGlobalBegin(da,localwork,INSERT_VALUES,global_out)); 4875f80ce2aSJacob Faibussowitsch CHKERRQ(DMLocalToGlobalEnd(da,localwork,INSERT_VALUES,global_out)); 488c4762a1bSJed Brown 489c4762a1bSJed Brown /* Print debugging information if desired */ 490c4762a1bSJed Brown if (appctx->debug) { 4915f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(appctx->comm,"RHS function vector\n")); 4925f80ce2aSJacob Faibussowitsch CHKERRQ(VecView(global_out,PETSC_VIEWER_STDOUT_WORLD)); 493c4762a1bSJed Brown } 494c4762a1bSJed Brown 495c4762a1bSJed Brown return 0; 496c4762a1bSJed Brown } 497c4762a1bSJed Brown /* --------------------------------------------------------------------- */ 498c4762a1bSJed Brown /* 499c4762a1bSJed Brown RHSJacobian - User-provided routine to compute the Jacobian of 500c4762a1bSJed Brown the nonlinear right-hand-side function of the ODE. 501c4762a1bSJed Brown 502c4762a1bSJed Brown Input Parameters: 503c4762a1bSJed Brown ts - the TS context 504c4762a1bSJed Brown t - current time 505c4762a1bSJed Brown global_in - global input vector 506c4762a1bSJed Brown dummy - optional user-defined context, as set by TSetRHSJacobian() 507c4762a1bSJed Brown 508c4762a1bSJed Brown Output Parameters: 509c4762a1bSJed Brown AA - Jacobian matrix 510c4762a1bSJed Brown BB - optionally different preconditioning matrix 511c4762a1bSJed Brown str - flag indicating matrix structure 512c4762a1bSJed Brown 513c4762a1bSJed Brown Notes: 514c4762a1bSJed Brown RHSJacobian computes entries for the locally owned part of the Jacobian. 515c4762a1bSJed Brown - Currently, all PETSc parallel matrix formats are partitioned by 516c4762a1bSJed Brown contiguous chunks of rows across the processors. 517c4762a1bSJed Brown - Each processor needs to insert only elements that it owns 518c4762a1bSJed Brown locally (but any non-local elements will be sent to the 519c4762a1bSJed Brown appropriate processor during matrix assembly). 520c4762a1bSJed Brown - Always specify global row and columns of matrix entries when 521c4762a1bSJed Brown using MatSetValues(). 522c4762a1bSJed Brown - Here, we set all entries for a particular row at once. 523c4762a1bSJed Brown - Note that MatSetValues() uses 0-based row and column numbers 524c4762a1bSJed Brown in Fortran as well as in C. 525c4762a1bSJed Brown */ 526c4762a1bSJed Brown PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec global_in,Mat AA,Mat BB,void *ctx) 527c4762a1bSJed Brown { 528c4762a1bSJed Brown AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */ 529c4762a1bSJed Brown Vec local_in = appctx->u_local; /* local ghosted input vector */ 530c4762a1bSJed Brown DM da = appctx->da; /* distributed array */ 531c4762a1bSJed Brown PetscScalar v[3],sc; 532c4762a1bSJed Brown const PetscScalar *localptr; 533c4762a1bSJed Brown PetscInt i,mstart,mend,mstarts,mends,idx[3],is; 534c4762a1bSJed Brown 535c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 536c4762a1bSJed Brown Get ready for local Jacobian computations 537c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 538c4762a1bSJed Brown /* 539c4762a1bSJed Brown Scatter ghost points to local vector, using the 2-step process 540c4762a1bSJed Brown DMGlobalToLocalBegin(), DMGlobalToLocalEnd(). 541c4762a1bSJed Brown By placing code between these two statements, computations can be 542c4762a1bSJed Brown done while messages are in transition. 543c4762a1bSJed Brown */ 5445f80ce2aSJacob Faibussowitsch CHKERRQ(DMGlobalToLocalBegin(da,global_in,INSERT_VALUES,local_in)); 5455f80ce2aSJacob Faibussowitsch CHKERRQ(DMGlobalToLocalEnd(da,global_in,INSERT_VALUES,local_in)); 546c4762a1bSJed Brown 547c4762a1bSJed Brown /* 548c4762a1bSJed Brown Get pointer to vector data 549c4762a1bSJed Brown */ 5505f80ce2aSJacob Faibussowitsch CHKERRQ(VecGetArrayRead(local_in,&localptr)); 551c4762a1bSJed Brown 552c4762a1bSJed Brown /* 553c4762a1bSJed Brown Get starting and ending locally owned rows of the matrix 554c4762a1bSJed Brown */ 5555f80ce2aSJacob Faibussowitsch CHKERRQ(MatGetOwnershipRange(BB,&mstarts,&mends)); 556c4762a1bSJed Brown mstart = mstarts; mend = mends; 557c4762a1bSJed Brown 558c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 559c4762a1bSJed Brown Compute entries for the locally owned part of the Jacobian. 560c4762a1bSJed Brown - Currently, all PETSc parallel matrix formats are partitioned by 561c4762a1bSJed Brown contiguous chunks of rows across the processors. 562c4762a1bSJed Brown - Each processor needs to insert only elements that it owns 563c4762a1bSJed Brown locally (but any non-local elements will be sent to the 564c4762a1bSJed Brown appropriate processor during matrix assembly). 565c4762a1bSJed Brown - Here, we set all entries for a particular row at once. 566c4762a1bSJed Brown - We can set matrix entries either using either 567c4762a1bSJed Brown MatSetValuesLocal() or MatSetValues(). 568c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 569c4762a1bSJed Brown 570c4762a1bSJed Brown /* 571c4762a1bSJed Brown Set matrix rows corresponding to boundary data 572c4762a1bSJed Brown */ 573c4762a1bSJed Brown if (mstart == 0) { 574c4762a1bSJed Brown v[0] = 0.0; 5755f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetValues(BB,1,&mstart,1,&mstart,v,INSERT_VALUES)); 576c4762a1bSJed Brown mstart++; 577c4762a1bSJed Brown } 578c4762a1bSJed Brown if (mend == appctx->m) { 579c4762a1bSJed Brown mend--; 580c4762a1bSJed Brown v[0] = 0.0; 5815f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetValues(BB,1,&mend,1,&mend,v,INSERT_VALUES)); 582c4762a1bSJed Brown } 583c4762a1bSJed Brown 584c4762a1bSJed Brown /* 585c4762a1bSJed Brown Set matrix rows corresponding to interior data. We construct the 586c4762a1bSJed Brown matrix one row at a time. 587c4762a1bSJed Brown */ 588c4762a1bSJed Brown sc = 1.0/(appctx->h*appctx->h*2.0*(1.0+t)*(1.0+t)); 589c4762a1bSJed Brown for (i=mstart; i<mend; i++) { 590c4762a1bSJed Brown idx[0] = i-1; idx[1] = i; idx[2] = i+1; 591c4762a1bSJed Brown is = i - mstart + 1; 592c4762a1bSJed Brown v[0] = sc*localptr[is]; 593c4762a1bSJed Brown v[1] = sc*(localptr[is+1] + localptr[is-1] - 4.0*localptr[is]); 594c4762a1bSJed Brown v[2] = sc*localptr[is]; 5955f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetValues(BB,1,&i,3,idx,v,INSERT_VALUES)); 596c4762a1bSJed Brown } 597c4762a1bSJed Brown 598c4762a1bSJed Brown /* 599c4762a1bSJed Brown Restore vector 600c4762a1bSJed Brown */ 6015f80ce2aSJacob Faibussowitsch CHKERRQ(VecRestoreArrayRead(local_in,&localptr)); 602c4762a1bSJed Brown 603c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 604c4762a1bSJed Brown Complete the matrix assembly process and set some options 605c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 606c4762a1bSJed Brown /* 607c4762a1bSJed Brown Assemble matrix, using the 2-step process: 608c4762a1bSJed Brown MatAssemblyBegin(), MatAssemblyEnd() 609c4762a1bSJed Brown Computations can be done while messages are in transition 610c4762a1bSJed Brown by placing code between these two statements. 611c4762a1bSJed Brown */ 6125f80ce2aSJacob Faibussowitsch CHKERRQ(MatAssemblyBegin(BB,MAT_FINAL_ASSEMBLY)); 6135f80ce2aSJacob Faibussowitsch CHKERRQ(MatAssemblyEnd(BB,MAT_FINAL_ASSEMBLY)); 614c4762a1bSJed Brown if (BB != AA) { 6155f80ce2aSJacob Faibussowitsch CHKERRQ(MatAssemblyBegin(AA,MAT_FINAL_ASSEMBLY)); 6165f80ce2aSJacob Faibussowitsch CHKERRQ(MatAssemblyEnd(AA,MAT_FINAL_ASSEMBLY)); 617c4762a1bSJed Brown } 618c4762a1bSJed Brown 619c4762a1bSJed Brown /* 620c4762a1bSJed Brown Set and option to indicate that we will never add a new nonzero location 621c4762a1bSJed Brown to the matrix. If we do, it will generate an error. 622c4762a1bSJed Brown */ 6235f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetOption(BB,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE)); 624c4762a1bSJed Brown 625c4762a1bSJed Brown return 0; 626c4762a1bSJed Brown } 627c4762a1bSJed Brown 628c4762a1bSJed Brown /*TEST 629c4762a1bSJed Brown 630c4762a1bSJed Brown test: 631c4762a1bSJed Brown args: -nox -ts_dt 10 -mymonitor 632c4762a1bSJed Brown nsize: 2 633c4762a1bSJed Brown requires: !single 634c4762a1bSJed Brown 635c4762a1bSJed Brown test: 636c4762a1bSJed Brown suffix: tut_1 637c4762a1bSJed Brown nsize: 1 638c4762a1bSJed Brown args: -ts_max_steps 10 -ts_monitor 639c4762a1bSJed Brown 640c4762a1bSJed Brown test: 641c4762a1bSJed Brown suffix: tut_2 642c4762a1bSJed Brown nsize: 4 643c4762a1bSJed Brown args: -ts_max_steps 10 -ts_monitor -snes_monitor -ksp_monitor 644c4762a1bSJed Brown 645c4762a1bSJed Brown test: 646c4762a1bSJed Brown suffix: tut_3 647c4762a1bSJed Brown nsize: 4 648c4762a1bSJed Brown args: ./ex2 -ts_max_steps 10 -ts_monitor -M 128 649c4762a1bSJed Brown 650c4762a1bSJed Brown TEST*/ 651