Actual source code: ex2.c
petsc-3.5.4 2015-05-23
1: /*
2: Formatted test for TS routines.
4: Solves U_t=F(t,u)
5: Where:
7: [2*u1+u2
8: F(t,u)= [u1+2*u2+u3
9: [ u2+2*u3
10: We can compare the solutions from euler, beuler and SUNDIALS to
11: see what is the difference.
13: */
15: static char help[] = "Solves a nonlinear ODE. \n\n";
17: #include <petscts.h>
18: #include <petscpc.h>
20: extern PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*);
21: extern PetscErrorCode RHSJacobian(TS,PetscReal,Vec,Mat,Mat,void*);
22: extern PetscErrorCode Monitor(TS,PetscInt,PetscReal,Vec,void*);
23: extern PetscErrorCode Initial(Vec,void*);
25: extern PetscReal solx(PetscReal);
26: extern PetscReal soly(PetscReal);
27: extern PetscReal solz(PetscReal);
31: int main(int argc,char **argv)
32: {
34: PetscInt time_steps = 100,steps;
35: PetscMPIInt size;
36: Vec global;
37: PetscReal dt,ftime;
38: TS ts;
39: Mat A = 0;
41: PetscInitialize(&argc,&argv,(char*)0,help);
42: MPI_Comm_size(PETSC_COMM_WORLD,&size);
44: PetscOptionsGetInt(NULL,"-time",&time_steps,NULL);
46: /* set initial conditions */
47: VecCreate(PETSC_COMM_WORLD,&global);
48: VecSetSizes(global,PETSC_DECIDE,3);
49: VecSetFromOptions(global);
50: Initial(global,NULL);
52: /* make timestep context */
53: TSCreate(PETSC_COMM_WORLD,&ts);
54: TSSetProblemType(ts,TS_NONLINEAR);
55: TSMonitorSet(ts,Monitor,NULL,NULL);
57: dt = 0.1;
59: /*
60: The user provides the RHS and Jacobian
61: */
62: TSSetRHSFunction(ts,NULL,RHSFunction,NULL);
63: MatCreate(PETSC_COMM_WORLD,&A);
64: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,3,3);
65: MatSetFromOptions(A);
66: MatSetUp(A);
67: RHSJacobian(ts,0.0,global,A,A,NULL);
68: TSSetRHSJacobian(ts,A,A,RHSJacobian,NULL);
70: TSSetFromOptions(ts);
72: TSSetInitialTimeStep(ts,0.0,dt);
73: TSSetDuration(ts,time_steps,1);
74: TSSetSolution(ts,global);
76: TSSolve(ts,global);
77: TSGetSolveTime(ts,&ftime);
78: TSGetTimeStepNumber(ts,&steps);
81: /* free the memories */
83: TSDestroy(&ts);
84: VecDestroy(&global);
85: MatDestroy(&A);
87: PetscFinalize();
88: return 0;
89: }
91: /* -------------------------------------------------------------------*/
94: /* this test problem has initial values (1,1,1). */
95: PetscErrorCode Initial(Vec global,void *ctx)
96: {
97: PetscScalar *localptr;
98: PetscInt i,mybase,myend,locsize;
101: /* determine starting point of each processor */
102: VecGetOwnershipRange(global,&mybase,&myend);
103: VecGetLocalSize(global,&locsize);
105: /* Initialize the array */
106: VecGetArray(global,&localptr);
107: for (i=0; i<locsize; i++) localptr[i] = 1.0;
109: if (mybase == 0) localptr[0]=1.0;
111: VecRestoreArray(global,&localptr);
112: return 0;
113: }
117: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec global,void *ctx)
118: {
119: VecScatter scatter;
120: IS from,to;
121: PetscInt i,n,*idx;
122: Vec tmp_vec;
124: PetscScalar *tmp;
126: /* Get the size of the vector */
127: VecGetSize(global,&n);
129: /* Set the index sets */
130: PetscMalloc1(n,&idx);
131: for (i=0; i<n; i++) idx[i]=i;
133: /* Create local sequential vectors */
134: VecCreateSeq(PETSC_COMM_SELF,n,&tmp_vec);
136: /* Create scatter context */
137: ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from);
138: ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to);
139: VecScatterCreate(global,from,tmp_vec,to,&scatter);
140: VecScatterBegin(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
141: VecScatterEnd(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
143: VecGetArray(tmp_vec,&tmp);
144: PetscPrintf(PETSC_COMM_WORLD,"At t =%14.6e u = %14.6e %14.6e %14.6e \n",
145: time,PetscRealPart(tmp[0]),PetscRealPart(tmp[1]),PetscRealPart(tmp[2]));
146: PetscPrintf(PETSC_COMM_WORLD,"At t =%14.6e errors = %14.6e %14.6e %14.6e \n",
147: time,PetscRealPart(tmp[0]-solx(time)),PetscRealPart(tmp[1]-soly(time)),PetscRealPart(tmp[2]-solz(time)));
148: VecRestoreArray(tmp_vec,&tmp);
149: VecScatterDestroy(&scatter);
150: ISDestroy(&from);
151: ISDestroy(&to);
152: PetscFree(idx);
153: VecDestroy(&tmp_vec);
154: return 0;
155: }
159: PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
160: {
161: PetscScalar *inptr,*outptr;
162: PetscInt i,n,*idx;
164: IS from,to;
165: VecScatter scatter;
166: Vec tmp_in,tmp_out;
168: /* Get the length of parallel vector */
169: VecGetSize(globalin,&n);
171: /* Set the index sets */
172: PetscMalloc1(n,&idx);
173: for (i=0; i<n; i++) idx[i]=i;
175: /* Create local sequential vectors */
176: VecCreateSeq(PETSC_COMM_SELF,n,&tmp_in);
177: VecDuplicate(tmp_in,&tmp_out);
179: /* Create scatter context */
180: ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from);
181: ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to);
182: VecScatterCreate(globalin,from,tmp_in,to,&scatter);
183: VecScatterBegin(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
184: VecScatterEnd(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
185: VecScatterDestroy(&scatter);
187: /*Extract income array */
188: VecGetArray(tmp_in,&inptr);
190: /* Extract outcome array*/
191: VecGetArray(tmp_out,&outptr);
193: outptr[0] = 2.0*inptr[0]+inptr[1];
194: outptr[1] = inptr[0]+2.0*inptr[1]+inptr[2];
195: outptr[2] = inptr[1]+2.0*inptr[2];
197: VecRestoreArray(tmp_in,&inptr);
198: VecRestoreArray(tmp_out,&outptr);
200: VecScatterCreate(tmp_out,from,globalout,to,&scatter);
201: VecScatterBegin(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
202: VecScatterEnd(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
204: /* Destroy idx aand scatter */
205: ISDestroy(&from);
206: ISDestroy(&to);
207: VecScatterDestroy(&scatter);
208: VecDestroy(&tmp_in);
209: VecDestroy(&tmp_out);
210: PetscFree(idx);
211: return 0;
212: }
216: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec x,Mat A,Mat BB,void *ctx)
217: {
218: PetscScalar v[3],*tmp;
219: PetscInt idx[3],i;
222: idx[0]=0; idx[1]=1; idx[2]=2;
223: VecGetArray(x,&tmp);
225: i = 0;
226: v[0] = 2.0; v[1] = 1.0; v[2] = 0.0;
227: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
229: i = 1;
230: v[0] = 1.0; v[1] = 2.0; v[2] = 1.0;
231: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
233: i = 2;
234: v[0] = 0.0; v[1] = 1.0; v[2] = 2.0;
235: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
237: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
238: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
240: VecRestoreArray(x,&tmp);
242: return 0;
243: }
245: /*
246: The exact solutions
247: */
248: PetscReal solx(PetscReal t)
249: {
250: return PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/(2.0*PetscSqrtReal(2.0)) +
251: PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/(2.0*PetscSqrtReal(2.0));
252: }
254: PetscReal soly(PetscReal t)
255: {
256: return PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/PetscSqrtReal(2.0) +
257: PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/PetscSqrtReal(2.0);
258: }
260: PetscReal solz(PetscReal t)
261: {
262: return PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/(2.0*PetscSqrtReal(2.0)) +
263: PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/(2.0*PetscSqrtReal(2.0));
264: }