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