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