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 PetscCall(PetscStrcpy(options->filename, "")); 38 39 ierr = PetscOptionsBegin(comm, "", "Harmonic Oscillator Options", "DMPLEX");PetscCall(ierr); 40 PetscCall(PetscOptionsInt("-output_step", "Number of time steps between output", "ex4.c", options->ostep, &options->ostep, PETSC_NULL)); 41 PetscCall(PetscOptionsInt("-dim", "The topological mesh dimension", "ex4.c", options->dim, &options->dim, NULL)); 42 PetscCall(PetscOptionsBool("-monitor", "Flag to use the TS monitor", "ex4.c", options->monitor, &options->monitor, NULL)); 43 PetscCall(PetscOptionsBool("-error", "Flag to print the error", "ex4.c", options->error, &options->error, NULL)); 44 PetscCall(PetscOptionsBool("-simplex", "The flag for simplices or tensor cells", "ex4.c", options->simplex, &options->simplex, NULL)); 45 PetscCall(PetscOptionsString("-mesh", "Name of the mesh filename if any", "ex4.c", options->filename, options->filename, sizeof(options->filename), NULL)); 46 PetscCall(PetscOptionsInt("-particles_per_cell", "Number of particles per cell", "ex4.c", options->particlesPerCell, &options->particlesPerCell, NULL)); 47 PetscCall(PetscOptionsReal("-omega", "Oscillator frequency", "ex4.c", options->omega, &options->omega, PETSC_NULL)); 48 ierr = PetscOptionsEnd();PetscCall(ierr); 49 50 PetscFunctionReturn(0); 51 52 } 53 54 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user) 55 { 56 PetscFunctionBeginUser; 57 PetscCall(DMCreate(comm, dm)); 58 PetscCall(DMSetType(*dm, DMPLEX)); 59 PetscCall(DMSetFromOptions(*dm)); 60 PetscCall(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 PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject) dmSw), &rnd)); 76 PetscCall(PetscRandomSetInterval(rnd, -1.0, 1.0)); 77 PetscCall(PetscRandomSetFromOptions(rnd)); 78 79 PetscCall(DMGetApplicationContext(dmSw, &user)); 80 Np = user->particlesPerCell; 81 PetscCall(DMSwarmGetCellDM(dmSw, &dm)); 82 PetscCall(DMGetDimension(dm, &dim)); 83 PetscCall(DMPlexIsSimplex(dm, &simplex)); 84 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 85 user->numberOfCells = cEnd - cStart; 86 PetscCall(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 PetscCall(DMSwarmGetField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords)); 89 for (c = cStart; c < cEnd; ++c) { 90 if (Np == 1) { 91 PetscCall(DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL)); 92 for (d = 0; d < dim; ++d) coords[c*dim+d] = centroid[d]; 93 } else { 94 PetscCall(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 PetscCall(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 PetscCall(DMSwarmRestoreField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords)); 110 PetscCall(PetscFree5(centroid, xi0, v0, J, invJ)); 111 PetscCall(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 PetscCall(DMGetApplicationContext(dmSw, &user)); 125 Np = user->particlesPerCell; 126 PetscCall(DMSwarmGetCellDM(dmSw, &dm)); 127 PetscCall(DMGetDimension(dm, &dim)); 128 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 129 PetscCall(DMSwarmGetField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords)); 130 PetscCall(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 PetscCall(VecRestoreArray(u, &initialConditions)); 139 PetscCall(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 PetscCall(DMGetDimension(dm, &dim)); 150 PetscCall(DMCreate(PetscObjectComm((PetscObject) dm), sw)); 151 PetscCall(DMSetType(*sw, DMSWARM)); 152 PetscCall(DMSetDimension(*sw, dim)); 153 Np = user->particlesPerCell; 154 PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC)); 155 PetscCall(DMSwarmSetCellDM(*sw, dm)); 156 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "kinematics", 2, PETSC_REAL)); 157 PetscCall(DMSwarmFinalizeFieldRegister(*sw)); 158 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 159 PetscCall(DMSwarmSetLocalSizes(*sw, (cEnd - cStart) * Np, 0)); 160 PetscCall(DMSetFromOptions(*sw)); 161 PetscCall(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 PetscCall(DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid)); 170 PetscCall(PetscObjectSetName((PetscObject) *sw, "Particles")); 171 PetscCall(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 PetscCall(PetscObjectGetComm((PetscObject) ts, &comm)); 187 if (!step) PetscCall(PetscPrintf(comm, "Time Step Part Energy Mod Energy\n")); 188 PetscCall(TSGetTimeStep(ts, &dt)); 189 PetscCall(VecGetArrayRead(U, &u)); 190 PetscCall(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 PetscCall(PetscPrintf(comm, "%.6lf %4D %4D %10.4lf %10.4lf\n", t, step, p, (double) E, (double) mE)); 198 } 199 PetscCall(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 PetscCall(TSGetDM(ts, &dm)); 211 PetscCall(DMGetApplicationContext(dm, &user)); 212 PetscCall(SetInitialCoordinates(dm)); 213 PetscCall(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 PetscCall(PetscObjectGetComm((PetscObject) ts, &comm)); 229 PetscCall(TSGetDM(ts, &sdm)); 230 PetscCall(DMGetApplicationContext(sdm, &user)); 231 omega = user->omega; 232 PetscCall(DMGetDimension(sdm, &dim)); 233 PetscCall(TSGetSolveTime(ts, &t)); 234 PetscCall(VecGetArray(E, &e)); 235 PetscCall(VecGetArrayRead(U, &u)); 236 PetscCall(VecGetLocalSize(U, &Np)); 237 Np /= 2; 238 PetscCall(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) PetscCall(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 PetscCall(DMSwarmRestoreField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords)); 251 252 PetscCall(VecRestoreArrayRead(U, &u)); 253 PetscCall(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 PetscCall(VecGetArray(Xres, &xres)); 266 PetscCall(VecGetArrayRead(V, &v)); 267 PetscCall(VecGetLocalSize(Xres, &Np)); 268 for (p = 0; p < Np; ++p) { 269 xres[p] = v[p]; 270 } 271 PetscCall(VecRestoreArrayRead(V, &v)); 272 PetscCall(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 PetscCall(VecGetArray(Vres, &vres)); 285 PetscCall(VecGetArrayRead(X, &x)); 286 PetscCall(VecGetLocalSize(Vres, &Np)); 287 for (p = 0; p < Np; ++p) { 288 vres[p] = -PetscSqr(user->omega)*x[p]; 289 } 290 PetscCall(VecRestoreArrayRead(X, &x)); 291 PetscCall(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 PetscCall(TSGetDM(ts, &dm)); 308 PetscCall(VecGetArray(G, &g)); 309 PetscCall(VecGetArrayRead(U, &u)); 310 PetscCall(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 PetscCall(VecRestoreArrayRead(U, &u)); 317 PetscCall(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 PetscCall(VecGetArrayRead(U, &u)); 335 PetscCall(VecGetLocalSize(U, &Np)); 336 Np /= 2; 337 PetscCall(MatGetOwnershipRange(J, &m, &n)); 338 for (i = 0; i < Np; ++i) { 339 const PetscInt rows[2] = {2*i, 2*i+1}; 340 PetscCall(MatSetValues(J, 2, rows, 2, rows, vals, INSERT_VALUES)); 341 } 342 PetscCall(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY)); 343 PetscCall(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 PetscCall(VecGetArrayRead(U, &u)); 370 PetscCall(VecGetLocalSize(U, &Np)); 371 Np /= 2; 372 PetscCall(MatGetOwnershipRange(S, &m, &n)); 373 for (i = 0; i < Np; ++i) { 374 const PetscInt rows[2] = {2*i, 2*i+1}; 375 PetscCall(MatSetValues(S, 2, rows, 2, rows, vals, INSERT_VALUES)); 376 } 377 PetscCall(MatAssemblyBegin(S,MAT_FINAL_ASSEMBLY)); 378 PetscCall(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 PetscCall(TSGetDM(ts, &dm)); 393 394 /* Define F */ 395 PetscCall(VecGetArrayRead(U, &u)); 396 PetscCall(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 PetscCall(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 PetscCall(TSGetDM(ts, &dm)); 416 PetscCall(VecGetArrayRead(U, &u)); 417 PetscCall(VecGetLocalSize(U, &Np)); 418 Np /= 2; 419 /*Define gradF*/ 420 PetscCall(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 PetscCall(VecRestoreArrayRead(U, &u)); 426 PetscCall(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 443 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 444 Initialize program 445 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 446 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 447 comm = PETSC_COMM_WORLD; 448 PetscCall(ProcessOptions(comm, &user)); 449 450 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 451 Create Particle-Mesh 452 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 453 PetscCall(CreateMesh(comm, &dm, &user)); 454 PetscCall(CreateParticles(dm, &sw, &user)); 455 PetscCall(DMSetApplicationContext(sw, &user)); 456 457 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 458 Setup timestepping solver 459 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 460 PetscCall(TSCreate(comm, &ts)); 461 PetscCall(TSSetProblemType(ts,TS_NONLINEAR)); 462 463 PetscCall(TSSetDM(ts, sw)); 464 PetscCall(TSSetMaxTime(ts, 0.1)); 465 PetscCall(TSSetTimeStep(ts, 0.00001)); 466 PetscCall(TSSetMaxSteps(ts, 100)); 467 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 468 if (user.monitor) PetscCall(TSMonitorSet(ts, Monitor, &user, NULL)); 469 470 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 471 Prepare to solve 472 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 473 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "kinematics", &u)); 474 PetscCall(VecGetLocalSize(u, &n)); 475 PetscCall(TSSetFromOptions(ts)); 476 PetscCall(TSSetComputeInitialCondition(ts, InitializeSolve)); 477 PetscCall(TSSetComputeExactError(ts, ComputeError)); 478 479 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 480 Define function F(U, Udot , x , t) = G(U, x, t) 481 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 482 483 /* - - - - - - - Basic Symplectic - - - - - - - - - - - - - - - - - - - - - -*/ 484 PetscCall(ISCreateStride(comm, n/2, 0, 2, &is1)); 485 PetscCall(ISCreateStride(comm, n/2, 1, 2, &is2)); 486 PetscCall(TSRHSSplitSetIS(ts, "position", is1)); 487 PetscCall(TSRHSSplitSetIS(ts, "momentum", is2)); 488 PetscCall(ISDestroy(&is1)); 489 PetscCall(ISDestroy(&is2)); 490 PetscCall(TSRHSSplitSetRHSFunction(ts, "position", NULL, RHSFunction1, &user)); 491 PetscCall(TSRHSSplitSetRHSFunction(ts, "momentum", NULL, RHSFunction2, &user)); 492 493 /* - - - - - - - Theta (Implicit Midpoint) - - - - - - - - - - - - - - - - -*/ 494 495 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, &user)); 496 497 PetscCall(MatCreate(PETSC_COMM_WORLD,&J)); 498 PetscCall(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n)); 499 PetscCall(MatSetFromOptions(J)); 500 PetscCall(MatSetUp(J)); 501 PetscCall(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY)); 502 PetscCall(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY)); 503 PetscCall(TSSetRHSJacobian(ts,J,J,RHSJacobian,&user)); 504 505 /* - - - - - - - Discrete Gradients - - - - - - - - - - - - - - - - - - - - */ 506 PetscCall(TSDiscGradSetFormulation(ts, Sfunc, Ffunc, gradFfunc, &user)); 507 508 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 509 Solve 510 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 511 PetscCall(TSComputeInitialCondition(ts, u)); 512 PetscCall(TSSolve(ts, u)); 513 514 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 515 Clean up workspace 516 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 517 PetscCall(MatDestroy(&J)); 518 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "kinematics", &u)); 519 PetscCall(TSDestroy(&ts)); 520 PetscCall(DMDestroy(&sw)); 521 PetscCall(DMDestroy(&dm)); 522 PetscCall(PetscFinalize()); 523 return 0; 524 } 525 526 /*TEST 527 528 build: 529 requires: triangle !single !complex 530 test: 531 suffix: 1 532 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 533 534 test: 535 suffix: 2 536 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 537 538 test: 539 suffix: 3 540 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 541 542 test: 543 suffix: 4 544 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 545 546 test: 547 suffix: 5 548 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 549 550 test: 551 suffix: 6 552 args: -dm_plex_box_faces 1,1 -ts_type discgrad -monitor -output_step 50 -error -ts_convergence_estimate -convest_num_refine 2 -dm_view 553 554 test: 555 suffix: 7 556 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 557 558 TEST*/ 559