1 static char help[] = "Vlasov example of many particles orbiting around a several massive points.\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 PetscInt particlesPerCircle; /* The number of partices per circle */ 14 PetscReal momentTol; /* Tolerance for checking moment conservation */ 15 PetscBool monitor; /* Flag for using the TS monitor */ 16 PetscBool error; /* Flag for printing the error */ 17 PetscInt ostep; /* print the energy at each ostep time steps */ 18 PetscReal center[6]; /* Centers of the two orbits */ 19 PetscReal radius[2]; /* Radii of the two orbits */ 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->particlesPerCircle = 1; 32 options->momentTol = 100.0*PETSC_MACHINE_EPSILON; 33 options->ostep = 100; 34 35 ierr = PetscStrcpy(options->filename, "");CHKERRQ(ierr); 36 37 ierr = PetscOptionsBegin(comm, "", "Vlasov Options", "DMPLEX");CHKERRQ(ierr); 38 ierr = PetscOptionsInt("-output_step", "Number of time steps between output", "ex4.c", options->ostep, &options->ostep, PETSC_NULL);CHKERRQ(ierr); 39 ierr = PetscOptionsInt("-dim", "The topological mesh dimension", "ex4.c", options->dim, &options->dim, NULL);CHKERRQ(ierr); 40 ierr = PetscOptionsBool("-monitor", "Flag to use the TS monitor", "ex4.c", options->monitor, &options->monitor, NULL);CHKERRQ(ierr); 41 ierr = PetscOptionsBool("-error", "Flag to print the error", "ex4.c", options->error, &options->error, NULL);CHKERRQ(ierr); 42 ierr = PetscOptionsBool("-simplex", "The flag for simplices or tensor cells", "ex4.c", options->simplex, &options->simplex, NULL);CHKERRQ(ierr); 43 ierr = PetscOptionsString("-mesh", "Name of the mesh filename if any", "ex4.c", options->filename, options->filename, PETSC_MAX_PATH_LEN, NULL);CHKERRQ(ierr); 44 ierr = PetscOptionsInt("-particles_per_circle", "Number of particles per circle", "ex4.c", options->particlesPerCircle, &options->particlesPerCircle, NULL);CHKERRQ(ierr); 45 ierr = PetscOptionsEnd();CHKERRQ(ierr); 46 47 PetscFunctionReturn(0); 48 } 49 50 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user) 51 { 52 PetscBool flg; 53 PetscErrorCode ierr; 54 55 PetscFunctionBeginUser; 56 ierr = PetscStrcmp(user->filename, "", &flg);CHKERRQ(ierr); 57 if (flg) { 58 ierr = DMPlexCreateBoxMesh(comm, user->dim, user->simplex, NULL, NULL, NULL, NULL, PETSC_TRUE, dm);CHKERRQ(ierr); 59 } else { 60 ierr = DMPlexCreateFromFile(comm, user->filename, PETSC_TRUE, dm);CHKERRQ(ierr); 61 ierr = DMGetDimension(*dm, &user->dim);CHKERRQ(ierr); 62 } 63 { 64 DM distributedMesh = NULL; 65 66 ierr = DMPlexDistribute(*dm, 0, NULL, &distributedMesh);CHKERRQ(ierr); 67 if (distributedMesh) { 68 ierr = DMDestroy(dm);CHKERRQ(ierr); 69 *dm = distributedMesh; 70 } 71 } 72 ierr = DMLocalizeCoordinates(*dm);CHKERRQ(ierr); /* needed for periodic */ 73 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 74 ierr = PetscObjectSetName((PetscObject) *dm, "Mesh");CHKERRQ(ierr); 75 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 76 PetscFunctionReturn(0); 77 } 78 79 static PetscErrorCode orbit(AppCtx *ctx, PetscInt c, PetscInt p, PetscReal t, PetscReal x[], PetscReal v[]) 80 { 81 const PetscInt Np = ctx->particlesPerCircle; 82 const PetscReal r = ctx->radius[c]; 83 const PetscReal omega = PetscSqrtReal(1000./r)/r; 84 const PetscReal t0 = (2.*PETSC_PI*p)/(Np*omega); 85 const PetscInt dim = 2; 86 87 PetscFunctionBeginUser; 88 if (x) { 89 x[0] = r*PetscCosReal(omega*(t + t0)) + ctx->center[c*dim + 0]; 90 x[1] = r*PetscSinReal(omega*(t + t0)) + ctx->center[c*dim + 1]; 91 } 92 if (v) { 93 v[0] = -r*omega*PetscSinReal(omega*(t + t0)); 94 v[1] = r*omega*PetscCosReal(omega*(t + t0)); 95 } 96 PetscFunctionReturn(0); 97 } 98 99 static PetscErrorCode force(AppCtx *ctx, PetscInt c, const PetscReal x[], PetscReal force[]) 100 { 101 const PetscReal r = ctx->radius[c]; 102 const PetscReal omega = PetscSqrtReal(1000./r)/r; 103 const PetscInt dim = 2; 104 PetscInt d; 105 106 PetscFunctionBeginUser; 107 for (d = 0; d < dim; ++d) force[d] = -PetscSqr(omega)*(x[d] - ctx->center[c*dim + d]); 108 PetscFunctionReturn(0); 109 } 110 111 static PetscReal energy(AppCtx *ctx, PetscInt c) 112 { 113 const PetscReal r = ctx->radius[c]; 114 const PetscReal omega = PetscSqrtReal(1000./r)/r; 115 116 return 0.5 * omega * r; 117 } 118 119 static PetscErrorCode SetInitialCoordinates(DM dmSw) 120 { 121 DM dm; 122 AppCtx *ctx; 123 Vec coordinates; 124 PetscSF cellSF = NULL; 125 PetscReal *coords; 126 PetscInt *cellid; 127 const PetscInt *found; 128 const PetscSFNode *cells; 129 PetscInt dim, d, c, Np, p; 130 PetscMPIInt rank; 131 PetscErrorCode ierr; 132 133 PetscFunctionBeginUser; 134 ierr = DMGetApplicationContext(dmSw, (void **) &ctx);CHKERRQ(ierr); 135 Np = ctx->particlesPerCircle; 136 ierr = DMSwarmGetCellDM(dmSw, &dm);CHKERRQ(ierr); 137 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 138 ierr = DMSwarmGetField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);CHKERRQ(ierr); 139 for (c = 0; c < 2; ++c) { 140 for (d = 0; d < dim; ++d) ctx->center[c*dim+d] = (!c && !d) ? 3.0 : 0.0; 141 ctx->radius[c] = 3.*c+1.; 142 for (p = 0; p < Np; ++p) { 143 const PetscInt n = c*Np + p; 144 145 ierr = orbit(ctx, c, p, 0.0, &coords[n*dim], NULL);CHKERRQ(ierr); 146 } 147 } 148 ierr = DMSwarmRestoreField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);CHKERRQ(ierr); 149 150 ierr = DMSwarmCreateGlobalVectorFromField(dmSw, DMSwarmPICField_coor, &coordinates);CHKERRQ(ierr); 151 ierr = DMLocatePoints(dm, coordinates, DM_POINTLOCATION_NONE, &cellSF);CHKERRQ(ierr); 152 ierr = DMSwarmDestroyGlobalVectorFromField(dmSw, DMSwarmPICField_coor, &coordinates);CHKERRQ(ierr); 153 154 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dmSw), &rank);CHKERRQ(ierr); 155 ierr = DMSwarmGetField(dmSw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);CHKERRQ(ierr); 156 ierr = PetscSFGetGraph(cellSF, NULL, &Np, &found, &cells);CHKERRQ(ierr); 157 for (p = 0; p < Np; ++p) { 158 const PetscInt part = found ? found[p] : p; 159 160 if (cells[p].rank != rank) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %D not found in the mesh", part); 161 cellid[part] = cells[p].index; 162 } 163 ierr = DMSwarmRestoreField(dmSw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);CHKERRQ(ierr); 164 ierr = PetscSFDestroy(&cellSF);CHKERRQ(ierr); 165 PetscFunctionReturn(0); 166 } 167 168 static PetscErrorCode SetInitialConditions(DM dmSw, Vec u) 169 { 170 DM dm; 171 AppCtx *ctx; 172 PetscScalar *initialConditions; 173 PetscInt dim, cStart, cEnd, c, Np, p; 174 PetscErrorCode ierr; 175 176 PetscFunctionBeginUser; 177 ierr = DMGetApplicationContext(dmSw, (void **) &ctx);CHKERRQ(ierr); 178 Np = ctx->particlesPerCircle; 179 ierr = DMSwarmGetCellDM(dmSw, &dm);CHKERRQ(ierr); 180 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 181 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 182 ierr = VecGetArray(u, &initialConditions);CHKERRQ(ierr); 183 for (c = 0; c < 2; ++c) { 184 for (p = 0; p < Np; ++p) { 185 const PetscInt n = c*Np + p; 186 187 ierr = orbit(ctx, c, p, 0.0, &initialConditions[(n*2 + 0)*dim], &initialConditions[(n*2 + 1)*dim]);CHKERRQ(ierr); 188 } 189 } 190 ierr = VecRestoreArray(u, &initialConditions);CHKERRQ(ierr); 191 PetscFunctionReturn(0); 192 } 193 194 static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user) 195 { 196 PetscInt dim, Np = user->particlesPerCircle; 197 PetscErrorCode ierr; 198 199 PetscFunctionBeginUser; 200 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 201 ierr = DMCreate(PetscObjectComm((PetscObject) dm), sw);CHKERRQ(ierr); 202 ierr = DMSetType(*sw, DMSWARM);CHKERRQ(ierr); 203 ierr = DMSetDimension(*sw, dim);CHKERRQ(ierr); 204 205 ierr = DMSwarmSetType(*sw, DMSWARM_PIC);CHKERRQ(ierr); 206 ierr = DMSwarmSetCellDM(*sw, dm);CHKERRQ(ierr); 207 ierr = DMSwarmRegisterPetscDatatypeField(*sw, "kinematics", 2*dim, PETSC_REAL);CHKERRQ(ierr); 208 ierr = DMSwarmFinalizeFieldRegister(*sw);CHKERRQ(ierr); 209 ierr = DMSwarmSetLocalSizes(*sw, 2 * Np, 0);CHKERRQ(ierr); 210 ierr = DMSetFromOptions(*sw);CHKERRQ(ierr); 211 ierr = PetscObjectSetName((PetscObject) *sw, "Particles");CHKERRQ(ierr); 212 ierr = DMViewFromOptions(*sw, NULL, "-sw_view");CHKERRQ(ierr); 213 PetscFunctionReturn(0); 214 } 215 216 /* Create particle RHS Functions for gravity with G = 1 for simplification */ 217 static PetscErrorCode RHSFunction1(TS ts, PetscReal t, Vec V, Vec Xres, void *ctx) 218 { 219 DM dm; 220 const PetscScalar *v; 221 PetscScalar *xres; 222 PetscInt Np, p, dim, d; 223 PetscErrorCode ierr; 224 225 PetscFunctionBeginUser; 226 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 227 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 228 ierr = VecGetLocalSize(Xres, &Np);CHKERRQ(ierr); 229 Np /= dim; 230 ierr = VecGetArray(Xres, &xres);CHKERRQ(ierr); 231 ierr = VecGetArrayRead(V, &v);CHKERRQ(ierr); 232 for (p = 0; p < Np; ++p) { 233 for(d = 0; d < dim; ++d) { 234 xres[p*dim+d] = v[p*dim+d]; 235 } 236 } 237 ierr = VecRestoreArrayRead(V,& v);CHKERRQ(ierr); 238 ierr = VecRestoreArray(Xres, &xres);CHKERRQ(ierr); 239 PetscFunctionReturn(0); 240 } 241 242 static PetscErrorCode RHSFunction2(TS ts, PetscReal t, Vec X, Vec Vres, void *user) 243 { 244 AppCtx *ctx = (AppCtx *) user; 245 DM dm; 246 const PetscScalar *x; 247 PetscScalar *vres; 248 PetscInt Np, p, dim; 249 PetscErrorCode ierr; 250 251 252 PetscFunctionBeginUser; 253 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 254 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 255 ierr = VecGetLocalSize(Vres, &Np);CHKERRQ(ierr); 256 Np /= dim; 257 ierr = VecGetArray(Vres, &vres);CHKERRQ(ierr); 258 ierr = VecGetArrayRead(X, &x);CHKERRQ(ierr); 259 for (p = 0; p < Np; ++p) { 260 const PetscInt c = p / ctx->particlesPerCircle; 261 262 ierr = force(ctx, c, &x[p*dim], &vres[p*dim]);CHKERRQ(ierr); 263 } 264 ierr = VecRestoreArray(Vres, &vres);CHKERRQ(ierr); 265 ierr = VecRestoreArrayRead(X, &x);CHKERRQ(ierr); 266 PetscFunctionReturn(0); 267 } 268 269 static PetscErrorCode RHSFunctionParticles(TS ts, PetscReal t , Vec U, Vec R, void *user) 270 { 271 AppCtx *ctx = (AppCtx *) user; 272 DM dm; 273 const PetscScalar *u; 274 PetscScalar *r; 275 PetscInt Np, p, dim, d; 276 PetscErrorCode ierr; 277 278 PetscFunctionBeginUser; 279 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 280 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 281 ierr = VecGetLocalSize(U, &Np);CHKERRQ(ierr); 282 Np /= 2*dim; 283 ierr = VecGetArrayRead(U, &u);CHKERRQ(ierr); 284 ierr = VecGetArray(R, &r);CHKERRQ(ierr); 285 for (p = 0; p < Np; ++p) { 286 const PetscInt c = p / ctx->particlesPerCircle; 287 288 for (d = 0; d < dim; ++d) r[(p*2 + 0)*dim + d] = u[(p*2 + 1)*dim + d]; 289 ierr = force(ctx, c, &u[(p*2 + 0)*dim], &r[(p*2 + 1)*dim]);CHKERRQ(ierr); 290 } 291 ierr = VecRestoreArrayRead(U, &u);CHKERRQ(ierr); 292 ierr = VecRestoreArray(R, &r);CHKERRQ(ierr); 293 PetscFunctionReturn(0); 294 } 295 296 static PetscErrorCode InitializeSolve(TS ts, Vec u) 297 { 298 DM dm; 299 AppCtx *user; 300 PetscErrorCode ierr; 301 302 PetscFunctionBeginUser; 303 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 304 ierr = DMGetApplicationContext(dm, (void **) &user);CHKERRQ(ierr); 305 ierr = SetInitialCoordinates(dm);CHKERRQ(ierr); 306 ierr = SetInitialConditions(dm, u);CHKERRQ(ierr); 307 PetscFunctionReturn(0); 308 } 309 310 static PetscErrorCode ComputeError(TS ts, Vec U, Vec E) 311 { 312 MPI_Comm comm; 313 DM sdm; 314 AppCtx *ctx; 315 const PetscScalar *u, *coords; 316 PetscScalar *e; 317 PetscReal t; 318 PetscInt dim, Np, p, c; 319 PetscErrorCode ierr; 320 321 PetscFunctionBeginUser; 322 ierr = PetscObjectGetComm((PetscObject) ts, &comm);CHKERRQ(ierr); 323 ierr = TSGetDM(ts, &sdm);CHKERRQ(ierr); 324 ierr = DMGetApplicationContext(sdm, (void **) &ctx);CHKERRQ(ierr); 325 ierr = DMGetDimension(sdm, &dim);CHKERRQ(ierr); 326 ierr = TSGetSolveTime(ts, &t);CHKERRQ(ierr); 327 ierr = VecGetArray(E, &e);CHKERRQ(ierr); 328 ierr = VecGetArrayRead(U, &u);CHKERRQ(ierr); 329 ierr = VecGetLocalSize(U, &Np);CHKERRQ(ierr); 330 Np /= 2*dim*2; 331 ierr = DMSwarmGetField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);CHKERRQ(ierr); 332 for (c = 0; c < 2; ++c) { 333 for (p = 0; p < Np; ++p) { 334 const PetscInt n = c*Np + p; 335 const PetscScalar *x = &u[(n*2+0)*dim]; 336 const PetscScalar *v = &u[(n*2+1)*dim]; 337 PetscReal xe[3], ve[3]; 338 PetscInt d; 339 340 ierr = orbit(ctx, c, p, t, xe, ve);CHKERRQ(ierr); 341 for (d = 0; d < dim; ++d) { 342 e[(p*2+0)*dim+d] = x[d] - xe[d]; 343 e[(p*2+1)*dim+d] = v[d] - ve[d]; 344 } 345 if (ctx->error) {ierr = PetscPrintf(comm, "p%D error [%.2g %.2g] sol [(%.6lf %.6lf) (%.6lf %.6lf)] exact [(%.6lf %.6lf) (%.6lf %.6lf)] energy/exact energy %g / %g\n", p, (double) DMPlex_NormD_Internal(dim, &e[(p*2+0)*dim]), (double) DMPlex_NormD_Internal(dim, &e[(p*2+1)*dim]), (double) x[0], (double) x[1], (double) v[0], (double) v[1], (double) xe[0], (double) xe[1], (double) ve[0], (double) ve[1], 0.5*DMPlex_NormD_Internal(dim, v), (double) energy(ctx, c));} 346 } 347 } 348 ierr = DMSwarmRestoreField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);CHKERRQ(ierr); 349 ierr = VecRestoreArrayRead(U, &u);CHKERRQ(ierr); 350 ierr = VecRestoreArray(E, &e);CHKERRQ(ierr); 351 PetscFunctionReturn(0); 352 } 353 354 int main(int argc, char **argv) 355 { 356 TS ts; 357 TSConvergedReason reason; 358 DM dm, sw; 359 Vec u; 360 IS is1, is2; 361 PetscInt *idx1, *idx2; 362 MPI_Comm comm; 363 AppCtx user; 364 const PetscScalar *endVals; 365 PetscReal ftime = .1; 366 PetscInt locSize, dim, d, Np, p, c, steps; 367 PetscErrorCode ierr; 368 369 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 370 comm = PETSC_COMM_WORLD; 371 ierr = ProcessOptions(comm, &user);CHKERRQ(ierr); 372 373 ierr = CreateMesh(comm, &dm, &user);CHKERRQ(ierr); 374 ierr = CreateParticles(dm, &sw, &user);CHKERRQ(ierr); 375 ierr = DMSetApplicationContext(sw, &user);CHKERRQ(ierr); 376 377 ierr = TSCreate(comm, &ts);CHKERRQ(ierr); 378 ierr = TSSetType(ts, TSBASICSYMPLECTIC);CHKERRQ(ierr); 379 ierr = TSSetDM(ts, sw);CHKERRQ(ierr); 380 ierr = TSSetMaxTime(ts, ftime);CHKERRQ(ierr); 381 ierr = TSSetTimeStep(ts, 0.0001);CHKERRQ(ierr); 382 ierr = TSSetMaxSteps(ts, 10);CHKERRQ(ierr); 383 ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr); 384 ierr = TSSetTime(ts, 0.0);CHKERRQ(ierr); 385 ierr = TSSetRHSFunction(ts, NULL, RHSFunctionParticles, &user);CHKERRQ(ierr); 386 387 ierr = DMSwarmCreateGlobalVectorFromField(sw, "kinematics", &u);CHKERRQ(ierr); 388 ierr = DMGetDimension(sw, &dim);CHKERRQ(ierr); 389 ierr = VecGetLocalSize(u, &locSize);CHKERRQ(ierr); 390 Np = locSize/(2*dim); 391 ierr = PetscMalloc1(locSize/2, &idx1);CHKERRQ(ierr); 392 ierr = PetscMalloc1(locSize/2, &idx2);CHKERRQ(ierr); 393 for (p = 0; p < Np; ++p) { 394 for (d = 0; d < dim; ++d) { 395 idx1[p*dim+d] = (p*2+0)*dim + d; 396 idx2[p*dim+d] = (p*2+1)*dim + d; 397 } 398 } 399 ierr = ISCreateGeneral(comm, locSize/2, idx1, PETSC_OWN_POINTER, &is1);CHKERRQ(ierr); 400 ierr = ISCreateGeneral(comm, locSize/2, idx2, PETSC_OWN_POINTER, &is2);CHKERRQ(ierr); 401 ierr = TSRHSSplitSetIS(ts, "position", is1);CHKERRQ(ierr); 402 ierr = TSRHSSplitSetIS(ts, "momentum", is2);CHKERRQ(ierr); 403 ierr = ISDestroy(&is1);CHKERRQ(ierr); 404 ierr = ISDestroy(&is2);CHKERRQ(ierr); 405 ierr = TSRHSSplitSetRHSFunction(ts,"position",NULL,RHSFunction1,&user);CHKERRQ(ierr); 406 ierr = TSRHSSplitSetRHSFunction(ts,"momentum",NULL,RHSFunction2,&user);CHKERRQ(ierr); 407 408 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 409 ierr = TSSetComputeInitialCondition(ts, InitializeSolve);CHKERRQ(ierr); 410 ierr = TSSetComputeExactError(ts, ComputeError);CHKERRQ(ierr); 411 ierr = TSComputeInitialCondition(ts, u);CHKERRQ(ierr); 412 ierr = VecViewFromOptions(u, NULL, "-init_view");CHKERRQ(ierr); 413 ierr = TSSolve(ts, u);CHKERRQ(ierr); 414 ierr = TSGetSolveTime(ts, &ftime);CHKERRQ(ierr); 415 ierr = TSGetConvergedReason(ts, &reason);CHKERRQ(ierr); 416 ierr = TSGetStepNumber(ts, &steps);CHKERRQ(ierr); 417 ierr = PetscPrintf(comm,"%s at time %g after %D steps\n", TSConvergedReasons[reason], (double) ftime, steps);CHKERRQ(ierr); 418 419 ierr = VecGetArrayRead(u, &endVals);CHKERRQ(ierr); 420 for (c = 0; c < 2; ++c) { 421 for (p = 0; p < Np/2; ++p) { 422 const PetscInt n = c*(Np/2) + p; 423 const PetscReal norm = DMPlex_NormD_Internal(dim, &endVals[(n*2 + 1)*dim]); 424 ierr = PetscPrintf(comm, "Particle %D initial Energy: %g Final Energy: %g\n", p, (double) (0.5*(1000./(3*c+1.))), (double) 0.5*PetscSqr(norm));CHKERRQ(ierr); 425 } 426 } 427 ierr = VecRestoreArrayRead(u, &endVals);CHKERRQ(ierr); 428 ierr = DMSwarmDestroyGlobalVectorFromField(sw, "kinematics", &u);CHKERRQ(ierr); 429 ierr = TSDestroy(&ts);CHKERRQ(ierr); 430 ierr = DMDestroy(&sw);CHKERRQ(ierr); 431 ierr = DMDestroy(&dm);CHKERRQ(ierr); 432 ierr = PetscFinalize(); 433 return ierr; 434 } 435 436 437 /*TEST 438 439 build: 440 requires: triangle !single !complex 441 test: 442 suffix: bsi1 443 args: -dim 2 -dm_plex_box_faces 1,1 -dm_view -sw_view -particles_per_circle 5 -ts_basicsymplectic_type 1 -ts_max_time 0.1 -ts_dt 0.001 -ts_monitor_sp_swarm -ts_convergence_estimate -convest_num_refine 2 444 test: 445 suffix: bsi2 446 args: -dim 2 -dm_plex_box_faces 1,1 -dm_view -sw_view -particles_per_circle 5 -ts_basicsymplectic_type 2 -ts_max_time 0.1 -ts_dt 0.001 -ts_monitor_sp_swarm -ts_convergence_estimate -convest_num_refine 2 447 test: 448 suffix: euler 449 args: -dim 2 -dm_plex_box_faces 1,1 -dm_view -sw_view -particles_per_circle 5 -ts_type euler -ts_max_time 0.1 -ts_dt 0.001 -ts_monitor_sp_swarm -ts_convergence_estimate -convest_num_refine 2 450 451 TEST*/ 452