Actual source code: ex1.c
petsc-3.14.6 2021-03-30
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 dim; /* The topological mesh dimension */
12: PetscBool simplex; /* Flag for simplices or tensor cells */
13: char mshNam[PETSC_MAX_PATH_LEN]; /* Name of the mesh filename if any */
14: PetscInt nbrVerEdge; /* Number of vertices per edge if unit square/cube generated */
15: } AppCtx;
17: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
18: {
22: options->dim = 2;
23: options->simplex = PETSC_TRUE;
24: PetscStrcpy(options->mshNam, "");
26: PetscOptionsBegin(comm, "", "Meshing Adaptation Options", "DMPLEX");
27: PetscOptionsInt("-dim", "The topological mesh dimension", "ex1.c", options->dim, &options->dim, NULL);
28: PetscOptionsBool("-simplex", "The flag for simplices or tensor cells", "ex1.c", options->simplex, &options->simplex, NULL);
29: PetscOptionsString("-msh", "Name of the mesh filename if any", "ex1.c", options->mshNam, options->mshNam, sizeof(options->mshNam), NULL);
30: PetscOptionsInt("-nbrVerEdge", "Number of vertices per edge if unit square/cube generated", "ex1.c", options->nbrVerEdge, &options->nbrVerEdge, NULL);
31: PetscOptionsEnd();
33: return(0);
34: }
36: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
37: {
38: PetscBool flag;
42: PetscStrcmp(user->mshNam, "", &flag);
43: if (flag) {
44: PetscInt faces[3];
46: faces[0] = user->nbrVerEdge-1; faces[1] = user->nbrVerEdge-1; faces[2] = user->nbrVerEdge-1;
47: DMPlexCreateBoxMesh(comm, user->dim, user->simplex, faces, NULL, NULL, NULL, PETSC_TRUE, dm);
48: } else {
49: DMPlexCreateFromFile(comm, user->mshNam, PETSC_TRUE, dm);
50: DMGetDimension(*dm, &user->dim);
51: }
52: {
53: DM distributedMesh = NULL;
55: /* Distribute mesh over processes */
56: DMPlexDistribute(*dm, 0, NULL, &distributedMesh);
57: if (distributedMesh) {
58: DMDestroy(dm);
59: *dm = distributedMesh;
60: }
61: }
62: DMSetFromOptions(*dm);
63: return(0);
64: }
66: static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
67: {
68: PetscInt d;
69: u[0] = 0.0;
70: for (d = 0; d < dim; ++d) u[0] += x[d];
71: return 0;
72: }
74: static void identity(PetscInt dim, PetscInt Nf, PetscInt NfAux,
75: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
76: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
77: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
78: {
79: g0[0] = 1.0;
80: }
82: static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
83: {
84: PetscDS prob;
85: PetscFE fe;
86: PetscQuadrature quad;
87: PetscScalar *vals;
88: PetscReal *v0, *J, *invJ, detJ, *coords, *xi0;
89: PetscInt *cellid;
90: const PetscReal *qpoints;
91: PetscInt Ncell, c, Nq, q, dim;
92: PetscErrorCode ierr;
95: DMGetDimension(dm, &dim);
96: PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, 1, user->simplex, NULL, -1, &fe);
97: DMGetDS(dm, &prob);
98: PetscDSSetDiscretization(prob, 0, (PetscObject) fe);
99: PetscFEDestroy(&fe);
100: PetscDSSetJacobian(prob, 0, 0, identity, NULL, NULL, NULL);
101: DMPlexGetHeightStratum(dm, 0, NULL, &Ncell);
102: PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
103: PetscFEGetQuadrature(fe, &quad);
104: PetscQuadratureGetData(quad, NULL, NULL, &Nq, &qpoints, NULL);
106: DMCreate(PetscObjectComm((PetscObject) dm), sw);
107: DMSetType(*sw, DMSWARM);
108: DMSetDimension(*sw, dim);
110: DMSwarmSetType(*sw, DMSWARM_PIC);
111: DMSwarmSetCellDM(*sw, dm);
112: DMSwarmRegisterPetscDatatypeField(*sw, "f_q", 1, PETSC_SCALAR);
113: DMSwarmFinalizeFieldRegister(*sw);
114: DMSwarmSetLocalSizes(*sw, Ncell * Nq, 0);
115: DMSetFromOptions(*sw);
117: PetscMalloc4(dim, &xi0, dim, &v0, dim*dim, &J, dim*dim, &invJ);
118: for (c = 0; c < dim; c++) xi0[c] = -1.;
119: DMSwarmGetField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
120: DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
121: DMSwarmGetField(*sw, "f_q", NULL, NULL, (void **) &vals);
122: for (c = 0; c < Ncell; ++c) {
123: for (q = 0; q < Nq; ++q) {
124: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
125: cellid[c*Nq + q] = c;
126: CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], &coords[(c*Nq + q)*dim]);
127: linear(dim, 0.0, &coords[(c*Nq + q)*dim], 1, &vals[c*Nq + q], NULL);
128: }
129: }
130: DMSwarmRestoreField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
131: DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
132: DMSwarmRestoreField(*sw, "f_q", NULL, NULL, (void **) &vals);
133: PetscFree4(xi0, v0, J, invJ);
134: return(0);
135: }
137: static PetscErrorCode TestL2Projection(DM dm, DM sw, AppCtx *user)
138: {
139: PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *);
140: KSP ksp;
141: Mat mass;
142: Vec u, rhs, uproj;
143: PetscReal error;
144: PetscErrorCode ierr;
147: funcs[0] = linear;
149: DMSwarmCreateGlobalVectorFromField(sw, "f_q", &u);
150: VecViewFromOptions(u, NULL, "-f_view");
151: DMGetGlobalVector(dm, &rhs);
152: DMCreateMassMatrix(sw, dm, &mass);
153: MatMult(mass, u, rhs);
154: MatDestroy(&mass);
155: VecDestroy(&u);
157: DMGetGlobalVector(dm, &uproj);
158: DMCreateMatrix(dm, &mass);
159: DMPlexSNESComputeJacobianFEM(dm, uproj, mass, mass, user);
160: MatViewFromOptions(mass, NULL, "-mass_mat_view");
161: KSPCreate(PETSC_COMM_WORLD, &ksp);
162: KSPSetOperators(ksp, mass, mass);
163: KSPSetFromOptions(ksp);
164: KSPSolve(ksp, rhs, uproj);
165: KSPDestroy(&ksp);
166: PetscObjectSetName((PetscObject) uproj, "Full Projection");
167: VecViewFromOptions(uproj, NULL, "-proj_vec_view");
168: DMComputeL2Diff(dm, 0.0, funcs, NULL, uproj, &error);
169: PetscPrintf(PETSC_COMM_WORLD, "Projected L2 Error: %g\n", (double) error);
170: MatDestroy(&mass);
171: DMRestoreGlobalVector(dm, &rhs);
172: DMRestoreGlobalVector(dm, &uproj);
173: return(0);
174: }
176: int main (int argc, char * argv[]) {
177: MPI_Comm comm;
178: DM dm, sw;
179: AppCtx user;
182: PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
183: comm = PETSC_COMM_WORLD;
184: ProcessOptions(comm, &user);
185: CreateMesh(comm, &dm, &user);
186: CreateParticles(dm, &sw, &user);
187: PetscObjectSetName((PetscObject) dm, "Mesh");
188: DMViewFromOptions(dm, NULL, "-dm_view");
189: DMViewFromOptions(sw, NULL, "-sw_view");
190: TestL2Projection(dm, sw, &user);
191: DMDestroy(&dm);
192: DMDestroy(&sw);
193: PetscFinalize();
194: return 0;
195: }
197: /*TEST
199: test:
200: suffix: proj_0
201: requires: pragmatic
202: TODO: broken
203: args: -dim 2 -nbrVerEdge 3 -dm_plex_separate_marker 0 -dm_view -sw_view -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu
204: test:
205: suffix: proj_1
206: requires: pragmatic
207: TODO: broken
208: args: -dim 2 -simplex 0 -nbrVerEdge 3 -dm_plex_separate_marker 0 -dm_view -sw_view -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu
209: test:
210: suffix: proj_2
211: requires: pragmatic
212: TODO: broken
213: args: -dim 3 -nbrVerEdge 3 -dm_view -sw_view -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu
214: test:
215: suffix: proj_3
216: requires: pragmatic
217: TODO: broken
218: args: -dim 2 -simplex 0 -nbrVerEdge 3 -dm_plex_separate_marker 0 -dm_view -sw_view -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu
220: TEST*/