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 PetscErrorCode ierr; 52 53 PetscFunctionBeginUser; 54 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 55 ierr = VecGetArray(Vres,&vres);CHKERRQ(ierr); 56 vres[0] = -user->omega*user->omega*x[0]; 57 ierr = VecRestoreArray(Vres,&vres);CHKERRQ(ierr); 58 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 59 PetscFunctionReturn(0); 60 } 61 62 static PetscErrorCode RHSFunction1(TS ts,PetscReal t,Vec V,Vec Xres,void *ctx) 63 { 64 const PetscScalar *v; 65 PetscScalar *xres; 66 PetscErrorCode ierr; 67 68 PetscFunctionBeginUser; 69 ierr = VecGetArray(Xres,&xres);CHKERRQ(ierr); 70 ierr = VecGetArrayRead(V,&v);CHKERRQ(ierr); 71 xres[0] = v[0]; 72 ierr = VecRestoreArrayRead(V,&v);CHKERRQ(ierr); 73 ierr = VecRestoreArray(Xres,&xres);CHKERRQ(ierr); 74 PetscFunctionReturn(0); 75 } 76 77 static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec R,void *ctx) 78 { 79 User user = (User)ctx; 80 const PetscScalar *u; 81 PetscScalar *r; 82 PetscErrorCode ierr; 83 84 PetscFunctionBeginUser; 85 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 86 ierr = VecGetArray(R,&r);CHKERRQ(ierr); 87 r[0] = u[1]; 88 r[1] = -user->omega*user->omega*u[0]; 89 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 90 ierr = VecRestoreArray(R,&r);CHKERRQ(ierr); 91 PetscFunctionReturn(0); 92 } 93 94 /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */ 95 static PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec U,void *ctx) 96 { 97 PetscErrorCode ierr; 98 const PetscScalar *u; 99 PetscReal dt; 100 PetscScalar energy,menergy; 101 User user = (User)ctx; 102 103 PetscFunctionBeginUser; 104 if (step%user->nts == 0) { 105 ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr); 106 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 107 menergy = (u[1]*u[1]+user->omega*user->omega*u[0]*u[0]-user->omega*user->omega*dt*u[0]*u[1])/2.; 108 energy = (u[1]*u[1]+user->omega*user->omega*u[0]*u[0])/2.; 109 ierr = PetscPrintf(PETSC_COMM_WORLD,"At time %.6lf, Energy = %8g, Modified Energy = %8g\n",t,(double)energy,(double)menergy);CHKERRQ(ierr); 110 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 111 } 112 PetscFunctionReturn(0); 113 } 114 115 int main(int argc,char **argv) 116 { 117 TS ts; /* nonlinear solver */ 118 Vec U; /* solution, residual vectors */ 119 IS is1,is2; 120 PetscInt nindices[1]; 121 PetscReal ftime = 0.1; 122 PetscBool monitor = PETSC_FALSE; 123 PetscScalar *u_ptr; 124 PetscMPIInt size; 125 struct _n_User user; 126 PetscErrorCode ierr; 127 128 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 129 Initialize program 130 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 131 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; 132 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 133 if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 134 135 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 136 Set runtime options 137 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 138 user.omega = 64.; 139 user.nts = 100; 140 ierr = PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);CHKERRQ(ierr); 141 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Physical parameters",NULL);CHKERRQ(ierr); 142 ierr = PetscOptionsReal("-omega","parameter","<64>",user.omega,&user.omega,PETSC_NULL);CHKERRQ(ierr); 143 ierr = PetscOptionsInt("-next_output","time steps for next output point","<100>",user.nts,&user.nts,PETSC_NULL);CHKERRQ(ierr); 144 ierr = PetscOptionsEnd();CHKERRQ(ierr); 145 146 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 147 Create necessary matrix and vectors, solve same ODE on every process 148 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 149 ierr = VecCreateSeq(PETSC_COMM_SELF,2,&U);CHKERRQ(ierr); 150 nindices[0] = 0; 151 ierr = ISCreateGeneral(PETSC_COMM_SELF,1,nindices,PETSC_COPY_VALUES,&is1);CHKERRQ(ierr); 152 nindices[0] = 1; 153 ierr = ISCreateGeneral(PETSC_COMM_SELF,1,nindices,PETSC_COPY_VALUES,&is2);CHKERRQ(ierr); 154 155 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 156 Create timestepping solver context 157 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 158 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 159 ierr = TSSetType(ts,TSBASICSYMPLECTIC);CHKERRQ(ierr); 160 ierr = TSRHSSplitSetIS(ts,"position",is1);CHKERRQ(ierr); 161 ierr = TSRHSSplitSetIS(ts,"momentum",is2);CHKERRQ(ierr); 162 ierr = TSRHSSplitSetRHSFunction(ts,"position",NULL,RHSFunction1,&user);CHKERRQ(ierr); 163 ierr = TSRHSSplitSetRHSFunction(ts,"momentum",NULL,RHSFunction2,&user);CHKERRQ(ierr); 164 ierr = TSSetRHSFunction(ts,NULL,RHSFunction,&user);CHKERRQ(ierr); 165 166 ierr = TSSetMaxTime(ts,ftime);CHKERRQ(ierr); 167 ierr = TSSetTimeStep(ts,0.0001);CHKERRQ(ierr); 168 ierr = TSSetMaxSteps(ts,1000);CHKERRQ(ierr); 169 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr); 170 if (monitor) { 171 ierr = TSMonitorSet(ts,Monitor,&user,NULL);CHKERRQ(ierr); 172 } 173 174 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 175 Set initial conditions 176 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 177 ierr = VecGetArray(U,&u_ptr);CHKERRQ(ierr); 178 u_ptr[0] = 0.2; 179 u_ptr[1] = 0.0; 180 ierr = VecRestoreArray(U,&u_ptr);CHKERRQ(ierr); 181 182 ierr = TSSetTime(ts,0.0);CHKERRQ(ierr); 183 184 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 185 Set runtime options 186 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 187 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 188 189 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 190 Solve nonlinear system 191 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 192 ierr = TSSolve(ts,U);CHKERRQ(ierr); 193 ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr); 194 ierr = VecView(U,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 195 196 ierr = 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));CHKERRQ(ierr); 197 198 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 199 Free work space. All PETSc objects should be destroyed when they 200 are no longer needed. 201 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 202 ierr = VecDestroy(&U);CHKERRQ(ierr); 203 ierr = TSDestroy(&ts);CHKERRQ(ierr); 204 ierr = ISDestroy(&is1);CHKERRQ(ierr); 205 ierr = ISDestroy(&is2);CHKERRQ(ierr); 206 ierr = PetscFinalize(); 207 return ierr; 208 } 209 210 /*TEST 211 build: 212 requires: !single !complex 213 214 test: 215 args: -ts_basicsymplectic_type 1 -monitor 216 217 test: 218 suffix: 2 219 args: -ts_basicsymplectic_type 2 -monitor 220 221 TEST*/ 222