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 PetscCall(VecGetArrayRead(U, &u)); 25 PetscCall(VecGetArray(F, &f)); 26 f[0] = u[0] / ctx->a; 27 f[1] = u[0] * PetscCosScalar(t * ctx->b); 28 PetscCall(VecRestoreArrayRead(U, &u)); 29 PetscCall(VecRestoreArray(F, &f)); 30 PetscFunctionReturn(PETSC_SUCCESS); 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] = u[0] / ctx->a; 42 PetscCall(VecRestoreArrayRead(U, &u)); 43 PetscCall(VecRestoreArray(F, &f)); 44 PetscFunctionReturn(PETSC_SUCCESS); 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] = u[0] * PetscCosScalar(t * ctx->b); 56 PetscCall(VecRestoreArrayRead(U, &u)); 57 PetscCall(VecRestoreArray(F, &f)); 58 PetscFunctionReturn(PETSC_SUCCESS); 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] = 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 PetscCall(VecRestoreArray(U, &u)); 73 PetscFunctionReturn(PETSC_SUCCESS); 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 PetscReal error, tt; 90 91 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 92 Initialize program 93 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 94 PetscFunctionBeginUser; 95 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 96 PetscCallMPI(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 PetscCall(PetscMalloc1(1, &indicess)); 103 indicess[0] = 0; 104 PetscCall(PetscMalloc1(1, &indicesf)); 105 indicesf[0] = 1; 106 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 1, indicess, PETSC_COPY_VALUES, &iss)); 107 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 1, indicesf, PETSC_COPY_VALUES, &isf)); 108 109 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 110 Create necessary vector 111 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 112 PetscCall(VecCreate(PETSC_COMM_WORLD, &U)); 113 PetscCall(VecSetSizes(U, n, PETSC_DETERMINE)); 114 PetscCall(VecSetFromOptions(U)); 115 PetscCall(VecDuplicate(U, &Utrue)); 116 PetscCall(VecCopy(U, Utrue)); 117 118 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 119 Set runtime options 120 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 121 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ODE options", ""); 122 { 123 ctx.a = 2.0; 124 ctx.b = 25.0; 125 PetscCall(PetscOptionsReal("-a", "", "", ctx.a, &ctx.a, NULL)); 126 PetscCall(PetscOptionsReal("-b", "", "", ctx.b, &ctx.b, NULL)); 127 ctx.Tf = 2; 128 ctx.dt = 0.01; 129 PetscCall(PetscOptionsReal("-Tf", "", "", ctx.Tf, &ctx.Tf, NULL)); 130 PetscCall(PetscOptionsReal("-dt", "", "", ctx.dt, &ctx.dt, NULL)); 131 } 132 PetscOptionsEnd(); 133 134 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 135 Initialize the solution 136 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 137 PetscCall(VecGetArray(U, &u)); 138 u[0] = 1.0; 139 u[1] = ctx.a / (1 + ctx.a * ctx.a * ctx.b * ctx.b); 140 PetscCall(VecRestoreArray(U, &u)); 141 142 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 143 Create timestepping solver context 144 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 145 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 146 PetscCall(TSSetType(ts, TSMPRK)); 147 148 PetscCall(TSSetRHSFunction(ts, NULL, (TSRHSFunctionFn *)RHSFunction, &ctx)); 149 PetscCall(TSRHSSplitSetIS(ts, "slow", iss)); 150 PetscCall(TSRHSSplitSetIS(ts, "fast", isf)); 151 PetscCall(TSRHSSplitSetRHSFunction(ts, "slow", NULL, (TSRHSFunctionFn *)RHSFunctionslow, &ctx)); 152 PetscCall(TSRHSSplitSetRHSFunction(ts, "fast", NULL, (TSRHSFunctionFn *)RHSFunctionfast, &ctx)); 153 154 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 155 Set initial conditions 156 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 157 PetscCall(TSSetSolution(ts, U)); 158 159 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 160 Set solver options 161 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 162 PetscCall(TSSetMaxTime(ts, ctx.Tf)); 163 PetscCall(TSSetTimeStep(ts, ctx.dt)); 164 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 165 PetscCall(TSSetFromOptions(ts)); 166 167 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 168 Solve linear system 169 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 170 PetscCall(TSSolve(ts, U)); 171 PetscCall(VecView(U, PETSC_VIEWER_STDOUT_WORLD)); 172 173 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 174 Check the error of the Petsc solution 175 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 176 PetscCall(TSGetTime(ts, &tt)); 177 PetscCall(sol_true(tt, Utrue, &ctx)); 178 PetscCall(VecAXPY(Utrue, -1.0, U)); 179 PetscCall(VecNorm(Utrue, NORM_2, &error)); 180 181 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 182 Print norm2 error 183 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 184 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "l2 error norm = %g\n", (double)error)); 185 186 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 187 Free work space. All PETSc objects should be destroyed when they are no longer needed. 188 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 189 PetscCall(VecDestroy(&U)); 190 PetscCall(VecDestroy(&Utrue)); 191 PetscCall(TSDestroy(&ts)); 192 PetscCall(ISDestroy(&iss)); 193 PetscCall(ISDestroy(&isf)); 194 PetscCall(PetscFree(indicess)); 195 PetscCall(PetscFree(indicesf)); 196 PetscCall(PetscFinalize()); 197 return 0; 198 } 199 200 /*TEST 201 build: 202 requires: !complex 203 204 test: 205 206 TEST*/ 207