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