1 static char help[] = "Constant velocity check with 1st-order generalized-alpha.\n";
2
3 #include <petscts.h>
4
5 typedef struct {
6 PetscReal v0; /* constant velocity */
7 PetscReal u0; /* initial condition */
8 PetscReal radius; /* spectral radius of integrator */
9 } UserParams;
10
Exact(PetscReal t,PetscReal v0,PetscReal u0,PetscScalar * ut)11 static void Exact(PetscReal t, PetscReal v0, PetscReal u0, PetscScalar *ut)
12 {
13 if (ut) *ut = u0 + v0 * t;
14 }
15
Residual(TS ts,PetscReal t,Vec U,Vec V,Vec R,PetscCtx ctx)16 PetscErrorCode Residual(TS ts, PetscReal t, Vec U, Vec V, Vec R, PetscCtx ctx)
17 {
18 UserParams *user = (UserParams *)ctx;
19 const PetscScalar *v;
20 PetscScalar *r;
21
22 PetscFunctionBeginUser;
23 PetscCall(VecGetArrayRead(V, &v));
24 PetscCall(VecGetArrayWrite(R, &r));
25
26 r[0] = v[0] - user->v0;
27
28 PetscCall(VecRestoreArrayRead(V, &v));
29 PetscCall(VecRestoreArrayWrite(R, &r));
30 PetscFunctionReturn(PETSC_SUCCESS);
31 }
32
Tangent(TS ts,PetscReal t,Vec U,Vec V,PetscReal shiftV,Mat J,Mat P,PetscCtx ctx)33 PetscErrorCode Tangent(TS ts, PetscReal t, Vec U, Vec V, PetscReal shiftV, Mat J, Mat P, PetscCtx ctx)
34 {
35 PetscReal T = 0;
36
37 PetscFunctionBeginUser;
38 T = shiftV;
39
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 }
49
main(int argc,char * argv[])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};
63
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");
68
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();
74
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));
81
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));
90
91 PetscCall(VecCreateSeq(PETSC_COMM_SELF, 1, &U));
92 PetscCall(VecGetArrayWrite(U, &u));
93 u[0] = user.u0;
94 PetscCall(VecRestoreArrayWrite(U, &u));
95
96 PetscCall(TSSetSolution(ts, U));
97 PetscCall(TSSetFromOptions(ts));
98 PetscCall(TSSolve(ts, NULL));
99
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));
106
107 PetscCall(VecDestroy(&U));
108 PetscCall(TSDestroy(&ts));
109 PetscCall(PetscFinalize());
110 return 0;
111 }
112
113 /*TEST
114
115 test:
116 suffix: a
117 args: -radius 0.0
118 requires: !single
119
120 test:
121 suffix: b
122 args: -radius 0.5
123 requires: !single
124
125 test:
126 suffix: c
127 args: -radius 1.0
128 requires: !single
129
130 TEST*/
131