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