Actual source code: ex1.c

  1: static char help[] = "Tests projection with DMSwarm using general particle shapes.\n";

  3: #include <petsc/private/petscfeimpl.h>
  4: #include <petscdmswarm.h>
  5: #include <petscdmplex.h>
  6: #include <petscds.h>
  7: #include <petscksp.h>

  9: typedef struct {
 10:   PetscReal L[3]; /* Dimensions of the mesh bounding box */
 11: } AppCtx;

 13: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
 14: {
 15:   PetscReal low[3], high[3];
 16:   PetscInt  dim;

 18:   PetscFunctionBeginUser;
 19:   PetscCall(DMCreate(comm, dm));
 20:   PetscCall(DMSetType(*dm, DMPLEX));
 21:   PetscCall(DMSetFromOptions(*dm));
 22:   PetscCall(DMGetBoundingBox(*dm, low, high));
 23:   PetscCall(DMGetDimension(*dm, &dim));
 24:   for (PetscInt d = 0; d < dim; ++d) user->L[d] = high[d] - low[d];
 25:   PetscCall(PetscObjectSetName((PetscObject)*dm, "Mesh"));
 26:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
 27:   PetscFunctionReturn(PETSC_SUCCESS);
 28: }

 30: static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *a_ctx)
 31: {
 32:   AppCtx  *ctx = (AppCtx *)a_ctx;
 33:   PetscInt d;
 34:   u[0] = 0.0;
 35:   for (d = 0; d < dim; ++d) u[0] += x[d] / ctx->L[d];
 36:   return PETSC_SUCCESS;
 37: }

 39: static void identity(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
 40: {
 41:   g0[0] = 1.0;
 42: }

 44: static PetscErrorCode CreateFEM(DM dm, AppCtx *user)
 45: {
 46:   PetscFE   fe;
 47:   PetscDS   ds;
 48:   PetscBool simplex;
 49:   PetscInt  dim;

 51:   PetscFunctionBeginUser;
 52:   PetscCall(DMGetDimension(dm, &dim));
 53:   PetscCall(DMPlexIsSimplex(dm, &simplex));
 54:   PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, 1, simplex, NULL, -1, &fe));
 55:   PetscCall(PetscObjectSetName((PetscObject)fe, "fe"));
 56:   PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
 57:   PetscCall(DMCreateDS(dm));
 58:   PetscCall(PetscFEDestroy(&fe));
 59:   /* Setup to form mass matrix */
 60:   PetscCall(DMGetDS(dm, &ds));
 61:   PetscCall(PetscDSSetJacobian(ds, 0, 0, identity, NULL, NULL, NULL));
 62:   PetscFunctionReturn(PETSC_SUCCESS);
 63: }

 65: static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
 66: {
 67:   PetscDS          prob;
 68:   PetscFE          fe;
 69:   PetscQuadrature  quad;
 70:   PetscScalar     *vals;
 71:   PetscReal       *v0, *J, *invJ, detJ, *coords, *xi0;
 72:   PetscInt        *cellid;
 73:   const PetscReal *qpoints, *qweights;
 74:   PetscInt         cStart, cEnd, c, Nq, q, dim;

 76:   PetscFunctionBeginUser;
 77:   PetscCall(DMGetDimension(dm, &dim));
 78:   PetscCall(DMGetDS(dm, &prob));
 79:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
 80:   PetscCall(PetscDSGetDiscretization(prob, 0, (PetscObject *)&fe));
 81:   PetscCall(PetscFEGetQuadrature(fe, &quad));
 82:   PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, &qpoints, &qweights));
 83:   PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw));
 84:   PetscCall(DMSetType(*sw, DMSWARM));
 85:   PetscCall(DMSetDimension(*sw, dim));
 86:   PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC));
 87:   PetscCall(DMSwarmSetCellDM(*sw, dm));
 88:   PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "f_q", 1, PETSC_REAL));
 89:   PetscCall(DMSwarmFinalizeFieldRegister(*sw));
 90:   PetscCall(DMSwarmSetLocalSizes(*sw, (cEnd - cStart) * Nq, 0));
 91:   PetscCall(DMSetFromOptions(*sw));
 92:   PetscCall(PetscMalloc4(dim, &xi0, dim, &v0, dim * dim, &J, dim * dim, &invJ));
 93:   for (c = 0; c < dim; c++) xi0[c] = -1.;
 94:   PetscCall(DMSwarmGetField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
 95:   PetscCall(DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid));
 96:   PetscCall(DMSwarmGetField(*sw, "f_q", NULL, NULL, (void **)&vals));
 97:   for (c = cStart; c < cEnd; ++c) {
 98:     for (q = 0; q < Nq; ++q) {
 99:       PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
100:       cellid[c * Nq + q] = c;
101:       CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q * dim], &coords[(c * Nq + q) * dim]);
102:       PetscCall(linear(dim, 0.0, &coords[(c * Nq + q) * dim], 1, &vals[c * Nq + q], (void *)user));
103:       vals[c * Nq + q] *= qweights[q] * detJ;
104:     }
105:   }
106:   PetscCall(DMSwarmRestoreField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
107:   PetscCall(DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid));
108:   PetscCall(DMSwarmRestoreField(*sw, "f_q", NULL, NULL, (void **)&vals));
109:   PetscCall(PetscFree4(xi0, v0, J, invJ));
110:   PetscCall(DMSwarmMigrate(*sw, PETSC_FALSE));
111:   PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
112:   PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view"));
113:   PetscFunctionReturn(PETSC_SUCCESS);
114: }

116: static PetscErrorCode TestL2Projection(DM dm, DM sw, AppCtx *user)
117: {
118:   PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
119:   KSP       ksp;
120:   Mat       mass;
121:   Vec       u, rhs, uproj;
122:   void     *ctxs[1];
123:   PetscReal error;

125:   PetscFunctionBeginUser;
126:   ctxs[0]  = (void *)user;
127:   funcs[0] = linear;

129:   PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "f_q", &u));
130:   PetscCall(VecViewFromOptions(u, NULL, "-f_view"));
131:   PetscCall(DMGetGlobalVector(dm, &rhs));
132:   PetscCall(DMCreateMassMatrix(sw, dm, &mass));
133:   PetscCall(MatMultTranspose(mass, u, rhs));
134:   PetscCall(VecViewFromOptions(rhs, NULL, "-rhs_view"));
135:   PetscCall(MatDestroy(&mass));
136:   PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "f_q", &u));
137:   PetscCall(DMGetGlobalVector(dm, &uproj));
138:   PetscCall(DMCreateMatrix(dm, &mass));
139:   PetscCall(DMPlexSNESComputeJacobianFEM(dm, uproj, mass, mass, user));
140:   PetscCall(MatViewFromOptions(mass, NULL, "-mass_mat_view"));
141:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
142:   PetscCall(KSPSetOperators(ksp, mass, mass));
143:   PetscCall(KSPSetFromOptions(ksp));
144:   PetscCall(KSPSolve(ksp, rhs, uproj));
145:   PetscCall(KSPDestroy(&ksp));
146:   PetscCall(PetscObjectSetName((PetscObject)uproj, "Full Projection"));
147:   PetscCall(VecViewFromOptions(uproj, NULL, "-proj_vec_view"));
148:   PetscCall(DMComputeL2Diff(dm, 0.0, funcs, ctxs, uproj, &error));
149:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Projected L2 Error: %g\n", (double)error));
150:   PetscCall(MatDestroy(&mass));
151:   PetscCall(DMRestoreGlobalVector(dm, &rhs));
152:   PetscCall(DMRestoreGlobalVector(dm, &uproj));
153:   PetscFunctionReturn(PETSC_SUCCESS);
154: }

156: int main(int argc, char *argv[])
157: {
158:   MPI_Comm comm;
159:   DM       dm, sw;
160:   AppCtx   user;

162:   PetscFunctionBeginUser;
163:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
164:   comm = PETSC_COMM_WORLD;
165:   PetscCall(CreateMesh(comm, &dm, &user));
166:   PetscCall(CreateFEM(dm, &user));
167:   PetscCall(CreateParticles(dm, &sw, &user));
168:   PetscCall(TestL2Projection(dm, sw, &user));
169:   PetscCall(DMDestroy(&dm));
170:   PetscCall(DMDestroy(&sw));
171:   PetscCall(PetscFinalize());
172:   return 0;
173: }

175: /*TEST

177:   build:
178:     requires: double !complex
179:   test:
180:     suffix: proj_0
181:     requires: triangle
182:     args: -dm_plex_box_faces 2,2\
183:           -dm_view -sw_view\
184:           -petscspace_degree 1 -petscfe_default_quadrature_order 1\
185:           -pc_type lu
186:   test:
187:     suffix: proj_1
188:     requires: triangle
189:     args: -dm_plex_simplex 0 -dm_plex_box_faces 2,2\
190:           -dm_view -sw_view\
191:           -petscspace_degree 1 -petscfe_default_quadrature_order 1\
192:           -pc_type lu
193:   test:
194:     suffix: proj_2
195:     requires: ctetgen
196:     args: -dm_plex_dim 3 -dm_plex_box_faces 2,2,2\
197:           -dm_view -sw_view\
198:           -petscspace_degree 1 -petscfe_default_quadrature_order 1\
199:           -pc_type lu
200:   test:
201:     suffix: proj_3
202:     requires: ctetgen
203:     args: -dm_plex_dim 3 -dm_plex_simplex 0\
204:           -dm_plex_box_faces 2,2,2\
205:           -dm_view -sw_view\
206:           -petscspace_degree 1 -petscfe_default_quadrature_order 1\
207:           -pc_type lu

209: TEST*/