Actual source code: ex4.c

petsc-3.13.4 2020-08-01
Report Typos and Errors
  1: static char help[] = "Example of simple hamiltonian system (harmonic oscillator) with particles and a basic symplectic integrator\n";

  3:  #include <petscdmplex.h>
  4:  #include <petsc/private/dmpleximpl.h>
  5:  #include <petsc/private/petscfeimpl.h>
  6:  #include <petscdmswarm.h>
  7:  #include <petscts.h>

  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;

 21: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 22: {

 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;

 35:   PetscStrcpy(options->filename, "");

 37:   PetscOptionsBegin(comm, "", "Harmonic Oscillator Options", "DMPLEX");
 38:   PetscOptionsInt("-output_step", "Number of time steps between output", "ex4.c", options->ostep, &options->ostep, PETSC_NULL);
 39:   PetscOptionsInt("-dim", "The topological mesh dimension", "ex4.c", options->dim, &options->dim, NULL);
 40:   PetscOptionsBool("-monitor", "Flag to use the TS monitor", "ex4.c", options->monitor, &options->monitor, NULL);
 41:   PetscOptionsBool("-error", "Flag to print the error", "ex4.c", options->error, &options->error, NULL);
 42:   PetscOptionsBool("-simplex", "The flag for simplices or tensor cells", "ex4.c", options->simplex, &options->simplex, NULL);
 43:   PetscOptionsString("-mesh", "Name of the mesh filename if any", "ex4.c", options->filename, options->filename, PETSC_MAX_PATH_LEN, NULL);
 44:   PetscOptionsInt("-particles_per_cell", "Number of particles per cell", "ex4.c", options->particlesPerCell, &options->particlesPerCell, NULL);
 45:   PetscOptionsReal("-omega", "Oscillator frequency", "ex4.c", options->omega, &options->omega, PETSC_NULL);
 46:   PetscOptionsEnd();

 48:   return(0);
 49: }

 51: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
 52: {
 53:   PetscBool      flg;

 57:   PetscStrcmp(user->filename, "", &flg);
 58:   if (flg) {
 59:     DMPlexCreateBoxMesh(comm, user->dim, user->simplex, NULL, NULL, NULL, NULL, PETSC_TRUE, dm);
 60:   } else {
 61:     DMPlexCreateFromFile(comm, user->filename, PETSC_TRUE, dm);
 62:     DMGetDimension(*dm, &user->dim);
 63:   }
 64:   {
 65:     DM distributedMesh = NULL;

 67:     DMPlexDistribute(*dm, 0, NULL, &distributedMesh);
 68:     if (distributedMesh) {
 69:       DMDestroy(dm);
 70:       *dm  = distributedMesh;
 71:     }
 72:   }
 73:   DMLocalizeCoordinates(*dm); /* needed for periodic */
 74:   DMSetFromOptions(*dm);
 75:   PetscObjectSetName((PetscObject) *dm, "Mesh");
 76:   DMViewFromOptions(*dm, NULL, "-dm_view");
 77:   return(0);
 78: }

 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;

 91:   PetscRandomCreate(PetscObjectComm((PetscObject) dmSw), &rnd);
 92:   PetscRandomSetInterval(rnd, -1.0, 1.0);
 93:   PetscRandomSetFromOptions(rnd);

 95:   DMGetApplicationContext(dmSw, (void **) &user);
 96:   simplex = user->simplex;
 97:   Np   = user->particlesPerCell;
 98:   DMSwarmGetCellDM(dmSw, &dm);
 99:   DMGetDimension(dm, &dim);
100:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
101:   PetscMalloc5(dim, &centroid, dim, &xi0, dim, &v0, dim*dim, &J, dim*dim, &invJ);
102:   for (d = 0; d < dim; ++d) xi0[d] = -1.0;
103:   DMSwarmGetField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
104:   for (c = cStart; c < cEnd; ++c) {
105:     if (Np == 1) {
106:       DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL);
107:       for (d = 0; d < dim; ++d) coords[c*dim+d] = centroid[d];
108:     } else {
109:       DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ); /* affine */
110:       for (p = 0; p < Np; ++p) {
111:         const PetscInt n   = c*Np + p;
112:         PetscReal      sum = 0.0, refcoords[3];

114:         for (d = 0; d < dim; ++d) {
115:           PetscRandomGetValueReal(rnd, &refcoords[d]);
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:   DMSwarmRestoreField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
124:   PetscFree5(centroid, xi0, v0, J, invJ);
125:   PetscRandomDestroy(&rnd);
126:   return(0);
127: }

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;

139:   DMGetApplicationContext(dmSw, (void **) &user);
140:   Np   = user->particlesPerCell;
141:   DMSwarmGetCellDM(dmSw, &dm);
142:   DMGetDimension(dm, &dim);
143:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
144:   DMSwarmGetField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
145:   VecGetArray(u, &initialConditions);
146:   for (c = cStart; c < cEnd; ++c) {
147:     for (p = 0; p < Np; ++p) {
148:       const PetscInt n = c*Np + p;

150:       initialConditions[n*2+0] = DMPlex_NormD_Internal(dim, &coords[n*dim]);
151:       initialConditions[n*2+1] = 0.0;
152:     }
153:   }
154:   VecRestoreArray(u, &initialConditions);
155:   DMSwarmRestoreField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
156:   return(0);
157: }

159: static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
160: {
161:   PetscInt      *cellid;
162:   PetscInt       dim, cStart, cEnd, c, Np = user->particlesPerCell, p;

166:   DMGetDimension(dm, &dim);
167:   DMCreate(PetscObjectComm((PetscObject) dm), sw);
168:   DMSetType(*sw, DMSWARM);
169:   DMSetDimension(*sw, dim);

171:   DMSwarmSetType(*sw, DMSWARM_PIC);
172:   DMSwarmSetCellDM(*sw, dm);
173:   DMSwarmRegisterPetscDatatypeField(*sw, "kinematics", 2, PETSC_REAL);
174:   DMSwarmFinalizeFieldRegister(*sw);
175:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
176:   DMSwarmSetLocalSizes(*sw, (cEnd - cStart) * Np, 0);
177:   DMSetFromOptions(*sw);
178:   DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
179:   for (c = cStart; c < cEnd; ++c) {
180:     for (p = 0; p < Np; ++p) {
181:       const PetscInt n = c*Np + p;

183:       cellid[n] = c;
184:     }
185:   }
186:   DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
187:   PetscObjectSetName((PetscObject) *sw, "Particles");
188:   DMViewFromOptions(*sw, NULL, "-sw_view");
189:   return(0);
190: }

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;

201:   VecGetArray(Xres, &xres);
202:   VecGetArrayRead(V, &v);
203:   VecGetLocalSize(Xres, &Np);
204:   for (p = 0; p < Np; ++p) {
205:      xres[p] = v[p];
206:   }
207:   VecRestoreArrayRead(V, &v);
208:   VecRestoreArray(Xres, &xres);
209:   return(0);
210: }

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;

221:   VecGetArray(Vres, &vres);
222:   VecGetArrayRead(X, &x);
223:   VecGetLocalSize(Vres, &Np);
224:   for (p = 0; p < Np; ++p) {
225:     vres[p] = -PetscSqr(user->omega)*x[p];
226:   }
227:   VecRestoreArrayRead(X, &x);
228:   VecRestoreArray(Vres, &vres);
229:   return(0);
230: }

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;

242:   TSGetDM(ts, &dm);
243:   VecGetArray(R, &r);
244:   VecGetArrayRead(U, &u);
245:   VecGetLocalSize(U, &Np);
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:   VecRestoreArrayRead(U, &u);
252:   VecRestoreArray(R, &r);
253:   return(0);
254: }

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;

267:   if (step%user->ostep == 0) {
268:     PetscObjectGetComm((PetscObject) ts, &comm);
269:     if (!step) {PetscPrintf(comm, "Time     Step Part     Energy Mod Energy\n");}
270:     TSGetTimeStep(ts, &dt);
271:     VecGetArrayRead(U, &u);
272:     VecGetLocalSize(U, &Np);
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);

280:       PetscPrintf(comm, "%.6lf %4D %4D %10.4lf %10.4lf\n", t, step, p, (double) E, (double) mE);
281:     }
282:     VecRestoreArrayRead(U, &u);
283:   }
284:   return(0);
285: }

287: static PetscErrorCode InitializeSolve(TS ts, Vec u)
288: {
289:   DM             dm;
290:   AppCtx        *user;

294:   TSGetDM(ts, &dm);
295:   DMGetApplicationContext(dm, (void **) &user);
296:   SetInitialCoordinates(dm);
297:   SetInitialConditions(dm, u);
298:   return(0);
299: }

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;

313:   PetscObjectGetComm((PetscObject) ts, &comm);
314:   TSGetDM(ts, &sdm);
315:   DMGetApplicationContext(sdm, (void **) &user);
316:   omega = user->omega;
317:   DMGetDimension(sdm, &dim);
318:   TSGetSolveTime(ts, &t);
319:   VecGetArray(E, &e);
320:   VecGetArrayRead(U, &u);
321:   VecGetLocalSize(U, &Np);
322:   Np  /= 2;
323:   DMSwarmGetField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
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);

331:     if (user->error) {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:   DMSwarmRestoreField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
336:   VecRestoreArrayRead(U, &u);
337:   VecRestoreArray(E, &e);
338:   return(0);
339: }

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;

352:   PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
353:   comm = PETSC_COMM_WORLD;
354:   ProcessOptions(comm, &user);

356:   CreateMesh(comm, &dm, &user);
357:   CreateParticles(dm, &sw, &user);
358:   DMSetApplicationContext(sw, &user);

360:   TSCreate(comm, &ts);
361:   TSSetType(ts, TSBASICSYMPLECTIC);
362:   TSSetDM(ts, sw);
363:   TSSetMaxTime(ts, 0.1);
364:   TSSetTimeStep(ts, 0.00001);
365:   TSSetMaxSteps(ts, 100);
366:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
367:   if (user.monitor) {TSMonitorSet(ts, Monitor, &user, NULL);}
368:   TSSetRHSFunction(ts, NULL, RHSFunctionParticles, &user);

370:   DMSwarmCreateGlobalVectorFromField(sw, "kinematics", &u);
371:   VecGetLocalSize(u, &n);
372:   ISCreateStride(comm, n/2, 0, 2, &is1);
373:   ISCreateStride(comm, n/2, 1, 2, &is2);
374:   TSRHSSplitSetIS(ts, "position", is1);
375:   TSRHSSplitSetIS(ts, "momentum", is2);
376:   ISDestroy(&is1);
377:   ISDestroy(&is2);
378:   TSRHSSplitSetRHSFunction(ts, "position", NULL, RHSFunction1, &user);
379:   TSRHSSplitSetRHSFunction(ts, "momentum", NULL, RHSFunction2, &user);

381:   TSSetFromOptions(ts);
382:   TSSetComputeInitialCondition(ts, InitializeSolve);
383:   TSSetComputeExactError(ts, ComputeError);
384:   TSComputeInitialCondition(ts, u);
385:   TSSolve(ts, u);
386:   DMSwarmDestroyGlobalVectorFromField(sw, "kinematics", &u);

388:   TSDestroy(&ts);
389:   DMDestroy(&sw);
390:   DMDestroy(&dm);
391:   PetscFinalize();
392:   return ierr;
393: }

395: /*TEST

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

412: TEST*/