Actual source code: ex2.c

  1: /*
  2:        Formatted test for TS routines.

  4:           Solves U_t=F(t,u)
  5:           Where:
  6:           
  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



 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:   MatStructure   A_structure;
 40:   Mat            A = 0;
 41: 
 42:   PetscInitialize(&argc,&argv,(char*)0,help);
 43:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 44: 
 45:   PetscOptionsGetInt(PETSC_NULL,"-time",&time_steps,PETSC_NULL);
 46: 
 47:   /* set initial conditions */
 48:   VecCreate(PETSC_COMM_WORLD,&global);
 49:   VecSetSizes(global,PETSC_DECIDE,3);
 50:   VecSetFromOptions(global);
 51:   Initial(global,PETSC_NULL);
 52: 
 53:   /* make timestep context */
 54:   TSCreate(PETSC_COMM_WORLD,&ts);
 55:   TSSetProblemType(ts,TS_NONLINEAR);
 56:   TSMonitorSet(ts,Monitor,PETSC_NULL,PETSC_NULL);

 58:   dt = 0.1;

 60:   /*
 61:     The user provides the RHS and Jacobian
 62:   */
 63:   TSSetRHSFunction(ts,RHSFunction,NULL);
 64:   MatCreate(PETSC_COMM_WORLD,&A);
 65:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,3,3);
 66:   MatSetFromOptions(A);
 67:   RHSJacobian(ts,0.0,global,&A,&A,&A_structure,NULL);
 68:   TSSetRHSJacobian(ts,A,A,RHSJacobian,NULL);
 69: 
 70:   TSSetFromOptions(ts);

 72:   TSSetInitialTimeStep(ts,0.0,dt);
 73:   TSSetDuration(ts,time_steps,1);
 74:   TSSetSolution(ts,global);


 77:   TSSetUp(ts);
 78:   TSStep(ts,&steps,&ftime);


 81:   /* free the memories */
 82: 
 83:   TSDestroy(ts);
 84:   VecDestroy(global);
 85:   ierr= 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++) {
108:     localptr[i] = 1.0;
109:   }
110: 
111:   if (mybase == 0) localptr[0]=1.0;

113:   VecRestoreArray(global,&localptr);
114:   return 0;
115: }

119: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec global,void *ctx)
120: {
121:   VecScatter     scatter;
122:   IS             from,to;
123:   PetscInt       i,n,*idx;
124:   Vec            tmp_vec;
126:   PetscScalar    *tmp;

128:   /* Get the size of the vector */
129:   VecGetSize(global,&n);

131:   /* Set the index sets */
132:   PetscMalloc(n*sizeof(PetscInt),&idx);
133:   for(i=0; i<n; i++) idx[i]=i;
134: 
135:   /* Create local sequential vectors */
136:   VecCreateSeq(PETSC_COMM_SELF,n,&tmp_vec);

138:   /* Create scatter context */
139:   ISCreateGeneral(PETSC_COMM_SELF,n,idx,&from);
140:   ISCreateGeneral(PETSC_COMM_SELF,n,idx,&to);
141:   VecScatterCreate(global,from,tmp_vec,to,&scatter);
142:   VecScatterBegin(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
143:   VecScatterEnd(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);

145:   VecGetArray(tmp_vec,&tmp);
146:   PetscPrintf(PETSC_COMM_WORLD,"At t =%14.6e u = %14.6e  %14.6e  %14.6e \n",
147:                      time,PetscRealPart(tmp[0]),PetscRealPart(tmp[1]),PetscRealPart(tmp[2]));
148:   PetscPrintf(PETSC_COMM_WORLD,"At t =%14.6e errors = %14.6e  %14.6e  %14.6e \n",
149:                      time,PetscRealPart(tmp[0]-solx(time)),PetscRealPart(tmp[1]-soly(time)),PetscRealPart(tmp[2]-solz(time)));
150:   VecRestoreArray(tmp_vec,&tmp);
151:   VecScatterDestroy(scatter);
152:   ISDestroy(from);
153:   ISDestroy(to);
154:   PetscFree(idx);
155:   VecDestroy(tmp_vec);
156:   return 0;
157: }

161: PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
162: {
163:   PetscScalar    *inptr,*outptr;
164:   PetscInt       i,n,*idx;
166:   IS             from,to;
167:   VecScatter     scatter;
168:   Vec            tmp_in,tmp_out;

170:   /* Get the length of parallel vector */
171:   VecGetSize(globalin,&n);

173:   /* Set the index sets */
174:   PetscMalloc(n*sizeof(PetscInt),&idx);
175:   for(i=0; i<n; i++) idx[i]=i;
176: 
177:   /* Create local sequential vectors */
178:   VecCreateSeq(PETSC_COMM_SELF,n,&tmp_in);
179:   VecDuplicate(tmp_in,&tmp_out);

181:   /* Create scatter context */
182:   ISCreateGeneral(PETSC_COMM_SELF,n,idx,&from);
183:   ISCreateGeneral(PETSC_COMM_SELF,n,idx,&to);
184:   VecScatterCreate(globalin,from,tmp_in,to,&scatter);
185:   VecScatterBegin(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
186:   VecScatterEnd(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
187:   VecScatterDestroy(scatter);

189:   /*Extract income array */
190:   VecGetArray(tmp_in,&inptr);

192:   /* Extract outcome array*/
193:   VecGetArray(tmp_out,&outptr);

195:   outptr[0] = 2.0*inptr[0]+inptr[1];
196:   outptr[1] = inptr[0]+2.0*inptr[1]+inptr[2];
197:   outptr[2] = inptr[1]+2.0*inptr[2];

199:   VecRestoreArray(tmp_in,&inptr);
200:   VecRestoreArray(tmp_out,&outptr);

202:   VecScatterCreate(tmp_out,from,globalout,to,&scatter);
203:   VecScatterBegin(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
204:   VecScatterEnd(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);

206:   /* Destroy idx aand scatter */
207:   ISDestroy(from);
208:   ISDestroy(to);
209:   VecScatterDestroy(scatter);
210:   VecDestroy(tmp_in);
211:   VecDestroy(tmp_out);
212:   PetscFree(idx);
213:   return 0;
214: }

218: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec x,Mat *AA,Mat *BB,MatStructure *str,void *ctx)
219: {
220:   Mat            A = *AA;
221:   PetscScalar    v[3],*tmp;
222:   PetscInt       idx[3],i;
224: 
225:   *str = SAME_NONZERO_PATTERN;

227:   idx[0]=0; idx[1]=1; idx[2]=2;
228:   VecGetArray(x,&tmp);

230:   i = 0;
231:   v[0] = 2.0; v[1] = 1.0; v[2] = 0.0;
232:   MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);

234:   i = 1;
235:   v[0] = 1.0; v[1] = 2.0; v[2] = 1.0;
236:   MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
237: 
238:   i = 2;
239:   v[0]= 0.0; v[1] = 1.0; v[2] = 2.0;
240:   MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);

242:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
243:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

245:   VecRestoreArray(x,&tmp);

247:   return 0;
248: }

250: /*
251:       The exact solutions 
252: */
253: PetscReal solx(PetscReal t)
254: {
255:   return exp((2.0 - sqrt(2.0))*t)/2.0 - exp((2.0 - sqrt(2.0))*t)/(2.0*sqrt(2.0)) +
256:          exp((2.0 + sqrt(2.0))*t)/2.0 + exp((2.0 + sqrt(2.0))*t)/(2.0*sqrt(2.0));
257: }

259: PetscReal soly(PetscReal t)
260: {
261:   return exp((2.0 - sqrt(2.0))*t)/2.0 - exp((2.0 - sqrt(2.0))*t)/sqrt(2.0) +
262:          exp((2.0 + sqrt(2.0))*t)/2.0 + exp((2.0 + sqrt(2.0))*t)/sqrt(2.0);
263: }
264: 
265: PetscReal solz(PetscReal t)
266: {
267:   return exp((2.0 - sqrt(2.0))*t)/2.0 - exp((2.0 - sqrt(2.0))*t)/(2.0*sqrt(2.0)) +
268:          exp((2.0 + sqrt(2.0))*t)/2.0 + exp((2.0 + sqrt(2.0))*t)/(2.0*sqrt(2.0));
269: }