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