1 static char help[] = "Example of simple hamiltonian system (harmonic oscillator) with particles and a basic symplectic integrator\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 PetscReal omega; /* Oscillation frequency omega */ 14 PetscInt particlesPerCell; /* The number of partices per cell */ 15 PetscReal momentTol; /* Tolerance for checking moment conservation */ 16 PetscBool monitor; /* Flag for using the TS monitor */ 17 PetscBool error; /* Flag for printing the error */ 18 PetscInt ostep; /* print the energy at each ostep time steps */ 19 } AppCtx; 20 21 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 22 { 23 PetscErrorCode ierr; 24 25 PetscFunctionBeginUser; 26 options->dim = 2; 27 options->simplex = PETSC_TRUE; 28 options->monitor = PETSC_FALSE; 29 options->error = PETSC_FALSE; 30 options->particlesPerCell = 1; 31 options->momentTol = 100.0*PETSC_MACHINE_EPSILON; 32 options->omega = 64.0; 33 options->ostep = 100; 34 35 ierr = PetscStrcpy(options->filename, "");CHKERRQ(ierr); 36 37 ierr = PetscOptionsBegin(comm, "", "Harmonic Oscillator 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_cell", "Number of particles per cell", "ex4.c", options->particlesPerCell, &options->particlesPerCell, NULL);CHKERRQ(ierr); 45 ierr = PetscOptionsReal("-omega", "Oscillator frequency", "ex4.c", options->omega, &options->omega, PETSC_NULL);CHKERRQ(ierr); 46 ierr = PetscOptionsEnd();CHKERRQ(ierr); 47 48 PetscFunctionReturn(0); 49 } 50 51 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user) 52 { 53 PetscBool flg; 54 PetscErrorCode ierr; 55 56 PetscFunctionBeginUser; 57 ierr = PetscStrcmp(user->filename, "", &flg);CHKERRQ(ierr); 58 if (flg) { 59 ierr = DMPlexCreateBoxMesh(comm, user->dim, user->simplex, NULL, NULL, NULL, NULL, PETSC_TRUE, dm);CHKERRQ(ierr); 60 } else { 61 ierr = DMPlexCreateFromFile(comm, user->filename, PETSC_TRUE, dm);CHKERRQ(ierr); 62 ierr = DMGetDimension(*dm, &user->dim);CHKERRQ(ierr); 63 } 64 { 65 DM distributedMesh = NULL; 66 67 ierr = DMPlexDistribute(*dm, 0, NULL, &distributedMesh);CHKERRQ(ierr); 68 if (distributedMesh) { 69 ierr = DMDestroy(dm);CHKERRQ(ierr); 70 *dm = distributedMesh; 71 } 72 } 73 ierr = DMLocalizeCoordinates(*dm);CHKERRQ(ierr); /* needed for periodic */ 74 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 75 ierr = PetscObjectSetName((PetscObject) *dm, "Mesh");CHKERRQ(ierr); 76 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 77 PetscFunctionReturn(0); 78 } 79 80 static PetscErrorCode SetInitialCoordinates(DM dmSw) 81 { 82 DM dm; 83 AppCtx *user; 84 PetscRandom rnd; 85 PetscBool simplex; 86 PetscReal *centroid, *coords, *xi0, *v0, *J, *invJ, detJ; 87 PetscInt dim, d, cStart, cEnd, c, Np, p; 88 PetscErrorCode ierr; 89 90 PetscFunctionBeginUser; 91 ierr = PetscRandomCreate(PetscObjectComm((PetscObject) dmSw), &rnd);CHKERRQ(ierr); 92 ierr = PetscRandomSetInterval(rnd, -1.0, 1.0);CHKERRQ(ierr); 93 ierr = PetscRandomSetFromOptions(rnd);CHKERRQ(ierr); 94 95 ierr = DMGetApplicationContext(dmSw, (void **) &user);CHKERRQ(ierr); 96 simplex = user->simplex; 97 Np = user->particlesPerCell; 98 ierr = DMSwarmGetCellDM(dmSw, &dm);CHKERRQ(ierr); 99 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 100 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 101 ierr = PetscMalloc5(dim, ¢roid, dim, &xi0, dim, &v0, dim*dim, &J, dim*dim, &invJ);CHKERRQ(ierr); 102 for (d = 0; d < dim; ++d) xi0[d] = -1.0; 103 ierr = DMSwarmGetField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);CHKERRQ(ierr); 104 for (c = cStart; c < cEnd; ++c) { 105 if (Np == 1) { 106 ierr = DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL);CHKERRQ(ierr); 107 for (d = 0; d < dim; ++d) coords[c*dim+d] = centroid[d]; 108 } else { 109 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); /* affine */ 110 for (p = 0; p < Np; ++p) { 111 const PetscInt n = c*Np + p; 112 PetscReal sum = 0.0, refcoords[3]; 113 114 for (d = 0; d < dim; ++d) { 115 ierr = PetscRandomGetValueReal(rnd, &refcoords[d]);CHKERRQ(ierr); 116 sum += refcoords[d]; 117 } 118 if (simplex && sum > 0.0) for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim)*sum; 119 CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n*dim]); 120 } 121 } 122 } 123 ierr = DMSwarmRestoreField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);CHKERRQ(ierr); 124 ierr = PetscFree5(centroid, xi0, v0, J, invJ);CHKERRQ(ierr); 125 ierr = PetscRandomDestroy(&rnd);CHKERRQ(ierr); 126 PetscFunctionReturn(0); 127 } 128 129 static PetscErrorCode SetInitialConditions(DM dmSw, Vec u) 130 { 131 DM dm; 132 AppCtx *user; 133 PetscReal *coords; 134 PetscScalar *initialConditions; 135 PetscInt dim, cStart, cEnd, c, Np, p; 136 PetscErrorCode ierr; 137 138 PetscFunctionBeginUser; 139 ierr = DMGetApplicationContext(dmSw, (void **) &user);CHKERRQ(ierr); 140 Np = user->particlesPerCell; 141 ierr = DMSwarmGetCellDM(dmSw, &dm);CHKERRQ(ierr); 142 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 143 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 144 ierr = DMSwarmGetField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);CHKERRQ(ierr); 145 ierr = VecGetArray(u, &initialConditions);CHKERRQ(ierr); 146 for (c = cStart; c < cEnd; ++c) { 147 for (p = 0; p < Np; ++p) { 148 const PetscInt n = c*Np + p; 149 150 initialConditions[n*2+0] = DMPlex_NormD_Internal(dim, &coords[n*dim]); 151 initialConditions[n*2+1] = 0.0; 152 } 153 } 154 ierr = VecRestoreArray(u, &initialConditions);CHKERRQ(ierr); 155 ierr = DMSwarmRestoreField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);CHKERRQ(ierr); 156 PetscFunctionReturn(0); 157 } 158 159 static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user) 160 { 161 PetscInt *cellid; 162 PetscInt dim, cStart, cEnd, c, Np = user->particlesPerCell, p; 163 PetscErrorCode ierr; 164 165 PetscFunctionBeginUser; 166 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 167 ierr = DMCreate(PetscObjectComm((PetscObject) dm), sw);CHKERRQ(ierr); 168 ierr = DMSetType(*sw, DMSWARM);CHKERRQ(ierr); 169 ierr = DMSetDimension(*sw, dim);CHKERRQ(ierr); 170 171 ierr = DMSwarmSetType(*sw, DMSWARM_PIC);CHKERRQ(ierr); 172 ierr = DMSwarmSetCellDM(*sw, dm);CHKERRQ(ierr); 173 ierr = DMSwarmRegisterPetscDatatypeField(*sw, "kinematics", 2, PETSC_REAL);CHKERRQ(ierr); 174 ierr = DMSwarmFinalizeFieldRegister(*sw);CHKERRQ(ierr); 175 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 176 ierr = DMSwarmSetLocalSizes(*sw, (cEnd - cStart) * Np, 0);CHKERRQ(ierr); 177 ierr = DMSetFromOptions(*sw);CHKERRQ(ierr); 178 ierr = DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);CHKERRQ(ierr); 179 for (c = cStart; c < cEnd; ++c) { 180 for (p = 0; p < Np; ++p) { 181 const PetscInt n = c*Np + p; 182 183 cellid[n] = c; 184 } 185 } 186 ierr = DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);CHKERRQ(ierr); 187 ierr = PetscObjectSetName((PetscObject) *sw, "Particles");CHKERRQ(ierr); 188 ierr = DMViewFromOptions(*sw, NULL, "-sw_view");CHKERRQ(ierr); 189 PetscFunctionReturn(0); 190 } 191 192 /* Create particle RHS Functions */ 193 static PetscErrorCode RHSFunction1(TS ts, PetscReal t, Vec V, Vec Xres, void *ctx) 194 { 195 const PetscScalar *v; 196 PetscScalar *xres; 197 PetscInt Np, p; 198 PetscErrorCode ierr; 199 200 PetscFunctionBeginUser; 201 ierr = VecGetArray(Xres, &xres);CHKERRQ(ierr); 202 ierr = VecGetArrayRead(V, &v);CHKERRQ(ierr); 203 ierr = VecGetLocalSize(Xres, &Np);CHKERRQ(ierr); 204 for (p = 0; p < Np; ++p) { 205 xres[p] = v[p]; 206 } 207 ierr = VecRestoreArrayRead(V, &v);CHKERRQ(ierr); 208 ierr = VecRestoreArray(Xres, &xres);CHKERRQ(ierr); 209 PetscFunctionReturn(0); 210 } 211 212 static PetscErrorCode RHSFunction2(TS ts, PetscReal t, Vec X, Vec Vres, void *ctx) 213 { 214 AppCtx *user = (AppCtx *)ctx; 215 const PetscScalar *x; 216 PetscScalar *vres; 217 PetscInt Np, p; 218 PetscErrorCode ierr; 219 220 PetscFunctionBeginUser; 221 ierr = VecGetArray(Vres, &vres);CHKERRQ(ierr); 222 ierr = VecGetArrayRead(X, &x);CHKERRQ(ierr); 223 ierr = VecGetLocalSize(Vres, &Np);CHKERRQ(ierr); 224 for (p = 0; p < Np; ++p) { 225 vres[p] = -PetscSqr(user->omega)*x[p]; 226 } 227 ierr = VecRestoreArrayRead(X, &x);CHKERRQ(ierr); 228 ierr = VecRestoreArray(Vres, &vres);CHKERRQ(ierr); 229 PetscFunctionReturn(0); 230 } 231 232 static PetscErrorCode RHSFunctionParticles(TS ts, PetscReal t, Vec U, Vec R, void *ctx) 233 { 234 AppCtx *user = (AppCtx *) ctx; 235 DM dm; 236 const PetscScalar *u; 237 PetscScalar *r; 238 PetscInt Np, p; 239 PetscErrorCode ierr; 240 241 PetscFunctionBeginUser; 242 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 243 ierr = VecGetArray(R, &r);CHKERRQ(ierr); 244 ierr = VecGetArrayRead(U, &u);CHKERRQ(ierr); 245 ierr = VecGetLocalSize(U, &Np);CHKERRQ(ierr); 246 Np /= 2; 247 for (p = 0; p < Np; ++p) { 248 r[p*2+0] = u[p*2+1]; 249 r[p*2+1] = -PetscSqr(user->omega)*u[p*2+0]; 250 } 251 ierr = VecRestoreArrayRead(U, &u);CHKERRQ(ierr); 252 ierr = VecRestoreArray(R, &r);CHKERRQ(ierr); 253 PetscFunctionReturn(0); 254 } 255 256 static PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx) 257 { 258 AppCtx *user = (AppCtx *) ctx; 259 const PetscReal omega = user->omega; 260 const PetscScalar *u; 261 MPI_Comm comm; 262 PetscReal dt; 263 PetscInt Np, p; 264 PetscErrorCode ierr; 265 266 PetscFunctionBeginUser; 267 if (step%user->ostep == 0) { 268 ierr = PetscObjectGetComm((PetscObject) ts, &comm);CHKERRQ(ierr); 269 if (!step) {ierr = PetscPrintf(comm, "Time Step Part Energy Mod Energy\n");CHKERRQ(ierr);} 270 ierr = TSGetTimeStep(ts, &dt);CHKERRQ(ierr); 271 ierr = VecGetArrayRead(U, &u);CHKERRQ(ierr); 272 ierr = VecGetLocalSize(U, &Np);CHKERRQ(ierr); 273 Np /= 2; 274 for (p = 0; p < Np; ++p) { 275 const PetscReal x = PetscRealPart(u[p*2+0]); 276 const PetscReal v = PetscRealPart(u[p*2+1]); 277 const PetscReal E = 0.5*(v*v + PetscSqr(omega)*x*x); 278 const PetscReal mE = 0.5*(v*v + PetscSqr(omega)*x*x - PetscSqr(omega)*dt*x*v); 279 280 ierr = PetscPrintf(comm, "%.6lf %4D %4D %10.4lf %10.4lf\n", t, step, p, (double) E, (double) mE);CHKERRQ(ierr); 281 } 282 ierr = VecRestoreArrayRead(U, &u);CHKERRQ(ierr); 283 } 284 PetscFunctionReturn(0); 285 } 286 287 static PetscErrorCode InitializeSolve(TS ts, Vec u) 288 { 289 DM dm; 290 AppCtx *user; 291 PetscErrorCode ierr; 292 293 PetscFunctionBeginUser; 294 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 295 ierr = DMGetApplicationContext(dm, (void **) &user);CHKERRQ(ierr); 296 ierr = SetInitialCoordinates(dm);CHKERRQ(ierr); 297 ierr = SetInitialConditions(dm, u);CHKERRQ(ierr); 298 PetscFunctionReturn(0); 299 } 300 301 static PetscErrorCode ComputeError(TS ts, Vec U, Vec E) 302 { 303 MPI_Comm comm; 304 DM sdm; 305 AppCtx *user; 306 const PetscScalar *u, *coords; 307 PetscScalar *e; 308 PetscReal t, omega; 309 PetscInt dim, Np, p; 310 PetscErrorCode ierr; 311 312 PetscFunctionBeginUser; 313 ierr = PetscObjectGetComm((PetscObject) ts, &comm);CHKERRQ(ierr); 314 ierr = TSGetDM(ts, &sdm);CHKERRQ(ierr); 315 ierr = DMGetApplicationContext(sdm, (void **) &user);CHKERRQ(ierr); 316 omega = user->omega; 317 ierr = DMGetDimension(sdm, &dim);CHKERRQ(ierr); 318 ierr = TSGetSolveTime(ts, &t);CHKERRQ(ierr); 319 ierr = VecGetArray(E, &e);CHKERRQ(ierr); 320 ierr = VecGetArrayRead(U, &u);CHKERRQ(ierr); 321 ierr = VecGetLocalSize(U, &Np);CHKERRQ(ierr); 322 Np /= 2; 323 ierr = DMSwarmGetField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);CHKERRQ(ierr); 324 for (p = 0; p < Np; ++p) { 325 const PetscReal x = PetscRealPart(u[p*2+0]); 326 const PetscReal v = PetscRealPart(u[p*2+1]); 327 const PetscReal x0 = DMPlex_NormD_Internal(dim, &coords[p*dim]); 328 const PetscReal ex = x0*PetscCosReal(omega*t); 329 const PetscReal ev = -x0*omega*PetscSinReal(omega*t); 330 331 if (user->error) {ierr = PetscPrintf(comm, "p%D error [%.2g %.2g] sol [%.6lf %.6lf] exact [%.6lf %.6lf] energy/exact energy %g / %g\n", p, (double) PetscAbsReal(x-ex), (double) PetscAbsReal(v-ev), (double) x, (double) v, (double) ex, (double) ev, 0.5*(v*v + PetscSqr(omega)*x*x), (double) 0.5*PetscSqr(omega*x0));} 332 e[p*2+0] = x - ex; 333 e[p*2+1] = v - ev; 334 } 335 ierr = DMSwarmRestoreField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);CHKERRQ(ierr); 336 ierr = VecRestoreArrayRead(U, &u);CHKERRQ(ierr); 337 ierr = VecRestoreArray(E, &e);CHKERRQ(ierr); 338 PetscFunctionReturn(0); 339 } 340 341 int main(int argc,char **argv) 342 { 343 TS ts; /* nonlinear solver */ 344 DM dm, sw; /* Mesh and particle managers */ 345 Vec u; /* swarm vector */ 346 IS is1, is2; 347 PetscInt n; 348 MPI_Comm comm; 349 AppCtx user; 350 PetscErrorCode ierr; 351 352 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 353 comm = PETSC_COMM_WORLD; 354 ierr = ProcessOptions(comm, &user);CHKERRQ(ierr); 355 356 ierr = CreateMesh(comm, &dm, &user);CHKERRQ(ierr); 357 ierr = CreateParticles(dm, &sw, &user);CHKERRQ(ierr); 358 ierr = DMSetApplicationContext(sw, &user);CHKERRQ(ierr); 359 360 ierr = TSCreate(comm, &ts);CHKERRQ(ierr); 361 ierr = TSSetType(ts, TSBASICSYMPLECTIC);CHKERRQ(ierr); 362 ierr = TSSetDM(ts, sw);CHKERRQ(ierr); 363 ierr = TSSetMaxTime(ts, 0.1);CHKERRQ(ierr); 364 ierr = TSSetTimeStep(ts, 0.00001);CHKERRQ(ierr); 365 ierr = TSSetMaxSteps(ts, 100);CHKERRQ(ierr); 366 ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr); 367 if (user.monitor) {ierr = TSMonitorSet(ts, Monitor, &user, NULL);CHKERRQ(ierr);} 368 ierr = TSSetRHSFunction(ts, NULL, RHSFunctionParticles, &user);CHKERRQ(ierr); 369 370 ierr = DMSwarmCreateGlobalVectorFromField(sw, "kinematics", &u);CHKERRQ(ierr); 371 ierr = VecGetLocalSize(u, &n);CHKERRQ(ierr); 372 ierr = ISCreateStride(comm, n/2, 0, 2, &is1);CHKERRQ(ierr); 373 ierr = ISCreateStride(comm, n/2, 1, 2, &is2);CHKERRQ(ierr); 374 ierr = TSRHSSplitSetIS(ts, "position", is1);CHKERRQ(ierr); 375 ierr = TSRHSSplitSetIS(ts, "momentum", is2);CHKERRQ(ierr); 376 ierr = ISDestroy(&is1);CHKERRQ(ierr); 377 ierr = ISDestroy(&is2);CHKERRQ(ierr); 378 ierr = TSRHSSplitSetRHSFunction(ts, "position", NULL, RHSFunction1, &user);CHKERRQ(ierr); 379 ierr = TSRHSSplitSetRHSFunction(ts, "momentum", NULL, RHSFunction2, &user);CHKERRQ(ierr); 380 381 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 382 ierr = TSSetComputeInitialCondition(ts, InitializeSolve);CHKERRQ(ierr); 383 ierr = TSSetComputeExactError(ts, ComputeError);CHKERRQ(ierr); 384 ierr = TSComputeInitialCondition(ts, u);CHKERRQ(ierr); 385 ierr = TSSolve(ts, u);CHKERRQ(ierr); 386 ierr = DMSwarmDestroyGlobalVectorFromField(sw, "kinematics", &u);CHKERRQ(ierr); 387 388 ierr = TSDestroy(&ts);CHKERRQ(ierr); 389 ierr = DMDestroy(&sw);CHKERRQ(ierr); 390 ierr = DMDestroy(&dm);CHKERRQ(ierr); 391 ierr = PetscFinalize(); 392 return ierr; 393 } 394 395 /*TEST 396 397 build: 398 requires: triangle !single !complex 399 test: 400 suffix: 1 401 args: -dm_plex_box_faces 1,1 -ts_basicsymplectic_type 1 -ts_convergence_estimate -convest_num_refine 2 -dm_view -sw_view -monitor -output_step 50 -error 402 test: 403 suffix: 2 404 args: -dm_plex_box_faces 1,1 -ts_basicsymplectic_type 2 -ts_convergence_estimate -convest_num_refine 2 -dm_view -sw_view -monitor -output_step 50 -error 405 test: 406 suffix: 3 407 args: -dm_plex_box_faces 1,1 -ts_basicsymplectic_type 3 -ts_convergence_estimate -convest_num_refine 2 -ts_dt 0.0001 -dm_view -sw_view -monitor -output_step 50 -error 408 test: 409 suffix: 4 410 args: -dm_plex_box_faces 1,1 -ts_basicsymplectic_type 4 -ts_convergence_estimate -convest_num_refine 2 -ts_dt 0.0001 -dm_view -sw_view -monitor -output_step 50 -error 411 412 TEST*/ 413