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