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