1 static char help[] = "Basic problem for multi-rate method.\n"; 2 3 /*F 4 5 \begin{eqnarray} 6 ys' = -sin(a*t)\\ 7 yf' = bcos(b*t)ys-sin(b*t)sin(a*t)\\ 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 PetscCall(VecGetArrayRead(U,&u)); 25 PetscCall(VecGetArray(F,&f)); 26 f[0] = -PetscSinScalar(ctx->a*t); 27 f[1] = ctx->b*PetscCosScalar(ctx->b*t)*u[0]-PetscSinScalar(ctx->b*t)*PetscSinScalar(ctx->a*t); 28 PetscCall(VecRestoreArrayRead(U,&u)); 29 PetscCall(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 PetscCall(VecGetArrayRead(U,&u)); 40 PetscCall(VecGetArray(F,&f)); 41 f[0] = -PetscSinScalar(ctx->a*t); 42 PetscCall(VecRestoreArrayRead(U,&u)); 43 PetscCall(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 PetscCall(VecGetArrayRead(U,&u)); 54 PetscCall(VecGetArray(F,&f)); 55 f[0] = ctx->b*PetscCosScalar(ctx->b*t)*u[0]-PetscSinScalar(ctx->b*t)*PetscSinScalar(ctx->a*t); 56 PetscCall(VecRestoreArrayRead(U,&u)); 57 PetscCall(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 PetscCall(VecGetArray(U,&u)); 70 u[0] = PetscCosScalar(ctx->a*t)/ctx->a; 71 u[1] = PetscSinScalar(ctx->b*t)*PetscCosScalar(ctx->a*t)/ctx->a; 72 PetscCall(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 PetscMPIInt size; 82 AppCtx ctx; 83 PetscScalar *u; 84 IS iss; 85 IS isf; 86 PetscInt *indicess; 87 PetscInt *indicesf; 88 PetscInt n=2; 89 PetscScalar error,tt; 90 91 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 92 Initialize program 93 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 94 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 95 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 96 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Only for sequential runs"); 97 98 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 99 Create index for slow part and fast part 100 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 101 PetscCall(PetscMalloc1(1,&indicess)); 102 indicess[0]=0; 103 PetscCall(PetscMalloc1(1,&indicesf)); 104 indicesf[0]=1; 105 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,1,indicess,PETSC_COPY_VALUES,&iss)); 106 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,1,indicesf,PETSC_COPY_VALUES,&isf)); 107 108 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 109 Create necesary vector 110 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 111 PetscCall(VecCreate(PETSC_COMM_WORLD,&U)); 112 PetscCall(VecSetSizes(U,n,PETSC_DETERMINE)); 113 PetscCall(VecSetFromOptions(U)); 114 PetscCall(VecDuplicate(U,&Utrue)); 115 PetscCall(VecCopy(U,Utrue)); 116 117 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 118 Set runtime options 119 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 120 PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ODE options",""); 121 { 122 ctx.a = 1.0; 123 ctx.b = 25.0; 124 PetscCall(PetscOptionsScalar("-a","","",ctx.a,&ctx.a,NULL)); 125 PetscCall(PetscOptionsScalar("-b","","",ctx.b,&ctx.b,NULL)); 126 ctx.Tf = 5.0; 127 ctx.dt = 0.01; 128 PetscCall(PetscOptionsScalar("-Tf","","",ctx.Tf,&ctx.Tf,NULL)); 129 PetscCall(PetscOptionsScalar("-dt","","",ctx.dt,&ctx.dt,NULL)); 130 } 131 PetscOptionsEnd(); 132 133 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 134 Initialize the solution 135 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 136 PetscCall(VecGetArray(U,&u)); 137 u[0] = 1.0/ctx.a; 138 u[1] = 0.0; 139 PetscCall(VecRestoreArray(U,&u)); 140 141 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 142 Create timestepping solver context 143 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 144 PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 145 PetscCall(TSSetType(ts,TSMPRK)); 146 147 PetscCall(TSSetRHSFunction(ts,NULL,(TSRHSFunction)RHSFunction,&ctx)); 148 PetscCall(TSRHSSplitSetIS(ts,"slow",iss)); 149 PetscCall(TSRHSSplitSetIS(ts,"fast",isf)); 150 PetscCall(TSRHSSplitSetRHSFunction(ts,"slow",NULL,(TSRHSFunction)RHSFunctionslow,&ctx)); 151 PetscCall(TSRHSSplitSetRHSFunction(ts,"fast",NULL,(TSRHSFunction)RHSFunctionfast,&ctx)); 152 153 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 154 Set initial conditions 155 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 156 PetscCall(TSSetSolution(ts,U)); 157 158 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 159 Set solver options 160 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 161 PetscCall(TSSetMaxTime(ts,ctx.Tf)); 162 PetscCall(TSSetTimeStep(ts,ctx.dt)); 163 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP)); 164 PetscCall(TSSetFromOptions(ts)); 165 166 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 167 Solve linear system 168 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 169 PetscCall(TSSolve(ts,U)); 170 PetscCall(VecView(U,PETSC_VIEWER_STDOUT_WORLD)); 171 172 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 173 Check the error of the Petsc solution 174 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 175 PetscCall(TSGetTime(ts,&tt)); 176 PetscCall(sol_true(tt,Utrue,&ctx)); 177 PetscCall(VecAXPY(Utrue,-1.0,U)); 178 PetscCall(VecNorm(Utrue,NORM_2,&error)); 179 180 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 181 Print norm2 error 182 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 183 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"l2 error norm: %g\n", (double)PetscAbsScalar(error))); 184 185 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 186 Free work space. All PETSc objects should be destroyed when they are no longer needed. 187 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 188 PetscCall(VecDestroy(&U)); 189 PetscCall(TSDestroy(&ts)); 190 PetscCall(VecDestroy(&Utrue)); 191 PetscCall(ISDestroy(&iss)); 192 PetscCall(ISDestroy(&isf)); 193 PetscCall(PetscFree(indicess)); 194 PetscCall(PetscFree(indicesf)); 195 PetscCall(PetscFinalize()); 196 return 0; 197 } 198 199 /*TEST 200 build: 201 requires: !complex 202 203 test: 204 205 TEST*/ 206