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 PetscFunctionBeginUser; 39 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 40 PetscCall(PetscOptionsGetInt(NULL,NULL,"-time",&time_steps,NULL)); 41 42 /* set initial conditions */ 43 PetscCall(VecCreate(PETSC_COMM_WORLD,&global)); 44 PetscCall(VecSetSizes(global,PETSC_DECIDE,3)); 45 PetscCall(VecSetFromOptions(global)); 46 PetscCall(Initial(global,NULL)); 47 48 /* make timestep context */ 49 PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 50 PetscCall(TSSetProblemType(ts,TS_NONLINEAR)); 51 PetscCall(TSMonitorSet(ts,Monitor,NULL,NULL)); 52 dt = 0.001; 53 54 /* 55 The user provides the RHS and Jacobian 56 */ 57 PetscCall(TSSetRHSFunction(ts,NULL,RHSFunction,NULL)); 58 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 59 PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,3,3)); 60 PetscCall(MatSetFromOptions(A)); 61 PetscCall(MatSetUp(A)); 62 PetscCall(RHSJacobian(ts,0.0,global,A,A,NULL)); 63 PetscCall(TSSetRHSJacobian(ts,A,A,RHSJacobian,NULL)); 64 65 PetscCall(MatCreateShell(PETSC_COMM_WORLD,3,3,3,3,NULL,&S)); 66 PetscCall(MatShellSetOperation(S,MATOP_MULT,(void (*)(void))MyMatMult)); 67 PetscCall(TSSetRHSJacobian(ts,S,A,RHSJacobian,NULL)); 68 69 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP)); 70 PetscCall(TSSetFromOptions(ts)); 71 72 PetscCall(TSSetTimeStep(ts,dt)); 73 PetscCall(TSSetMaxSteps(ts,time_steps)); 74 PetscCall(TSSetMaxTime(ts,1)); 75 PetscCall(TSSetSolution(ts,global)); 76 77 PetscCall(TSSolve(ts,global)); 78 PetscCall(TSGetSolveTime(ts,&ftime)); 79 PetscCall(TSGetStepNumber(ts,&steps)); 80 81 /* free the memories */ 82 83 PetscCall(TSDestroy(&ts)); 84 PetscCall(VecDestroy(&global)); 85 PetscCall(MatDestroy(&A)); 86 PetscCall(MatDestroy(&S)); 87 88 PetscCall(PetscFinalize()); 89 return 0; 90 } 91 92 PetscErrorCode MyMatMult(Mat S,Vec x,Vec y) 93 { 94 const PetscScalar *inptr; 95 PetscScalar *outptr; 96 97 PetscFunctionBegin; 98 PetscCall(VecGetArrayRead(x,&inptr)); 99 PetscCall(VecGetArrayWrite(y,&outptr)); 100 101 outptr[0] = 2.0*inptr[0]+inptr[1]; 102 outptr[1] = inptr[0]+2.0*inptr[1]+inptr[2]; 103 outptr[2] = inptr[1]+2.0*inptr[2]; 104 105 PetscCall(VecRestoreArrayRead(x,&inptr)); 106 PetscCall(VecRestoreArrayWrite(y,&outptr)); 107 PetscFunctionReturn(0); 108 } 109 110 /* this test problem has initial values (1,1,1). */ 111 PetscErrorCode Initial(Vec global,void *ctx) 112 { 113 PetscScalar *localptr; 114 PetscInt i,mybase,myend,locsize; 115 116 /* determine starting point of each processor */ 117 PetscCall(VecGetOwnershipRange(global,&mybase,&myend)); 118 PetscCall(VecGetLocalSize(global,&locsize)); 119 120 /* Initialize the array */ 121 PetscCall(VecGetArrayWrite(global,&localptr)); 122 for (i=0; i<locsize; i++) localptr[i] = 1.0; 123 124 if (mybase == 0) localptr[0]=1.0; 125 126 PetscCall(VecRestoreArrayWrite(global,&localptr)); 127 return 0; 128 } 129 130 PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec global,void *ctx) 131 { 132 VecScatter scatter; 133 IS from,to; 134 PetscInt i,n,*idx; 135 Vec tmp_vec; 136 const PetscScalar *tmp; 137 138 /* Get the size of the vector */ 139 PetscCall(VecGetSize(global,&n)); 140 141 /* Set the index sets */ 142 PetscCall(PetscMalloc1(n,&idx)); 143 for (i=0; i<n; i++) idx[i]=i; 144 145 /* Create local sequential vectors */ 146 PetscCall(VecCreateSeq(PETSC_COMM_SELF,n,&tmp_vec)); 147 148 /* Create scatter context */ 149 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from)); 150 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to)); 151 PetscCall(VecScatterCreate(global,from,tmp_vec,to,&scatter)); 152 PetscCall(VecScatterBegin(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD)); 153 PetscCall(VecScatterEnd(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD)); 154 155 PetscCall(VecGetArrayRead(tmp_vec,&tmp)); 156 PetscCall(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]))); 158 PetscCall(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)))); 160 PetscCall(VecRestoreArrayRead(tmp_vec,&tmp)); 161 PetscCall(VecScatterDestroy(&scatter)); 162 PetscCall(ISDestroy(&from)); 163 PetscCall(ISDestroy(&to)); 164 PetscCall(PetscFree(idx)); 165 PetscCall(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 PetscCall(VecGetSize(globalin,&n)); 180 181 /* Set the index sets */ 182 PetscCall(PetscMalloc1(n,&idx)); 183 for (i=0; i<n; i++) idx[i]=i; 184 185 /* Create local sequential vectors */ 186 PetscCall(VecCreateSeq(PETSC_COMM_SELF,n,&tmp_in)); 187 PetscCall(VecDuplicate(tmp_in,&tmp_out)); 188 189 /* Create scatter context */ 190 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from)); 191 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to)); 192 PetscCall(VecScatterCreate(globalin,from,tmp_in,to,&scatter)); 193 PetscCall(VecScatterBegin(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD)); 194 PetscCall(VecScatterEnd(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD)); 195 PetscCall(VecScatterDestroy(&scatter)); 196 197 /*Extract income array */ 198 PetscCall(VecGetArrayRead(tmp_in,&inptr)); 199 200 /* Extract outcome array*/ 201 PetscCall(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 PetscCall(VecRestoreArrayRead(tmp_in,&inptr)); 208 PetscCall(VecRestoreArrayWrite(tmp_out,&outptr)); 209 210 PetscCall(VecScatterCreate(tmp_out,from,globalout,to,&scatter)); 211 PetscCall(VecScatterBegin(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD)); 212 PetscCall(VecScatterEnd(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD)); 213 214 /* Destroy idx aand scatter */ 215 PetscCall(ISDestroy(&from)); 216 PetscCall(ISDestroy(&to)); 217 PetscCall(VecScatterDestroy(&scatter)); 218 PetscCall(VecDestroy(&tmp_in)); 219 PetscCall(VecDestroy(&tmp_out)); 220 PetscCall(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 PetscCall(VecGetArrayRead(x,&tmp)); 232 233 i = 0; 234 v[0] = 2.0; v[1] = 1.0; v[2] = 0.0; 235 PetscCall(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 PetscCall(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 PetscCall(MatSetValues(BB,1,&i,3,idx,v,INSERT_VALUES)); 244 245 PetscCall(MatAssemblyBegin(BB,MAT_FINAL_ASSEMBLY)); 246 PetscCall(MatAssemblyEnd(BB,MAT_FINAL_ASSEMBLY)); 247 248 if (A != BB) { 249 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 250 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 251 } 252 PetscCall(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