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