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: KSP ksp;
13: PC pc;
14: Mat M_p, PM_p = NULL;
15: Vec f, rho, rhs;
16: PetscInt dim, Nc = 1, timestep = 0, i, faces[3];
17: PetscInt Np = 10, p, field = 0, zero = 0, bs;
18: PetscReal time = 0.0, norm, energy_0, energy_1;
19: PetscReal lo[3], hi[3], h[3];
20: PetscBool removePoints = PETSC_TRUE;
21: PetscReal *wq, *coords;
22: PetscDataType dtype;
23: PetscBool is_bjac;
25: PetscFunctionBeginUser;
26: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
27: /* Create a mesh */
28: PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
29: PetscCall(DMSetType(dm, DMPLEX));
30: PetscCall(DMSetFromOptions(dm));
31: PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
33: PetscCall(DMGetDimension(dm, &dim));
34: i = dim;
35: PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-dm_plex_box_faces", faces, &i, NULL));
36: PetscCall(PetscOptionsGetInt(NULL, NULL, "-np", &Np, NULL));
37: PetscCall(DMGetBoundingBox(dm, lo, hi));
38: for (i = 0; i < dim; i++) {
39: h[i] = (hi[i] - lo[i]) / faces[i];
40: 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]));
41: }
43: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, Nc, PETSC_FALSE, "", PETSC_DECIDE, &fe));
44: PetscCall(PetscFESetFromOptions(fe));
45: PetscCall(PetscObjectSetName((PetscObject)fe, "fe"));
46: PetscCall(DMSetField(dm, field, NULL, (PetscObject)fe));
47: PetscCall(DMCreateDS(dm));
48: PetscCall(PetscFEDestroy(&fe));
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 (p = 0, energy_0 = 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 particles to field */
74: /* This gives M f = \int_\Omega \phi f, which looks like a rhs for a PDE */
75: PetscCall(DMCreateMassMatrix(sw, dm, &M_p));
76: PetscCall(DMCreateGlobalVector(dm, &rho));
77: PetscCall(PetscObjectSetName((PetscObject)rho, "rho"));
79: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f));
80: PetscCall(PetscObjectSetName((PetscObject)f, "weights"));
81: PetscCall(MatMultTranspose(M_p, f, rho));
83: /* Visualize mesh field */
84: PetscCall(DMSetOutputSequenceNumber(dm, timestep, time));
85: PetscCall(VecViewFromOptions(rho, NULL, "-rho_view"));
87: /* Project field to particles */
88: /* This gives f_p = M_p^+ M f */
89: PetscCall(DMCreateGlobalVector(dm, &rhs));
90: PetscCall(VecCopy(rho, rhs)); /* Identity: M^1 M rho */
92: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
93: PetscCall(KSPSetOptionsPrefix(ksp, "ftop_"));
94: PetscCall(KSPSetFromOptions(ksp));
95: PetscCall(KSPGetPC(ksp, &pc));
96: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCBJACOBI, &is_bjac));
97: if (is_bjac) {
98: PetscCall(DMSwarmCreateMassMatrixSquare(sw, dm, &PM_p));
99: PetscCall(KSPSetOperators(ksp, M_p, PM_p));
100: } else {
101: PetscCall(KSPSetOperators(ksp, M_p, M_p));
102: }
103: PetscCall(KSPSolveTranspose(ksp, rhs, f));
104: PetscCall(KSPDestroy(&ksp));
105: PetscCall(VecDestroy(&rhs));
107: /* Visualize particle field */
108: PetscCall(DMSetOutputSequenceNumber(sw, timestep, time));
109: PetscCall(VecViewFromOptions(f, NULL, "-weights_view"));
110: PetscCall(VecNorm(f, NORM_1, &norm));
111: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f));
113: /* compute energy */
114: PetscCall(DMSwarmGetField(sw, "w_q", &bs, &dtype, (void **)&wq));
115: PetscCall(DMSwarmGetField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords));
116: for (p = 0, energy_1 = 0; p < Np; p++) energy_1 += wq[p] * (PetscSqr(coords[p * 2 + 0]) + PetscSqr(coords[p * 2 + 1]));
117: PetscCall(DMSwarmRestoreField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords));
118: PetscCall(DMSwarmRestoreField(sw, "w_q", &bs, &dtype, (void **)&wq));
119: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Total number = %20.12e. energy = %20.12e error = %20.12e\n", (double)norm, (double)energy_0, (double)((energy_1 - energy_0) / energy_0)));
120: /* Cleanup */
121: PetscCall(MatDestroy(&M_p));
122: PetscCall(MatDestroy(&PM_p));
123: PetscCall(VecDestroy(&rho));
124: PetscCall(DMDestroy(&sw));
125: PetscCall(DMDestroy(&dm));
126: PetscCall(PetscFinalize());
127: return 0;
128: }
130: /*TEST
132: build:
133: requires: !complex
135: test:
136: suffix: 0
137: requires: double triangle
138: args: -dm_plex_simplex 0 -dm_plex_box_faces 4,2 -np 50 -dm_plex_box_lower -2.0,0.0 -dm_plex_box_upper 2.0,2.0 -petscspace_degree 2 -ftop_ksp_type lsqr -ftop_pc_type none -dm_view -swarm_view -ftop_ksp_rtol 1.e-14
139: filter: grep -v DM_ | grep -v atomic
141: test:
142: suffix: bjacobi
143: requires: double triangle
144: args: -dm_plex_simplex 0 -dm_plex_box_faces 4,2 -np 50 -dm_plex_box_lower -2.0,0.0 -dm_plex_box_upper 2.0,2.0 -petscspace_degree 2 -dm_plex_hash_location -ftop_ksp_type lsqr -ftop_pc_type bjacobi -ftop_sub_pc_type lu -ftop_sub_pc_factor_shift_type nonzero -dm_view -swarm_view -ftop_ksp_rtol 1.e-14
145: filter: grep -v DM_ | grep -v atomic
147: TEST*/