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 PetscFunctionBegin; 45 Exact(t,user->Omega,user->Xi,user->u0,user->v0,&u,&v); 46 CHKERRQ(VecGetArray(X,&x)); 47 x[0] = (PetscScalar)u; 48 CHKERRQ(VecRestoreArray(X,&x)); 49 PetscFunctionReturn(0); 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 PetscFunctionBegin; 60 CHKERRQ(VecGetArrayRead(U,&u)); 61 CHKERRQ(VecGetArrayRead(A,&a)); 62 CHKERRQ(VecGetArrayWrite(R,&r)); 63 64 r[0] = a[0] + (Omega*Omega)*u[0]; 65 66 CHKERRQ(VecRestoreArrayRead(U,&u)); 67 CHKERRQ(VecRestoreArrayRead(A,&a)); 68 CHKERRQ(VecRestoreArrayWrite(R,&r)); 69 CHKERRQ(VecAssemblyBegin(R)); 70 CHKERRQ(VecAssemblyEnd(R)); 71 PetscFunctionReturn(0); 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 PetscFunctionBegin; 81 82 T = shiftA + (Omega*Omega); 83 84 CHKERRQ(MatSetValue(P,0,0,T,INSERT_VALUES)); 85 CHKERRQ(MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY)); 86 CHKERRQ(MatAssemblyEnd (P,MAT_FINAL_ASSEMBLY)); 87 if (J != P) { 88 CHKERRQ(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY)); 89 CHKERRQ(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY)); 90 } 91 PetscFunctionReturn(0); 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 PetscFunctionBegin; 102 CHKERRQ(VecGetArrayRead(U,&u)); 103 CHKERRQ(VecGetArrayRead(V,&v)); 104 CHKERRQ(VecGetArrayRead(A,&a)); 105 CHKERRQ(VecGetArrayWrite(R,&r)); 106 107 r[0] = a[0] + (2*Xi*Omega)*v[0] + (Omega*Omega)*u[0]; 108 109 CHKERRQ(VecRestoreArrayRead(U,&u)); 110 CHKERRQ(VecRestoreArrayRead(V,&v)); 111 CHKERRQ(VecRestoreArrayRead(A,&a)); 112 CHKERRQ(VecRestoreArrayWrite(R,&r)); 113 CHKERRQ(VecAssemblyBegin(R)); 114 CHKERRQ(VecAssemblyEnd(R)); 115 PetscFunctionReturn(0); 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 PetscFunctionBegin; 125 126 T = shiftA + shiftV * (2*Xi*Omega) + (Omega*Omega); 127 128 CHKERRQ(MatSetValue(P,0,0,T,INSERT_VALUES)); 129 CHKERRQ(MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY)); 130 CHKERRQ(MatAssemblyEnd (P,MAT_FINAL_ASSEMBLY)); 131 if (J != P) { 132 CHKERRQ(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY)); 133 CHKERRQ(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY)); 134 } 135 PetscFunctionReturn(0); 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 PetscErrorCode ierr; 148 149 CHKERRQ(PetscInitialize(&argc,&argv,NULL,help)); 150 CHKERRMPI(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 ierr = PetscOptionsBegin(PETSC_COMM_SELF,"","ex43 options","");CHKERRQ(ierr); 154 CHKERRQ(PetscOptionsReal("-frequency","Natual frequency",__FILE__,user.Omega,&user.Omega,NULL)); 155 CHKERRQ(PetscOptionsReal("-damping","Damping coefficient",__FILE__,user.Xi,&user.Xi,NULL)); 156 CHKERRQ(PetscOptionsReal("-initial_u","Initial displacement",__FILE__,user.u0,&user.u0,NULL)); 157 CHKERRQ(PetscOptionsReal("-initial_v","Initial velocity",__FILE__,user.v0,&user.v0,NULL)); 158 ierr = PetscOptionsEnd();CHKERRQ(ierr); 159 160 CHKERRQ(TSCreate(PETSC_COMM_SELF,&ts)); 161 CHKERRQ(TSSetType(ts,TSALPHA2)); 162 CHKERRQ(TSSetMaxTime(ts,5*(2*PETSC_PI))); 163 CHKERRQ(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 164 CHKERRQ(TSSetTimeStep(ts,0.01)); 165 166 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,1,&R)); 167 CHKERRQ(VecSetUp(R)); 168 CHKERRQ(MatCreateSeqDense(PETSC_COMM_SELF,1,1,NULL,&J)); 169 CHKERRQ(MatSetUp(J)); 170 if (user.Xi) { 171 CHKERRQ(TSSetI2Function(ts,R,Residual2,&user)); 172 CHKERRQ(TSSetI2Jacobian(ts,J,J,Tangent2,&user)); 173 } else { 174 CHKERRQ(TSSetIFunction(ts,R,Residual1,&user)); 175 CHKERRQ(TSSetIJacobian(ts,J,J,Tangent1,&user)); 176 } 177 CHKERRQ(VecDestroy(&R)); 178 CHKERRQ(MatDestroy(&J)); 179 CHKERRQ(TSSetSolutionFunction(ts,Solution,&user)); 180 181 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,1,&U)); 182 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,1,&V)); 183 CHKERRQ(VecGetArrayWrite(U,&u)); 184 CHKERRQ(VecGetArrayWrite(V,&v)); 185 u[0] = user.u0; 186 v[0] = user.v0; 187 CHKERRQ(VecRestoreArrayWrite(U,&u)); 188 CHKERRQ(VecRestoreArrayWrite(V,&v)); 189 190 CHKERRQ(TS2SetSolution(ts,U,V)); 191 CHKERRQ(TSSetFromOptions(ts)); 192 CHKERRQ(TSSolve(ts,NULL)); 193 194 CHKERRQ(VecDestroy(&U)); 195 CHKERRQ(VecDestroy(&V)); 196 CHKERRQ(TSDestroy(&ts)); 197 CHKERRQ(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