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 11 static void Exact(PetscReal t, PetscReal v0, PetscReal u0, PetscScalar *ut) 12 { 13 if (ut) *ut = u0 + v0 * t; 14 } 15 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; 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 33 PetscErrorCode Tangent(TS ts, PetscReal t, Vec U, Vec V, PetscReal shiftV, Mat J, Mat P, void *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 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