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