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