1 static char help[] = "A fast-slow system for testing ARKIMEX.\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 /* 13 This example demonstrates how to use ARKIMEX for solving a fast-slow system. The system is partitioned additively and component-wise at the same time. 14 ys stands for the slow component and yf stands for the fast component. On the RHS for yf, only the term -\frac{-2.0+yf^2-\cos(5.0*t)}{2.0*yf} is treated implicitly while the rest is treated explicitly. 15 */ 16 17 #include <petscts.h> 18 19 typedef struct { 20 PetscReal Tf, dt; 21 } AppCtx; 22 23 static PetscErrorCode RHSFunctionslow(TS ts, PetscReal t, Vec U, Vec F, AppCtx *ctx) 24 { 25 const PetscScalar *u; 26 PetscScalar *f; 27 28 PetscFunctionBegin; 29 PetscCall(VecGetArrayRead(U, &u)); 30 PetscCall(VecGetArray(F, &f)); 31 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]); 32 PetscCall(VecRestoreArrayRead(U, &u)); 33 PetscCall(VecRestoreArray(F, &f)); 34 PetscFunctionReturn(PETSC_SUCCESS); 35 } 36 37 static PetscErrorCode RHSFunctionfast(TS ts, PetscReal t, Vec U, Vec F, AppCtx *ctx) 38 { 39 const PetscScalar *u; 40 PetscScalar *f; 41 42 PetscFunctionBegin; 43 PetscCall(VecGetArrayRead(U, &u)); 44 PetscCall(VecGetArray(F, &f)); 45 f[0] = 0.05 * (-1.0 + u[0] * u[0] - PetscCosScalar(t)) / (2.0 * u[0]) - 5.0 * PetscSinScalar(5.0 * t) / (2.0 * u[1]); 46 PetscCall(VecRestoreArrayRead(U, &u)); 47 PetscCall(VecRestoreArray(F, &f)); 48 PetscFunctionReturn(PETSC_SUCCESS); 49 } 50 51 static PetscErrorCode IFunctionfast(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, AppCtx *ctx) 52 { 53 const PetscScalar *u, *udot; 54 PetscScalar *f; 55 56 PetscFunctionBegin; 57 PetscCall(VecGetArrayRead(U, &u)); 58 PetscCall(VecGetArrayRead(Udot, &udot)); 59 PetscCall(VecGetArray(F, &f)); 60 f[0] = udot[0] + (-2.0 + u[1] * u[1] - PetscCosScalar(5.0 * t)) / (2.0 * u[1]); 61 PetscCall(VecRestoreArrayRead(Udot, &udot)); 62 PetscCall(VecRestoreArrayRead(U, &u)); 63 PetscCall(VecRestoreArray(F, &f)); 64 PetscFunctionReturn(PETSC_SUCCESS); 65 } 66 67 static PetscErrorCode IJacobianfast(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, AppCtx *ctx) 68 { 69 PetscInt rowcol[] = {0}; 70 const PetscScalar *u; 71 PetscScalar J[1][1]; 72 73 PetscFunctionBeginUser; 74 PetscCall(VecGetArrayRead(U, &u)); 75 J[0][0] = a + (2.0 + PetscCosScalar(5.0 * t)) / (2.0 * u[1] * u[1]) + 0.5; 76 PetscCall(MatSetValues(B, 1, rowcol, 1, rowcol, &J[0][0], INSERT_VALUES)); 77 PetscCall(VecRestoreArrayRead(U, &u)); 78 79 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 80 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 81 if (A != B) { 82 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 83 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 84 } 85 PetscFunctionReturn(PETSC_SUCCESS); 86 } 87 88 /* 89 Define the analytic solution for check method easily 90 */ 91 static PetscErrorCode sol_true(PetscReal t, Vec U) 92 { 93 PetscScalar *u; 94 95 PetscFunctionBegin; 96 PetscCall(VecGetArray(U, &u)); 97 u[0] = PetscSqrtScalar(1.0 + PetscCosScalar(t)); 98 u[1] = PetscSqrtScalar(2.0 + PetscCosScalar(5.0 * t)); 99 PetscCall(VecRestoreArray(U, &u)); 100 PetscFunctionReturn(PETSC_SUCCESS); 101 } 102 103 int main(int argc, char **argv) 104 { 105 TS ts; /* ODE integrator */ 106 Vec U; /* solution will be stored here */ 107 Vec Utrue; 108 Mat A; 109 PetscMPIInt size; 110 AppCtx ctx; 111 PetscScalar *u; 112 IS iss; 113 IS isf; 114 PetscInt *indicess; 115 PetscInt *indicesf; 116 PetscInt n = 2; 117 PetscScalar error, tt; 118 119 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 120 Initialize program 121 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 122 PetscFunctionBeginUser; 123 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 124 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 125 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs"); 126 127 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 128 Create index for slow part and fast part 129 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 130 PetscCall(PetscMalloc1(1, &indicess)); 131 indicess[0] = 0; 132 PetscCall(PetscMalloc1(1, &indicesf)); 133 indicesf[0] = 1; 134 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 1, indicess, PETSC_COPY_VALUES, &iss)); 135 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 1, indicesf, PETSC_COPY_VALUES, &isf)); 136 137 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 138 Create necessary vector 139 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 140 PetscCall(VecCreate(PETSC_COMM_WORLD, &U)); 141 PetscCall(VecSetSizes(U, n, PETSC_DETERMINE)); 142 PetscCall(VecSetFromOptions(U)); 143 PetscCall(VecDuplicate(U, &Utrue)); 144 PetscCall(VecCopy(U, Utrue)); 145 146 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 147 Set runtime options 148 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 149 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ODE options", ""); 150 ctx.Tf = 0.3; 151 ctx.dt = 0.01; 152 PetscCall(PetscOptionsScalar("-Tf", "", "", ctx.Tf, &ctx.Tf, NULL)); 153 PetscCall(PetscOptionsScalar("-dt", "", "", ctx.dt, &ctx.dt, NULL)); 154 PetscOptionsEnd(); 155 156 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 157 Initialize the solution 158 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 159 PetscCall(VecGetArray(U, &u)); 160 u[0] = PetscSqrtScalar(2.0); 161 u[1] = PetscSqrtScalar(3.0); 162 PetscCall(VecRestoreArray(U, &u)); 163 164 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 165 Create timestepping solver context 166 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 167 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 168 PetscCall(TSSetType(ts, TSARKIMEX)); 169 PetscCall(TSARKIMEXSetFastSlowSplit(ts, PETSC_TRUE)); 170 171 PetscCall(TSRHSSplitSetIS(ts, "slow", iss)); 172 PetscCall(TSRHSSplitSetIS(ts, "fast", isf)); 173 PetscCall(TSRHSSplitSetRHSFunction(ts, "slow", NULL, (TSRHSFunctionFn *)RHSFunctionslow, &ctx)); 174 PetscCall(TSRHSSplitSetRHSFunction(ts, "fast", NULL, (TSRHSFunctionFn *)RHSFunctionfast, &ctx)); 175 PetscCall(TSRHSSplitSetIFunction(ts, "fast", NULL, (TSIFunctionFn *)IFunctionfast, &ctx)); 176 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 177 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 1, 1)); 178 PetscCall(MatSetFromOptions(A)); 179 PetscCall(MatSetUp(A)); 180 PetscCall(TSRHSSplitSetIJacobian(ts, "fast", A, A, (TSIJacobianFn *)IJacobianfast, &ctx)); 181 182 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 183 Set initial conditions 184 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 185 PetscCall(TSSetSolution(ts, U)); 186 187 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 188 Set solver options 189 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 190 PetscCall(TSSetMaxTime(ts, ctx.Tf)); 191 PetscCall(TSSetTimeStep(ts, ctx.dt)); 192 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 193 PetscCall(TSSetFromOptions(ts)); 194 195 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 196 Solve linear system 197 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 198 PetscCall(TSSolve(ts, U)); 199 PetscCall(VecView(U, PETSC_VIEWER_STDOUT_WORLD)); 200 201 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 202 Check the error of the Petsc solution 203 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 204 PetscCall(TSGetTime(ts, &tt)); 205 PetscCall(sol_true(tt, Utrue)); 206 PetscCall(VecAXPY(Utrue, -1.0, U)); 207 PetscCall(VecNorm(Utrue, NORM_2, &error)); 208 209 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 210 Print norm2 error 211 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 212 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "l2 error norm: %g\n", (double)PetscAbsScalar(error))); 213 214 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 215 Free work space. All PETSc objects should be destroyed when they are no longer needed. 216 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 217 PetscCall(VecDestroy(&U)); 218 PetscCall(TSDestroy(&ts)); 219 PetscCall(VecDestroy(&Utrue)); 220 PetscCall(ISDestroy(&iss)); 221 PetscCall(ISDestroy(&isf)); 222 PetscCall(PetscFree(indicess)); 223 PetscCall(PetscFree(indicesf)); 224 PetscCall(MatDestroy(&A)); 225 PetscCall(PetscFinalize()); 226 return 0; 227 } 228 229 /*TEST 230 build: 231 requires: !complex 232 233 test: 234 235 test: 236 suffix: 2 237 args: -ts_arkimex_initial_guess_extrapolate 1 238 output_file: output/ex3fastslowsplit_1.out 239 240 TEST*/ 241