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