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