xref: /petsc/src/ts/tutorials/ex2.c (revision b122ec5aa1bd4469eb4e0673542fb7de3f411254)
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,&copyptr));
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,&copyptr));
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