1 static char help[] = "Constant acceleration check with 2nd-order generalized-alpha.\n";
2
3 #include <petscts.h>
4
5 typedef struct {
6 PetscReal a0; /* constant acceleration */
7 PetscReal u0, v0; /* initial conditions */
8 PetscReal radius; /* spectral radius of integrator */
9 } UserParams;
10
Exact(PetscReal t,PetscReal a0,PetscReal u0,PetscReal v0,PetscScalar * ut,PetscScalar * vt)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 }
16
Residual(TS ts,PetscReal t,Vec U,Vec V,Vec A,Vec R,PetscCtx ctx)17 PetscErrorCode Residual(TS ts, PetscReal t, Vec U, Vec V, Vec A, Vec R, PetscCtx ctx)
18 {
19 UserParams *user = (UserParams *)ctx;
20 const PetscScalar *a;
21 PetscScalar *r;
22
23 PetscFunctionBeginUser;
24 PetscCall(VecGetArrayRead(A, &a));
25 PetscCall(VecGetArrayWrite(R, &r));
26
27 r[0] = a[0] - user->a0;
28
29 PetscCall(VecRestoreArrayRead(A, &a));
30 PetscCall(VecRestoreArrayWrite(R, &r));
31 PetscFunctionReturn(PETSC_SUCCESS);
32 }
33
Tangent(TS ts,PetscReal t,Vec U,Vec V,Vec A,PetscReal shiftV,PetscReal shiftA,Mat J,Mat P,PetscCtx ctx)34 PetscErrorCode Tangent(TS ts, PetscReal t, Vec U, Vec V, Vec A, PetscReal shiftV, PetscReal shiftA, Mat J, Mat P, PetscCtx ctx)
35 {
36 PetscReal T = 0;
37
38 PetscFunctionBeginUser;
39 T = shiftA;
40
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 }
50
main(int argc,char * argv[])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};
64
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");
69
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();
76
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));
83
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));
92
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));
101
102 PetscCall(TS2SetSolution(ts, U, V));
103 PetscCall(TSSetFromOptions(ts));
104 PetscCall(TSSolve(ts, NULL));
105
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));
117
118 PetscCall(VecDestroy(&U));
119 PetscCall(VecDestroy(&V));
120 PetscCall(TSDestroy(&ts));
121 PetscCall(PetscFinalize());
122 return 0;
123 }
124
125 /*TEST
126
127 test:
128 suffix: a
129 args: -radius 0.0
130 requires: !single
131
132 test:
133 suffix: b
134 args: -radius 0.5
135 requires: !single
136
137 test:
138 suffix: c
139 args: -radius 1.0
140 requires: !single
141
142 TEST*/
143