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