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

  3: #include <petscts.h>

  5: typedef struct {
  6:   PetscReal v0;     /* constant velocity */
  7:   PetscReal u0;     /* initial condition */
  8:   PetscReal radius; /* spectral radius of integrator */
  9: } UserParams;

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

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

 22:   PetscFunctionBeginUser;
 23:   PetscCall(VecGetArrayRead(V, &v));
 24:   PetscCall(VecGetArrayWrite(R, &r));

 26:   r[0] = v[0] - user->v0;

 28:   PetscCall(VecRestoreArrayRead(V, &v));
 29:   PetscCall(VecRestoreArrayWrite(R, &r));
 30:   PetscFunctionReturn(PETSC_SUCCESS);
 31: }

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

 37:   PetscFunctionBeginUser;
 38:   T = shiftV;

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

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

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

 69:   PetscOptionsBegin(PETSC_COMM_SELF, "", "ex81 options", "");
 70:   PetscCall(PetscOptionsReal("-velocity", "Constant velocity", __FILE__, user.v0, &user.v0, NULL));
 71:   PetscCall(PetscOptionsReal("-initial_u", "Initial displacement", __FILE__, user.u0, &user.u0, NULL));
 72:   PetscCall(PetscOptionsReal("-radius", "Spectral radius", __FILE__, user.radius, &user.radius, NULL));
 73:   PetscOptionsEnd();

 75:   PetscCall(TSCreate(PETSC_COMM_SELF, &ts));
 76:   PetscCall(TSSetType(ts, TSALPHA));
 77:   PetscCall(TSSetMaxTime(ts, t_final));
 78:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
 79:   PetscCall(TSSetTimeStep(ts, t_final / n_step));
 80:   PetscCall(TSAlphaSetRadius(ts, user.radius));

 82:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, 1, &R));
 83:   PetscCall(VecSetUp(R));
 84:   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, 1, 1, NULL, &J));
 85:   PetscCall(MatSetUp(J));
 86:   PetscCall(TSSetIFunction(ts, R, Residual, &user));
 87:   PetscCall(TSSetIJacobian(ts, J, J, Tangent, &user));
 88:   PetscCall(VecDestroy(&R));
 89:   PetscCall(MatDestroy(&J));

 91:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, 1, &U));
 92:   PetscCall(VecGetArrayWrite(U, &u));
 93:   u[0] = user.u0;
 94:   PetscCall(VecRestoreArrayWrite(U, &u));

 96:   PetscCall(TSSetSolution(ts, U));
 97:   PetscCall(TSSetFromOptions(ts));
 98:   PetscCall(TSSolve(ts, NULL));

100:   PetscCall(VecGetArray(U, &u));
101:   Exact(t_final, user.v0, user.u0, &u_exact);
102:   u_err = PetscAbsScalar(u[0] - u_exact);
103:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "u(t=%g) = %g (error = %g)\n", (double)t_final, (double)PetscRealPart(u[0]), (double)u_err));
104:   PetscCheck(u_err < atol, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Inexact displacement.");
105:   PetscCall(VecRestoreArray(U, &u));

107:   PetscCall(VecDestroy(&U));
108:   PetscCall(TSDestroy(&ts));
109:   PetscCall(PetscFinalize());
110:   return 0;
111: }

113: /*TEST

115:     test:
116:       suffix: a
117:       args: -radius 0.0
118:       requires: !single

120:     test:
121:       suffix: b
122:       args: -radius 0.5
123:       requires: !single

125:     test:
126:       suffix: c
127:       args: -radius 1.0
128:       requires: !single

130: TEST*/