1 static char help[] = "Example program demonstrating projection between particle and finite element spaces\n\n";
2
3 #include <petscdmplex.h>
4 #include <petscds.h>
5 #include <petscdmswarm.h>
6 #include <petscksp.h>
7
main(int argc,char ** argv)8 int main(int argc, char **argv)
9 {
10 DM dm, sw;
11 PetscFE fe;
12 Vec u_f;
13 DMPolytopeType ct;
14 PetscInt dim, Nc = 1, faces[3];
15 PetscInt Np = 10, field = 0, zero = 0, bs, cStart;
16 PetscReal energy_0 = 0, energy_1 = 0;
17 PetscReal lo[3], hi[3], h[3];
18 PetscBool removePoints = PETSC_TRUE;
19 PetscReal *wq, *coords;
20 PetscDataType dtype;
21
22 PetscFunctionBeginUser;
23 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
24 /* Create a mesh */
25 PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
26 PetscCall(DMSetType(dm, DMPLEX));
27 PetscCall(DMSetFromOptions(dm));
28 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
29
30 PetscCall(DMGetDimension(dm, &dim));
31 bs = dim;
32 PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-dm_plex_box_faces", faces, &bs, NULL));
33 PetscCall(PetscOptionsGetInt(NULL, NULL, "-np", &Np, NULL));
34 PetscCall(DMGetBoundingBox(dm, lo, hi));
35 for (PetscInt i = 0; i < dim; ++i) {
36 h[i] = (hi[i] - lo[i]) / faces[i];
37 PetscCall(PetscPrintf(PETSC_COMM_SELF, " lo = %g hi = %g n = %" PetscInt_FMT " h = %g\n", (double)lo[i], (double)hi[i], faces[i], (double)h[i]));
38 }
39 // Create FE space
40 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
41 PetscCall(DMPlexGetCellType(dm, cStart, &ct));
42 PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, Nc, ct, NULL, PETSC_DECIDE, &fe));
43 PetscCall(PetscFESetFromOptions(fe));
44 PetscCall(PetscObjectSetName((PetscObject)fe, "fe"));
45 PetscCall(DMSetField(dm, field, NULL, (PetscObject)fe));
46 PetscCall(DMCreateDS(dm));
47 PetscCall(PetscFEDestroy(&fe));
48 PetscCall(DMCreateGlobalVector(dm, &u_f));
49 // Create particle swarm
50 PetscCall(DMCreate(PETSC_COMM_SELF, &sw));
51 PetscCall(DMSetType(sw, DMSWARM));
52 PetscCall(DMSetDimension(sw, dim));
53 PetscCall(DMSwarmSetType(sw, DMSWARM_PIC));
54 PetscCall(DMSwarmSetCellDM(sw, dm));
55 PetscCall(DMSwarmRegisterPetscDatatypeField(sw, "w_q", Nc, PETSC_SCALAR));
56 PetscCall(DMSwarmFinalizeFieldRegister(sw));
57 PetscCall(DMSwarmSetLocalSizes(sw, Np, zero));
58 PetscCall(DMSetFromOptions(sw));
59 PetscCall(DMSwarmGetField(sw, "w_q", &bs, &dtype, (void **)&wq));
60 PetscCall(DMSwarmGetField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords));
61 for (PetscInt p = 0; p < Np; ++p) {
62 coords[p * 2 + 0] = -PetscCosReal((PetscReal)(p + 1) / (PetscReal)(Np + 1) * PETSC_PI);
63 coords[p * 2 + 1] = PetscSinReal((PetscReal)(p + 1) / (PetscReal)(Np + 1) * PETSC_PI);
64 wq[p] = 1.0;
65 energy_0 += wq[p] * (PetscSqr(coords[p * 2 + 0]) + PetscSqr(coords[p * 2 + 1]));
66 }
67 PetscCall(DMSwarmRestoreField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords));
68 PetscCall(DMSwarmRestoreField(sw, "w_q", &bs, &dtype, (void **)&wq));
69 PetscCall(DMSwarmMigrate(sw, removePoints));
70 PetscCall(PetscObjectSetName((PetscObject)sw, "Particle Grid"));
71 PetscCall(DMViewFromOptions(sw, NULL, "-swarm_view"));
72
73 // Project between particles and continuum field
74 const char *fieldnames[1] = {"w_q"};
75 Vec fields[1] = {u_f};
76 PetscCall(DMSwarmProjectFields(sw, NULL, 1, fieldnames, fields, SCATTER_FORWARD));
77 PetscCall(DMSwarmProjectFields(sw, NULL, 1, fieldnames, fields, SCATTER_REVERSE));
78
79 // Compute energy
80 PetscCall(DMSwarmGetField(sw, "w_q", &bs, &dtype, (void **)&wq));
81 PetscCall(DMSwarmGetField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords));
82 for (PetscInt p = 0; p < Np; ++p) energy_1 += wq[p] * (PetscSqr(coords[p * 2 + 0]) + PetscSqr(coords[p * 2 + 1]));
83 PetscCall(DMSwarmRestoreField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords));
84 PetscCall(DMSwarmRestoreField(sw, "w_q", &bs, &dtype, (void **)&wq));
85 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Energy = %20.12e error = %20.12e\n", (double)energy_0, (double)((energy_1 - energy_0) / energy_0)));
86
87 // Cleanup
88 PetscCall(VecDestroy(&u_f));
89 PetscCall(DMDestroy(&sw));
90 PetscCall(DMDestroy(&dm));
91 PetscCall(PetscFinalize());
92 return 0;
93 }
94
95 /*TEST
96
97 build:
98 requires: !complex
99
100 test:
101 suffix: 0
102 requires: double triangle
103 args: -dm_plex_simplex 0 -dm_plex_box_faces 4,2 -dm_plex_box_lower -2.0,0.0 -dm_plex_box_upper 2.0,2.0 \
104 -np 50 -petscspace_degree 2 \
105 -ptof_ksp_type cg -ptof_pc_type ilu -ptof_ksp_rtol 1.e-14 \
106 -ftop_ksp_type lsqr -ftop_pc_type none -ftop_ksp_rtol 1.e-14 \
107 -dm_view -swarm_view
108 filter: grep -v DM_ | grep -v atomic
109
110 test:
111 suffix: bjacobi
112 requires: double triangle
113 args: -dm_plex_simplex 0 -dm_plex_box_faces 4,2 -dm_plex_box_lower -2.0,0.0 -dm_plex_box_upper 2.0,2.0 \
114 -np 50 -petscspace_degree 2 -dm_plex_hash_location \
115 -ptof_ksp_type cg -ptof_pc_type ilu -ptof_ksp_rtol 1.e-14 \
116 -ftop_ksp_type lsqr -ftop_pc_type bjacobi -ftop_sub_pc_type lu -ftop_sub_pc_factor_shift_type nonzero \
117 -dm_view -swarm_view -ftop_ksp_rtol 1.e-14
118 filter: grep -v DM_ | grep -v atomic
119
120 TEST*/
121