1 static char help[] = "Comparing basic symplectic, theta and discrete gradients interators on a simple hamiltonian system (harmonic oscillator) with particles\n"; 2 3 #include <petscdmplex.h> 4 #include <petsc/private/dmpleximpl.h> /* For norm */ 5 #include <petsc/private/petscfeimpl.h> /* For CoordinatesRefToReal() */ 6 #include <petscdmswarm.h> 7 #include <petscts.h> 8 9 typedef struct { 10 PetscInt dim; /* The topological mesh dimension */ 11 PetscBool simplex; /* Flag for simplices or tensor cells */ 12 char filename[PETSC_MAX_PATH_LEN]; /* Name of the mesh filename if any */ 13 PetscReal omega; /* Oscillation frequency omega */ 14 PetscInt particlesPerCell; /* The number of partices per cell */ 15 PetscInt numberOfCells; /* Number of cells in mesh */ 16 PetscReal momentTol; /* Tolerance for checking moment conservation */ 17 PetscBool monitor; /* Flag for using the TS monitor */ 18 PetscBool error; /* Flag for printing the error */ 19 PetscInt ostep; /* print the energy at each ostep time steps */ 20 } AppCtx; 21 22 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 23 { 24 PetscErrorCode ierr; 25 26 PetscFunctionBeginUser; 27 options->dim = 2; 28 options->simplex = PETSC_TRUE; 29 options->monitor = PETSC_FALSE; 30 options->error = PETSC_FALSE; 31 options->particlesPerCell = 1; 32 options->numberOfCells = 2; 33 options->momentTol = 100.0*PETSC_MACHINE_EPSILON; 34 options->omega = 64.0; 35 options->ostep = 100; 36 37 CHKERRQ(PetscStrcpy(options->filename, "")); 38 39 ierr = PetscOptionsBegin(comm, "", "Harmonic Oscillator Options", "DMPLEX");CHKERRQ(ierr); 40 CHKERRQ(PetscOptionsInt("-output_step", "Number of time steps between output", "ex4.c", options->ostep, &options->ostep, PETSC_NULL)); 41 CHKERRQ(PetscOptionsInt("-dim", "The topological mesh dimension", "ex4.c", options->dim, &options->dim, NULL)); 42 CHKERRQ(PetscOptionsBool("-monitor", "Flag to use the TS monitor", "ex4.c", options->monitor, &options->monitor, NULL)); 43 CHKERRQ(PetscOptionsBool("-error", "Flag to print the error", "ex4.c", options->error, &options->error, NULL)); 44 CHKERRQ(PetscOptionsBool("-simplex", "The flag for simplices or tensor cells", "ex4.c", options->simplex, &options->simplex, NULL)); 45 CHKERRQ(PetscOptionsString("-mesh", "Name of the mesh filename if any", "ex4.c", options->filename, options->filename, sizeof(options->filename), NULL)); 46 CHKERRQ(PetscOptionsInt("-particles_per_cell", "Number of particles per cell", "ex4.c", options->particlesPerCell, &options->particlesPerCell, NULL)); 47 CHKERRQ(PetscOptionsReal("-omega", "Oscillator frequency", "ex4.c", options->omega, &options->omega, PETSC_NULL)); 48 ierr = PetscOptionsEnd();CHKERRQ(ierr); 49 50 PetscFunctionReturn(0); 51 52 } 53 54 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user) 55 { 56 PetscFunctionBeginUser; 57 CHKERRQ(DMCreate(comm, dm)); 58 CHKERRQ(DMSetType(*dm, DMPLEX)); 59 CHKERRQ(DMSetFromOptions(*dm)); 60 CHKERRQ(DMViewFromOptions(*dm, NULL, "-dm_view")); 61 62 PetscFunctionReturn(0); 63 } 64 65 static PetscErrorCode SetInitialCoordinates(DM dmSw) 66 { 67 DM dm; 68 AppCtx *user; 69 PetscRandom rnd; 70 PetscBool simplex; 71 PetscReal *centroid, *coords, *xi0, *v0, *J, *invJ, detJ; 72 PetscInt dim, d, cStart, cEnd, c, Np, p; 73 74 PetscFunctionBeginUser; 75 CHKERRQ(PetscRandomCreate(PetscObjectComm((PetscObject) dmSw), &rnd)); 76 CHKERRQ(PetscRandomSetInterval(rnd, -1.0, 1.0)); 77 CHKERRQ(PetscRandomSetFromOptions(rnd)); 78 79 CHKERRQ(DMGetApplicationContext(dmSw, &user)); 80 Np = user->particlesPerCell; 81 CHKERRQ(DMSwarmGetCellDM(dmSw, &dm)); 82 CHKERRQ(DMGetDimension(dm, &dim)); 83 CHKERRQ(DMPlexIsSimplex(dm, &simplex)); 84 CHKERRQ(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 85 user->numberOfCells = cEnd - cStart; 86 CHKERRQ(PetscMalloc5(dim, ¢roid, dim, &xi0, dim, &v0, dim*dim, &J, dim*dim, &invJ)); 87 for (d = 0; d < dim; ++d) xi0[d] = -1.0; 88 CHKERRQ(DMSwarmGetField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords)); 89 for (c = cStart; c < cEnd; ++c) { 90 if (Np == 1) { 91 CHKERRQ(DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL)); 92 for (d = 0; d < dim; ++d) coords[c*dim+d] = centroid[d]; 93 } else { 94 CHKERRQ(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ)); /* affine */ 95 for (p = 0; p < Np; ++p) { 96 const PetscInt n = c*Np + p; 97 PetscReal sum = 0.0, refcoords[3]; 98 99 for (d = 0; d < dim; ++d) { 100 CHKERRQ(PetscRandomGetValueReal(rnd, &refcoords[d])); 101 sum += refcoords[d]; 102 } 103 if (simplex && sum > 0.0) for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim)*sum; 104 CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n*dim]); 105 } 106 } 107 } 108 109 CHKERRQ(DMSwarmRestoreField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords)); 110 CHKERRQ(PetscFree5(centroid, xi0, v0, J, invJ)); 111 CHKERRQ(PetscRandomDestroy(&rnd)); 112 PetscFunctionReturn(0); 113 } 114 115 static PetscErrorCode SetInitialConditions(DM dmSw, Vec u) 116 { 117 DM dm; 118 AppCtx *user; 119 PetscReal *coords; 120 PetscScalar *initialConditions; 121 PetscInt dim, cStart, cEnd, c, Np, p; 122 123 PetscFunctionBeginUser; 124 CHKERRQ(DMGetApplicationContext(dmSw, &user)); 125 Np = user->particlesPerCell; 126 CHKERRQ(DMSwarmGetCellDM(dmSw, &dm)); 127 CHKERRQ(DMGetDimension(dm, &dim)); 128 CHKERRQ(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 129 CHKERRQ(DMSwarmGetField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords)); 130 CHKERRQ(VecGetArray(u, &initialConditions)); 131 for (c = cStart; c < cEnd; ++c) { 132 for (p = 0; p < Np; ++p) { 133 const PetscInt n = c*Np + p; 134 initialConditions[n*2+0] = DMPlex_NormD_Internal(dim, &coords[n*dim]); 135 initialConditions[n*2+1] = 0.0; 136 } 137 } 138 CHKERRQ(VecRestoreArray(u, &initialConditions)); 139 CHKERRQ(DMSwarmRestoreField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords)); 140 PetscFunctionReturn(0); 141 } 142 143 static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user) 144 { 145 PetscInt *cellid; 146 PetscInt dim, cStart, cEnd, c, Np, p; 147 148 PetscFunctionBeginUser; 149 CHKERRQ(DMGetDimension(dm, &dim)); 150 CHKERRQ(DMCreate(PetscObjectComm((PetscObject) dm), sw)); 151 CHKERRQ(DMSetType(*sw, DMSWARM)); 152 CHKERRQ(DMSetDimension(*sw, dim)); 153 Np = user->particlesPerCell; 154 CHKERRQ(DMSwarmSetType(*sw, DMSWARM_PIC)); 155 CHKERRQ(DMSwarmSetCellDM(*sw, dm)); 156 CHKERRQ(DMSwarmRegisterPetscDatatypeField(*sw, "kinematics", 2, PETSC_REAL)); 157 CHKERRQ(DMSwarmFinalizeFieldRegister(*sw)); 158 CHKERRQ(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 159 CHKERRQ(DMSwarmSetLocalSizes(*sw, (cEnd - cStart) * Np, 0)); 160 CHKERRQ(DMSetFromOptions(*sw)); 161 CHKERRQ(DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid)); 162 for (c = cStart; c < cEnd; ++c) { 163 for (p = 0; p < Np; ++p) { 164 const PetscInt n = c*Np + p; 165 166 cellid[n] = c; 167 } 168 } 169 CHKERRQ(DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid)); 170 CHKERRQ(PetscObjectSetName((PetscObject) *sw, "Particles")); 171 CHKERRQ(DMViewFromOptions(*sw, NULL, "-sw_view")); 172 PetscFunctionReturn(0); 173 } 174 175 static PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx) 176 { 177 AppCtx *user = (AppCtx *) ctx; 178 const PetscReal omega = user->omega; 179 const PetscScalar *u; 180 MPI_Comm comm; 181 PetscReal dt; 182 PetscInt Np, p; 183 184 PetscFunctionBeginUser; 185 if (step%user->ostep == 0) { 186 CHKERRQ(PetscObjectGetComm((PetscObject) ts, &comm)); 187 if (!step) CHKERRQ(PetscPrintf(comm, "Time Step Part Energy Mod Energy\n")); 188 CHKERRQ(TSGetTimeStep(ts, &dt)); 189 CHKERRQ(VecGetArrayRead(U, &u)); 190 CHKERRQ(VecGetLocalSize(U, &Np)); 191 Np /= 2; 192 for (p = 0; p < Np; ++p) { 193 const PetscReal x = PetscRealPart(u[p*2+0]); 194 const PetscReal v = PetscRealPart(u[p*2+1]); 195 const PetscReal E = 0.5*(v*v + PetscSqr(omega)*x*x); 196 const PetscReal mE = 0.5*(v*v + PetscSqr(omega)*x*x - PetscSqr(omega)*dt*x*v); 197 CHKERRQ(PetscPrintf(comm, "%.6lf %4D %4D %10.4lf %10.4lf\n", t, step, p, (double) E, (double) mE)); 198 } 199 CHKERRQ(VecRestoreArrayRead(U, &u)); 200 } 201 PetscFunctionReturn(0); 202 } 203 204 static PetscErrorCode InitializeSolve(TS ts, Vec u) 205 { 206 DM dm; 207 AppCtx *user; 208 209 PetscFunctionBeginUser; 210 CHKERRQ(TSGetDM(ts, &dm)); 211 CHKERRQ(DMGetApplicationContext(dm, &user)); 212 CHKERRQ(SetInitialCoordinates(dm)); 213 CHKERRQ(SetInitialConditions(dm, u)); 214 PetscFunctionReturn(0); 215 } 216 217 static PetscErrorCode ComputeError(TS ts, Vec U, Vec E) 218 { 219 MPI_Comm comm; 220 DM sdm; 221 AppCtx *user; 222 const PetscScalar *u, *coords; 223 PetscScalar *e; 224 PetscReal t, omega; 225 PetscInt dim, Np, p; 226 227 PetscFunctionBeginUser; 228 CHKERRQ(PetscObjectGetComm((PetscObject) ts, &comm)); 229 CHKERRQ(TSGetDM(ts, &sdm)); 230 CHKERRQ(DMGetApplicationContext(sdm, &user)); 231 omega = user->omega; 232 CHKERRQ(DMGetDimension(sdm, &dim)); 233 CHKERRQ(TSGetSolveTime(ts, &t)); 234 CHKERRQ(VecGetArray(E, &e)); 235 CHKERRQ(VecGetArrayRead(U, &u)); 236 CHKERRQ(VecGetLocalSize(U, &Np)); 237 Np /= 2; 238 CHKERRQ(DMSwarmGetField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords)); 239 for (p = 0; p < Np; ++p) { 240 const PetscReal x = PetscRealPart(u[p*2+0]); 241 const PetscReal v = PetscRealPart(u[p*2+1]); 242 const PetscReal x0 = DMPlex_NormD_Internal(dim, &coords[p*dim]); 243 const PetscReal ex = x0*PetscCosReal(omega*t); 244 const PetscReal ev = -x0*omega*PetscSinReal(omega*t); 245 246 if (user->error) CHKERRQ(PetscPrintf(comm, "p%D error [%.2g %.2g] sol [%.6lf %.6lf] exact [%.6lf %.6lf] energy/exact energy %g / %g\n", p, (double) PetscAbsReal(x-ex), (double) PetscAbsReal(v-ev), (double) x, (double) v, (double) ex, (double) ev, 0.5*(v*v + PetscSqr(omega)*x*x), (double) 0.5*PetscSqr(omega*x0))); 247 e[p*2+0] = x - ex; 248 e[p*2+1] = v - ev; 249 } 250 CHKERRQ(DMSwarmRestoreField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords)); 251 252 CHKERRQ(VecRestoreArrayRead(U, &u)); 253 CHKERRQ(VecRestoreArray(E, &e)); 254 PetscFunctionReturn(0); 255 } 256 257 /*---------------------Create particle RHS Functions--------------------------*/ 258 static PetscErrorCode RHSFunction1(TS ts, PetscReal t, Vec V, Vec Xres, void *ctx) 259 { 260 const PetscScalar *v; 261 PetscScalar *xres; 262 PetscInt Np, p; 263 264 PetscFunctionBeginUser; 265 CHKERRQ(VecGetArray(Xres, &xres)); 266 CHKERRQ(VecGetArrayRead(V, &v)); 267 CHKERRQ(VecGetLocalSize(Xres, &Np)); 268 for (p = 0; p < Np; ++p) { 269 xres[p] = v[p]; 270 } 271 CHKERRQ(VecRestoreArrayRead(V, &v)); 272 CHKERRQ(VecRestoreArray(Xres, &xres)); 273 PetscFunctionReturn(0); 274 } 275 276 static PetscErrorCode RHSFunction2(TS ts, PetscReal t, Vec X, Vec Vres, void *ctx) 277 { 278 AppCtx *user = (AppCtx *)ctx; 279 const PetscScalar *x; 280 PetscScalar *vres; 281 PetscInt Np, p; 282 283 PetscFunctionBeginUser; 284 CHKERRQ(VecGetArray(Vres, &vres)); 285 CHKERRQ(VecGetArrayRead(X, &x)); 286 CHKERRQ(VecGetLocalSize(Vres, &Np)); 287 for (p = 0; p < Np; ++p) { 288 vres[p] = -PetscSqr(user->omega)*x[p]; 289 } 290 CHKERRQ(VecRestoreArrayRead(X, &x)); 291 CHKERRQ(VecRestoreArray(Vres, &vres)); 292 PetscFunctionReturn(0); 293 } 294 295 /*----------------------------------------------------------------------------*/ 296 297 /*--------------------Define RHSFunction, RHSJacobian (Theta)-----------------*/ 298 static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec G, void *ctx) 299 { 300 AppCtx *user = (AppCtx *) ctx; 301 DM dm; 302 const PetscScalar *u; 303 PetscScalar *g; 304 PetscInt Np, p; 305 306 PetscFunctionBeginUser; 307 CHKERRQ(TSGetDM(ts, &dm)); 308 CHKERRQ(VecGetArray(G, &g)); 309 CHKERRQ(VecGetArrayRead(U, &u)); 310 CHKERRQ(VecGetLocalSize(U, &Np)); 311 Np /= 2; 312 for (p = 0; p < Np; ++p) { 313 g[p*2+0] = u[p*2+1]; 314 g[p*2+1] = -PetscSqr(user->omega)*u[p*2+0]; 315 } 316 CHKERRQ(VecRestoreArrayRead(U, &u)); 317 CHKERRQ(VecRestoreArray(G, &g)); 318 PetscFunctionReturn(0); 319 } 320 321 /*Ji = dFi/dxj 322 J= (0 1) 323 (-w^2 0) 324 */ 325 static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U , Mat J, Mat P, void *ctx) 326 { 327 AppCtx *user = (AppCtx *) ctx; 328 PetscInt Np; 329 PetscInt i, m, n; 330 const PetscScalar *u; 331 PetscScalar vals[4] = {0., 1., -PetscSqr(user->omega), 0.}; 332 333 PetscFunctionBeginUser; 334 CHKERRQ(VecGetArrayRead(U, &u)); 335 CHKERRQ(VecGetLocalSize(U, &Np)); 336 Np /= 2; 337 CHKERRQ(MatGetOwnershipRange(J, &m, &n)); 338 for (i = 0; i < Np; ++i) { 339 const PetscInt rows[2] = {2*i, 2*i+1}; 340 CHKERRQ(MatSetValues(J, 2, rows, 2, rows, vals, INSERT_VALUES)); 341 } 342 CHKERRQ(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY)); 343 CHKERRQ(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY)); 344 PetscFunctionReturn(0); 345 346 } 347 348 /*----------------------------------------------------------------------------*/ 349 350 /*----------------Define S, F, G Functions (Discrete Gradients)---------------*/ 351 /* 352 u_t = S * gradF 353 --or-- 354 u_t = S * G 355 356 + Sfunc - constructor for the S matrix from the formulation 357 . Ffunc - functional F from the formulation 358 - Gfunc - constructor for the gradient of F from the formulation 359 */ 360 361 PetscErrorCode Sfunc(TS ts, PetscReal t, Vec U, Mat S, void *ctx) 362 { 363 PetscInt Np; 364 PetscInt i, m, n; 365 const PetscScalar *u; 366 PetscScalar vals[4] = {0., 1., -1, 0.}; 367 368 PetscFunctionBeginUser; 369 CHKERRQ(VecGetArrayRead(U, &u)); 370 CHKERRQ(VecGetLocalSize(U, &Np)); 371 Np /= 2; 372 CHKERRQ(MatGetOwnershipRange(S, &m, &n)); 373 for (i = 0; i < Np; ++i) { 374 const PetscInt rows[2] = {2*i, 2*i+1}; 375 CHKERRQ(MatSetValues(S, 2, rows, 2, rows, vals, INSERT_VALUES)); 376 } 377 CHKERRQ(MatAssemblyBegin(S,MAT_FINAL_ASSEMBLY)); 378 CHKERRQ(MatAssemblyEnd(S,MAT_FINAL_ASSEMBLY)); 379 PetscFunctionReturn(0); 380 381 } 382 383 PetscErrorCode Ffunc(TS ts, PetscReal t, Vec U, PetscScalar *F, void *ctx) 384 { 385 AppCtx *user = (AppCtx *) ctx; 386 DM dm; 387 const PetscScalar *u; 388 PetscInt Np; 389 PetscInt p; 390 391 PetscFunctionBeginUser; 392 CHKERRQ(TSGetDM(ts, &dm)); 393 394 /* Define F */ 395 CHKERRQ(VecGetArrayRead(U, &u)); 396 CHKERRQ(VecGetLocalSize(U, &Np)); 397 Np /= 2; 398 for (p = 0; p < Np; ++p) { 399 *F += 0.5*PetscSqr(user->omega)*PetscSqr(u[p*2+0]) + 0.5*PetscSqr(u[p*2+1]); 400 } 401 CHKERRQ(VecRestoreArrayRead(U, &u)); 402 PetscFunctionReturn(0); 403 } 404 405 PetscErrorCode gradFfunc(TS ts, PetscReal t, Vec U, Vec gradF, void *ctx) 406 { 407 AppCtx *user = (AppCtx *) ctx; 408 DM dm; 409 const PetscScalar *u; 410 PetscScalar *g; 411 PetscInt Np; 412 PetscInt p; 413 414 PetscFunctionBeginUser; 415 CHKERRQ(TSGetDM(ts, &dm)); 416 CHKERRQ(VecGetArrayRead(U, &u)); 417 CHKERRQ(VecGetLocalSize(U, &Np)); 418 Np /= 2; 419 /*Define gradF*/ 420 CHKERRQ(VecGetArray(gradF, &g)); 421 for (p = 0; p < Np; ++p) { 422 g[p*2+0] = PetscSqr(user->omega)*u[p*2+0]; /*dF/dx*/ 423 g[p*2+1] = u[p*2+1]; /*dF/dv*/ 424 } 425 CHKERRQ(VecRestoreArrayRead(U, &u)); 426 CHKERRQ(VecRestoreArray(gradF, &g)); 427 PetscFunctionReturn(0); 428 } 429 430 /*----------------------------------------------------------------------------*/ 431 432 int main(int argc,char **argv) 433 { 434 TS ts; /* nonlinear solver */ 435 DM dm, sw; /* Mesh and particle managers */ 436 Vec u; /* swarm vector */ 437 PetscInt n; /* swarm vector size */ 438 IS is1, is2; 439 MPI_Comm comm; 440 Mat J; /* Jacobian matrix */ 441 AppCtx user; 442 PetscErrorCode ierr; 443 444 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 445 Initialize program 446 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 447 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 448 comm = PETSC_COMM_WORLD; 449 CHKERRQ(ProcessOptions(comm, &user)); 450 451 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 452 Create Particle-Mesh 453 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 454 CHKERRQ(CreateMesh(comm, &dm, &user)); 455 CHKERRQ(CreateParticles(dm, &sw, &user)); 456 CHKERRQ(DMSetApplicationContext(sw, &user)); 457 458 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 459 Setup timestepping solver 460 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 461 CHKERRQ(TSCreate(comm, &ts)); 462 CHKERRQ(TSSetProblemType(ts,TS_NONLINEAR)); 463 464 CHKERRQ(TSSetDM(ts, sw)); 465 CHKERRQ(TSSetMaxTime(ts, 0.1)); 466 CHKERRQ(TSSetTimeStep(ts, 0.00001)); 467 CHKERRQ(TSSetMaxSteps(ts, 100)); 468 CHKERRQ(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 469 if (user.monitor) CHKERRQ(TSMonitorSet(ts, Monitor, &user, NULL)); 470 471 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 472 Prepare to solve 473 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 474 CHKERRQ(DMSwarmCreateGlobalVectorFromField(sw, "kinematics", &u)); 475 CHKERRQ(VecGetLocalSize(u, &n)); 476 CHKERRQ(TSSetFromOptions(ts)); 477 CHKERRQ(TSSetComputeInitialCondition(ts, InitializeSolve)); 478 CHKERRQ(TSSetComputeExactError(ts, ComputeError)); 479 480 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 481 Define function F(U, Udot , x , t) = G(U, x, t) 482 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 483 484 /* - - - - - - - Basic Symplectic - - - - - - - - - - - - - - - - - - - - - -*/ 485 CHKERRQ(ISCreateStride(comm, n/2, 0, 2, &is1)); 486 CHKERRQ(ISCreateStride(comm, n/2, 1, 2, &is2)); 487 CHKERRQ(TSRHSSplitSetIS(ts, "position", is1)); 488 CHKERRQ(TSRHSSplitSetIS(ts, "momentum", is2)); 489 CHKERRQ(ISDestroy(&is1)); 490 CHKERRQ(ISDestroy(&is2)); 491 CHKERRQ(TSRHSSplitSetRHSFunction(ts, "position", NULL, RHSFunction1, &user)); 492 CHKERRQ(TSRHSSplitSetRHSFunction(ts, "momentum", NULL, RHSFunction2, &user)); 493 494 /* - - - - - - - Theta (Implicit Midpoint) - - - - - - - - - - - - - - - - -*/ 495 496 CHKERRQ(TSSetRHSFunction(ts, NULL, RHSFunction, &user)); 497 498 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&J)); 499 CHKERRQ(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n)); 500 CHKERRQ(MatSetFromOptions(J)); 501 CHKERRQ(MatSetUp(J)); 502 CHKERRQ(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY)); 503 CHKERRQ(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY)); 504 CHKERRQ(TSSetRHSJacobian(ts,J,J,RHSJacobian,&user)); 505 506 /* - - - - - - - Discrete Gradients - - - - - - - - - - - - - - - - - - - - */ 507 CHKERRQ(TSDiscGradSetFormulation(ts, Sfunc, Ffunc, gradFfunc, &user)); 508 509 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 510 Solve 511 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 512 CHKERRQ(TSComputeInitialCondition(ts, u)); 513 CHKERRQ(TSSolve(ts, u)); 514 515 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 516 Clean up workspace 517 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 518 CHKERRQ(MatDestroy(&J)); 519 CHKERRQ(DMSwarmDestroyGlobalVectorFromField(sw, "kinematics", &u)); 520 CHKERRQ(TSDestroy(&ts)); 521 CHKERRQ(DMDestroy(&sw)); 522 CHKERRQ(DMDestroy(&dm)); 523 ierr = PetscFinalize(); 524 return ierr; 525 } 526 527 /*TEST 528 529 build: 530 requires: triangle !single !complex 531 test: 532 suffix: 1 533 args: -dm_plex_box_faces 1,1 -ts_type basicsymplectic -ts_basicsymplectic_type 1 -ts_convergence_estimate -convest_num_refine 2 -dm_view -monitor -output_step 50 -error 534 535 test: 536 suffix: 2 537 args: -dm_plex_box_faces 1,1 -ts_type basicsymplectic -ts_basicsymplectic_type 2 -ts_convergence_estimate -convest_num_refine 2 -dm_view -monitor -output_step 50 -error 538 539 test: 540 suffix: 3 541 args: -dm_plex_box_faces 1,1 -ts_type basicsymplectic -ts_basicsymplectic_type 3 -ts_convergence_estimate -convest_num_refine 2 -ts_dt 0.0001 -dm_view -monitor -output_step 50 -error 542 543 test: 544 suffix: 4 545 args: -dm_plex_box_faces 1,1 -ts_type basicsymplectic -ts_basicsymplectic_type 4 -ts_convergence_estimate -convest_num_refine 2 -ts_dt 0.0001 -dm_view -monitor -output_step 50 -error 546 547 test: 548 suffix: 5 549 args: -dm_plex_box_faces 1,1 -ts_type theta -ts_theta_theta 0.5 -monitor -output_step 50 -error -ts_convergence_estimate -convest_num_refine 2 -dm_view 550 551 test: 552 suffix: 6 553 args: -dm_plex_box_faces 1,1 -ts_type discgrad -monitor -output_step 50 -error -ts_convergence_estimate -convest_num_refine 2 -dm_view 554 555 test: 556 suffix: 7 557 args: -dm_plex_box_faces 1,1 -ts_type discgrad -ts_discgrad_gonzalez -monitor -output_step 50 -error -ts_convergence_estimate -convest_num_refine 2 -dm_view 558 559 TEST*/ 560