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: {
19: DMCreate(comm, dm);
20: DMSetType(*dm, DMPLEX);
21: DMSetFromOptions(*dm);
22: DMViewFromOptions(*dm, NULL, "-dm_view");
23: return(0);
24: }
26: static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
27: {
28: PetscInt d;
29: u[0] = 0.0;
30: for (d = 0; d < dim; ++d) u[0] += x[d];
31: return 0;
32: }
34: static void identity(PetscInt dim, PetscInt Nf, PetscInt NfAux,
35: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
36: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
37: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
38: {
39: g0[0] = 1.0;
40: }
42: static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
43: {
44: PetscDS prob;
45: PetscFE fe;
46: PetscQuadrature quad;
47: PetscScalar *vals;
48: PetscReal *v0, *J, *invJ, detJ, *coords, *xi0;
49: PetscInt *cellid;
50: const PetscReal *qpoints;
51: PetscInt Ncell, c, Nq, q, dim;
52: PetscBool simplex;
53: PetscErrorCode ierr;
56: DMGetDimension(dm, &dim);
57: DMPlexIsSimplex(dm, &simplex);
58: PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, 1, simplex, NULL, -1, &fe);
59: DMGetDS(dm, &prob);
60: PetscDSSetDiscretization(prob, 0, (PetscObject) fe);
61: PetscFEDestroy(&fe);
62: PetscDSSetJacobian(prob, 0, 0, identity, NULL, NULL, NULL);
63: DMPlexGetHeightStratum(dm, 0, NULL, &Ncell);
64: PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
65: PetscFEGetQuadrature(fe, &quad);
66: PetscQuadratureGetData(quad, NULL, NULL, &Nq, &qpoints, NULL);
68: DMCreate(PetscObjectComm((PetscObject) dm), sw);
69: DMSetType(*sw, DMSWARM);
70: DMSetDimension(*sw, dim);
72: DMSwarmSetType(*sw, DMSWARM_PIC);
73: DMSwarmSetCellDM(*sw, dm);
74: DMSwarmRegisterPetscDatatypeField(*sw, "f_q", 1, PETSC_SCALAR);
75: DMSwarmFinalizeFieldRegister(*sw);
76: DMSwarmSetLocalSizes(*sw, Ncell * Nq, 0);
77: DMSetFromOptions(*sw);
79: PetscMalloc4(dim, &xi0, dim, &v0, dim*dim, &J, dim*dim, &invJ);
80: for (c = 0; c < dim; c++) xi0[c] = -1.;
81: DMSwarmGetField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
82: DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
83: DMSwarmGetField(*sw, "f_q", NULL, NULL, (void **) &vals);
84: for (c = 0; c < Ncell; ++c) {
85: for (q = 0; q < Nq; ++q) {
86: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
87: cellid[c*Nq + q] = c;
88: CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], &coords[(c*Nq + q)*dim]);
89: linear(dim, 0.0, &coords[(c*Nq + q)*dim], 1, &vals[c*Nq + q], NULL);
90: }
91: }
92: DMSwarmRestoreField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
93: DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
94: DMSwarmRestoreField(*sw, "f_q", NULL, NULL, (void **) &vals);
95: PetscFree4(xi0, v0, J, invJ);
96: return(0);
97: }
99: static PetscErrorCode TestL2Projection(DM dm, DM sw, AppCtx *user)
100: {
101: PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *);
102: KSP ksp;
103: Mat mass;
104: Vec u, rhs, uproj;
105: PetscReal error;
106: PetscErrorCode ierr;
109: funcs[0] = linear;
111: DMSwarmCreateGlobalVectorFromField(sw, "f_q", &u);
112: VecViewFromOptions(u, NULL, "-f_view");
113: DMGetGlobalVector(dm, &rhs);
114: DMCreateMassMatrix(sw, dm, &mass);
115: MatMult(mass, u, rhs);
116: MatDestroy(&mass);
117: VecDestroy(&u);
119: DMGetGlobalVector(dm, &uproj);
120: DMCreateMatrix(dm, &mass);
121: DMPlexSNESComputeJacobianFEM(dm, uproj, mass, mass, user);
122: MatViewFromOptions(mass, NULL, "-mass_mat_view");
123: KSPCreate(PETSC_COMM_WORLD, &ksp);
124: KSPSetOperators(ksp, mass, mass);
125: KSPSetFromOptions(ksp);
126: KSPSolve(ksp, rhs, uproj);
127: KSPDestroy(&ksp);
128: PetscObjectSetName((PetscObject) uproj, "Full Projection");
129: VecViewFromOptions(uproj, NULL, "-proj_vec_view");
130: DMComputeL2Diff(dm, 0.0, funcs, NULL, uproj, &error);
131: PetscPrintf(PETSC_COMM_WORLD, "Projected L2 Error: %g\n", (double) error);
132: MatDestroy(&mass);
133: DMRestoreGlobalVector(dm, &rhs);
134: DMRestoreGlobalVector(dm, &uproj);
135: return(0);
136: }
138: int main (int argc, char * argv[]) {
139: MPI_Comm comm;
140: DM dm, sw;
141: AppCtx user;
144: PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
145: comm = PETSC_COMM_WORLD;
146: CreateMesh(comm, &dm, &user);
147: CreateParticles(dm, &sw, &user);
148: PetscObjectSetName((PetscObject) dm, "Mesh");
149: DMViewFromOptions(dm, NULL, "-dm_view");
150: DMViewFromOptions(sw, NULL, "-sw_view");
151: TestL2Projection(dm, sw, &user);
152: DMDestroy(&dm);
153: DMDestroy(&sw);
154: PetscFinalize();
155: return 0;
156: }
158: /*TEST
160: test:
161: suffix: proj_0
162: requires: pragmatic
163: TODO: broken
164: args: -dm_plex_separate_marker 0 -dm_view -sw_view -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu
165: test:
166: suffix: proj_1
167: requires: pragmatic
168: TODO: broken
169: 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
170: test:
171: suffix: proj_2
172: requires: pragmatic
173: TODO: broken
174: 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
175: test:
176: suffix: proj_3
177: requires: pragmatic
178: TODO: broken
179: 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
181: TEST*/