Actual source code: ex1.c
1: static char help[] = "Example program demonstrating projection between particle and finite element spaces\n\n";
3: #include <petscdmplex.h>
4: #include <petscds.h>
5: #include <petscdmswarm.h>
6: #include <petscksp.h>
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;
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"));
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"));
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));
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)));
87: // Cleanup
88: PetscCall(VecDestroy(&u_f));
89: PetscCall(DMDestroy(&sw));
90: PetscCall(DMDestroy(&dm));
91: PetscCall(PetscFinalize());
92: return 0;
93: }
95: /*TEST
97: build:
98: requires: !complex
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
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
120: TEST*/