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