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 PetscFunctionBeginUser; 89 *resize = PETSC_TRUE; 90 PetscFunctionReturn(PETSC_SUCCESS); 91 } 92 93 PetscErrorCode Monitor(TS ts, PetscInt n, PetscReal t, Vec x, void *ctx) 94 { 95 const PetscScalar *a; 96 PetscScalar *store = (PetscScalar *)ctx; 97 98 PetscFunctionBeginUser; 99 PetscCall(VecGetArrayRead(x, &a)); 100 if (n < 10) store[n] = a[0]; 101 PetscCall(VecRestoreArrayRead(x, &a)); 102 PetscFunctionReturn(PETSC_SUCCESS); 103 } 104 105 int main(int argc, char **argv) 106 { 107 TS ts; 108 Vec x; 109 Mat A; 110 PetscInt order = 2; 111 PetscScalar results[2][10]; 112 /* I would like to use 0 here, but linux-gcc-complex-opt-32bit 113 errors with arkimex with 1.e-18 errors */ 114 PetscReal tol = PETSC_MACHINE_EPSILON; 115 116 PetscFunctionBeginUser; 117 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 118 PetscCall(PetscOptionsGetInt(NULL, NULL, "-order", &order, NULL)); 119 PetscCall(PetscOptionsGetReal(NULL, NULL, "-tol", &tol, NULL)); 120 121 for (PetscInt i = 0; i < 2; i++) { 122 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 123 PetscCall(TSSetProblemType(ts, TS_LINEAR)); 124 125 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, &order)); 126 127 PetscCall(CreateMat(1, &A)); 128 PetscCall(TSSetRHSJacobian(ts, A, A, RHSJacobian, NULL)); 129 PetscCall(MatDestroy(&A)); 130 131 PetscCall(CreateVec(1, &x)); 132 PetscCall(TSSetSolution(ts, x)); 133 PetscCall(VecDestroy(&x)); 134 135 for (PetscInt j = 0; j < 10; j++) results[i][j] = 0; 136 PetscCall(TSMonitorSet(ts, Monitor, results[i], NULL)); 137 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 138 if (i) PetscCall(TSSetResize(ts, TransferSetUp, Transfer, NULL)); 139 PetscCall(TSSetTime(ts, 0)); 140 PetscCall(TSSetTimeStep(ts, 1. / 4.)); 141 PetscCall(TSSetMaxSteps(ts, 10)); 142 PetscCall(TSSetFromOptions(ts)); 143 144 PetscCall(TSSolve(ts, NULL)); 145 146 PetscCall(TSDestroy(&ts)); 147 } 148 149 /* Dump errors if any */ 150 PetscBool flg = PETSC_FALSE; 151 for (PetscInt i = 0; i < 10; i++) { 152 PetscReal err = PetscAbsScalar(results[0][i] - results[1][i]); 153 if (err > tol) { 154 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error step %" PetscInt_FMT ": %g\n", i, (double)err)); 155 flg = PETSC_TRUE; 156 } 157 } 158 if (flg) { 159 PetscCall(PetscScalarView(10, results[0], PETSC_VIEWER_STDOUT_WORLD)); 160 PetscCall(PetscScalarView(10, results[1], PETSC_VIEWER_STDOUT_WORLD)); 161 } 162 PetscCall(PetscFinalize()); 163 return 0; 164 } 165 166 /*TEST 167 168 test: 169 suffix: bdf 170 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 171 output_file: output/ex17.out 172 173 test: 174 suffix: expl 175 args: -ts_adapt_wnormtype infinity -ts_type {{euler rk ssp}} -order 6 -ts_adapt_type {{none basic dsp}} 176 output_file: output/ex17.out 177 178 test: 179 suffix: impl 180 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 181 output_file: output/ex17.out 182 183 TEST*/ 184