1 static char help[] = "Basic problem for multi-rate method.\n"; 2 3 /*F 4 5 \begin{eqnarray} 6 ys' = \frac{ys}{a}\\ 7 yf' = ys*cos(bt)\\ 8 \end{eqnarray} 9 10 F*/ 11 12 #include <petscts.h> 13 14 typedef struct { 15 PetscReal a,b,Tf,dt; 16 } AppCtx; 17 18 static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,AppCtx *ctx) 19 { 20 PetscErrorCode ierr; 21 const PetscScalar *u; 22 PetscScalar *f; 23 24 PetscFunctionBegin; 25 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 26 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 27 f[0] = u[0]/ctx->a; 28 f[1] = u[0]*PetscCosScalar(t*ctx->b); 29 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 30 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 31 PetscFunctionReturn(0); 32 } 33 34 static PetscErrorCode RHSFunctionslow(TS ts,PetscReal t,Vec U,Vec F,AppCtx *ctx) 35 { 36 PetscErrorCode ierr; 37 const PetscScalar *u; 38 PetscScalar *f; 39 40 PetscFunctionBegin; 41 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 42 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 43 f[0] = u[0]/ctx->a; 44 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 45 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 46 PetscFunctionReturn(0); 47 } 48 49 static PetscErrorCode RHSFunctionfast(TS ts,PetscReal t,Vec U,Vec F,AppCtx *ctx) 50 { 51 PetscErrorCode ierr; 52 const PetscScalar *u; 53 PetscScalar *f; 54 55 PetscFunctionBegin; 56 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 57 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 58 f[0] = u[0]*PetscCosScalar(t*ctx->b); 59 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 60 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 61 PetscFunctionReturn(0); 62 } 63 64 /* 65 Define the analytic solution for check method easily 66 */ 67 static PetscErrorCode sol_true(PetscReal t, Vec U,AppCtx *ctx) 68 { 69 PetscErrorCode ierr; 70 PetscScalar *u; 71 72 PetscFunctionBegin; 73 ierr = VecGetArray(U,&u);CHKERRQ(ierr); 74 u[0] = PetscExpScalar(t/ctx->a); 75 u[1] = (ctx->a*PetscCosScalar(ctx->b*t)+ctx->a*ctx->a*ctx->b*PetscSinScalar(ctx->b*t))*PetscExpScalar(t/ctx->a)/(1+ctx->a*ctx->a*ctx->b*ctx->b); 76 ierr = VecRestoreArray(U,&u);CHKERRQ(ierr); 77 PetscFunctionReturn(0); 78 } 79 80 int main(int argc,char **argv) 81 { 82 TS ts; /* ODE integrator */ 83 Vec U; /* solution will be stored here */ 84 Vec Utrue; 85 PetscErrorCode ierr; 86 PetscMPIInt size; 87 AppCtx ctx; 88 PetscScalar *u; 89 IS iss; 90 IS isf; 91 PetscInt *indicess; 92 PetscInt *indicesf; 93 PetscInt n=2; 94 PetscReal error,tt; 95 96 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 97 Initialize program 98 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 99 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 100 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 101 if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs"); 102 103 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 104 Create index for slow part and fast part 105 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 106 ierr = PetscMalloc1(1,&indicess);CHKERRQ(ierr); 107 indicess[0] = 0; 108 ierr = PetscMalloc1(1,&indicesf);CHKERRQ(ierr); 109 indicesf[0] = 1; 110 ierr = ISCreateGeneral(PETSC_COMM_SELF,1,indicess,PETSC_COPY_VALUES,&iss);CHKERRQ(ierr); 111 ierr = ISCreateGeneral(PETSC_COMM_SELF,1,indicesf,PETSC_COPY_VALUES,&isf);CHKERRQ(ierr); 112 113 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 114 Create necesary vector 115 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 116 ierr = VecCreate(PETSC_COMM_WORLD,&U);CHKERRQ(ierr); 117 ierr = VecSetSizes(U,n,PETSC_DETERMINE);CHKERRQ(ierr); 118 ierr = VecSetFromOptions(U);CHKERRQ(ierr); 119 ierr = VecDuplicate(U,&Utrue);CHKERRQ(ierr); 120 ierr = VecCopy(U,Utrue);CHKERRQ(ierr); 121 122 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 123 Set runtime options 124 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 125 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ODE options","");CHKERRQ(ierr); 126 { 127 ctx.a = 2.0; 128 ctx.b = 25.0; 129 ierr = PetscOptionsReal("-a","","",ctx.a,&ctx.a,NULL);CHKERRQ(ierr); 130 ierr = PetscOptionsReal("-b","","",ctx.b,&ctx.b,NULL);CHKERRQ(ierr); 131 ctx.Tf = 2; 132 ctx.dt = 0.01; 133 ierr = PetscOptionsReal("-Tf","","",ctx.Tf,&ctx.Tf,NULL);CHKERRQ(ierr); 134 ierr = PetscOptionsReal("-dt","","",ctx.dt,&ctx.dt,NULL);CHKERRQ(ierr); 135 } 136 ierr = PetscOptionsEnd();CHKERRQ(ierr); 137 138 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 139 Initialize the solution 140 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 141 ierr = VecGetArray(U,&u);CHKERRQ(ierr); 142 u[0] = 1.0; 143 u[1] = ctx.a/(1+ctx.a*ctx.a*ctx.b*ctx.b); 144 ierr = VecRestoreArray(U,&u);CHKERRQ(ierr); 145 146 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 147 Create timestepping solver context 148 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 149 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 150 ierr = TSSetType(ts,TSMPRK);CHKERRQ(ierr); 151 152 ierr = TSSetRHSFunction(ts,NULL,(TSRHSFunction)RHSFunction,&ctx);CHKERRQ(ierr); 153 ierr = TSRHSSplitSetIS(ts,"slow",iss);CHKERRQ(ierr); 154 ierr = TSRHSSplitSetIS(ts,"fast",isf);CHKERRQ(ierr); 155 ierr = TSRHSSplitSetRHSFunction(ts,"slow",NULL,(TSRHSFunction)RHSFunctionslow,&ctx);CHKERRQ(ierr); 156 ierr = TSRHSSplitSetRHSFunction(ts,"fast",NULL,(TSRHSFunction)RHSFunctionfast,&ctx);CHKERRQ(ierr); 157 158 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 159 Set initial conditions 160 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 161 ierr = TSSetSolution(ts,U);CHKERRQ(ierr); 162 163 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 164 Set solver options 165 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 166 ierr = TSSetMaxTime(ts,ctx.Tf);CHKERRQ(ierr); 167 ierr = TSSetTimeStep(ts,ctx.dt);CHKERRQ(ierr); 168 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr); 169 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 170 171 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 172 Solve linear system 173 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 174 ierr = TSSolve(ts,U);CHKERRQ(ierr); 175 ierr = VecView(U,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 176 177 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 178 Check the error of the Petsc solution 179 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 180 ierr = TSGetTime(ts,&tt);CHKERRQ(ierr); 181 ierr = sol_true(tt,Utrue,&ctx);CHKERRQ(ierr); 182 ierr = VecAXPY(Utrue,-1.0,U);CHKERRQ(ierr); 183 ierr = VecNorm(Utrue,NORM_2,&error); 184 185 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 186 Print norm2 error 187 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 188 ierr = PetscPrintf(PETSC_COMM_WORLD,"l2 error norm = %g\n", error);CHKERRQ(ierr); 189 190 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 191 Free work space. All PETSc objects should be destroyed when they are no longer needed. 192 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 193 ierr = VecDestroy(&U);CHKERRQ(ierr); 194 ierr = VecDestroy(&Utrue);CHKERRQ(ierr); 195 ierr = TSDestroy(&ts);CHKERRQ(ierr); 196 ierr = ISDestroy(&iss);CHKERRQ(ierr); 197 ierr = ISDestroy(&isf);CHKERRQ(ierr); 198 ierr = PetscFree(indicess);CHKERRQ(ierr); 199 ierr = PetscFree(indicesf);CHKERRQ(ierr); 200 ierr = PetscFinalize(); 201 return ierr; 202 } 203 204 /*TEST 205 build: 206 requires: !complex 207 208 test: 209 210 TEST*/ 211