Actual source code: ex36.c
1: static char help[] = "First example in homogenization book\n\n";
3: #include <petscsnes.h>
4: #include <petscdmplex.h>
5: #include <petscds.h>
6: #include <petscconvest.h>
7: #include <petscbag.h>
9: /*
10: To control the refinement use -dm_plex_box_faces <n> or -dm_refine <k>, or both
12: To see the exact and computed solutions
14: -compare_view draw -draw_size 500,500 -draw_pause -1
16: To see the delay in convergence of the discretization use
18: -snes_convergence_estimate -convest_num_refine 7 -convest_monitor
20: and to see the proper rate use
22: -dm_refine 5 -snes_convergence_estimate -convest_num_refine 2 -convest_monitor
23: */
25: typedef enum {
26: MOD_CONSTANT,
27: MOD_OSCILLATORY,
28: MOD_TANH,
29: NUM_MOD_TYPES
30: } ModType;
31: const char *modTypes[NUM_MOD_TYPES + 1] = {"constant", "oscillatory", "tanh", "unknown"};
33: /* Constants */
34: enum {
35: EPSILON,
36: NUM_CONSTANTS
37: };
39: typedef struct {
40: PetscReal epsilon; /* Wavelength of fine scale oscillation */
41: } Parameter;
43: typedef struct {
44: PetscBag bag; /* Holds problem parameters */
45: ModType modType; /* Model type */
46: } AppCtx;
48: static PetscErrorCode trig_homogeneous_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
49: {
50: PetscInt d;
51: *u = 1.0;
52: for (d = 0; d < dim; ++d) *u *= PetscSinReal(2.0 * PETSC_PI * x[d]);
53: return PETSC_SUCCESS;
54: }
56: static PetscErrorCode oscillatory_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
57: {
58: Parameter *param = (Parameter *)ctx;
59: const PetscReal eps = param->epsilon;
61: u[0] = x[0] - x[0] * x[0] + (eps / (2. * PETSC_PI)) * (0.5 - x[0]) * PetscSinReal(2. * PETSC_PI * x[0] / eps) + PetscSqr(eps / (2. * PETSC_PI)) * (1. - PetscCosReal(2. * PETSC_PI * x[0] / eps));
62: return PETSC_SUCCESS;
63: }
65: static void f0_trig_homogeneous_u(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, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
66: {
67: PetscInt d;
68: for (d = 0; d < dim; ++d) {
69: PetscScalar v = 1.;
70: for (PetscInt e = 0; e < dim; e++) {
71: if (e == d) {
72: v *= -4.0 * PetscSqr(PETSC_PI) * PetscSinReal(2.0 * PETSC_PI * x[d]);
73: } else {
74: v *= PetscSinReal(2.0 * PETSC_PI * x[d]);
75: }
76: }
77: f0[0] += v;
78: }
79: }
81: static void f1_u(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, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
82: {
83: PetscInt d;
84: for (d = 0; d < dim; ++d) f1[d] = u_x[d];
85: }
87: static void g3_uu(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 g3[])
88: {
89: PetscInt d;
90: for (d = 0; d < dim; ++d) g3[d * dim + d] = 1.0;
91: }
93: static void f0_oscillatory_u(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, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
94: {
95: f0[0] = -1.;
96: }
98: static void f1_oscillatory_u(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, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
99: {
100: const PetscReal eps = PetscRealPart(constants[EPSILON]);
102: f1[0] = u_x[0] / (2. + PetscCosReal(2. * PETSC_PI * x[0] / eps));
103: }
105: static void g3_oscillatory_uu(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 g3[])
106: {
107: const PetscReal eps = PetscRealPart(constants[EPSILON]);
109: g3[0] = 1. / (2. + PetscCosReal(2. * PETSC_PI * x[0] / eps));
110: }
112: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
113: {
114: PetscInt mod;
116: PetscFunctionBeginUser;
117: options->modType = MOD_CONSTANT;
118: PetscOptionsBegin(comm, "", "Homogenization Problem Options", "DMPLEX");
119: mod = options->modType;
120: PetscCall(PetscOptionsEList("-mod_type", "The model type", "ex36.c", modTypes, NUM_MOD_TYPES, modTypes[options->modType], &mod, NULL));
121: options->modType = (ModType)mod;
122: PetscOptionsEnd();
123: PetscFunctionReturn(PETSC_SUCCESS);
124: }
126: static PetscErrorCode SetupParameters(MPI_Comm comm, AppCtx *user)
127: {
128: PetscBag bag;
129: Parameter *p;
131: PetscFunctionBeginUser;
132: PetscCall(PetscBagCreate(comm, sizeof(Parameter), &user->bag));
133: PetscCall(PetscBagGetData(user->bag, (void **)&p));
134: PetscCall(PetscBagSetName(user->bag, "par", "Homogenization parameters"));
135: bag = user->bag;
136: PetscCall(PetscBagRegisterReal(bag, &p->epsilon, 1.0, "epsilon", "Wavelength of fine scale oscillation"));
137: PetscFunctionReturn(PETSC_SUCCESS);
138: }
140: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
141: {
142: PetscFunctionBeginUser;
143: PetscCall(DMCreate(comm, dm));
144: PetscCall(DMSetType(*dm, DMPLEX));
145: PetscCall(DMSetFromOptions(*dm));
146: PetscCall(DMSetApplicationContext(*dm, user));
147: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
148: PetscFunctionReturn(PETSC_SUCCESS);
149: }
151: static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
152: {
153: PetscDS ds;
154: DMLabel label;
155: PetscSimplePointFn *ex;
156: const PetscInt id = 1;
157: void *ctx;
159: PetscFunctionBeginUser;
160: PetscCall(DMGetDS(dm, &ds));
161: PetscCall(PetscBagGetData(user->bag, (void **)&ctx));
162: switch (user->modType) {
163: case MOD_CONSTANT:
164: PetscCall(PetscDSSetResidual(ds, 0, f0_trig_homogeneous_u, f1_u));
165: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu));
166: PetscCall(DMGetLabel(dm, "marker", &label));
167: ex = trig_homogeneous_u;
168: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))ex, NULL, ctx, NULL));
169: break;
170: case MOD_OSCILLATORY:
171: PetscCall(PetscDSSetResidual(ds, 0, f0_oscillatory_u, f1_oscillatory_u));
172: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_oscillatory_uu));
173: PetscCall(DMGetLabel(dm, "marker", &label));
174: ex = oscillatory_u;
175: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))ex, NULL, ctx, NULL));
176: break;
177: default:
178: SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Unsupported model type: %s (%d)", modTypes[PetscMin(user->modType, NUM_MOD_TYPES)], user->modType);
179: }
180: PetscCall(PetscDSSetExactSolution(ds, 0, ex, ctx));
181: /* Setup constants */
182: {
183: Parameter *param;
184: PetscScalar constants[NUM_CONSTANTS];
186: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
188: constants[EPSILON] = param->epsilon;
189: PetscCall(PetscDSSetConstants(ds, NUM_CONSTANTS, constants));
190: }
191: PetscFunctionReturn(PETSC_SUCCESS);
192: }
194: static PetscErrorCode SetupDiscretization(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user)
195: {
196: DM cdm = dm;
197: PetscFE fe;
198: PetscBool simplex;
199: PetscInt dim;
200: char prefix[PETSC_MAX_PATH_LEN];
202: PetscFunctionBeginUser;
203: PetscCall(DMGetDimension(dm, &dim));
204: PetscCall(DMPlexIsSimplex(dm, &simplex));
205: /* Create finite element */
206: PetscCall(PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name));
207: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, name ? prefix : NULL, -1, &fe));
208: PetscCall(PetscObjectSetName((PetscObject)fe, name));
209: /* Set discretization and boundary conditions for each mesh */
210: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
211: PetscCall(DMCreateDS(dm));
212: PetscCall((*setup)(dm, user));
213: while (cdm) {
214: PetscCall(DMCopyDisc(dm, cdm));
215: PetscCall(DMGetCoarseDM(cdm, &cdm));
216: }
217: PetscCall(PetscFEDestroy(&fe));
218: PetscFunctionReturn(PETSC_SUCCESS);
219: }
221: static PetscErrorCode CompareView(Vec u)
222: {
223: DM dm;
224: Vec v[2], lv[2], exact;
225: PetscOptions options;
226: PetscViewer viewer;
227: PetscViewerFormat format;
228: PetscBool flg;
229: PetscInt i;
230: const char *name, *prefix;
232: PetscFunctionBeginUser;
233: PetscCall(VecGetDM(u, &dm));
234: PetscCall(PetscObjectGetOptions((PetscObject)dm, &options));
235: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
236: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)dm), options, prefix, "-compare_view", &viewer, &format, &flg));
237: if (flg) {
238: PetscCall(DMGetGlobalVector(dm, &exact));
239: PetscCall(DMComputeExactSolution(dm, 0.0, exact, NULL));
240: v[0] = u;
241: v[1] = exact;
242: for (i = 0; i < 2; ++i) {
243: PetscCall(DMGetLocalVector(dm, &lv[i]));
244: PetscCall(PetscObjectGetName((PetscObject)v[i], &name));
245: PetscCall(PetscObjectSetName((PetscObject)lv[i], name));
246: PetscCall(DMGlobalToLocalBegin(dm, v[i], INSERT_VALUES, lv[i]));
247: PetscCall(DMGlobalToLocalEnd(dm, v[i], INSERT_VALUES, lv[i]));
248: PetscCall(DMPlexInsertBoundaryValues(dm, PETSC_TRUE, lv[i], 0., NULL, NULL, NULL));
249: }
250: PetscCall(DMPlexVecView1D(dm, 2, lv, viewer));
251: for (i = 0; i < 2; ++i) PetscCall(DMRestoreLocalVector(dm, &lv[i]));
252: PetscCall(DMRestoreGlobalVector(dm, &exact));
253: PetscCall(PetscViewerDestroy(&viewer));
254: }
255: PetscFunctionReturn(PETSC_SUCCESS);
256: }
258: typedef struct {
259: Mat Mcoarse; /* Mass matrix on the coarse space */
260: Mat Mfine; /* Mass matrix on the fine space */
261: Mat Ifine; /* Interpolator from coarse to fine */
262: Vec Iscale; /* Scaling vector for restriction */
263: KSP kspCoarse; /* Solver for the coarse mass matrix */
264: Vec tmpfine; /* Temporary vector in the fine space */
265: Vec tmpcoarse; /* Temporary vector in the coarse space */
266: } ProjStruct;
268: static PetscErrorCode DestroyCoarseProjection(Mat Pi)
269: {
270: ProjStruct *ctx;
272: PetscFunctionBegin;
273: PetscCall(MatShellGetContext(Pi, (void **)&ctx));
274: PetscCall(MatDestroy(&ctx->Mcoarse));
275: PetscCall(MatDestroy(&ctx->Mfine));
276: PetscCall(MatDestroy(&ctx->Ifine));
277: PetscCall(VecDestroy(&ctx->Iscale));
278: PetscCall(KSPDestroy(&ctx->kspCoarse));
279: PetscCall(VecDestroy(&ctx->tmpcoarse));
280: PetscCall(VecDestroy(&ctx->tmpfine));
281: PetscCall(PetscFree(ctx));
282: PetscCall(MatShellSetContext(Pi, NULL));
283: PetscFunctionReturn(PETSC_SUCCESS);
284: }
286: static PetscErrorCode CoarseProjection(Mat Pi, Vec x, Vec y)
287: {
288: ProjStruct *ctx;
290: PetscFunctionBegin;
291: PetscCall(MatShellGetContext(Pi, (void **)&ctx));
292: PetscCall(MatMult(ctx->Mfine, x, ctx->tmpfine));
293: PetscCall(PetscObjectSetName((PetscObject)ctx->tmpfine, "Fine DG RHS"));
294: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ctx->tmpfine, "fine_dg_"));
295: PetscCall(VecViewFromOptions(ctx->tmpfine, NULL, "-rhs_view"));
296: PetscCall(MatMultTranspose(ctx->Ifine, ctx->tmpfine, ctx->tmpcoarse));
297: PetscCall(PetscObjectSetName((PetscObject)ctx->tmpcoarse, "Coarse DG RHS"));
298: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ctx->tmpcoarse, "coarse_dg_"));
299: PetscCall(VecViewFromOptions(ctx->tmpcoarse, NULL, "-rhs_view"));
300: PetscCall(VecPointwiseMult(ctx->tmpcoarse, ctx->Iscale, ctx->tmpcoarse));
301: PetscCall(VecViewFromOptions(ctx->tmpcoarse, NULL, "-rhs_view"));
302: PetscCall(KSPSolve(ctx->kspCoarse, ctx->tmpcoarse, y));
303: PetscFunctionReturn(PETSC_SUCCESS);
304: }
306: static PetscErrorCode CreateCoarseProjection(DM dmc, DM dmf, Mat *Pi)
307: {
308: ProjStruct *ctx;
309: PetscInt m, n, M, N;
311: PetscFunctionBegin;
312: PetscCall(PetscMalloc1(1, &ctx));
313: PetscCall(DMCreateGlobalVector(dmc, &ctx->tmpcoarse));
314: PetscCall(DMCreateGlobalVector(dmf, &ctx->tmpfine));
315: PetscCall(VecGetLocalSize(ctx->tmpcoarse, &m));
316: PetscCall(VecGetSize(ctx->tmpcoarse, &M));
317: PetscCall(VecGetLocalSize(ctx->tmpfine, &n));
318: PetscCall(VecGetSize(ctx->tmpfine, &N));
319: PetscCall(DMCreateMassMatrix(dmc, dmc, &ctx->Mcoarse));
320: PetscCall(DMCreateMassMatrix(dmf, dmf, &ctx->Mfine));
321: PetscCall(DMCreateInterpolation(dmc, dmf, &ctx->Ifine, &ctx->Iscale));
322: PetscCall(KSPCreate(PetscObjectComm((PetscObject)dmc), &ctx->kspCoarse));
323: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ctx->kspCoarse, "coarse_"));
324: PetscCall(KSPSetOperators(ctx->kspCoarse, ctx->Mcoarse, ctx->Mcoarse));
325: PetscCall(KSPSetFromOptions(ctx->kspCoarse));
326: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)dmc), m, n, M, N, ctx, Pi));
327: PetscCall(MatShellSetOperation(*Pi, MATOP_DESTROY, (void (*)(void))DestroyCoarseProjection));
328: PetscCall(MatShellSetOperation(*Pi, MATOP_MULT, (void (*)(void))CoarseProjection));
329: PetscFunctionReturn(PETSC_SUCCESS);
330: }
332: typedef struct {
333: Mat Ifdg; /* Embed the fine space into the DG version */
334: Mat Pi; /* The L_2 stable projection to the DG coarse space */
335: Vec tmpc; /* A temporary vector in the DG coarse space */
336: Vec tmpf; /* A temporary vector in the DG fine space */
337: } QuasiInterp;
339: static PetscErrorCode DestroyQuasiInterpolator(Mat P)
340: {
341: QuasiInterp *ctx;
343: PetscFunctionBegin;
344: PetscCall(MatShellGetContext(P, (void **)&ctx));
345: PetscCall(MatDestroy(&ctx->Ifdg));
346: PetscCall(MatDestroy(&ctx->Pi));
347: PetscCall(VecDestroy(&ctx->tmpc));
348: PetscCall(VecDestroy(&ctx->tmpf));
349: PetscCall(PetscFree(ctx));
350: PetscCall(MatShellSetContext(P, NULL));
351: PetscFunctionReturn(PETSC_SUCCESS);
352: }
354: static PetscErrorCode QuasiInterpolate(Mat P, Vec x, Vec y)
355: {
356: QuasiInterp *ctx;
357: DM dmcdg, dmc;
358: Vec ly;
360: PetscFunctionBegin;
361: PetscCall(MatShellGetContext(P, (void **)&ctx));
362: PetscCall(MatMult(ctx->Ifdg, x, ctx->tmpf));
364: PetscCall(PetscObjectSetName((PetscObject)ctx->tmpf, "Fine DG Potential"));
365: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ctx->tmpf, "fine_dg_"));
366: PetscCall(VecViewFromOptions(ctx->tmpf, NULL, "-vec_view"));
367: PetscCall(MatMult(ctx->Pi, x, ctx->tmpc));
369: PetscCall(PetscObjectSetName((PetscObject)ctx->tmpc, "Coarse DG Potential"));
370: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ctx->tmpc, "coarse_dg_"));
371: PetscCall(VecViewFromOptions(ctx->tmpc, NULL, "-vec_view"));
372: PetscCall(VecGetDM(ctx->tmpc, &dmcdg));
374: PetscCall(VecGetDM(y, &dmc));
375: PetscCall(DMGetLocalVector(dmc, &ly));
376: PetscCall(DMPlexComputeClementInterpolant(dmcdg, ctx->tmpc, ly));
377: PetscCall(DMLocalToGlobalBegin(dmc, ly, INSERT_VALUES, y));
378: PetscCall(DMLocalToGlobalEnd(dmc, ly, INSERT_VALUES, y));
379: PetscCall(DMRestoreLocalVector(dmc, &ly));
380: PetscFunctionReturn(PETSC_SUCCESS);
381: }
383: static PetscErrorCode CreateQuasiInterpolator(DM dmc, DM dmf, Mat *P)
384: {
385: QuasiInterp *ctx;
386: DM dmcdg, dmfdg;
387: PetscFE fe;
388: Vec x, y;
389: DMPolytopeType ct;
390: PetscInt dim, cStart, m, n, M, N;
392: PetscFunctionBegin;
393: PetscCall(PetscCalloc1(1, &ctx));
394: PetscCall(DMGetGlobalVector(dmc, &x));
395: PetscCall(DMGetGlobalVector(dmf, &y));
396: PetscCall(VecGetLocalSize(x, &m));
397: PetscCall(VecGetSize(x, &M));
398: PetscCall(VecGetLocalSize(y, &n));
399: PetscCall(VecGetSize(y, &N));
400: PetscCall(DMRestoreGlobalVector(dmc, &x));
401: PetscCall(DMRestoreGlobalVector(dmf, &y));
403: PetscCall(DMClone(dmf, &dmfdg));
404: PetscCall(DMGetDimension(dmfdg, &dim));
405: PetscCall(DMPlexGetHeightStratum(dmfdg, 0, &cStart, NULL));
406: PetscCall(DMPlexGetCellType(dmfdg, cStart, &ct));
407: PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 1, ct, "fine_dg_", PETSC_DETERMINE, &fe));
408: PetscCall(DMSetField(dmfdg, 0, NULL, (PetscObject)fe));
409: PetscCall(PetscFEDestroy(&fe));
410: PetscCall(DMCreateDS(dmfdg));
411: PetscCall(DMCreateInterpolation(dmf, dmfdg, &ctx->Ifdg, NULL));
412: PetscCall(DMCreateGlobalVector(dmfdg, &ctx->tmpf));
413: PetscCall(DMDestroy(&dmfdg));
415: PetscCall(DMClone(dmc, &dmcdg));
416: PetscCall(DMGetDimension(dmcdg, &dim));
417: PetscCall(DMPlexGetHeightStratum(dmcdg, 0, &cStart, NULL));
418: PetscCall(DMPlexGetCellType(dmcdg, cStart, &ct));
419: PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 1, ct, "coarse_dg_", PETSC_DETERMINE, &fe));
420: PetscCall(DMSetField(dmcdg, 0, NULL, (PetscObject)fe));
421: PetscCall(PetscFEDestroy(&fe));
422: PetscCall(DMCreateDS(dmcdg));
424: PetscCall(CreateCoarseProjection(dmcdg, dmf, &ctx->Pi));
425: PetscCall(DMCreateGlobalVector(dmcdg, &ctx->tmpc));
426: PetscCall(DMDestroy(&dmcdg));
428: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)dmc), m, n, M, N, ctx, P));
429: PetscCall(MatShellSetOperation(*P, MATOP_DESTROY, (void (*)(void))DestroyQuasiInterpolator));
430: PetscCall(MatShellSetOperation(*P, MATOP_MULT, (void (*)(void))QuasiInterpolate));
431: PetscFunctionReturn(PETSC_SUCCESS);
432: }
434: static PetscErrorCode CoarseTest(DM dm, Vec u, AppCtx *user)
435: {
436: DM dmc;
437: Mat P; /* The quasi-interpolator to the coarse space */
438: Vec uc;
440: PetscFunctionBegin;
441: if (user->modType == MOD_CONSTANT) PetscFunctionReturn(PETSC_SUCCESS);
442: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &dmc));
443: PetscCall(DMSetType(dmc, DMPLEX));
444: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)dmc, "coarse_"));
445: PetscCall(DMSetApplicationContext(dmc, user));
446: PetscCall(DMSetFromOptions(dmc));
447: PetscCall(DMViewFromOptions(dmc, NULL, "-dm_view"));
449: PetscCall(SetupDiscretization(dmc, "potential", SetupPrimalProblem, user));
450: PetscCall(DMCreateGlobalVector(dmc, &uc));
451: PetscCall(PetscObjectSetName((PetscObject)uc, "potential"));
452: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)uc, "coarse_"));
454: PetscCall(CreateQuasiInterpolator(dmc, dm, &P));
455: #if 1
456: PetscCall(MatMult(P, u, uc));
457: #else
458: {
459: Mat In;
460: Vec sc;
462: PetscCall(DMCreateInterpolation(dmc, dm, &In, &sc));
463: PetscCall(MatMultTranspose(In, u, uc));
464: PetscCall(VecPointwiseMult(uc, sc, uc));
465: PetscCall(MatDestroy(&In));
466: PetscCall(VecDestroy(&sc));
467: }
468: #endif
469: PetscCall(CompareView(uc));
471: PetscCall(MatDestroy(&P));
472: PetscCall(VecDestroy(&uc));
473: PetscCall(DMDestroy(&dmc));
474: PetscFunctionReturn(PETSC_SUCCESS);
475: }
477: int main(int argc, char **argv)
478: {
479: DM dm; /* Problem specification */
480: SNES snes; /* Nonlinear solver */
481: Vec u; /* Solutions */
482: AppCtx user; /* User-defined work context */
484: PetscFunctionBeginUser;
485: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
486: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
487: PetscCall(SetupParameters(PETSC_COMM_WORLD, &user));
488: /* Primal system */
489: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
490: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
491: PetscCall(SNESSetDM(snes, dm));
492: PetscCall(SetupDiscretization(dm, "potential", SetupPrimalProblem, &user));
493: PetscCall(DMCreateGlobalVector(dm, &u));
494: PetscCall(VecSet(u, 0.0));
495: PetscCall(PetscObjectSetName((PetscObject)u, "potential"));
496: PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user));
497: PetscCall(SNESSetFromOptions(snes));
498: PetscCall(DMSNESCheckFromOptions(snes, u));
499: PetscCall(SNESSolve(snes, NULL, u));
500: PetscCall(SNESGetSolution(snes, &u));
501: PetscCall(VecViewFromOptions(u, NULL, "-potential_view"));
502: PetscCall(CompareView(u));
503: /* Looking at a coarse problem */
504: PetscCall(CoarseTest(dm, u, &user));
505: /* Cleanup */
506: PetscCall(VecDestroy(&u));
507: PetscCall(SNESDestroy(&snes));
508: PetscCall(DMDestroy(&dm));
509: PetscCall(PetscBagDestroy(&user.bag));
510: PetscCall(PetscFinalize());
511: return 0;
512: }
514: /*TEST
516: test:
517: suffix: 1d_p1_constant
518: args: -dm_plex_dim 1 -dm_plex_box_faces 4 -potential_petscspace_degree 1 -dmsnes_check
520: test:
521: suffix: 1d_p1_constant_conv
522: args: -dm_plex_dim 1 -dm_plex_box_faces 4 -potential_petscspace_degree 1 \
523: -snes_convergence_estimate -convest_num_refine 2
525: test:
526: suffix: 1d_p1_oscillatory
527: args: -mod_type oscillatory -epsilon 0.03125 \
528: -dm_plex_dim 1 -dm_plex_box_faces 4 -potential_petscspace_degree 1 -dm_refine 2 -dmsnes_check \
529: -coarse_dm_plex_dim 1 -coarse_dm_plex_box_faces 4 -coarse_dm_plex_hash_location \
530: -fine_dg_petscspace_degree 1 -fine_dg_petscdualspace_lagrange_continuity 0 \
531: -coarse_dg_petscspace_degree 1 -coarse_dg_petscdualspace_lagrange_continuity 0
533: TEST*/