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*/