Actual source code: ex1.c

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

  3: #include <petsc/private/dmswarmimpl.h>
  4: #include <petsc/private/petscfeimpl.h>

  6: #include <petscdmplex.h>
  7: #include <petscds.h>
  8: #include <petscksp.h>

 10: typedef struct {
 11:   PetscInt dummy;
 12: } AppCtx;

 14: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
 15: {
 17:   DMCreate(comm, dm);
 18:   DMSetType(*dm, DMPLEX);
 19:   DMSetFromOptions(*dm);
 20:   DMViewFromOptions(*dm, NULL, "-dm_view");
 21:   return 0;
 22: }

 24: static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 25: {
 26:   PetscInt d;
 27:   u[0] = 0.0;
 28:   for (d = 0; d < dim; ++d) u[0] += x[d];
 29:   return 0;
 30: }

 32: 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[])
 33: {
 34:   g0[0] = 1.0;
 35: }

 37: static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
 38: {
 39:   PetscDS          prob;
 40:   PetscFE          fe;
 41:   PetscQuadrature  quad;
 42:   PetscScalar     *vals;
 43:   PetscReal       *v0, *J, *invJ, detJ, *coords, *xi0;
 44:   PetscInt        *cellid;
 45:   const PetscReal *qpoints;
 46:   PetscInt         Ncell, c, Nq, q, dim;
 47:   PetscBool        simplex;

 50:   DMGetDimension(dm, &dim);
 51:   DMPlexIsSimplex(dm, &simplex);
 52:   PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, 1, simplex, NULL, -1, &fe);
 53:   DMGetDS(dm, &prob);
 54:   PetscDSSetDiscretization(prob, 0, (PetscObject)fe);
 55:   PetscFEDestroy(&fe);
 56:   PetscDSSetJacobian(prob, 0, 0, identity, NULL, NULL, NULL);
 57:   DMPlexGetHeightStratum(dm, 0, NULL, &Ncell);
 58:   PetscDSGetDiscretization(prob, 0, (PetscObject *)&fe);
 59:   PetscFEGetQuadrature(fe, &quad);
 60:   PetscQuadratureGetData(quad, NULL, NULL, &Nq, &qpoints, NULL);

 62:   DMCreate(PetscObjectComm((PetscObject)dm), sw);
 63:   DMSetType(*sw, DMSWARM);
 64:   DMSetDimension(*sw, dim);

 66:   DMSwarmSetType(*sw, DMSWARM_PIC);
 67:   DMSwarmSetCellDM(*sw, dm);
 68:   DMSwarmRegisterPetscDatatypeField(*sw, "f_q", 1, PETSC_SCALAR);
 69:   DMSwarmFinalizeFieldRegister(*sw);
 70:   DMSwarmSetLocalSizes(*sw, Ncell * Nq, 0);
 71:   DMSetFromOptions(*sw);

 73:   PetscMalloc4(dim, &xi0, dim, &v0, dim * dim, &J, dim * dim, &invJ);
 74:   for (c = 0; c < dim; c++) xi0[c] = -1.;
 75:   DMSwarmGetField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
 76:   DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid);
 77:   DMSwarmGetField(*sw, "f_q", NULL, NULL, (void **)&vals);
 78:   for (c = 0; c < Ncell; ++c) {
 79:     for (q = 0; q < Nq; ++q) {
 80:       DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
 81:       cellid[c * Nq + q] = c;
 82:       CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q * dim], &coords[(c * Nq + q) * dim]);
 83:       linear(dim, 0.0, &coords[(c * Nq + q) * dim], 1, &vals[c * Nq + q], NULL);
 84:     }
 85:   }
 86:   DMSwarmRestoreField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
 87:   DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid);
 88:   DMSwarmRestoreField(*sw, "f_q", NULL, NULL, (void **)&vals);
 89:   PetscFree4(xi0, v0, J, invJ);
 90:   return 0;
 91: }

 93: static PetscErrorCode TestL2Projection(DM dm, DM sw, AppCtx *user)
 94: {
 95:   PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
 96:   KSP       ksp;
 97:   Mat       mass;
 98:   Vec       u, rhs, uproj;
 99:   PetscReal error;

102:   funcs[0] = linear;

104:   DMSwarmCreateGlobalVectorFromField(sw, "f_q", &u);
105:   VecViewFromOptions(u, NULL, "-f_view");
106:   DMGetGlobalVector(dm, &rhs);
107:   DMCreateMassMatrix(sw, dm, &mass);
108:   MatMult(mass, u, rhs);
109:   MatDestroy(&mass);
110:   VecDestroy(&u);

112:   DMGetGlobalVector(dm, &uproj);
113:   DMCreateMatrix(dm, &mass);
114:   DMPlexSNESComputeJacobianFEM(dm, uproj, mass, mass, user);
115:   MatViewFromOptions(mass, NULL, "-mass_mat_view");
116:   KSPCreate(PETSC_COMM_WORLD, &ksp);
117:   KSPSetOperators(ksp, mass, mass);
118:   KSPSetFromOptions(ksp);
119:   KSPSolve(ksp, rhs, uproj);
120:   KSPDestroy(&ksp);
121:   PetscObjectSetName((PetscObject)uproj, "Full Projection");
122:   VecViewFromOptions(uproj, NULL, "-proj_vec_view");
123:   DMComputeL2Diff(dm, 0.0, funcs, NULL, uproj, &error);
124:   PetscPrintf(PETSC_COMM_WORLD, "Projected L2 Error: %g\n", (double)error);
125:   MatDestroy(&mass);
126:   DMRestoreGlobalVector(dm, &rhs);
127:   DMRestoreGlobalVector(dm, &uproj);
128:   return 0;
129: }

131: int main(int argc, char *argv[])
132: {
133:   MPI_Comm comm;
134:   DM       dm, sw;
135:   AppCtx   user;

138:   PetscInitialize(&argc, &argv, NULL, help);
139:   comm = PETSC_COMM_WORLD;
140:   CreateMesh(comm, &dm, &user);
141:   CreateParticles(dm, &sw, &user);
142:   PetscObjectSetName((PetscObject)dm, "Mesh");
143:   DMViewFromOptions(dm, NULL, "-dm_view");
144:   DMViewFromOptions(sw, NULL, "-sw_view");
145:   TestL2Projection(dm, sw, &user);
146:   DMDestroy(&dm);
147:   DMDestroy(&sw);
148:   PetscFinalize();
149:   return 0;
150: }

152: /*TEST

154:   test:
155:     suffix: proj_0
156:     requires: pragmatic
157:     TODO: broken
158:     args: -dm_plex_separate_marker 0 -dm_view -sw_view -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu -dm_adaptor pragmatic
159:   test:
160:     suffix: proj_1
161:     requires: pragmatic
162:     TODO: broken
163:     args: -dm_plex_simplex 0 -dm_plex_separate_marker 0 -dm_view -sw_view -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu -dm_adaptor pragmatic
164:   test:
165:     suffix: proj_2
166:     requires: pragmatic
167:     TODO: broken
168:     args: -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_view -sw_view -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu -dm_adaptor pragmatic
169:   test:
170:     suffix: proj_3
171:     requires: pragmatic
172:     TODO: broken
173:     args: -dm_plex_simplex 0 -dm_plex_separate_marker 0 -dm_view -sw_view -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu -dm_adaptor pragmatic

175: TEST*/