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;
26: PetscInitialize(&argc, &argv, NULL, help);
27: /* Create a mesh */
28: DMCreate(PETSC_COMM_WORLD, &dm);
29: DMSetType(dm, DMPLEX);
30: DMSetFromOptions(dm);
31: DMViewFromOptions(dm, NULL, "-dm_view");
33: DMGetDimension(dm, &dim);
34: i = dim;
35: PetscOptionsGetIntArray(NULL, NULL, "-dm_plex_box_faces", faces, &i, NULL);
36: PetscOptionsGetInt(NULL, NULL, "-np", &Np, NULL);
37: DMGetBoundingBox(dm, lo, hi);
38: for (i = 0; i < dim; i++) {
39: h[i] = (hi[i] - lo[i]) / faces[i];
40: 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: PetscFECreateDefault(PETSC_COMM_SELF, dim, Nc, PETSC_FALSE, "", PETSC_DECIDE, &fe);
44: PetscFESetFromOptions(fe);
45: PetscObjectSetName((PetscObject)fe, "fe");
46: DMSetField(dm, field, NULL, (PetscObject)fe);
47: DMCreateDS(dm);
48: PetscFEDestroy(&fe);
49: /* Create particle swarm */
50: DMCreate(PETSC_COMM_SELF, &sw);
51: DMSetType(sw, DMSWARM);
52: DMSetDimension(sw, dim);
53: DMSwarmSetType(sw, DMSWARM_PIC);
54: DMSwarmSetCellDM(sw, dm);
55: DMSwarmRegisterPetscDatatypeField(sw, "w_q", Nc, PETSC_SCALAR);
56: DMSwarmFinalizeFieldRegister(sw);
57: DMSwarmSetLocalSizes(sw, Np, zero);
58: DMSetFromOptions(sw);
59: DMSwarmGetField(sw, "w_q", &bs, &dtype, (void **)&wq);
60: 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: DMSwarmRestoreField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords);
68: DMSwarmRestoreField(sw, "w_q", &bs, &dtype, (void **)&wq);
69: DMSwarmMigrate(sw, removePoints);
70: PetscObjectSetName((PetscObject)sw, "Particle Grid");
71: 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: DMCreateMassMatrix(sw, dm, &M_p);
76: DMCreateGlobalVector(dm, &rho);
77: PetscObjectSetName((PetscObject)rho, "rho");
79: DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f);
80: PetscObjectSetName((PetscObject)f, "weights");
81: MatMultTranspose(M_p, f, rho);
83: /* Visualize mesh field */
84: DMSetOutputSequenceNumber(dm, timestep, time);
85: VecViewFromOptions(rho, NULL, "-rho_view");
87: /* Project field to particles */
88: /* This gives f_p = M_p^+ M f */
89: DMCreateGlobalVector(dm, &rhs);
90: VecCopy(rho, rhs); /* Identity: M^1 M rho */
92: KSPCreate(PETSC_COMM_WORLD, &ksp);
93: KSPSetOptionsPrefix(ksp, "ftop_");
94: KSPSetFromOptions(ksp);
95: KSPGetPC(ksp, &pc);
96: PetscObjectTypeCompare((PetscObject)pc, PCBJACOBI, &is_bjac);
97: if (is_bjac) {
98: DMSwarmCreateMassMatrixSquare(sw, dm, &PM_p);
99: KSPSetOperators(ksp, M_p, PM_p);
100: } else {
101: KSPSetOperators(ksp, M_p, M_p);
102: }
103: KSPSolveTranspose(ksp, rhs, f);
104: KSPDestroy(&ksp);
105: VecDestroy(&rhs);
107: /* Visualize particle field */
108: DMSetOutputSequenceNumber(sw, timestep, time);
109: VecViewFromOptions(f, NULL, "-weights_view");
110: VecNorm(f, NORM_1, &norm);
111: DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f);
113: /* compute energy */
114: DMSwarmGetField(sw, "w_q", &bs, &dtype, (void **)&wq);
115: 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: DMSwarmRestoreField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords);
118: DMSwarmRestoreField(sw, "w_q", &bs, &dtype, (void **)&wq);
119: 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: MatDestroy(&M_p);
122: MatDestroy(&PM_p);
123: VecDestroy(&rho);
124: DMDestroy(&sw);
125: DMDestroy(&dm);
126: 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*/