1 static const char help[] = "Time-dependent PDE in 1d. Simplified from ex15.c for illustrating how to solve DAEs. \n";
2 /*
3 u_t = uxx
4 0 < x < 1;
5 At t=0: u(x) = exp(c*r*r*r), if r=PetscSqrtReal((x-.5)*(x-.5)) < .125
6 u(x) = 0.0 if r >= .125
7
8
9 Boundary conditions:
10 Dirichlet BC:
11 At x=0, x=1, u = 0.0
12
13 Neumann BC:
14 At x=0, x=1: du(x,t)/dx = 0
15
16 mpiexec -n 2 ./ex17 -da_grid_x 40 -ts_max_steps 2 -snes_monitor -ksp_monitor
17 ./ex17 -da_grid_x 40 -monitor_solution
18 ./ex17 -da_grid_x 100 -ts_type theta -ts_theta_theta 0.5 # Midpoint is not L-stable
19 ./ex17 -jac_type fd_coloring -da_grid_x 500 -boundary 1
20 ./ex17 -da_grid_x 100 -ts_type gl -ts_adapt_type none -ts_max_steps 2
21 */
22
23 #include <petscdm.h>
24 #include <petscdmda.h>
25 #include <petscts.h>
26
27 typedef enum {JACOBIAN_ANALYTIC,JACOBIAN_FD_COLORING,JACOBIAN_FD_FULL} JacobianType;
28 static const char *const JacobianTypes[] = {"analytic","fd_coloring","fd_full","JacobianType","fd_",0};
29
30 /*
31 User-defined data structures and routines
32 */
33 typedef struct {
34 PetscReal c;
35 PetscInt boundary; /* Type of boundary condition */
36 PetscBool viewJacobian;
37 } AppCtx;
38
39 static PetscErrorCode FormIFunction(TS,PetscReal,Vec,Vec,Vec,void*);
40 static PetscErrorCode FormIJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*);
41 static PetscErrorCode FormInitialSolution(TS,Vec,void*);
42
main(int argc,char ** argv)43 int main(int argc,char **argv)
44 {
45 TS ts; /* nonlinear solver */
46 Vec u; /* solution, residual vectors */
47 Mat J; /* Jacobian matrix */
48 PetscInt nsteps;
49 PetscReal vmin,vmax,norm;
50 PetscErrorCode ierr;
51 DM da;
52 PetscReal ftime,dt;
53 AppCtx user; /* user-defined work context */
54 JacobianType jacType;
55
56 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
57
58 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
59 Create distributed array (DMDA) to manage parallel grid and vectors
60 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
61 ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,11,1,1,NULL,&da);CHKERRQ(ierr);
62 ierr = DMSetFromOptions(da);CHKERRQ(ierr);
63 ierr = DMSetUp(da);CHKERRQ(ierr);
64
65 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
66 Extract global vectors from DMDA;
67 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
68 ierr = DMCreateGlobalVector(da,&u);CHKERRQ(ierr);
69
70 /* Initialize user application context */
71 user.c = -30.0;
72 user.boundary = 0; /* 0: Dirichlet BC; 1: Neumann BC */
73 user.viewJacobian = PETSC_FALSE;
74
75 ierr = PetscOptionsGetInt(NULL,NULL,"-boundary",&user.boundary,NULL);CHKERRQ(ierr);
76 ierr = PetscOptionsHasName(NULL,NULL,"-viewJacobian",&user.viewJacobian);CHKERRQ(ierr);
77
78 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
79 Create timestepping solver context
80 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
81 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
82 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr);
83 ierr = TSSetType(ts,TSTHETA);CHKERRQ(ierr);
84 ierr = TSThetaSetTheta(ts,1.0);CHKERRQ(ierr); /* Make the Theta method behave like backward Euler */
85 ierr = TSSetIFunction(ts,NULL,FormIFunction,&user);CHKERRQ(ierr);
86
87 ierr = DMSetMatType(da,MATAIJ);CHKERRQ(ierr);
88 ierr = DMCreateMatrix(da,&J);CHKERRQ(ierr);
89 jacType = JACOBIAN_ANALYTIC; /* use user-provide Jacobian */
90
91 ierr = TSSetDM(ts,da);CHKERRQ(ierr); /* Use TSGetDM() to access. Setting here allows easy use of geometric multigrid. */
92
93 ftime = 1.0;
94 ierr = TSSetMaxTime(ts,ftime);CHKERRQ(ierr);
95 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr);
96
97 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
98 Set initial conditions
99 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
100 ierr = FormInitialSolution(ts,u,&user);CHKERRQ(ierr);
101 ierr = TSSetSolution(ts,u);CHKERRQ(ierr);
102 dt = .01;
103 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr);
104
105
106 /* Use slow fd Jacobian or fast fd Jacobian with colorings.
107 Note: this requirs snes which is not created until TSSetUp()/TSSetFromOptions() is called */
108 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Options for Jacobian evaluation",NULL);CHKERRQ(ierr);
109 ierr = PetscOptionsEnum("-jac_type","Type of Jacobian","",JacobianTypes,(PetscEnum)jacType,(PetscEnum*)&jacType,0);CHKERRQ(ierr);
110 ierr = PetscOptionsEnd();CHKERRQ(ierr);
111 if (jacType == JACOBIAN_ANALYTIC) {
112 ierr = TSSetIJacobian(ts,J,J,FormIJacobian,&user);CHKERRQ(ierr);
113 } else if (jacType == JACOBIAN_FD_COLORING) {
114 SNES snes;
115 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr);
116 ierr = SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,0);CHKERRQ(ierr);
117 } else if (jacType == JACOBIAN_FD_FULL) {
118 SNES snes;
119 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr);
120 ierr = SNESSetJacobian(snes,J,J,SNESComputeJacobianDefault,&user);CHKERRQ(ierr);
121 }
122
123 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124 Set runtime options
125 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
126 ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
127
128 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129 Integrate ODE system
130 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
131 ierr = TSSolve(ts,u);CHKERRQ(ierr);
132
133 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134 Compute diagnostics of the solution
135 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
136 ierr = VecNorm(u,NORM_1,&norm);CHKERRQ(ierr);
137 ierr = VecMax(u,NULL,&vmax);CHKERRQ(ierr);
138 ierr = VecMin(u,NULL,&vmin);CHKERRQ(ierr);
139 ierr = TSGetStepNumber(ts,&nsteps);CHKERRQ(ierr);
140 ierr = TSGetTime(ts,&ftime);CHKERRQ(ierr);
141 ierr = PetscPrintf(PETSC_COMM_WORLD,"timestep %D: time %g, solution norm %g, max %g, min %g\n",nsteps,(double)ftime,(double)norm,(double)vmax,(double)vmin);CHKERRQ(ierr);
142
143 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144 Free work space.
145 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
146 ierr = MatDestroy(&J);CHKERRQ(ierr);
147 ierr = VecDestroy(&u);CHKERRQ(ierr);
148 ierr = TSDestroy(&ts);CHKERRQ(ierr);
149 ierr = DMDestroy(&da);CHKERRQ(ierr);
150 ierr = PetscFinalize();
151 return ierr;
152 }
153 /* ------------------------------------------------------------------- */
FormIFunction(TS ts,PetscReal ftime,Vec U,Vec Udot,Vec F,void * ptr)154 static PetscErrorCode FormIFunction(TS ts,PetscReal ftime,Vec U,Vec Udot,Vec F,void *ptr)
155 {
156 AppCtx *user=(AppCtx*)ptr;
157 DM da;
158 PetscErrorCode ierr;
159 PetscInt i,Mx,xs,xm;
160 PetscReal hx,sx;
161 PetscScalar *u,*udot,*f;
162 Vec localU;
163
164 PetscFunctionBeginUser;
165 ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
166 ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr);
167 ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr);
168
169 hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx);
170
171 /*
172 Scatter ghost points to local vector,using the 2-step process
173 DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
174 By placing code between these two statements, computations can be
175 done while messages are in transition.
176 */
177 ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);
178 ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);
179
180 /* Get pointers to vector data */
181 ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr);
182 ierr = DMDAVecGetArrayRead(da,Udot,&udot);CHKERRQ(ierr);
183 ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr);
184
185 /* Get local grid boundaries */
186 ierr = DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr);
187
188 /* Compute function over the locally owned part of the grid */
189 for (i=xs; i<xs+xm; i++) {
190 if (user->boundary == 0) { /* Dirichlet BC */
191 if (i == 0 || i == Mx-1) f[i] = u[i]; /* F = U */
192 else f[i] = udot[i] + (2.*u[i] - u[i-1] - u[i+1])*sx;
193 } else { /* Neumann BC */
194 if (i == 0) f[i] = u[0] - u[1];
195 else if (i == Mx-1) f[i] = u[i] - u[i-1];
196 else f[i] = udot[i] + (2.*u[i] - u[i-1] - u[i+1])*sx;
197 }
198 }
199
200 /* Restore vectors */
201 ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr);
202 ierr = DMDAVecRestoreArrayRead(da,Udot,&udot);CHKERRQ(ierr);
203 ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr);
204 ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr);
205 PetscFunctionReturn(0);
206 }
207
208 /* --------------------------------------------------------------------- */
209 /*
210 IJacobian - Compute IJacobian = dF/dU + a dF/dUdot
211 */
FormIJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat J,Mat Jpre,void * ctx)212 PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat J,Mat Jpre,void *ctx)
213 {
214 PetscErrorCode ierr;
215 PetscInt i,rstart,rend,Mx;
216 PetscReal hx,sx;
217 AppCtx *user = (AppCtx*)ctx;
218 DM da;
219 MatStencil col[3],row;
220 PetscInt nc;
221 PetscScalar vals[3];
222
223 PetscFunctionBeginUser;
224 ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
225 ierr = MatGetOwnershipRange(Jpre,&rstart,&rend);CHKERRQ(ierr);
226 ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr);
227 hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx);
228 for (i=rstart; i<rend; i++) {
229 nc = 0;
230 row.i = i;
231 if (user->boundary == 0 && (i == 0 || i == Mx-1)) {
232 col[nc].i = i; vals[nc++] = 1.0;
233 } else if (user->boundary > 0 && i == 0) { /* Left Neumann */
234 col[nc].i = i; vals[nc++] = 1.0;
235 col[nc].i = i+1; vals[nc++] = -1.0;
236 } else if (user->boundary > 0 && i == Mx-1) { /* Right Neumann */
237 col[nc].i = i-1; vals[nc++] = -1.0;
238 col[nc].i = i; vals[nc++] = 1.0;
239 } else { /* Interior */
240 col[nc].i = i-1; vals[nc++] = -1.0*sx;
241 col[nc].i = i; vals[nc++] = 2.0*sx + a;
242 col[nc].i = i+1; vals[nc++] = -1.0*sx;
243 }
244 ierr = MatSetValuesStencil(Jpre,1,&row,nc,col,vals,INSERT_VALUES);CHKERRQ(ierr);
245 }
246
247 ierr = MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
248 ierr = MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
249 if (J != Jpre) {
250 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
251 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
252 }
253 if (user->viewJacobian) {
254 ierr = PetscPrintf(PETSC_COMM_WORLD,"Jpre:\n");CHKERRQ(ierr);
255 ierr = MatView(Jpre,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
256 }
257 PetscFunctionReturn(0);
258 }
259
260 /* ------------------------------------------------------------------- */
FormInitialSolution(TS ts,Vec U,void * ptr)261 PetscErrorCode FormInitialSolution(TS ts,Vec U,void *ptr)
262 {
263 AppCtx *user=(AppCtx*)ptr;
264 PetscReal c =user->c;
265 DM da;
266 PetscErrorCode ierr;
267 PetscInt i,xs,xm,Mx;
268 PetscScalar *u;
269 PetscReal hx,x,r;
270
271 PetscFunctionBeginUser;
272 ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
273 ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr);
274
275 hx = 1.0/(PetscReal)(Mx-1);
276
277 /* Get pointers to vector data */
278 ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr);
279
280 /* Get local grid boundaries */
281 ierr = DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr);
282
283 /* Compute function over the locally owned part of the grid */
284 for (i=xs; i<xs+xm; i++) {
285 x = i*hx;
286 r = PetscSqrtReal((x-.5)*(x-.5));
287 if (r < .125) u[i] = PetscExpReal(c*r*r*r);
288 else u[i] = 0.0;
289 }
290
291 /* Restore vectors */
292 ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr);
293 PetscFunctionReturn(0);
294 }
295
296 /*TEST
297
298 test:
299 requires: !single
300 args: -da_grid_x 40 -ts_max_steps 2 -snes_monitor_short -ksp_monitor_short -ts_monitor
301
302 test:
303 suffix: 2
304 requires: !single
305 args: -da_grid_x 100 -ts_type theta -ts_theta_theta 0.5
306
307
308 TEST*/
309
310