1 2 static char help[] = "Solves the motion of spring.\n\ 3 Input parameters include:\n"; 4 5 /* 6 Concepts: TS^Separable Hamiltonian problems 7 Concepts: TS^Symplectic intergrators 8 Processors: 1 9 */ 10 /* ------------------------------------------------------------------------ 11 12 This program solves the motion of spring by Hooke's law 13 x' = f(t,v) = v 14 v' = g(t,x) = -omega^2*x 15 on the interval 0 <= t <= 0.1, with the initial conditions 16 x(0) = 0.2, x'(0) = v(0) = 0, 17 and 18 omega = 64. 19 The exact solution is 20 x(t) = A*sin(t*omega) + B*cos(t*omega) 21 where A and B are constants that can be determined from the initial conditions. 22 In this case, B=0.2, A=0. 23 24 Notes: 25 This code demonstrates the TS solver interface to solve a separable Hamiltonian 26 system, which can be split into two subsystems involving two coupling components, 27 named generailized momentum and generailized position respectively. 28 Using a symplectic intergrator can preserve energy 29 E = (v^2+omega^2*x^2-omega^2*h*v*x)/2 30 ------------------------------------------------------------------------- */ 31 32 #include <petscts.h> 33 #include <petscvec.h> 34 35 typedef struct _n_User *User; 36 struct _n_User { 37 PetscReal omega; 38 PetscInt nts; /* print the energy at each nts time steps */ 39 }; 40 41 /* 42 User-defined routines. 43 The first RHS function provides f(t,x), the residual for the generalized momentum, 44 and the second one provides g(t,v), the residual for the generalized position. 45 */ 46 static PetscErrorCode RHSFunction2(TS ts,PetscReal t,Vec X,Vec Vres,void *ctx) 47 { 48 User user = (User)ctx; 49 const PetscScalar *x; 50 PetscScalar *vres; 51 52 PetscFunctionBeginUser; 53 PetscCall(VecGetArrayRead(X,&x)); 54 PetscCall(VecGetArray(Vres,&vres)); 55 vres[0] = -user->omega*user->omega*x[0]; 56 PetscCall(VecRestoreArray(Vres,&vres)); 57 PetscCall(VecRestoreArrayRead(X,&x)); 58 PetscFunctionReturn(0); 59 } 60 61 static PetscErrorCode RHSFunction1(TS ts,PetscReal t,Vec V,Vec Xres,void *ctx) 62 { 63 const PetscScalar *v; 64 PetscScalar *xres; 65 66 PetscFunctionBeginUser; 67 PetscCall(VecGetArray(Xres,&xres)); 68 PetscCall(VecGetArrayRead(V,&v)); 69 xres[0] = v[0]; 70 PetscCall(VecRestoreArrayRead(V,&v)); 71 PetscCall(VecRestoreArray(Xres,&xres)); 72 PetscFunctionReturn(0); 73 } 74 75 static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec R,void *ctx) 76 { 77 User user = (User)ctx; 78 const PetscScalar *u; 79 PetscScalar *r; 80 81 PetscFunctionBeginUser; 82 PetscCall(VecGetArrayRead(U,&u)); 83 PetscCall(VecGetArray(R,&r)); 84 r[0] = u[1]; 85 r[1] = -user->omega*user->omega*u[0]; 86 PetscCall(VecRestoreArrayRead(U,&u)); 87 PetscCall(VecRestoreArray(R,&r)); 88 PetscFunctionReturn(0); 89 } 90 91 /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */ 92 static PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec U,void *ctx) 93 { 94 const PetscScalar *u; 95 PetscReal dt; 96 PetscScalar energy,menergy; 97 User user = (User)ctx; 98 99 PetscFunctionBeginUser; 100 if (step%user->nts == 0) { 101 PetscCall(TSGetTimeStep(ts,&dt)); 102 PetscCall(VecGetArrayRead(U,&u)); 103 menergy = (u[1]*u[1]+user->omega*user->omega*u[0]*u[0]-user->omega*user->omega*dt*u[0]*u[1])/2.; 104 energy = (u[1]*u[1]+user->omega*user->omega*u[0]*u[0])/2.; 105 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"At time %.6lf, Energy = %8g, Modified Energy = %8g\n",t,(double)energy,(double)menergy)); 106 PetscCall(VecRestoreArrayRead(U,&u)); 107 } 108 PetscFunctionReturn(0); 109 } 110 111 int main(int argc,char **argv) 112 { 113 TS ts; /* nonlinear solver */ 114 Vec U; /* solution, residual vectors */ 115 IS is1,is2; 116 PetscInt nindices[1]; 117 PetscReal ftime = 0.1; 118 PetscBool monitor = PETSC_FALSE; 119 PetscScalar *u_ptr; 120 PetscMPIInt size; 121 struct _n_User user; 122 PetscErrorCode ierr; 123 124 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 125 Initialize program 126 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 127 PetscCall(PetscInitialize(&argc,&argv,NULL,help)); 128 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 129 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 130 131 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 132 Set runtime options 133 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 134 user.omega = 64.; 135 user.nts = 100; 136 PetscCall(PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL)); 137 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Physical parameters",NULL);PetscCall(ierr); 138 PetscCall(PetscOptionsReal("-omega","parameter","<64>",user.omega,&user.omega,PETSC_NULL)); 139 PetscCall(PetscOptionsInt("-next_output","time steps for next output point","<100>",user.nts,&user.nts,PETSC_NULL)); 140 ierr = PetscOptionsEnd();PetscCall(ierr); 141 142 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 143 Create necessary matrix and vectors, solve same ODE on every process 144 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 145 PetscCall(VecCreateSeq(PETSC_COMM_SELF,2,&U)); 146 nindices[0] = 0; 147 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,1,nindices,PETSC_COPY_VALUES,&is1)); 148 nindices[0] = 1; 149 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,1,nindices,PETSC_COPY_VALUES,&is2)); 150 151 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 152 Create timestepping solver context 153 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 154 PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 155 PetscCall(TSSetType(ts,TSBASICSYMPLECTIC)); 156 PetscCall(TSRHSSplitSetIS(ts,"position",is1)); 157 PetscCall(TSRHSSplitSetIS(ts,"momentum",is2)); 158 PetscCall(TSRHSSplitSetRHSFunction(ts,"position",NULL,RHSFunction1,&user)); 159 PetscCall(TSRHSSplitSetRHSFunction(ts,"momentum",NULL,RHSFunction2,&user)); 160 PetscCall(TSSetRHSFunction(ts,NULL,RHSFunction,&user)); 161 162 PetscCall(TSSetMaxTime(ts,ftime)); 163 PetscCall(TSSetTimeStep(ts,0.0001)); 164 PetscCall(TSSetMaxSteps(ts,1000)); 165 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP)); 166 if (monitor) { 167 PetscCall(TSMonitorSet(ts,Monitor,&user,NULL)); 168 } 169 170 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 171 Set initial conditions 172 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 173 PetscCall(VecGetArray(U,&u_ptr)); 174 u_ptr[0] = 0.2; 175 u_ptr[1] = 0.0; 176 PetscCall(VecRestoreArray(U,&u_ptr)); 177 178 PetscCall(TSSetTime(ts,0.0)); 179 180 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 181 Set runtime options 182 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 183 PetscCall(TSSetFromOptions(ts)); 184 185 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 186 Solve nonlinear system 187 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 188 PetscCall(TSSolve(ts,U)); 189 PetscCall(TSGetSolveTime(ts,&ftime)); 190 PetscCall(VecView(U,PETSC_VIEWER_STDOUT_WORLD)); 191 192 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"The exact solution at time %.6lf is [%g %g]\n",(double)ftime,(double)0.2*PetscCosReal(user.omega*ftime),(double)-0.2*user.omega*PetscSinReal(user.omega*ftime))); 193 194 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 195 Free work space. All PETSc objects should be destroyed when they 196 are no longer needed. 197 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 198 PetscCall(VecDestroy(&U)); 199 PetscCall(TSDestroy(&ts)); 200 PetscCall(ISDestroy(&is1)); 201 PetscCall(ISDestroy(&is2)); 202 PetscCall(PetscFinalize()); 203 return 0; 204 } 205 206 /*TEST 207 build: 208 requires: !single !complex 209 210 test: 211 args: -ts_basicsymplectic_type 1 -monitor 212 213 test: 214 suffix: 2 215 args: -ts_basicsymplectic_type 2 -monitor 216 217 TEST*/ 218