1: static char help[] = "Constant acceleration check with 2nd-order generalized-alpha.\n";

  3: #include <petscts.h>

  5: typedef struct {
  6:   PetscReal a0;     /* constant acceleration  */
  7:   PetscReal u0, v0; /* initial conditions */
  8:   PetscReal radius; /* spectral radius of integrator */
  9: } UserParams;

 11: static void Exact(PetscReal t, PetscReal a0, PetscReal u0, PetscReal v0, PetscScalar *ut, PetscScalar *vt)
 12: {
 13:   if (ut) *ut = u0 + v0 * t + 0.5 * a0 * t * t;
 14:   if (vt) *vt = v0 + a0 * t;
 15: }

 17: PetscErrorCode Residual(TS ts, PetscReal t, Vec U, Vec V, Vec A, Vec R, void *ctx)
 18: {
 19:   UserParams        *user = (UserParams *)ctx;
 20:   const PetscScalar *a;
 21:   PetscScalar       *r;

 23:   PetscFunctionBeginUser;
 24:   PetscCall(VecGetArrayRead(A, &a));
 25:   PetscCall(VecGetArrayWrite(R, &r));

 27:   r[0] = a[0] - user->a0;

 29:   PetscCall(VecRestoreArrayRead(A, &a));
 30:   PetscCall(VecRestoreArrayWrite(R, &r));
 31:   PetscFunctionReturn(PETSC_SUCCESS);
 32: }

 34: PetscErrorCode Tangent(TS ts, PetscReal t, Vec U, Vec V, Vec A, PetscReal shiftV, PetscReal shiftA, Mat J, Mat P, void *ctx)
 35: {
 36:   PetscReal T = 0;

 38:   PetscFunctionBeginUser;
 39:   T = shiftA;

 41:   PetscCall(MatSetValue(P, 0, 0, T, INSERT_VALUES));
 42:   PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
 43:   PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
 44:   if (J != P) {
 45:     PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
 46:     PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
 47:   }
 48:   PetscFunctionReturn(PETSC_SUCCESS);
 49: }

 51: int main(int argc, char *argv[])
 52: {
 53:   PetscMPIInt  size;
 54:   TS           ts;
 55:   Vec          R;
 56:   Mat          J;
 57:   Vec          U, V;
 58:   PetscScalar *u, *v, u_exact, v_exact;
 59:   PetscReal    u_err, v_err;
 60:   PetscReal    atol    = 1e-15;
 61:   PetscReal    t_final = 3.0;
 62:   PetscInt     n_step  = 8;
 63:   UserParams   user    = {/*a0=*/1, /*u0=*/1, /*v0=*/0, /*radius=*/0.0};

 65:   PetscFunctionBeginUser;
 66:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 67:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 68:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");

 70:   PetscOptionsBegin(PETSC_COMM_SELF, "", "ex80 options", "");
 71:   PetscCall(PetscOptionsReal("-acceleration", "Constant acceleration", __FILE__, user.a0, &user.a0, NULL));
 72:   PetscCall(PetscOptionsReal("-initial_u", "Initial displacement", __FILE__, user.u0, &user.u0, NULL));
 73:   PetscCall(PetscOptionsReal("-initial_v", "Initial velocity", __FILE__, user.v0, &user.v0, NULL));
 74:   PetscCall(PetscOptionsReal("-radius", "Spectral radius", __FILE__, user.radius, &user.radius, NULL));
 75:   PetscOptionsEnd();

 77:   PetscCall(TSCreate(PETSC_COMM_SELF, &ts));
 78:   PetscCall(TSSetType(ts, TSALPHA2));
 79:   PetscCall(TSSetMaxTime(ts, t_final));
 80:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
 81:   PetscCall(TSSetTimeStep(ts, t_final / n_step));
 82:   PetscCall(TSAlpha2SetRadius(ts, user.radius));

 84:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, 1, &R));
 85:   PetscCall(VecSetUp(R));
 86:   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, 1, 1, NULL, &J));
 87:   PetscCall(MatSetUp(J));
 88:   PetscCall(TSSetI2Function(ts, R, Residual, &user));
 89:   PetscCall(TSSetI2Jacobian(ts, J, J, Tangent, &user));
 90:   PetscCall(VecDestroy(&R));
 91:   PetscCall(MatDestroy(&J));

 93:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, 1, &U));
 94:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, 1, &V));
 95:   PetscCall(VecGetArrayWrite(U, &u));
 96:   PetscCall(VecGetArrayWrite(V, &v));
 97:   u[0] = user.u0;
 98:   v[0] = user.v0;
 99:   PetscCall(VecRestoreArrayWrite(U, &u));
100:   PetscCall(VecRestoreArrayWrite(V, &v));

102:   PetscCall(TS2SetSolution(ts, U, V));
103:   PetscCall(TSSetFromOptions(ts));
104:   PetscCall(TSSolve(ts, NULL));

106:   PetscCall(VecGetArray(U, &u));
107:   PetscCall(VecGetArray(V, &v));
108:   Exact(t_final, user.a0, user.u0, user.v0, &u_exact, &v_exact);
109:   u_err = PetscAbsScalar(u[0] - u_exact);
110:   v_err = PetscAbsScalar(v[0] - v_exact);
111:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "u(t=%g) = %g (error = %g)\n", (double)t_final, (double)PetscRealPart(u[0]), (double)u_err));
112:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "v(t=%g) = %g (error = %g)\n", (double)t_final, (double)PetscRealPart(v[0]), (double)v_err));
113:   PetscCheck(u_err < atol, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Inexact displacement.");
114:   PetscCheck(v_err < atol, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Inexact velocity.");
115:   PetscCall(VecRestoreArray(U, &u));
116:   PetscCall(VecRestoreArray(V, &v));

118:   PetscCall(VecDestroy(&U));
119:   PetscCall(VecDestroy(&V));
120:   PetscCall(TSDestroy(&ts));
121:   PetscCall(PetscFinalize());
122:   return 0;
123: }

125: /*TEST

127:     test:
128:       suffix: a
129:       args: -radius 0.0
130:       requires: !single

132:     test:
133:       suffix: b
134:       args: -radius 0.5
135:       requires: !single

137:     test:
138:       suffix: c
139:       args: -radius 1.0
140:       requires: !single

142: TEST*/