1 /* 2 Formatted test for TS routines. 3 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. 12 13 */ 14 15 static char help[] = "Solves a linear ODE. \n\n"; 16 17 #include <petscts.h> 18 #include <petscpc.h> 19 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*); 24 extern PetscErrorCode MyMatMult(Mat,Vec,Vec); 25 26 extern PetscReal solx(PetscReal); 27 extern PetscReal soly(PetscReal); 28 extern PetscReal solz(PetscReal); 29 30 int main(int argc,char **argv) 31 { 32 PetscInt time_steps = 100,steps; 33 Vec global; 34 PetscReal dt,ftime; 35 TS ts; 36 Mat A = 0,S; 37 38 CHKERRQ(PetscInitialize(&argc,&argv,(char*)0,help)); 39 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-time",&time_steps,NULL)); 40 41 /* set initial conditions */ 42 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&global)); 43 CHKERRQ(VecSetSizes(global,PETSC_DECIDE,3)); 44 CHKERRQ(VecSetFromOptions(global)); 45 CHKERRQ(Initial(global,NULL)); 46 47 /* make timestep context */ 48 CHKERRQ(TSCreate(PETSC_COMM_WORLD,&ts)); 49 CHKERRQ(TSSetProblemType(ts,TS_NONLINEAR)); 50 CHKERRQ(TSMonitorSet(ts,Monitor,NULL,NULL)); 51 dt = 0.001; 52 53 /* 54 The user provides the RHS and Jacobian 55 */ 56 CHKERRQ(TSSetRHSFunction(ts,NULL,RHSFunction,NULL)); 57 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&A)); 58 CHKERRQ(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,3,3)); 59 CHKERRQ(MatSetFromOptions(A)); 60 CHKERRQ(MatSetUp(A)); 61 CHKERRQ(RHSJacobian(ts,0.0,global,A,A,NULL)); 62 CHKERRQ(TSSetRHSJacobian(ts,A,A,RHSJacobian,NULL)); 63 64 CHKERRQ(MatCreateShell(PETSC_COMM_WORLD,3,3,3,3,NULL,&S)); 65 CHKERRQ(MatShellSetOperation(S,MATOP_MULT,(void (*)(void))MyMatMult)); 66 CHKERRQ(TSSetRHSJacobian(ts,S,A,RHSJacobian,NULL)); 67 68 CHKERRQ(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP)); 69 CHKERRQ(TSSetFromOptions(ts)); 70 71 CHKERRQ(TSSetTimeStep(ts,dt)); 72 CHKERRQ(TSSetMaxSteps(ts,time_steps)); 73 CHKERRQ(TSSetMaxTime(ts,1)); 74 CHKERRQ(TSSetSolution(ts,global)); 75 76 CHKERRQ(TSSolve(ts,global)); 77 CHKERRQ(TSGetSolveTime(ts,&ftime)); 78 CHKERRQ(TSGetStepNumber(ts,&steps)); 79 80 /* free the memories */ 81 82 CHKERRQ(TSDestroy(&ts)); 83 CHKERRQ(VecDestroy(&global)); 84 CHKERRQ(MatDestroy(&A)); 85 CHKERRQ(MatDestroy(&S)); 86 87 CHKERRQ(PetscFinalize()); 88 return 0; 89 } 90 91 PetscErrorCode MyMatMult(Mat S,Vec x,Vec y) 92 { 93 const PetscScalar *inptr; 94 PetscScalar *outptr; 95 96 PetscFunctionBegin; 97 CHKERRQ(VecGetArrayRead(x,&inptr)); 98 CHKERRQ(VecGetArrayWrite(y,&outptr)); 99 100 outptr[0] = 2.0*inptr[0]+inptr[1]; 101 outptr[1] = inptr[0]+2.0*inptr[1]+inptr[2]; 102 outptr[2] = inptr[1]+2.0*inptr[2]; 103 104 CHKERRQ(VecRestoreArrayRead(x,&inptr)); 105 CHKERRQ(VecRestoreArrayWrite(y,&outptr)); 106 PetscFunctionReturn(0); 107 } 108 109 /* this test problem has initial values (1,1,1). */ 110 PetscErrorCode Initial(Vec global,void *ctx) 111 { 112 PetscScalar *localptr; 113 PetscInt i,mybase,myend,locsize; 114 115 /* determine starting point of each processor */ 116 CHKERRQ(VecGetOwnershipRange(global,&mybase,&myend)); 117 CHKERRQ(VecGetLocalSize(global,&locsize)); 118 119 /* Initialize the array */ 120 CHKERRQ(VecGetArrayWrite(global,&localptr)); 121 for (i=0; i<locsize; i++) localptr[i] = 1.0; 122 123 if (mybase == 0) localptr[0]=1.0; 124 125 CHKERRQ(VecRestoreArrayWrite(global,&localptr)); 126 return 0; 127 } 128 129 PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec global,void *ctx) 130 { 131 VecScatter scatter; 132 IS from,to; 133 PetscInt i,n,*idx; 134 Vec tmp_vec; 135 PetscErrorCode ierr; 136 const PetscScalar *tmp; 137 138 /* Get the size of the vector */ 139 CHKERRQ(VecGetSize(global,&n)); 140 141 /* Set the index sets */ 142 CHKERRQ(PetscMalloc1(n,&idx)); 143 for (i=0; i<n; i++) idx[i]=i; 144 145 /* Create local sequential vectors */ 146 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,n,&tmp_vec)); 147 148 /* Create scatter context */ 149 CHKERRQ(ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from)); 150 CHKERRQ(ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to)); 151 CHKERRQ(VecScatterCreate(global,from,tmp_vec,to,&scatter)); 152 CHKERRQ(VecScatterBegin(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD)); 153 CHKERRQ(VecScatterEnd(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD)); 154 155 CHKERRQ(VecGetArrayRead(tmp_vec,&tmp)); 156 ierr = PetscPrintf(PETSC_COMM_WORLD,"At t =%14.6e u = %14.6e %14.6e %14.6e \n", 157 (double)time,(double)PetscRealPart(tmp[0]),(double)PetscRealPart(tmp[1]),(double)PetscRealPart(tmp[2]));CHKERRQ(ierr); 158 ierr = PetscPrintf(PETSC_COMM_WORLD,"At t =%14.6e errors = %14.6e %14.6e %14.6e \n", 159 (double)time,(double)PetscRealPart(tmp[0]-solx(time)),(double)PetscRealPart(tmp[1]-soly(time)),(double)PetscRealPart(tmp[2]-solz(time)));CHKERRQ(ierr); 160 CHKERRQ(VecRestoreArrayRead(tmp_vec,&tmp)); 161 CHKERRQ(VecScatterDestroy(&scatter)); 162 CHKERRQ(ISDestroy(&from)); 163 CHKERRQ(ISDestroy(&to)); 164 CHKERRQ(PetscFree(idx)); 165 CHKERRQ(VecDestroy(&tmp_vec)); 166 return 0; 167 } 168 169 PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx) 170 { 171 PetscScalar *outptr; 172 const PetscScalar *inptr; 173 PetscInt i,n,*idx; 174 IS from,to; 175 VecScatter scatter; 176 Vec tmp_in,tmp_out; 177 178 /* Get the length of parallel vector */ 179 CHKERRQ(VecGetSize(globalin,&n)); 180 181 /* Set the index sets */ 182 CHKERRQ(PetscMalloc1(n,&idx)); 183 for (i=0; i<n; i++) idx[i]=i; 184 185 /* Create local sequential vectors */ 186 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,n,&tmp_in)); 187 CHKERRQ(VecDuplicate(tmp_in,&tmp_out)); 188 189 /* Create scatter context */ 190 CHKERRQ(ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from)); 191 CHKERRQ(ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to)); 192 CHKERRQ(VecScatterCreate(globalin,from,tmp_in,to,&scatter)); 193 CHKERRQ(VecScatterBegin(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD)); 194 CHKERRQ(VecScatterEnd(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD)); 195 CHKERRQ(VecScatterDestroy(&scatter)); 196 197 /*Extract income array */ 198 CHKERRQ(VecGetArrayRead(tmp_in,&inptr)); 199 200 /* Extract outcome array*/ 201 CHKERRQ(VecGetArrayWrite(tmp_out,&outptr)); 202 203 outptr[0] = 2.0*inptr[0]+inptr[1]; 204 outptr[1] = inptr[0]+2.0*inptr[1]+inptr[2]; 205 outptr[2] = inptr[1]+2.0*inptr[2]; 206 207 CHKERRQ(VecRestoreArrayRead(tmp_in,&inptr)); 208 CHKERRQ(VecRestoreArrayWrite(tmp_out,&outptr)); 209 210 CHKERRQ(VecScatterCreate(tmp_out,from,globalout,to,&scatter)); 211 CHKERRQ(VecScatterBegin(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD)); 212 CHKERRQ(VecScatterEnd(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD)); 213 214 /* Destroy idx aand scatter */ 215 CHKERRQ(ISDestroy(&from)); 216 CHKERRQ(ISDestroy(&to)); 217 CHKERRQ(VecScatterDestroy(&scatter)); 218 CHKERRQ(VecDestroy(&tmp_in)); 219 CHKERRQ(VecDestroy(&tmp_out)); 220 CHKERRQ(PetscFree(idx)); 221 return 0; 222 } 223 224 PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec x,Mat A,Mat BB,void *ctx) 225 { 226 PetscScalar v[3]; 227 const PetscScalar *tmp; 228 PetscInt idx[3],i; 229 230 idx[0]=0; idx[1]=1; idx[2]=2; 231 CHKERRQ(VecGetArrayRead(x,&tmp)); 232 233 i = 0; 234 v[0] = 2.0; v[1] = 1.0; v[2] = 0.0; 235 CHKERRQ(MatSetValues(BB,1,&i,3,idx,v,INSERT_VALUES)); 236 237 i = 1; 238 v[0] = 1.0; v[1] = 2.0; v[2] = 1.0; 239 CHKERRQ(MatSetValues(BB,1,&i,3,idx,v,INSERT_VALUES)); 240 241 i = 2; 242 v[0] = 0.0; v[1] = 1.0; v[2] = 2.0; 243 CHKERRQ(MatSetValues(BB,1,&i,3,idx,v,INSERT_VALUES)); 244 245 CHKERRQ(MatAssemblyBegin(BB,MAT_FINAL_ASSEMBLY)); 246 CHKERRQ(MatAssemblyEnd(BB,MAT_FINAL_ASSEMBLY)); 247 248 if (A != BB) { 249 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 250 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 251 } 252 CHKERRQ(VecRestoreArrayRead(x,&tmp)); 253 254 return 0; 255 } 256 257 /* 258 The exact solutions 259 */ 260 PetscReal solx(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 } 265 266 PetscReal soly(PetscReal t) 267 { 268 return PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/PetscSqrtReal(2.0) + 269 PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/PetscSqrtReal(2.0); 270 } 271 272 PetscReal solz(PetscReal t) 273 { 274 return PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/(2.0*PetscSqrtReal(2.0)) + 275 PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/(2.0*PetscSqrtReal(2.0)); 276 } 277 278 /*TEST 279 280 test: 281 suffix: euler 282 args: -ts_type euler 283 requires: !single 284 285 test: 286 suffix: beuler 287 args: -ts_type beuler 288 requires: !single 289 290 TEST*/ 291