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