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,
 33:                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 34:                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 35:                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
 36: {
 37:   g0[0] = 1.0;
 38: }

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

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

 65:   DMCreate(PetscObjectComm((PetscObject) dm), sw);
 66:   DMSetType(*sw, DMSWARM);
 67:   DMSetDimension(*sw, dim);

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

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

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

105:   funcs[0] = linear;

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

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

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

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

153: /*TEST

155:   test:
156:     suffix: proj_0
157:     requires: pragmatic
158:     TODO: broken
159:     args: -dm_plex_separate_marker 0 -dm_view -sw_view -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu -dm_adaptor pragmatic
160:   test:
161:     suffix: proj_1
162:     requires: pragmatic
163:     TODO: broken
164:     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
165:   test:
166:     suffix: proj_2
167:     requires: pragmatic
168:     TODO: broken
169:     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
170:   test:
171:     suffix: proj_3
172:     requires: pragmatic
173:     TODO: broken
174:     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

176: TEST*/