1 static char help[] = "Single-DOF oscillator formulated as a second-order system.\n"; 2 3 #include <petscts.h> 4 5 typedef struct { 6 PetscReal Omega; /* natural frequency */ 7 PetscReal Xi; /* damping coefficient */ 8 PetscReal u0, v0; /* initial conditions */ 9 } UserParams; 10 11 static void Exact(PetscReal t, PetscReal omega, PetscReal xi, PetscReal u0, PetscReal v0, PetscReal *ut, PetscReal *vt) 12 { 13 PetscReal u, v; 14 if (xi < 1) { 15 PetscReal a = xi * omega; 16 PetscReal w = PetscSqrtReal(1 - xi * xi) * omega; 17 PetscReal C1 = (v0 + a * u0) / w; 18 PetscReal C2 = u0; 19 u = PetscExpReal(-a * t) * (C1 * PetscSinReal(w * t) + C2 * PetscCosReal(w * t)); 20 v = (-a * PetscExpReal(-a * t) * (C1 * PetscSinReal(w * t) + C2 * PetscCosReal(w * t)) + w * PetscExpReal(-a * t) * (C1 * PetscCosReal(w * t) - C2 * PetscSinReal(w * t))); 21 } else if (xi > 1) { 22 PetscReal w = PetscSqrtReal(xi * xi - 1) * omega; 23 PetscReal C1 = (w * u0 + xi * u0 + v0) / (2 * w); 24 PetscReal C2 = (w * u0 - xi * u0 - v0) / (2 * w); 25 u = C1 * PetscExpReal((-xi + w) * t) + C2 * PetscExpReal((-xi - w) * t); 26 v = C1 * (-xi + w) * PetscExpReal((-xi + w) * t) + C2 * (-xi - w) * PetscExpReal((-xi - w) * t); 27 } else { 28 PetscReal a = xi * omega; 29 PetscReal C1 = v0 + a * u0; 30 PetscReal C2 = u0; 31 u = (C1 * t + C2) * PetscExpReal(-a * t); 32 v = (C1 - a * (C1 * t + C2)) * PetscExpReal(-a * t); 33 } 34 if (ut) *ut = u; 35 if (vt) *vt = v; 36 } 37 38 PetscErrorCode Solution(TS ts, PetscReal t, Vec X, void *ctx) 39 { 40 UserParams *user = (UserParams *)ctx; 41 PetscReal u, v; 42 PetscScalar *x; 43 44 PetscFunctionBeginUser; 45 Exact(t, user->Omega, user->Xi, user->u0, user->v0, &u, &v); 46 PetscCall(VecGetArray(X, &x)); 47 x[0] = (PetscScalar)u; 48 PetscCall(VecRestoreArray(X, &x)); 49 PetscFunctionReturn(PETSC_SUCCESS); 50 } 51 52 PetscErrorCode Residual1(TS ts, PetscReal t, Vec U, Vec A, Vec R, void *ctx) 53 { 54 UserParams *user = (UserParams *)ctx; 55 PetscReal Omega = user->Omega; 56 const PetscScalar *u, *a; 57 PetscScalar *r; 58 59 PetscFunctionBeginUser; 60 PetscCall(VecGetArrayRead(U, &u)); 61 PetscCall(VecGetArrayRead(A, &a)); 62 PetscCall(VecGetArrayWrite(R, &r)); 63 64 r[0] = a[0] + (Omega * Omega) * u[0]; 65 66 PetscCall(VecRestoreArrayRead(U, &u)); 67 PetscCall(VecRestoreArrayRead(A, &a)); 68 PetscCall(VecRestoreArrayWrite(R, &r)); 69 PetscCall(VecAssemblyBegin(R)); 70 PetscCall(VecAssemblyEnd(R)); 71 PetscFunctionReturn(PETSC_SUCCESS); 72 } 73 74 PetscErrorCode Tangent1(TS ts, PetscReal t, Vec U, Vec A, PetscReal shiftA, Mat J, Mat P, void *ctx) 75 { 76 UserParams *user = (UserParams *)ctx; 77 PetscReal Omega = user->Omega; 78 PetscReal T = 0; 79 80 PetscFunctionBeginUser; 81 82 T = shiftA + (Omega * Omega); 83 84 PetscCall(MatSetValue(P, 0, 0, T, INSERT_VALUES)); 85 PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY)); 86 PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY)); 87 if (J != P) { 88 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 89 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 90 } 91 PetscFunctionReturn(PETSC_SUCCESS); 92 } 93 94 PetscErrorCode Residual2(TS ts, PetscReal t, Vec U, Vec V, Vec A, Vec R, void *ctx) 95 { 96 UserParams *user = (UserParams *)ctx; 97 PetscReal Omega = user->Omega, Xi = user->Xi; 98 const PetscScalar *u, *v, *a; 99 PetscScalar *r; 100 101 PetscFunctionBeginUser; 102 PetscCall(VecGetArrayRead(U, &u)); 103 PetscCall(VecGetArrayRead(V, &v)); 104 PetscCall(VecGetArrayRead(A, &a)); 105 PetscCall(VecGetArrayWrite(R, &r)); 106 107 r[0] = a[0] + (2 * Xi * Omega) * v[0] + (Omega * Omega) * u[0]; 108 109 PetscCall(VecRestoreArrayRead(U, &u)); 110 PetscCall(VecRestoreArrayRead(V, &v)); 111 PetscCall(VecRestoreArrayRead(A, &a)); 112 PetscCall(VecRestoreArrayWrite(R, &r)); 113 PetscCall(VecAssemblyBegin(R)); 114 PetscCall(VecAssemblyEnd(R)); 115 PetscFunctionReturn(PETSC_SUCCESS); 116 } 117 118 PetscErrorCode Tangent2(TS ts, PetscReal t, Vec U, Vec V, Vec A, PetscReal shiftV, PetscReal shiftA, Mat J, Mat P, void *ctx) 119 { 120 UserParams *user = (UserParams *)ctx; 121 PetscReal Omega = user->Omega, Xi = user->Xi; 122 PetscReal T = 0; 123 124 PetscFunctionBeginUser; 125 126 T = shiftA + shiftV * (2 * Xi * Omega) + (Omega * Omega); 127 128 PetscCall(MatSetValue(P, 0, 0, T, INSERT_VALUES)); 129 PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY)); 130 PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY)); 131 if (J != P) { 132 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 133 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 134 } 135 PetscFunctionReturn(PETSC_SUCCESS); 136 } 137 138 int main(int argc, char *argv[]) 139 { 140 PetscMPIInt size; 141 TS ts; 142 Vec R; 143 Mat J; 144 Vec U, V; 145 PetscScalar *u, *v; 146 UserParams user = {/*Omega=*/1, /*Xi=*/0, /*u0=*/1, /*,v0=*/0}; 147 148 PetscFunctionBeginUser; 149 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 150 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 151 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs"); 152 153 PetscOptionsBegin(PETSC_COMM_SELF, "", "ex43 options", ""); 154 PetscCall(PetscOptionsReal("-frequency", "Natual frequency", __FILE__, user.Omega, &user.Omega, NULL)); 155 PetscCall(PetscOptionsReal("-damping", "Damping coefficient", __FILE__, user.Xi, &user.Xi, NULL)); 156 PetscCall(PetscOptionsReal("-initial_u", "Initial displacement", __FILE__, user.u0, &user.u0, NULL)); 157 PetscCall(PetscOptionsReal("-initial_v", "Initial velocity", __FILE__, user.v0, &user.v0, NULL)); 158 PetscOptionsEnd(); 159 160 PetscCall(TSCreate(PETSC_COMM_SELF, &ts)); 161 PetscCall(TSSetType(ts, TSALPHA2)); 162 PetscCall(TSSetMaxTime(ts, 5 * (2 * PETSC_PI))); 163 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 164 PetscCall(TSSetTimeStep(ts, 0.01)); 165 166 PetscCall(VecCreateSeq(PETSC_COMM_SELF, 1, &R)); 167 PetscCall(VecSetUp(R)); 168 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, 1, 1, NULL, &J)); 169 PetscCall(MatSetUp(J)); 170 if (user.Xi) { 171 PetscCall(TSSetI2Function(ts, R, Residual2, &user)); 172 PetscCall(TSSetI2Jacobian(ts, J, J, Tangent2, &user)); 173 } else { 174 PetscCall(TSSetIFunction(ts, R, Residual1, &user)); 175 PetscCall(TSSetIJacobian(ts, J, J, Tangent1, &user)); 176 } 177 PetscCall(VecDestroy(&R)); 178 PetscCall(MatDestroy(&J)); 179 PetscCall(TSSetSolutionFunction(ts, Solution, &user)); 180 181 PetscCall(VecCreateSeq(PETSC_COMM_SELF, 1, &U)); 182 PetscCall(VecCreateSeq(PETSC_COMM_SELF, 1, &V)); 183 PetscCall(VecGetArrayWrite(U, &u)); 184 PetscCall(VecGetArrayWrite(V, &v)); 185 u[0] = user.u0; 186 v[0] = user.v0; 187 PetscCall(VecRestoreArrayWrite(U, &u)); 188 PetscCall(VecRestoreArrayWrite(V, &v)); 189 190 PetscCall(TS2SetSolution(ts, U, V)); 191 PetscCall(TSSetFromOptions(ts)); 192 PetscCall(TSSolve(ts, NULL)); 193 194 PetscCall(VecDestroy(&U)); 195 PetscCall(VecDestroy(&V)); 196 PetscCall(TSDestroy(&ts)); 197 PetscCall(PetscFinalize()); 198 return 0; 199 } 200 201 /*TEST 202 203 test: 204 suffix: a 205 args: -ts_max_steps 10 -ts_view 206 requires: !single 207 208 test: 209 suffix: b 210 args: -ts_max_steps 10 -ts_rtol 0 -ts_atol 1e-5 -ts_adapt_type basic -ts_adapt_monitor 211 requires: !single 212 213 TEST*/ 214