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