1 static char help[] = "Testing integrators on the simple harmonic oscillator\n"; 2 3 /* 4 In order to check long time behavior, you can give 5 6 -ts_max_steps 10000 -ts_max_time 10.0 -output_step 1000 7 8 To look at the long time behavior for different integrators, we can use the nergy monitor. Below we see that Euler is almost 100 times worse at conserving energy than the symplectic integrator of the same order. 9 10 make -f ./gmakefile test search="dm_impls_swarm_tests-ex4_bsi_1" 11 EXTRA_OPTIONS="-ts_max_steps 10000 -ts_max_time 10.0 -output_step 1000 -ts_type euler" | grep error 12 13 energy/exact energy 398.236 / 396.608 (0.4104399231) 14 energy/exact energy 1579.52 / 1573.06 (0.4104399231) 15 energy/exact energy 397.421 / 396.608 (0.2050098479) 16 energy/exact energy 1576.29 / 1573.06 (0.2050098479) 17 energy/exact energy 397.014 / 396.608 (0.1024524454) 18 energy/exact energy 1574.68 / 1573.06 (0.1024524454) 19 20 make -f ./gmakefile test search="dm_impls_swarm_tests-ex4_bsi_1" 21 EXTRA_OPTIONS="-ts_max_steps 10000 -ts_max_time 10.0 -output_step 1000" | grep error 22 23 energy/exact energy 396.579 / 396.608 (0.0074080434) 24 energy/exact energy 1572.95 / 1573.06 (0.0074080434) 25 energy/exact energy 396.593 / 396.608 (0.0037040885) 26 energy/exact energy 1573.01 / 1573.06 (0.0037040886) 27 energy/exact energy 396.601 / 396.608 (0.0018520613) 28 energy/exact energy 1573.04 / 1573.06 (0.0018520613) 29 30 We can look at third order integrators in the same way, but we need to use more steps. 31 32 make -f ./gmakefile test search="dm_impls_swarm_tests-ex4_bsi_1" 33 EXTRA_OPTIONS="-ts_max_steps 1000000 -ts_max_time 1000.0 -output_step 100000 -ts_type rk -ts_adapt_type none" | grep error 34 35 energy/exact energy 396.608 / 396.608 (0.0000013981) 36 energy/exact energy 1573.06 / 1573.06 (0.0000013981) 37 energy/exact energy 396.608 / 396.608 (0.0000001747) 38 energy/exact energy 1573.06 / 1573.06 (0.0000001748) 39 energy/exact energy 396.608 / 396.608 (0.0000000218) 40 energy/exact energy 1573.06 / 1573.06 (0.0000000218) 41 42 make -f ./gmakefile test search="dm_impls_swarm_tests-ex4_bsi_3" 43 EXTRA_OPTIONS="-ts_max_steps 1000000 -ts_max_time 1000.0 -output_step 100000 -ts_adapt_type none" | grep error 44 45 energy/exact energy 396.608 / 396.608 (0.0000000007) 46 energy/exact energy 1573.06 / 1573.06 (0.0000000007) 47 energy/exact energy 396.608 / 396.608 (0.0000000001) 48 energy/exact energy 1573.06 / 1573.06 (0.0000000001) 49 energy/exact energy 396.608 / 396.608 (0.0000000000) 50 energy/exact energy 1573.06 / 1573.06 (0.0000000000) 51 */ 52 53 #include <petscts.h> 54 #include <petscdmplex.h> 55 #include <petscdmswarm.h> 56 #include <petsc/private/dmpleximpl.h> /* For norm and dot */ 57 58 typedef struct { 59 PetscReal omega; /* Oscillation frequency omega */ 60 PetscBool error; /* Flag for printing the error */ 61 PetscInt ostep; /* print the energy at each ostep time steps */ 62 } AppCtx; 63 64 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 65 { 66 PetscFunctionBeginUser; 67 options->omega = 64.0; 68 options->error = PETSC_FALSE; 69 options->ostep = 100; 70 71 PetscOptionsBegin(comm, "", "Harmonic Oscillator Options", "DMSWARM"); 72 PetscCall(PetscOptionsReal("-omega", "Oscillator frequency", "ex4.c", options->omega, &options->omega, PETSC_NULL)); 73 PetscCall(PetscOptionsBool("-error", "Flag to print the error", "ex4.c", options->error, &options->error, NULL)); 74 PetscCall(PetscOptionsInt("-output_step", "Number of time steps between output", "ex4.c", options->ostep, &options->ostep, PETSC_NULL)); 75 PetscOptionsEnd(); 76 PetscFunctionReturn(0); 77 } 78 79 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 80 { 81 PetscFunctionBeginUser; 82 PetscCall(DMCreate(comm, dm)); 83 PetscCall(DMSetType(*dm, DMPLEX)); 84 PetscCall(DMSetFromOptions(*dm)); 85 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 86 PetscFunctionReturn(0); 87 } 88 89 static PetscErrorCode CreateSwarm(DM dm, AppCtx *user, DM *sw) 90 { 91 PetscReal v0[1] = {0.}; // Initialize velocities to 0 92 PetscInt dim; 93 94 PetscFunctionBeginUser; 95 PetscCall(DMGetDimension(dm, &dim)); 96 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw)); 97 PetscCall(DMSetType(*sw, DMSWARM)); 98 PetscCall(DMSetDimension(*sw, dim)); 99 PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC)); 100 PetscCall(DMSwarmSetCellDM(*sw, dm)); 101 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR)); 102 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "velocity", dim, PETSC_REAL)); 103 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "species", 1, PETSC_INT)); 104 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "initCoordinates", dim, PETSC_REAL)); 105 PetscCall(DMSwarmFinalizeFieldRegister(*sw)); 106 PetscCall(DMSwarmComputeLocalSizeFromOptions(*sw)); 107 PetscCall(DMSwarmInitializeCoordinates(*sw)); 108 PetscCall(DMSwarmInitializeVelocitiesFromOptions(*sw, v0)); 109 PetscCall(DMSetFromOptions(*sw)); 110 PetscCall(DMSetApplicationContext(*sw, user)); 111 PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles")); 112 PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view")); 113 { 114 Vec gc, gc0; 115 116 PetscCall(DMSwarmCreateGlobalVectorFromField(*sw, DMSwarmPICField_coor, &gc)); 117 PetscCall(DMSwarmCreateGlobalVectorFromField(*sw, "initCoordinates", &gc0)); 118 PetscCall(VecCopy(gc, gc0)); 119 PetscCall(DMSwarmDestroyGlobalVectorFromField(*sw, DMSwarmPICField_coor, &gc)); 120 PetscCall(DMSwarmDestroyGlobalVectorFromField(*sw, "initCoordinates", &gc0)); 121 } 122 PetscFunctionReturn(0); 123 } 124 125 static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec G, void *ctx) 126 { 127 const PetscReal omega = ((AppCtx *)ctx)->omega; 128 DM sw; 129 const PetscScalar *u; 130 PetscScalar *g; 131 PetscInt dim, d, Np, p; 132 133 PetscFunctionBeginUser; 134 PetscCall(TSGetDM(ts, &sw)); 135 PetscCall(DMGetDimension(sw, &dim)); 136 PetscCall(VecGetLocalSize(U, &Np)); 137 PetscCall(VecGetArray(G, &g)); 138 PetscCall(VecGetArrayRead(U, &u)); 139 Np /= 2 * dim; 140 for (p = 0; p < Np; ++p) { 141 for (d = 0; d < dim; ++d) { 142 g[(p * 2 + 0) * dim + d] = u[(p * 2 + 1) * dim + d]; 143 g[(p * 2 + 1) * dim + d] = -PetscSqr(omega) * u[(p * 2 + 0) * dim + d]; 144 } 145 } 146 PetscCall(VecRestoreArrayRead(U, &u)); 147 PetscCall(VecRestoreArray(G, &g)); 148 PetscFunctionReturn(0); 149 } 150 151 /* J_{ij} = dF_i/dx_j 152 J_p = ( 0 1) 153 (-w^2 0) 154 */ 155 static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat J, Mat P, void *ctx) 156 { 157 PetscScalar vals[4] = {0., 1., -PetscSqr(((AppCtx *)ctx)->omega), 0.}; 158 DM sw; 159 PetscInt dim, d, Np, p, rStart; 160 161 PetscFunctionBeginUser; 162 PetscCall(TSGetDM(ts, &sw)); 163 PetscCall(DMGetDimension(sw, &dim)); 164 PetscCall(VecGetLocalSize(U, &Np)); 165 PetscCall(MatGetOwnershipRange(J, &rStart, NULL)); 166 Np /= 2 * dim; 167 for (p = 0; p < Np; ++p) { 168 for (d = 0; d < dim; ++d) { 169 const PetscInt rows[2] = {(p * 2 + 0) * dim + d + rStart, (p * 2 + 1) * dim + d + rStart}; 170 PetscCall(MatSetValues(J, 2, rows, 2, rows, vals, INSERT_VALUES)); 171 } 172 } 173 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 174 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 175 PetscFunctionReturn(0); 176 } 177 178 static PetscErrorCode RHSFunctionX(TS ts, PetscReal t, Vec V, Vec Xres, void *ctx) 179 { 180 const PetscScalar *v; 181 PetscScalar *xres; 182 PetscInt Np, p; 183 184 PetscFunctionBeginUser; 185 PetscCall(VecGetLocalSize(Xres, &Np)); 186 PetscCall(VecGetArrayRead(V, &v)); 187 PetscCall(VecGetArray(Xres, &xres)); 188 for (p = 0; p < Np; ++p) xres[p] = v[p]; 189 PetscCall(VecRestoreArrayRead(V, &v)); 190 PetscCall(VecRestoreArray(Xres, &xres)); 191 PetscFunctionReturn(0); 192 } 193 194 static PetscErrorCode RHSFunctionV(TS ts, PetscReal t, Vec X, Vec Vres, void *ctx) 195 { 196 const PetscReal omega = ((AppCtx *)ctx)->omega; 197 const PetscScalar *x; 198 PetscScalar *vres; 199 PetscInt Np, p; 200 201 PetscFunctionBeginUser; 202 PetscCall(VecGetArray(Vres, &vres)); 203 PetscCall(VecGetArrayRead(X, &x)); 204 PetscCall(VecGetLocalSize(Vres, &Np)); 205 for (p = 0; p < Np; ++p) vres[p] = -PetscSqr(omega) * x[p]; 206 PetscCall(VecRestoreArrayRead(X, &x)); 207 PetscCall(VecRestoreArray(Vres, &vres)); 208 PetscFunctionReturn(0); 209 } 210 211 PetscErrorCode RHSJacobianS(TS ts, PetscReal t, Vec U, Mat S, void *ctx) 212 { 213 PetscScalar vals[4] = {0., 1., -1., 0.}; 214 DM sw; 215 PetscInt dim, d, Np, p, rStart; 216 217 PetscFunctionBeginUser; 218 PetscCall(TSGetDM(ts, &sw)); 219 PetscCall(DMGetDimension(sw, &dim)); 220 PetscCall(VecGetLocalSize(U, &Np)); 221 PetscCall(MatGetOwnershipRange(S, &rStart, NULL)); 222 Np /= 2 * dim; 223 for (p = 0; p < Np; ++p) { 224 for (d = 0; d < dim; ++d) { 225 const PetscInt rows[2] = {(p * 2 + 0) * dim + d + rStart, (p * 2 + 1) * dim + d + rStart}; 226 PetscCall(MatSetValues(S, 2, rows, 2, rows, vals, INSERT_VALUES)); 227 } 228 } 229 PetscCall(MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY)); 230 PetscCall(MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY)); 231 PetscFunctionReturn(0); 232 } 233 234 PetscErrorCode RHSObjectiveF(TS ts, PetscReal t, Vec U, PetscScalar *F, void *ctx) 235 { 236 const PetscReal omega = ((AppCtx *)ctx)->omega; 237 DM sw; 238 const PetscScalar *u; 239 PetscInt dim, Np, p; 240 241 PetscFunctionBeginUser; 242 PetscCall(TSGetDM(ts, &sw)); 243 PetscCall(DMGetDimension(sw, &dim)); 244 PetscCall(VecGetArrayRead(U, &u)); 245 PetscCall(VecGetLocalSize(U, &Np)); 246 Np /= 2 * dim; 247 for (p = 0; p < Np; ++p) { 248 const PetscReal x2 = DMPlex_DotRealD_Internal(dim, &u[(p * 2 + 0) * dim], &u[(p * 2 + 0) * dim]); 249 const PetscReal v2 = DMPlex_DotRealD_Internal(dim, &u[(p * 2 + 1) * dim], &u[(p * 2 + 1) * dim]); 250 const PetscReal E = 0.5 * (v2 + PetscSqr(omega) * x2); 251 252 *F += E; 253 } 254 PetscCall(VecRestoreArrayRead(U, &u)); 255 PetscFunctionReturn(0); 256 } 257 258 /* dF/dx = omega^2 x dF/dv = v */ 259 PetscErrorCode RHSFunctionG(TS ts, PetscReal t, Vec U, Vec G, void *ctx) 260 { 261 const PetscReal omega = ((AppCtx *)ctx)->omega; 262 DM sw; 263 const PetscScalar *u; 264 PetscScalar *g; 265 PetscInt dim, d, Np, p; 266 267 PetscFunctionBeginUser; 268 PetscCall(TSGetDM(ts, &sw)); 269 PetscCall(DMGetDimension(sw, &dim)); 270 PetscCall(VecGetArray(G, &g)); 271 PetscCall(VecGetArrayRead(U, &u)); 272 PetscCall(VecGetLocalSize(U, &Np)); 273 Np /= 2 * dim; 274 for (p = 0; p < Np; ++p) { 275 for (d = 0; d < dim; ++d) { 276 g[(p * 2 + 0) * dim + d] = PetscSqr(omega) * u[(p * 2 + 0) * dim + d]; 277 g[(p * 2 + 1) * dim + d] = u[(p * 2 + 1) * dim + d]; 278 } 279 } 280 PetscCall(VecRestoreArrayRead(U, &u)); 281 PetscCall(VecRestoreArray(G, &g)); 282 PetscFunctionReturn(0); 283 } 284 285 static PetscErrorCode CreateSolution(TS ts) 286 { 287 DM sw; 288 Vec u; 289 PetscInt dim, Np; 290 291 PetscFunctionBegin; 292 PetscCall(TSGetDM(ts, &sw)); 293 PetscCall(DMGetDimension(sw, &dim)); 294 PetscCall(DMSwarmGetLocalSize(sw, &Np)); 295 PetscCall(VecCreate(PETSC_COMM_WORLD, &u)); 296 PetscCall(VecSetBlockSize(u, dim)); 297 PetscCall(VecSetSizes(u, 2 * Np * dim, PETSC_DECIDE)); 298 PetscCall(VecSetUp(u)); 299 PetscCall(TSSetSolution(ts, u)); 300 PetscCall(VecDestroy(&u)); 301 PetscFunctionReturn(0); 302 } 303 304 static PetscErrorCode SetProblem(TS ts) 305 { 306 AppCtx *user; 307 DM sw; 308 309 PetscFunctionBegin; 310 PetscCall(TSGetDM(ts, &sw)); 311 PetscCall(DMGetApplicationContext(sw, (void **)&user)); 312 // Define unified system for (X, V) 313 { 314 Mat J; 315 PetscInt dim, Np; 316 317 PetscCall(DMGetDimension(sw, &dim)); 318 PetscCall(DMSwarmGetLocalSize(sw, &Np)); 319 PetscCall(MatCreate(PETSC_COMM_WORLD, &J)); 320 PetscCall(MatSetSizes(J, 2 * Np * dim, 2 * Np * dim, PETSC_DECIDE, PETSC_DECIDE)); 321 PetscCall(MatSetBlockSize(J, 2 * dim)); 322 PetscCall(MatSetFromOptions(J)); 323 PetscCall(MatSetUp(J)); 324 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, user)); 325 PetscCall(TSSetRHSJacobian(ts, J, J, RHSJacobian, user)); 326 PetscCall(MatDestroy(&J)); 327 } 328 // Define split system for X and V 329 { 330 IS isx, isv, istmp; 331 const PetscInt *idx; 332 PetscInt dim, Np; 333 334 PetscCall(DMGetDimension(sw, &dim)); 335 PetscCall(DMSwarmGetLocalSize(sw, &Np)); 336 PetscCall(ISCreateStride(PETSC_COMM_SELF, Np, 0, 2, &istmp)); 337 PetscCall(ISGetIndices(istmp, &idx)); 338 PetscCall(ISCreateBlock(PETSC_COMM_SELF, dim, Np, idx, PETSC_COPY_VALUES, &isx)); 339 PetscCall(ISRestoreIndices(istmp, &idx)); 340 PetscCall(ISDestroy(&istmp)); 341 PetscCall(ISCreateStride(PETSC_COMM_SELF, Np, 1, 2, &istmp)); 342 PetscCall(ISGetIndices(istmp, &idx)); 343 PetscCall(ISCreateBlock(PETSC_COMM_SELF, dim, Np, idx, PETSC_COPY_VALUES, &isv)); 344 PetscCall(ISRestoreIndices(istmp, &idx)); 345 PetscCall(ISDestroy(&istmp)); 346 PetscCall(TSRHSSplitSetIS(ts, "position", isx)); 347 PetscCall(TSRHSSplitSetIS(ts, "momentum", isv)); 348 PetscCall(ISDestroy(&isx)); 349 PetscCall(ISDestroy(&isv)); 350 PetscCall(TSRHSSplitSetRHSFunction(ts, "position", NULL, RHSFunctionX, user)); 351 PetscCall(TSRHSSplitSetRHSFunction(ts, "momentum", NULL, RHSFunctionV, user)); 352 } 353 // Define symplectic formulation U_t = S . G, where G = grad F 354 { 355 PetscCall(TSDiscGradSetFormulation(ts, RHSJacobianS, RHSObjectiveF, RHSFunctionG, user)); 356 } 357 PetscFunctionReturn(0); 358 } 359 360 static PetscErrorCode InitializeSolve(TS ts, Vec u) 361 { 362 DM sw; 363 Vec gc, gc0; 364 IS isx, isv; 365 AppCtx *user; 366 367 PetscFunctionBeginUser; 368 PetscCall(TSGetDM(ts, &sw)); 369 PetscCall(DMGetApplicationContext(sw, &user)); 370 { 371 PetscReal v0[1] = {1.}; 372 373 PetscCall(DMSwarmInitializeCoordinates(sw)); 374 PetscCall(DMSwarmInitializeVelocitiesFromOptions(sw, v0)); 375 PetscCall(SetProblem(ts)); 376 } 377 PetscCall(TSRHSSplitGetIS(ts, "position", &isx)); 378 PetscCall(TSRHSSplitGetIS(ts, "momentum", &isv)); 379 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc)); 380 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "initCoordinates", &gc0)); 381 PetscCall(VecCopy(gc, gc0)); 382 PetscCall(VecISCopy(u, isx, SCATTER_FORWARD, gc)); 383 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc)); 384 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "initCoordinates", &gc0)); 385 PetscCall(VecISSet(u, isv, 0.)); 386 PetscFunctionReturn(0); 387 } 388 389 static PetscErrorCode ComputeError(TS ts, Vec U, Vec E) 390 { 391 MPI_Comm comm; 392 DM sw; 393 AppCtx *user; 394 const PetscScalar *u; 395 const PetscReal *coords; 396 PetscScalar *e; 397 PetscReal t; 398 PetscInt dim, d, Np, p; 399 400 PetscFunctionBeginUser; 401 PetscCall(PetscObjectGetComm((PetscObject)ts, &comm)); 402 PetscCall(TSGetDM(ts, &sw)); 403 PetscCall(DMGetApplicationContext(sw, &user)); 404 PetscCall(DMGetDimension(sw, &dim)); 405 PetscCall(TSGetSolveTime(ts, &t)); 406 PetscCall(VecGetArray(E, &e)); 407 PetscCall(VecGetArrayRead(U, &u)); 408 PetscCall(VecGetLocalSize(U, &Np)); 409 PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords)); 410 Np /= 2 * dim; 411 for (p = 0; p < Np; ++p) { 412 const PetscReal omega = user->omega; 413 const PetscReal ct = PetscCosReal(omega * t); 414 const PetscReal st = PetscSinReal(omega * t); 415 const PetscReal x0 = DMPlex_NormD_Internal(dim, &coords[p * dim]); 416 const PetscReal ex = x0 * ct; 417 const PetscReal ev = -x0 * omega * st; 418 const PetscReal x = DMPlex_DotRealD_Internal(dim, &u[(p * 2 + 0) * dim], &coords[p * dim]) / x0; 419 const PetscReal v = DMPlex_DotRealD_Internal(dim, &u[(p * 2 + 1) * dim], &coords[p * dim]) / x0; 420 421 if (user->error) { 422 const PetscReal en = 0.5 * (v * v + PetscSqr(omega) * x * x); 423 const PetscReal exen = 0.5 * PetscSqr(omega * x0); 424 PetscCall(PetscPrintf(comm, "p%" PetscInt_FMT " error [%.2g %.2g] sol [%.6lf %.6lf] exact [%.6lf %.6lf] energy/exact energy %g / %g (%.10lf%%)\n", p, (double)PetscAbsReal(x - ex), (double)PetscAbsReal(v - ev), (double)x, (double)v, (double)ex, (double)ev, (double)en, (double)exen, (double)(PetscAbsReal(exen - en) * 100. / exen))); 425 } 426 for (d = 0; d < dim; ++d) { 427 e[(p * 2 + 0) * dim + d] = u[(p * 2 + 0) * dim + d] - coords[p * dim + d] * ct; 428 e[(p * 2 + 1) * dim + d] = u[(p * 2 + 1) * dim + d] + coords[p * dim + d] * omega * st; 429 } 430 } 431 PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords)); 432 PetscCall(VecRestoreArrayRead(U, &u)); 433 PetscCall(VecRestoreArray(E, &e)); 434 PetscFunctionReturn(0); 435 } 436 437 static PetscErrorCode EnergyMonitor(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx) 438 { 439 const PetscReal omega = ((AppCtx *)ctx)->omega; 440 const PetscInt ostep = ((AppCtx *)ctx)->ostep; 441 DM sw; 442 const PetscScalar *u; 443 PetscReal dt; 444 PetscInt dim, Np, p; 445 MPI_Comm comm; 446 447 PetscFunctionBeginUser; 448 if (step % ostep == 0) { 449 PetscCall(PetscObjectGetComm((PetscObject)ts, &comm)); 450 PetscCall(TSGetDM(ts, &sw)); 451 PetscCall(TSGetTimeStep(ts, &dt)); 452 PetscCall(DMGetDimension(sw, &dim)); 453 PetscCall(VecGetArrayRead(U, &u)); 454 PetscCall(VecGetLocalSize(U, &Np)); 455 Np /= 2 * dim; 456 if (!step) PetscCall(PetscPrintf(comm, "Time Step Part Energy Mod Energy\n")); 457 for (p = 0; p < Np; ++p) { 458 const PetscReal x2 = DMPlex_DotRealD_Internal(dim, &u[(p * 2 + 0) * dim], &u[(p * 2 + 0) * dim]); 459 const PetscReal v2 = DMPlex_DotRealD_Internal(dim, &u[(p * 2 + 1) * dim], &u[(p * 2 + 1) * dim]); 460 const PetscReal E = 0.5 * (v2 + PetscSqr(omega) * x2); 461 const PetscReal mE = 0.5 * (v2 + PetscSqr(omega) * x2 - PetscSqr(omega) * dt * PetscSqrtReal(x2 * v2)); 462 463 PetscCall(PetscPrintf(comm, "%.6lf %4" PetscInt_FMT " %4" PetscInt_FMT " %10.4lf %10.4lf\n", (double)t, step, p, (double)E, (double)mE)); 464 } 465 PetscCall(VecRestoreArrayRead(U, &u)); 466 } 467 PetscFunctionReturn(0); 468 } 469 470 int main(int argc, char **argv) 471 { 472 DM dm, sw; 473 TS ts; 474 Vec u; 475 AppCtx user; 476 477 PetscFunctionBeginUser; 478 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 479 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 480 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 481 PetscCall(CreateSwarm(dm, &user, &sw)); 482 PetscCall(DMSetApplicationContext(sw, &user)); 483 484 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 485 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 486 PetscCall(TSSetDM(ts, sw)); 487 PetscCall(TSSetMaxTime(ts, 0.1)); 488 PetscCall(TSSetTimeStep(ts, 0.00001)); 489 PetscCall(TSSetMaxSteps(ts, 100)); 490 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 491 PetscCall(TSMonitorSet(ts, EnergyMonitor, &user, NULL)); 492 PetscCall(TSSetFromOptions(ts)); 493 PetscCall(TSSetComputeInitialCondition(ts, InitializeSolve)); 494 PetscCall(TSSetComputeExactError(ts, ComputeError)); 495 496 PetscCall(CreateSolution(ts)); 497 PetscCall(TSGetSolution(ts, &u)); 498 PetscCall(TSComputeInitialCondition(ts, u)); 499 PetscCall(TSSolve(ts, NULL)); 500 501 PetscCall(TSDestroy(&ts)); 502 PetscCall(DMDestroy(&sw)); 503 PetscCall(DMDestroy(&dm)); 504 PetscCall(PetscFinalize()); 505 return 0; 506 } 507 508 /*TEST 509 510 build: 511 requires: triangle !single !complex 512 513 testset: 514 args: -dm_plex_dim 1 -dm_plex_box_faces 1 -dm_plex_box_lower -1 -dm_plex_box_upper 1 \ 515 -dm_swarm_num_particles 2 -dm_swarm_coordinate_density constant \ 516 -ts_type basicsymplectic -ts_convergence_estimate -convest_num_refine 2 \ 517 -dm_view -output_step 50 -error 518 test: 519 suffix: bsi_1 520 args: -ts_basicsymplectic_type 1 521 test: 522 suffix: bsi_2 523 args: -ts_basicsymplectic_type 2 524 test: 525 suffix: bsi_3 526 args: -ts_basicsymplectic_type 3 527 test: 528 suffix: bsi_4 529 args: -ts_basicsymplectic_type 4 -ts_dt 0.0001 530 531 testset: 532 args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -dm_plex_box_lower -1,-1 -dm_plex_box_upper 1,1 \ 533 -dm_swarm_num_particles 2 -dm_swarm_coordinate_density constant \ 534 -ts_type basicsymplectic -ts_convergence_estimate -convest_num_refine 2 \ 535 -dm_view -output_step 50 -error 536 test: 537 suffix: bsi_2d_1 538 args: -ts_basicsymplectic_type 1 539 test: 540 suffix: bsi_2d_2 541 args: -ts_basicsymplectic_type 2 542 test: 543 suffix: bsi_2d_3 544 args: -ts_basicsymplectic_type 3 545 test: 546 suffix: bsi_2d_4 547 args: -ts_basicsymplectic_type 4 -ts_dt 0.0001 548 549 testset: 550 args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -dm_plex_box_lower -1,-1,-1 -dm_plex_box_upper 1,1,1 \ 551 -dm_swarm_num_particles 2 -dm_swarm_coordinate_density constant \ 552 -ts_type basicsymplectic -ts_convergence_estimate -convest_num_refine 2 \ 553 -dm_view -output_step 50 -error 554 test: 555 suffix: bsi_3d_1 556 args: -ts_basicsymplectic_type 1 557 test: 558 suffix: bsi_3d_2 559 args: -ts_basicsymplectic_type 2 560 test: 561 suffix: bsi_3d_3 562 args: -ts_basicsymplectic_type 3 563 test: 564 suffix: bsi_3d_4 565 args: -ts_basicsymplectic_type 4 -ts_dt 0.0001 566 567 testset: 568 args: -dm_swarm_num_particles 2 -dm_swarm_coordinate_density constant \ 569 -ts_type theta -ts_theta_theta 0.5 -ts_convergence_estimate -convest_num_refine 2 \ 570 -mat_type baij -ksp_error_if_not_converged -pc_type lu \ 571 -dm_view -output_step 50 -error 572 test: 573 suffix: im_1d 574 args: -dm_plex_dim 1 -dm_plex_box_faces 1 -dm_plex_box_lower -1 -dm_plex_box_upper 1 575 test: 576 suffix: im_2d 577 args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -dm_plex_box_lower -1,-1 -dm_plex_box_upper 1,1 578 test: 579 suffix: im_3d 580 args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1 -dm_plex_box_lower -1,-1,-1 -dm_plex_box_upper 1,1,1 581 582 testset: 583 args: -dm_swarm_num_particles 2 -dm_swarm_coordinate_density constant \ 584 -ts_type discgrad -ts_discgrad_gonzalez -ts_convergence_estimate -convest_num_refine 2 \ 585 -mat_type baij -ksp_error_if_not_converged -pc_type lu \ 586 -dm_view -output_step 50 -error 587 test: 588 suffix: dg_1d 589 args: -dm_plex_dim 1 -dm_plex_box_faces 1 -dm_plex_box_lower -1 -dm_plex_box_upper 1 590 test: 591 suffix: dg_2d 592 args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -dm_plex_box_lower -1,-1 -dm_plex_box_upper 1,1 593 test: 594 suffix: dg_3d 595 args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1 -dm_plex_box_lower -1,-1,-1 -dm_plex_box_upper 1,1,1 596 597 TEST*/ 598