1 static char help[] = "Solves the ODE du/dt = poly(t), u(0) = 0. Tests TSResize for varying size.\n\n"; 2 3 #include <petscts.h> 4 5 PetscScalar poly(PetscInt p, PetscReal t) 6 { 7 return p ? t * poly(p - 1, t) : 1.0; 8 } 9 10 PetscScalar dpoly(PetscInt p, PetscReal t) 11 { 12 return p > 0 ? (PetscReal)p * poly(p - 1, t) : 0.0; 13 } 14 15 PetscErrorCode CreateVec(PetscInt lsize, Vec *out) 16 { 17 Vec x; 18 19 PetscFunctionBeginUser; 20 PetscCall(VecCreate(PETSC_COMM_WORLD, &x)); 21 PetscCall(VecSetSizes(x, lsize, PETSC_DECIDE)); 22 PetscCall(VecSetFromOptions(x)); 23 PetscCall(VecSetUp(x)); 24 *out = x; 25 PetscFunctionReturn(PETSC_SUCCESS); 26 } 27 28 PetscErrorCode CreateMat(PetscInt lsize, Mat *out) 29 { 30 Mat A; 31 32 PetscFunctionBeginUser; 33 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 34 PetscCall(MatSetSizes(A, lsize, lsize, PETSC_DECIDE, PETSC_DECIDE)); 35 PetscCall(MatSetFromOptions(A)); 36 PetscCall(MatSetUp(A)); 37 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 38 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 39 /* ensure that the Jacobian matrix has diagonal entries since that is required by TS */ 40 PetscCall(MatShift(A, (PetscReal)1)); 41 PetscCall(MatShift(A, (PetscReal)-1)); 42 *out = A; 43 PetscFunctionReturn(PETSC_SUCCESS); 44 } 45 46 PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec x, Vec f, void *ctx) 47 { 48 PetscInt *order = (PetscInt *)ctx; 49 50 PetscFunctionBeginUser; 51 PetscCall(VecSet(f, dpoly(*order, t))); 52 PetscFunctionReturn(PETSC_SUCCESS); 53 } 54 55 PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec x, Mat A, Mat B, void *ctx) 56 { 57 PetscFunctionBeginUser; 58 PetscCall(MatZeroEntries(B)); 59 if (B != A) PetscCall(MatZeroEntries(A)); 60 PetscFunctionReturn(PETSC_SUCCESS); 61 } 62 63 PetscErrorCode Transfer(TS ts, PetscInt nv, Vec vecsin[], Vec vecsout[], void *ctx) 64 { 65 PetscInt n, nnew; 66 67 PetscFunctionBeginUser; 68 PetscAssert(nv > 0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Zero vectors"); 69 PetscCall(VecGetLocalSize(vecsin[0], &n)); 70 nnew = n == 2 ? 1 : 2; 71 for (PetscInt i = 0; i < nv; i++) { 72 const PetscScalar *vals; 73 74 PetscCall(CreateVec(nnew, &vecsout[i])); 75 PetscCall(VecGetArrayRead(vecsin[i], &vals)); 76 PetscCall(VecSet(vecsout[i], vals[0])); 77 PetscCall(VecRestoreArrayRead(vecsin[i], &vals)); 78 } 79 Mat A; 80 PetscCall(CreateMat(nnew, &A)); 81 PetscCall(TSSetRHSJacobian(ts, A, A, RHSJacobian, NULL)); 82 PetscCall(MatDestroy(&A)); 83 PetscFunctionReturn(PETSC_SUCCESS); 84 } 85 86 PetscErrorCode TransferSetUp(TS ts, PetscInt step, PetscReal time, Vec sol, PetscBool *resize, void *ctx) 87 { 88 PetscBool *alreadydone = (PetscBool *)ctx; 89 90 PetscFunctionBeginUser; 91 *alreadydone = (PetscBool) !(*alreadydone); 92 *resize = *alreadydone; 93 PetscFunctionReturn(PETSC_SUCCESS); 94 } 95 96 PetscErrorCode Monitor(TS ts, PetscInt n, PetscReal t, Vec x, void *ctx) 97 { 98 const PetscScalar *a; 99 PetscScalar *store = (PetscScalar *)ctx; 100 101 PetscFunctionBeginUser; 102 PetscCall(VecGetArrayRead(x, &a)); 103 if (n < 10) store[n] = a[0]; 104 PetscCall(VecRestoreArrayRead(x, &a)); 105 PetscFunctionReturn(PETSC_SUCCESS); 106 } 107 108 int main(int argc, char **argv) 109 { 110 TS ts; 111 Vec x; 112 Mat A; 113 PetscInt order = 2; 114 PetscScalar results[3][10]; 115 /* I would like to use 0 here, but arkimex errors with 1.e-14 discrepancy when using TSRResize without restart on some machines (mostly arm-osx) */ 116 PetscReal tol = PETSC_SMALL; 117 118 PetscFunctionBeginUser; 119 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 120 PetscCall(PetscOptionsGetInt(NULL, NULL, "-order", &order, NULL)); 121 PetscCall(PetscOptionsGetReal(NULL, NULL, "-tol", &tol, NULL)); 122 123 for (PetscInt i = 0; i < 3; i++) { 124 PetscBool alreadydone = PETSC_TRUE; 125 126 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 127 PetscCall(TSSetProblemType(ts, TS_LINEAR)); 128 129 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, &order)); 130 131 PetscCall(CreateMat(1, &A)); 132 PetscCall(TSSetRHSJacobian(ts, A, A, RHSJacobian, NULL)); 133 PetscCall(MatDestroy(&A)); 134 135 PetscCall(CreateVec(1, &x)); 136 PetscCall(TSSetSolution(ts, x)); 137 PetscCall(VecDestroy(&x)); 138 139 for (PetscInt j = 0; j < 10; j++) results[i][j] = 0; 140 PetscCall(TSMonitorSet(ts, Monitor, results[i], NULL)); 141 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 142 if (i) PetscCall(TSSetResize(ts, i == 1 ? PETSC_TRUE : PETSC_FALSE, TransferSetUp, Transfer, &alreadydone)); 143 PetscCall(TSSetTime(ts, 0)); 144 PetscCall(TSSetTimeStep(ts, 1. / 4.)); 145 PetscCall(TSSetMaxSteps(ts, 10)); 146 PetscCall(TSSetFromOptions(ts)); 147 148 PetscCall(TSSolve(ts, NULL)); 149 150 PetscCall(TSDestroy(&ts)); 151 } 152 153 /* Dump errors if any */ 154 PetscBool flg = PETSC_FALSE; 155 for (PetscInt i = 0; i < 10; i++) { 156 PetscReal err; 157 158 err = PetscAbsScalar(results[0][i] - results[1][i]); 159 if (err > tol) { 160 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error with restart for step %" PetscInt_FMT ": %g\n", i, (double)err)); 161 flg = PETSC_TRUE; 162 } 163 err = PetscAbsScalar(results[0][i] - results[2][i]); 164 if (err > tol) { 165 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error without restart for step %" PetscInt_FMT ": %g\n", i, (double)err)); 166 flg = PETSC_TRUE; 167 } 168 } 169 if (flg) { 170 PetscCall(PetscScalarView(10, results[0], PETSC_VIEWER_STDOUT_WORLD)); 171 PetscCall(PetscScalarView(10, results[1], PETSC_VIEWER_STDOUT_WORLD)); 172 PetscCall(PetscScalarView(10, results[2], PETSC_VIEWER_STDOUT_WORLD)); 173 } 174 PetscCall(PetscFinalize()); 175 return 0; 176 } 177 178 /*TEST 179 testset: 180 # using gemv gives larger error which failed error checking 181 args: -vec_mdot_use_gemv 0 -vec_maxpy_use_gemv 0 182 183 test: 184 suffix: bdf 185 args: -ts_adapt_wnormtype infinity -ts_type bdf -ts_bdf_order {{2 3 4 5 6}} -order 6 -ts_adapt_type {{none basic dsp}} -ksp_type preonly -pc_type lu 186 output_file: output/ex17.out 187 188 test: 189 suffix: expl 190 args: -ts_adapt_wnormtype infinity -ts_type {{euler rk ssp}} -order 6 -ts_adapt_type {{none basic dsp}} 191 output_file: output/ex17.out 192 193 test: 194 suffix: impl 195 args: -ts_adapt_wnormtype infinity -ts_type {{rosw beuler cn alpha theta arkimex}} -order 6 -ts_adapt_type {{none basic dsp}} -ksp_type preonly -pc_type lu 196 output_file: output/ex17.out 197 198 TEST*/ 199