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 T = shiftA + (Omega * Omega); 97 98 PetscCall(MatSetValue(P, 0, 0, T, INSERT_VALUES)); 99 PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY)); 100 PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY)); 101 if (J != P) { 102 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 103 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 104 } 105 PetscFunctionReturn(PETSC_SUCCESS); 106 } 107 108 PetscErrorCode Residual2(TS ts, PetscReal t, Vec U, Vec V, Vec A, Vec R, void *ctx) 109 { 110 UserParams *user = (UserParams *)ctx; 111 PetscReal Omega = user->Omega, Xi = user->Xi; 112 const PetscScalar *u, *v, *a; 113 PetscScalar *r; 114 115 PetscFunctionBeginUser; 116 PetscCall(VecGetArrayRead(U, &u)); 117 PetscCall(VecGetArrayRead(V, &v)); 118 PetscCall(VecGetArrayRead(A, &a)); 119 PetscCall(VecGetArrayWrite(R, &r)); 120 121 r[0] = a[0] + (2 * Xi * Omega) * v[0] + (Omega * Omega) * u[0]; 122 123 PetscCall(VecRestoreArrayRead(U, &u)); 124 PetscCall(VecRestoreArrayRead(V, &v)); 125 PetscCall(VecRestoreArrayRead(A, &a)); 126 PetscCall(VecRestoreArrayWrite(R, &r)); 127 PetscCall(VecAssemblyBegin(R)); 128 PetscCall(VecAssemblyEnd(R)); 129 PetscFunctionReturn(PETSC_SUCCESS); 130 } 131 132 PetscErrorCode Tangent2(TS ts, PetscReal t, Vec U, Vec V, Vec A, PetscReal shiftV, PetscReal shiftA, Mat J, Mat P, void *ctx) 133 { 134 UserParams *user = (UserParams *)ctx; 135 PetscReal Omega = user->Omega, Xi = user->Xi; 136 PetscReal T = 0; 137 138 PetscFunctionBeginUser; 139 T = shiftA + shiftV * (2 * Xi * Omega) + (Omega * Omega); 140 141 PetscCall(MatSetValue(P, 0, 0, T, INSERT_VALUES)); 142 PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY)); 143 PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY)); 144 if (J != P) { 145 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 146 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 147 } 148 PetscFunctionReturn(PETSC_SUCCESS); 149 } 150 151 int main(int argc, char *argv[]) 152 { 153 PetscMPIInt size; 154 TS ts; 155 Vec R; 156 Mat J; 157 Vec U, V; 158 PetscScalar *u, *v; 159 UserParams user = {/*Omega=*/1, /*Xi=*/0, /*u0=*/1, /*,v0=*/0, /*,use_pred=*/PETSC_FALSE}; 160 161 PetscFunctionBeginUser; 162 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 163 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 164 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs"); 165 166 PetscOptionsBegin(PETSC_COMM_SELF, "", "ex43 options", ""); 167 PetscCall(PetscOptionsReal("-frequency", "Natual frequency", __FILE__, user.Omega, &user.Omega, NULL)); 168 PetscCall(PetscOptionsReal("-damping", "Damping coefficient", __FILE__, user.Xi, &user.Xi, NULL)); 169 PetscCall(PetscOptionsReal("-initial_u", "Initial displacement", __FILE__, user.u0, &user.u0, NULL)); 170 PetscCall(PetscOptionsReal("-initial_v", "Initial velocity", __FILE__, user.v0, &user.v0, NULL)); 171 PetscCall(PetscOptionsBool("-use_pred", "Use a custom predictor", __FILE__, user.use_pred, &user.use_pred, NULL)); 172 PetscOptionsEnd(); 173 174 PetscCall(TSCreate(PETSC_COMM_SELF, &ts)); 175 PetscCall(TSSetType(ts, TSALPHA2)); 176 PetscCall(TSSetMaxTime(ts, 5 * (2 * PETSC_PI))); 177 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 178 PetscCall(TSSetTimeStep(ts, 0.01)); 179 180 PetscCall(VecCreateSeq(PETSC_COMM_SELF, 1, &R)); 181 PetscCall(VecSetUp(R)); 182 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, 1, 1, NULL, &J)); 183 PetscCall(MatSetUp(J)); 184 if (user.Xi) { 185 PetscCall(TSSetI2Function(ts, R, Residual2, &user)); 186 PetscCall(TSSetI2Jacobian(ts, J, J, Tangent2, &user)); 187 } else { 188 PetscCall(TSSetIFunction(ts, R, Residual1, &user)); 189 PetscCall(TSSetIJacobian(ts, J, J, Tangent1, &user)); 190 } 191 PetscCall(VecDestroy(&R)); 192 PetscCall(MatDestroy(&J)); 193 PetscCall(TSSetSolutionFunction(ts, Solution, &user)); 194 195 PetscCall(VecCreateSeq(PETSC_COMM_SELF, 1, &U)); 196 PetscCall(VecCreateSeq(PETSC_COMM_SELF, 1, &V)); 197 PetscCall(VecGetArrayWrite(U, &u)); 198 PetscCall(VecGetArrayWrite(V, &v)); 199 u[0] = user.u0; 200 v[0] = user.v0; 201 PetscCall(VecRestoreArrayWrite(U, &u)); 202 PetscCall(VecRestoreArrayWrite(V, &v)); 203 204 if (user.use_pred) PetscCall(TSAlpha2SetPredictor(ts, Predictor, NULL)); 205 206 PetscCall(TS2SetSolution(ts, U, V)); 207 PetscCall(TSSetFromOptions(ts)); 208 PetscCall(TSSolve(ts, NULL)); 209 210 PetscCall(VecDestroy(&U)); 211 PetscCall(VecDestroy(&V)); 212 PetscCall(TSDestroy(&ts)); 213 PetscCall(PetscFinalize()); 214 return 0; 215 } 216 217 /*TEST 218 219 test: 220 suffix: a 221 args: -ts_max_steps 10 -ts_view -use_pred 222 requires: !single 223 224 test: 225 suffix: b 226 args: -ts_max_steps 10 -ts_rtol 0 -ts_atol 1e-5 -ts_adapt_type basic -ts_adapt_monitor 227 requires: !single 228 229 TEST*/ 230