Actual source code: ex17.c

  1: static char help[] = "Solves the ODE du/dt = poly(t), u(0) = 0. Tests TSResize for varying size.\n\n";

  3: #include <petscts.h>

  5: PetscScalar poly(PetscInt p, PetscReal t)
  6: {
  7:   return p ? t * poly(p - 1, t) : 1.0;
  8: }

 10: PetscScalar dpoly(PetscInt p, PetscReal t)
 11: {
 12:   return p > 0 ? (PetscReal)p * poly(p - 1, t) : 0.0;
 13: }

 15: PetscErrorCode CreateVec(PetscInt lsize, Vec *out)
 16: {
 17:   Vec x;

 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: }

 28: PetscErrorCode CreateMat(PetscInt lsize, Mat *out)
 29: {
 30:   Mat A;

 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: }

 46: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec x, Vec f, void *ctx)
 47: {
 48:   PetscInt *order = (PetscInt *)ctx;

 50:   PetscFunctionBeginUser;
 51:   PetscCall(VecSet(f, dpoly(*order, t)));
 52:   PetscFunctionReturn(PETSC_SUCCESS);
 53: }

 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: }

 63: PetscErrorCode Transfer(TS ts, PetscInt nv, Vec vecsin[], Vec vecsout[], void *ctx)
 64: {
 65:   PetscInt n, nnew;

 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;

 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: }

 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: }

 93: PetscErrorCode Monitor(TS ts, PetscInt n, PetscReal t, Vec x, void *ctx)
 94: {
 95:   const PetscScalar *a;
 96:   PetscScalar       *store = (PetscScalar *)ctx;

 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: }

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  errors with arkimex with 1.e-18 errors, macOS clang requires an even larger tolerance */
113:   PetscReal tol = 10 * PETSC_MACHINE_EPSILON;

115:   PetscFunctionBeginUser;
116:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
117:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-order", &order, NULL));
118:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-tol", &tol, NULL));

120:   for (PetscInt i = 0; i < 2; i++) {
121:     PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
122:     PetscCall(TSSetProblemType(ts, TS_LINEAR));

124:     PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, &order));

126:     PetscCall(CreateMat(1, &A));
127:     PetscCall(TSSetRHSJacobian(ts, A, A, RHSJacobian, NULL));
128:     PetscCall(MatDestroy(&A));

130:     PetscCall(CreateVec(1, &x));
131:     PetscCall(TSSetSolution(ts, x));
132:     PetscCall(VecDestroy(&x));

134:     for (PetscInt j = 0; j < 10; j++) results[i][j] = 0;
135:     PetscCall(TSMonitorSet(ts, Monitor, results[i], NULL));
136:     PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
137:     if (i) PetscCall(TSSetResize(ts, TransferSetUp, Transfer, NULL));
138:     PetscCall(TSSetTime(ts, 0));
139:     PetscCall(TSSetTimeStep(ts, 1. / 4.));
140:     PetscCall(TSSetMaxSteps(ts, 10));
141:     PetscCall(TSSetFromOptions(ts));

143:     PetscCall(TSSolve(ts, NULL));

145:     PetscCall(TSDestroy(&ts));
146:   }

148:   /* Dump errors if any */
149:   PetscBool flg = PETSC_FALSE;
150:   for (PetscInt i = 0; i < 10; i++) {
151:     PetscReal err = PetscAbsScalar(results[0][i] - results[1][i]);
152:     if (err > tol) {
153:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error step %" PetscInt_FMT ": %g\n", i, (double)err));
154:       flg = PETSC_TRUE;
155:     }
156:   }
157:   if (flg) {
158:     PetscCall(PetscScalarView(10, results[0], PETSC_VIEWER_STDOUT_WORLD));
159:     PetscCall(PetscScalarView(10, results[1], PETSC_VIEWER_STDOUT_WORLD));
160:   }
161:   PetscCall(PetscFinalize());
162:   return 0;
163: }

165: /*TEST

167:     test:
168:       suffix: bdf
169:       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
170:       output_file: output/ex17.out

172:     test:
173:       suffix: expl
174:       args: -ts_adapt_wnormtype infinity -ts_type {{euler rk ssp}} -order 6 -ts_adapt_type {{none basic dsp}}
175:       output_file: output/ex17.out

177:     test:
178:       suffix: impl
179:       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
180:       output_file: output/ex17.out

182: TEST*/