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