Actual source code: ex2.c

  1: /*$Id: ex2.c,v 1.36 2001/08/07 21:31:27 bsmith Exp $*/
  2: /*
  3:        Formatted test for TS routines.

  5:           Solves U_t=F(t,u)
  6:           Where:
  7:           
  8:                   [2*u1+u2
  9:           F(t,u)= [u1+2*u2+u3
 10:                   [   u2+2*u3
 11:        We can compare the solutions from euler, beuler and PVODE to
 12:        see what is the difference.

 14: */

 16: static char help[] = "Solves a nonlinear ODE. \n\n";

 18:  #include petscsys.h
 19:  #include petscts.h
 20:  #include petscpc.h

 22: extern int RHSFunction(TS,PetscReal,Vec,Vec,void*);
 23: extern int RHSJacobian(TS,PetscReal,Vec,Mat*,Mat*,MatStructure *,void*);
 24: extern int Monitor(TS,int,PetscReal,Vec,void *);
 25: extern int Initial(Vec,void *);

 27: extern PetscReal solx(PetscReal);
 28: extern PetscReal soly(PetscReal);
 29: extern PetscReal solz(PetscReal);

 33: int main(int argc,char **argv)
 34: {
 35:   int           ierr,time_steps = 100,steps,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:   TSSetMonitor(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,PETSC_DECIDE,PETSC_DECIDE,3,3,&A);
 65:   MatSetFromOptions(A);
 66:   RHSJacobian(ts,0.0,global,&A,&A,&A_structure,NULL);
 67:   TSSetRHSJacobian(ts,A,A,RHSJacobian,NULL);
 68: 
 69:   TSSetFromOptions(ts);

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


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


 80:   /* free the memories */
 81: 
 82:   TSDestroy(ts);
 83:   VecDestroy(global);
 84:   ierr= MatDestroy(A);

 86:   PetscFinalize();
 87:   return 0;
 88: }

 90: /* -------------------------------------------------------------------*/
 93: /* this test problem has initial values (1,1,1).                      */
 94: int Initial(Vec global,void *ctx)
 95: {
 96:   PetscScalar *localptr;
 97:   int    i,mybase,myend,ierr,locsize;

 99:   /* determine starting point of each processor */
100:   VecGetOwnershipRange(global,&mybase,&myend);
101:   VecGetLocalSize(global,&locsize);

103:   /* Initialize the array */
104:   VecGetArray(global,&localptr);
105:   for (i=0; i<locsize; i++) {
106:     localptr[i] = 1.0;
107:   }
108: 
109:   if (mybase == 0) localptr[0]=1.0;

111:   VecRestoreArray(global,&localptr);
112:   return 0;
113: }

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

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

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

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

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: int RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
160: {
161:   PetscScalar  *inptr,*outptr;
162:   int          i,n,ierr,*idx;
163:   IS           from,to;
164:   VecScatter   scatter;
165:   Vec          tmp_in,tmp_out;

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

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

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

186:   /*Extract income array */
187:   VecGetArray(tmp_in,&inptr);

189:   /* Extract outcome array*/
190:   VecGetArray(tmp_out,&outptr);

192:   outptr[0] = 2.0*inptr[0]+inptr[1];
193:   outptr[1] = inptr[0]+2.0*inptr[1]+inptr[2];
194:   outptr[2] = inptr[1]+2.0*inptr[2];

196:   VecRestoreArray(tmp_in,&inptr);
197:   VecRestoreArray(tmp_out,&outptr);

199:   VecScatterCreate(tmp_out,from,globalout,to,&scatter);
200:   VecScatterBegin(tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD,scatter);
201:   VecScatterEnd(tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD,scatter);

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

215: int RHSJacobian(TS ts,PetscReal t,Vec x,Mat *AA,Mat *BB,MatStructure *str,void *ctx)
216: {
217:   Mat    A = *AA;
218:   PetscScalar v[3],*tmp;
219:   int    idx[3],i,ierr;
220: 
221:   *str = SAME_NONZERO_PATTERN;

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

226:   i = 0;
227:   v[0] = 2.0; v[1] = 1.0; v[2] = 0.0;
228:   MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);

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

238:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
239:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

241:   VecRestoreArray(x,&tmp);

243:   return 0;
244: }

246: /*
247:       The exact solutions 
248: */
249: PetscReal solx(PetscReal t)
250: {
251:   return exp((2.0 - sqrt(2.0))*t)/2.0 - exp((2.0 - sqrt(2.0))*t)/(2.0*sqrt(2.0)) +
252:          exp((2.0 + sqrt(2.0))*t)/2.0 + exp((2.0 + sqrt(2.0))*t)/(2.0*sqrt(2.0));
253: }

255: PetscReal soly(PetscReal t)
256: {
257:   return exp((2.0 - sqrt(2.0))*t)/2.0 - exp((2.0 - sqrt(2.0))*t)/sqrt(2.0) +
258:          exp((2.0 + sqrt(2.0))*t)/2.0 + exp((2.0 + sqrt(2.0))*t)/sqrt(2.0);
259: }
260: 
261: PetscReal solz(PetscReal t)
262: {
263:   return exp((2.0 - sqrt(2.0))*t)/2.0 - exp((2.0 - sqrt(2.0))*t)/(2.0*sqrt(2.0)) +
264:          exp((2.0 + sqrt(2.0))*t)/2.0 + exp((2.0 + sqrt(2.0))*t)/(2.0*sqrt(2.0));
265: }