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");PetscCall(ierr); 32 PetscCall(PetscOptionsInt("-output_step", "Number of time steps between output", "ex4.c", options->ostep, &options->ostep, PETSC_NULL)); 33 PetscCall(PetscOptionsBool("-monitor", "Flag to use the TS monitor", "ex4.c", options->monitor, &options->monitor, NULL)); 34 PetscCall(PetscOptionsBool("-error", "Flag to print the error", "ex4.c", options->error, &options->error, NULL)); 35 PetscCall(PetscOptionsInt("-particles_per_circle", "Number of particles per circle", "ex4.c", options->particlesPerCircle, &options->particlesPerCircle, NULL)); 36 ierr = PetscOptionsEnd();PetscCall(ierr); 37 38 PetscFunctionReturn(0); 39 } 40 41 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user) 42 { 43 PetscFunctionBeginUser; 44 PetscCall(DMCreate(comm, dm)); 45 PetscCall(DMSetType(*dm, DMPLEX)); 46 PetscCall(DMSetFromOptions(*dm)); 47 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 48 PetscCall(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 PetscCall(DMGetApplicationContext(dmSw, &ctx)); 107 Np = ctx->particlesPerCircle; 108 PetscCall(DMSwarmGetCellDM(dmSw, &dm)); 109 PetscCall(DMGetDimension(dm, &dim)); 110 PetscCall(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 PetscCall(orbit(ctx, c, p, 0.0, &coords[n*dim], NULL)); 118 } 119 } 120 PetscCall(DMSwarmRestoreField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords)); 121 122 PetscCall(DMSwarmCreateGlobalVectorFromField(dmSw, DMSwarmPICField_coor, &coordinates)); 123 PetscCall(DMLocatePoints(dm, coordinates, DM_POINTLOCATION_NONE, &cellSF)); 124 PetscCall(DMSwarmDestroyGlobalVectorFromField(dmSw, DMSwarmPICField_coor, &coordinates)); 125 126 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject) dmSw), &rank)); 127 PetscCall(DMSwarmGetField(dmSw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid)); 128 PetscCall(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 PetscCall(DMSwarmRestoreField(dmSw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid)); 136 PetscCall(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 PetscCall(DMGetApplicationContext(dmSw, &ctx)); 149 Np = ctx->particlesPerCircle; 150 PetscCall(DMSwarmGetCellDM(dmSw, &dm)); 151 PetscCall(DMGetDimension(dm, &dim)); 152 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 153 PetscCall(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 PetscCall(orbit(ctx, c, p, 0.0, &initialConditions[(n*2 + 0)*dim], &initialConditions[(n*2 + 1)*dim])); 159 } 160 } 161 PetscCall(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 PetscCall(DMGetDimension(dm, &dim)); 171 PetscCall(DMCreate(PetscObjectComm((PetscObject) dm), sw)); 172 PetscCall(DMSetType(*sw, DMSWARM)); 173 PetscCall(DMSetDimension(*sw, dim)); 174 175 PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC)); 176 PetscCall(DMSwarmSetCellDM(*sw, dm)); 177 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "kinematics", 2*dim, PETSC_REAL)); 178 PetscCall(DMSwarmFinalizeFieldRegister(*sw)); 179 PetscCall(DMSwarmSetLocalSizes(*sw, 2 * Np, 0)); 180 PetscCall(DMSetFromOptions(*sw)); 181 PetscCall(PetscObjectSetName((PetscObject) *sw, "Particles")); 182 PetscCall(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 PetscCall(VecGetLocalSize(Xres, &Np)); 197 Np /= dim; 198 PetscCall(VecGetArray(Xres, &xres)); 199 PetscCall(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 PetscCall(VecRestoreArrayRead(V,& v)); 206 PetscCall(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 PetscCall(VecGetLocalSize(Vres, &Np)); 221 Np /= dim; 222 PetscCall(VecGetArray(Vres, &vres)); 223 PetscCall(VecGetArrayRead(X, &x)); 224 for (p = 0; p < Np; ++p) { 225 const PetscInt c = p / ctx->particlesPerCircle; 226 227 PetscCall(force(ctx, c, &x[p*dim], &vres[p*dim])); 228 } 229 PetscCall(VecRestoreArray(Vres, &vres)); 230 PetscCall(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 PetscCall(TSGetDM(ts, &dm)); 244 PetscCall(DMGetDimension(dm, &dim)); 245 PetscCall(VecGetLocalSize(U, &Np)); 246 Np /= 2*dim; 247 PetscCall(VecGetArrayRead(U, &u)); 248 PetscCall(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 PetscCall(force(ctx, c, &u[(p*2 + 0)*dim], &r[(p*2 + 1)*dim])); 254 } 255 PetscCall(VecRestoreArrayRead(U, &u)); 256 PetscCall(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 PetscCall(TSGetDM(ts, &dm)); 267 PetscCall(DMGetApplicationContext(dm, &user)); 268 PetscCall(SetInitialCoordinates(dm)); 269 PetscCall(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 PetscCall(PetscObjectGetComm((PetscObject) ts, &comm)); 285 PetscCall(TSGetDM(ts, &sdm)); 286 PetscCall(DMGetApplicationContext(sdm, &ctx)); 287 PetscCall(DMGetDimension(sdm, &dim)); 288 PetscCall(TSGetSolveTime(ts, &t)); 289 PetscCall(VecGetArray(E, &e)); 290 PetscCall(VecGetArrayRead(U, &u)); 291 PetscCall(VecGetLocalSize(U, &Np)); 292 Np /= 2*dim*2; 293 PetscCall(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 PetscCall(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) PetscCall(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 PetscCall(DMSwarmRestoreField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords)); 311 PetscCall(VecRestoreArrayRead(U, &u)); 312 PetscCall(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 330 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 331 comm = PETSC_COMM_WORLD; 332 PetscCall(ProcessOptions(comm, &user)); 333 334 PetscCall(CreateMesh(comm, &dm, &user)); 335 PetscCall(CreateParticles(dm, &sw, &user)); 336 PetscCall(DMSetApplicationContext(sw, &user)); 337 338 PetscCall(TSCreate(comm, &ts)); 339 PetscCall(TSSetType(ts, TSBASICSYMPLECTIC)); 340 PetscCall(TSSetDM(ts, sw)); 341 PetscCall(TSSetMaxTime(ts, ftime)); 342 PetscCall(TSSetTimeStep(ts, 0.0001)); 343 PetscCall(TSSetMaxSteps(ts, 10)); 344 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 345 PetscCall(TSSetTime(ts, 0.0)); 346 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunctionParticles, &user)); 347 348 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "kinematics", &u)); 349 PetscCall(DMGetDimension(sw, &dim)); 350 PetscCall(VecGetLocalSize(u, &locSize)); 351 Np = locSize/(2*dim); 352 PetscCall(PetscMalloc1(locSize/2, &idx1)); 353 PetscCall(PetscMalloc1(locSize/2, &idx2)); 354 for (p = 0; p < Np; ++p) { 355 for (d = 0; d < dim; ++d) { 356 idx1[p*dim+d] = (p*2+0)*dim + d; 357 idx2[p*dim+d] = (p*2+1)*dim + d; 358 } 359 } 360 PetscCall(ISCreateGeneral(comm, locSize/2, idx1, PETSC_OWN_POINTER, &is1)); 361 PetscCall(ISCreateGeneral(comm, locSize/2, idx2, PETSC_OWN_POINTER, &is2)); 362 PetscCall(TSRHSSplitSetIS(ts, "position", is1)); 363 PetscCall(TSRHSSplitSetIS(ts, "momentum", is2)); 364 PetscCall(ISDestroy(&is1)); 365 PetscCall(ISDestroy(&is2)); 366 PetscCall(TSRHSSplitSetRHSFunction(ts,"position",NULL,RHSFunction1,&user)); 367 PetscCall(TSRHSSplitSetRHSFunction(ts,"momentum",NULL,RHSFunction2,&user)); 368 369 PetscCall(TSSetFromOptions(ts)); 370 PetscCall(TSSetComputeInitialCondition(ts, InitializeSolve)); 371 PetscCall(TSSetComputeExactError(ts, ComputeError)); 372 PetscCall(TSComputeInitialCondition(ts, u)); 373 PetscCall(VecViewFromOptions(u, NULL, "-init_view")); 374 PetscCall(TSSolve(ts, u)); 375 PetscCall(TSGetSolveTime(ts, &ftime)); 376 PetscCall(TSGetConvergedReason(ts, &reason)); 377 PetscCall(TSGetStepNumber(ts, &steps)); 378 PetscCall(PetscPrintf(comm,"%s at time %g after %D steps\n", TSConvergedReasons[reason], (double) ftime, steps)); 379 380 PetscCall(VecGetArrayRead(u, &endVals)); 381 for (c = 0; c < 2; ++c) { 382 for (p = 0; p < Np/2; ++p) { 383 const PetscInt n = c*(Np/2) + p; 384 const PetscReal norm = DMPlex_NormD_Internal(dim, &endVals[(n*2 + 1)*dim]); 385 PetscCall(PetscPrintf(comm, "Particle %D initial Energy: %g Final Energy: %g\n", p, (double) (0.5*(1000./(3*c+1.))), (double) 0.5*PetscSqr(norm))); 386 } 387 } 388 PetscCall(VecRestoreArrayRead(u, &endVals)); 389 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "kinematics", &u)); 390 PetscCall(TSDestroy(&ts)); 391 PetscCall(DMDestroy(&sw)); 392 PetscCall(DMDestroy(&dm)); 393 PetscCall(PetscFinalize()); 394 return 0; 395 } 396 397 /*TEST 398 399 build: 400 requires: triangle !single !complex 401 test: 402 suffix: bsi1 403 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 404 test: 405 suffix: bsi2 406 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 407 test: 408 suffix: euler 409 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 410 411 TEST*/ 412