static char help[] ="Solves a simple time-dependent linear PDE (the heat equation).\n\ Input parameters include:\n\ -m , where = number of grid points\n\ -time_dependent_rhs : Treat the problem as having a time-dependent right-hand side\n\ -debug : Activate debugging printouts\n\ -nox : Deactivate x-window graphics\n\n"; /* Concepts: TS^time-dependent linear problems Concepts: TS^heat equation Concepts: TS^diffusion equation Processors: n */ /* ------------------------------------------------------------------------ This program solves the one-dimensional heat equation (also called the diffusion equation), u_t = u_xx, on the domain 0 <= x <= 1, with the boundary conditions u(t,0) = 0, u(t,1) = 0, and the initial condition u(0,x) = sin(6*pi*x) + 3*sin(2*pi*x). This is a linear, second-order, parabolic equation. We discretize the right-hand side using finite differences with uniform grid spacing h: u_xx = (u_{i+1} - 2u_{i} + u_{i-1})/(h^2) We then demonstrate time evolution using the various TS methods by running the program via mpiexec -n ex3 -ts_type We compare the approximate solution with the exact solution, given by u_exact(x,t) = exp(-36*pi*pi*t) * sin(6*pi*x) + 3*exp(-4*pi*pi*t) * sin(2*pi*x) Notes: This code demonstrates the TS solver interface to two variants of linear problems, u_t = f(u,t), namely - time-dependent f: f(u,t) is a function of t - time-independent f: f(u,t) is simply f(u) The uniprocessor version of this code is ts/tutorials/ex3.c ------------------------------------------------------------------------- */ /* Include "petscdmda.h" so that we can use distributed arrays (DMDAs) to manage the parallel grid. Include "petscts.h" so that we can use TS solvers. Note that this file automatically includes: petscsys.h - base PETSc routines petscvec.h - vectors petscmat.h - matrices petscis.h - index sets petscksp.h - Krylov subspace methods petscviewer.h - viewers petscpc.h - preconditioners petscksp.h - linear solvers petscsnes.h - nonlinear solvers */ #include #include #include #include /* User-defined application context - contains data needed by the application-provided call-back routines. */ typedef struct { MPI_Comm comm; /* communicator */ DM da; /* distributed array data structure */ Vec localwork; /* local ghosted work vector */ Vec u_local; /* local ghosted approximate solution vector */ Vec solution; /* global exact solution vector */ PetscInt m; /* total number of grid points */ PetscReal h; /* mesh width h = 1/(m-1) */ PetscBool debug; /* flag (1 indicates activation of debugging printouts) */ PetscViewer viewer1,viewer2; /* viewers for the solution and error */ PetscReal norm_2,norm_max; /* error norms */ } AppCtx; /* User-defined routines */ extern PetscErrorCode InitialConditions(Vec,AppCtx*); extern PetscErrorCode RHSMatrixHeat(TS,PetscReal,Vec,Mat,Mat,void*); extern PetscErrorCode RHSFunctionHeat(TS,PetscReal,Vec,Vec,void*); extern PetscErrorCode Monitor(TS,PetscInt,PetscReal,Vec,void*); extern PetscErrorCode ExactSolution(PetscReal,Vec,AppCtx*); int main(int argc,char **argv) { AppCtx appctx; /* user-defined application context */ TS ts; /* timestepping context */ Mat A; /* matrix data structure */ Vec u; /* approximate solution vector */ PetscReal time_total_max = 1.0; /* default max total time */ PetscInt time_steps_max = 100; /* default max timesteps */ PetscDraw draw; /* drawing context */ PetscErrorCode ierr; PetscInt steps,m; PetscMPIInt size; PetscReal dt,ftime; PetscBool flg; TSProblemType tsproblem = TS_LINEAR; /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Initialize program and set problem parameters - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; appctx.comm = PETSC_COMM_WORLD; m = 60; ierr = PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);CHKERRQ(ierr); ierr = PetscOptionsHasName(NULL,NULL,"-debug",&appctx.debug);CHKERRQ(ierr); appctx.m = m; appctx.h = 1.0/(m-1.0); appctx.norm_2 = 0.0; appctx.norm_max = 0.0; ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"Solving a linear TS problem, number of processors = %d\n",size);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Create vector data structures - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* Create distributed array (DMDA) to manage parallel grid and vectors and to set up the ghost point communication pattern. There are M total grid values spread equally among all the processors. */ ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,m,1,1,NULL,&appctx.da);CHKERRQ(ierr); ierr = DMSetFromOptions(appctx.da);CHKERRQ(ierr); ierr = DMSetUp(appctx.da);CHKERRQ(ierr); /* Extract global and local vectors from DMDA; we use these to store the approximate solution. Then duplicate these for remaining vectors that have the same types. */ ierr = DMCreateGlobalVector(appctx.da,&u);CHKERRQ(ierr); ierr = DMCreateLocalVector(appctx.da,&appctx.u_local);CHKERRQ(ierr); /* Create local work vector for use in evaluating right-hand-side function; create global work vector for storing exact solution. */ ierr = VecDuplicate(appctx.u_local,&appctx.localwork);CHKERRQ(ierr); ierr = VecDuplicate(u,&appctx.solution);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Set up displays to show graphs of the solution and error - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"",80,380,400,160,&appctx.viewer1);CHKERRQ(ierr); ierr = PetscViewerDrawGetDraw(appctx.viewer1,0,&draw);CHKERRQ(ierr); ierr = PetscDrawSetDoubleBuffer(draw);CHKERRQ(ierr); ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"",80,0,400,160,&appctx.viewer2);CHKERRQ(ierr); ierr = PetscViewerDrawGetDraw(appctx.viewer2,0,&draw);CHKERRQ(ierr); ierr = PetscDrawSetDoubleBuffer(draw);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Create timestepping solver context - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); flg = PETSC_FALSE; ierr = PetscOptionsGetBool(NULL,NULL,"-nonlinear",&flg,NULL);CHKERRQ(ierr); ierr = TSSetProblemType(ts,flg ? TS_NONLINEAR : TS_LINEAR);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Set optional user-defined monitoring routine - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = TSMonitorSet(ts,Monitor,&appctx,NULL);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Create matrix data structure; set matrix evaluation routine. - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,m);CHKERRQ(ierr); ierr = MatSetFromOptions(A);CHKERRQ(ierr); ierr = MatSetUp(A);CHKERRQ(ierr); flg = PETSC_FALSE; ierr = PetscOptionsGetBool(NULL,NULL,"-time_dependent_rhs",&flg,NULL);CHKERRQ(ierr); if (flg) { /* For linear problems with a time-dependent f(u,t) in the equation u_t = f(u,t), the user provides the discretized right-hand-side as a time-dependent matrix. */ ierr = TSSetRHSFunction(ts,NULL,TSComputeRHSFunctionLinear,&appctx);CHKERRQ(ierr); ierr = TSSetRHSJacobian(ts,A,A,RHSMatrixHeat,&appctx);CHKERRQ(ierr); } else { /* For linear problems with a time-independent f(u) in the equation u_t = f(u), the user provides the discretized right-hand-side as a matrix only once, and then sets a null matrix evaluation routine. */ ierr = RHSMatrixHeat(ts,0.0,u,A,A,&appctx);CHKERRQ(ierr); ierr = TSSetRHSFunction(ts,NULL,TSComputeRHSFunctionLinear,&appctx);CHKERRQ(ierr); ierr = TSSetRHSJacobian(ts,A,A,TSComputeRHSJacobianConstant,&appctx);CHKERRQ(ierr); } if (tsproblem == TS_NONLINEAR) { SNES snes; ierr = TSSetRHSFunction(ts,NULL,RHSFunctionHeat,&appctx);CHKERRQ(ierr); ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); ierr = SNESSetJacobian(snes,NULL,NULL,SNESComputeJacobianDefault,NULL);CHKERRQ(ierr); } /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Set solution vector and initial timestep - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ dt = appctx.h*appctx.h/2.0; ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); ierr = TSSetSolution(ts,u);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Customize timestepping solver: - Set the solution method to be the Backward Euler method. - Set timestepping duration info Then set runtime options, which can override these defaults. For example, -ts_max_steps -ts_max_time to override the defaults set by TSSetMaxSteps()/TSSetMaxTime(). - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = TSSetMaxSteps(ts,time_steps_max);CHKERRQ(ierr); ierr = TSSetMaxTime(ts,time_total_max);CHKERRQ(ierr); ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); ierr = TSSetFromOptions(ts);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Solve the problem - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* Evaluate initial conditions */ ierr = InitialConditions(u,&appctx);CHKERRQ(ierr); /* Run the timestepping solver */ ierr = TSSolve(ts,u);CHKERRQ(ierr); ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr); ierr = TSGetStepNumber(ts,&steps);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - View timestepping solver info - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = PetscPrintf(PETSC_COMM_WORLD,"Total timesteps %D, Final time %g\n",steps,(double)ftime);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"Avg. error (2 norm) = %g Avg. error (max norm) = %g\n",(double)(appctx.norm_2/steps),(double)(appctx.norm_max/steps));CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Free work space. All PETSc objects should be destroyed when they are no longer needed. - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = TSDestroy(&ts);CHKERRQ(ierr); ierr = MatDestroy(&A);CHKERRQ(ierr); ierr = VecDestroy(&u);CHKERRQ(ierr); ierr = PetscViewerDestroy(&appctx.viewer1);CHKERRQ(ierr); ierr = PetscViewerDestroy(&appctx.viewer2);CHKERRQ(ierr); ierr = VecDestroy(&appctx.localwork);CHKERRQ(ierr); ierr = VecDestroy(&appctx.solution);CHKERRQ(ierr); ierr = VecDestroy(&appctx.u_local);CHKERRQ(ierr); ierr = DMDestroy(&appctx.da);CHKERRQ(ierr); /* Always call PetscFinalize() before exiting a program. This routine - finalizes the PETSc libraries as well as MPI - provides summary and diagnostic information if certain runtime options are chosen (e.g., -log_view). */ ierr = PetscFinalize(); return ierr; } /* --------------------------------------------------------------------- */ /* InitialConditions - Computes the solution at the initial time. Input Parameter: u - uninitialized solution vector (global) appctx - user-defined application context Output Parameter: u - vector with solution at initial time (global) */ PetscErrorCode InitialConditions(Vec u,AppCtx *appctx) { PetscScalar *u_localptr,h = appctx->h; PetscInt i,mybase,myend; PetscErrorCode ierr; /* Determine starting point of each processor's range of grid values. */ ierr = VecGetOwnershipRange(u,&mybase,&myend);CHKERRQ(ierr); /* Get a pointer to vector data. - For default PETSc vectors, VecGetArray() returns a pointer to the data array. Otherwise, the routine is implementation dependent. - You MUST call VecRestoreArray() when you no longer need access to the array. - Note that the Fortran interface to VecGetArray() differs from the C version. See the users manual for details. */ ierr = VecGetArray(u,&u_localptr);CHKERRQ(ierr); /* We initialize the solution array by simply writing the solution directly into the array locations. Alternatively, we could use VecSetValues() or VecSetValuesLocal(). */ for (i=mybase; idebug) { ierr = PetscPrintf(appctx->comm,"initial guess vector\n");CHKERRQ(ierr); ierr = VecView(u,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); } return 0; } /* --------------------------------------------------------------------- */ /* ExactSolution - Computes the exact solution at a given time. Input Parameters: t - current time solution - vector in which exact solution will be computed appctx - user-defined application context Output Parameter: solution - vector with the newly computed exact solution */ PetscErrorCode ExactSolution(PetscReal t,Vec solution,AppCtx *appctx) { PetscScalar *s_localptr,h = appctx->h,ex1,ex2,sc1,sc2; PetscInt i,mybase,myend; PetscErrorCode ierr; /* Determine starting and ending points of each processor's range of grid values */ ierr = VecGetOwnershipRange(solution,&mybase,&myend);CHKERRQ(ierr); /* Get a pointer to vector data. */ ierr = VecGetArray(solution,&s_localptr);CHKERRQ(ierr); /* Simply write the solution directly into the array locations. Alternatively, we culd use VecSetValues() or VecSetValuesLocal(). */ ex1 = PetscExpReal(-36.*PETSC_PI*PETSC_PI*t); ex2 = PetscExpReal(-4.*PETSC_PI*PETSC_PI*t); sc1 = PETSC_PI*6.*h; sc2 = PETSC_PI*2.*h; for (i=mybase; iviewer2);CHKERRQ(ierr); /* Compute the exact solution */ ierr = ExactSolution(time,appctx->solution,appctx);CHKERRQ(ierr); /* Print debugging information if desired */ if (appctx->debug) { ierr = PetscPrintf(appctx->comm,"Computed solution vector\n");CHKERRQ(ierr); ierr = VecView(u,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); ierr = PetscPrintf(appctx->comm,"Exact solution vector\n");CHKERRQ(ierr); ierr = VecView(appctx->solution,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); } /* Compute the 2-norm and max-norm of the error */ ierr = VecAXPY(appctx->solution,-1.0,u);CHKERRQ(ierr); ierr = VecNorm(appctx->solution,NORM_2,&norm_2);CHKERRQ(ierr); norm_2 = PetscSqrtReal(appctx->h)*norm_2; ierr = VecNorm(appctx->solution,NORM_MAX,&norm_max);CHKERRQ(ierr); if (norm_2 < 1e-14) norm_2 = 0; if (norm_max < 1e-14) norm_max = 0; /* PetscPrintf() causes only the first processor in this communicator to print the timestep information. */ ierr = PetscPrintf(appctx->comm,"Timestep %D: time = %g 2-norm error = %g max norm error = %g\n",step,(double)time,(double)norm_2,(double)norm_max);CHKERRQ(ierr); appctx->norm_2 += norm_2; appctx->norm_max += norm_max; /* View a graph of the error */ ierr = VecView(appctx->solution,appctx->viewer1);CHKERRQ(ierr); /* Print debugging information if desired */ if (appctx->debug) { ierr = PetscPrintf(appctx->comm,"Error vector\n");CHKERRQ(ierr); ierr = VecView(appctx->solution,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); } return 0; } /* --------------------------------------------------------------------- */ /* RHSMatrixHeat - User-provided routine to compute the right-hand-side matrix for the heat equation. Input Parameters: ts - the TS context t - current time global_in - global input vector dummy - optional user-defined context, as set by TSetRHSJacobian() Output Parameters: AA - Jacobian matrix BB - optionally different preconditioning matrix str - flag indicating matrix structure Notes: RHSMatrixHeat computes entries for the locally owned part of the system. - Currently, all PETSc parallel matrix formats are partitioned by contiguous chunks of rows across the processors. - Each processor needs to insert only elements that it owns locally (but any non-local elements will be sent to the appropriate processor during matrix assembly). - Always specify global row and columns of matrix entries when using MatSetValues(); we could alternatively use MatSetValuesLocal(). - Here, we set all entries for a particular row at once. - Note that MatSetValues() uses 0-based row and column numbers in Fortran as well as in C. */ PetscErrorCode RHSMatrixHeat(TS ts,PetscReal t,Vec X,Mat AA,Mat BB,void *ctx) { Mat A = AA; /* Jacobian matrix */ AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */ PetscErrorCode ierr; PetscInt i,mstart,mend,idx[3]; PetscScalar v[3],stwo = -2./(appctx->h*appctx->h),sone = -.5*stwo; /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Compute entries for the locally owned part of the matrix - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = MatGetOwnershipRange(A,&mstart,&mend);CHKERRQ(ierr); /* Set matrix rows corresponding to boundary data */ if (mstart == 0) { /* first processor only */ v[0] = 1.0; ierr = MatSetValues(A,1,&mstart,1,&mstart,v,INSERT_VALUES);CHKERRQ(ierr); mstart++; } if (mend == appctx->m) { /* last processor only */ mend--; v[0] = 1.0; ierr = MatSetValues(A,1,&mend,1,&mend,v,INSERT_VALUES);CHKERRQ(ierr); } /* Set matrix rows corresponding to interior data. We construct the matrix one row at a time. */ v[0] = sone; v[1] = stwo; v[2] = sone; for (i=mstart; i