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