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