1 static char help[] = "Basic problem for multi-rate method.\n"; 2 3 /*F 4 5 \begin{eqnarray} 6 ys' = -2.0*\frac{-1.0+ys^2.0-\cos(t)}{2.0*ys}+0.05*\frac{-2.0+yf^2-\cos(5.0*t)}{2.0*yf}-\frac{\sin(t)}{2.0*ys}\\ 7 yf' = 0.05*\frac{-1.0+ys^2-\cos(t)}{2.0*ys}-\frac{-2.0+yf^2-\cos(5.0*t)}{2.0*yf}-5.0*\frac{\sin(5.0*t)}{2.0*yf}\\ 8 \end{eqnarray} 9 10 F*/ 11 12 #include <petscts.h> 13 14 typedef struct { 15 PetscReal 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] = -2.0 * (-1.0 + u[0] * u[0] - PetscCosScalar(t)) / (2.0 * u[0]) + 0.05 * (-2.0 + u[1] * u[1] - PetscCosScalar(5.0 * t)) / (2.0 * u[1]) - PetscSinScalar(t) / (2.0 * u[0]); 26 f[1] = 0.05 * (-1.0 + u[0] * u[0] - PetscCosScalar(t)) / (2.0 * u[0]) - (-2.0 + u[1] * u[1] - PetscCosScalar(5.0 * t)) / (2.0 * u[1]) - 5.0 * PetscSinScalar(5.0 * t) / (2.0 * u[1]); 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] = -2.0 * (-1.0 + u[0] * u[0] - PetscCosScalar(t)) / (2.0 * u[0]) + 0.05 * (-2.0 + u[1] * u[1] - PetscCosScalar(5.0 * t)) / (2.0 * u[1]) - PetscSinScalar(t) / (2.0 * u[0]); 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] = 0.05 * (-1.0 + u[0] * u[0] - PetscCosScalar(t)) / (2.0 * u[0]) - (-2.0 + u[1] * u[1] - PetscCosScalar(5.0 * t)) / (2.0 * u[1]) - 5.0 * PetscSinScalar(5.0 * t) / (2.0 * u[1]); 53 PetscCall(VecRestoreArrayRead(U, &u)); 54 PetscCall(VecRestoreArray(F, &f)); 55 PetscFunctionReturn(0); 56 } 57 58 static PetscErrorCode sol_true(PetscReal t, Vec U) { 59 PetscScalar *u; 60 61 PetscFunctionBegin; 62 PetscCall(VecGetArray(U, &u)); 63 u[0] = PetscSqrtScalar(1.0 + PetscCosScalar(t)); 64 u[1] = PetscSqrtScalar(2.0 + PetscCosScalar(5.0 * t)); 65 PetscCall(VecRestoreArray(U, &u)); 66 PetscFunctionReturn(0); 67 } 68 69 int main(int argc, char **argv) { 70 TS ts; /* ODE integrator */ 71 Vec U; /* solution will be stored here */ 72 Vec Utrue; 73 PetscMPIInt size; 74 AppCtx ctx; 75 PetscScalar *u; 76 IS iss; 77 IS isf; 78 PetscInt *indicess; 79 PetscInt *indicesf; 80 PetscInt n = 2; 81 PetscReal error, tt; 82 83 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 84 Initialize program 85 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 86 PetscFunctionBeginUser; 87 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 88 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 89 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs"); 90 91 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 92 Create index for slow part and fast part 93 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 94 PetscCall(PetscMalloc1(1, &indicess)); 95 indicess[0] = 0; 96 PetscCall(PetscMalloc1(1, &indicesf)); 97 indicesf[0] = 1; 98 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 1, indicess, PETSC_COPY_VALUES, &iss)); 99 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 1, indicesf, PETSC_COPY_VALUES, &isf)); 100 101 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 102 Create necessary vector 103 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 104 PetscCall(VecCreate(PETSC_COMM_WORLD, &U)); 105 PetscCall(VecSetSizes(U, n, PETSC_DETERMINE)); 106 PetscCall(VecSetFromOptions(U)); 107 PetscCall(VecDuplicate(U, &Utrue)); 108 PetscCall(VecCopy(U, Utrue)); 109 110 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 111 Set initial condition 112 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 113 PetscCall(VecGetArray(U, &u)); 114 u[0] = PetscSqrtScalar(2.0); 115 u[1] = PetscSqrtScalar(3.0); 116 PetscCall(VecRestoreArray(U, &u)); 117 118 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 119 Create timestepping solver context 120 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 121 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 122 PetscCall(TSSetType(ts, TSMPRK)); 123 124 PetscCall(TSSetRHSFunction(ts, NULL, (TSRHSFunction)RHSFunction, &ctx)); 125 PetscCall(TSRHSSplitSetIS(ts, "slow", iss)); 126 PetscCall(TSRHSSplitSetIS(ts, "fast", isf)); 127 PetscCall(TSRHSSplitSetRHSFunction(ts, "slow", NULL, (TSRHSFunction)RHSFunctionslow, &ctx)); 128 PetscCall(TSRHSSplitSetRHSFunction(ts, "fast", NULL, (TSRHSFunction)RHSFunctionfast, &ctx)); 129 130 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 131 Set initial conditions 132 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 133 PetscCall(TSSetSolution(ts, U)); 134 135 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 136 Set solver options 137 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 138 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ODE options", ""); 139 { 140 ctx.Tf = 0.3; 141 ctx.dt = 0.01; 142 PetscCall(PetscOptionsScalar("-Tf", "", "", ctx.Tf, &ctx.Tf, NULL)); 143 PetscCall(PetscOptionsScalar("-dt", "", "", ctx.dt, &ctx.dt, NULL)); 144 } 145 PetscOptionsEnd(); 146 PetscCall(TSSetMaxTime(ts, ctx.Tf)); 147 PetscCall(TSSetTimeStep(ts, ctx.dt)); 148 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 149 PetscCall(TSSetFromOptions(ts)); 150 151 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 152 Solve linear system 153 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 154 PetscCall(TSSolve(ts, U)); 155 PetscCall(VecView(U, PETSC_VIEWER_STDOUT_WORLD)); 156 157 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 158 Check the error of the Petsc solution 159 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 160 PetscCall(TSGetTime(ts, &tt)); 161 PetscCall(sol_true(tt, Utrue)); 162 PetscCall(VecAXPY(Utrue, -1.0, U)); 163 PetscCall(VecNorm(Utrue, NORM_2, &error)); 164 165 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 166 Print norm2 error 167 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 168 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "l2 error norm: %g\n", (double)error)); 169 170 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 171 Free work space. All PETSc objects should be destroyed when they are no longer needed. 172 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 173 PetscCall(VecDestroy(&U)); 174 PetscCall(TSDestroy(&ts)); 175 PetscCall(VecDestroy(&Utrue)); 176 PetscCall(ISDestroy(&iss)); 177 PetscCall(ISDestroy(&isf)); 178 PetscCall(PetscFree(indicess)); 179 PetscCall(PetscFree(indicesf)); 180 PetscCall(PetscFinalize()); 181 return 0; 182 } 183 184 /*TEST 185 build: 186 requires: !complex 187 188 test: 189 190 TEST*/ 191