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